1/* copyright(c) 1991-2022 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 <https://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/* glibc
's intent is to support the IEC 559 math functionality, real
28 and complex. If the GCC (4.9 and later) predefined macros
29 specifying compiler intent are available, use them to determine
30 whether the overall intent is to support these features; otherwise,
31 presume an older compiler has intent to support these features and
32 define these macros by default. */
36/* wchar_t uses Unicode 10.0.0. Version 10.0 of the Unicode Standard is
37 synchronized with ISO/IEC 10646:2017, fifth edition, plus
38 the following additions from Amendment 1 to the fifth edition:
41 - 3 additional Zanabazar Square characters */
43 SUBROUTINE TAUOLA(MODE,KEYPOL)
44C *************************************
45C general tauola interface, should work in every case until
46C hepevt is OK, does not check if hepevt is 'clean
'
47C in particular will decay decayed taus...
48C only longitudinal spin effects are included.
49C in W decay v-a vertex is assumed
50C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
51C this is the hepevt class in old style. No d_h_ class pre-name
53 PARAMETER (NMXHEP=10000)
54 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
55 INTEGER nevhep,nhep,isthep,idhep,jmohep,
58 $ nevhep, ! serial number
59 $ nhep, ! number of particles
60 $ isthep(nmxhep), ! status code
61 $ idhep(nmxhep), ! particle ident KF
62 $ jmohep(2,nmxhep), ! parent particles
63 $ jdahep(2,nmxhep), ! childreen particles
64 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
65 $ vhep(4,nmxhep) ! vertex [mm]
66* ----------------------------------------------------------------------
69 $ qedrad(nmxhep) ! Photos flag
70* ----------------------------------------------------------------------
72 COMMON /TAUPOS/ NP1, NP2
73 REAL*4 PHOI(4),PHOF(4)
74 double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
75 COMMON / MOMDEC / Q1,Q2,P1,P2,P3,P4
76* tauola, photos and jetset overall switches
77 COMMON /LIBRA/ JAK1,JAK2,ITDKRC,IFPHOT,IFHADM,IFHADP
80 common /pseudocoup/ csc,ssc
83 COMMON / INOUT / INUT,IOUT
85C to switch tau polarization OFF in taus
86 DIMENSION POL1(4), POL2(4)
87 double precision POL1x(4), POL2x(4)
89 DATA POL1 /0.0,0.0,0.0,0.0/
90 DATA POL2 /0.0,0.0,0.0,0.0/
91 DATA PI /3.141592653589793238462643D0/
94 DIMENSION IMOTHER (20)
97C store daughter pointers
99 COMMON /ISONS_TAU/ISON(2)
103C ***********************
105 JAK1 = 0 ! decay mode first tau
106 JAK2 = 0 ! decay mode second tau
107 ITDKRC=1.0 ! switch of radiative corrections in decay
108 IFPHOT=1.0 ! PHOTOS switch
111 POL=1.0 ! tau polarization dipswitch must be 1 or 0
121C couplings of the 'pseudoscalar higgs
' as in CERN-TH/2003-166
123 xmtau=1.777 ! tau mass
124 xmh=120 ! higgs boson mass
125 betah=sqrt(1d0-4*xmtau**2/xmh**2)
128C write(*,*) ' scalar component=
',csc,' pseudo-scalar component=
',ssc
129.EQ.
IF (IFPHOT1) CALL PHOINI ! this if PHOTOS was not initialized earlier
130 CALL INIETC(JAK1,JAK2,ITDKRC,IFPHOT)
134C activation of pi0 and eta decays: (1) means on, (0) off
138 CALL TAUPI0(-1,1,ION)
140 WRITE(IOUT,7001) pol,psi,ION(1),ION(2),ION(3)
141.EQ.
ELSEIF(MODE0) THEN
142C ***********************
144C..... find tau-s and fill common block /TAUPOS/
145C this is to avoid LUND history fillings. This call is optional
146 CALL PHYFIX(NSTOP,NSTART)
147C clear mothers of the previous event
155C ... and to find mothers giving taus.
157C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0)
159.EQ..AND..EQ..AND.
IF(ABS(IDHEP(I))KFTAUISTHEP(I)1
160.GE..OR..LT.
$ (ISTHEP(I)125ISTHEP(I)120)) THEN
162.EQ.
DO WHILE (ABS(IDHEP(IMOTH))KFTAU) ! KEEP WALKING UP
163 IMOTH=JMOHEP(1,IMOTH)
165.EQ..OR.
IF (ISTHEP(IMOTH)3
166.GE..AND..LE.
$ (ISTHEP(IMOTH)120ISTHEP(IMOTH)125)) THEN
167 DO J=NSTART,NHEP ! WE HAVE WALKED INTO HARD RECORD
168.EQ..AND.
IF (IDHEP(J)IDHEP(IMOTH)
169.EQ..AND.
$ JMOHEP(1,J)IMOTH
170.EQ.
$ ISTHEP(J)2) THEN
180.EQ.
IF(JMOTHIMOTHER(II)) GOTO 9999
189C ... taus of every mother are treated in this main loop
198C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT..
200.EQ.
IF (IDHEP(JMOHEP(1,IM0))IDHEP(IM0)) IM0=JMOHEP(1,IM0)
203.EQ..OR.
IF (ISTHEP(I)3
204.GE..AND..LE.
$ (ISTHEP(I)120ISTHEP(I)125)) THEN ! HARD RECORD
208.EQ..OR..EQ.
DO WHILE (IDHEP(IMOTH)IDHEP(I)ABS(IDHEP(IMOTH))KFTAU) ! KEEP WALKING UP
209 IMOTH=JMOHEP(1,IMOTH)
211.EQ..OR..EQ..AND..EQ.
IF ((IMOTHIM0IMOTHIM)ISEL-1) THEN
214.EQ..OR..EQ..AND..EQ.
ELSEIF ((IMOTHIM0IMOTHIM)ISEL0) THEN
216.NE..AND..NE..AND..EQ.
ELSEIF ((IMOTHIM0IMOTHIM)ISEL0) THEN
226C ... we correct HEPEVT (fix developped with Catherine BISCARAT)
227.EQ.
c IF (JDAHEP(2,IM)0) THEN ! ID of second daughter was missing
229c DO I=JDAHEP(1,IM)+1,NHEP ! OK lets look for it
230.EQ..AND..EQ.
c IF (JMOHEP(1,I)IMISECU1) THEN ! we have found one
232.EQ..AND..NE.
c ELSEIF (JMOHEP(1,I)IMISECU1) THEN ! we have found one after there
233c JDAHEP(2,IM)=0 ! was something else, lets kill game
235.NE.
c IF (JMOHEP(1,I)IM) ISECU=0 ! other stuff starts
239C ... we check whether there are just two or more tau-likes
241.EQ..OR..EQ.
IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NCOUNT=NCOUNT+1
242.EQ..OR..EQ.
IF(IDHEP(I) KFTAUIDHEP(I) KFNUE) NCOUNT=NCOUNT+1
245C ... if there will be more we will come here again
249 DO I=MAX(NP1+1,ISON(1)),ISON(2)
251.EQ..OR..EQ.
IF(IDHEP(I)-KFTAUIDHEP(I)-KFNUE) NP1=I
254 DO I=MAX(NP2+1,ISON(1)),ISON(2)
256.EQ..OR..EQ.
IF(IDHEP(I) KFTAUIDHEP(I) KFNUE) NP2=I
259 P1(I)= PHEP(I,NP1) !momentum of tau+
260 P2(I)= PHEP(I,NP2) !momentum of tau-
268c.....include polarisation effect
271.EQ..OR..EQ..OR.
IF(IDHEP(IM)KFHIGGS(1)IDHEP(IM)KFHIGGS(2)
272.EQ.
$ IDHEP(IM)KFHIGGS(3)) THEN ! case of Higgs
273.LT.
IF(RRR(1)0.5) THEN
280.EQ..OR..EQ.
ELSEIF((IDHEP(IM)KFZ0)(IDHEP(IM)KFGAM)) THEN ! case of gamma/Z
281C there is no angular dependence in gamma/Z polarization
282C there is no s-dependence in gamma/Z polarization at all
283C there is even no Z polarization in any form
284C main reason is that nobody asked ...
285C but it is prepared and longitudinal correlations
286C can be included up to KORALZ standards
288 POLZ0=PLZAPX(.true.,IM,NP1,NP2)
289.LT.
IF(RRR(1)POLZ0) THEN
296.EQ.
ELSEIF(IDHEP(NP1)-IDHEP(NP2))THEN ! undef orig. only s-dep poss.
297 POLZ0=PLZAPX(.true.,IM,NP1,NP2)
298.LT.
IF(RRR(1)POLZ0) THEN
305.EQ.
ELSEIF(ABS(IDHEP(IM))KFHIGCH) THEN ! case of charged Higgs
308 ELSE ! case of W+ or W-
312c.....include polarisation effect
315.EQ..OR..EQ..OR.
IF(IDHEP(IM)KFHIGGS(1)IDHEP(IM)KFHIGGS(2)
316.EQ.
$ IDHEP(IM)KFHIGGS(3)) THEN
317.EQ..AND.
IF(IDHEP(NP1)-KFTAU
318.LE..OR..GT..AND.
$ (JDAHEP(1,NP1)NP1JDAHEP(1,NP1)NHEP)
319.EQ..AND.
$ IDHEP(NP2) KFTAU
320.LE..OR..GT.
$ (JDAHEP(1,NP2)NP2JDAHEP(1,NP2)NHEP)
322.EQ.
IF (IDHEP(IM)KFHIGGS(1)) THEN
324.EQ.
ELSEIF (IDHEP(IM)KFHIGGS(2)) THEN
326.EQ.
ELSEIF (IDHEP(IM)KFHIGGS(3)) THEN
329 WRITE(*,*) 'warning from tauola:
'
330 WRITE(*,*) 'i stop this run, wrong idhep(im)=
',IDHEP(IM)
333 CALL SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
334.EQ.
IF (IFPHOT1) CALL PHOTOS(IM) ! Bremsstrahlung in Higgs decay
335 ! AFTER adding taus !!
340.EQ..AND.
IF(IDHEP(NP1)-KFTAU
341.LE..OR..GT.
$ (JDAHEP(1,NP1)NP1JDAHEP(1,NP1)NHEP)) THEN
342C here check on if NP1 was not decayed should be verified
344.EQ.
IF (IFPHOT1) CALL PHOTOS(NP1)
348.EQ..AND.
IF(IDHEP(NP2) KFTAU
349.LE..OR..GT.
$ (JDAHEP(1,NP2)NP2JDAHEP(1,NP2)NHEP)) THEN
350C here check on if NP2 was not decayed should be added
352.EQ.
IF (IFPHOT1) CALL PHOTOS(NP2)
357.GT.
IF (NCOUNT0) GOTO 666
360.EQ.
ELSEIF(MODE1) THEN
361C ***********************
364 CALL DEKAY(100,POL1x)
368 7001 FORMAT(///1X,15(5H*****)
369 $ /,' *
', 25X,'*****tauola universal interface: ******
',9X,1H*,
370 $ /,' *
', 25X,'*****version 1.22, april 2009 (gfort)**
',9X,1H*,
371 $ /,' *
', 25X,'**authors: p. golonka, b. kersevan, ***
',9X,1H*,
372 $ /,' *
', 25X,'**t. pierzchala, e. richter-was, ******
',9X,1H*,
373 $ /,' *
', 25X,'****** z. was, m. worek ***************
',9X,1H*,
374 $ /,' *
', 25X,'**useful discussions, in particular ***
',9X,1H*,
375 $ /,' *
', 25X,'*with c. biscarat and s. slabospitzky**
',9X,1H*,
376 $ /,' *
', 25X,'****** are warmly acknowledged ********
',9X,1H*,
377 $ /,' *
', 25X,' ',9X,1H*,
378 $ /,' *
', 25X,'********** initialization ************
',9X,1H*,
379 $ /,' *
',F20.5,5X,'tau polarization switch must be 1 or 0
',9X,1H*,
380 $ /,' *
',F20.5,5X,'higs scalar/pseudo mix cern-th/2003-166
',9X,1H*,
381 $ /,' *
',I10, 15X,'pi0 decay switch must be 1 or 0
',9X,1H*,
382 $ /,' *
',I10, 15X,'eta decay switch must be 1 or 0
',9X,1H*,
383 $ /,' *
',I10, 15X,'k0s decay switch must be 1 or 0
',9X,1H*,
386 7002 FORMAT(///1X,15(5H*****)
387 $ /,' *
', 25X,'*****tauola universal interface: ******
',9X,1H*,
388 $ /,' *
', 25X,'*****version 1.22, april 2009 (gfort)**
',9X,1H*,
389 $ /,' *
', 25X,'**authors: p. golonka, b. kersevan, ***
',9X,1H*,
390 $ /,' *
', 25X,'**t. pierzchala, e. richter-was, ******
',9X,1H*,
391 $ /,' *
', 25X,'****** z. was, m. worek ***************
',9X,1H*,
392 $ /,' *
', 25X,'**useful discussions, in particular ***
',9X,1H*,
393 $ /,' *
', 25X,'*with c. biscarat and s. slabospitzky**
',9X,1H*,
394 $ /,' *
', 25X,'****** are warmly acknowledged ********
',9X,1H*,
395 $ /,' *
', 25X,'******
END OF MODULE OPERATION ********
',9X,1H*,
400 SUBROUTINE SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
402 REAL*8 HH1,HH2,wthiggs
403 DIMENSION POL1(4), POL2(4),HH1(4),HH2(4), RRR(1)
404! CALL DEXAY(1,POL1) ! Kept for tests
405! CALL DEXAY(2,POL2) ! Kept for tests
411 wt=wthiggs(IFPSEUDO,HH1,HH2)
412.GT.
IF (RRR(1)WT) GOTO 10
418 FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
420 common /pseudocoup/ csc,ssc
423 REAL*8 HH1(4),HH2(4),R(4,4),wthiggs
431 R(4,4)= 1D0 ! unpolarized part
432 R(3,3)=-1D0 ! longitudinal
437 R(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
438 R(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
439 R(1,2)=2*csc*ssc/(csc**2+ssc**2)
440 R(2,1)=-2*csc*ssc/(csc**2+ssc**2)
450 WTHIGGS=WTHIGGS+R(K,L)*HH1(K)*HH2(L)
456 FUNCTION PLZAPX(HOPEin,IM0,NP1,NP2)
457C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block.
458C the purpose of this routine is to calculate polarization of Z along tau direction.
459C this is highly non-trivial due to necessity of reading infromation from hard process
460C history in HEPEVT, which is often written not up to the gramatic rules.
461 REAL*8 PLZAP0,SVAR,COSTHE,sini,sfin,ZPROP2,
462 $ P1(4),P2(4),Q1(4),Q2(4),QQ(4),PH(4),PD1(4),PD2(4),
463 $ PQ1(4),PQ2(4),PB(4),PA(4)
467C this is the hepevt class in old style. No d_h_ class pre-name
469 PARAMETER (NMXHEP=10000)
470 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
471 INTEGER nevhep,nhep,isthep,idhep,jmohep,
474 $ nevhep, ! serial number
475 $ nhep, ! number of particles
476 $ isthep(nmxhep), ! status code
477 $ idhep(nmxhep), ! particle ident KF
478 $ jmohep(2,nmxhep), ! parent particles
479 $ jdahep(2,nmxhep), ! childreen particles
480 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
481 $ vhep(4,nmxhep) ! vertex [mm]
482* ----------------------------------------------------------------------
485 $ qedrad(nmxhep) ! Photos flag
486* ----------------------------------------------------------------------
489C(BPK)--> BROTHERS OF TAU ALREADY FOUND
491 COMMON /ISONS_TAU/ISON(2)
494C >> STEP 1: find where are particles in hepevent and pick them up
497C sometimes shade Z of Z is its mother ...
500C to protect against check on mother of beam particles.
502.EQ.
IF (IDHEP(IM0)IDHEP(IM00)) IM=JMOHEP(1,IM0)
505C find (host generator-level) incoming beam-bare-particles which form Z and co.
509C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS
511.EQ.
IF (ISTHEP(IM00)120) THEN
518C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first
' in hepevt.
519.EQ..AND..EQ.
IF (IMO10IMO20) THEN
522.EQ.
IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
526.EQ.
IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
529C case when it was like qq --> tau+tau- gammas and qq were NOT 'first
' in hepevt.
531.NE..AND..NE.
ELSEIF (IDHEP(IM)22IDHEP(IM)23) THEN
534.EQ.
IF (IDHEP(IMO1)IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
538.EQ.
IF (IDHEP(IMO2)IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
544C and check if it really happened
545.EQ.
IF (IMO10) HOPE=.FALSE.
546.EQ.
IF (IMO20) HOPE=.FALSE.
547.EQ.
IF (IMO1IMO2) HOPE=.FALSE.
551 Q1(I)= PHEP(I,NP1) !momentum of tau+
552 Q2(I)= PHEP(I,NP2) !momentum of tau-
555C corrections due to possible differences in 4-momentum of shadow vs true Z.
557.EQ..AND.
IF (IMJMOHEP(1,IM0)
558.EQ..OR..EQ.
$ (IDHEP(IM)22IDHEP(IM)23)) THEN
565 CALL BOSTDQ( 1,PA, Q1, Q1)
566 CALL BOSTDQ( 1,PA, Q2, Q2)
567 CALL BOSTDQ(-1,PB, Q1, Q1)
568 CALL BOSTDQ(-1,PB, Q2, Q2)
573 QQ(I)= Q1(I)+Q2(I) !momentum of Z
574 IF (HOPE) P1(I)=PHEP(I,IMO1) !momentum of beam1
575 IF (HOPE) P2(I)=PHEP(I,IMO2) !momentum of beam2
580! These momenta correspond to quarks, gluons photons or taus
583 IF (HOPE) IDFP1=IDHEP(IMO1)
584 IF (HOPE) IDFP2=IDHEP(IMO2)
586 SVAR=QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2
588C options which may be useful in some cases of two heavy boson production
589C need individual considerations. To be developed.
590C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
591C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
592 PLZAPX=0.5 ! pure gamma
596C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
598C let us define beginning and end of particles which are produced in parallel to Z
599C in principle following should work
601C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
605 INBR=IM ! OK, HARD RECORD Z/GAMMA
606.EQ.
IF (IFFULL1) INBR=NP1 ! OK, NO Z/GAMMA
607.EQ.
IF (IDHEP(JMOHEP(1,INBR))IDHEP(INBR)) INBR=JMOHEP(1,INBR) ! FORCE HARD RECORD
609.EQ..OR..EQ.
IF(NX10NX20) THEN
613.EQ.
IF(JMOHEP(1,INBR-K)JMOHEP(1,INBR)) THEN
622.EQ.
IF(JMOHEP(1,K)JMOHEP(1,INBR)) THEN
631C case of annihilation of two bosons is hopeles
632.GE..AND..GE.
IF (ABS(IDFP1)20ABS(IDFP2)20) HOPE=.FALSE.
633C case of annihilation of non-matching flavors is hopeless
634.LE..AND..LE..AND..NE.
IF (ABS(IDFP1)20ABS(IDFP2)20IDFP1+IDFP20)
637C options which may be useful in some cases of two heavy boson production
638C need individual considerations. To be developed.
639C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam
640C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z
641 PLZAPX=0.5 ! pure gamma
644.LT.
IF (ABS(IDFP1)20) IDE= IDFP1
645.LT.
IF (ABS(IDFP2)20) IDE=-IDFP2
649C >> STEP 3 we combine gluons, photons into incoming effective beams
652C in the following we will ignore the possibility of photon emission from taus
653C however at certain moment it will be necessary to take care of
665 IFLAV=min(ABS(IDFP1),ABS(IDFP2))
667*--------------------------------------------------------------------------
668c IFLAV=min(ABS(IDFP1),ABS(IDFP2))
669c that means that always quark or lepton i.e. process like
670c f g(gamma) --> f Z0 --> tau tau
671c we glue fermions to effective beams that is f f --> Z0 --> tau tau
672c with gamma/g emission from initial fermion.
673*---------------------------------------------------------------------------
675.GE.
IF (ABS(IDFP1)20) THEN
678.EQ.
IF (ABS(IDP)IFLAV) THEN
686.GE.
IF (ABS(IDFP2)20) THEN
689.EQ.
IF (ABS(IDP)IFLAV) THEN
697C if first beam was boson: gluing
699.GE.
IF (ABS(IDFP1)20) THEN
703 xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
704 $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
705 xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
706 $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
718C if second beam was boson: gluing
721.GE.
IF (ABS(IDFP2)20) THEN
725 xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
726 $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
727 xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
728 $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
745.EQ.
IF (IDHEP(JMOHEP(1,NP1))IDHEP(NP1)) NPH1=JMOHEP(1,NP1) ! SHOULD PUT US IN HARD REC.
746.EQ.
IF (IDHEP(JMOHEP(1,NP2))IDHEP(NP2)) NPH2=JMOHEP(1,NP2) ! SHOULD PUT US IN HARD REC.
750.NE..AND..NE..AND.
IF (ABS(IDHEP(K))IFLAVKIM
752.NE..AND..NE.
$ KNPH1KNPH2) THEN
754.EQ..AND..EQ.
IF(IDHEP(K)22IFFULL1) THEN
758 xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
759 $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
760 xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
761 $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
762 xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
763 $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
764 xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
765 $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)
768 sini=abs((PD1(4)+PD2(4)-PH(4))**2-(PD1(3)+PD2(3)-PH(3))**2
769 $ -(PD1(2)+PD2(2)-PH(2))**2-(PD1(1)+PD2(1)-PH(1))**2)
770 sfin=abs((PD1(4)+PD2(4) )**2-(PD1(3)+PD2(3) )**2
771 $ -(PD1(2)+PD2(2) )**2-(PD1(1)+PD2(1) )**2)
781 XM=MIN(XM1,XM2,XM3,XM4)
786.EQ.
ELSEIF (XM2XM) THEN
790.EQ.
ELSEIF (XM3XM) THEN
803 xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
804 $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
805 xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
806 $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
822C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in
823c >> effective outcoming taus
825C let us define beginning and end of particles which are produced in
830C find outcoming particles which come from Z
835C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER
836.EQ..OR..EQ.
IF (ABS(IDHEP(IM0))22abs(IDHEP(IM0))23) THEN
838.EQ.
IF(ABS(IDHEP(K))22) THEN
845 xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
846 $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
847 xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
848 $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)
867*------------------------------------------------------------------------
870C out of effective momenta we calculate COSTHE and later polarization
871 CALL ANGULU(PD1,PD2,Q1,Q2,COSTHE)
873 PLZAPX=PLZAP0(IDE,IDFQ1,SVAR,COSTHE)
876 SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE)
877 REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
878C take effective beam which is less massive, it should be irrelevant
879C but in case HEPEVT is particulary dirty may help.
880C this routine calculate reduced system transver and cosine of scattering
883 XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2)
884 XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2)
896C calculate space like part of P (in Z restframe)
902 XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2)
904 QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1)
906 QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2
909 PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1)
911 P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2
914 PXP =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
915 QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2)
916 PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4)
917 COSTHE=PXQT/PXP/QTXQT
921 FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0)
922C this function calculates probability for the helicity +1 +1 configuration
923C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
924 REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN
927.LT.
C >>>>> IF (IDE*IDF0) COSTHE=-COSTH0 ! this is probably not needed ID
928C >>>>> of first beam is used by T_GIVIZ0 including sign
931 CALL INITWK(IDE,IDF,SVAR)
933 CALL INITWK(-IDE,-IDF,SVAR)
935 PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0)
936 $ /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0))
940 FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB)
941C ----------------------------------------------------------------------
942C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME
943C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER
944C EXAMPLE OF THE METHOD APPLIED THERE
945C INPUT PARAMETERS ARE: SVAR -- transfer
946C COSTHE -- cosine of angle between tau+ and 1st beam
947C TA,TB -- helicity states of tau+ tau-
949C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
950C ----------------------------------------------------------------------
951 IMPLICIT REAL*8(A-H,O-Z)
952 COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
953 REAL*8 ENE ,AMIN,AMFIN
954 COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
955 & ,XUPGI ,XUPZI ,XUPGF ,XUPZF
956 & ,NDIAG0,NDIAGA,KEYA,KEYZ
957 & ,ITCE,JTCE,ITCF,JTCF,KOLOR
958 REAL*8 SS,POLN,T3E,QE,T3F,QF
959 & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
961C=====================================================================
962 COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
963 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
964C SWSQ = sin2 (theta Weinberg)
965C AMW,AMZ = W & Z boson masses respectively
966C AMH = the Higgs mass
967C AMTOP = the top mass
969 COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
970 COMPLEX*16 XUPZFP(2),XUPZIP(2)
971 COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
972 COMPLEX*16 PROPA,PROPZ
974 COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
975 COMPLEX*16 XTHING,XVE,XVF,XVEF
976 DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/
979 DATA SVAR0,COST0 /-5.D0,-6.D0/
980 DATA PI /3.141592653589793238462643D0/
981 DATA SEPS1,SEPS2 /0D0,0D0/
983C MEMORIZATION =========================================================
984.NE..OR..NE..OR..NE.
IF ( MODEMODE0SVARSVAR0COSTHECOST0
985.OR..NE.
$ IDE0IDE)THEN
993 SINTHE=SQRT(1.D0-COSTHE**2)
994 BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR))
995C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
996 XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2))
997 XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2))
998 XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2))
999 XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2))
1000C FINAL STATE VECTOR COUPLING
1001 XUPF =0.5D0*(XUPZF(1)+XUPZF(2))
1002 XUPI =0.5D0*(XUPZI(1)+XUPZI(2))
1006 PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ)
1007.EQ.
IF (KEYGSW0) PROPZ=0.D0
1010 REGULA= (3-2*I)*(3-2*J) + COSTHE
1011 REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR)
1012 APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA)
1013 AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA
1014 ABORN(I,J)=APHOT(I,J)+AZETT(I,J)
1015 APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM
1016 AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM
1017 ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J)
1022C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
1023C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
1024C* HELICITY CONSERVATION EXPLICITLY OBEYED
1032 FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0
1033 FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB)
1034 FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB)
1036 BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR
1039 BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
1045.LT.
IF(FUNT0.D0) FUNT=BORN
1048.GT.
IF (SVAR4D0*AMFIN**2) THEN
1049C PHASE SPACE THRESHOLD FACTOR
1050 THRESH=SQRT(1-4D0*AMFIN**2/SVAR)
1051 T_BORN= FUNT*SVAR**2*THRESH
1056C ZW HERE WAS AN ERROR 19. 05. 1989
1057! write(*,*) 'KKKK
',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
1058! write(*,*) 'KKKK X
',svar,costhe,TA,TB,T_BORN
1061 SUBROUTINE INITWK(IDEX,IDFX,SVAR)
1062! initialization routine coupling masses etc.
1063 IMPLICIT REAL*8 (A-H,O-Z)
1064 COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
1065 REAL*8 ENE ,AMIN,AMFIN
1066 COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
1067 & ,XUPGI ,XUPZI ,XUPGF ,XUPZF
1068 & ,NDIAG0,NDIAGA,KEYA,KEYZ
1069 & ,ITCE,JTCE,ITCF,JTCF,KOLOR
1070 REAL*8 SS,POLN,T3E,QE,T3F,QF
1071 & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
1072 COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
1073 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
1074C SWSQ = sin2 (theta Weinberg)
1075C AMW,AMZ = W & Z boson masses respectively
1076C AMH = the Higgs mass
1077C AMTOP = the top mass
1085.EQ.
IF (IDFX 15) then
1086 IDF=2 ! denotes tau +2 tau-
1087 AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
1088.EQ.
ELSEIF (IDFX-15) then
1089 IDF=-2 ! denotes tau -2 tau-
1090 AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
1092 WRITE(*,*) 'INITWK: WRONG IDFX
'
1096.EQ.
IF (IDEX 11) then !electron
1099.EQ.
ELSEIF (IDEX-11) then !positron
1102.EQ.
ELSEIF (IDEX 13) then !mu+
1105.EQ.
ELSEIF (IDEX-13) then !mu-
1108.EQ.
ELSEIF (IDEX 1) then !d
1111.EQ.
ELSEIF (IDEX- 1) then !d~
1114.EQ.
ELSEIF (IDEX 2) then !u
1117.EQ.
ELSEIF (IDEX- 2) then !u~
1120.EQ.
ELSEIF (IDEX 3) then !s
1123.EQ.
ELSEIF (IDEX- 3) then !s~
1126.EQ.
ELSEIF (IDEX 4) then !c
1129.EQ.
ELSEIF (IDEX- 4) then !c~
1132.EQ.
ELSEIF (IDEX 5) then !b
1135.EQ.
ELSEIF (IDEX- 5) then !b~
1138.EQ.
ELSEIF (IDEX 12) then !nu_e
1141.EQ.
ELSEIF (IDEX- 12) then !nu_e~
1144.EQ.
ELSEIF (IDEX 14) then !nu_mu
1147.EQ.
ELSEIF (IDEX- 14) then !nu_mu~
1150.EQ.
ELSEIF (IDEX 16) then !nu_tau
1153.EQ.
ELSEIF (IDEX- 16) then !nu_tau~
1158 WRITE(*,*) 'INITWK: WRONG IDEX
'
1162C ----------------------------------------------------------------------
1164C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
1167C ----------------------------------------------------------------------
1172 CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM)
1173 CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM)
1177 XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1178 XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1179 CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR)
1180 CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR)
1184 XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1185 XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1196 SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1197C ----------------------------------------------------------------------
1198C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
1199C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
1200C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
1201C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
1202C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
1203C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
1204C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
1206C called by : EVENTE, EVENTM, FUNTIH, .....
1207C ----------------------------------------------------------------------
1208 IMPLICIT REAL*8(A-H,O-Z)
1210.EQ..OR..GT.
IF(IDFERM0IABS(IDFERM)4) GOTO 901
1211.NE.
IF(IABS(IHELIC)1) GOTO 901
1213 IDTYPE =IABS(IDFERM)
1215 LEPQUA=INT(IDTYPE*0.4999999D0)
1216 IUPDOW=IDTYPE-2*LEPQUA-1
1217 CHARGE =(-IUPDOW+2D0/3D0*LEPQUA)*IC
1218 SIZO3 =0.25D0*(IC-IH)*(1-2*IUPDOW)
1220C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
1221C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
1223 901 PRINT *,' STOP IN GIVIZO: WRONG PARAMS.
'
1226 SUBROUTINE PHYFIX(NSTOP,NSTART)
1227 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
1229C NSTOP NSTART : when PHYTIA history ends and event starts.
1233.NE.
IF(K(I,1)21) THEN
1241 SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1242C ----------------------------------------------------------------------
1243C this subroutine fills one entry into the HEPEVT common
1244C and updates the information for affected mother entries
1246C written by Martin W. Gruenewald (91/01/28)
1248C called by : ZTOHEP,BTOHEP,DWLUxy
1249C ----------------------------------------------------------------------
1251C this is the hepevt class in old style. No d_h_ class pre-name
1252C this is the hepevt class in old style. No d_h_ class pre-name
1254 PARAMETER (NMXHEP=10000)
1255 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1256 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1259 $ nevhep, ! serial number
1260 $ nhep, ! number of particles
1261 $ isthep(nmxhep), ! status code
1262 $ idhep(nmxhep), ! particle ident KF
1263 $ jmohep(2,nmxhep), ! parent particles
1264 $ jdahep(2,nmxhep), ! childreen particles
1265 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1266 $ vhep(4,nmxhep) ! vertex [mm]
1267* ----------------------------------------------------------------------
1270 $ qedrad(nmxhep) ! Photos flag
1271* ----------------------------------------------------------------------
1282.GT.
ELSE IF (N0) THEN
1293.LE..OR..GT.
IF ((IHEP0)(IHEPNMXHEP)) RETURN
1300.LT.
IF(JMO10)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP
1302.LT.
IF(JMO20)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP
1309C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
1317 DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
1320C if there is a daughter at IHEP, mother entry at IP has decayed
1321.EQ.
IF(ISTHEP(IP)1)ISTHEP(IP)=2
1323C and daughter pointers of mother entry must be updated
1324.EQ.
IF(JDAHEP(1,IP)0)THEN
1328 JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP))
1337 FUNCTION IHEPDIM(DUM)
1338C this is the hepevt class in old style. No d_h_ class pre-name
1339C this is the hepevt class in old style. No d_h_ class pre-name
1341 PARAMETER (NMXHEP=10000)
1342 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1343 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1346 $ nevhep, ! serial number
1347 $ nhep, ! number of particles
1348 $ isthep(nmxhep), ! status code
1349 $ idhep(nmxhep), ! particle ident KF
1350 $ jmohep(2,nmxhep), ! parent particles
1351 $ jdahep(2,nmxhep), ! childreen particles
1352 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1353 $ vhep(4,nmxhep) ! vertex [mm]
1354* ----------------------------------------------------------------------
1357 $ qedrad(nmxhep) ! Photos flag
1358* ----------------------------------------------------------------------
1363 IMPLICIT REAL*8(A-H,O-Z)
1364 COMPLEX*16 CPRZ0,CPRZ0M
1367 CPRZ0=DCMPLX((S-AMZ**2),S/AMZ*GAMMZ)
1369 ZPROP2=(ABS(CPRZ0M))**2
1372 SUBROUTINE TAUPI0(MODE,JAK,ION)
1373C no initialization required. Must be called once after every:
1374C 1) CALL DEKAY(1+10,...)
1375C 2) CALL DEKAY(2+10,...)
1376C 3) CALL DEXAY(1,...)
1377C 4) CALL DEXAY(2,...)
1378C subroutine to decay originating from TAUOLA's taus:
1379c 1) etas(with
CALL taueta(jak))
1380c 2) later pi0
's from taus.
1381C 3) extensions to other applications possible.
1382C this routine belongs to >tauola universal interface<, but uses
1383C routines from >tauola< utilities as well. 25.08.2005
1384C this is the hepevt class in old style. No d_h_ class pre-name
1386 PARAMETER (NMXHEP=10000)
1387 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1388 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1391 $ nevhep, ! serial number
1392 $ nhep, ! number of particles
1393 $ isthep(nmxhep), ! status code
1394 $ idhep(nmxhep), ! particle ident KF
1395 $ jmohep(2,nmxhep), ! parent particles
1396 $ jdahep(2,nmxhep), ! childreen particles
1397 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1398 $ vhep(4,nmxhep) ! vertex [mm]
1399* ----------------------------------------------------------------------
1402 $ qedrad(nmxhep) ! Photos flag
1403* ----------------------------------------------------------------------
1406C position of taus, must be defined by host program:
1407 COMMON /TAUPOS/ NP1,NP2
1409 REAL PHOT1(4),PHOT2(4)
1410 REAL*8 R,X(4),Y(4),PI0(4)
1411 INTEGER JEZELI(3),ION(3)
1414.EQ.
IF (MODE-1) THEN
1420.EQ.
IF (JEZELI(1)0) RETURN
1421.EQ.
IF (JEZELI(2)1) CALL TAUETA(JAK)
1422.EQ.
IF (JEZELI(3)1) CALL TAUK0S(JAK)
1423C position of decaying particle:
1424.EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN
1429 nhepM=nhep ! to avoid infinite loop
1430 DO K=JDAHEP(1,NPS),nhepM ! we search for pi0's from tau till eor.
1431 IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k)
THEN
1436 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1445 CALL bostdq(-1,pi0,x,x)
1446 CALL bostdq(-1,pi0,y,y)
1452 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1453 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1458 SUBROUTINE taueta(JAK)
1459c
subroutine to decay etas
's from taus.
1460C this routine belongs to tauola universal interface, but uses
1461C routines from tauola utilities. Just flat phase space, but 4 channels.
1462C it is called at the beginning of SUBR. TAUPI0(JAK)
1463C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1464C this is the hepevt class in old style. No d_h_ class pre-name
1466 PARAMETER (NMXHEP=10000)
1467 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1468 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1471 $ nevhep, ! serial number
1472 $ nhep, ! number of particles
1473 $ isthep(nmxhep), ! status code
1474 $ idhep(nmxhep), ! particle ident KF
1475 $ jmohep(2,nmxhep), ! parent particles
1476 $ jdahep(2,nmxhep), ! childreen particles
1477 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1478 $ vhep(4,nmxhep) ! vertex [mm]
1479* ----------------------------------------------------------------------
1482 $ qedrad(nmxhep) ! Photos flag
1483* ----------------------------------------------------------------------
1485 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1486 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1487 * ,AMK,AMKZ,AMKST,GAMKST
1489 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1490 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1491 * ,AMK,AMKZ,AMKST,GAMKST
1493C position of taus, must be defined by host program:
1494 COMMON /TAUPOS/ NP1,NP2
1496 REAL RRR(1),BRSUM(3), RR(2)
1497 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1498 REAL*8 X(4), Y(4), Z(4)
1500 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1502 XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1503C position of decaying particle:
1504.EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN
1509 nhepM=nhep ! to avoid infinite loop
1510 DO K=JDAHEP(1,NPS),nhepM ! we search for etas's from tau till eor.
1511 IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k)
THEN
1515c eta cumulated branching ratios:
1517 brsum(2)=brsum(1)+0.319
1518 brsum(3)=brsum(2)+0.237
1521 IF (rrr(1).LT.brsum(1))
THEN
1523 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1532 CALL bostdq(-1,peta,x,x)
1533 CALL bostdq(-1,peta,y,y)
1539 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1540 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1542 IF(rrr(1).LT.brsum(2))
THEN
1549 ELSEIF(rrr(1).LT.brsum(3))
THEN
1566 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1569 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1570c weight for flat phase space
1571 wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1572 & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1574 IF (rr(2).GT.wt)
GOTO 7
1576 ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2
1580 x(4)=sqrt(ru**2+xm1**2)
1581 y(4)=sqrt(ru**2+xm2**2)
1586c generate momentum of that pair in rest frame of eta:
1587 ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1589 z(4)=sqrt(ru**2+am2**2)
1590c and boost first two decay products to rest frame of eta.
1591 CALL bostdq(-1,z,x,x)
1592 CALL bostdq(-1,z,y,y)
1593c redefine z(4) to 4-momentum of the last decay product:
1597 z(4)=sqrt(ru**2+xm3**2)
1598c boost all to lab and move to real*4; also masses
1599 CALL bostdq(-1,peta,x,x)
1600 CALL bostdq(-1,peta,y,y)
1601 CALL bostdq(-1,peta,z,z)
1611 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1612 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1613 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1620 SUBROUTINE tauk0s(JAK)
1621c
subroutine to decay k0s
's from taus.
1622C this routine belongs to tauola universal interface, but uses
1623C routines from tauola utilities. Just flat phase space, but 4 channels.
1624C it is called at the beginning of SUBR. TAUPI0(JAK)
1625C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005
1626C this is the hepevt class in old style. No d_h_ class pre-name
1628 PARAMETER (NMXHEP=10000)
1629 REAL*8 phep, vhep ! to be real*4/ *8 depending on host
1630 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1633 $ nevhep, ! serial number
1634 $ nhep, ! number of particles
1635 $ isthep(nmxhep), ! status code
1636 $ idhep(nmxhep), ! particle ident KF
1637 $ jmohep(2,nmxhep), ! parent particles
1638 $ jdahep(2,nmxhep), ! childreen particles
1639 $ phep(5,nmxhep), ! four-momentum, mass [GeV]
1640 $ vhep(4,nmxhep) ! vertex [mm]
1641* ----------------------------------------------------------------------
1644 $ qedrad(nmxhep) ! Photos flag
1645* ----------------------------------------------------------------------
1648 COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1649 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1650 * ,AMK,AMKZ,AMKST,GAMKST
1652 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1653 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1654 * ,AMK,AMKZ,AMKST,GAMKST
1656C position of taus, must be defined by host program:
1657 COMMON /TAUPOS/ NP1,NP2
1659 REAL RRR(1),BRSUM(3), RR(2)
1660 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1661 REAL*8 X(4), Y(4), Z(4)
1663 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1665 XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
1666C position of decaying particle:
1667.EQ..OR..EQ.
IF((KTO 1)(KTO11)) THEN
1672 nhepM=nhep ! to avoid infinite loop
1673 DO K=JDAHEP(1,NPS),nhepM ! we search for K0S's from tau till eor.
1674 IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k)
THEN
1680c k0s cumulated branching ratios:
1683 brsum(3)=brsum(2)+0.237
1686 IF(rrr(1).LT.brsum(1))
THEN
1691 ELSEIF(rrr(1).LT.brsum(2))
THEN
1704 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1706 r=sqrt(abs(r**2-xm1**2))
1715 CALL bostdq(-1,peta,x,x)
1716 CALL bostdq(-1,peta,y,y)
1725 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1726 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)