tauola-F/jetset-F/tauface-jetset.f

00001 
00002 
00003       SUBROUTINE TAUOLA(MODE,KEYPOL) 
00004 C     *************************************
00005 C general tauola interface, should work in every case until 
00006 C hepevt is OK, does not check if hepevt is 'clean'
00007 C in particular will decay decayed taus...
00008 C only longitudinal spin effects are included.
00009 C in W decay v-a vertex is assumed
00010 C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
00011 C this is the hepevt class in old style. No d_h_ class pre-name
00012       INTEGER NMXHEP
00013       PARAMETER (NMXHEP=10000)
00014       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
00015       INTEGER nevhep,nhep,isthep,idhep,jmohep,
00016      $        jdahep
00017       COMMON /hepevt/
00018      $      nevhep,               ! serial number
00019      $      nhep,                 ! number of particles
00020      $      isthep(nmxhep),   ! status code
00021      $      idhep(nmxhep),    ! particle ident KF
00022      $      jmohep(2,nmxhep), ! parent particles
00023      $      jdahep(2,nmxhep), ! childreen particles
00024      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
00025      $      vhep(4,nmxhep)    ! vertex [mm]
00026 * ----------------------------------------------------------------------
00027       LOGICAL qedrad
00028       COMMON /phoqed/ 
00029      $     qedrad(nmxhep)    ! Photos flag
00030 * ----------------------------------------------------------------------
00031       SAVE hepevt,phoqed
00032       COMMON /TAUPOS/ NP1, NP2 
00033       REAL*4 PHOI(4),PHOF(4)
00034       double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
00035       COMMON / MOMDEC / Q1,Q2,P1,P2,P3,P4
00036 * tauola, photos and jetset overall switches
00037       COMMON /LIBRA/ JAK1,JAK2,ITDKRC,IFPHOT,IFHADM,IFHADP
00038       REAL*4 RRR(1)
00039       LOGICAL IFPSEUDO
00040       common /pseudocoup/ csc,ssc
00041       REAL*4 csc,ssc
00042       save pseudocoup
00043       COMMON / INOUT / INUT,IOUT
00044 
00045 C to switch tau polarization OFF in taus 
00046       DIMENSION POL1(4), POL2(4)
00047       double precision POL1x(4), POL2x(4)
00048       INTEGER ION(3)
00049       DATA  POL1 /0.0,0.0,0.0,0.0/
00050       DATA  POL2 /0.0,0.0,0.0,0.0/
00051       DATA PI /3.141592653589793238462643D0/
00052 
00053 C store decay vertexes 
00054       DIMENSION IMOTHER (20)
00055       INTEGER KFHIGGS(3)
00056 
00057 C store daughter pointers
00058       INTEGER ISON
00059       COMMON /ISONS_TAU/ISON(2)
00060       SAVE /ISONS_TAU/
00061 
00062       IF(MODE.EQ.-1) THEN
00063 C     ***********************
00064 
00065          JAK1  =  0     ! decay mode first tau 
00066          JAK2  =  0     ! decay mode second tau
00067          ITDKRC=1.0     ! switch of radiative corrections in decay
00068          IFPHOT=1.0     ! PHOTOS switch
00069          IFHADM=1.0
00070          IFHADP=1.0
00071          POL=1.0        ! tau polarization dipswitch must be 1 or 0
00072 
00073          KFHIGGS(1) = 25
00074          KFHIGGS(2) = 35
00075          KFHIGGS(3) = 36
00076          KFHIGCH = 37
00077          KFZ0    = 23
00078          KFGAM   = 22
00079          KFTAU   = 15
00080          KFNUE   = 16
00081 C couplings of the 'pseudoscalar higgs' as in CERN-TH/2003-166
00082          psi=0.5*PI ! 0.15*PI
00083          xmtau=1.777 ! tau mass
00084          xmh=120     ! higgs boson mass
00085          betah=sqrt(1d0-4*xmtau**2/xmh**2)
00086          csc=cos(psi)*betah
00087          ssc=sin(psi)
00088 C         write(*,*) ' scalar component=',csc,' pseudo-scalar component=',ssc
00089          IF (IFPHOT.EQ.1) CALL  PHOINI  ! this if PHOTOS was not initialized earlier
00090          CALL  INIETC(JAK1,JAK2,ITDKRC,IFPHOT)
00091          CALL  INIMAS
00092          CALL  INIPHX(0.01d0)
00093          CALL  INITDK
00094 C activation of pi0 and eta decays: (1) means on,  (0) off
00095          ION(1)=0
00096          ION(2)=0
00097          ION(3)=0
00098          CALL  TAUPI0(-1,1,ION)
00099          CALL DEKAY(-1,POL1x)
00100            WRITE(IOUT,7001) pol,psi,ION(1),ION(2),ION(3)
00101       ELSEIF(MODE.EQ.0) THEN
00102 C     ***********************
00103 C
00104 C..... find tau-s and fill common block /TAUPOS/
00105 C this is to avoid LUND history fillings. This call is optional
00106       CALL PHYFIX(NSTOP,NSTART)
00107 C clear mothers of the previous event
00108       DO II=1,20
00109        IMOTHER(II)=0
00110       ENDDO
00111 
00112       DO II=1,2
00113          ISON(II)=0
00114       ENDDO
00115 C ... and to find mothers giving taus.
00116       NDEC    = 0
00117 C(BPK)--> LOOK FOR MOTHER, CHECK THAT IT IS NOT THE HISTORY ENTRY (E.G. MSTP(128)=0)
00118       DO I=NSTART,NHEP
00119         IF(ABS(IDHEP(I)).EQ.KFTAU.AND.ISTHEP(I).EQ.1.AND.
00120      $        (ISTHEP(I).GE.125.OR.ISTHEP(I).LT.120)) THEN
00121            IMOTH=JMOHEP(1,I)
00122            DO WHILE (ABS(IDHEP(IMOTH)).EQ.KFTAU) ! KEEP WALKING UP
00123               IMOTH=JMOHEP(1,IMOTH)
00124            ENDDO
00125            IF (ISTHEP(IMOTH).EQ.3.OR.
00126      $        (ISTHEP(IMOTH).GE.120.AND.ISTHEP(IMOTH).LE.125)) THEN 
00127               DO J=NSTART,NHEP  ! WE HAVE WALKED INTO HARD RECORD
00128                  IF (IDHEP(J).EQ.IDHEP(IMOTH).AND.
00129      $               JMOHEP(1,J).EQ.IMOTH.AND.
00130      $               ISTHEP(J).EQ.2) THEN
00131                     JMOTH=J
00132                     GOTO 66
00133                  ENDIF
00134               ENDDO
00135            ELSE
00136               JMOTH=IMOTH
00137            ENDIF
00138  66        CONTINUE
00139            DO II=1,NDEC
00140             IF(JMOTH.EQ.IMOTHER(II)) GOTO 9999
00141            ENDDO
00142 C(BPK)--<
00143            NDEC=NDEC+1
00144            IMOTHER(NDEC)= JMOTH
00145         ENDIF
00146  9999   CONTINUE
00147       ENDDO  
00148 
00149 C ... taus of every mother are treated in this main loop
00150       DO II=1,NDEC
00151          IM=IMOTHER(II)
00152          NCOUNT=0
00153          NP1=0
00154          NP2=0
00155 
00156 
00157 C(BPK)--> 
00158 C CORRECTING HEPEVT IS OUT OF QUESTION AT THIS POINT..
00159          IM0=IM
00160          IF (IDHEP(JMOHEP(1,IM0)).EQ.IDHEP(IM0)) IM0=JMOHEP(1,IM0)
00161          ISEL=-1
00162          DO I=NSTART,NHEP
00163             IF (ISTHEP(I).EQ.3.OR.
00164      $           (ISTHEP(I).GE.120.AND.ISTHEP(I).LE.125)) THEN ! HARD RECORD
00165                GOTO 76
00166             ENDIF
00167             IMOTH=JMOHEP(1,I)
00168             DO WHILE (IDHEP(IMOTH).EQ.IDHEP(I).OR.ABS(IDHEP(IMOTH)).EQ.KFTAU) ! KEEP WALKING UP
00169                IMOTH=JMOHEP(1,IMOTH)
00170             ENDDO
00171             IF ((IMOTH.EQ.IM0.OR.IMOTH.EQ.IM).AND.ISEL.EQ.-1) THEN
00172                ISON(1)=I
00173                ISEL=0
00174             ELSEIF ((IMOTH.EQ.IM0.OR.IMOTH.EQ.IM).AND.ISEL.EQ.0) THEN
00175                ISON(2)=I
00176             ELSEIF ((IMOTH.NE.IM0.AND.IMOTH.NE.IM).AND.ISEL.EQ.0) THEN
00177                ISEL=1
00178                GOTO 77
00179             ENDIF
00180  76         CONTINUE
00181          ENDDO
00182  77      CONTINUE
00183 C(BPK)--< 
00184 
00185 
00186 C ... we correct HEPEVT (fix developped with Catherine BISCARAT)
00187 c         IF (JDAHEP(2,IM).EQ.0) THEN                      ! ID of second daughter was missing
00188 c          ISECU=1
00189 c          DO I=JDAHEP(1,IM)+1,NHEP                        ! OK lets look for it
00190 c           IF (JMOHEP(1,I).EQ.IM.AND.ISECU.EQ.1) THEN     ! we have found one
00191 c              JDAHEP(2,IM)=I
00192 c           ELSEIF (JMOHEP(1,I).EQ.IM.AND.ISECU.NE.1) THEN ! we have found one after there
00193 c              JDAHEP(2,IM)=0                              ! was something else, lets kill game
00194 c           ENDIF
00195 c           IF (JMOHEP(1,I).NE.IM) ISECU=0                 ! other stuff starts
00196 c          ENDDO
00197 c         ENDIF
00198 
00199 C ... we check whether there are just two or more tau-likes 
00200          DO I=ISON(1),ISON(2)
00201             IF(IDHEP(I).EQ.-KFTAU.OR.IDHEP(I).EQ.-KFNUE) NCOUNT=NCOUNT+1
00202             IF(IDHEP(I).EQ. KFTAU.OR.IDHEP(I).EQ. KFNUE) NCOUNT=NCOUNT+1
00203          ENDDO
00204 
00205 C ... if there will be more we will come here again
00206  666     CONTINUE
00207 
00208 C(BPK)--> 
00209          DO I=MAX(NP1+1,ISON(1)),ISON(2) 
00210 C(BPK)--< 
00211             IF(IDHEP(I).EQ.-KFTAU.OR.IDHEP(I).EQ.-KFNUE) NP1=I
00212          ENDDO
00213 C(BPK)--> 
00214          DO I=MAX(NP2+1,ISON(1)),ISON(2)
00215 C(BPK)--< 
00216             IF(IDHEP(I).EQ. KFTAU.OR.IDHEP(I).EQ. KFNUE) NP2=I
00217          ENDDO
00218          DO I=1,4
00219             P1(I)= PHEP(I,NP1)    !momentum of tau+
00220             P2(I)= PHEP(I,NP2)    !momentum of tau-
00221             Q1(I)= P1(I)+P2(I)
00222          ENDDO
00223 
00224          POL1(3)=  0D0
00225          POL2(3)=  0D0
00226          
00227          IF(KEYPOL.EQ.1) THEN
00228 c.....include polarisation effect
00229          CALL RANMAR(RRR,1)
00230 
00231          IF(IDHEP(IM).EQ.KFHIGGS(1).OR.IDHEP(IM).EQ.KFHIGGS(2).OR.
00232      $    IDHEP(IM).EQ.KFHIGGS(3)) THEN   ! case of Higgs 
00233             IF(RRR(1).LT.0.5) THEN
00234                POL1(3)= POL
00235                POL2(3)=-POL
00236             ELSE     
00237                POL1(3)=-POL
00238                POL2(3)= POL
00239             ENDIF
00240          ELSEIF((IDHEP(IM).EQ.KFZ0).OR.(IDHEP(IM).EQ.KFGAM)) THEN ! case of gamma/Z 
00241 C there is no angular dependence in gamma/Z polarization 
00242 C there is no s-dependence in gamma/Z polarization at all
00243 C there is even no Z polarization in any form
00244 C main reason is that nobody asked ...
00245 C but it is prepared and longitudinal correlations 
00246 C can be included up to KORALZ standards
00247 
00248             POLZ0=PLZAPX(.true.,IM,NP1,NP2)
00249             IF(RRR(1).LT.POLZ0) THEN
00250                POL1(3)= POL
00251                POL2(3)= POL
00252             ELSE     
00253                POL1(3)=-POL
00254                POL2(3)=-POL
00255             ENDIF
00256          ELSEIF(IDHEP(NP1).EQ.-IDHEP(NP2))THEN ! undef orig. only s-dep poss.
00257             POLZ0=PLZAPX(.true.,IM,NP1,NP2)
00258             IF(RRR(1).LT.POLZ0) THEN
00259                POL1(3)= POL
00260                POL2(3)= POL
00261             ELSE     
00262                POL1(3)=-POL
00263                POL2(3)=-POL
00264             ENDIF
00265          ELSEIF(ABS(IDHEP(IM)).EQ.KFHIGCH) THEN ! case of charged Higgs 
00266             POL1(3)=  POL
00267             POL2(3)=  POL
00268          ELSE ! case of W+ or W-
00269             POL1(3)= -POL
00270             POL2(3)= -POL
00271          ENDIF
00272 c.....include polarisation effect
00273          ENDIF  
00274 
00275          IF(IDHEP(IM).EQ.KFHIGGS(1).OR.IDHEP(IM).EQ.KFHIGGS(2).OR.
00276      $   IDHEP(IM).EQ.KFHIGGS(3)) THEN
00277            IF(IDHEP(NP1).EQ.-KFTAU                           .AND.
00278      $       (JDAHEP(1,NP1).LE.NP1.OR.JDAHEP(1,NP1).GT.NHEP) .AND.
00279      $        IDHEP(NP2).EQ. KFTAU                            .AND.
00280      $       (JDAHEP(1,NP2).LE.NP2.OR.JDAHEP(1,NP2).GT.NHEP)
00281      $                                                       ) THEN
00282                IF     (IDHEP(IM).EQ.KFHIGGS(1)) THEN
00283                  IFPSEUDO= .FALSE.
00284                ELSEIF (IDHEP(IM).EQ.KFHIGGS(2)) THEN
00285                  IFPSEUDO= .FALSE.
00286                ELSEIF (IDHEP(IM).EQ.KFHIGGS(3)) THEN
00287                  IFPSEUDO= .TRUE.
00288                ELSE
00289                  WRITE(*,*) 'Warning from TAUOLA:'
00290                  WRITE(*,*) 'I stop this run, wrong IDHEP(IM)=',IDHEP(IM)
00291                  STOP
00292                ENDIF
00293                CALL SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
00294                IF (IFPHOT.EQ.1) CALL PHOTOS(IM)  ! Bremsstrahlung in Higgs decay
00295                                                  ! AFTER adding taus !!
00296 
00297 
00298            ENDIF
00299          ELSE
00300            IF(IDHEP(NP1).EQ.-KFTAU.AND.
00301      $       (JDAHEP(1,NP1).LE.NP1.OR.JDAHEP(1,NP1).GT.NHEP)) THEN
00302 C            here check on if NP1 was not decayed should be verified
00303              CALL DEXAY(1,POL1)
00304              IF (IFPHOT.EQ.1) CALL PHOTOS(NP1)
00305              CALL TAUPI0(0,1,ION)
00306            ENDIF
00307 
00308            IF(IDHEP(NP2).EQ. KFTAU.AND.
00309      $       (JDAHEP(1,NP2).LE.NP2.OR.JDAHEP(1,NP2).GT.NHEP)) THEN
00310 C            here check on if NP2 was not decayed should be added
00311              CALL DEXAY(2,POL2)
00312              IF (IFPHOT.EQ.1) CALL PHOTOS(NP2)
00313              CALL TAUPI0(0,2,ION)
00314            ENDIF
00315          ENDIF
00316          NCOUNT=NCOUNT-2
00317          IF (NCOUNT.GT.0) GOTO 666
00318       ENDDO
00319 
00320       ELSEIF(MODE.EQ.1) THEN
00321 C     ***********************
00322 C
00323       CALL DEXAY(100,POL1)
00324       CALL DEKAY(100,POL1x)
00325            WRITE(IOUT,7002)
00326       ENDIF
00327 C     *****
00328  7001 FORMAT(///1X,15(5H*****)
00329      $ /,' *',     25X,'*****TAUOLA UNIVERSAL INTERFACE: ******',9X,1H*,
00330      $ /,' *',     25X,'*****VERSION 1.22, April 2009 (gfort)**',9X,1H*,
00331      $ /,' *',     25X,'**AUTHORS: P. Golonka, B. Kersevan, ***',9X,1H*,
00332      $ /,' *',     25X,'**T. Pierzchala, E. Richter-Was, ******',9X,1H*,
00333      $ /,' *',     25X,'****** Z. Was, M. Worek ***************',9X,1H*,
00334      $ /,' *',     25X,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9X,1H*,
00335      $ /,' *',     25X,'*WITH C. Biscarat and S. Slabospitzky**',9X,1H*,
00336      $ /,' *',     25X,'****** are warmly acknowledged ********',9X,1H*,
00337      $ /,' *',     25X,'                                       ',9X,1H*,
00338      $ /,' *',     25X,'********** INITIALIZATION  ************',9X,1H*,
00339      $ /,' *',F20.5,5X,'tau polarization switch must be 1 or 0 ',9X,1H*,
00340      $ /,' *',F20.5,5X,'Higs scalar/pseudo mix CERN-TH/2003-166',9X,1H*,
00341      $ /,' *',I10, 15X,'PI0 decay switch must be 1 or 0        ',9X,1H*,
00342      $ /,' *',I10, 15X,'ETA decay switch must be 1 or 0        ',9X,1H*,
00343      $ /,' *',I10, 15X,'K0S decay switch must be 1 or 0        ',9X,1H*,
00344      $  /,1X,15(5H*****)/)
00345 
00346  7002 FORMAT(///1X,15(5H*****)
00347      $ /,' *',     25X,'*****TAUOLA UNIVERSAL INTERFACE: ******',9X,1H*,
00348      $ /,' *',     25X,'*****VERSION 1.22, April 2009 (gfort)**',9X,1H*,
00349      $ /,' *',     25X,'**AUTHORS: P. Golonka, B. Kersevan, ***',9X,1H*,
00350      $ /,' *',     25X,'**T. Pierzchala, E. Richter-Was, ******',9X,1H*,
00351      $ /,' *',     25X,'****** Z. Was, M. Worek ***************',9X,1H*,
00352      $ /,' *',     25X,'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9X,1H*,
00353      $ /,' *',     25X,'*WITH C. Biscarat and S. Slabospitzky**',9X,1H*,
00354      $ /,' *',     25X,'****** are warmly acknowledged ********',9X,1H*,
00355      $ /,' *',     25X,'****** END OF MODULE OPERATION ********',9X,1H*,
00356      $  /,1X,15(5H*****)/)
00357 
00358       END
00359 
00360       SUBROUTINE SPINHIGGS(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
00361       LOGICAL IFPSEUDO
00362       REAL*8 HH1,HH2,wthiggs
00363       DIMENSION POL1(4), POL2(4),HH1(4),HH2(4), RRR(1)
00364 !             CALL DEXAY(1,POL1)  ! Kept for tests
00365 !             CALL DEXAY(2,POL2)  ! Kept for tests
00366       INTEGER ION(3)
00367  10   CONTINUE
00368          CALL RANMAR(RRR,1)
00369          CALL DEKAY(1,HH1)
00370          CALL DEKAY(2,HH2)
00371          wt=wthiggs(IFPSEUDO,HH1,HH2)
00372       IF (RRR(1).GT.WT) GOTO 10
00373       CALL DEKAY(1+10,HH1)
00374       CALL TAUPI0(0,1,ION)
00375       CALL DEKAY(2+10,HH2)
00376       CALL TAUPI0(0,2,ION)
00377       END
00378       FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
00379       LOGICAL IFPSEUDO
00380       common /pseudocoup/ csc,ssc
00381       REAL*4 csc,ssc
00382       save pseudocoup
00383       REAL*8 HH1(4),HH2(4),R(4,4),wthiggs
00384       DO K=1,4
00385        DO L=1,4
00386         R(K,L)=0
00387        ENDDO
00388       ENDDO
00389       WTHIGGS=0D0
00390       
00391       R(4,4)= 1D0    !  unpolarized part
00392       R(3,3)=-1D0    !  longitudinal
00393                      !  other missing
00394       IF (IFPSEUDO) THEN
00395         R(1,1)=-1
00396         R(2,2)= -1
00397         R(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
00398         R(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
00399         R(1,2)=2*csc*ssc/(csc**2+ssc**2)
00400         R(2,1)=-2*csc*ssc/(csc**2+ssc**2)
00401       ELSE
00402         R(1,1)=1
00403         R(2,2)=1
00404       ENDIF
00405 
00406 
00407 
00408       DO K=1,4
00409        DO L=1,4
00410         WTHIGGS=WTHIGGS+R(K,L)*HH1(K)*HH2(L)
00411        ENDDO
00412       ENDDO
00413         WTHIGGS=WTHIGGS/4D0
00414       END
00415 
00416       FUNCTION PLZAPX(HOPEin,IM0,NP1,NP2)
00417 C     IM0 NP1 NP2 are the positions of Z/gamma tau tau in hepevt common block.
00418 C     the purpose of this routine is to calculate polarization of Z along tau direction.
00419 C     this is highly non-trivial due to necessity of reading infromation from hard process
00420 C     history in HEPEVT, which is often written not up to the gramatic rules.
00421       REAL*8 PLZAP0,SVAR,COSTHE,sini,sfin,ZPROP2,
00422      $       P1(4),P2(4),Q1(4),Q2(4),QQ(4),PH(4),PD1(4),PD2(4),
00423      $       PQ1(4),PQ2(4),PB(4),PA(4)
00424       REAL*4 PLZAPX
00425       INTEGER IM
00426       LOGICAL HOPE,HOPEin
00427 C this is the hepevt class in old style. No d_h_ class pre-name
00428       INTEGER NMXHEP
00429       PARAMETER (NMXHEP=10000)
00430       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
00431       INTEGER nevhep,nhep,isthep,idhep,jmohep,
00432      $        jdahep
00433       COMMON /hepevt/
00434      $      nevhep,               ! serial number
00435      $      nhep,                 ! number of particles
00436      $      isthep(nmxhep),   ! status code
00437      $      idhep(nmxhep),    ! particle ident KF
00438      $      jmohep(2,nmxhep), ! parent particles
00439      $      jdahep(2,nmxhep), ! childreen particles
00440      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
00441      $      vhep(4,nmxhep)    ! vertex [mm]
00442 * ----------------------------------------------------------------------
00443       LOGICAL qedrad
00444       COMMON /phoqed/ 
00445      $     qedrad(nmxhep)    ! Photos flag
00446 * ----------------------------------------------------------------------
00447       SAVE hepevt,phoqed
00448 
00449 C(BPK)--> BROTHERS OF TAU ALREADY FOUND
00450       INTEGER ISON
00451       COMMON /ISONS_TAU/ISON(2)
00452 C(BPK)--<
00453 C >>
00454 C >> STEP 1: find where are particles in hepevent and pick them up
00455 C >>
00456              HOPE=HOPEin
00457 C sometimes shade Z of Z is its mother ...
00458             IM=IM0
00459             IM00=JMOHEP(1,IM0)
00460 C to protect against check on mother of beam particles.
00461             IF (IM00.GT.0) THEN
00462               IF (IDHEP(IM0).EQ.IDHEP(IM00)) IM=JMOHEP(1,IM0)
00463             ENDIF                                                              
00464 C
00465 C find (host generator-level) incoming beam-bare-particles which form Z and co.
00466             IMO1=JMOHEP(1,IM)
00467             IMO2=JMOHEP(2,IM)
00468 
00469 C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS
00470             IM00=IMO1
00471             IF (ISTHEP(IM00).EQ.120) THEN
00472                IMO1=JMOHEP(1,IM00)
00473                IMO2=JMOHEP(2,IM00)
00474             ENDIF
00475 C(BPK)--<
00476 
00477             IFFULL=0
00478 C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt.
00479       IF (IMO1.EQ.0.AND.IMO2.EQ.0) THEN
00480          IMO1=JMOHEP(1,NP1)
00481 C(BPK)-->
00482          IF (IDHEP(IMO1).EQ.IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
00483 C(BPK)--<
00484          IMO2=JMOHEP(2,NP1)
00485 C(BPK)-->
00486          IF (IDHEP(IMO2).EQ.IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
00487 C(BPK)--<
00488          IFFULL=1
00489 C case when it was like  qq --> tau+tau- gammas and qq were NOT 'first' in hepevt.
00490 
00491       ELSEIF (IDHEP(IM).NE.22.AND.IDHEP(IM).NE.23) THEN
00492          IMO1=JMOHEP(1,NP1)
00493 C(BPK)-->
00494          IF (IDHEP(IMO1).EQ.IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
00495 C(BPK)--<
00496          IMO2=JMOHEP(2,NP1)
00497 C(BPK)-->
00498          IF (IDHEP(IMO2).EQ.IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
00499 C(BPK)--<
00500          IFFULL=1
00501       ENDIF
00502 
00503             
00504 C and check if it really happened
00505             IF (IMO1.EQ.0) HOPE=.FALSE.
00506             IF (IMO2.EQ.0) HOPE=.FALSE.
00507             IF (IMO1.EQ.IMO2) HOPE=.FALSE.
00508 
00509 C 
00510          DO I=1,4
00511             Q1(I)= PHEP(I,NP1)              !momentum of tau+
00512             Q2(I)= PHEP(I,NP2)              !momentum of tau-
00513          ENDDO
00514 
00515 C corrections due to possible differences in 4-momentum of shadow vs true Z.
00516 C(BPK)-->
00517       IF (IM.EQ.JMOHEP(1,IM0).AND.
00518      $     (IDHEP(IM).EQ.22.OR.IDHEP(IM).EQ.23)) THEN
00519          DO K=1,4
00520             PB(K)=PHEP(K,IM)
00521             PA(K)=PHEP(K,IM0)
00522          ENDDO
00523 C(BPK)--<
00524 
00525                CALL BOSTDQ( 1,PA, Q1, Q1)
00526                CALL BOSTDQ( 1,PA, Q2, Q2)
00527                CALL BOSTDQ(-1,PB, Q1, Q1)
00528                CALL BOSTDQ(-1,PB, Q2, Q2)
00529 
00530               ENDIF
00531 
00532          DO I=1,4                                                                                   
00533             QQ(I)= Q1(I)+Q2(I)              !momentum of Z
00534             IF (HOPE) P1(I)=PHEP(I,IMO1)    !momentum of beam1
00535             IF (HOPE) P2(I)=PHEP(I,IMO2)    !momentum of beam2
00536             PH(I)=0D0
00537             PD1(I)=0D0
00538             PD2(I)=0D0
00539          ENDDO
00540 !        These momenta correspond to  quarks, gluons photons or taus
00541                    IDFQ1=IDHEP(NP1)
00542                    IDFQ2=IDHEP(NP2)
00543          IF (HOPE) IDFP1=IDHEP(IMO1)
00544          IF (HOPE) IDFP2=IDHEP(IMO2)
00545 
00546          SVAR=QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2
00547          IF (.NOT.HOPE) THEN
00548 C options which may be useful in some cases of two heavy boson production
00549 C need individual considerations. To be developed.
00550 C           PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0)  ! gamma/Z mixture as if produced from e beam
00551 C           PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0)  ! pure Z
00552            PLZAPX=0.5                         ! pure gamma
00553            RETURN
00554          ENDIF
00555 C >>
00556 C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
00557 C >>
00558 C let us define beginning and end of particles which are produced in parallel to Z
00559 C in principle following should work
00560 
00561 C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
00562          NX1=JDAHEP(1,IM00)
00563          NX2=JDAHEP(2,IM00)
00564 C but ...
00565          INBR=IM ! OK, HARD RECORD Z/GAMMA 
00566          IF (IFFULL.EQ.1) INBR=NP1  ! OK, NO Z/GAMMA
00567          IF (IDHEP(JMOHEP(1,INBR)).EQ.IDHEP(INBR)) INBR=JMOHEP(1,INBR) ! FORCE HARD RECORD
00568 C(BPK)--<
00569          IF(NX1.EQ.0.OR.NX2.EQ.0) THEN
00570            NX1=INBR
00571            NX2=INBR
00572            DO K=1,INBR-1
00573              IF(JMOHEP(1,INBR-K).EQ.JMOHEP(1,INBR)) THEN
00574               NX1=INBR-K
00575              ELSE
00576                 GOTO 7
00577              ENDIF
00578            ENDDO
00579  7         CONTINUE
00580 
00581            DO K=INBR+1,NHEP
00582              IF(JMOHEP(1,K).EQ.JMOHEP(1,INBR)) THEN
00583               NX2=K
00584              ELSE
00585                 GOTO 8
00586              ENDIF
00587            ENDDO
00588  8         CONTINUE
00589          ENDIF
00590 
00591 C case of annihilation of two bosons is hopeles
00592          IF (ABS(IDFP1).GE.20.AND.ABS(IDFP2).GE.20) HOPE=.FALSE.
00593 C case of annihilation of non-matching flavors is hopeless
00594          IF (ABS(IDFP1).LE.20.AND.ABS(IDFP2).LE.20.AND.IDFP1+IDFP2.NE.0)
00595      $       HOPE=.FALSE.
00596          IF (.NOT.HOPE) THEN
00597 C options which may be useful in some cases of two heavy boson production
00598 C need individual considerations. To be developed.
00599 C           PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0)  ! gamma/Z mixture as if produced from e beam
00600 C           PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0)  ! pure Z
00601            PLZAPX=0.5                         ! pure gamma
00602            RETURN
00603          ENDIF
00604          IF (ABS(IDFP1).LT.20) IDE= IDFP1
00605          IF (ABS(IDFP2).LT.20) IDE=-IDFP2
00606 
00607 
00608 C >>
00609 C >> STEP 3  we combine gluons, photons  into incoming effective beams
00610 C >>
00611 
00612 C in the following we will ignore the possibility of photon emission from taus
00613 C however at certain moment it will be necessary to take care of
00614 
00615            DO L=1,4
00616             PD1(L)=P1(L)
00617             PD2(L)=P2(L)
00618            ENDDO 
00619 
00620                 DO L=1,4
00621                 PQ1(L)=Q1(L)
00622                 PQ2(L)=Q2(L)
00623                 ENDDO 
00624 
00625          IFLAV=min(ABS(IDFP1),ABS(IDFP2))
00626 
00627 *--------------------------------------------------------------------------
00628 c                  IFLAV=min(ABS(IDFP1),ABS(IDFP2))
00629 c that means that always quark or lepton i.e. process like
00630 c      f g(gamma) --> f Z0 --> tau tau
00631 c we glue  fermions to effective  beams that is f f --> Z0 --> tau tau
00632 c  with  gamma/g emission from initial fermion.
00633 *---------------------------------------------------------------------------
00634  
00635          IF (ABS(IDFP1).GE.20) THEN
00636            DO k=NX1,NX2
00637              IDP=IDHEP(k)
00638              IF (ABS(IDP).EQ.IFLAV) THEN
00639                DO L=1,4
00640                  PD1(L)=-PHEP(L,K)
00641                ENDDO
00642              ENDIF
00643            ENDDO
00644          ENDIF
00645 
00646          IF (ABS(IDFP2).GE.20) THEN
00647            DO k=NX1,NX2
00648              IDP=IDHEP(k)
00649              IF (ABS(IDP).EQ.IFLAV) THEN
00650                DO L=1,4
00651                  PD2(L)=-PHEP(L,K)
00652                ENDDO
00653              ENDIF
00654            ENDDO
00655          ENDIF
00656 
00657 C if first beam was boson: gluing
00658 
00659          IF (ABS(IDFP1).GE.20) THEN
00660            DO L=1,4
00661              PH(L)=P1(L)
00662            ENDDO
00663            xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
00664      $            -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
00665            xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
00666      $            -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
00667           IF (XM1.LT.XM2) THEN
00668              DO L=1,4
00669                PD1(L)=PD1(L)+P1(L)
00670              ENDDO
00671            ELSE
00672              DO L=1,4
00673                PD2(L)=PD2(L)+P1(L)
00674              ENDDO
00675            ENDIF
00676          ENDIF
00677 
00678 C if second beam was boson: gluing 
00679 
00680 
00681          IF (ABS(IDFP2).GE.20) THEN
00682            DO L=1,4
00683              PH(L)=P2(L)
00684            ENDDO
00685            xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
00686      $            -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
00687            xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
00688      $            -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
00689            IF (XM1.LT.XM2) THEN
00690              DO L=1,4
00691                PD1(L)=PD1(L)+P2(L)
00692              ENDDO
00693            ELSE
00694              DO L=1,4
00695                PD2(L)=PD2(L)+P2(L)
00696              ENDDO
00697            ENDIF
00698          ENDIF
00699 
00700 C now spectators ...
00701 
00702 C(BPK)-->
00703       NPH1=NP1
00704       NPH2=NP2
00705       IF (IDHEP(JMOHEP(1,NP1)).EQ.IDHEP(NP1)) NPH1=JMOHEP(1,NP1) ! SHOULD PUT US IN HARD REC.
00706       IF (IDHEP(JMOHEP(1,NP2)).EQ.IDHEP(NP2)) NPH2=JMOHEP(1,NP2) ! SHOULD PUT US IN HARD REC.
00707 C(BPK)--<
00708 
00709          DO k=NX1,NX2
00710          IF (ABS(IDHEP(K)).NE.IFLAV.AND.K.NE.IM.AND.
00711 C(BPK)-->
00712      $       K.NE.NPH1.AND.K.NE.NPH2) THEN 
00713 C(BPK)--<
00714           IF(IDHEP(K).EQ.22.AND.IFFULL.EQ.1) THEN
00715             DO L=1,4
00716               PH(L)=PHEP(L,K)
00717             ENDDO
00718             xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
00719      $             -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
00720             xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
00721      $             -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
00722            xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
00723      $            -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
00724            xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
00725      $            -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)
00726 
00727   
00728             sini=abs((PD1(4)+PD2(4)-PH(4))**2-(PD1(3)+PD2(3)-PH(3))**2
00729      $              -(PD1(2)+PD2(2)-PH(2))**2-(PD1(1)+PD2(1)-PH(1))**2)
00730             sfin=abs((PD1(4)+PD2(4)      )**2-(PD1(3)+PD2(3)      )**2
00731      $              -(PD1(2)+PD2(2)      )**2-(PD1(1)+PD2(1)      )**2)
00732 
00733            FACINI=ZPROP2(sini)
00734            FACFIN=ZPROP2(sfin)
00735 
00736            XM1=XM1/FACINI
00737            XM2=XM2/FACINI
00738            XM3=XM3/FACFIN
00739            XM4=XM4/FACFIN
00740 
00741            XM=MIN(XM1,XM2,XM3,XM4)
00742                   IF      (XM1.EQ.XM) THEN 
00743                      DO L=1,4
00744                        PD1(L)=PD1(L)-PH(L)
00745                      ENDDO
00746                   ELSEIF   (XM2.EQ.XM) THEN 
00747                      DO L=1,4
00748                        PD2(L)=PD2(L)-PH(L)
00749                      ENDDO
00750                   ELSEIF   (XM3.EQ.XM) THEN 
00751                      DO L=1,4
00752                         Q1(L)=PQ1(L)+PH(L)
00753                      ENDDO
00754                   ELSE
00755                      DO L=1,4
00756                         Q2(L)=PQ2(L)+PH(L)
00757                      ENDDO
00758                   ENDIF
00759            ELSE
00760             DO L=1,4
00761               PH(L)=PHEP(L,K)
00762             ENDDO
00763             xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
00764      $             -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
00765             xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
00766      $             -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
00767             IF (XM1.LT.XM2) THEN
00768               DO L=1,4
00769                 PD1(L)=PD1(L)-PH(L)
00770               ENDDO
00771             ELSE
00772               DO L=1,4
00773                 PD2(L)=PD2(L)-PH(L)
00774               ENDDO
00775             ENDIF
00776            ENDIF
00777           ENDIF
00778          ENDDO
00779 
00780 
00781 C >>
00782 C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in 
00783 c >>          effective outcoming taus
00784 C >>
00785 C let us define beginning and end of particles which are produced in 
00786 c  parallel to tau
00787 
00788 
00789 
00790 C find outcoming particles which come from Z   
00791 
00792        
00793  
00794 
00795 C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER       
00796       IF (ABS(IDHEP(IM0)).EQ.22.OR.abs(IDHEP(IM0)).EQ.23) THEN
00797          DO K=ISON(1),ISON(2)
00798             IF(ABS(IDHEP(K)).EQ.22) THEN
00799 C(BPK)--< 
00800 
00801               do l=1,4
00802               ph(l)=phep(l,k)
00803               enddo
00804 
00805            xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
00806      $            -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
00807            xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
00808      $            -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)  
00809 
00810            XM=MIN(XM3,XM4) 
00811 
00812                   IF   (XM3.EQ.XM) THEN 
00813                      DO L=1,4
00814                         Q1(L)=PQ1(L)+PH(L)
00815                      ENDDO
00816                   ELSE
00817                      DO L=1,4
00818                         Q2(L)=PQ2(L)+PH(L)
00819                      ENDDO
00820                   ENDIF
00821             endif
00822           enddo
00823           ENDIF
00824 
00825 
00826 
00827 *------------------------------------------------------------------------
00828 
00829 
00830 C out of effective momenta we calculate COSTHE and later polarization
00831       CALL ANGULU(PD1,PD2,Q1,Q2,COSTHE)
00832      
00833       PLZAPX=PLZAP0(IDE,IDFQ1,SVAR,COSTHE)
00834       END
00835 
00836       SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE)
00837       REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
00838 C take effective beam which is less massive, it should be irrelevant
00839 C but in case HEPEVT is particulary dirty may help.
00840 C this routine calculate reduced system transver and cosine of scattering 
00841 C angle.
00842 
00843       XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2)
00844       XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2)
00845       IF (XM1.LT.XM2) THEN
00846         SIGN=1D0
00847         DO K=1,4
00848           P(K)=PD1(K)
00849         ENDDO
00850       ELSE
00851         SIGN=-1D0
00852         DO K=1,4
00853           P(K)=PD2(K)
00854         ENDDO
00855       ENDIF
00856 C calculate space like part of P (in Z restframe)
00857       DO K=1,4
00858        QQ(K)=Q1(k)+Q2(K)
00859        QT(K)=Q1(K)-Q2(K)
00860       ENDDO
00861 
00862        XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2)
00863 
00864        QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1)
00865       DO K=1,4
00866        QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2
00867       ENDDO
00868 
00869        PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1)
00870       DO K=1,4
00871        P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2
00872       ENDDO
00873 C calculate costhe
00874        PXP  =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
00875        QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2)
00876        PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4)
00877        COSTHE=PXQT/PXP/QTXQT
00878        COSTHE=COSTHE*SIGN
00879       END
00880 
00881       FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0)
00882 C this function calculates probability for the helicity +1 +1 configuration
00883 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
00884       REAL*8 PLZAP0,SVAR,COSTHE,COSTH0,T_BORN
00885 
00886       COSTHE=COSTH0
00887 C >>>>>      IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
00888 C >>>>>      of first beam is used by T_GIVIZ0 including sign
00889 
00890       IF (IDF.GT.0) THEN
00891         CALL INITWK(IDE,IDF,SVAR)
00892       ELSE
00893         CALL INITWK(-IDE,-IDF,SVAR)
00894       ENDIF
00895       PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0)
00896      $  /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0))
00897 
00898 !      PLZAP0=0.5
00899       END
00900       FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB)
00901 C ----------------------------------------------------------------------
00902 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME         
00903 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER       
00904 C EXAMPLE OF THE METHOD APPLIED THERE                               
00905 C INPUT PARAMETERS ARE: SVAR    -- transfer
00906 C                       COSTHE  -- cosine of angle between tau+ and 1st beam
00907 C                       TA,TB   -- helicity states of tau+ tau-
00908 C
00909 C     called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
00910 C ----------------------------------------------------------------------
00911       IMPLICIT REAL*8(A-H,O-Z)
00912       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
00913       REAL*8              ENE ,AMIN,AMFIN
00914       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
00915      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
00916      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
00917      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
00918       REAL*8             SS,POLN,T3E,QE,T3F,QF
00919      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
00920       REAL*8            SEPS1,SEPS2
00921 C=====================================================================
00922       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00923       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00924 C     SWSQ        = sin2 (theta Weinberg)
00925 C     AMW,AMZ     = W & Z boson masses respectively
00926 C     AMH         = the Higgs mass
00927 C     AMTOP       = the top mass
00928 C     GAMMZ       = Z0 width
00929       COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
00930       COMPLEX*16 XUPZFP(2),XUPZIP(2)
00931       COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
00932       COMPLEX*16 PROPA,PROPZ
00933       COMPLEX*16 XR,XI
00934       COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
00935       COMPLEX*16 XTHING,XVE,XVF,XVEF
00936       DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/
00937       DATA MODE0 /-5/
00938       DATA IDE0 /-55/
00939       DATA SVAR0,COST0 /-5.D0,-6.D0/
00940       DATA PI /3.141592653589793238462643D0/
00941       DATA SEPS1,SEPS2 /0D0,0D0/
00942 C
00943 C MEMORIZATION =========================================================
00944       IF ( MODE.NE.MODE0.OR.SVAR.NE.SVAR0.OR.COSTHE.NE.COST0
00945      $    .OR.IDE0.NE.IDE)THEN
00946 C
00947         KEYGSW=1
00948 C ** PROPAGATORS
00949         IDE0=IDE
00950         MODE0=MODE
00951         SVAR0=SVAR
00952         COST0=COSTHE
00953         SINTHE=SQRT(1.D0-COSTHE**2)
00954         BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR))
00955 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
00956         XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2))
00957         XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2))
00958         XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2))
00959         XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2))
00960 C FINAL STATE VECTOR COUPLING
00961         XUPF     =0.5D0*(XUPZF(1)+XUPZF(2))
00962         XUPI     =0.5D0*(XUPZI(1)+XUPZI(2))
00963         XTHING   =0D0
00964 
00965         PROPA =1D0/SVAR
00966         PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ)
00967         IF (KEYGSW.EQ.0) PROPZ=0.D0
00968         DO 50 I=1,2
00969          DO 50 J=1,2
00970           REGULA= (3-2*I)*(3-2*J) + COSTHE
00971           REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR)
00972           APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA)
00973           AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA
00974           ABORN(I,J)=APHOT(I,J)+AZETT(I,J)
00975           APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM
00976           AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM
00977           ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J)
00978    50   CONTINUE
00979       ENDIF
00980 C
00981 C******************
00982 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
00983 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
00984 C* HELICITY CONSERVATION EXPLICITLY OBEYED
00985       POLAR1=  (SEPS1)
00986       POLAR2= (-SEPS2)
00987       BORN=0D0
00988       DO 150 I=1,2
00989        HELIC= 3-2*I
00990        DO 150 J=1,2
00991         HELIT=3-2*J
00992         FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0
00993         FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB)
00994         FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB)
00995 
00996         BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR
00997 C      MASS TERM IN BORN
00998         IF (MODE.GE.1) THEN
00999          BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
01000         ENDIF
01001 
01002   150 CONTINUE
01003 C************
01004       FUNT=BORN
01005       IF(FUNT.LT.0.D0)  FUNT=BORN
01006 
01007 C
01008       IF (SVAR.GT.4D0*AMFIN**2) THEN
01009 C PHASE SPACE THRESHOLD FACTOR
01010         THRESH=SQRT(1-4D0*AMFIN**2/SVAR)
01011         T_BORN= FUNT*SVAR**2*THRESH
01012       ELSE
01013         THRESH=0.D0
01014         T_BORN=0.D0
01015       ENDIF
01016 C ZW HERE WAS AN ERROR 19. 05. 1989
01017 !      write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
01018 !      write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
01019       END
01020 
01021       SUBROUTINE INITWK(IDEX,IDFX,SVAR)
01022 ! initialization routine coupling masses etc.
01023       IMPLICIT REAL*8 (A-H,O-Z)
01024       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
01025       REAL*8              ENE ,AMIN,AMFIN
01026       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
01027      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
01028      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
01029      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
01030       REAL*8             SS,POLN,T3E,QE,T3F,QF
01031      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
01032       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
01033       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
01034 C     SWSQ        = sin2 (theta Weinberg)
01035 C     AMW,AMZ     = W & Z boson masses respectively
01036 C     AMH         = the Higgs mass
01037 C     AMTOP       = the top mass
01038 C     GAMMZ       = Z0 width
01039 C
01040       ENE=SQRT(SVAR)/2
01041       AMIN=0.511D-3
01042       SWSQ=0.23147
01043       AMZ=91.1882
01044       GAMMZ=2.4952
01045       IF     (IDFX.EQ. 15) then       
01046         IDF=2  ! denotes tau +2 tau-
01047         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
01048       ELSEIF (IDFX.EQ.-15) then
01049         IDF=-2  ! denotes tau -2 tau-
01050         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
01051       ELSE
01052         WRITE(*,*) 'INITWK: WRONG IDFX'
01053         STOP
01054       ENDIF
01055 
01056       IF     (IDEX.EQ. 11) then      !electron
01057         IDE= 2
01058         AMIN=0.511D-3
01059       ELSEIF (IDEX.EQ.-11) then      !positron
01060         IDE=-2
01061         AMIN=0.511D-3
01062       ELSEIF (IDEX.EQ. 13) then      !mu+
01063         IDE= 2
01064         AMIN=0.105659
01065       ELSEIF (IDEX.EQ.-13) then      !mu-
01066         IDE=-2
01067         AMIN=0.105659
01068       ELSEIF (IDEX.EQ.  1) then      !d
01069         IDE= 4
01070         AMIN=0.05
01071       ELSEIF (IDEX.EQ.- 1) then      !d~
01072         IDE=-4
01073         AMIN=0.05
01074       ELSEIF (IDEX.EQ.  2) then      !u
01075         IDE= 3
01076         AMIN=0.02
01077       ELSEIF (IDEX.EQ.- 2) then      !u~
01078         IDE=-3
01079         AMIN=0.02
01080       ELSEIF (IDEX.EQ.  3) then      !s
01081         IDE= 4
01082         AMIN=0.3
01083       ELSEIF (IDEX.EQ.- 3) then      !s~
01084         IDE=-4
01085         AMIN=0.3
01086       ELSEIF (IDEX.EQ.  4) then      !c
01087         IDE= 3
01088         AMIN=1.3
01089       ELSEIF (IDEX.EQ.- 4) then      !c~
01090         IDE=-3
01091         AMIN=1.3
01092       ELSEIF (IDEX.EQ.  5) then      !b
01093         IDE= 4
01094         AMIN=4.5
01095       ELSEIF (IDEX.EQ.- 5) then      !b~
01096         IDE=-4
01097         AMIN=4.5
01098       ELSEIF (IDEX.EQ.  12) then     !nu_e
01099         IDE= 1
01100         AMIN=0.1D-3
01101       ELSEIF (IDEX.EQ.- 12) then     !nu_e~
01102         IDE=-1
01103         AMIN=0.1D-3
01104       ELSEIF (IDEX.EQ.  14) then     !nu_mu
01105         IDE= 1
01106         AMIN=0.1D-3
01107       ELSEIF (IDEX.EQ.- 14) then     !nu_mu~
01108         IDE=-1
01109         AMIN=0.1D-3
01110       ELSEIF (IDEX.EQ.  16) then     !nu_tau
01111         IDE= 1
01112         AMIN=0.1D-3
01113       ELSEIF (IDEX.EQ.- 16) then     !nu_tau~
01114         IDE=-1
01115         AMIN=0.1D-3
01116 
01117       ELSE
01118         WRITE(*,*) 'INITWK: WRONG IDEX'
01119         STOP
01120       ENDIF
01121 
01122 C ----------------------------------------------------------------------
01123 C
01124 C     INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
01125 C
01126 C     called by : KORALZ
01127 C ----------------------------------------------------------------------
01128       ITCE=IDE/IABS(IDE)
01129       JTCE=(1-ITCE)/2
01130       ITCF=IDF/IABS(IDF)
01131       JTCF=(1-ITCF)/2
01132       CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM)
01133       CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM)
01134       XUPGI(1)=QE
01135       XUPGI(2)=QE
01136       T3E    = AIZOL+AIZOR
01137       XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01138       XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01139       CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR)
01140       CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR)
01141       XUPGF(1)=QF
01142       XUPGF(2)=QF
01143       T3F    =  AIZOL+AIZOR
01144       XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01145       XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01146 C
01147       NDIAG0=2
01148       NDIAGA=11
01149       KEYA  = 1
01150       KEYZ  = 1
01151 C
01152 C
01153       RETURN
01154       END
01155 
01156       SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
01157 C ----------------------------------------------------------------------
01158 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
01159 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
01160 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
01161 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
01162 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
01163 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
01164 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
01165 C
01166 C     called by : EVENTE, EVENTM, FUNTIH, .....
01167 C ----------------------------------------------------------------------
01168       IMPLICIT REAL*8(A-H,O-Z)
01169 C
01170       IF(IDFERM.EQ.0.OR.IABS(IDFERM).GT.4) GOTO 901
01171       IF(IABS(IHELIC).NE.1)                GOTO 901
01172       IH  =IHELIC
01173       IDTYPE =IABS(IDFERM)
01174       IC  =IDFERM/IDTYPE
01175       LEPQUA=INT(IDTYPE*0.4999999D0)
01176       IUPDOW=IDTYPE-2*LEPQUA-1
01177       CHARGE  =(-IUPDOW+2D0/3D0*LEPQUA)*IC
01178       SIZO3   =0.25D0*(IC-IH)*(1-2*IUPDOW)
01179       KOLOR=1+2*LEPQUA
01180 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
01181 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01182       RETURN
01183  901  PRINT *,' STOP IN GIVIZO: WRONG PARAMS.'
01184       STOP
01185       END
01186       SUBROUTINE PHYFIX(NSTOP,NSTART)
01187       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) 
01188       SAVE /LUJETS/ 
01189 C NSTOP NSTART : when PHYTIA history ends and event starts.
01190       NSTOP=0
01191       NSTART=1
01192       DO I=1, N
01193        IF(K(I,1).NE.21) THEN
01194            NSTOP = I-1
01195            NSTART= I
01196            GOTO 500
01197        ENDIF
01198       ENDDO
01199  500  CONTINUE
01200       END
01201       SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
01202 C ----------------------------------------------------------------------
01203 C this subroutine fills one entry into the HEPEVT common
01204 C and updates the information for affected mother entries
01205 C
01206 C written by Martin W. Gruenewald (91/01/28)
01207 C
01208 C     called by : ZTOHEP,BTOHEP,DWLUxy
01209 C ----------------------------------------------------------------------
01210 C
01211 C this is the hepevt class in old style. No d_h_ class pre-name
01212 C this is the hepevt class in old style. No d_h_ class pre-name
01213       INTEGER NMXHEP
01214       PARAMETER (NMXHEP=10000)
01215       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01216       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01217      $        jdahep
01218       COMMON /hepevt/
01219      $      nevhep,               ! serial number
01220      $      nhep,                 ! number of particles
01221      $      isthep(nmxhep),   ! status code
01222      $      idhep(nmxhep),    ! particle ident KF
01223      $      jmohep(2,nmxhep), ! parent particles
01224      $      jdahep(2,nmxhep), ! childreen particles
01225      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01226      $      vhep(4,nmxhep)    ! vertex [mm]
01227 * ----------------------------------------------------------------------
01228       LOGICAL qedrad
01229       COMMON /phoqed/ 
01230      $     qedrad(nmxhep)    ! Photos flag
01231 * ----------------------------------------------------------------------
01232       SAVE hepevt,phoqed
01233       LOGICAL PHFLAG
01234 C
01235       REAL*4  P4(4)
01236 C
01237 C check address mode
01238       IF (N.EQ.0) THEN
01239 C
01240 C append mode
01241         IHEP=NHEP+1
01242       ELSE IF (N.GT.0) THEN
01243 C
01244 C absolute position
01245         IHEP=N
01246       ELSE
01247 C
01248 C relative position
01249         IHEP=NHEP+N
01250       END IF
01251 C
01252 C check on IHEP
01253       IF ((IHEP.LE.0).OR.(IHEP.GT.NMXHEP)) RETURN
01254 C
01255 C add entry
01256       NHEP=IHEP
01257       ISTHEP(IHEP)=IST
01258       IDHEP(IHEP)=ID
01259       JMOHEP(1,IHEP)=JMO1
01260       IF(JMO1.LT.0)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP
01261       JMOHEP(2,IHEP)=JMO2
01262       IF(JMO2.LT.0)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP
01263       JDAHEP(1,IHEP)=JDA1
01264       JDAHEP(2,IHEP)=JDA2
01265 C
01266       DO I=1,4
01267         PHEP(I,IHEP)=P4(I)
01268 C
01269 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
01270         VHEP(I,IHEP)=0.0
01271       END DO
01272       PHEP(5,IHEP)=PINV
01273 C FLAG FOR PHOTOS...
01274       QEDRAD(IHEP)=PHFLAG
01275 C
01276 C update process:
01277       DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
01278         IF(IP.GT.0)THEN
01279 C
01280 C if there is a daughter at IHEP, mother entry at IP has decayed
01281           IF(ISTHEP(IP).EQ.1)ISTHEP(IP)=2
01282 C
01283 C and daughter pointers of mother entry must be updated
01284           IF(JDAHEP(1,IP).EQ.0)THEN
01285             JDAHEP(1,IP)=IHEP
01286             JDAHEP(2,IP)=IHEP
01287           ELSE
01288             JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP))
01289           END IF
01290         END IF
01291       END DO
01292 C
01293       RETURN
01294       END
01295 
01296 
01297       FUNCTION IHEPDIM(DUM) 
01298 C this is the hepevt class in old style. No d_h_ class pre-name
01299 C this is the hepevt class in old style. No d_h_ class pre-name
01300       INTEGER NMXHEP
01301       PARAMETER (NMXHEP=10000)
01302       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01303       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01304      $        jdahep
01305       COMMON /hepevt/
01306      $      nevhep,               ! serial number
01307      $      nhep,                 ! number of particles
01308      $      isthep(nmxhep),   ! status code
01309      $      idhep(nmxhep),    ! particle ident KF
01310      $      jmohep(2,nmxhep), ! parent particles
01311      $      jdahep(2,nmxhep), ! childreen particles
01312      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01313      $      vhep(4,nmxhep)    ! vertex [mm]
01314 * ----------------------------------------------------------------------
01315       LOGICAL qedrad
01316       COMMON /phoqed/ 
01317      $     qedrad(nmxhep)    ! Photos flag
01318 * ----------------------------------------------------------------------
01319       SAVE hepevt,phoqed
01320       IHEPDIM=NHEP
01321       END
01322       FUNCTION ZPROP2(S)
01323       IMPLICIT REAL*8(A-H,O-Z)
01324       COMPLEX*16 CPRZ0,CPRZ0M
01325       AMZ=91.1882
01326       GAMMZ=2.49
01327       CPRZ0=DCMPLX((S-AMZ**2),S/AMZ*GAMMZ)
01328       CPRZ0M=1/CPRZ0
01329       ZPROP2=(ABS(CPRZ0M))**2
01330       END
01331 
01332       SUBROUTINE TAUPI0(MODE,JAK,ION)
01333 C no initialization required. Must be called once after every:
01334 C   1)    CALL DEKAY(1+10,...)
01335 C   2)    CALL DEKAY(2+10,...)
01336 C   3)    CALL DEXAY(1,...)
01337 C   4)    CALL DEXAY(2,...)
01338 C subroutine to decay originating from TAUOLA's taus: 
01339 C 1) etas (with CALL TAUETA(JAK))
01340 C 2) later pi0's from taus.
01341 C 3) extensions to other applications possible. 
01342 C this routine belongs to >tauola universal interface<, but uses 
01343 C routines from >tauola< utilities as well.  25.08.2005      
01344 C this is the hepevt class in old style. No d_h_ class pre-name
01345       INTEGER NMXHEP
01346       PARAMETER (NMXHEP=10000)
01347       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01348       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01349      $        jdahep
01350       COMMON /hepevt/
01351      $      nevhep,               ! serial number
01352      $      nhep,                 ! number of particles
01353      $      isthep(nmxhep),   ! status code
01354      $      idhep(nmxhep),    ! particle ident KF
01355      $      jmohep(2,nmxhep), ! parent particles
01356      $      jdahep(2,nmxhep), ! childreen particles
01357      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01358      $      vhep(4,nmxhep)    ! vertex [mm]
01359 * ----------------------------------------------------------------------
01360       LOGICAL qedrad
01361       COMMON /phoqed/ 
01362      $     qedrad(nmxhep)    ! Photos flag
01363 * ----------------------------------------------------------------------
01364       SAVE hepevt,phoqed
01365 
01366 C position of taus, must be defined by host program:
01367       COMMON /TAUPOS/ NP1,NP2
01368 c
01369       REAL  PHOT1(4),PHOT2(4)
01370       REAL*8  R,X(4),Y(4),PI0(4)
01371       INTEGER JEZELI(3),ION(3)
01372       DATA JEZELI /0,0,0/
01373       SAVE JEZELI
01374       IF (MODE.EQ.-1) THEN
01375         JEZELI(1)=ION(1)
01376         JEZELI(2)=ION(2)
01377         JEZELI(3)=ION(3)
01378         RETURN
01379       ENDIF
01380       IF (JEZELI(1).EQ.0) RETURN
01381       IF (JEZELI(2).EQ.1) CALL TAUETA(JAK)
01382       IF (JEZELI(3).EQ.1) CALL TAUK0S(JAK)
01383 C position of decaying particle:
01384       IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
01385         NPS=NP1
01386       ELSE
01387         NPS=NP2
01388       ENDIF
01389       nhepM=nhep                ! to avoid infinite loop
01390       DO K=JDAHEP(1,NPS),nhepM  ! we search for pi0's from tau till eor.
01391        IF (IDHEP(K).EQ.111.AND.JDAHEP(1,K).LE.K) THEN ! IF we found pi0
01392         DO L=1,4
01393           PI0(L)= phep(L,K)
01394         ENDDO
01395 ! random 3 vector on the sphere, masless
01396         R=SQRT(PI0(4)**2-PI0(3)**2-PI0(2)**2-PI0(1)**2)/2D0
01397         CALL SPHERD(R,X)
01398         X(4)=R
01399         Y(4)=R
01400         
01401         Y(1)=-X(1)
01402         Y(2)=-X(2)
01403         Y(3)=-X(3)
01404 ! boost to lab and to real*4
01405         CALL bostdq(-1,PI0,X,X)
01406         CALL bostdq(-1,PI0,Y,Y)
01407         DO L=1,4
01408          PHOT1(L)=X(L)
01409          PHOT2(L)=Y(L)
01410         ENDDO
01411 C to hepevt
01412         CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
01413         CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
01414        ENDIF
01415       ENDDO
01416 C
01417       END
01418       SUBROUTINE TAUETA(JAK)
01419 C subroutine to decay etas's from taus. 
01420 C this routine belongs to tauola universal interface, but uses 
01421 C routines from tauola utilities. Just flat phase space, but 4 channels.
01422 C it is called at the beginning of SUBR. TAUPI0(JAK)
01423 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005    
01424 C this is the hepevt class in old style. No d_h_ class pre-name
01425       INTEGER NMXHEP
01426       PARAMETER (NMXHEP=10000)
01427       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01428       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01429      $        jdahep
01430       COMMON /hepevt/
01431      $      nevhep,               ! serial number
01432      $      nhep,                 ! number of particles
01433      $      isthep(nmxhep),   ! status code
01434      $      idhep(nmxhep),    ! particle ident KF
01435      $      jmohep(2,nmxhep), ! parent particles
01436      $      jdahep(2,nmxhep), ! childreen particles
01437      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01438      $      vhep(4,nmxhep)    ! vertex [mm]
01439 * ----------------------------------------------------------------------
01440       LOGICAL qedrad
01441       COMMON /phoqed/ 
01442      $     qedrad(nmxhep)    ! Photos flag
01443 * ----------------------------------------------------------------------
01444       SAVE hepevt,phoqed
01445       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01446      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01447      *                 ,AMK,AMKZ,AMKST,GAMKST
01448 *
01449       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01450      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01451      *                 ,AMK,AMKZ,AMKST,GAMKST
01452 
01453 C position of taus, must be defined by host program:
01454       COMMON /TAUPOS/ NP1,NP2
01455 c
01456       REAL  RRR(1),BRSUM(3), RR(2)
01457       REAL  PHOT1(4),PHOT2(4),PHOT3(4)
01458       REAL*8    X(4),    Y(4),    Z(4)
01459       REAL                                YM1,YM2,YM3
01460       REAL*8  R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
01461       REAL*8 a,b,c
01462       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
01463 C position of decaying particle:
01464       IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
01465         NPS=NP1
01466       ELSE
01467         NPS=NP2
01468       ENDIF
01469       nhepM=nhep                ! to avoid infinite loop              
01470       DO K=JDAHEP(1,NPS),nhepM  ! we search for etas's from tau till eor.
01471        IF (IDHEP(K).EQ.221.AND.JDAHEP(1,K).LE.K) THEN ! IF we found eta
01472         DO L=1,4
01473           PETA(L)= phep(L,K)  ! eta 4 momentum
01474         ENDDO
01475 C       eta cumulated branching ratios:
01476         BRSUM(1)=0.389  ! gamma gamma
01477         BRSUM(2)=BRSUM(1)+0.319  ! 3 pi0
01478         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
01479         CALL RANMAR(RRR,1) 
01480         
01481         IF (RRR(1).LT.BRSUM(1)) THEN ! gamma gamma channel exactly like pi0
01482 ! random 3 vector on the sphere, masless   
01483          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
01484          CALL SPHERD(R,X) 
01485          X(4)=R
01486          Y(4)=R
01487         
01488          Y(1)=-X(1)
01489          Y(2)=-X(2)
01490          Y(3)=-X(3)
01491 ! boost to lab and to real*4
01492          CALL bostdq(-1,PETA,X,X)  
01493          CALL bostdq(-1,PETA,Y,Y)
01494          DO L=1,4
01495           PHOT1(L)=X(L)
01496           PHOT2(L)=Y(L)
01497          ENDDO
01498 C to hepevt
01499          CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
01500          CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
01501         ELSE ! 3 body channels
01502          IF(RRR(1).LT.BRSUM(2)) THEN  ! 3 pi0
01503           ID1= 111
01504           ID2= 111
01505           ID3= 111
01506           XM1=AMPIZ ! masses
01507           XM2=AMPIZ
01508           XM3=AMPIZ
01509          ELSEIF(RRR(1).LT.BRSUM(3)) THEN ! pi+ pi- pi0
01510           ID1= 211
01511           ID2=-211
01512           ID3= 111
01513           XM1=AMPI ! masses
01514           XM2=AMPI
01515           XM3=AMPIZ
01516          ELSE                            ! pi+ pi- gamma 
01517           ID1= 211
01518           ID2=-211
01519           ID3=  22
01520           XM1=AMPI ! masses
01521           XM2=AMPI
01522           XM3=0.0
01523          ENDIF
01524  7       CONTINUE  ! we generate mass of the first pair:
01525           CALL RANMAR(RR,2)
01526           R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)
01527           AMIN=XM1+XM2
01528           AMAX=R-XM3
01529           AM2=SQRT(AMIN**2+RR(1)*(AMAX**2-AMIN**2))
01530 C         weight for flat phase space
01531           WT=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)
01532      &      *XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)
01533      &           /R**2                    /AM2**2
01534          IF (RR(2).GT.WT) GOTO 7
01535 
01536          RU=XLAM(1D0*AM2**2,1D0*XM1**2,1D0*XM2**2)/AM2/2  ! momenta of the
01537                                               ! first two products
01538                                               ! in the rest frame of that pair
01539          CALL SPHERD(RU,X)
01540          X(4)=SQRT(RU**2+XM1**2)
01541          Y(4)=SQRT(RU**2+XM2**2)
01542         
01543          Y(1)=-X(1)
01544          Y(2)=-X(2)
01545          Y(3)=-X(3)
01546 C generate momentum of that pair in rest frame of eta:
01547          RU=XLAM(1D0*R**2,1D0*AM2**2,1D0*XM3**2)/R/2
01548          CALL SPHERD(RU,Z)
01549          Z(4)=SQRT(RU**2+AM2**2)
01550 C and boost first two decay products to rest frame of eta.
01551          CALL bostdq(-1,Z,X,X)
01552          CALL bostdq(-1,Z,Y,Y)
01553 C redefine Z(4) to 4-momentum of the last decay product: 
01554          Z(1)=-Z(1)
01555          Z(2)=-Z(2)
01556          Z(3)=-Z(3)
01557          Z(4)=SQRT(RU**2+XM3**2)
01558 C boost all to lab and move to real*4; also masses
01559          CALL bostdq(-1,PETA,X,X)
01560          CALL bostdq(-1,PETA,Y,Y)
01561          CALL bostdq(-1,PETA,Z,Z)
01562          DO L=1,4
01563           PHOT1(L)=X(L)
01564           PHOT2(L)=Y(L)
01565           PHOT3(L)=Z(L)
01566          ENDDO
01567          YM1=XM1
01568          YM2=XM2
01569          YM3=XM3
01570 C to hepevt
01571          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
01572          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
01573          CALL FILHEP(0,1,ID3,K,K,0,0,PHOT3,YM3,.TRUE.)
01574         ENDIF
01575 
01576        ENDIF
01577       ENDDO
01578 C
01579       END
01580       SUBROUTINE TAUK0S(JAK)
01581 C subroutine to decay K0S's from taus. 
01582 C this routine belongs to tauola universal interface, but uses 
01583 C routines from tauola utilities. Just flat phase space, but 4 channels.
01584 C it is called at the beginning of SUBR. TAUPI0(JAK)
01585 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005    
01586 C this is the hepevt class in old style. No d_h_ class pre-name
01587       INTEGER NMXHEP
01588       PARAMETER (NMXHEP=10000)
01589       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01590       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01591      $        jdahep
01592       COMMON /hepevt/
01593      $      nevhep,               ! serial number
01594      $      nhep,                 ! number of particles
01595      $      isthep(nmxhep),   ! status code
01596      $      idhep(nmxhep),    ! particle ident KF
01597      $      jmohep(2,nmxhep), ! parent particles
01598      $      jdahep(2,nmxhep), ! childreen particles
01599      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01600      $      vhep(4,nmxhep)    ! vertex [mm]
01601 * ----------------------------------------------------------------------
01602       LOGICAL qedrad
01603       COMMON /phoqed/ 
01604      $     qedrad(nmxhep)    ! Photos flag
01605 * ----------------------------------------------------------------------
01606       SAVE hepevt,phoqed
01607 
01608       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01609      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01610      *                 ,AMK,AMKZ,AMKST,GAMKST
01611 *
01612       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01613      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01614      *                 ,AMK,AMKZ,AMKST,GAMKST
01615 
01616 C position of taus, must be defined by host program:
01617       COMMON /TAUPOS/ NP1,NP2
01618 c
01619       REAL  RRR(1),BRSUM(3), RR(2)
01620       REAL  PHOT1(4),PHOT2(4),PHOT3(4)
01621       REAL*8    X(4),    Y(4),    Z(4)
01622       REAL                                YM1,YM2,YM3
01623       REAL*8  R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
01624       REAL*8 a,b,c
01625       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
01626 C position of decaying particle:
01627       IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
01628         NPS=NP1
01629       ELSE
01630         NPS=NP2
01631       ENDIF
01632       nhepM=nhep                ! to avoid infinite loop              
01633       DO K=JDAHEP(1,NPS),nhepM  ! we search for K0S's from tau till eor.
01634        IF (IDHEP(K).EQ.310.AND.JDAHEP(1,K).LE.K) THEN ! IF we found K0S
01635 
01636       
01637         DO L=1,4
01638           PETA(L)= phep(L,K)  ! K0S 4 momentum  (this is cloned from eta decay)
01639         ENDDO
01640 C       K0S cumulated branching ratios:
01641         BRSUM(1)=0.313  ! 2 PI0
01642         BRSUM(2)=1.0 ! BRSUM(1)+0.319  ! Pi+ PI-
01643         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
01644         CALL RANMAR(RRR,1) 
01645 
01646          IF(RRR(1).LT.BRSUM(1)) THEN  ! 2 pi0
01647           ID1= 111
01648           ID2= 111
01649           XM1=AMPIZ ! masses
01650           XM2=AMPIZ
01651          ELSEIF(RRR(1).LT.BRSUM(2)) THEN ! pi+ pi- 
01652           ID1= 211
01653           ID2=-211
01654           XM1=AMPI ! masses
01655           XM2=AMPI
01656          ELSE                            ! gamma gamma unused !!!
01657           ID1= 22
01658           ID2= 22
01659           XM1= 0.0 ! masses
01660           XM2= 0.0
01661          ENDIF
01662         
01663 ! random 3 vector on the sphere, of equal mass !!  
01664          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
01665          R4=R
01666          R=SQRT(ABS(R**2-XM1**2))
01667          CALL SPHERD(R,X) 
01668          X(4)=R4
01669          Y(4)=R4
01670         
01671          Y(1)=-X(1)
01672          Y(2)=-X(2)
01673          Y(3)=-X(3)
01674 ! boost to lab and to real*4
01675          CALL bostdq(-1,PETA,X,X)  
01676          CALL bostdq(-1,PETA,Y,Y)
01677          DO L=1,4
01678           PHOT1(L)=X(L)
01679           PHOT2(L)=Y(L)
01680          ENDDO
01681 
01682          YM1=XM1
01683          YM2=XM2
01684 C to hepevt
01685          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
01686          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
01687 
01688 C
01689        ENDIF
01690       ENDDO
01691 
01692       END
Generated on Sun Oct 20 20:24:10 2013 for C++InterfacetoTauola by  doxygen 1.6.3