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

00001       SUBROUTINE TAUOLA(MODE,KEYPOL) 
00002 C     *************************************
00003 C general tauola interface, should work in every case until 
00004 C hepevt is OK, does not check if hepevt is 'clean'
00005 C in particular will decay decayed taus...
00006 C only longitudinal spin effects are included.
00007 C in W decay v-a vertex is assumed
00008 C date: 12 DEC 1998. date: 21 June 1999. date: 24 Jan 2001 date: 24 Aug 2001
00009 C this is the hepevt class in old style. No d_h_ class pre-name
00010       INTEGER NMXHEP
00011       PARAMETER (NMXHEP=4000)
00012       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
00013       INTEGER nevhep,nhep,isthep,idhep,jmohep,
00014      $        jdahep
00015       COMMON /hepevt/
00016      $      nevhep,               ! serial number
00017      $      nhep,                 ! number of particles
00018      $      isthep(nmxhep),   ! status code
00019      $      idhep(nmxhep),    ! particle ident KF
00020      $      jmohep(2,nmxhep), ! parent particles
00021      $      jdahep(2,nmxhep), ! childreen particles
00022      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
00023      $      vhep(4,nmxhep)    ! vertex [mm]
00024 * ----------------------------------------------------------------------
00025       LOGICAL qedrad
00026       COMMON /phoqed/ 
00027      $     qedrad(nmxhep)    ! Photos flag
00028 * ----------------------------------------------------------------------
00029       SAVE hepevt,phoqed
00030 
00031 
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,1) means on,  (-1,0) off
00095          ION(1)=1
00096          ION(2)=1
00097          ION(3)=1
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.21, September 2005******',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.21, September2005 ******',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 PLZAPX,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       INTEGER IM
00425       LOGICAL HOPE,HOPEin
00426 C this is the hepevt class in old style. No d_h_ class pre-name
00427       INTEGER NMXHEP
00428       PARAMETER (NMXHEP=4000)
00429       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
00430       INTEGER nevhep,nhep,isthep,idhep,jmohep,
00431      $        jdahep
00432       COMMON /hepevt/
00433      $      nevhep,               ! serial number
00434      $      nhep,                 ! number of particles
00435      $      isthep(nmxhep),   ! status code
00436      $      idhep(nmxhep),    ! particle ident KF
00437      $      jmohep(2,nmxhep), ! parent particles
00438      $      jdahep(2,nmxhep), ! childreen particles
00439      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
00440      $      vhep(4,nmxhep)    ! vertex [mm]
00441 * ----------------------------------------------------------------------
00442       LOGICAL qedrad
00443       COMMON /phoqed/ 
00444      $     qedrad(nmxhep)    ! Photos flag
00445 * ----------------------------------------------------------------------
00446       SAVE hepevt,phoqed
00447 
00448 
00449 
00450 C(BPK)--> BROTHERS OF TAU ALREADY FOUND
00451       INTEGER ISON
00452       COMMON /ISONS_TAU/ISON(2)
00453 C(BPK)--<
00454 C >>
00455 C >> STEP 1: find where are particles in hepevent and pick them up
00456 C >>
00457              HOPE=HOPEin
00458 C sometimes shade Z of Z is its mother ...
00459             IM=IM0
00460             IM00=JMOHEP(1,IM0)
00461 C to protect against check on mother of beam particles.
00462             IF (IM00.GT.0) THEN
00463               IF (IDHEP(IM0).EQ.IDHEP(IM00)) IM=JMOHEP(1,IM0)
00464             ENDIF                                                              
00465 C
00466 C find (host generator-level) incoming beam-bare-particles which form Z and co.
00467             IMO1=JMOHEP(1,IM)
00468             IMO2=JMOHEP(2,IM)
00469 
00470 C(BPK)--> IN HERWIG THE POINTER MIGHT BE TO HARD CMS
00471             IM00=IMO1
00472             IF (ISTHEP(IM00).EQ.120) THEN
00473                IMO1=JMOHEP(1,IM00)
00474                IMO2=JMOHEP(2,IM00)
00475             ENDIF
00476 C(BPK)--<
00477 
00478             IFFULL=0
00479 C case when it was like e+e- --> tau+tau- gammas and e+e- were 'first' in hepevt.
00480       IF (IMO1.EQ.0.AND.IMO2.EQ.0) THEN
00481          IMO1=JMOHEP(1,NP1)
00482 C(BPK)-->
00483          IF (IDHEP(IMO1).EQ.IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
00484 C(BPK)--<
00485          IMO2=JMOHEP(2,NP1)
00486 C(BPK)-->
00487          IF (IDHEP(IMO2).EQ.IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
00488 C(BPK)--<
00489          IFFULL=1
00490 C case when it was like  qq --> tau+tau- gammas and qq were NOT 'first' in hepevt.
00491 
00492       ELSEIF (IDHEP(IM).NE.22.AND.IDHEP(IM).NE.23) THEN
00493          IMO1=JMOHEP(1,NP1)
00494 C(BPK)-->
00495          IF (IDHEP(IMO1).EQ.IDHEP(NP1)) IMO1=JMOHEP(1,IMO1) ! PROTECT AGAINST COPIES
00496 C(BPK)--<
00497          IMO2=JMOHEP(2,NP1)
00498 C(BPK)-->
00499          IF (IDHEP(IMO2).EQ.IDHEP(NP2)) IMO2=JMOHEP(1,IMO2) ! PROTECT AGAINST COPIES
00500 C(BPK)--<
00501          IFFULL=1
00502       ENDIF
00503 
00504             
00505 C and check if it really happened
00506             IF (IMO1.EQ.0) HOPE=.FALSE.
00507             IF (IMO2.EQ.0) HOPE=.FALSE.
00508             IF (IMO1.EQ.IMO2) HOPE=.FALSE.
00509 
00510 C 
00511          DO I=1,4
00512             Q1(I)= PHEP(I,NP1)              !momentum of tau+
00513             Q2(I)= PHEP(I,NP2)              !momentum of tau-
00514          ENDDO
00515 
00516 C corrections due to possible differences in 4-momentum of shadow vs true Z.
00517 C(BPK)-->
00518       IF (IM.EQ.JMOHEP(1,IM0).AND.
00519      $     (IDHEP(IM).EQ.22.OR.IDHEP(IM).EQ.23)) THEN
00520          DO K=1,4
00521             PB(K)=PHEP(K,IM)
00522             PA(K)=PHEP(K,IM0)
00523          ENDDO
00524 C(BPK)--<
00525 
00526                CALL BOSTDQ( 1,PA, Q1, Q1)
00527                CALL BOSTDQ( 1,PA, Q2, Q2)
00528                CALL BOSTDQ(-1,PB, Q1, Q1)
00529                CALL BOSTDQ(-1,PB, Q2, Q2)
00530 
00531               ENDIF
00532 
00533          DO I=1,4                                                                                   
00534             QQ(I)= Q1(I)+Q2(I)              !momentum of Z
00535             IF (HOPE) P1(I)=PHEP(I,IMO1)    !momentum of beam1
00536             IF (HOPE) P2(I)=PHEP(I,IMO2)    !momentum of beam2
00537             PH(I)=0D0
00538             PD1(I)=0D0
00539             PD2(I)=0D0
00540          ENDDO
00541 !        These momenta correspond to  quarks, gluons photons or taus
00542                    IDFQ1=IDHEP(NP1)
00543                    IDFQ2=IDHEP(NP2)
00544          IF (HOPE) IDFP1=IDHEP(IMO1)
00545          IF (HOPE) IDFP2=IDHEP(IMO2)
00546 
00547          SVAR=QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2
00548          IF (.NOT.HOPE) THEN
00549 C options which may be useful in some cases of two heavy boson production
00550 C need individual considerations. To be developed.
00551 C           PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0)  ! gamma/Z mixture as if produced from e beam
00552 C           PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0)  ! pure Z
00553            PLZAPX=0.5                         ! pure gamma
00554            RETURN
00555          ENDIF
00556 C >>
00557 C >> STEP 2 look for brothers of Z which have to be included in effective incoming particles
00558 C >>
00559 C let us define beginning and end of particles which are produced in parallel to Z
00560 C in principle following should work
00561 
00562 C(BPK)--> ACCOMMODATE FOR HERWIG - IM00 POINTS TO BEAM PARTICLE OR HARD CMS
00563          NX1=JDAHEP(1,IM00)
00564          NX2=JDAHEP(2,IM00)
00565 C but ...
00566          INBR=IM ! OK, HARD RECORD Z/GAMMA 
00567          IF (IFFULL.EQ.1) INBR=NP1  ! OK, NO Z/GAMMA
00568          IF (IDHEP(JMOHEP(1,INBR)).EQ.IDHEP(INBR)) INBR=JMOHEP(1,INBR) ! FORCE HARD RECORD
00569 C(BPK)--<
00570          IF(NX1.EQ.0.OR.NX2.EQ.0) THEN
00571            NX1=INBR
00572            NX2=INBR
00573            DO K=1,INBR-1
00574              IF(JMOHEP(1,INBR-K).EQ.JMOHEP(1,INBR)) THEN
00575               NX1=INBR-K
00576              ELSE
00577                 GOTO 7
00578              ENDIF
00579            ENDDO
00580  7         CONTINUE
00581 
00582            DO K=INBR+1,NHEP
00583              IF(JMOHEP(1,K).EQ.JMOHEP(1,INBR)) THEN
00584               NX2=K
00585              ELSE
00586                 GOTO 8
00587              ENDIF
00588            ENDDO
00589  8         CONTINUE
00590          ENDIF
00591 
00592 C case of annihilation of two bosons is hopeles
00593          IF (ABS(IDFP1).GE.20.AND.ABS(IDFP2).GE.20) HOPE=.FALSE.
00594 C case of annihilation of non-matching flavors is hopeless
00595          IF (ABS(IDFP1).LE.20.AND.ABS(IDFP2).LE.20.AND.IDFP1+IDFP2.NE.0)
00596      $       HOPE=.FALSE.
00597          IF (.NOT.HOPE) THEN
00598 C options which may be useful in some cases of two heavy boson production
00599 C need individual considerations. To be developed.
00600 C           PLZAPX=PLZAP0(11,IDFQ1,SVAR,0D0)  ! gamma/Z mixture as if produced from e beam
00601 C           PLZAPX=PLZAP0(12,IDFQ1,SVAR,0D0)  ! pure Z
00602            PLZAPX=0.5                         ! pure gamma
00603            RETURN
00604          ENDIF
00605          IF (ABS(IDFP1).LT.20) IDE= IDFP1
00606          IF (ABS(IDFP2).LT.20) IDE=-IDFP2
00607 
00608 
00609 C >>
00610 C >> STEP 3  we combine gluons, photons  into incoming effective beams
00611 C >>
00612 
00613 C in the following we will ignore the possibility of photon emission from taus
00614 C however at certain moment it will be necessary to take care of
00615 
00616            DO L=1,4
00617             PD1(L)=P1(L)
00618             PD2(L)=P2(L)
00619            ENDDO 
00620 
00621                 DO L=1,4
00622                 PQ1(L)=Q1(L)
00623                 PQ2(L)=Q2(L)
00624                 ENDDO 
00625 
00626          IFLAV=min(ABS(IDFP1),ABS(IDFP2))
00627 
00628 *--------------------------------------------------------------------------
00629 c                  IFLAV=min(ABS(IDFP1),ABS(IDFP2))
00630 c that means that always quark or lepton i.e. process like
00631 c      f g(gamma) --> f Z0 --> tau tau
00632 c we glue  fermions to effective  beams that is f f --> Z0 --> tau tau
00633 c  with  gamma/g emission from initial fermion.
00634 *---------------------------------------------------------------------------
00635  
00636          IF (ABS(IDFP1).GE.20) THEN
00637            DO k=NX1,NX2
00638              IDP=IDHEP(k)
00639              IF (ABS(IDP).EQ.IFLAV) THEN
00640                DO L=1,4
00641                  PD1(L)=-PHEP(L,K)
00642                ENDDO
00643              ENDIF
00644            ENDDO
00645          ENDIF
00646 
00647          IF (ABS(IDFP2).GE.20) THEN
00648            DO k=NX1,NX2
00649              IDP=IDHEP(k)
00650              IF (ABS(IDP).EQ.IFLAV) THEN
00651                DO L=1,4
00652                  PD2(L)=-PHEP(L,K)
00653                ENDDO
00654              ENDIF
00655            ENDDO
00656          ENDIF
00657 
00658 C if first beam was boson: gluing
00659 
00660          IF (ABS(IDFP1).GE.20) THEN
00661            DO L=1,4
00662              PH(L)=P1(L)
00663            ENDDO
00664            xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
00665      $            -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
00666            xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
00667      $            -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
00668           IF (XM1.LT.XM2) THEN
00669              DO L=1,4
00670                PD1(L)=PD1(L)+P1(L)
00671              ENDDO
00672            ELSE
00673              DO L=1,4
00674                PD2(L)=PD2(L)+P1(L)
00675              ENDDO
00676            ENDIF
00677          ENDIF
00678 
00679 C if second beam was boson: gluing 
00680 
00681 
00682          IF (ABS(IDFP2).GE.20) THEN
00683            DO L=1,4
00684              PH(L)=P2(L)
00685            ENDDO
00686            xm1=abs((PD1(4)+PH(4))**2-(PD1(3)+PH(3))**2
00687      $            -(PD1(2)+PH(2))**2-(PD1(1)+PH(1))**2)
00688            xm2=abs((PD2(4)+PH(4))**2-(PD2(3)+PH(3))**2
00689      $            -(PD2(2)+PH(2))**2-(PD2(1)+PH(1))**2)
00690            IF (XM1.LT.XM2) THEN
00691              DO L=1,4
00692                PD1(L)=PD1(L)+P2(L)
00693              ENDDO
00694            ELSE
00695              DO L=1,4
00696                PD2(L)=PD2(L)+P2(L)
00697              ENDDO
00698            ENDIF
00699          ENDIF
00700 
00701 C now spectators ...
00702 
00703 C(BPK)-->
00704       NPH1=NP1
00705       NPH2=NP2
00706       IF (IDHEP(JMOHEP(1,NP1)).EQ.IDHEP(NP1)) NPH1=JMOHEP(1,NP1) ! SHOULD PUT US IN HARD REC.
00707       IF (IDHEP(JMOHEP(1,NP2)).EQ.IDHEP(NP2)) NPH2=JMOHEP(1,NP2) ! SHOULD PUT US IN HARD REC.
00708 C(BPK)--<
00709 
00710          DO k=NX1,NX2
00711          IF (ABS(IDHEP(K)).NE.IFLAV.AND.K.NE.IM.AND.
00712 C(BPK)-->
00713      $       K.NE.NPH1.AND.K.NE.NPH2) THEN 
00714 C(BPK)--<
00715           IF(IDHEP(K).EQ.22.AND.IFFULL.EQ.1) THEN
00716             DO L=1,4
00717               PH(L)=PHEP(L,K)
00718             ENDDO
00719             xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
00720      $             -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
00721             xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
00722      $             -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
00723            xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
00724      $            -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
00725            xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
00726      $            -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)
00727 
00728   
00729             sini=abs((PD1(4)+PD2(4)-PH(4))**2-(PD1(3)+PD2(3)-PH(3))**2
00730      $              -(PD1(2)+PD2(2)-PH(2))**2-(PD1(1)+PD2(1)-PH(1))**2)
00731             sfin=abs((PD1(4)+PD2(4)      )**2-(PD1(3)+PD2(3)      )**2
00732      $              -(PD1(2)+PD2(2)      )**2-(PD1(1)+PD2(1)      )**2)
00733 
00734            FACINI=ZPROP2(sini)
00735            FACFIN=ZPROP2(sfin)
00736 
00737            XM1=XM1/FACINI
00738            XM2=XM2/FACINI
00739            XM3=XM3/FACFIN
00740            XM4=XM4/FACFIN
00741 
00742            XM=MIN(XM1,XM2,XM3,XM4)
00743                   IF      (XM1.EQ.XM) THEN 
00744                      DO L=1,4
00745                        PD1(L)=PD1(L)-PH(L)
00746                      ENDDO
00747                   ELSEIF   (XM2.EQ.XM) THEN 
00748                      DO L=1,4
00749                        PD2(L)=PD2(L)-PH(L)
00750                      ENDDO
00751                   ELSEIF   (XM3.EQ.XM) THEN 
00752                      DO L=1,4
00753                         Q1(L)=PQ1(L)+PH(L)
00754                      ENDDO
00755                   ELSE
00756                      DO L=1,4
00757                         Q2(L)=PQ2(L)+PH(L)
00758                      ENDDO
00759                   ENDIF
00760            ELSE
00761             DO L=1,4
00762               PH(L)=PHEP(L,K)
00763             ENDDO
00764             xm1=abs((PD1(4)-PH(4))**2-(PD1(3)-PH(3))**2
00765      $             -(PD1(2)-PH(2))**2-(PD1(1)-PH(1))**2)
00766             xm2=abs((PD2(4)-PH(4))**2-(PD2(3)-PH(3))**2
00767      $             -(PD2(2)-PH(2))**2-(PD2(1)-PH(1))**2)
00768             IF (XM1.LT.XM2) THEN
00769               DO L=1,4
00770                 PD1(L)=PD1(L)-PH(L)
00771               ENDDO
00772             ELSE
00773               DO L=1,4
00774                 PD2(L)=PD2(L)-PH(L)
00775               ENDDO
00776             ENDIF
00777            ENDIF
00778           ENDIF
00779          ENDDO
00780 
00781 
00782 C >>
00783 C >> STEP 4 look for brothers of tau (sons of Z!) which have to be included in 
00784 c >>          effective outcoming taus
00785 C >>
00786 C let us define beginning and end of particles which are produced in 
00787 c  parallel to tau
00788 
00789 
00790 
00791 C find outcoming particles which come from Z   
00792 
00793        
00794  
00795 
00796 C(BPK)--> OK, IT WOULD HAVE TO BE ALONG TAUS IN HARD RECORD WITH THE SAME MOTHER       
00797       IF (ABS(IDHEP(IM0)).EQ.22.OR.abs(IDHEP(IM0)).EQ.23) THEN
00798          DO K=ISON(1),ISON(2)
00799             IF(ABS(IDHEP(K)).EQ.22) THEN
00800 C(BPK)--< 
00801 
00802               do l=1,4
00803               ph(l)=phep(l,k)
00804               enddo
00805 
00806            xm3=abs((PQ1(4)+PH(4))**2-(PQ1(3)+PH(3))**2
00807      $            -(PQ1(2)+PH(2))**2-(PQ1(1)+PH(1))**2)
00808            xm4=abs((PQ2(4)+PH(4))**2-(PQ2(3)+PH(3))**2
00809      $            -(PQ2(2)+PH(2))**2-(PQ2(1)+PH(1))**2)  
00810 
00811            XM=MIN(XM3,XM4) 
00812 
00813                   IF   (XM3.EQ.XM) THEN 
00814                      DO L=1,4
00815                         Q1(L)=PQ1(L)+PH(L)
00816                      ENDDO
00817                   ELSE
00818                      DO L=1,4
00819                         Q2(L)=PQ2(L)+PH(L)
00820                      ENDDO
00821                   ENDIF
00822             endif
00823           enddo
00824           ENDIF
00825 
00826 
00827 
00828 *------------------------------------------------------------------------
00829 
00830 
00831 C out of effective momenta we calculate COSTHE and later polarization
00832       CALL ANGULU(PD1,PD2,Q1,Q2,COSTHE)
00833      
00834       PLZAPX=PLZAP0(IDE,IDFQ1,SVAR,COSTHE)
00835       END
00836 
00837       SUBROUTINE ANGULU(PD1,PD2,Q1,Q2,COSTHE)
00838       REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
00839 C take effective beam which is less massive, it should be irrelevant
00840 C but in case HEPEVT is particulary dirty may help.
00841 C this routine calculate reduced system transver and cosine of scattering 
00842 C angle.
00843 
00844       XM1=ABS(PD1(4)**2-PD1(3)**2-PD1(2)**2-PD1(1)**2)
00845       XM2=ABS(PD2(4)**2-PD2(3)**2-PD2(2)**2-PD2(1)**2)
00846       IF (XM1.LT.XM2) THEN
00847         SIGN=1D0
00848         DO K=1,4
00849           P(K)=PD1(K)
00850         ENDDO
00851       ELSE
00852         SIGN=-1D0
00853         DO K=1,4
00854           P(K)=PD2(K)
00855         ENDDO
00856       ENDIF
00857 C calculate space like part of P (in Z restframe)
00858       DO K=1,4
00859        QQ(K)=Q1(k)+Q2(K)
00860        QT(K)=Q1(K)-Q2(K)
00861       ENDDO
00862 
00863        XMQQ=SQRT(QQ(4)**2-QQ(3)**2-QQ(2)**2-QQ(1)**2)
00864 
00865        QTXQQ=QT(4)*QQ(4)-QT(3)*QQ(3)-QT(2)*QQ(2)-QT(1)*QQ(1)
00866       DO K=1,4
00867        QT(K)=QT(K)-QQ(K)*QTXQQ/XMQQ**2
00868       ENDDO
00869 
00870        PXQQ=P(4)*QQ(4)-P(3)*QQ(3)-P(2)*QQ(2)-P(1)*QQ(1)
00871       DO K=1,4
00872        P(K)=P(K)-QQ(K)*PXQQ/XMQQ**2
00873       ENDDO
00874 C calculate costhe
00875        PXP  =SQRT(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
00876        QTXQT=SQRT(QT(3)**2+QT(2)**2+QT(1)**2-QT(4)**2)
00877        PXQT =P(3)*QT(3)+P(2)*QT(2)+P(1)*QT(1)-P(4)*QT(4)
00878        COSTHE=PXQT/PXP/QTXQT
00879        COSTHE=COSTHE*SIGN
00880       END
00881 
00882       FUNCTION PLZAP0(IDE,IDF,SVAR,COSTH0)
00883 C this function calculates probability for the helicity +1 +1 configuration
00884 C of taus for given Z/gamma transfer and COSTH0 cosine of scattering angle
00885       REAL*8 PLZAP0,SVAR,COSTHE,COSTH0
00886 
00887       COSTHE=COSTH0
00888 C >>>>>      IF (IDE*IDF.LT.0) COSTHE=-COSTH0 ! this is probably not needed ID
00889 C >>>>>      of first beam is used by T_GIVIZ0 including sign
00890 
00891       IF (IDF.GT.0) THEN
00892         CALL INITWK(IDE,IDF,SVAR)
00893       ELSE
00894         CALL INITWK(-IDE,-IDF,SVAR)
00895       ENDIF
00896       PLZAP0=T_BORN(0,SVAR,COSTHE,1D0,1D0)
00897      $  /(T_BORN(0,SVAR,COSTHE,1D0,1D0)+T_BORN(0,SVAR,COSTHE,-1D0,-1D0))
00898 
00899 !      PLZAP0=0.5
00900       END
00901       FUNCTION T_BORN(MODE,SVAR,COSTHE,TA,TB)
00902 C ----------------------------------------------------------------------
00903 C THIS ROUTINE PROVIDES BORN CROSS SECTION. IT HAS THE SAME         
00904 C STRUCTURE AS FUNTIS AND FUNTIH, THUS CAN BE USED AS SIMPLER       
00905 C EXAMPLE OF THE METHOD APPLIED THERE                               
00906 C INPUT PARAMETERS ARE: SVAR    -- transfer
00907 C                       COSTHE  -- cosine of angle between tau+ and 1st beam
00908 C                       TA,TB   -- helicity states of tau+ tau-
00909 C
00910 C     called by : BORNY, BORAS, BORNV, WAGA, WEIGHT
00911 C ----------------------------------------------------------------------
00912       IMPLICIT REAL*8(A-H,O-Z)
00913       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
00914       REAL*8              ENE ,AMIN,AMFIN
00915       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
00916      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
00917      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
00918      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
00919       REAL*8             SS,POLN,T3E,QE,T3F,QF
00920      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
00921       REAL*8            SEPS1,SEPS2
00922 C=====================================================================
00923       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00924       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
00925 C     SWSQ        = sin2 (theta Weinberg)
00926 C     AMW,AMZ     = W & Z boson masses respectively
00927 C     AMH         = the Higgs mass
00928 C     AMTOP       = the top mass
00929 C     GAMMZ       = Z0 width
00930       COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
00931       COMPLEX*16 XUPZFP(2),XUPZIP(2)
00932       COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
00933       COMPLEX*16 PROPA,PROPZ
00934       COMPLEX*16 XR,XI
00935       COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
00936       COMPLEX*16 XTHING,XVE,XVF,XVEF
00937       DATA XI/(0.D0,1.D0)/,XR/(1.D0,0.D0)/
00938       DATA MODE0 /-5/
00939       DATA IDE0 /-55/
00940       DATA SVAR0,COST0 /-5.D0,-6.D0/
00941       DATA PI /3.141592653589793238462643D0/
00942       DATA SEPS1,SEPS2 /0D0,0D0/
00943 C
00944 C MEMORIZATION =========================================================
00945       IF ( MODE.NE.MODE0.OR.SVAR.NE.SVAR0.OR.COSTHE.NE.COST0
00946      $    .OR.IDE0.NE.IDE)THEN
00947 C
00948         KEYGSW=1
00949 C ** PROPAGATORS
00950         IDE0=IDE
00951         MODE0=MODE
00952         SVAR0=SVAR
00953         COST0=COSTHE
00954         SINTHE=SQRT(1.D0-COSTHE**2)
00955         BETA=SQRT(MAX(0D0,1D0-4D0*AMFIN**2/SVAR))
00956 C I MULTIPLY AXIAL COUPLING BY BETA FACTOR.
00957         XUPZFP(1)=0.5D0*(XUPZF(1)+XUPZF(2))+0.5*BETA*(XUPZF(1)-XUPZF(2))
00958         XUPZFP(2)=0.5D0*(XUPZF(1)+XUPZF(2))-0.5*BETA*(XUPZF(1)-XUPZF(2))
00959         XUPZIP(1)=0.5D0*(XUPZI(1)+XUPZI(2))+0.5*(XUPZI(1)-XUPZI(2))
00960         XUPZIP(2)=0.5D0*(XUPZI(1)+XUPZI(2))-0.5*(XUPZI(1)-XUPZI(2))
00961 C FINAL STATE VECTOR COUPLING
00962         XUPF     =0.5D0*(XUPZF(1)+XUPZF(2))
00963         XUPI     =0.5D0*(XUPZI(1)+XUPZI(2))
00964         XTHING   =0D0
00965 
00966         PROPA =1D0/SVAR
00967         PROPZ =1D0/DCMPLX(SVAR-AMZ**2,SVAR/AMZ*GAMMZ)
00968         IF (KEYGSW.EQ.0) PROPZ=0.D0
00969         DO 50 I=1,2
00970          DO 50 J=1,2
00971           REGULA= (3-2*I)*(3-2*J) + COSTHE
00972           REGULM=-(3-2*I)*(3-2*J) * SINTHE *2.D0*AMFIN/SQRT(SVAR)
00973           APHOT(I,J)=PROPA*(XUPGI(I)*XUPGF(J)*REGULA)
00974           AZETT(I,J)=PROPZ*(XUPZIP(I)*XUPZFP(J)+XTHING)*REGULA
00975           ABORN(I,J)=APHOT(I,J)+AZETT(I,J)
00976           APHOTM(I,J)=PROPA*DCMPLX(0D0,1D0)*XUPGI(I)*XUPGF(J)*REGULM
00977           AZETTM(I,J)=PROPZ*DCMPLX(0D0,1D0)*(XUPZIP(I)*XUPF+XTHING)*REGULM
00978           ABORNM(I,J)=APHOTM(I,J)+AZETTM(I,J)
00979    50   CONTINUE
00980       ENDIF
00981 C
00982 C******************
00983 C* IN CALCULATING CROSS SECTION ONLY DIAGONAL ELEMENTS
00984 C* OF THE SPIN DENSITY MATRICES ENTER (LONGITUD. POL. ONLY.)
00985 C* HELICITY CONSERVATION EXPLICITLY OBEYED
00986       POLAR1=  (SEPS1)
00987       POLAR2= (-SEPS2)
00988       BORN=0D0
00989       DO 150 I=1,2
00990        HELIC= 3-2*I
00991        DO 150 J=1,2
00992         HELIT=3-2*J
00993         FACTOR=KOLOR*(1D0+HELIC*POLAR1)*(1D0-HELIC*POLAR2)/4D0
00994         FACTOM=FACTOR*(1+HELIT*TA)*(1-HELIT*TB)
00995         FACTOR=FACTOR*(1+HELIT*TA)*(1+HELIT*TB)
00996 
00997         BORN=BORN+CDABS(ABORN(I,J))**2*FACTOR
00998 C      MASS TERM IN BORN
00999         IF (MODE.GE.1) THEN
01000          BORN=BORN+CDABS(ABORNM(I,J))**2*FACTOM
01001         ENDIF
01002 
01003   150 CONTINUE
01004 C************
01005       FUNT=BORN
01006       IF(FUNT.LT.0.D0)  FUNT=BORN
01007 
01008 C
01009       IF (SVAR.GT.4D0*AMFIN**2) THEN
01010 C PHASE SPACE THRESHOLD FACTOR
01011         THRESH=SQRT(1-4D0*AMFIN**2/SVAR)
01012         T_BORN= FUNT*SVAR**2*THRESH
01013       ELSE
01014         THRESH=0.D0
01015         T_BORN=0.D0
01016       ENDIF
01017 C ZW HERE WAS AN ERROR 19. 05. 1989
01018 !      write(*,*) 'KKKK ',PROPA,PROPZ,XUPGI,XUPGF,XUPZI,XUPZF
01019 !      write(*,*) 'KKKK X',svar,costhe,TA,TB,T_BORN
01020       END
01021 
01022       SUBROUTINE INITWK(IDEX,IDFX,SVAR)
01023 ! initialization routine coupling masses etc.
01024       IMPLICIT REAL*8 (A-H,O-Z)
01025       COMMON / T_BEAMPM / ENE ,AMIN,AMFIN,IDE,IDF
01026       REAL*8              ENE ,AMIN,AMFIN
01027       COMMON / T_GAUSPM /SS,POLN,T3E,QE,T3F,QF
01028      &                  ,XUPGI   ,XUPZI   ,XUPGF   ,XUPZF
01029      &                  ,NDIAG0,NDIAGA,KEYA,KEYZ
01030      &                  ,ITCE,JTCE,ITCF,JTCF,KOLOR
01031       REAL*8             SS,POLN,T3E,QE,T3F,QF
01032      &                  ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
01033       COMMON / T_GSWPRM /SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
01034       REAL*8             SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
01035 C     SWSQ        = sin2 (theta Weinberg)
01036 C     AMW,AMZ     = W & Z boson masses respectively
01037 C     AMH         = the Higgs mass
01038 C     AMTOP       = the top mass
01039 C     GAMMZ       = Z0 width
01040 C
01041       ENE=SQRT(SVAR)/2
01042       AMIN=0.511D-3
01043       SWSQ=0.23147
01044       AMZ=91.1882
01045       GAMMZ=2.4952
01046       IF     (IDFX.EQ. 15) then       
01047         IDF=2  ! denotes tau +2 tau-
01048         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
01049       ELSEIF (IDFX.EQ.-15) then
01050         IDF=-2  ! denotes tau -2 tau-
01051         AMFIN=1.77703 !this mass is irrelevant if small, used in ME only
01052       ELSE
01053         WRITE(*,*) 'INITWK: WRONG IDFX'
01054         STOP
01055       ENDIF
01056 
01057       IF     (IDEX.EQ. 11) then      !electron
01058         IDE= 2
01059         AMIN=0.511D-3
01060       ELSEIF (IDEX.EQ.-11) then      !positron
01061         IDE=-2
01062         AMIN=0.511D-3
01063       ELSEIF (IDEX.EQ. 13) then      !mu+
01064         IDE= 2
01065         AMIN=0.105659
01066       ELSEIF (IDEX.EQ.-13) then      !mu-
01067         IDE=-2
01068         AMIN=0.105659
01069       ELSEIF (IDEX.EQ.  1) then      !d
01070         IDE= 4
01071         AMIN=0.05
01072       ELSEIF (IDEX.EQ.- 1) then      !d~
01073         IDE=-4
01074         AMIN=0.05
01075       ELSEIF (IDEX.EQ.  2) then      !u
01076         IDE= 3
01077         AMIN=0.02
01078       ELSEIF (IDEX.EQ.- 2) then      !u~
01079         IDE=-3
01080         AMIN=0.02
01081       ELSEIF (IDEX.EQ.  3) then      !s
01082         IDE= 4
01083         AMIN=0.3
01084       ELSEIF (IDEX.EQ.- 3) then      !s~
01085         IDE=-4
01086         AMIN=0.3
01087       ELSEIF (IDEX.EQ.  4) then      !c
01088         IDE= 3
01089         AMIN=1.3
01090       ELSEIF (IDEX.EQ.- 4) then      !c~
01091         IDE=-3
01092         AMIN=1.3
01093       ELSEIF (IDEX.EQ.  5) then      !b
01094         IDE= 4
01095         AMIN=4.5
01096       ELSEIF (IDEX.EQ.- 5) then      !b~
01097         IDE=-4
01098         AMIN=4.5
01099       ELSEIF (IDEX.EQ.  12) then     !nu_e
01100         IDE= 1
01101         AMIN=0.1D-3
01102       ELSEIF (IDEX.EQ.- 12) then     !nu_e~
01103         IDE=-1
01104         AMIN=0.1D-3
01105       ELSEIF (IDEX.EQ.  14) then     !nu_mu
01106         IDE= 1
01107         AMIN=0.1D-3
01108       ELSEIF (IDEX.EQ.- 14) then     !nu_mu~
01109         IDE=-1
01110         AMIN=0.1D-3
01111       ELSEIF (IDEX.EQ.  16) then     !nu_tau
01112         IDE= 1
01113         AMIN=0.1D-3
01114       ELSEIF (IDEX.EQ.- 16) then     !nu_tau~
01115         IDE=-1
01116         AMIN=0.1D-3
01117 
01118       ELSE
01119         WRITE(*,*) 'INITWK: WRONG IDEX'
01120         STOP
01121       ENDIF
01122 
01123 C ----------------------------------------------------------------------
01124 C
01125 C     INITIALISATION OF COUPLING CONSTANTS AND FERMION-GAMMA / Z0 VERTEX
01126 C
01127 C     called by : KORALZ
01128 C ----------------------------------------------------------------------
01129       ITCE=IDE/IABS(IDE)
01130       JTCE=(1-ITCE)/2
01131       ITCF=IDF/IABS(IDF)
01132       JTCF=(1-ITCF)/2
01133       CALL T_GIVIZO( IDE, 1,AIZOR,QE,KDUMM)
01134       CALL T_GIVIZO( IDE,-1,AIZOL,QE,KDUMM)
01135       XUPGI(1)=QE
01136       XUPGI(2)=QE
01137       T3E    = AIZOL+AIZOR
01138       XUPZI(1)=(AIZOR-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01139       XUPZI(2)=(AIZOL-QE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01140       CALL T_GIVIZO( IDF, 1,AIZOR,QF,KOLOR)
01141       CALL T_GIVIZO( IDF,-1,AIZOL,QF,KOLOR)
01142       XUPGF(1)=QF
01143       XUPGF(2)=QF
01144       T3F    =  AIZOL+AIZOR
01145       XUPZF(1)=(AIZOR-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01146       XUPZF(2)=(AIZOL-QF*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01147 C
01148       NDIAG0=2
01149       NDIAGA=11
01150       KEYA  = 1
01151       KEYZ  = 1
01152 C
01153 C
01154       RETURN
01155       END
01156 
01157       SUBROUTINE T_GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
01158 C ----------------------------------------------------------------------
01159 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
01160 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
01161 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
01162 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
01163 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
01164 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
01165 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
01166 C
01167 C     called by : EVENTE, EVENTM, FUNTIH, .....
01168 C ----------------------------------------------------------------------
01169       IMPLICIT REAL*8(A-H,O-Z)
01170 C
01171       IF(IDFERM.EQ.0.OR.IABS(IDFERM).GT.4) GOTO 901
01172       IF(IABS(IHELIC).NE.1)                GOTO 901
01173       IH  =IHELIC
01174       IDTYPE =IABS(IDFERM)
01175       IC  =IDFERM/IDTYPE
01176       LEPQUA=INT(IDTYPE*0.4999999D0)
01177       IUPDOW=IDTYPE-2*LEPQUA-1
01178       CHARGE  =(-IUPDOW+2D0/3D0*LEPQUA)*IC
01179       SIZO3   =0.25D0*(IC-IH)*(1-2*IUPDOW)
01180       KOLOR=1+2*LEPQUA
01181 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
01182 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
01183       RETURN
01184  901  PRINT *,' STOP IN GIVIZO: WRONG PARAMS.'
01185       STOP
01186       END
01187       SUBROUTINE PHYFIX(NSTOP,NSTART)
01188       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) 
01189       SAVE /LUJETS/ 
01190 C NSTOP NSTART : when PHYTIA history ends and event starts.
01191       NSTOP=0
01192       NSTART=1
01193       DO I=1, N
01194        IF(K(I,1).NE.21) THEN
01195            NSTOP = I-1
01196            NSTART= I
01197            GOTO 500
01198        ENDIF
01199       ENDDO
01200  500  CONTINUE
01201       END
01202       SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
01203 C ----------------------------------------------------------------------
01204 C this subroutine fills one entry into the HEPEVT common
01205 C and updates the information for affected mother entries
01206 C
01207 C written by Martin W. Gruenewald (91/01/28)
01208 C
01209 C     called by : ZTOHEP,BTOHEP,DWLUxy
01210 C ----------------------------------------------------------------------
01211 C
01212 C this is the hepevt class in old style. No d_h_ class pre-name
01213 C this is the hepevt class in old style. No d_h_ class pre-name
01214       INTEGER NMXHEP
01215       PARAMETER (NMXHEP=4000)
01216       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01217       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01218      $        jdahep
01219       COMMON /hepevt/
01220      $      nevhep,               ! serial number
01221      $      nhep,                 ! number of particles
01222      $      isthep(nmxhep),   ! status code
01223      $      idhep(nmxhep),    ! particle ident KF
01224      $      jmohep(2,nmxhep), ! parent particles
01225      $      jdahep(2,nmxhep), ! childreen particles
01226      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01227      $      vhep(4,nmxhep)    ! vertex [mm]
01228 * ----------------------------------------------------------------------
01229       LOGICAL qedrad
01230       COMMON /phoqed/ 
01231      $     qedrad(nmxhep)    ! Photos flag
01232 * ----------------------------------------------------------------------
01233       SAVE hepevt,phoqed
01234 
01235 
01236       LOGICAL PHFLAG
01237 C
01238       REAL*4  P4(4)
01239 C
01240 C check address mode
01241       IF (N.EQ.0) THEN
01242 C
01243 C append mode
01244         IHEP=NHEP+1
01245       ELSE IF (N.GT.0) THEN
01246 C
01247 C absolute position
01248         IHEP=N
01249       ELSE
01250 C
01251 C relative position
01252         IHEP=NHEP+N
01253       END IF
01254 C
01255 C check on IHEP
01256       IF ((IHEP.LE.0).OR.(IHEP.GT.NMXHEP)) RETURN
01257 C
01258 C add entry
01259       NHEP=IHEP
01260       ISTHEP(IHEP)=IST
01261       IDHEP(IHEP)=ID
01262       JMOHEP(1,IHEP)=JMO1
01263       IF(JMO1.LT.0)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP
01264       JMOHEP(2,IHEP)=JMO2
01265       IF(JMO2.LT.0)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP
01266       JDAHEP(1,IHEP)=JDA1
01267       JDAHEP(2,IHEP)=JDA2
01268 C
01269       DO I=1,4
01270         PHEP(I,IHEP)=P4(I)
01271 C
01272 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
01273         VHEP(I,IHEP)=0.0
01274       END DO
01275       PHEP(5,IHEP)=PINV
01276 C FLAG FOR PHOTOS...
01277       QEDRAD(IHEP)=PHFLAG
01278 C
01279 C update process:
01280       DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
01281         IF(IP.GT.0)THEN
01282 C
01283 C if there is a daughter at IHEP, mother entry at IP has decayed
01284           IF(ISTHEP(IP).EQ.1)ISTHEP(IP)=2
01285 C
01286 C and daughter pointers of mother entry must be updated
01287           IF(JDAHEP(1,IP).EQ.0)THEN
01288             JDAHEP(1,IP)=IHEP
01289             JDAHEP(2,IP)=IHEP
01290           ELSE
01291             JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP))
01292           END IF
01293         END IF
01294       END DO
01295 C
01296       RETURN
01297       END
01298 
01299 
01300       FUNCTION IHEPDIM(DUM) 
01301 C this is the hepevt class in old style. No d_h_ class pre-name
01302 C this is the hepevt class in old style. No d_h_ class pre-name
01303       INTEGER NMXHEP
01304       PARAMETER (NMXHEP=4000)
01305       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01306       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01307      $        jdahep
01308       COMMON /hepevt/
01309      $      nevhep,               ! serial number
01310      $      nhep,                 ! number of particles
01311      $      isthep(nmxhep),   ! status code
01312      $      idhep(nmxhep),    ! particle ident KF
01313      $      jmohep(2,nmxhep), ! parent particles
01314      $      jdahep(2,nmxhep), ! childreen particles
01315      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01316      $      vhep(4,nmxhep)    ! vertex [mm]
01317 * ----------------------------------------------------------------------
01318       LOGICAL qedrad
01319       COMMON /phoqed/ 
01320      $     qedrad(nmxhep)    ! Photos flag
01321 * ----------------------------------------------------------------------
01322       SAVE hepevt,phoqed
01323 
01324 
01325       IHEPDIM=NHEP
01326       END
01327       FUNCTION ZPROP2(S)
01328       IMPLICIT REAL*8(A-H,O-Z)
01329       COMPLEX*16 CPRZ0,CPRZ0M
01330       AMZ=91.1882
01331       GAMMZ=2.49
01332       CPRZ0=DCMPLX((S-AMZ**2),S/AMZ*GAMMZ)
01333       CPRZ0M=1/CPRZ0
01334       ZPROP2=(ABS(CPRZ0M))**2
01335       END
01336 
01337       SUBROUTINE TAUPI0(MODE,JAK,ION)
01338 C no initialization required. Must be called once after every:
01339 C   1)    CALL DEKAY(1+10,...)
01340 C   2)    CALL DEKAY(2+10,...)
01341 C   3)    CALL DEXAY(1,...)
01342 C   4)    CALL DEXAY(2,...)
01343 C subroutine to decay originating from TAUOLA's taus: 
01344 C 1) etas (with CALL TAUETA(JAK))
01345 C 2) later pi0's from taus.
01346 C 3) extensions to other applications possible. 
01347 C this routine belongs to >tauola universal interface<, but uses 
01348 C routines from >tauola< utilities as well.  25.08.2005      
01349 C this is the hepevt class in old style. No d_h_ class pre-name
01350       INTEGER NMXHEP
01351       PARAMETER (NMXHEP=4000)
01352       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01353       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01354      $        jdahep
01355       COMMON /hepevt/
01356      $      nevhep,               ! serial number
01357      $      nhep,                 ! number of particles
01358      $      isthep(nmxhep),   ! status code
01359      $      idhep(nmxhep),    ! particle ident KF
01360      $      jmohep(2,nmxhep), ! parent particles
01361      $      jdahep(2,nmxhep), ! childreen particles
01362      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01363      $      vhep(4,nmxhep)    ! vertex [mm]
01364 * ----------------------------------------------------------------------
01365       LOGICAL qedrad
01366       COMMON /phoqed/ 
01367      $     qedrad(nmxhep)    ! Photos flag
01368 * ----------------------------------------------------------------------
01369       SAVE hepevt,phoqed
01370 
01371 
01372 
01373 C position of taus, must be defined by host program:
01374       COMMON /TAUPOS/ NP1,NP2
01375 c
01376       REAL  PHOT1(4),PHOT2(4)
01377       REAL*8  R,X(4),Y(4),PI0(4)
01378       INTEGER JEZELI(3),ION(3)
01379       DATA JEZELI /0,0,0/
01380       SAVE JEZELI
01381       IF (MODE.EQ.-1) THEN
01382         JEZELI(1)=ION(1)
01383         JEZELI(2)=ION(2)
01384         JEZELI(3)=ION(3)
01385         RETURN
01386       ENDIF
01387       IF (JEZELI(1).EQ.0) RETURN
01388       IF (JEZELI(2).EQ.1) CALL TAUETA(JAK)
01389       IF (JEZELI(3).EQ.1) CALL TAUK0S(JAK)
01390 C position of decaying particle:
01391       IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
01392         NPS=NP1
01393       ELSE
01394         NPS=NP2
01395       ENDIF
01396       nhepM=nhep                ! to avoid infinite loop
01397       DO K=JDAHEP(1,NPS),nhepM  ! we search for pi0's from tau till eor.
01398        IF (IDHEP(K).EQ.111.AND.JDAHEP(1,K).LE.K) THEN ! IF we found pi0
01399         DO L=1,4
01400           PI0(L)= phep(L,K)
01401         ENDDO
01402 ! random 3 vector on the sphere, masless
01403         R=SQRT(PI0(4)**2-PI0(3)**2-PI0(2)**2-PI0(1)**2)/2D0
01404         CALL SPHERD(R,X)
01405         X(4)=R
01406         Y(4)=R
01407         
01408         Y(1)=-X(1)
01409         Y(2)=-X(2)
01410         Y(3)=-X(3)
01411 ! boost to lab and to real*4
01412         CALL bostdq(-1,PI0,X,X)
01413         CALL bostdq(-1,PI0,Y,Y)
01414         DO L=1,4
01415          PHOT1(L)=X(L)
01416          PHOT2(L)=Y(L)
01417         ENDDO
01418 C to hepevt
01419         CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
01420         CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
01421        ENDIF
01422       ENDDO
01423 C
01424       END
01425       SUBROUTINE TAUETA(JAK)
01426 C subroutine to decay etas's from taus. 
01427 C this routine belongs to tauola universal interface, but uses 
01428 C routines from tauola utilities. Just flat phase space, but 4 channels.
01429 C it is called at the beginning of SUBR. TAUPI0(JAK)
01430 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005    
01431 C this is the hepevt class in old style. No d_h_ class pre-name
01432       INTEGER NMXHEP
01433       PARAMETER (NMXHEP=4000)
01434       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01435       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01436      $        jdahep
01437       COMMON /hepevt/
01438      $      nevhep,               ! serial number
01439      $      nhep,                 ! number of particles
01440      $      isthep(nmxhep),   ! status code
01441      $      idhep(nmxhep),    ! particle ident KF
01442      $      jmohep(2,nmxhep), ! parent particles
01443      $      jdahep(2,nmxhep), ! childreen particles
01444      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01445      $      vhep(4,nmxhep)    ! vertex [mm]
01446 * ----------------------------------------------------------------------
01447       LOGICAL qedrad
01448       COMMON /phoqed/ 
01449      $     qedrad(nmxhep)    ! Photos flag
01450 * ----------------------------------------------------------------------
01451       SAVE hepevt,phoqed
01452 
01453 
01454       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01455      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01456      *                 ,AMK,AMKZ,AMKST,GAMKST
01457 *
01458       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01459      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01460      *                 ,AMK,AMKZ,AMKST,GAMKST
01461 
01462 C position of taus, must be defined by host program:
01463       COMMON /TAUPOS/ NP1,NP2
01464 c
01465       REAL  RRR(1),BRSUM(3), RR(2)
01466       REAL  PHOT1(4),PHOT2(4),PHOT3(4)
01467       REAL*8    X(4),    Y(4),    Z(4)
01468       REAL                                YM1,YM2,YM3
01469       REAL*8  R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
01470       REAL*8 a,b,c
01471       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
01472 C position of decaying particle:
01473       IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
01474         NPS=NP1
01475       ELSE
01476         NPS=NP2
01477       ENDIF
01478       nhepM=nhep                ! to avoid infinite loop              
01479       DO K=JDAHEP(1,NPS),nhepM  ! we search for etas's from tau till eor.
01480        IF (IDHEP(K).EQ.221.AND.JDAHEP(1,K).LE.K) THEN ! IF we found eta
01481         DO L=1,4
01482           PETA(L)= phep(L,K)  ! eta 4 momentum
01483         ENDDO
01484 C       eta cumulated branching ratios:
01485         BRSUM(1)=0.389  ! gamma gamma
01486         BRSUM(2)=BRSUM(1)+0.319  ! 3 pi0
01487         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
01488         CALL RANMAR(RRR,1) 
01489         
01490         IF (RRR(1).LT.BRSUM(1)) THEN ! gamma gamma channel exactly like pi0
01491 ! random 3 vector on the sphere, masless   
01492          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
01493          CALL SPHERD(R,X) 
01494          X(4)=R
01495          Y(4)=R
01496         
01497          Y(1)=-X(1)
01498          Y(2)=-X(2)
01499          Y(3)=-X(3)
01500 ! boost to lab and to real*4
01501          CALL bostdq(-1,PETA,X,X)  
01502          CALL bostdq(-1,PETA,Y,Y)
01503          DO L=1,4
01504           PHOT1(L)=X(L)
01505           PHOT2(L)=Y(L)
01506          ENDDO
01507 C to hepevt
01508          CALL FILHEP(0,1,22,K,K,0,0,PHOT1,0.0,.TRUE.)
01509          CALL FILHEP(0,1,22,K,K,0,0,PHOT2,0.0,.TRUE.)
01510         ELSE ! 3 body channels
01511          IF(RRR(1).LT.BRSUM(2)) THEN  ! 3 pi0
01512           ID1= 111
01513           ID2= 111
01514           ID3= 111
01515           XM1=AMPIZ ! masses
01516           XM2=AMPIZ
01517           XM3=AMPIZ
01518          ELSEIF(RRR(1).LT.BRSUM(3)) THEN ! pi+ pi- pi0
01519           ID1= 211
01520           ID2=-211
01521           ID3= 111
01522           XM1=AMPI ! masses
01523           XM2=AMPI
01524           XM3=AMPIZ
01525          ELSE                            ! pi+ pi- gamma 
01526           ID1= 211
01527           ID2=-211
01528           ID3=  22
01529           XM1=AMPI ! masses
01530           XM2=AMPI
01531           XM3=0.0
01532          ENDIF
01533  7       CONTINUE  ! we generate mass of the first pair:
01534           CALL RANMAR(RR,2)
01535           R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)
01536           AMIN=XM1+XM2
01537           AMAX=R-XM3
01538           AM2=SQRT(AMIN**2+RR(1)*(AMAX**2-AMIN**2))
01539 C         weight for flat phase space
01540           WT=XLAM(R**2,AM2**2,XM3**2)*XLAM(AM2**2,XM1**2,XM2**2)
01541      &           /R**2                    /AM2**2
01542          IF (RR(2).GT.WT) GOTO 7
01543 
01544          RU=XLAM(AM2**2,XM1**2,XM2**2)/AM2/2  ! momenta of the
01545                                               ! first two products
01546                                               ! in the rest frame of that pair
01547          CALL SPHERD(RU,X)
01548          X(4)=SQRT(RU**2+XM1**2)
01549          Y(4)=SQRT(RU**2+XM2**2)
01550         
01551          Y(1)=-X(1)
01552          Y(2)=-X(2)
01553          Y(3)=-X(3)
01554 C generate momentum of that pair in rest frame of eta:
01555          RU=XLAM(R**2,AM2**2,XM3**2)/R/2
01556          CALL SPHERD(RU,Z)
01557          Z(4)=SQRT(RU**2+AM2**2)
01558 C and boost first two decay products to rest frame of eta.
01559          CALL bostdq(-1,Z,X,X)
01560          CALL bostdq(-1,Z,Y,Y)
01561 C redefine Z(4) to 4-momentum of the last decay product: 
01562          Z(1)=-Z(1)
01563          Z(2)=-Z(2)
01564          Z(3)=-Z(3)
01565          Z(4)=SQRT(RU**2+XM3**2)
01566 C boost all to lab and move to real*4; also masses
01567          CALL bostdq(-1,PETA,X,X)
01568          CALL bostdq(-1,PETA,Y,Y)
01569          CALL bostdq(-1,PETA,Z,Z)
01570          DO L=1,4
01571           PHOT1(L)=X(L)
01572           PHOT2(L)=Y(L)
01573           PHOT3(L)=Z(L)
01574          ENDDO
01575          YM1=XM1
01576          YM2=XM2
01577          YM3=XM3
01578 C to hepevt
01579          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
01580          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
01581          CALL FILHEP(0,1,ID3,K,K,0,0,PHOT3,YM3,.TRUE.)
01582         ENDIF
01583 
01584        ENDIF
01585       ENDDO
01586 C
01587       END
01588       SUBROUTINE TAUK0S(JAK)
01589 C subroutine to decay K0S's from taus. 
01590 C this routine belongs to tauola universal interface, but uses 
01591 C routines from tauola utilities. Just flat phase space, but 4 channels.
01592 C it is called at the beginning of SUBR. TAUPI0(JAK)
01593 C and as far as hepevt search it is basically the same as TAUPI0.  25.08.2005    
01594 C this is the hepevt class in old style. No d_h_ class pre-name
01595       INTEGER NMXHEP
01596       PARAMETER (NMXHEP=4000)
01597       REAL*8  phep,  vhep ! to be real*4/ *8  depending on host
01598       INTEGER nevhep,nhep,isthep,idhep,jmohep,
01599      $        jdahep
01600       COMMON /hepevt/
01601      $      nevhep,               ! serial number
01602      $      nhep,                 ! number of particles
01603      $      isthep(nmxhep),   ! status code
01604      $      idhep(nmxhep),    ! particle ident KF
01605      $      jmohep(2,nmxhep), ! parent particles
01606      $      jdahep(2,nmxhep), ! childreen particles
01607      $      phep(5,nmxhep),   ! four-momentum, mass [GeV]
01608      $      vhep(4,nmxhep)    ! vertex [mm]
01609 * ----------------------------------------------------------------------
01610       LOGICAL qedrad
01611       COMMON /phoqed/ 
01612      $     qedrad(nmxhep)    ! Photos flag
01613 * ----------------------------------------------------------------------
01614       SAVE hepevt,phoqed
01615 
01616 
01617 
01618       COMMON / PARMAS / AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01619      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01620      *                 ,AMK,AMKZ,AMKST,GAMKST
01621 *
01622       REAL*4            AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
01623      *                 ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
01624      *                 ,AMK,AMKZ,AMKST,GAMKST
01625 
01626 C position of taus, must be defined by host program:
01627       COMMON /TAUPOS/ NP1,NP2
01628 c
01629       REAL  RRR(1),BRSUM(3), RR(2)
01630       REAL  PHOT1(4),PHOT2(4),PHOT3(4)
01631       REAL*8    X(4),    Y(4),    Z(4)
01632       REAL                                YM1,YM2,YM3
01633       REAL*8  R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
01634       REAL*8 a,b,c
01635       XLAM(a,b,c)=SQRT(ABS((a-b-c)**2-4.0*b*c))
01636 C position of decaying particle:
01637       IF((KTO.EQ. 1).OR.(KTO.EQ.11)) THEN
01638         NPS=NP1
01639       ELSE
01640         NPS=NP2
01641       ENDIF
01642       nhepM=nhep                ! to avoid infinite loop              
01643       DO K=JDAHEP(1,NPS),nhepM  ! we search for K0S's from tau till eor.
01644        IF (IDHEP(K).EQ.310.AND.JDAHEP(1,K).LE.K) THEN ! IF we found K0S
01645 
01646       
01647         DO L=1,4
01648           PETA(L)= phep(L,K)  ! K0S 4 momentum  (this is cloned from eta decay)
01649         ENDDO
01650 C       K0S cumulated branching ratios:
01651         BRSUM(1)=0.313  ! 2 PI0
01652         BRSUM(2)=1.0 ! BRSUM(1)+0.319  ! Pi+ PI-
01653         BRSUM(3)=BRSUM(2)+0.237  ! pi+ pi- pi0 rest is thus pi+pi-gamma
01654         CALL RANMAR(RRR,1) 
01655 
01656          IF(RRR(1).LT.BRSUM(1)) THEN  ! 2 pi0
01657           ID1= 111
01658           ID2= 111
01659           XM1=AMPIZ ! masses
01660           XM2=AMPIZ
01661          ELSEIF(RRR(1).LT.BRSUM(2)) THEN ! pi+ pi- 
01662           ID1= 211
01663           ID2=-211
01664           XM1=AMPI ! masses
01665           XM2=AMPI
01666          ELSE                            ! gamma gamma unused !!!
01667           ID1= 22
01668           ID2= 22
01669           XM1= 0.0 ! masses
01670           XM2= 0.0
01671          ENDIF
01672         
01673 ! random 3 vector on the sphere, of equal mass !!  
01674          R=SQRT(PETA(4)**2-PETA(3)**2-PETA(2)**2-PETA(1)**2)/2D0
01675          R4=R
01676          R=SQRT(ABS(R**2-XM1**2))
01677          CALL SPHERD(R,X) 
01678          X(4)=R4
01679          Y(4)=R4
01680         
01681          Y(1)=-X(1)
01682          Y(2)=-X(2)
01683          Y(3)=-X(3)
01684 ! boost to lab and to real*4
01685          CALL bostdq(-1,PETA,X,X)  
01686          CALL bostdq(-1,PETA,Y,Y)
01687          DO L=1,4
01688           PHOT1(L)=X(L)
01689           PHOT2(L)=Y(L)
01690          ENDDO
01691 
01692          YM1=XM1
01693          YM2=XM2
01694 C to hepevt
01695          CALL FILHEP(0,1,ID1,K,K,0,0,PHOT1,YM1,.TRUE.)
01696          CALL FILHEP(0,1,ID2,K,K,0,0,PHOT2,YM2,.TRUE.)
01697 
01698 C
01699        ENDIF
01700       ENDDO
01701 
01702       END
Generated on Sun Oct 20 20:24:10 2013 for C++InterfacetoTauola by  doxygen 1.6.3