CDECK ID>, HWHBGF. *CMZ :- -26/04/91 11.11.55 by Bryan Webber *-- Author : Giovanni Abbiendi & Luca Stanco C----------------------------------------------------------------------- SUBROUTINE HWHBGF C----------------------------------------------------------------------- C Order Alpha_s processes in charged lepton-hadron collisions C C Process code IPROC has to be set in the Main Program C the following codes IPROC may be selected C C 9100 : NC BOSON-GLUON FUSION C 9100+IQK (IQK=1,...,6) : produced flavour is IQK C 9107 : produced J/psi + gluon C C 9110 : NC QCD COMPTON C 9110+IQK (IQK=1,...,12) : struck parton is IQK C C 9130 : NC order alpha_s processes (9100+9110) C C Select maximum and minimum generated flavour when IQK=0 C setting IFLMIN and IFLMAX in the Main Program C (allowed values from 1 to 6), default are 1 and 5 C allowing d,u,s,c,b,dbar,ubar,sbar,cbar,bbar C C CHARGED CURRENT Boson-Gluon Fusion processes C 9141 : CC s cbar (c sbar) C 9142 : CC b cbar (c bbar) C 9143 : CC s tbar (t cbar) C 9144 : CC b tbar (t bbar) C C other inputs : Q2MIN,Q2MAX,YBMIN,YBMAX,PTMIN,EMMIN,EMMAX C when IPROC=(1)9107 : as above but Q2WWMN, Q2WWMX substitute C Q2MIN and Q2MAX (EPA is used); ZJMAX cut C C Add 10000 to suppress soft remnant fragmentation C C Mean EVWGT = cross section in nanoBarn C C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP, & ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,FSIGMA(18), & SIGSUM,PROB,PRAN,PVRT(4),X INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,LEPFIN,ID1,ID2,I,IDD LOGICAL CHARGD,INCLUD(18),INSIDE(18) EXTERNAL HWR SAVE LEPFIN,ID1,ID2,FSIGMA,SIGSUM COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF, & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL, & IPROO,CHARGD,INCLUD,INSIDE C---Initialization IF (FSTWGT) THEN C---LEP = 1 FOR LEPTONS, -1 FOR ANTILEPTONS LEP=ZERO IF (IDHW(1).GE.121.AND.IDHW(1).LE.126) THEN LEP=ONE ELSEIF (IDHW(1).GE.127.AND.IDHW(1).LE.132) THEN LEP=-ONE ENDIF IF (LEP.EQ.ZERO) CALL HWWARN('HWHBGF',500,*999) IPROO=MOD(IPROC,100)/10 IF (IPROO.EQ.0.OR.IPROO.EQ.4) THEN IQK=MOD(IPROC,10) IFL=IQK IF (IQK.EQ.7) IFL=164 CHARGD=IPROO.EQ.4 ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN IQK=MOD(IPROC,100)-10 IFL=IQK+6 CHARGD=.FALSE. ELSEIF (IPROO.EQ.3) THEN IQK=0 IFL=0 CHARGD=.FALSE. ELSE CALL HWWARN('HWHBGF',501,*999) ENDIF C LEPFIN = IDHW(1) IF(CHARGD) THEN LEPFIN = IDHW(1)+1 IF (IQK.EQ.1) THEN IFLAVU=4 IFLAVD=3 ID1 = 3 ID2 = 10 ELSEIF (IQK.EQ.2) THEN IFLAVU=4 IFLAVD=5 ID1 = 5 ID2 = 10 ELSEIF (IQK.EQ.3) THEN IFLAVU=6 IFLAVD=3 ID1 = 3 ID2 =12 ELSE IFLAVU=6 IFLAVD=5 ID1 = 5 ID2 =12 ENDIF IF (LEP.EQ.-1.) THEN IDD=ID1 ID1=ID2-6 ID2=IDD+6 ENDIF ENDIF C IF (IQK.EQ.0) THEN DO I=1,18 INCLUD(I)=.TRUE. ENDDO IMIN=1 IMAX=18 DO I=1,6 IF (I.LT.IFLMIN.OR.I.GT.IFLMAX) INCLUD(I)=.FALSE. ENDDO DO I=7,18 IF (I.LE.12) THEN IF (I-6.LT.IFLMIN.OR.I-6.GT.IFLMAX) INCLUD(I)=.FALSE. ELSE IF (I-12.LT.IFLMIN.OR.I-12.GT.IFLMAX) INCLUD(I)=.FALSE. ENDIF ENDDO IF (IPROO.EQ.0) THEN DO I=7,18 INCLUD(I)=.FALSE. ENDDO IMIN=IFLMIN IMAX=IFLMAX ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN DO I=1,6 INCLUD(I)=.FALSE. ENDDO IMIN=IFLMIN+6 IMAX=IFLMAX+12 ELSEIF (IPROO.EQ.3) THEN IMIN=IFLMIN IMAX=IFLMAX+12 ENDIF ELSEIF (IQK.NE.0 .AND. (.NOT.CHARGD)) THEN DO I=1,18 INCLUD(I)=.FALSE. ENDDO IF (IFL.LE.18) THEN INCLUD(IFL)=.TRUE. IMIN=IFL IMAX=IFL ELSEIF (IFL.EQ.164) THEN INCLUD(7)=.TRUE. IMIN=7 IMAX=7 ENDIF ENDIF ENDIF C---End of initialization IF(GENEV) THEN IF (.NOT.CHARGD) THEN IF (IQK.EQ.0) THEN PRAN= SIGSUM * HWR() PROB=ZERO DO 10 IFL=IMIN,IMAX IF (.NOT.INSIDE(IFL)) GOTO 10 PROB=PROB+FSIGMA(IFL) IF (PROB.GE.PRAN) GOTO 20 10 CONTINUE ENDIF C---at this point the subprocess has been selected (IFL) 20 CONTINUE IF (IFL.LE.6) THEN C---Boson-Gluon Fusion event IDHW(NHEP+1)=IDHW(1) IDHW(NHEP+2)=13 IDHW(NHEP+3)=15 IDHW(NHEP+4)=LEPFIN IDHW(NHEP+5)=IFL IDHW(NHEP+6)=IFL+6 ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN C---QCD_Compton event IDHW(NHEP+1)=IDHW(1) IDHW(NHEP+2)=IFL-6 IDHW(NHEP+3)=15 IDHW(NHEP+4)=LEPFIN IDHW(NHEP+5)=IFL-6 IDHW(NHEP+6)=13 ELSEIF (IFL.EQ.164) THEN C---gamma+gluon-->J/Psi+gluon IDHW(NHEP+1)=IDHW(1) IDHW(NHEP+2)=13 IDHW(NHEP+3)=15 IDHW(NHEP+4)=LEPFIN IDHW(NHEP+5)=164 IDHW(NHEP+6)=13 ELSE CALL HWWARN('HWHBGF',503,*999) ENDIF ELSE C---Charged current event of specified flavours IDHW(NHEP+1)=IDHW(1) IDHW(NHEP+2)=13 IDHW(NHEP+3)=15 IDHW(NHEP+4)=LEPFIN IDHW(NHEP+5)=ID1 IDHW(NHEP+6)=ID2 ENDIF C DO 1 I=NHEP+1,NHEP+6 1 IDHEP(I)=IDPDG(IDHW(I)) C C---Codes common for all processes ISTHEP(NHEP+1)=111 ISTHEP(NHEP+2)=112 ISTHEP(NHEP+3)=110 ISTHEP(NHEP+4)=113 ISTHEP(NHEP+5)=114 ISTHEP(NHEP+6)=114 C DO I=NHEP+1,NHEP+6 JMOHEP(1,I)=NHEP+3 JDAHEP(1,I)=0 ENDDO C---Incoming lepton JMOHEP(2,NHEP+1)=NHEP+4 JDAHEP(2,NHEP+1)=NHEP+4 C---Hard Process C.M. JMOHEP(1,NHEP+3)=NHEP+1 JMOHEP(2,NHEP+3)=NHEP+2 JDAHEP(1,NHEP+3)=NHEP+4 JDAHEP(2,NHEP+3)=NHEP+6 C---Outgoing lepton JMOHEP(2,NHEP+4)=NHEP+1 JDAHEP(2,NHEP+4)=NHEP+1 C IF (IFL.LE.6 .OR. CHARGD) THEN C---Codes for boson-gluon fusion processes C--- Incoming gluon JMOHEP(2,NHEP+2)=NHEP+6 JDAHEP(2,NHEP+2)=NHEP+5 C--- Outgoing quark JMOHEP(2,NHEP+5)=NHEP+2 JDAHEP(2,NHEP+5)=NHEP+6 C--- Outgoing antiquark JMOHEP(2,NHEP+6)=NHEP+5 JDAHEP(2,NHEP+6)=NHEP+2 ELSEIF (IFL.GE.7 .AND. IFL.LE.12) THEN C---Codes for V+q --> q+g C--- Incoming quark JMOHEP(2,NHEP+2)=NHEP+5 JDAHEP(2,NHEP+2)=NHEP+6 C--- Outgoing quark JMOHEP(2,NHEP+5)=NHEP+6 JDAHEP(2,NHEP+5)=NHEP+2 C--- Outgoing gluon JMOHEP(2,NHEP+6)=NHEP+2 JDAHEP(2,NHEP+6)=NHEP+5 ELSEIF (IFL.GE.13 .AND. IFL.LE.18) THEN C---Codes for V+qbar --> qbar+g C--- Incoming antiquark JMOHEP(2,NHEP+2)=NHEP+6 JDAHEP(2,NHEP+2)=NHEP+5 C--- Outgoing antiquark JMOHEP(2,NHEP+5)=NHEP+2 JDAHEP(2,NHEP+5)=NHEP+6 C--- Outgoing gluon JMOHEP(2,NHEP+6)=NHEP+5 JDAHEP(2,NHEP+6)=NHEP+2 ELSEIF (IFL.EQ.164) THEN C---Codes for Gamma+gluon --> J/Psi+gluon C--- Incoming gluon JMOHEP(2,NHEP+2)=NHEP+6 JDAHEP(2,NHEP+2)=NHEP+6 C--- Outgoing J/Psi JMOHEP(2,NHEP+5)=NHEP+1 JDAHEP(2,NHEP+5)=NHEP+1 C--- Outgoing gluon JMOHEP(2,NHEP+6)=NHEP+2 JDAHEP(2,NHEP+6)=NHEP+2 ENDIF C---Computation of momenta in Laboratory frame of reference CALL HWHBKI NHEP=NHEP+6 C Decide which quark radiated and assign production vertices IF (IFL.LE.6) THEN C Boson-Gluon fusion case IF (1-Z.LT.HWR()) THEN C Gluon splitting to quark CALL HWVZRO(4,VHEP(1,NHEP-1)) CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT) CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP)) CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4)) ELSE C Gluon splitting to antiquark CALL HWVZRO(4,VHEP(1,NHEP)) CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP-1),PVRT) CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP-1)) CALL HWVEQU(4,VHEP(1,NHEP-1),VHEP(1,NHEP-4)) ENDIF ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN C QCD Compton case X=1/(1+SHAT/Q2) IF (1.LT.HWR()*(1+(1-X-Z)**2+6*X*(1-X)*Z*(1-Z))) THEN C Incoming quark radiated the gluon CALL HWVZRO(4,VHEP(1,NHEP-1)) CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT) CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP)) CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4)) ELSE C Outgoing quark radiated the gluon CALL HWVZRO(4,VHEP(1,NHEP-4)) CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP),PVRT) CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP)) CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1)) ENDIF ENDIF C---HERWIG gets confused if lepton momentum is different from beam C momentum, which it can be if incoming hadron has negative virtuality C As a temporary fix, simply copy the momentum. C Momentum conservation somehow gets taken care of HWBGEN! call hwvequ(5,phep(1,1),phep(1,nhep-5)) ELSE EVWGT=ZERO C---generation of the 5 variables Y,Q2,SHAT,Z,PHI and Jacobian computation C---in the largest phase space avalaible for selected processes and C---filling of logical vector INSIDE to tag contributing ones CALL HWHBRN (*999) C---calculate differential cross section corresponding to the chosen C---variables and the weight for MC generation IF (IQK.EQ.0) THEN C---many subprocesses included DO I=1,18 FSIGMA(I)=ZERO ENDDO SIGSUM=ZERO DO I=IMIN,IMAX IF (INSIDE(I)) THEN IFL=I DSIGMA=ZERO CALL HWHBSG FSIGMA(I)=DSIGMA SIGSUM=SIGSUM+DSIGMA ENDIF ENDDO EVWGT=SIGSUM * AJACOB ELSE C---only one subprocess included CALL HWHBSG EVWGT= DSIGMA * AJACOB ENDIF IF (EVWGT.LT.ZERO) EVWGT=ZERO ENDIF 999 END CDECK ID>, HWHBKI. *CMZ :- -26/04/91 13.19.32 by Federico Carminati *-- Author : Giovanni Abbiendi & Luca Stanco C---------------------------------------------------------------------- SUBROUTINE HWHBKI C---------------------------------------------------------------------- C gives the fourmomenta in the laboratory system for the particles C of the hard 2-->3 subprocess, to match with HERWIG routines of C jet evolution. C---------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWUECM,HWUPCM,HWUSQR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB, & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT, & PGAMMA(5),SG,MF1,MF2,EP,PP,EL,PL,E1,E2,Q1,COSBET,SINBET,COSTHE, & SINTHE,SINAZI,COSAZI,ROTAZI(3,3),EGAM,A,PPROT,MREMIN,PGAM,PEP(5), & COSPHI,SINPHI,ROT(3,3),EPROT,PROTON(5),MPART INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,I,IHAD,J,IS,ICMF LOGICAL CHARGD,INCLUD(18),INSIDE(18) EXTERNAL HWUECM,HWUPCM,HWUSQR COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF, & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL, & IPROO,CHARGD,INCLUD,INSIDE C IHAD=2 IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD) C---Set masses IF (CHARGD) THEN MPART=ZERO MF1=RMASS(IDHW(NHEP+5)) MF2=RMASS(IDHW(NHEP+6)) MREMIN=MP ELSE IS = IFL IF (IFL.EQ.164) IS=IQK MPART=ZERO IF (IFL.GE.7.AND.IFL.LE.18) MPART=RMASS(IFL-6) MF1=MFIN1(IS) MF2=MFIN2(IS) MREMIN = MREMIF(IS) ENDIF C---Calculation of kinematical variables for the generated event C in the center of mass frame of the incoming boson and parton C with parton along +z EGAM = HWUECM (SHAT, -Q2, MPART**2) PGAM = SQRT( EGAM**2 + Q2 ) EP = RSHAT-EGAM PP = PGAM A = (W2+Q2-MP**2)/TWO PPROT = (A*PGAM-EGAM*SQRT(A**2+MP**2*Q2))/Q2 IF (PPROT.LT.ZERO) CALL HWWARN('HWHBKI',101,*999) EPROT = SQRT(PPROT**2+MP**2) IF ((EPROT+PPROT).LT.(EP+PP)) CALL HWWARN('HWHBKI',102,*999) EL = ( PGAM / PPROT * SMA - Q2 ) / TWO + / (EGAM + PGAM / PPROT * EPROT) IF (EL.GT.ME) THEN PL = SQRT ( EL**2 - ME**2 ) ELSE CALL HWWARN ('HWHBKI',103,*999) ENDIF COSBET = (TWO * EPROT * EL - SMA) / (TWO * PPROT * PL) IF ( ABS(COSBET) .GE. ONE ) THEN COSBET = SIGN (ONE,COSBET) SINBET = ZERO ELSE SINBET = SQRT (ONE - COSBET**2) ENDIF SG = ME**2 + MPART**2 + Q2 + TWO * RSHAT * EL IF (SG.LE.(RSHAT+ML)**2 .OR. SG.GE.(RS-MREMIN)**2) + CALL HWWARN ('HWHBKI',104,*999) Q1 = HWUPCM( RSHAT, MF1, MF2) E1 = SQRT(Q1**2+MF1**2) E2 = SQRT(Q1**2+MF2**2) IF (Q1 .GT. ZERO) THEN COSTHE=(TWO*EP*E1 - Z*(SHAT+Q2))/(TWO*PP*Q1) IF (ABS(COSTHE) .GT. ONE) THEN COSTHE=SIGN(ONE,COSTHE) SINTHE=ZERO ELSE SINTHE=SQRT(ONE-COSTHE**2) ENDIF ELSE COSTHE=ZERO SINTHE=ONE ENDIF C---Initial lepton PHEP(1,NHEP+1)=PL*SINBET PHEP(2,NHEP+1)=ZERO PHEP(3,NHEP+1)=PL*COSBET PHEP(4,NHEP+1)=EL PHEP(5,NHEP+1)=RMASS(IDHW(1)) C---Initial Hadron PROTON(1)=ZERO PROTON(2)=ZERO PROTON(3)=PPROT PROTON(4)=EPROT CALL HWUMAS (PROTON) C---Initial parton PHEP(1,NHEP+2)=ZERO PHEP(2,NHEP+2)=ZERO PHEP(3,NHEP+2)=PP PHEP(4,NHEP+2)=EP PHEP(5,NHEP+2)=MPART C---HARD SUBPROCESS 2-->3 CENTRE OF MASS PHEP(1,NHEP+3)=PHEP(1,NHEP+1)+PHEP(1,NHEP+2) PHEP(2,NHEP+3)=PHEP(2,NHEP+1)+PHEP(2,NHEP+2) PHEP(3,NHEP+3)=PHEP(3,NHEP+1)+PHEP(3,NHEP+2) PHEP(4,NHEP+3)=PHEP(4,NHEP+1)+PHEP(4,NHEP+2) CALL HWUMAS ( PHEP(1,NHEP+3) ) C---Virtual boson PGAMMA(1)=ZERO PGAMMA(2)=ZERO PGAMMA(3)=-PGAM PGAMMA(4)=EGAM PGAMMA(5)=HWUSQR(Q2) C---Scattered lepton PHEP(1,NHEP+4)=PHEP(1,NHEP+1)-PGAMMA(1) PHEP(2,NHEP+4)=PHEP(2,NHEP+1)-PGAMMA(2) PHEP(3,NHEP+4)=PHEP(3,NHEP+1)-PGAMMA(3) PHEP(4,NHEP+4)=PHEP(4,NHEP+1)-PGAMMA(4) PHEP(5,NHEP+4)=RMASS(IDHW(1)) IF (CHARGD) PHEP(5,NHEP+4)=ZERO C---First Final parton: quark (or J/psi) in Boson-Gluon Fusion C--- quark or antiquark in QCD Compton PHEP(1,NHEP+5)=Q1*SINTHE*COS(PHI) PHEP(2,NHEP+5)=Q1*SINTHE*SIN(PHI) PHEP(3,NHEP+5)=Q1*COSTHE PHEP(4,NHEP+5)=E1 PHEP(5,NHEP+5)=MF1 C---Second Final parton: antiquark in Boson-Gluon Fusion C--- gluon in QCD Compton PHEP(1,NHEP+6)=-PHEP(1,NHEP+5) PHEP(2,NHEP+6)=-PHEP(2,NHEP+5) PHEP(3,NHEP+6)=-PHEP(3,NHEP+5) PHEP(4,NHEP+6)=E2 PHEP(5,NHEP+6)=MF2 C---Boost to lepton-hadron CM frame PEP(1) = PHEP(1,NHEP+1) PEP(2) = PHEP(2,NHEP+1) PEP(3) = PHEP(3,NHEP+1) + PPROT PEP(4) = PHEP(4,NHEP+1) + EPROT CALL HWUMAS (PEP) DO I=1,6 CALL HWULOF (PEP,PHEP(1,NHEP+I),PHEP(1,NHEP+I)) ENDDO CALL HWULOF (PEP,PROTON,PROTON) CALL HWULOF (PEP,PGAMMA,PGAMMA) C---Rotation around y-axis to align lepton beam with z-axis COSPHI = PHEP(3,NHEP+1) / & SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 ) SINPHI = PHEP(1,NHEP+1) / & SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 ) DO I=1,3 DO J=1,3 ROT(I,J)=ZERO ENDDO ENDDO ROT(1,1) = COSPHI ROT(1,3) = -SINPHI ROT(2,2) = ONE ROT(3,1) = SINPHI ROT(3,3) = COSPHI DO I=1,6 CALL HWUROF (ROT,PHEP(1,NHEP+I),PHEP(1,NHEP+I)) ENDDO CALL HWUROF (ROT,PROTON,PROTON) CALL HWUROF (ROT,PGAMMA,PGAMMA) C---Boost to the LAB frame ICMF=3 DO I=1,6 CALL HWULOB (PHEP(1,ICMF),PHEP(1,NHEP+I),PHEP(1,NHEP+I)) ENDDO CALL HWULOB (PHEP(1,ICMF),PROTON,PROTON) CALL HWULOB (PHEP(1,ICMF),PGAMMA,PGAMMA) C---Random azimuthal rotation CALL HWRAZM (ONE,COSAZI,SINAZI) DO I=1,3 DO J=1,3 ROTAZI(I,J)=ZERO ENDDO ENDDO ROTAZI(1,1) = COSAZI ROTAZI(1,2) = SINAZI ROTAZI(2,1) = -SINAZI ROTAZI(2,2) = COSAZI ROTAZI(3,3) = ONE DO I=1,6 CALL HWUROF (ROTAZI,PHEP(1,NHEP+I),PHEP(1,NHEP+I)) ENDDO CALL HWUROF (ROTAZI,PROTON,PROTON) CALL HWUROF (ROTAZI,PGAMMA,PGAMMA) 999 END CDECK ID>, HWHBRN. *CMZ :- -03/07/95 19.02.12 by Giovanni Abbiendi *-- Author : Giovanni Abbiendi & Luca Stanco C----------------------------------------------------------------------- SUBROUTINE HWHBRN (*) C---------------------------------------------------------------------- C Returns a point in the phase space (Y,Q2,SHAT,Z,PHI) and the C corresponding Jacobian factor AJACOB C Fill the logical vector INSIDE to tag contributing subprocesses C to the cross-section C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWRUNI,HWR,HWUPCM,LEP,Y,Q2,SHAT,Z,PHI,AJACOB, & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT, & MF1,MF2,YMIN,YMAX,YJAC,Q2INF,Q2SUP,Q2JAC,EMW2,ZMIN,ZMAX,ZJAC, & GAMMA2,LAMBDA,PHIJAC,ZINT,ZLMIN,ZL,EMW,TMIN,TMAX,EMLMIN,EMLMAX, & SHMIN,EMMIF(18),EMMAF(18),WMIF(18),WMIN,MREMIN,YMIF(18),Q1CM(18), & Q2MAF(18),EMMAWF(18),ZMIF(18),ZMAF(18),PLMAX,PINC,SHINF,SHSUP, & SHJAC,CTHLIM,Q1,DETDSH,SRY,SRY0,SRY1 INTEGER IQK,IFLAVU,IFLAVD,I,IMIN,IMAX,IFL,IPROO,IHAD,NTRY,DEBUG LOGICAL CHARGD,INCLUD(18),INSIDE(18) EXTERNAL HWRUNI,HWR,HWUPCM SAVE EMLMIN,EMLMAX,EMMIF,EMMAF,MREMIN,MF1,MF2,YMIF, & YMIN,YMAX,WMIN,WMIF COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF, & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL, & IPROO,CHARGD,INCLUD,INSIDE EQUIVALENCE (EMW,RMASS(198)) C IHAD=2 IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD) C---Initialization IF (FSTWGT.OR.IHAD.NE.2) THEN ME = RMASS(IDHW(1)) MP = RMASS(IDHW(IHAD)) RS = PHEP(5,3) SMA = RS**2-ME**2-MP**2 PINC = HWUPCM(RS,ME,MP) C---Charged current IF (CHARGD) THEN ML=RMASS(IDHW(1)+1) YMAX = ONE - TWO*ML*MP / SMA YMAX = MIN(YMAX,YBMAX) MREMIN=MP IF (LEP.EQ.ONE) THEN MF1=RMASS(IFLAVD) MF2=RMASS(IFLAVU) ELSE MF1=RMASS(IFLAVU) MF2=RMASS(IFLAVD) ENDIF SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 + + TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2) EMLMIN=MAX(EMMIN,SQRT(SHMIN)) EMLMAX=MIN(EMMAX,RS-ML-MREMIN) DEBUG=1 IF (EMLMIN.GT.EMLMAX) GOTO 888 WMIN=EMLMIN+MREMIN PLMAX=HWUPCM(RS,ML,WMIN) YMIN = ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+ + PINC*PLMAX)/SMA YMIN = MAX(YMIN,YBMIN) DEBUG=2 IF (YMIN.GT.YMAX) GOTO 888 ELSE C---Neutral current ML = ME YMAX = ONE - TWO*ML*MP / SMA YMAX = MIN(YMAX,YBMAX) DO I=1,18 YMIF(I)=ZERO EMMIF(I)=ZERO EMMAF(I)=ZERO WMIF(I)=ZERO IF (I.LE.8) THEN C---Boson-Gluon Fusion (also J/Psi) and QCD Compton with struck u or d MREMIF(I)=MP IF (I.LE.6) THEN MFIN1(I)=RMASS(I) MFIN2(I)=RMASS(I+6) ELSE MFIN1(I)=RMASS(I-6) MFIN2(I)=ZERO ENDIF ELSE C---QCD Compton with struck non-valence parton MREMIF(I)=MP+RMASS(I-6) MFIN1(I)=RMASS(I-6) MFIN2(I)=ZERO ENDIF ENDDO IF (IFL.EQ.164) THEN C---J/Psi MFIN1(7)=RMASS(164) MFIN2(7)=ZERO ENDIF C---y boundaries for different flavours and processes DO 100 I=IMIN,IMAX IF (INCLUD(I)) THEN MF1=MFIN1(I) MF2=MFIN2(I) MREMIN=MREMIF(I) SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 + + TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2) EMMIF(I) = MAX(EMMIN,SQRT(SHMIN)) EMMAF(I) = MIN(EMMAX,RS-ML-MREMIN) IF (EMMIF(I).GT.EMMAF(I)) THEN INCLUD(I)=.FALSE. CALL HWWARN('HWHBRN',3,*999) GOTO 100 ENDIF WMIF(I) = EMMIF(I)+MREMIF(I) WMIN = WMIF(I) PLMAX = HWUPCM(RS,ML,WMIN) YMIF(I)=ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+ + PINC*PLMAX)/SMA IF (YMIF(I).GT.YMAX) THEN INCLUD(I)=.FALSE. CALL HWWARN('HWHBRN',4,*999) GOTO 100 ENDIF ENDIF 100 CONTINUE C---considering the largest boundaries EMLMIN=EMMIF(IMIN) EMLMAX=EMMAF(IMIN) IF (IPROO.EQ.3) THEN EMLMIN=MIN(EMMIF(IMIN),EMMIF(IMIN+6)) EMLMAX=MAX(EMMAF(IMIN),EMMAF(IMIN+6)) ENDIF DEBUG=3 IF (EMLMIN.GT.EMLMAX) GOTO 888 YMIN=YMIF(IMIN) IF (IPROO.EQ.3) YMIN=MIN(YMIF(IMIN),YMIF(IMIN+6)) YMIN = MAX(YMIN,YBMIN) DEBUG=4 IF (YMIN.GT.YMAX) GOTO 888 WMIN = WMIF(IMIN) MREMIN = MREMIF(IMIN) MF1=MFIN1(IMIN) MF2=MFIN2(IMIN) IF (IPROO.EQ.3) THEN WMIN = MIN(WMIF(IMIN),WMIF(IMIN+6)) MREMIN = MIN(MREMIF(IMIN),MREMIF(IMIN+6)) ENDIF ENDIF ENDIF C---Random generation in largest phase space Y=ZERO Q2=ZERO SHAT=ZERO Z=ZERO PHI=ZERO AJACOB=ZERO C---y generation IF (.NOT.CHARGD) THEN IF (IFL.LE.5.OR.(IFL.GE.7.AND.IFL.LE.18)) THEN SRY0 = SQRT(YMIN) SRY1 = SQRT(YMAX) SRY = HWRUNI(0,SRY0,SRY1) Y = SRY**2 YJAC = TWO*SRY*(SRY1-SRY0) ELSEIF (IFL.EQ.6) THEN Y = SQRT(HWRUNI(0,YMIN**2,YMAX**2)) YJAC = HALF * (YMAX**2-YMIN**2) / Y ELSEIF (IFL.EQ.164) THEN C---in J/psi photoproduction Y and Q2 are given by the Equivalent Photon C Approximation 10 NTRY=0 20 NTRY=NTRY+1 IF (NTRY.GT.NETRY) CALL HWWARN('HWHBRN',50,*10) Y = (YMIN/YMAX)**HWR()*YMAX IF (ONE+(ONE-Y)**2.LT.TWO*HWR()) GOTO 20 YJAC=(TWO*LOG(YMAX/YMIN)-TWO*(YMAX-YMIN) & +HALF*(YMAX**2-YMIN**2)) ENDIF ELSE IF (IPRO.EQ.5) THEN Y = EXP(HWRUNI(0,LOG(YMIN),LOG(YMAX))) YJAC = Y * LOG(YMAX/YMIN) ELSE Y = HWRUNI(0,YMIN,YMAX) YJAC = YMAX - YMIN ENDIF ENDIF C---Q**2 generation Q2INF = ME**2*Y**2 / (ONE-Y) Q2SUP = MP**2 + SMA*Y - WMIN**2 IF (IFL.EQ.164) THEN Q2INF = MAX(Q2INF,Q2WWMN) Q2SUP = MIN(Q2SUP,Q2WWMX) ELSE Q2INF = MAX(Q2INF,Q2MIN) Q2SUP = MIN(Q2SUP,Q2MAX) ENDIF DEBUG=5 IF (Q2INF .GT. Q2SUP) GOTO 888 C IF (.NOT.CHARGD) THEN IF (IFL.EQ.164) THEN Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP))) Q2JAC = LOG(Q2SUP/Q2INF) ELSEIF (Q2INF.LT.RMASS(4)**2) THEN Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP))) Q2JAC = Q2 * LOG(Q2SUP/Q2INF) ELSE Q2 = Q2INF*Q2SUP/HWRUNI(0,Q2INF,Q2SUP) Q2JAC = Q2**2 * (Q2SUP-Q2INF)/(Q2SUP*Q2INF) ENDIF ELSE EMW2=EMW**2 Q2=(Q2INF+EMW2)*(Q2SUP+EMW2)/(HWRUNI(0,Q2INF,Q2SUP)+EMW2)-EMW2 Q2JAC=(Q2+EMW2)**2*(Q2SUP-Q2INF)/((Q2SUP+EMW2)*(Q2INF+EMW2)) ENDIF W2 = MP**2 + SMA*Y - Q2 C---s_hat generation SHINF = EMLMIN **2 SHSUP = (MIN(SQRT(W2)-MREMIN,EMLMAX))**2 DEBUG=6 IF (SHINF .GT. SHSUP) GOTO 888 C IF (IPRO.EQ.91) THEN IF (.NOT.CHARGD) THEN SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP) SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF) ELSE SHAT = EXP(HWRUNI(0,LOG(SHINF),LOG(SHSUP))) SHJAC = SHAT*(LOG(SHSUP/SHINF)) ENDIF ELSE EMW2=EMW**2 IF (SHINF.GT.EMW2+10*GAMW*EMW) THEN SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP) SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF) ELSEIF (SHSUP.LT.EMW2-10*EMW*GAMW) THEN SHAT = HWRUNI(0,SHINF,SHSUP) SHJAC = SHSUP-SHINF ELSE TMIN=ATAN((SHINF-EMW2)/(GAMW*EMW)) TMAX=ATAN((SHSUP-EMW2)/(GAMW*EMW)) SHAT = GAMW*EMW*TAN(HWRUNI(0,TMIN,TMAX))+EMW2 SHJAC=((SHAT-EMW2)**2+(GAMW*EMW)**2)/(GAMW*EMW)*(TMAX-TMIN) ENDIF ENDIF DETDSH = ONE/SMA/Y SHJAC=SHJAC*DETDSH RSHAT = SQRT (SHAT) C--- z generation ZMIN = 10E10 ZMAX = -ONE IF (.NOT.CHARGD) THEN DO I=1,18 Q1CM(I) = ZERO ZMIF(I) = ZERO ZMAF(I) = ZERO ENDDO DO 150 I=IMIN,IMAX IF (INCLUD(I)) THEN Q1CM(I) = HWUPCM( RSHAT, MFIN1(I), MFIN2(I) ) IF (Q1CM(I) .LT. PTMIN) THEN ZMAF(I)=-ONE GOTO 150 ENDIF CTHLIM = SQRT(ONE - (PTMIN / Q1CM(I))**2) GAMMA2 = SHAT + MFIN1(I)**2 - MFIN2(I)**2 LAMBDA = (SHAT-MFIN1(I)**2-MFIN2(I)**2)**2 - + 4.D0*MFIN1(I)**2*MFIN2(I)**2 ZMIF(I) = (GAMMA2 - SQRT(LAMBDA)*CTHLIM)/TWO/SHAT ZMIF(I) = MAX(ZMIF(I),ZERO) ZMAF(I) = (GAMMA2 + SQRT(LAMBDA)*CTHLIM)/TWO/SHAT ZMAF(I) = MIN(ZMAF(I),ONE) ZMIN = MIN( ZMIN, ZMIF(I) ) ZMAX = MAX( ZMAX, ZMAF(I) ) ENDIF 150 CONTINUE IF (IFL.EQ.164) ZMAX=MIN(ZMAX,ZJMAX) ELSE Q1 = HWUPCM(RSHAT,MF1,MF2) DEBUG=7 IF (Q1.LT.PTMIN) GOTO 888 CTHLIM = SQRT(ONE-(PTMIN/Q1)**2) GAMMA2 = SHAT+MF1**2-MF2**2 LAMBDA = (SHAT-MF1**2-MF2**2)**2-4.D0*MF1**2*MF2**2 ZMIN = (GAMMA2-SQRT(LAMBDA)*CTHLIM)/TWO/SHAT ZMIN = MAX(ZMIN,1D-6) ZMAX = (GAMMA2+SQRT(LAMBDA)*CTHLIM)/TWO/SHAT ZMAX = MIN(ZMAX,ONE-1D-6) ENDIF DEBUG=8 IF (ZMIN .GT. ZMAX) GOTO 888 ZLMIN = LOG(ZMIN/(ONE-ZMIN)) ZINT = LOG(ZMAX/(ONE-ZMAX)) - LOG(ZMIN/(ONE-ZMIN)) ZL = ZLMIN+HWR()*ZINT Z = EXP(ZL)/(ONE+EXP(ZL)) ZJAC = Z*(ONE-Z)*ZINT C DEBUG=9 IF ((Y.LT.YMIN.OR.Y.GT.YMAX).OR.(Q2.LT.Q2INF.OR.Q2.GT.Q2SUP).OR. + (SHAT.LT.SHINF.OR.SHAT.GT.SHSUP).OR.(Z.LT.ZMIN.OR.Z.GT.ZMAX)) + GOTO 888 C---Phi generation PHI = HWRUNI(0,ZERO,2*PIFAC) PHIJAC = 2 * PIFAC IF (IFL.EQ.164) PHIJAC=ONE C AJACOB = YJAC * Q2JAC * SHJAC * ZJAC * PHIJAC C IF (IQK.NE.0.OR.IPRO.EQ.5) GOTO 999 C---contributing subprocesses: filling of logical vector INSIDE DO I=1,18 INSIDE(I)=.FALSE. Q2MAF(I)=ZERO EMMAWF(I)=ZERO ENDDO DO 200 I=IMIN,IMAX IF (INCLUD(I)) THEN IF ( Y.LT.YMIF(I) ) GOTO 200 C Q2MAF(I) = MP**2 + SMA*Y - WMIF(I)**2 Q2MAF(I) = MIN( Q2MAF(I), Q2MAX) IF (Q2INF .GT. Q2MAF(I)) GOTO 200 IF (Q2.LT.Q2INF .OR. Q2.GT.Q2MAF(I)) GOTO 200 C EMMAWF(I) = SQRT(W2) - MREMIF(I) EMMAWF(I) = MIN( EMMAWF(I), EMLMAX ) C IF (EMMIF(I) .GT. EMMAWF(I)) GOTO 200 IF (SHAT.LT.EMMIF(I)**2.OR.SHAT.GT.EMMAWF(I)**2) GOTO 200 C IF (ZMIF(I) .GT. ZMAF(I)) GOTO 200 IF (Z.LT.ZMIF(I) .OR. Z.GT.ZMAF(I)) GOTO 200 INSIDE(I)=.TRUE. ENDIF 200 CONTINUE 999 RETURN 888 EVWGT=ZERO C---UNCOMMENT THIS LINE TO GET A DEBUGGING WARNING FOR NO PHASE-SPACE C CALL HWWARN('HWHBRN',DEBUG,*777) 777 RETURN 1 END