1 /* copyright(c) 1991-2012 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public
6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public
15 license along with the gnu c library;
if not, see
16 <http://www.gnu.org/licenses/>. */
19 /* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27 /* we
do support the iec 559 math functionality,
real and complex. */
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
32 /* we
do not support c11 <threads.h>. */
34 c this file is created by hand from taumain.f
35 c actions: remove routines: taudem dectes taufil filhep
36 c add: inietc will not necesarily work fine ...
38 c rename iniphy to iniphx
40 SUBROUTINE inietc(jakk1,jakk2,itd,ifpho)
42 COMMON / taurad / xk0dec,itdkrc
43 DOUBLE PRECISION xk0dec
44 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
45 COMMON /phoact/ ifphot
47 c kto=1 will denote tau+, thus :: idff=-15
51 c radiative correction switch in tau --> e(mu) decays
53 c switches of tau+ tau- decay modes
56 c photos activation switch
60 SUBROUTINE tralo4(KTOS,PHOI,PHOF,AM)
64 COMMON / momdec / q1,q2,p1,p2,p3,p4
66 double precision q1(4),q2(4),p1(4),p2(4),p3(4),p4(4),p1qq(4),p2qq(4)
67 double precision pin(4),pout(4),pbst(4),pbs1(4),qq(4),pi
68 double precision thet,phi,exe
69 REAL*4 phoi(4),phof(4)
71 DATA pi /3.141592653589793238462643d0/
73 $ (phoi(4)**2-phoi(3)**2-phoi(2)**2-phoi(1)**2))
85 ELSEIF(idtra.EQ.2)
THEN
90 ELSEIF(idtra.EQ.3)
THEN
104 CALL bostdq(1,qq,pbst,pbst)
105 CALL bostdq(1,qq,p1,p1qq)
106 CALL bostdq(1,qq,p2,p2qq)
108 pbs1(3)=sqrt(pbst(3)**2+pbst(2)**2+pbst(1)**2)
111 exe=(pbs1(4)+pbs1(3))/dsqrt(pbs1(4)**2-pbs1(3)**2)
112 c for ktos=1 boost is antiparallel to 4-momentum of p2.
113 c restframes of tau+ tau- and
'first' frame of
'higgs' are all connected
114 c by boosts along z axis
115 IF(ktos.EQ.1) exe=(pbs1(4)-pbs1(3))/dsqrt(pbs1(4)**2-pbs1(3)**2)
116 CALL bostd3(exe,pin,pout)
118 c once in z/gamma/higgs rest frame we control further kinematics by p2qq for ktos=1,2
119 thet=acos(p2qq(3)/sqrt(p2qq(3)**2+p2qq(2)**2+p2qq(1)**2))
121 phi=acos(p2qq(1)/sqrt(p2qq(2)**2+p2qq(1)**2))
122 IF(p2qq(2).LT.0d0) phi=2*pi-phi
124 CALL rotpox(thet,phi,pout)
125 CALL bostdq(-1,qq,pout,pout)
132 SUBROUTINE choice(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
133 $ amrx,gamrx,amra,gamra,amrb,gamrb)
134 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
135 * ,ampiz,ampi,amro,gamro,ama1,gama1
136 * ,amk,amkz,amkst,gamkst
138 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
139 * ,ampiz,ampi,amro,gamro,ama1,gama1
140 * ,amk,amkz,amkst,gamkst
146 c xxxxa correspond to s2 channel
156 ELSEIF(mnum.EQ.1)
THEN
165 ELSEIF(mnum.EQ.2)
THEN
174 ELSEIF(mnum.EQ.3)
THEN
183 ELSEIF(mnum.EQ.4)
THEN
192 ELSEIF(mnum.EQ.5)
THEN
201 ELSEIF(mnum.EQ.6)
THEN
210 ELSEIF(mnum.EQ.7)
THEN
219 ELSEIF(mnum.EQ.8)
THEN
228 ELSEIF(mnum.EQ.101)
THEN
237 ELSEIF(mnum.EQ.102)
THEN
257 IF (rr.LE.prob1)
THEN
259 ELSEIF(rr.LE.(prob1+prob2))
THEN
274 prob3=1.0-prob1-prob2
277 * ----------------------------------------------------------------------
278 * initialisation of tau decay parameters and routines
281 * ----------------------------------------------------------------------
283 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
284 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
285 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
286 * ,ampiz,ampi,amro,gamro,ama1,gama1
287 * ,amk,amkz,amkst,gamkst
289 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
290 * ,ampiz,ampi,amro,gamro,ama1,gama1
291 * ,amk,amkz,amkst,gamkst
292 COMMON / taubra / gamprt(30),jlist(30),nchan
293 COMMON / taukle / bra1,brk0,brk0b,brks
294 REAL*4 bra1,brk0,brk0b,brks
295 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
296 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
298 CHARACTER names(nmode)*31
299 CHARACTER oldnames(7)*31
302 $ bxinit =
'(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
307 * list of branching ratios
308 cam normalised to e nu nutau channel
309 cam enu munu pinu rhonu a1nu knu k*nu pi
310 cam
DATA jlist / 1, 2, 3, 4, 5, 6, 7,
311 *am
DATA gamprt /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,o.o811,0.616
315 * conventions of particles names
316 * k-,p-,k+, k0,p-,kb, k-,p0,k0
317 * 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
318 * p0,p0,k-, k-,p-,p+, p-,kb,p0
319 * 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
324 dimension nopik(6,nmode),npik(nmode)
325 *am outgoing multiplicity and flavors of multi-pion /multi-k modes
334 DATA nopik / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
335 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
336 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
337 3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
338 4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
339 5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
340 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
341 c ajwmod fix sign bug, 2/22/99
342 7 -3,-4, 0, 0, 0, 0 /
343 * list of branching ratios
348 IF(i.EQ. 1) gamprt(i) =0.1800
349 IF(i.EQ. 2) gamprt(i) =0.1751
350 IF(i.EQ. 3) gamprt(i) =0.1110
351 IF(i.EQ. 4) gamprt(i) =0.2515
352 IF(i.EQ. 5) gamprt(i) =0.1790
353 IF(i.EQ. 6) gamprt(i) =0.0071
354 IF(i.EQ. 7) gamprt(i) =0.0134
355 IF(i.EQ. 8) gamprt(i) =0.0450
356 IF(i.EQ. 9) gamprt(i) =0.0100
357 IF(i.EQ.10) gamprt(i) =0.0009
358 IF(i.EQ.11) gamprt(i) =0.0004
359 IF(i.EQ.12) gamprt(i) =0.0003
360 IF(i.EQ.13) gamprt(i) =0.0005
361 IF(i.EQ.14) gamprt(i) =0.0015
362 IF(i.EQ.15) gamprt(i) =0.0015
363 IF(i.EQ.16) gamprt(i) =0.0015
364 IF(i.EQ.17) gamprt(i) =0.0005
365 IF(i.EQ.18) gamprt(i) =0.0050
366 IF(i.EQ.19) gamprt(i) =0.0055
367 IF(i.EQ.20) gamprt(i) =0.0017
368 IF(i.EQ.21) gamprt(i) =0.0013
369 IF(i.EQ.22) gamprt(i) =0.0010
370 IF(i.EQ. 1) oldnames(i)=
' TAU- --> E- '
371 IF(i.EQ. 2) oldnames(i)=
' TAU- --> MU- '
372 IF(i.EQ. 3) oldnames(i)=
' TAU- --> PI- '
373 IF(i.EQ. 4) oldnames(i)=
' TAU- --> PI-, PI0 '
374 IF(i.EQ. 5) oldnames(i)=
' TAU- --> A1- (two subch) '
375 IF(i.EQ. 6) oldnames(i)=
' TAU- --> K- '
376 IF(i.EQ. 7) oldnames(i)=
' TAU- --> K*- (two subch) '
377 IF(i.EQ. 8) names(i-7)=
' TAU- --> 2PI-, PI0, PI+ '
378 IF(i.EQ. 9) names(i-7)=
' TAU- --> 3PI0, PI- '
379 IF(i.EQ.10) names(i-7)=
' TAU- --> 2PI-, PI+, 2PI0 '
380 IF(i.EQ.11) names(i-7)=
' TAU- --> 3PI-, 2PI+, '
381 IF(i.EQ.12) names(i-7)=
' TAU- --> 3PI-, 2PI+, PI0 '
382 IF(i.EQ.13) names(i-7)=
' TAU- --> 2PI-, PI+, 3PI0 '
383 IF(i.EQ.14) names(i-7)=
' TAU- --> K-, PI-, K+ '
384 IF(i.EQ.15) names(i-7)=
' TAU- --> K0, PI-, K0B '
385 IF(i.EQ.16) names(i-7)=
' TAU- --> K-, K0, PI0 '
386 IF(i.EQ.17) names(i-7)=
' TAU- --> PI0 PI0 K- '
387 IF(i.EQ.18) names(i-7)=
' TAU- --> K- PI- PI+ '
388 IF(i.EQ.19) names(i-7)=
' TAU- --> PI- K0B PI0 '
389 IF(i.EQ.20) names(i-7)=
' TAU- --> ETA PI- PI0 '
390 IF(i.EQ.21) names(i-7)=
' TAU- --> PI- PI0 GAM '
391 IF(i.EQ.22) names(i-7)=
' TAU- --> K- K0 '
400 idffin(j,i)=nopik(j,i)
405 * --- coefficients to fix ratio of:
406 * --- a1 3charged/ a1 1charged 2 neutrals matrix elements(masless lim.)
407 * --- probability of k0 to be ks
408 * --- probability of k0b to be ks
409 * --- ratio of coefficients for k*--> k0 pi-
410 * --- all coefficents should be in the range(0.0,1.0)
411 * --- they meaning is probability of the first choice only
IF one
412 * --- neglects mass-phase space effects
426 * zw 13.04.89 here was an error
427 scabib = sqrt(1.-ccabib**2)
429 gamel = gfermi**2*amtau**5/(192*pi**3)
435 FUNCTION dcdmas(IDENT)
436 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
437 * ,ampiz,ampi,amro,gamro,ama1,gama1
438 * ,amk,amkz,amkst,gamkst
440 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
441 * ,ampiz,ampi,amro,gamro,ama1,gama1
442 * ,amk,amkz,amkst,gamkst
443 IF (ident.EQ. 1)
THEN
445 ELSEIF (ident.EQ.-1)
THEN
447 ELSEIF (ident.EQ. 2)
THEN
449 ELSEIF (ident.EQ.-2)
THEN
451 ELSEIF (ident.EQ. 3)
THEN
453 ELSEIF (ident.EQ.-3)
THEN
455 ELSEIF (ident.EQ. 4)
THEN
457 ELSEIF (ident.EQ.-4)
THEN
459 ELSEIF (ident.EQ. 8)
THEN
461 ELSEIF (ident.EQ.-8)
THEN
463 ELSEIF (ident.EQ. 9)
THEN
465 ELSEIF (ident.EQ.-9)
THEN
468 print *,
'STOP IN APKMAS, WRONG IDENT=',ident
473 FUNCTION lunpik(ID,ISGN)
474 COMMON / taukle / bra1,brk0,brk0b,brks
475 REAL*4 bra1,brk0,brk0b,brks
478 IF (ident.EQ. 1)
THEN
480 ELSEIF (ident.EQ.-1)
THEN
482 ELSEIF (ident.EQ. 2)
THEN
484 ELSEIF (ident.EQ.-2)
THEN
486 ELSEIF (ident.EQ. 3)
THEN
488 ELSEIF (ident.EQ.-3)
THEN
490 ELSEIF (ident.EQ. 4)
THEN
492 * k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
494 IF (xio(1).GT.brk0)
THEN
499 ELSEIF (ident.EQ.-4)
THEN
501 * k0b--> k0_long(is 130) / k0_short(is 310) = 1/1
503 IF (xio(1).GT.brk0b)
THEN
508 ELSEIF (ident.EQ. 8)
THEN
510 ELSEIF (ident.EQ.-8)
THEN
512 ELSEIF (ident.EQ. 9)
THEN
514 ELSEIF (ident.EQ.-9)
THEN
517 print *,
'STOP IN IPKDEF, WRONG IDENT=',ident
525 SUBROUTINE taurdf(KTO)
526 c this routine can be called before any tau+ or tau- event is generated
527 c it can be used to generate tau+ and tau- samples of different
529 COMMON / taukle / bra1,brk0,brk0b,brks
530 REAL*4 bra1,brk0,brk0b,brks
531 COMMON / taubra / gamprt(30),jlist(30),nchan
534 c ajwmod: set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
541 c ajwmod: set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
550 SUBROUTINE iniphx(XK00)
551 * ----------------------------------------------------------------------
552 * initialisation of parameters
553 * used in qed and/or gsw routines
554 * ----------------------------------------------------------------------
555 COMMON / qedprm /alfinv,alfpi,xk0
556 REAL*8 alfinv,alfpi,xk0
559 pi8 = 4.d0*datan(1.d0)
561 alfpi = 1d0/(alfinv*pi8)
566 c ----------------------------------------------------------------------
567 c initialisation of masses
570 c ----------------------------------------------------------------------
571 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
572 * ,ampiz,ampi,amro,gamro,ama1,gama1
573 * ,amk,amkz,amkst,gamkst
575 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
576 * ,ampiz,ampi,amro,gamro,ama1,gama1
577 * ,amk,amkz,amkst,gamkst
579 c in-coming / out-going fermion masses
581 c --- tau mass must be the same as in the host program, what-so-ever
589 * masses used in tau decays
603 c in-coming / out-going fermion masses
608 c masses used in tau decays cleo settings
623 subroutine bostdq(idir,vv,pp,q)
624 * *******************************
625 c boost along arbitrary vector v(see eg. j.d. jacson, classical
627 c four-vector pp is boosted from an actual frame to the rest frame
628 c of the four-vector v(for idir=1) or back(for idir=-1).
629 c q is a resulting four-vector.
630 c note: v must be time-like, pp may be arbitrary.
632 c written by: wieslaw placzek date: 22.07.1994
633 c last update: 3/29/95 by: m.s.
635 implicit DOUBLE PRECISION (a-h,o-z)
637 DOUBLE PRECISION v(4),p(4),q(4),pp(4),vv(4)
643 amv=(v(4)**2-v(1)**2-v(2)**2-v(3)**2)
645 write(6,*)
'bosstv: warning amv**2=',amv
649 q(4)=( p(1)*v(1)+p(2)*v(2)+p(3)*v(3)+p(4)*v(4))/amv
650 wsp =(q(4)+p(4))/(v(4)+amv)
651 elseif (idir.eq.1)
then
652 q(4)=(-p(1)*v(1)-p(2)*v(2)-p(3)*v(3)+p(4)*v(4))/amv
653 wsp =-(q(4)+p(4))/(v(4)+amv)
655 write(nout,*)
' >>> boostv: wrong value of idir = ',idir