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>. */
36 c note that the routines are not like in cpc deck this is historical
37 c=======================================================================
38 c====================== dectes : test of tau decay library===========
39 c====================== ktory = 1 :
INTERFACE of koral-z
TYPE ==========
40 c====================== ktory = 2 :
INTERFACE of koral-b
TYPE =========
41 c=======================================================================
42 c
COMMON /pawc/ blan(10000)
43 COMMON / / blan(10000)
45 COMMON / inout / inut,iout
51 OPEN(iout,file=
"./tauola.output")
52 OPEN(inut,file=
"./dane.dat")
58 SUBROUTINE dectes(KTORY)
59 c ************************
61 DOUBLE PRECISION hh(4)
62 c switches for tauola;
63 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
66 COMMON / inout / inut,iout
67 c lund
TYPE identifier for a1
69 c /ptau/ is used in routine tralo4
71 COMMON / taurad / xk0dec,itdkrc
74 c special switch for tests of dgamma/dq**2 in a1 decay
75 c keya1=1 constant width of a1 and rho
76 c keya1=2 free choice of rho propagator(defined in
function fpik)
77 c and free choice of a1 mass and width.
function g(Q**2)
78 c(see formula 3.48 in comp. phys. comm. 64 (1991) 275)
79 c hard coded both in monte carlo and in testing distribution.
80 c keya1=3
function g(Q**2) hardcoded in the Monte Carlo
81 c(it is timy to calculate
82 c testing distribution.
83 c-----------------------------------------------------------------------
85 c-----------------------------------------------------------------------
86 c======================================
94 READ( ninp,3000) testit
95 WRITE(nout,3000) testit
96 READ( ninp,3001) kat1,kat2,kat3,kat4,kat5,kat6
97 READ( ninp,3002) nevt,jak1,jak2,itdkrc
98 READ( ninp,3003) ptau,xk0dec
100 c======================================
102 WRITE(nout,
'(6A6/6I6)')
103 $
'KAT1',
'KAT2',
'KAT3',
'KAT4',
'KAT5',
'KAT6',
104 $ kat1 , kat2 , kat3 , kat4 , kat5 , kat6
105 WRITE(nout,
'(4A12/4I12)')
106 $
'NEVT',
'JAK1',
'JAK2',
'ITDKRC',
107 $ nevt, jak1 , jak2 , itdkrc
108 WRITE(nout,
'(2A12/2F12.6)')
111 c======================================
115 c lund identifier(for tau+) -15
121 c kto=1 denotes tau defined by idff(i.e. tau+)
122 c kto=2 denotes the opposite(i.e. tau-)
125 print *,
'for the sake of these tests KTO has to be 2'
126 print *,
'to change tau- to tau+ change IDFF from -15 to 15'
129 c tau polarization in its restframe;
133 c tau momentum in gev;
135 c number of events to be generated;
138 print *,
'NEVTES= ',nevtes
139 WRITE(iout,7011) keya1
142 WRITE(iout,7001) jak,idff,pol(3),ptau
144 WRITE(iout,7004) jak,idff,pol(3),ptau
146 c initialisation of tau decay package tauola
147 c ******************************************
158 c-----------------------------------------------------------------------
160 c-----------------------------------------------------------------------
164 c reslu initialise the lund record
171 CALL dekay(kto+10,hh)
183 IF(ipri.EQ.1) print *,
' event no: ',nev,
' NEVTES: ',nevtes
186 c-----------------------------------------------------------------------
188 c-----------------------------------------------------------------------
195 7001
FORMAT(//4(/1x,15(5h=====))
196 $ /,
' ', 19x,
' TEST OF RAD. CORR IN ELECTRON DECAY ',9x,1h ,
197 $ /,
' ', 19x,
' TESTS OF TAU DECAY ROUTINES ',9x,1h ,
198 $ /,
' ', 19x,
' INTERFACE OF THE KORAL-Z TYPE ',9x,1h ,
199 $ 2(/,1x,15(5h=====)),
200 $ /,5x ,
'JAK =',i7 ,
' KEY DEFINING DECAY TYPE ',9x,1h ,
201 $ /,5x ,
'IDFF =',i7 ,
' LUND IDENTIFIER FOR FIRST TAU ',9x,1h ,
202 $ /,5x ,
'POL(3)=',f7.2,
' THIRD COMPONENT OF TAU POLARIZ. ',9x,1h ,
203 $ /,5x ,
'PTAU =',f7.2,
' THIRD COMPONENT OF TAU MOM. GEV ',9x,1h ,
204 $ 2(/,1x,15(5h=====))/)
205 7002
FORMAT(///1x,
'===== EVENT NO.',i4,1x,5h=====)
206 7003
FORMAT(5x,
'POLARIMETRIC VECTOR: ',
207 $ 7x,
'HH(1)',7x,
'HH(2)',7x,
'HH(3)',7x,
'HH(4)',
208 $ /, 5x,
' ', 4(1x,f11.8) )
209 7004
FORMAT(//4(/1x,15(5h=====))
210 $ /,
' ', 19x,
' TEST OF RAD. CORR IN ELECTRON DECAY ',9x,1h ,
211 $ /,
' ', 19x,
' TESTS OF TAU DECAY ROUTINES ',9x,1h ,
212 $ /,
' ', 19x,
' INTERFACE OF THE KORAL-B TYPE ',9x,1h ,
213 $ 2(/,1x,15(5h=====)),
214 $ /,5x ,
'JAK =',i7 ,
' KEY DEFINING DECAY TYPE ',9x,1h ,
215 $ /,5x ,
'IDFF =',i7 ,
' LUND IDENTIFIER FOR FIRST TAU ',9x,1h ,
216 $ /,5x ,
'POL(3)=',f7.2,
' THIRD COMPONENT OF TAU POLARIZ. ',9x,1h ,
217 $ /,5x ,
'PTAU =',f7.2,
' THIRD COMPONENT OF TAU MOM. GEV ',9x,1h ,
218 $ 2(/,1x,15(5h=====))/)
219 7011
FORMAT(///1x,
'===== TYPE OF CURRENT',i4,1x,5h=====)
221 SUBROUTINE choice(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
222 $ amrx,gamrx,amra,gamra,amrb,gamrb)
223 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
224 * ,ampiz,ampi,amro,gamro,ama1,gama1
225 * ,amk,amkz,amkst,gamkst
227 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
228 * ,ampiz,ampi,amro,gamro,ama1,gama1
229 * ,amk,amkz,amkst,gamkst
235 c xxxxa correspond to s2 channel
245 ELSEIF(mnum.EQ.1)
THEN
254 ELSEIF(mnum.EQ.2)
THEN
263 ELSEIF(mnum.EQ.3)
THEN
272 ELSEIF(mnum.EQ.4)
THEN
281 ELSEIF(mnum.EQ.5)
THEN
290 ELSEIF(mnum.EQ.6)
THEN
299 ELSEIF(mnum.EQ.7)
THEN
308 ELSEIF(mnum.EQ.8)
THEN
317 ELSEIF(mnum.EQ.101)
THEN
326 ELSEIF(mnum.EQ.102)
THEN
346 IF (rr.LE.prob1)
THEN
348 ELSEIF(rr.LE.(prob1+prob2))
THEN
363 prob3=1.0-prob1-prob2
366 * ----------------------------------------------------------------------
367 * initialisation of tau decay parameters and routines
370 * ----------------------------------------------------------------------
372 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
373 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
374 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
375 * ,ampiz,ampi,amro,gamro,ama1,gama1
376 * ,amk,amkz,amkst,gamkst
378 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
379 * ,ampiz,ampi,amro,gamro,ama1,gama1
380 * ,amk,amkz,amkst,gamkst
381 COMMON / taubra / gamprt(30),jlist(30),nchan
382 COMMON / taukle / bra1,brk0,brk0b,brks
383 REAL*4 bra1,brk0,brk0b,brks
384 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
385 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
387 CHARACTER names(nmode)*31
388 CHARACTER oldnames(7)*31
391 $ bxinit =
'(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
396 * list of branching ratios
397 cam normalised to e nu nutau channel
398 cam enu munu pinu rhonu a1nu knu k*nu pi
399 cam
DATA jlist / 1, 2, 3, 4, 5, 6, 7,
400 *am
DATA gamprt /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,o.o811,0.616
404 * conventions of particles names
405 * k-,p-,k+, k0,p-,kb, k-,p0,k0
406 * 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
407 * p0,p0,k-, k-,p-,p+, p-,kb,p0
408 * 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
413 dimension nopik(6,nmode),npik(nmode)
414 *am outgoing multiplicity and flavors of multi-pion /multi-k modes
423 DATA nopik / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
424 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
425 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
426 3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
427 4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
428 5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
429 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
430 c ajwmod fix sign bug, 2/22/99
431 7 -3,-4, 0, 0, 0, 0 /
432 * list of branching ratios
437 IF(i.EQ. 1) gamprt(i) =0.1800
438 IF(i.EQ. 2) gamprt(i) =0.1751
439 IF(i.EQ. 3) gamprt(i) =0.1110
440 IF(i.EQ. 4) gamprt(i) =0.2515
441 IF(i.EQ. 5) gamprt(i) =0.1790
442 IF(i.EQ. 6) gamprt(i) =0.0071
443 IF(i.EQ. 7) gamprt(i) =0.0134
444 IF(i.EQ. 8) gamprt(i) =0.0450
445 IF(i.EQ. 9) gamprt(i) =0.0100
446 IF(i.EQ.10) gamprt(i) =0.0009
447 IF(i.EQ.11) gamprt(i) =0.0004
448 IF(i.EQ.12) gamprt(i) =0.0003
449 IF(i.EQ.13) gamprt(i) =0.0005
450 IF(i.EQ.14) gamprt(i) =0.0015
451 IF(i.EQ.15) gamprt(i) =0.0015
452 IF(i.EQ.16) gamprt(i) =0.0015
453 IF(i.EQ.17) gamprt(i) =0.0005
454 IF(i.EQ.18) gamprt(i) =0.0050
455 IF(i.EQ.19) gamprt(i) =0.0055
456 IF(i.EQ.20) gamprt(i) =0.0017
457 IF(i.EQ.21) gamprt(i) =0.0013
458 IF(i.EQ.22) gamprt(i) =0.0010
459 IF(i.EQ. 1) oldnames(i)=
' TAU- --> E- '
460 IF(i.EQ. 2) oldnames(i)=
' TAU- --> MU- '
461 IF(i.EQ. 3) oldnames(i)=
' TAU- --> PI- '
462 IF(i.EQ. 4) oldnames(i)=
' TAU- --> PI-, PI0 '
463 IF(i.EQ. 5) oldnames(i)=
' TAU- --> A1- (two subch) '
464 IF(i.EQ. 6) oldnames(i)=
' TAU- --> K- '
465 IF(i.EQ. 7) oldnames(i)=
' TAU- --> K*- (two subch) '
466 IF(i.EQ. 8) names(i-7)=
' TAU- --> 2PI-, PI0, PI+ '
467 IF(i.EQ. 9) names(i-7)=
' TAU- --> 3PI0, PI- '
468 IF(i.EQ.10) names(i-7)=
' TAU- --> 2PI-, PI+, 2PI0 '
469 IF(i.EQ.11) names(i-7)=
' TAU- --> 3PI-, 2PI+, '
470 IF(i.EQ.12) names(i-7)=
' TAU- --> 3PI-, 2PI+, PI0 '
471 IF(i.EQ.13) names(i-7)=
' TAU- --> 2PI-, PI+, 3PI0 '
472 IF(i.EQ.14) names(i-7)=
' TAU- --> K-, PI-, K+ '
473 IF(i.EQ.15) names(i-7)=
' TAU- --> K0, PI-, K0B '
474 IF(i.EQ.16) names(i-7)=
' TAU- --> K-, K0, PI0 '
475 IF(i.EQ.17) names(i-7)=
' TAU- --> PI0 PI0 K- '
476 IF(i.EQ.18) names(i-7)=
' TAU- --> K- PI- PI+ '
477 IF(i.EQ.19) names(i-7)=
' TAU- --> PI- K0B PI0 '
478 IF(i.EQ.20) names(i-7)=
' TAU- --> ETA PI- PI0 '
479 IF(i.EQ.21) names(i-7)=
' TAU- --> PI- PI0 GAM '
480 IF(i.EQ.22) names(i-7)=
' TAU- --> K- K0 '
489 idffin(j,i)=nopik(j,i)
494 * --- coefficients to fix ratio of:
495 * --- a1 3charged/ a1 1charged 2 neutrals matrix elements(masless lim.)
496 * --- probability of k0 to be ks
497 * --- probability of k0b to be ks
498 * --- ratio of coefficients for k*--> k0 pi-
499 * --- all coefficents should be in the range(0.0,1.0)
500 * --- they meaning is probability of the first choice only
IF one
501 * --- neglects mass-phase space effects
515 * zw 13.04.89 here was an error
516 scabib = sqrt(1.-ccabib**2)
518 gamel = gfermi**2*amtau**5/(192*pi**3)
520 * CALL dexay(-1,pol1)
524 FUNCTION dcdmas(IDENT)
525 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
526 * ,ampiz,ampi,amro,gamro,ama1,gama1
527 * ,amk,amkz,amkst,gamkst
529 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
530 * ,ampiz,ampi,amro,gamro,ama1,gama1
531 * ,amk,amkz,amkst,gamkst
532 IF (ident.EQ. 1)
THEN
534 ELSEIF (ident.EQ.-1)
THEN
536 ELSEIF (ident.EQ. 2)
THEN
538 ELSEIF (ident.EQ.-2)
THEN
540 ELSEIF (ident.EQ. 3)
THEN
542 ELSEIF (ident.EQ.-3)
THEN
544 ELSEIF (ident.EQ. 4)
THEN
546 ELSEIF (ident.EQ.-4)
THEN
548 ELSEIF (ident.EQ. 8)
THEN
550 ELSEIF (ident.EQ.-8)
THEN
552 ELSEIF (ident.EQ. 9)
THEN
554 ELSEIF (ident.EQ.-9)
THEN
557 print *,
'STOP IN APKMAS, WRONG IDENT=',ident
562 FUNCTION lunpik(ID,ISGN)
563 COMMON / taukle / bra1,brk0,brk0b,brks
564 REAL*4 bra1,brk0,brk0b,brks
567 IF (ident.EQ. 1)
THEN
569 ELSEIF (ident.EQ.-1)
THEN
571 ELSEIF (ident.EQ. 2)
THEN
573 ELSEIF (ident.EQ.-2)
THEN
575 ELSEIF (ident.EQ. 3)
THEN
577 ELSEIF (ident.EQ.-3)
THEN
579 ELSEIF (ident.EQ. 4)
THEN
581 * k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
583 IF (xio(1).GT.brk0)
THEN
588 ELSEIF (ident.EQ.-4)
THEN
590 * k0b--> k0_long(is 130) / k0_short(is 310) = 1/1
592 IF (xio(1).GT.brk0b)
THEN
597 ELSEIF (ident.EQ. 8)
THEN
599 ELSEIF (ident.EQ.-8)
THEN
601 ELSEIF (ident.EQ. 9)
THEN
603 ELSEIF (ident.EQ.-9)
THEN
606 print *,
'STOP IN IPKDEF, WRONG IDENT=',ident
614 SUBROUTINE taurdf(KTO)
615 c this routine can be called before any tau+ or tau- event is generated
616 c it can be used to generate tau+ and tau- samples of different
618 COMMON / taukle / bra1,brk0,brk0b,brks
619 REAL*4 bra1,brk0,brk0b,brks
620 COMMON / taubra / gamprt(30),jlist(30),nchan
623 c ajwmod: set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
630 c ajwmod: set the brs for(a1+ -> rho+ pi0) and(k*+ -> k0 pi+)
639 SUBROUTINE iniphy(XK00)
640 * ----------------------------------------------------------------------
641 * initialisation of parameters
642 * used in qed and/or gsw routines
643 * ----------------------------------------------------------------------
644 COMMON / qedprm /alfinv,alfpi,xk0
645 REAL*8 alfinv,alfpi,xk0
648 pi8 = 4.d0*datan(1.d0)
650 alfpi = 1d0/(alfinv*pi8)
655 c ----------------------------------------------------------------------
656 c initialisation of masses
659 c ----------------------------------------------------------------------
660 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
661 * ,ampiz,ampi,amro,gamro,ama1,gama1
662 * ,amk,amkz,amkst,gamkst
664 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
665 * ,ampiz,ampi,amro,gamro,ama1,gama1
666 * ,amk,amkz,amkst,gamkst
668 c in-coming / out-going fermion masses
670 c --- let us update tau mass ...
678 * masses used in tau decays
692 c in-coming / out-going fermion masses
697 c masses used in tau decays cleo settings
714 c subsitute of tau production generator
716 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
717 * ,ampiz,ampi,amro,gamro,ama1,gama1
718 * ,amk,amkz,amkst,gamkst
720 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
721 * ,ampiz,ampi,amro,gamro,ama1,gama1
722 * ,amk,amkz,amkst,gamkst
724 c positions of taus in the lund
common block
725 c it will be used by tauola output routines.
726 COMMON /taupos / npa,npb
727 dimension xpb1(4),xpb2(4),aqf1(4),aqf2(4)
729 c --- defining dummy events momenta
739 CALL tralo4(1,aqf1,aqf1,am)
740 CALL tralo4(2,aqf2,aqf2,am)
741 c --- beams momenta and identifiers
742 kfb1= 11*idff/iabs(idff)
743 kfb2=-11*idff/iabs(idff)
747 $ xpb1(3)= aqf1(4)*aqf1(3)/abs(aqf1(3))
751 $ xpb2(3)= aqf2(4)*aqf2(3)/abs(aqf2(3))
752 c --- position of first and second tau in lund
common
755 c --- fill to lund
COMMON
756 CALL filhep( 1,3, kfb1,0,0,0,0,xpb1, amel,.true.)
757 CALL filhep( 2,3, kfb2,0,0,0,0,xpb2, amel,.true.)
758 CALL filhep(npa,1, idff,1,2,0,0,aqf1,amtau,.true.)
759 CALL filhep(npb,1,-idff,1,2,0,0,aqf2,amtau,.true.)
761 SUBROUTINE tralo4(KTO,P,Q,AM)
762 c **************************
763 c subsitute of tralo4
766 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
767 * ,ampiz,ampi,amro,gamro,ama1,gama1
768 * ,amk,amkz,amkst,gamkst
770 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
771 * ,ampiz,ampi,amro,gamro,ama1,gama1
772 * ,amk,amkz,amkst,gamkst
775 etau=sqrt(ptau**2+amtau**2)
776 exe=(etau+ptau)/amtau
777 IF(kto.EQ.2) exe=(etau-ptau)/amtau
779 c ======================================================================
780 c
END of the test job
781 c ======================================================================
783 SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
784 c ----------------------------------------------------------------------
785 c this
subroutine fills one entry into the hepevt common
786 c and updates the information for affected mother entries
788 c written by martin w. gruenewald(91/01/28)
790 c called by : ztohep,btohep,dwluxy
791 c ----------------------------------------------------------------------
793 c this is the hepevt class in old style. no d_h_ class pre-name
795 parameter(nmxhep=10000)
797 INTEGER nevhep,nhep,isthep,idhep,jmohep,
808 * ----------------------------------------------------------------------
812 * ----------------------------------------------------------------------
814 c parameter(nmxhep=2000)
815 c common/hepevt/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
816 c &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
818 c common/phoqed/qedrad(nmxhep)
830 ELSE IF (n.GT.0)
THEN
841 IF ((ihep.LE.0).OR.(ihep.GT.nmxhep))
RETURN
848 IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
850 IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
857 c koral-b and koral-z
do not provide vertex and/or lifetime informations
865 DO ip=jmohep(1,ihep),jmohep(2,ihep)
868 c
if there is a daughter at ihep, mother entry at ip has decayed
869 IF(isthep(ip).EQ.1)isthep(ip)=2
871 c and daughter pointers of mother entry must be updated
872 IF(jdahep(1,ip).EQ.0)
THEN
876 jdahep(2,ip)=max(ihep,jdahep(2,ip))