SUBROUTINE TAUOLA(MODE,KEYPOL) C ************************************* C general tauola interface, should work in every case until C hepevt is OK, does not check if hepevt is 'clean' C in particular will decay decayed taus... C only longitudinal spin effects are included. C in W decay v-a vertex is assumed C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001 #include "../../include/HEPEVT.h" COMMON /TAUPOS/ NP1, NP2 REAL*4 PHOI(4),PHOF(4) double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4) COMMON / MOMDEC / Q1,Q2,P1,P2,P3,P4 * tauola, photos and jetset overall switches COMMON /LIBRA/ JAK1,JAK2,ITDKRC,IFPHOT,IFHADM,IFHADP REAL*4 RRR(1) LOGICAL IFPSEUDO common /pseudocoup/ csc,ssc REAL*4 csc,ssc save pseudocoup COMMON / INOUT / INUT,IOUT C to switch tau polarization OFF in taus DIMENSION POL1(4), POL2(4) double precision POL1x(4), POL2x(4) INTEGER ION(3) DATA POL1 /0.0,0.0,0.0,0.0/ DATA POL2 /0.0,0.0,0.0,0.0/ DATA PI /3.141592653589793238462643D0/ C store decay vertexes DIMENSION IMOTHER (20) INTEGER KFHIGGS(3) C store daughter pointers INTEGER ISON COMMON /ISONS_TAU/ISON(2) SAVE /ISONS_TAU/ IF(MODE.EQ.-1) THEN C *********************** JAK1 = 0 ! decay mode first tau JAK2 = 0 ! decay mode second tau ITDKRC=1.0 ! switch of radiative corrections in decay IFPHOT=1.0 ! PHOTOS switch IFHADM=1.0 IFHADP=1.0 POL=1.0 ! tau polarization dipswitch must be 1 or 0 KFHIGGS(1) = 25 KFHIGGS(2) = 35 KFHIGGS(3) = 36 KFHIGCH = 37 KFZ0 = 23 KFGAM = 22 KFTAU = 15 KFNUE = 16 C couplings of the 'pseudoscalar higgs' as in CERN-TH/2003-166 psi=0.5*PI ! 0.15*PI xmtau=1.777 ! tau mass xmh=120 ! higgs boson mass betah=sqrt(1d0-4*xmtau**2/xmh**2) csc=cos(psi)*betah ssc=sin(psi) C write(*,*) ' scalar component=',csc,' pseudo-scalar component=',ssc IF (IFPHOT.EQ.1) CALL PHOINI ! this if PHOTOS was not initialized earlier CALL INIETC(JAK1,JAK2,ITDKRC,IFPHOT) CALL INIMAS CALL INIPHX(0.01d0) CALL INITDK C activation of pi0 and eta decays: (-1,1) means on, (-1,0) off ION(1)=1 ION(2)=1 ION(3)=1 CALL TAUPI0(-1,1,ION) CALL DEKAY(-1,POL1x) WRITE(IOUT,7001) pol,psi,ION(1),ION(2),ION(3) ELSEIF(MODE.EQ.0) THEN C *********************** C C..... find tau-s and fill common block /TAUPOS/ C this is to avoid LUND history fillings. This call is optional CALL PHYFIX(NSTOP,NSTART) C clear mothers of the previous event DO II=1,20 IMOTHER(II)=0 ENDDO DO II=1,2 ISON(II)=0 ENDDO C ... and to find mothers giving taus. NDEC = 0 C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0) DO I=NSTART,NHEP IF(ABS(IDHEP(I)).EQ.KFTAU.AND.ISTHEP(I).EQ.1.AND. $ (ISTHEP(I).GE.125.OR.ISTHEP(I).LT.120)) THEN IMOTH=JMOHEP(1,I) DO WHILE (ABS(IDHEP(IMOTH)).EQ.KFTAU) ! KEEP WALKING UP IMOTH=JMOHEP(1,IMOTH) ENDDO IF (ISTHEP(IMOTH).EQ.3.OR. $ (ISTHEP(IMOTH).GE.120.AND.ISTHEP(IMOTH).LE.125)) THEN DO J=NSTART,NHEP ! WE HAVE WALKED INTO HARD RECORD IF (IDHEP(J).EQ.IDHEP(IMOTH).AND. $ JMOHEP(1,J).EQ.IMOTH.AND. $ ISTHEP(J).EQ.2) THEN JMOTH=J GOTO 66 ENDIF ENDDO ELSE JMOTH=IMOTH ENDIF 66 CONTINUE DO II=1,NDEC IF(JMOTH.EQ.IMOTHER(II)) GOTO 9999 ENDDO C(BPK)--< NDEC=NDEC+1 IMOTHER(NDEC)= JMOTH ENDIF 9999 CONTINUE ENDDO C ... taus of every mother are treated in this main loop DO II=1,NDEC IM=IMOTHER(II) NCOUNT=0 NP1=0 NP2=0 C(BPK)--> C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT.. IM0=IM IF (IDHEP(JMOHEP(1,IM0)).EQ.IDHEP(IM0)) IM0=JMOHEP(1,IM0) ISEL=-1 DO I=NSTART,NHEP IF (ISTHEP(I).EQ.3.OR. $ (ISTHEP(I).GE.120.AND.ISTHEP(I).LE.125)) THEN ! HARD RECORD GOTO 76 ENDIF IMOTH=JMOHEP(1,I) DO WHILE (IDHEP(IMOTH).EQ.IDHEP(I).OR.ABS(IDHEP(IMOTH)).EQ.KFTAU) ! KEEP WALKING UP IMOTH=JMOHEP(1,IMOTH) ENDDO IF ((IMOTH.EQ.IM0.OR.IMOTH.EQ.IM).AND.ISEL.EQ.-1) THEN ISON(1)=I ISEL=0 ELSEIF ((IMOTH.EQ.IM0.OR.IMOTH.EQ.IM).AND.ISEL.EQ.0) THEN ISON(2)=I ELSEIF ((IMOTH.NE.IM0.AND.IMOTH.NE.IM).AND.ISEL.EQ.0) THEN ISEL=1 GOTO 77 ENDIF 76 CONTINUE ENDDO 77 CONTINUE C(BPK)--< C ... we correct HEPEVT (fix developped with Catherine BISCARAT) c IF (JDAHEP(2,IM).EQ.0) THEN ! ID of second daughter was missing c ISECU=1 c DO I=JDAHEP(1,IM)+1,NHEP ! OK lets look for it c IF (JMOHEP(1,I).EQ.IM.AND.ISECU.EQ.1) THEN ! we have found one c JDAHEP(2,IM)=I c ELSEIF (JMOHEP(1,I).EQ.IM.AND.ISECU.NE.1) THEN ! we have found one after there c JDAHEP(2,IM)=0 ! was something else, lets kill game c ENDIF c IF (JMOHEP(1,I).NE.IM) ISECU=0 ! other stuff starts c ENDDO c ENDIF C ... we check whether there are just two or more tau-likes DO I=ISON(1),ISON(2) IF(IDHEP(I).EQ.-KFTAU.OR.IDHEP(I).EQ.-KFNUE) NCOUNT=NCOUNT+1 IF(IDHEP(I).EQ. KFTAU.OR.IDHEP(I).EQ. KFNUE) NCOUNT=NCOUNT+1 ENDDO C ... if there will be more we will come here again 666 CONTINUE C(BPK)--> DO I=MAX(NP1+1,ISON(1)),ISON(2) C(BPK)--< IF(IDHEP(I).EQ.-KFTAU.OR.IDHEP(I).EQ.-KFNUE) NP1=I ENDDO C(BPK)--> DO I=MAX(NP2+1,ISON(1)),ISON(2) C(BPK)--< IF(IDHEP(I).EQ. KFTAU.OR.IDHEP(I).EQ. KFNUE) NP2=I ENDDO DO I=1,4 P1(I)= PHEP(I,NP1) !momentum of tau+ P2(I)= PHEP(I,NP2) !momentum of tau- Q1(I)= P1(I)+P2(I) ENDDO POL1(3)= 0D0 POL2(3)= 0D0 IF(KEYPOL.EQ.1) THEN c.....include polarisation effect CALL RANMAR(RRR,1) IF(IDHEP(IM).EQ.KFHIGGS(1).OR.IDHEP(IM).EQ.KFHIGGS(2).OR. $ IDHEP(IM).EQ.KFHIGGS(3)) THEN ! case of Higgs IF(RRR(1).LT.0.5) THEN POL1(3)= POL POL2(3)=-POL ELSE POL1(3)=-POL POL2(3)= POL ENDIF ELSEIF((IDHEP(IM).EQ.KFZ0).OR.(IDHEP(IM).EQ.KFGAM)) THEN ! case of gamma/Z C there is no angular dependence in gamma/Z polarization C there is no s-dependence in gamma/Z polarization at all C there is even no Z polarization in any form C main reason is that nobody asked ... C but it is prepared and longitudinal correlations C can be included up to KORALZ standards POLZ0=PLZAPX(.true.,IM,NP1,NP2) IF(RRR(1).LT.POLZ0) THEN POL1(3)= POL POL2(3)= POL ELSE POL1(3)=-POL POL2(3)=-POL ENDIF ELSEIF(IDHEP(NP1).EQ.-IDHEP(NP2))THEN ! undef orig. only s-dep poss. POLZ0=PLZAPX(.true.,IM,NP1,NP2) IF(RRR(1).LT.POLZ0) THEN POL1(3)= POL POL2(3)= POL ELSE POL1(3)=-POL POL2(3)=-POL ENDIF ELSEIF(ABS(IDHEP(IM)).EQ.KFHIGCH) THEN ! case of charged Higgs POL1(3)= POL POL2(3)= POL ELSE ! case of W+ or W- POL1(3)= -POL POL2(3)= -POL ENDIF c.....include polarisation effect ENDIF IF(IDHEP(IM).EQ.KFHIGGS(1).OR.IDHEP(IM).EQ.KFHIGGS(2).OR. $ IDHEP(IM).EQ.KFHIGGS(3)) THEN IF(IDHEP(NP1).EQ.-KFTAU .AND. $ (JDAHEP(1,NP1).LE.NP1.OR.JDAHEP(1,NP1).GT.NHEP) .AND. $ IDHEP(NP2).EQ. KFTAU .AND. $ (JDAHEP(1,NP2).LE.NP2.OR.JDAHEP(1,NP2).GT.NHEP) $ ) THEN IF (IDHEP(IM).EQ.KFHIGGS(1)) THEN IFPSEUDO= .FALSE. ELSEIF (IDHEP(IM).EQ.KFHIGGS(2)) THEN IFPSEUDO= .FALSE. ELSEIF (IDHEP(IM).EQ.KFHIGGS(3)) THEN IFPSEUDO= .TRUE. ELSE WRITE(*,*) 'Warning from TAUOLA:' WRITE(*,*) 'I stop this run, wrong IDHEP(IM)=',IDHEP(IM) STOP ENDIF CALL SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2) IF (IFPHOT.EQ.1) CALL PHOTOS(IM) ! Bremsstrahlung in Higgs decay ! AFTER adding taus !! ENDIF ELSE IF(IDHEP(NP1).EQ.-KFTAU.AND. $ (JDAHEP(1,NP1).LE.NP1.OR.JDAHEP(1,NP1).GT.NHEP)) THEN C here check on if NP1 was not decayed should be verified CALL DEXAY(1,POL1) IF (IFPHOT.EQ.1) CALL PHOTOS(NP1) CALL TAUPI0(0,1,ION) ENDIF IF(IDHEP(NP2).EQ. KFTAU.AND. $ (JDAHEP(1,NP2).LE.NP2.OR.JDAHEP(1,NP2).GT.NHEP)) THEN C here check on if NP2 was not decayed should be added CALL DEXAY(2,POL2) IF (IFPHOT.EQ.1) CALL PHOTOS(NP2) CALL TAUPI0(0,2,ION) ENDIF ENDIF NCOUNT=NCOUNT-2 IF (NCOUNT.GT.0) GOTO 666 ENDDO ELSEIF(MODE.EQ.1) THEN C *********************** C CALL DEXAY(100,POL1) CALL DEKAY(100,POL1x) WRITE(IOUT,7002) ENDIF C ***** 7001 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'*****TAUOLA UNIVERSAL INTERFACE: ******',9X,1H*, $ /,' *', 25X,'*****VERSION 1.22, April 2009 (gfort)**',9X,1H*, $ /,' *', 25X,'**AUTHORS: P. Golonka, B. Kersevan, ***',9X,1H*, $ /,' *', 25X,'**T. Pierzchala, E. Richter-Was, ******',9X,1H*, $ /,' *', 25X,'****** Z. Was, M. Worek ***************',9X,1H*, $ /,' *', 25X,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9X,1H*, $ /,' *', 25X,'*WITH C. Biscarat and S. Slabospitzky**',9X,1H*, $ /,' *', 25X,'****** are warmly acknowledged ********',9X,1H*, $ /,' *', 25X,' ',9X,1H*, $ /,' *', 25X,'********** INITIALIZATION ************',9X,1H*, $ /,' *',F20.5,5X,'tau polarization switch must be 1 or 0 ',9X,1H*, $ /,' *',F20.5,5X,'Higs scalar/pseudo mix CERN-TH/2003-166',9X,1H*, $ /,' *',I10, 15X,'PI0 decay switch must be 1 or 0 ',9X,1H*, $ /,' *',I10, 15X,'ETA decay switch must be 1 or 0 ',9X,1H*, $ /,' *',I10, 15X,'K0S decay switch must be 1 or 0 ',9X,1H*, $ /,1X,15(5H*****)/) 7002 FORMAT(///1X,15(5H*****) $ /,' *', 25X,'*****TAUOLA UNIVERSAL INTERFACE: ******',9X,1H*, $ /,' *', 25X,'*****VERSION 1.22, April 2009 (gfort)**',9X,1H*, $ /,' *', 25X,'**AUTHORS: P. Golonka, B. Kersevan, ***',9X,1H*, $ /,' *', 25X,'**T. Pierzchala, E. Richter-Was, ******',9X,1H*, $ /,' *', 25X,'****** Z. Was, M. Worek ***************',9X,1H*, $ /,' *', 25X,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9X,1H*, $ /,' *', 25X,'*WITH C. Biscarat and S. Slabospitzky**',9X,1H*, $ /,' *', 25X,'****** are warmly acknowledged ********',9X,1H*, $ /,' *', 25X,'****** END OF MODULE OPERATION ********',9X,1H*, $ /,1X,15(5H*****)/) END SUBROUTINE SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2) LOGICAL IFPSEUDO REAL*8 HH1,HH2,wthiggs DIMENSION POL1(4), POL2(4),HH1(4),HH2(4), RRR(1) ! CALL DEXAY(1,POL1) ! Kept for tests ! CALL DEXAY(2,POL2) ! Kept for tests INTEGER ION(3) 10 CONTINUE CALL RANMAR(RRR,1) CALL DEKAY(1,HH1) CALL DEKAY(2,HH2) wt=wthiggs(IFPSEUDO,HH1,HH2) IF (RRR(1).GT.WT) GOTO 10 CALL DEKAY(1+10,HH1) CALL TAUPI0(0,1,ION) CALL DEKAY(2+10,HH2) CALL TAUPI0(0,2,ION) END FUNCTION wthiggs(IFPSEUDO,HH1,HH2) LOGICAL IFPSEUDO common /pseudocoup/ csc,ssc REAL*4 csc,ssc save pseudocoup REAL*8 HH1(4),HH2(4),R(4,4),wthiggs DO K=1,4 DO L=1,4 R(K,L)=0 ENDDO ENDDO WTHIGGS=0D0 R(4,4)= 1D0 ! unpolarized part R(3,3)=-1D0 ! longitudinal ! other missing IF (IFPSEUDO) THEN R(1,1)=-1 R(2,2)= -1 R(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2) R(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2) R(1,2)=2*csc*ssc/(csc**2+ssc**2) R(2,1)=-2*csc*ssc/(csc**2+ssc**2) ELSE R(1,1)=1 R(2,2)=1 ENDIF DO K=1,4 DO L=1,4 WTHIGGS=WTHIGGS+R(K,L)*HH1(K)*HH2(L) ENDDO ENDDO WTHIGGS=WTHIGGS/4D0 END FUNCTION PLZAPX(HOPEin,IM0,NP1,NP2) C IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block. C the purpose of this routine is to calculate polarization of Z along tau direction. C this is highly non-trivial due to necessity of reading infromation from hard process C history in HEPEVT, which is often written not up to the gramatic rules. REAL*8 PLZAP0,SVAR,COSTHE,sini,sfin,ZPROP2, $ P1(4),P2(4),Q1(4),Q2(4),QQ(4),PH(4),PD1(4),PD2(4), $ PQ1(4),PQ2(4),PB(4),PA(4) REAL*4 PLZAPX INTEGER IM LOGICAL HOPE,HOPEin #include "../../include/HEPEVT.h" C(BPK)--> BROTHERS OF TAU ALREADY FOUND INTEGER ISON COMMON /ISONS_TAU/ISON(2) C(BPK)--< C >> C >> STEP 1: find where are particles in hepevent and pick them up C >> HOPE=HOPEin C sometimes shade Z of Z is its mother ... IM=IM0 IM00=JMOHEP(1,IM0) C to protect against check on mother of beam particles. IF (IM00.GT.0) THEN IF (IDHEP(IM0).EQ.IDHEP(IM00)) IM=JMOHEP(1,IM0) ENDIF C C find (host generator-level) incoming beam-bare-particles which form Z and co. IMO1=JMOHEP(1,IM) IMO2=JMOHEP(2,IM) C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS IM00=IMO1 IF (ISTHEP(IM00).EQ.120) THEN IMO1=JMOHEP(1,IM00) IMO2=JMOHEP(2,IM00) ENDIF C(BPK)--< IFFULL=0 C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt. IF (IMO1.EQ.0.AND.IMO2.EQ.0) THEN IMO1=JMOHEP(1,NP1) C(BPK)--> IF (IDHEP(IMO1).EQ.IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES C(BPK)--< IMO2=JMOHEP(2,NP1) C(BPK)--> IF (IDHEP(IMO2).EQ.IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES C(BPK)--< IFFULL=1 C case when it was like qq --> tau+tau- gammas and qq were NOT 'first' in hepevt. ELSEIF (IDHEP(IM).NE.22.AND.IDHEP(IM).NE.23) THEN IMO1=JMOHEP(1,NP1) C(BPK)--> IF (IDHEP(IMO1).EQ.IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES C(BPK)--< IMO2=JMOHEP(2,NP1) C(BPK)--> IF (IDHEP(IMO2).EQ.IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES C(BPK)--< IFFULL=1 ENDIF C and check if it really happened IF (IMO1.EQ.0) HOPE=.FALSE. IF (IMO2.EQ.0) HOPE=.FALSE. IF (IMO1.EQ.IMO2) HOPE=.FALSE. C DO I=1,4 Q1(I)= PHEP(I,NP1) !momentum of tau+ Q2(I)= PHEP(I,NP2) !momentum of tau- ENDDO C corrections due to possible differences in 4-momentum of shadow vs true Z. C(BPK)--> IF (IM.EQ.JMOHEP(1,IM0).AND. $ (IDHEP(IM).EQ.22.OR.IDHEP(IM).EQ.23)) THEN DO K=1,4 PB(K)=PHEP(K,IM) PA(K)=PHEP(K,IM0) ENDDO C(BPK)--< CALL BOSTDQ( 1,PA, Q1, Q1) CALL BOSTDQ( 1,PA, Q2, Q2) CALL BOSTDQ(-1,PB, Q1, Q1) CALL BOSTDQ(-1,PB, Q2, Q2) ENDIF DO I=1,4 QQ(I)= Q1(I)+Q2(I) !momentum of Z IF (HOPE) P1(I)=PHEP(I,IMO1) !momentum of beam1 IF (HOPE) P2(I)=PHEP(I,IMO2) !momentum of beam2 PH(I)=0D0 PD1(I)=0D0 PD2(I)=0D0 ENDDO ! These momenta correspond to quarks, gluons photons or taus IDFQ1=IDHEP(NP1) IDFQ2=IDHEP(NP2) IF (HOPE) IDFP1=IDHEP(IMO1) IF (HOPE) IDFP2=IDHEP(IMO2) SVAR=QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2 IF (.NOT.HOPE) THEN C options which may be useful in some cases of two heavy boson production C need individual considerations. To be developed. C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z PLZAPX=0.5 ! pure gamma RETURN ENDIF C >> C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles C >> C let us define beginning and end of particles which are produced in parallel to Z C in principle following should work C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS NX1=JDAHEP(1,IM00) NX2=JDAHEP(2,IM00) C but ... INBR=IM ! OK, HARD RECORD Z/GAMMA IF (IFFULL.EQ.1) INBR=NP1 ! OK, NO Z/GAMMA IF (IDHEP(JMOHEP(1,INBR)).EQ.IDHEP(INBR)) INBR=JMOHEP(1,INBR) ! FORCE HARD RECORD C(BPK)--< IF(NX1.EQ.0.OR.NX2.EQ.0) THEN NX1=INBR NX2=INBR DO K=1,INBR-1 IF(JMOHEP(1,INBR-K).EQ.JMOHEP(1,INBR)) THEN NX1=INBR-K ELSE GOTO 7 ENDIF ENDDO 7 CONTINUE DO K=INBR+1,NHEP IF(JMOHEP(1,K).EQ.JMOHEP(1,INBR)) THEN NX2=K ELSE GOTO 8 ENDIF ENDDO 8 CONTINUE ENDIF C case of annihilation of two bosons is hopeles IF (ABS(IDFP1).GE.20.AND.ABS(IDFP2).GE.20) HOPE=.FALSE. C case of annihilation of non-matching flavors is hopeless IF (ABS(IDFP1).LE.20.AND.ABS(IDFP2).LE.20.AND.IDFP1+IDFP2.NE.0) $ HOPE=.FALSE. IF (.NOT.HOPE) THEN C options which may be useful in some cases of two heavy boson production C need individual considerations. To be developed. C PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0) ! gamma/Z mixture as if produced from e beam C PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0) ! pure Z PLZAPX=0.5 ! pure gamma RETURN ENDIF IF (ABS(IDFP1).LT.20) IDE= IDFP1 IF (ABS(IDFP2).LT.20) IDE=-IDFP2 C >> C >> STEP 3 we combine gluons, photons into incoming effective beams C >> C in the following we will ignore the possibility of photon emission from taus C however at certain moment it will be necessary to take care of DO L=1,4 PD1(L)=P1(L) PD2(L)=P2(L) ENDDO DO L=1,4 PQ1(L)=Q1(L) PQ2(L)=Q2(L) ENDDO IFLAV=min(ABS(IDFP1),ABS(IDFP2)) *-------------------------------------------------------------------------- c IFLAV=min(ABS(IDFP1),ABS(IDFP2)) c that means that always quark or lepton i.e. process like c f g(gamma) --> f Z0 --> tau tau c we glue fermions to effective beams that is f f --> Z0 --> tau tau c with gamma/g emission from initial fermion. *--------------------------------------------------------------------------- IF (ABS(IDFP1).GE.20) THEN DO k=NX1,NX2 IDP=IDHEP(k) IF (ABS(IDP).EQ.IFLAV) THEN DO L=1,4 PD1(L)=-PHEP(L,K) ENDDO ENDIF ENDDO ENDIF IF (ABS(IDFP2).GE.20) THEN DO k=NX1,NX2 IDP=IDHEP(k) IF (ABS(IDP).EQ.IFLAV) THEN DO L=1,4 PD2(L)=-PHEP(L,K) ENDDO ENDIF ENDDO ENDIF C if first beam was boson: gluing IF (ABS(IDFP1).GE.20) THEN DO L=1,4 PH(L)=P1(L) ENDDO xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2 $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2) xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2 $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2) IF (XM1.LT.XM2) THEN DO L=1,4 PD1(L)=PD1(L)+P1(L) ENDDO ELSE DO L=1,4 PD2(L)=PD2(L)+P1(L) ENDDO ENDIF ENDIF C if second beam was boson: gluing IF (ABS(IDFP2).GE.20) THEN DO L=1,4 PH(L)=P2(L) ENDDO xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2 $ -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2) xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2 $ -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2) IF (XM1.LT.XM2) THEN DO L=1,4 PD1(L)=PD1(L)+P2(L) ENDDO ELSE DO L=1,4 PD2(L)=PD2(L)+P2(L) ENDDO ENDIF ENDIF C now spectators ... C(BPK)--> NPH1=NP1 NPH2=NP2 IF (IDHEP(JMOHEP(1,NP1)).EQ.IDHEP(NP1)) NPH1=JMOHEP(1,NP1) ! SHOULD PUT US IN HARD REC. IF (IDHEP(JMOHEP(1,NP2)).EQ.IDHEP(NP2)) NPH2=JMOHEP(1,NP2) ! SHOULD PUT US IN HARD REC. C(BPK)--< DO k=NX1,NX2 IF (ABS(IDHEP(K)).NE.IFLAV.AND.K.NE.IM.AND. C(BPK)--> $ K.NE.NPH1.AND.K.NE.NPH2) THEN C(BPK)--< IF(IDHEP(K).EQ.22.AND.IFFULL.EQ.1) THEN DO L=1,4 PH(L)=PHEP(L,K) ENDDO xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2 $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2) xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2 $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2) xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2 $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2) xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2 $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2) sini=abs((PD1(4)+PD2(4)-PH(4))**2-(PD1(3)+PD2(3)-PH(3))**2 $ -(PD1(2)+PD2(2)-PH(2))**2-(PD1(1)+PD2(1)-PH(1))**2) sfin=abs((PD1(4)+PD2(4) )**2-(PD1(3)+PD2(3) )**2 $ -(PD1(2)+PD2(2) )**2-(PD1(1)+PD2(1) )**2) FACINI=ZPROP2(sini) FACFIN=ZPROP2(sfin) XM1=XM1/FACINI XM2=XM2/FACINI XM3=XM3/FACFIN XM4=XM4/FACFIN XM=MIN(XM1,XM2,XM3,XM4) IF (XM1.EQ.XM) THEN DO L=1,4 PD1(L)=PD1(L)-PH(L) ENDDO ELSEIF (XM2.EQ.XM) THEN DO L=1,4 PD2(L)=PD2(L)-PH(L) ENDDO ELSEIF (XM3.EQ.XM) THEN DO L=1,4 Q1(L)=PQ1(L)+PH(L) ENDDO ELSE DO L=1,4 Q2(L)=PQ2(L)+PH(L) ENDDO ENDIF ELSE DO L=1,4 PH(L)=PHEP(L,K) ENDDO xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2 $ -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2) xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2 $ -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2) IF (XM1.LT.XM2) THEN DO L=1,4 PD1(L)=PD1(L)-PH(L) ENDDO ELSE DO L=1,4 PD2(L)=PD2(L)-PH(L) ENDDO ENDIF ENDIF ENDIF ENDDO C >> C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in c >> effective outcoming taus C >> C let us define beginning and end of particles which are produced in c parallel to tau C find outcoming particles which come from Z C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER IF (ABS(IDHEP(IM0)).EQ.22.OR.abs(IDHEP(IM0)).EQ.23) THEN DO K=ISON(1),ISON(2) IF(ABS(IDHEP(K)).EQ.22) THEN C(BPK)--< do l=1,4 ph(l)=phep(l,k) enddo xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2 $ -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2) xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2 $ -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2) XM=MIN(XM3,XM4) IF (XM3.EQ.XM) THEN DO L=1,4 Q1(L)=PQ1(L)+PH(L) ENDDO ELSE DO L=1,4 Q2(L)=PQ2(L)+PH(L) ENDDO ENDIF endif enddo ENDIF *------------------------------------------------------------------------ C out of effective momenta we calculate COSTHE and later polarization CALL ANGULU(PD1,PD2,Q1,Q2,COSTHE) PLZAPX=PLZAP0(IDE,IDFQ1,SVAR,COSTHE) END SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE) REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4) C take effective beam which is less massive, it should be irrelevant C but in case HEPEVT is particulary dirty may help. C this routine calculate reduced system transver and cosine of scattering C angle. XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2) XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2) IF (XM1.LT.XM2) THEN SIGN=1D0 DO K=1,4 P(K)=PD1(K) ENDDO ELSE SIGN=-1D0 DO K=1,4 P(K)=PD2(K) ENDDO ENDIF C calculate space like part of P (in Z restframe) DO K=1,4 QQ(K)=Q1(k)+Q2(K) QT(K)=Q1(K)-Q2(K) ENDDO XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2) QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1) DO K=1,4 QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2 ENDDO PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1) DO K=1,4 P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2 ENDDO C calculate costhe PXP =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2) QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2) PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4) COSTHE=PXQT/PXP/QTXQT COSTHE=COSTHE*SIGN END FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0) C this function calculates probability for the helicity +1 +1 configuration C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN COSTHE=COSTH0 C >>>>> IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID C >>>>> of first beam is used by T_GIVIZ0 including sign IF (IDF.GT.0) THEN CALL INITWK(IDE,IDF,SVAR) ELSE CALL INITWK(-IDE,-IDF,SVAR) ENDIF PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0) $ /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0)) ! PLZAP0=0.5 END FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB) C ---------------------------------------------------------------------- C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER C EXAMPLE OF THE METHOD APPLIED THERE C INPUT PARAMETERS ARE: SVAR -- transfer C COSTHE -- cosine of angle between tau+ and 1st beam C TA,TB -- helicity states of tau+ tau- C C called by : BORNY, BORAS, BORNV, WAGA, WEIGHT C ---------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF REAL*8 ENE ,AMIN,AMFIN COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF & ,XUPGI ,XUPZI ,XUPGF ,XUPZF & ,NDIAG0,NDIAGA,KEYA,KEYZ & ,ITCE,JTCE,ITCF,JTCF,KOLOR REAL*8 SS,POLN,T3E,QE,T3F,QF & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2) REAL*8 SEPS1,SEPS2 C===================================================================== COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ C SWSQ = sin2 (theta Weinberg) C AMW,AMZ = W & Z boson masses respectively C AMH = the Higgs mass C AMTOP = the top mass C GAMMZ = Z0 width COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2) COMPLEX*16 XUPZFP(2),XUPZIP(2) COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2) COMPLEX*16 PROPA,PROPZ COMPLEX*16 XR,XI COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF COMPLEX*16 XTHING,XVE,XVF,XVEF DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/ DATA MODE0 /-5/ DATA IDE0 /-55/ DATA SVAR0,COST0 /-5.D0,-6.D0/ DATA PI /3.141592653589793238462643D0/ DATA SEPS1,SEPS2 /0D0,0D0/ C C MEMORIZATION ========================================================= IF ( MODE.NE.MODE0.OR.SVAR.NE.SVAR0.OR.COSTHE.NE.COST0 $ .OR.IDE0.NE.IDE)THEN C KEYGSW=1 C ** PROPAGATORS IDE0=IDE MODE0=MODE SVAR0=SVAR COST0=COSTHE SINTHE=SQRT(1.D0-COSTHE**2) BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR)) C I MULTIPLY AXIAL COUPLING BY BETA FACTOR. XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2)) XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2)) XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2)) XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2)) C FINAL STATE VECTOR COUPLING XUPF =0.5D0*(XUPZF(1)+XUPZF(2)) XUPI =0.5D0*(XUPZI(1)+XUPZI(2)) XTHING =0D0 PROPA =1D0/SVAR PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ) IF (KEYGSW.EQ.0) PROPZ=0.D0 DO 50 I=1,2 DO 50 J=1,2 REGULA= (3-2*I)*(3-2*J) + COSTHE REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR) APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA) AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA ABORN(I,J)=APHOT(I,J)+AZETT(I,J) APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J) 50 CONTINUE ENDIF C C****************** C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.) C* HELICITY CONSERVATION EXPLICITLY OBEYED POLAR1= (SEPS1) POLAR2= (-SEPS2) BORN=0D0 DO 150 I=1,2 HELIC= 3-2*I DO 150 J=1,2 HELIT=3-2*J FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0 FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB) FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB) BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR C MASS TERM IN BORN IF (MODE.GE.1) THEN BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM ENDIF 150 CONTINUE C************ FUNT=BORN IF(FUNT.LT.0.D0) FUNT=BORN C IF (SVAR.GT.4D0*AMFIN**2) THEN C PHASE SPACE THRESHOLD FACTOR THRESH=SQRT(1-4D0*AMFIN**2/SVAR) T_BORN= FUNT*SVAR**2*THRESH ELSE THRESH=0.D0 T_BORN=0.D0 ENDIF C ZW HERE WAS AN ERROR 19. 05. 1989 ! write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF ! write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN END SUBROUTINE INITWK(IDEX,IDFX,SVAR) ! initialization routine coupling masses etc. IMPLICIT REAL*8 (A-H,O-Z) COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF REAL*8 ENE ,AMIN,AMFIN COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF & ,XUPGI ,XUPZI ,XUPGF ,XUPZF & ,NDIAG0,NDIAGA,KEYA,KEYZ & ,ITCE,JTCE,ITCF,JTCF,KOLOR REAL*8 SS,POLN,T3E,QE,T3F,QF & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2) COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ C SWSQ = sin2 (theta Weinberg) C AMW,AMZ = W & Z boson masses respectively C AMH = the Higgs mass C AMTOP = the top mass C GAMMZ = Z0 width C ENE=SQRT(SVAR)/2 AMIN=0.511D-3 SWSQ=0.23147 AMZ=91.1882 GAMMZ=2.4952 IF (IDFX.EQ. 15) then IDF=2 ! denotes tau +2 tau- AMFIN=1.77703 !this mass is irrelevant if small, used in ME only ELSEIF (IDFX.EQ.-15) then IDF=-2 ! denotes tau -2 tau- AMFIN=1.77703 !this mass is irrelevant if small, used in ME only ELSE WRITE(*,*) 'INITWK: WRONG IDFX' STOP ENDIF IF (IDEX.EQ. 11) then !electron IDE= 2 AMIN=0.511D-3 ELSEIF (IDEX.EQ.-11) then !positron IDE=-2 AMIN=0.511D-3 ELSEIF (IDEX.EQ. 13) then !mu+ IDE= 2 AMIN=0.105659 ELSEIF (IDEX.EQ.-13) then !mu- IDE=-2 AMIN=0.105659 ELSEIF (IDEX.EQ. 1) then !d IDE= 4 AMIN=0.05 ELSEIF (IDEX.EQ.- 1) then !d~ IDE=-4 AMIN=0.05 ELSEIF (IDEX.EQ. 2) then !u IDE= 3 AMIN=0.02 ELSEIF (IDEX.EQ.- 2) then !u~ IDE=-3 AMIN=0.02 ELSEIF (IDEX.EQ. 3) then !s IDE= 4 AMIN=0.3 ELSEIF (IDEX.EQ.- 3) then !s~ IDE=-4 AMIN=0.3 ELSEIF (IDEX.EQ. 4) then !c IDE= 3 AMIN=1.3 ELSEIF (IDEX.EQ.- 4) then !c~ IDE=-3 AMIN=1.3 ELSEIF (IDEX.EQ. 5) then !b IDE= 4 AMIN=4.5 ELSEIF (IDEX.EQ.- 5) then !b~ IDE=-4 AMIN=4.5 ELSEIF (IDEX.EQ. 12) then !nu_e IDE= 1 AMIN=0.1D-3 ELSEIF (IDEX.EQ.- 12) then !nu_e~ IDE=-1 AMIN=0.1D-3 ELSEIF (IDEX.EQ. 14) then !nu_mu IDE= 1 AMIN=0.1D-3 ELSEIF (IDEX.EQ.- 14) then !nu_mu~ IDE=-1 AMIN=0.1D-3 ELSEIF (IDEX.EQ. 16) then !nu_tau IDE= 1 AMIN=0.1D-3 ELSEIF (IDEX.EQ.- 16) then !nu_tau~ IDE=-1 AMIN=0.1D-3 ELSE WRITE(*,*) 'INITWK: WRONG IDEX' STOP ENDIF C ---------------------------------------------------------------------- C C INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX C C called by : KORALZ C ---------------------------------------------------------------------- ITCE=IDE/IABS(IDE) JTCE=(1-ITCE)/2 ITCF=IDF/IABS(IDF) JTCF=(1-ITCF)/2 CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM) CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM) XUPGI(1)=QE XUPGI(2)=QE T3E = AIZOL+AIZOR XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR) CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR) XUPGF(1)=QF XUPGF(2)=QF T3F = AIZOL+AIZOR XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ)) XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ)) C NDIAG0=2 NDIAGA=11 KEYA = 1 KEYZ = 1 C C RETURN END SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR) C ---------------------------------------------------------------------- C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY) C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF) C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS C C called by : EVENTE, EVENTM, FUNTIH, ..... C ---------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) C IF(IDFERM.EQ.0.OR.IABS(IDFERM).GT.4) GOTO 901 IF(IABS(IHELIC).NE.1) GOTO 901 IH =IHELIC IDTYPE =IABS(IDFERM) IC =IDFERM/IDTYPE LEPQUA=INT(IDTYPE*0.4999999D0) IUPDOW=IDTYPE-2*LEPQUA-1 CHARGE =(-IUPDOW+2D0/3D0*LEPQUA)*IC SIZO3 =0.25D0*(IC-IH)*(1-2*IUPDOW) KOLOR=1+2*LEPQUA C** NOTE THAT CONVENTIONALY Z0 COUPLING IS C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ)) RETURN 901 PRINT *,' STOP IN GIVIZO: WRONG PARAMS.' STOP END #include "../../include/phyfix.h" SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG) C ---------------------------------------------------------------------- C this subroutine fills one entry into the HEPEVT common C and updates the information for affected mother entries C C written by Martin W. Gruenewald (91/01/28) C C called by : ZTOHEP,BTOHEP,DWLUxy C ---------------------------------------------------------------------- C C this is the hepevt class in old style. No d_h_ class pre-name #include "../../include/HEPEVT.h" LOGICAL PHFLAG C REAL*4 P4(4) C C check address mode IF (N.EQ.0) THEN C C append mode IHEP=NHEP+1 ELSE IF (N.GT.0) THEN C C absolute position IHEP=N ELSE C C relative position IHEP=NHEP+N END IF C C check on IHEP IF ((IHEP.LE.0).OR.(IHEP.GT.NMXHEP)) RETURN C C add entry NHEP=IHEP ISTHEP(IHEP)=IST IDHEP(IHEP)=ID JMOHEP(1,IHEP)=JMO1 IF(JMO1.LT.0)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP JMOHEP(2,IHEP)=JMO2 IF(JMO2.LT.0)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP JDAHEP(1,IHEP)=JDA1 JDAHEP(2,IHEP)=JDA2 C DO I=1,4 PHEP(I,IHEP)=P4(I) C C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations VHEP(I,IHEP)=0.0 END DO PHEP(5,IHEP)=PINV C FLAG FOR PHOTOS... QEDRAD(IHEP)=PHFLAG C C update process: DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP) IF(IP.GT.0)THEN C C if there is a daughter at IHEP, mother entry at IP has decayed IF(ISTHEP(IP).EQ.1)ISTHEP(IP)=2 C C and daughter pointers of mother entry must be updated IF(JDAHEP(1,IP).EQ.0)THEN JDAHEP(1,IP)=IHEP JDAHEP(2,IP)=IHEP ELSE JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP)) END IF END IF END DO C RETURN END FUNCTION IHEPDIM(DUM) C this is the hepevt class in old style. No d_h_ class pre-name #include "../../include/HEPEVT.h" IHEPDIM=NHEP END FUNCTION ZPROP2(S) IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 CPRZ0,CPRZ0M AMZ=91.1882 GAMMZ=2.49 CPRZ0=DCMPLX((S-AMZ**2),S/AMZ*GAMMZ) CPRZ0M=1/CPRZ0 ZPROP2=(ABS(CPRZ0M))**2 END SUBROUTINE TAUPI0(MODE,JAK,ION) C no initialization required. Must be called once after every: C 1) CALL DEKAY(1+10,...) C 2) CALL DEKAY(2+10,...) C 3) CALL DEXAY(1,...) C 4) CALL DEXAY(2,...) C subroutine to decay originating from TAUOLA's taus: C 1) etas (with CALL TAUETA(JAK)) C 2) later pi0's from taus. C 3) extensions to other applications possible. C this routine belongs to >tauola universal interface<, but uses C routines from >tauola< utilities as well. 25.08.2005 #include "../../include/HEPEVT.h" C position of taus, must be defined by host program: COMMON /TAUPOS/ NP1,NP2 c REAL PHOT1(4),PHOT2(4) REAL*8 R,X(4),Y(4),PI0(4) INTEGER JEZELI(3),ION(3) DATA JEZELI /0,0,0/ SAVE JEZELI IF (MODE.EQ.-1) THEN JEZELI(1)=ION(1) JEZELI(2)=ION(2) JEZELI(3)=ION(3) RETURN ENDIF IF (JEZELI(1).EQ.0) RETURN IF (JEZELI(2).EQ.1) CALL TAUETA(JAK) IF (JEZELI(3).EQ.1) CALL TAUK0S(JAK) C position of decaying particle: IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN NPS=NP1 ELSE NPS=NP2 ENDIF nhepM=nhep ! to avoid infinite loop DO K=JDAHEP(1,NPS),nhepM ! we search for pi0's from tau till eor. IF (IDHEP(K).EQ.111.AND.JDAHEP(1,K).LE.K) THEN ! IF we found pi0 DO L=1,4 PI0(L)= phep(L,K) ENDDO ! random 3 vector on the sphere, masless R=SQRT(PI0(4)**2-PI0(3)**2-PI0(2)**2-PI0(1)**2)/2D0 CALL SPHERD(R,X) X(4)=R Y(4)=R Y(1)=-X(1) Y(2)=-X(2) Y(3)=-X(3) ! boost to lab and to real*4 CALL bostdq(-1,PI0,X,X) CALL bostdq(-1,PI0,Y,Y) DO L=1,4 PHOT1(L)=X(L) PHOT2(L)=Y(L) ENDDO C to hepevt CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.) CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.) ENDIF ENDDO C END SUBROUTINE TAUETA(JAK) C subroutine to decay etas's from taus. C this routine belongs to tauola universal interface, but uses C routines from tauola utilities. Just flat phase space, but 4 channels. C it is called at the beginning of SUBR. TAUPI0(JAK) C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005 #include "../../include/HEPEVT.h" COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 * ,AMK,AMKZ,AMKST,GAMKST * REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 * ,AMK,AMKZ,AMKST,GAMKST C position of taus, must be defined by host program: COMMON /TAUPOS/ NP1,NP2 c REAL RRR(1),BRSUM(3), RR(2) REAL PHOT1(4),PHOT2(4),PHOT3(4) REAL*8 X(4), Y(4), Z(4) REAL YM1,YM2,YM3 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM REAL*8 a,b,c XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c)) C position of decaying particle: IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN NPS=NP1 ELSE NPS=NP2 ENDIF nhepM=nhep ! to avoid infinite loop DO K=JDAHEP(1,NPS),nhepM ! we search for etas's from tau till eor. IF (IDHEP(K).EQ.221.AND.JDAHEP(1,K).LE.K) THEN ! IF we found eta DO L=1,4 PETA(L)= phep(L,K) ! eta 4 momentum ENDDO C eta cumulated branching ratios: BRSUM(1)=0.389 ! gamma gamma BRSUM(2)=BRSUM(1)+0.319 ! 3 pi0 BRSUM(3)=BRSUM(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma CALL RANMAR(RRR,1) IF (RRR(1).LT.BRSUM(1)) THEN ! gamma gamma channel exactly like pi0 ! random 3 vector on the sphere, masless R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0 CALL SPHERD(R,X) X(4)=R Y(4)=R Y(1)=-X(1) Y(2)=-X(2) Y(3)=-X(3) ! boost to lab and to real*4 CALL bostdq(-1,PETA,X,X) CALL bostdq(-1,PETA,Y,Y) DO L=1,4 PHOT1(L)=X(L) PHOT2(L)=Y(L) ENDDO C to hepevt CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.) CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.) ELSE ! 3 body channels IF(RRR(1).LT.BRSUM(2)) THEN ! 3 pi0 ID1= 111 ID2= 111 ID3= 111 XM1=AMPIZ ! masses XM2=AMPIZ XM3=AMPIZ ELSEIF(RRR(1).LT.BRSUM(3)) THEN ! pi+ pi- pi0 ID1= 211 ID2=-211 ID3= 111 XM1=AMPI ! masses XM2=AMPI XM3=AMPIZ ELSE ! pi+ pi- gamma ID1= 211 ID2=-211 ID3= 22 XM1=AMPI ! masses XM2=AMPI XM3=0.0 ENDIF 7 CONTINUE ! we generate mass of the first pair: CALL RANMAR(RR,2) R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2) AMIN=XM1+XM2 AMAX=R-XM3 AM2=SQRT(AMIN**2+RR(1)*(AMAX**2-AMIN**2)) C weight for flat phase space WT=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2) & *XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2) & /R**2 /AM2**2 IF (RR(2).GT.WT) GOTO 7 RU=XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)/AM2/2 ! momenta of the ! first two products ! in the rest frame of that pair CALL SPHERD(RU,X) X(4)=SQRT(RU**2+XM1**2) Y(4)=SQRT(RU**2+XM2**2) Y(1)=-X(1) Y(2)=-X(2) Y(3)=-X(3) C generate momentum of that pair in rest frame of eta: RU=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)/R/2 CALL SPHERD(RU,Z) Z(4)=SQRT(RU**2+AM2**2) C and boost first two decay products to rest frame of eta. CALL bostdq(-1,Z,X,X) CALL bostdq(-1,Z,Y,Y) C redefine Z(4) to 4-momentum of the last decay product: Z(1)=-Z(1) Z(2)=-Z(2) Z(3)=-Z(3) Z(4)=SQRT(RU**2+XM3**2) C boost all to lab and move to real*4; also masses CALL bostdq(-1,PETA,X,X) CALL bostdq(-1,PETA,Y,Y) CALL bostdq(-1,PETA,Z,Z) DO L=1,4 PHOT1(L)=X(L) PHOT2(L)=Y(L) PHOT3(L)=Z(L) ENDDO YM1=XM1 YM2=XM2 YM3=XM3 C to hepevt CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.) CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.) CALL FILHEP(0,1,ID3,K,K,0,0,PHOT3,YM3,.TRUE.) ENDIF ENDIF ENDDO C END SUBROUTINE TAUK0S(JAK) C subroutine to decay K0S's from taus. C this routine belongs to tauola universal interface, but uses C routines from tauola utilities. Just flat phase space, but 4 channels. C it is called at the beginning of SUBR. TAUPI0(JAK) C and as far as hepevt search it is basically the same as TAUPI0. 25.08.2005 #include "../../include/HEPEVT.h" COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 * ,AMK,AMKZ,AMKST,GAMKST * REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1 * ,AMK,AMKZ,AMKST,GAMKST C position of taus, must be defined by host program: COMMON /TAUPOS/ NP1,NP2 c REAL RRR(1),BRSUM(3), RR(2) REAL PHOT1(4),PHOT2(4),PHOT3(4) REAL*8 X(4), Y(4), Z(4) REAL YM1,YM2,YM3 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM REAL*8 a,b,c XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c)) C position of decaying particle: IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN NPS=NP1 ELSE NPS=NP2 ENDIF nhepM=nhep ! to avoid infinite loop DO K=JDAHEP(1,NPS),nhepM ! we search for K0S's from tau till eor. IF (IDHEP(K).EQ.310.AND.JDAHEP(1,K).LE.K) THEN ! IF we found K0S DO L=1,4 PETA(L)= phep(L,K) ! K0S 4 momentum (this is cloned from eta decay) ENDDO C K0S cumulated branching ratios: BRSUM(1)=0.313 ! 2 PI0 BRSUM(2)=1.0 ! BRSUM(1)+0.319 ! Pi+ PI- BRSUM(3)=BRSUM(2)+0.237 ! pi+ pi- pi0 rest is thus pi+pi-gamma CALL RANMAR(RRR,1) IF(RRR(1).LT.BRSUM(1)) THEN ! 2 pi0 ID1= 111 ID2= 111 XM1=AMPIZ ! masses XM2=AMPIZ ELSEIF(RRR(1).LT.BRSUM(2)) THEN ! pi+ pi- ID1= 211 ID2=-211 XM1=AMPI ! masses XM2=AMPI ELSE ! gamma gamma unused !!! ID1= 22 ID2= 22 XM1= 0.0 ! masses XM2= 0.0 ENDIF ! random 3 vector on the sphere, of equal mass !! R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0 R4=R R=SQRT(ABS(R**2-XM1**2)) CALL SPHERD(R,X) X(4)=R4 Y(4)=R4 Y(1)=-X(1) Y(2)=-X(2) Y(3)=-X(3) ! boost to lab and to real*4 CALL bostdq(-1,PETA,X,X) CALL bostdq(-1,PETA,Y,Y) DO L=1,4 PHOT1(L)=X(L) PHOT2(L)=Y(L) ENDDO YM1=XM1 YM2=XM2 C to hepevt CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.) CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.) C ENDIF ENDDO END