4 *CMZ :- -26/04/91 11.11.55 by Bryan Webber
6 *-- Author : Giovanni Abbiendi & Luca Stanco
8 C-----------------------------------------------------------------------
12 C-----------------------------------------------------------------------
14 C Order Alpha_s processes in charged lepton-hadron collisions
18 C Process code IPROC has to be set in the Main Program
20 C the following codes IPROC may be selected
24 C 9100 : NC BOSON-GLUON FUSION
26 C 9100+IQK (IQK=1,...,6) : produced flavour is IQK
28 C 9107 : produced J/psi + gluon
32 C 9110 : NC QCD COMPTON
34 C 9110+IQK (IQK=1,...,12) : struck parton is IQK
38 C 9130 : NC order alpha_s processes (9100+9110)
42 C Select maximum and minimum generated flavour when IQK=0
44 C setting IFLMIN and IFLMAX in the Main Program
46 C (allowed values from 1 to 6), default are 1 and 5
48 C allowing d,u,s,c,b,dbar,ubar,sbar,cbar,bbar
52 C CHARGED CURRENT Boson-Gluon Fusion processes
54 C 9141 : CC s cbar (c sbar)
56 C 9142 : CC b cbar (c bbar)
58 C 9143 : CC s tbar (t cbar)
60 C 9144 : CC b tbar (t bbar)
64 C other inputs : Q2MIN,Q2MAX,YBMIN,YBMAX,PTMIN,EMMIN,EMMAX
66 C when IPROC=(1)9107 : as above but Q2WWMN, Q2WWMX substitute
68 C Q2MIN and Q2MAX (EPA is used); ZJMAX cut
72 C Add 10000 to suppress soft remnant fragmentation
76 C Mean EVWGT = cross section in nanoBarn
80 C-----------------------------------------------------------------------
82 INCLUDE 'HERWIG61.INC'
84 DOUBLE PRECISION HWR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,
86 & ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,FSIGMA(18),
88 & SIGSUM,PROB,PRAN,PVRT(4),X
90 INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,LEPFIN,ID1,ID2,I,IDD
92 LOGICAL CHARGD,INCLUD(18),INSIDE(18)
96 SAVE LEPFIN,ID1,ID2,FSIGMA,SIGSUM
98 COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
100 & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
102 & IPROO,CHARGD,INCLUD,INSIDE
108 C---LEP = 1 FOR LEPTONS, -1 FOR ANTILEPTONS
112 IF (IDHW(1).GE.121.AND.IDHW(1).LE.126) THEN
116 ELSEIF (IDHW(1).GE.127.AND.IDHW(1).LE.132) THEN
122 IF (LEP.EQ.ZERO) CALL HWWARN('HWHBGF',500,*999)
124 IPROO=MOD(IPROC,100)/10
126 IF (IPROO.EQ.0.OR.IPROO.EQ.4) THEN
132 IF (IQK.EQ.7) IFL=164
136 ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN
138 IQK=MOD(IPROC,100)-10
144 ELSEIF (IPROO.EQ.3) THEN
154 CALL HWWARN('HWHBGF',501,*999)
176 ELSEIF (IQK.EQ.2) THEN
186 ELSEIF (IQK.EQ.3) THEN
236 IF (I.LT.IFLMIN.OR.I.GT.IFLMAX) INCLUD(I)=.FALSE.
244 IF (I-6.LT.IFLMIN.OR.I-6.GT.IFLMAX) INCLUD(I)=.FALSE.
248 IF (I-12.LT.IFLMIN.OR.I-12.GT.IFLMAX) INCLUD(I)=.FALSE.
266 ELSEIF (IPROO.EQ.1.OR.IPROO.EQ.2) THEN
278 ELSEIF (IPROO.EQ.3) THEN
286 ELSEIF (IQK.NE.0 .AND. (.NOT.CHARGD)) THEN
302 ELSEIF (IFL.EQ.164) THEN
316 C---End of initialization
320 IF (.NOT.CHARGD) THEN
330 IF (.NOT.INSIDE(IFL)) GOTO 10
332 PROB=PROB+FSIGMA(IFL)
334 IF (PROB.GE.PRAN) GOTO 20
340 C---at this point the subprocess has been selected (IFL)
346 C---Boson-Gluon Fusion event
360 ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN
362 C---QCD_Compton event
376 ELSEIF (IFL.EQ.164) THEN
378 C---gamma+gluon-->J/Psi+gluon
394 CALL HWWARN('HWHBGF',503,*999)
400 C---Charged current event of specified flavours
420 1 IDHEP(I)=IDPDG(IDHW(I))
424 C---Codes common for all processes
450 JMOHEP(2,NHEP+1)=NHEP+4
452 JDAHEP(2,NHEP+1)=NHEP+4
454 C---Hard Process C.M.
456 JMOHEP(1,NHEP+3)=NHEP+1
458 JMOHEP(2,NHEP+3)=NHEP+2
460 JDAHEP(1,NHEP+3)=NHEP+4
462 JDAHEP(2,NHEP+3)=NHEP+6
466 JMOHEP(2,NHEP+4)=NHEP+1
468 JDAHEP(2,NHEP+4)=NHEP+1
472 IF (IFL.LE.6 .OR. CHARGD) THEN
474 C---Codes for boson-gluon fusion processes
478 JMOHEP(2,NHEP+2)=NHEP+6
480 JDAHEP(2,NHEP+2)=NHEP+5
484 JMOHEP(2,NHEP+5)=NHEP+2
486 JDAHEP(2,NHEP+5)=NHEP+6
488 C--- Outgoing antiquark
490 JMOHEP(2,NHEP+6)=NHEP+5
492 JDAHEP(2,NHEP+6)=NHEP+2
494 ELSEIF (IFL.GE.7 .AND. IFL.LE.12) THEN
496 C---Codes for V+q --> q+g
500 JMOHEP(2,NHEP+2)=NHEP+5
502 JDAHEP(2,NHEP+2)=NHEP+6
506 JMOHEP(2,NHEP+5)=NHEP+6
508 JDAHEP(2,NHEP+5)=NHEP+2
512 JMOHEP(2,NHEP+6)=NHEP+2
514 JDAHEP(2,NHEP+6)=NHEP+5
516 ELSEIF (IFL.GE.13 .AND. IFL.LE.18) THEN
518 C---Codes for V+qbar --> qbar+g
520 C--- Incoming antiquark
522 JMOHEP(2,NHEP+2)=NHEP+6
524 JDAHEP(2,NHEP+2)=NHEP+5
526 C--- Outgoing antiquark
528 JMOHEP(2,NHEP+5)=NHEP+2
530 JDAHEP(2,NHEP+5)=NHEP+6
534 JMOHEP(2,NHEP+6)=NHEP+5
536 JDAHEP(2,NHEP+6)=NHEP+2
538 ELSEIF (IFL.EQ.164) THEN
540 C---Codes for Gamma+gluon --> J/Psi+gluon
544 JMOHEP(2,NHEP+2)=NHEP+6
546 JDAHEP(2,NHEP+2)=NHEP+6
550 JMOHEP(2,NHEP+5)=NHEP+1
552 JDAHEP(2,NHEP+5)=NHEP+1
556 JMOHEP(2,NHEP+6)=NHEP+2
558 JDAHEP(2,NHEP+6)=NHEP+2
562 C---Computation of momenta in Laboratory frame of reference
568 C Decide which quark radiated and assign production vertices
572 C Boson-Gluon fusion case
574 IF (1-Z.LT.HWR()) THEN
576 C Gluon splitting to quark
578 CALL HWVZRO(4,VHEP(1,NHEP-1))
580 CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT)
582 CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP))
584 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4))
588 C Gluon splitting to antiquark
590 CALL HWVZRO(4,VHEP(1,NHEP))
592 CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP-1),PVRT)
594 CALL HWUDKL(IFL,PVRT,VHEP(1,NHEP-1))
596 CALL HWVEQU(4,VHEP(1,NHEP-1),VHEP(1,NHEP-4))
600 ELSEIF (IFL.GE.7.AND.IFL.LE.18) THEN
606 IF (1.LT.HWR()*(1+(1-X-Z)**2+6*X*(1-X)*Z*(1-Z))) THEN
608 C Incoming quark radiated the gluon
610 CALL HWVZRO(4,VHEP(1,NHEP-1))
612 CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT)
614 CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP))
616 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4))
620 C Outgoing quark radiated the gluon
622 CALL HWVZRO(4,VHEP(1,NHEP-4))
624 CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP),PVRT)
626 CALL HWUDKL(IFL-6,PVRT,VHEP(1,NHEP))
628 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
634 C---HERWIG gets confused if lepton momentum is different from beam
636 C momentum, which it can be if incoming hadron has negative virtuality
638 C As a temporary fix, simply copy the momentum.
640 C Momentum conservation somehow gets taken care of HWBGEN!
642 call hwvequ(5,phep(1,1),phep(1,nhep-5))
648 C---generation of the 5 variables Y,Q2,SHAT,Z,PHI and Jacobian computation
650 C---in the largest phase space avalaible for selected processes and
652 C---filling of logical vector INSIDE to tag contributing ones
656 C---calculate differential cross section corresponding to the chosen
658 C---variables and the weight for MC generation
662 C---many subprocesses included
690 EVWGT=SIGSUM * AJACOB
694 C---only one subprocess included
698 EVWGT= DSIGMA * AJACOB
702 IF (EVWGT.LT.ZERO) EVWGT=ZERO
710 *CMZ :- -26/04/91 13.19.32 by Federico Carminati
712 *-- Author : Giovanni Abbiendi & Luca Stanco
714 C----------------------------------------------------------------------
718 C----------------------------------------------------------------------
720 C gives the fourmomenta in the laboratory system for the particles
722 C of the hard 2-->3 subprocess, to match with HERWIG routines of
726 C----------------------------------------------------------------------
728 INCLUDE 'HERWIG61.INC'
730 DOUBLE PRECISION HWUECM,HWUPCM,HWUSQR,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,
732 & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,
734 & PGAMMA(5),SG,MF1,MF2,EP,PP,EL,PL,E1,E2,Q1,COSBET,SINBET,COSTHE,
736 & SINTHE,SINAZI,COSAZI,ROTAZI(3,3),EGAM,A,PPROT,MREMIN,PGAM,PEP(5),
738 & COSPHI,SINPHI,ROT(3,3),EPROT,PROTON(5),MPART
740 INTEGER IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,IPROO,I,IHAD,J,IS,ICMF
742 LOGICAL CHARGD,INCLUD(18),INSIDE(18)
744 EXTERNAL HWUECM,HWUPCM,HWUSQR
746 COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
748 & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
750 & IPROO,CHARGD,INCLUD,INSIDE
756 IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
764 MF1=RMASS(IDHW(NHEP+5))
766 MF2=RMASS(IDHW(NHEP+6))
774 IF (IFL.EQ.164) IS=IQK
778 IF (IFL.GE.7.AND.IFL.LE.18) MPART=RMASS(IFL-6)
788 C---Calculation of kinematical variables for the generated event
790 C in the center of mass frame of the incoming boson and parton
792 C with parton along +z
794 EGAM = HWUECM (SHAT, -Q2, MPART**2)
796 PGAM = SQRT( EGAM**2 + Q2 )
802 A = (W2+Q2-MP**2)/TWO
804 PPROT = (A*PGAM-EGAM*SQRT(A**2+MP**2*Q2))/Q2
806 IF (PPROT.LT.ZERO) CALL HWWARN('HWHBKI',101,*999)
808 EPROT = SQRT(PPROT**2+MP**2)
810 IF ((EPROT+PPROT).LT.(EP+PP)) CALL HWWARN('HWHBKI',102,*999)
812 EL = ( PGAM / PPROT * SMA - Q2 ) / TWO
814 + / (EGAM + PGAM / PPROT * EPROT)
818 PL = SQRT ( EL**2 - ME**2 )
822 CALL HWWARN ('HWHBKI',103,*999)
826 COSBET = (TWO * EPROT * EL - SMA) / (TWO * PPROT * PL)
828 IF ( ABS(COSBET) .GE. ONE ) THEN
830 COSBET = SIGN (ONE,COSBET)
836 SINBET = SQRT (ONE - COSBET**2)
840 SG = ME**2 + MPART**2 + Q2 + TWO * RSHAT * EL
842 IF (SG.LE.(RSHAT+ML)**2 .OR. SG.GE.(RS-MREMIN)**2)
844 + CALL HWWARN ('HWHBKI',104,*999)
846 Q1 = HWUPCM( RSHAT, MF1, MF2)
848 E1 = SQRT(Q1**2+MF1**2)
850 E2 = SQRT(Q1**2+MF2**2)
852 IF (Q1 .GT. ZERO) THEN
854 COSTHE=(TWO*EP*E1 - Z*(SHAT+Q2))/(TWO*PP*Q1)
856 IF (ABS(COSTHE) .GT. ONE) THEN
858 COSTHE=SIGN(ONE,COSTHE)
864 SINTHE=SQRT(ONE-COSTHE**2)
878 PHEP(1,NHEP+1)=PL*SINBET
882 PHEP(3,NHEP+1)=PL*COSBET
886 PHEP(5,NHEP+1)=RMASS(IDHW(1))
912 C---HARD SUBPROCESS 2-->3 CENTRE OF MASS
914 PHEP(1,NHEP+3)=PHEP(1,NHEP+1)+PHEP(1,NHEP+2)
916 PHEP(2,NHEP+3)=PHEP(2,NHEP+1)+PHEP(2,NHEP+2)
918 PHEP(3,NHEP+3)=PHEP(3,NHEP+1)+PHEP(3,NHEP+2)
920 PHEP(4,NHEP+3)=PHEP(4,NHEP+1)+PHEP(4,NHEP+2)
922 CALL HWUMAS ( PHEP(1,NHEP+3) )
938 PHEP(1,NHEP+4)=PHEP(1,NHEP+1)-PGAMMA(1)
940 PHEP(2,NHEP+4)=PHEP(2,NHEP+1)-PGAMMA(2)
942 PHEP(3,NHEP+4)=PHEP(3,NHEP+1)-PGAMMA(3)
944 PHEP(4,NHEP+4)=PHEP(4,NHEP+1)-PGAMMA(4)
946 PHEP(5,NHEP+4)=RMASS(IDHW(1))
948 IF (CHARGD) PHEP(5,NHEP+4)=ZERO
950 C---First Final parton: quark (or J/psi) in Boson-Gluon Fusion
952 C--- quark or antiquark in QCD Compton
954 PHEP(1,NHEP+5)=Q1*SINTHE*COS(PHI)
956 PHEP(2,NHEP+5)=Q1*SINTHE*SIN(PHI)
958 PHEP(3,NHEP+5)=Q1*COSTHE
964 C---Second Final parton: antiquark in Boson-Gluon Fusion
966 C--- gluon in QCD Compton
968 PHEP(1,NHEP+6)=-PHEP(1,NHEP+5)
970 PHEP(2,NHEP+6)=-PHEP(2,NHEP+5)
972 PHEP(3,NHEP+6)=-PHEP(3,NHEP+5)
978 C---Boost to lepton-hadron CM frame
980 PEP(1) = PHEP(1,NHEP+1)
982 PEP(2) = PHEP(2,NHEP+1)
984 PEP(3) = PHEP(3,NHEP+1) + PPROT
986 PEP(4) = PHEP(4,NHEP+1) + EPROT
992 CALL HWULOF (PEP,PHEP(1,NHEP+I),PHEP(1,NHEP+I))
996 CALL HWULOF (PEP,PROTON,PROTON)
998 CALL HWULOF (PEP,PGAMMA,PGAMMA)
1000 C---Rotation around y-axis to align lepton beam with z-axis
1002 COSPHI = PHEP(3,NHEP+1) /
1004 & SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 )
1006 SINPHI = PHEP(1,NHEP+1) /
1008 & SQRT( PHEP(1,NHEP+1)**2 + PHEP(3,NHEP+1)**2 )
1032 CALL HWUROF (ROT,PHEP(1,NHEP+I),PHEP(1,NHEP+I))
1036 CALL HWUROF (ROT,PROTON,PROTON)
1038 CALL HWUROF (ROT,PGAMMA,PGAMMA)
1040 C---Boost to the LAB frame
1046 CALL HWULOB (PHEP(1,ICMF),PHEP(1,NHEP+I),PHEP(1,NHEP+I))
1050 CALL HWULOB (PHEP(1,ICMF),PROTON,PROTON)
1052 CALL HWULOB (PHEP(1,ICMF),PGAMMA,PGAMMA)
1054 C---Random azimuthal rotation
1056 CALL HWRAZM (ONE,COSAZI,SINAZI)
1068 ROTAZI(1,1) = COSAZI
1070 ROTAZI(1,2) = SINAZI
1072 ROTAZI(2,1) = -SINAZI
1074 ROTAZI(2,2) = COSAZI
1080 CALL HWUROF (ROTAZI,PHEP(1,NHEP+I),PHEP(1,NHEP+I))
1084 CALL HWUROF (ROTAZI,PROTON,PROTON)
1086 CALL HWUROF (ROTAZI,PGAMMA,PGAMMA)
1092 *CMZ :- -03/07/95 19.02.12 by Giovanni Abbiendi
1094 *-- Author : Giovanni Abbiendi & Luca Stanco
1096 C-----------------------------------------------------------------------
1098 SUBROUTINE HWHBRN (*)
1100 C----------------------------------------------------------------------
1102 C Returns a point in the phase space (Y,Q2,SHAT,Z,PHI) and the
1104 C corresponding Jacobian factor AJACOB
1106 C Fill the logical vector INSIDE to tag contributing subprocesses
1108 C to the cross-section
1110 C-----------------------------------------------------------------------
1112 INCLUDE 'HERWIG61.INC'
1114 DOUBLE PRECISION HWRUNI,HWR,HWUPCM,LEP,Y,Q2,SHAT,Z,PHI,AJACOB,
1116 & DSIGMA,ME,MP,ML,MREMIF(18),MFIN1(18),MFIN2(18),RS,SMA,W2,RSHAT,
1118 & MF1,MF2,YMIN,YMAX,YJAC,Q2INF,Q2SUP,Q2JAC,EMW2,ZMIN,ZMAX,ZJAC,
1120 & GAMMA2,LAMBDA,PHIJAC,ZINT,ZLMIN,ZL,EMW,TMIN,TMAX,EMLMIN,EMLMAX,
1122 & SHMIN,EMMIF(18),EMMAF(18),WMIF(18),WMIN,MREMIN,YMIF(18),Q1CM(18),
1124 & Q2MAF(18),EMMAWF(18),ZMIF(18),ZMAF(18),PLMAX,PINC,SHINF,SHSUP,
1126 & SHJAC,CTHLIM,Q1,DETDSH,SRY,SRY0,SRY1
1128 INTEGER IQK,IFLAVU,IFLAVD,I,IMIN,IMAX,IFL,IPROO,IHAD,NTRY,DEBUG
1130 LOGICAL CHARGD,INCLUD(18),INSIDE(18)
1132 EXTERNAL HWRUNI,HWR,HWUPCM
1134 SAVE EMLMIN,EMLMAX,EMMIF,EMMAF,MREMIN,MF1,MF2,YMIF,
1136 & YMIN,YMAX,WMIN,WMIF
1138 COMMON /HWAREA/ LEP,Y,Q2,SHAT,Z,PHI,AJACOB,DSIGMA,ME,MP,ML,MREMIF,
1140 & MFIN1,MFIN2,RS,SMA,W2,RSHAT,IQK,IFLAVU,IFLAVD,IMIN,IMAX,IFL,
1142 & IPROO,CHARGD,INCLUD,INSIDE
1144 EQUIVALENCE (EMW,RMASS(198))
1150 IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
1154 IF (FSTWGT.OR.IHAD.NE.2) THEN
1158 MP = RMASS(IDHW(IHAD))
1162 SMA = RS**2-ME**2-MP**2
1164 PINC = HWUPCM(RS,ME,MP)
1172 YMAX = ONE - TWO*ML*MP / SMA
1174 YMAX = MIN(YMAX,YBMAX)
1178 IF (LEP.EQ.ONE) THEN
1192 SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 +
1194 + TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2)
1196 EMLMIN=MAX(EMMIN,SQRT(SHMIN))
1198 EMLMAX=MIN(EMMAX,RS-ML-MREMIN)
1202 IF (EMLMIN.GT.EMLMAX) GOTO 888
1206 PLMAX=HWUPCM(RS,ML,WMIN)
1208 YMIN = ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+
1212 YMIN = MAX(YMIN,YBMIN)
1216 IF (YMIN.GT.YMAX) GOTO 888
1224 YMAX = ONE - TWO*ML*MP / SMA
1226 YMAX = MIN(YMAX,YBMAX)
1240 C---Boson-Gluon Fusion (also J/Psi) and QCD Compton with struck u or d
1260 C---QCD Compton with struck non-valence parton
1262 MREMIF(I)=MP+RMASS(I-6)
1272 IF (IFL.EQ.164) THEN
1282 C---y boundaries for different flavours and processes
1294 SHMIN = MF1**2+MF2**2 + TWO * PTMIN**2 +
1296 + TWO * SQRT(PTMIN**2+MF1**2) * SQRT(PTMIN**2+MF2**2)
1298 EMMIF(I) = MAX(EMMIN,SQRT(SHMIN))
1300 EMMAF(I) = MIN(EMMAX,RS-ML-MREMIN)
1302 IF (EMMIF(I).GT.EMMAF(I)) THEN
1306 CALL HWWARN('HWHBRN',3,*999)
1312 WMIF(I) = EMMIF(I)+MREMIF(I)
1316 PLMAX = HWUPCM(RS,ML,WMIN)
1318 YMIF(I)=ONE-TWO*(SQRT(PINC**2+MP**2)*SQRT(PLMAX**2+ML**2)+
1322 IF (YMIF(I).GT.YMAX) THEN
1326 CALL HWWARN('HWHBRN',4,*999)
1336 C---considering the largest boundaries
1342 IF (IPROO.EQ.3) THEN
1344 EMLMIN=MIN(EMMIF(IMIN),EMMIF(IMIN+6))
1346 EMLMAX=MAX(EMMAF(IMIN),EMMAF(IMIN+6))
1352 IF (EMLMIN.GT.EMLMAX) GOTO 888
1356 IF (IPROO.EQ.3) YMIN=MIN(YMIF(IMIN),YMIF(IMIN+6))
1358 YMIN = MAX(YMIN,YBMIN)
1362 IF (YMIN.GT.YMAX) GOTO 888
1366 MREMIN = MREMIF(IMIN)
1372 IF (IPROO.EQ.3) THEN
1374 WMIN = MIN(WMIF(IMIN),WMIF(IMIN+6))
1376 MREMIN = MIN(MREMIF(IMIN),MREMIF(IMIN+6))
1384 C---Random generation in largest phase space
1400 IF (.NOT.CHARGD) THEN
1402 IF (IFL.LE.5.OR.(IFL.GE.7.AND.IFL.LE.18)) THEN
1408 SRY = HWRUNI(0,SRY0,SRY1)
1412 YJAC = TWO*SRY*(SRY1-SRY0)
1414 ELSEIF (IFL.EQ.6) THEN
1416 Y = SQRT(HWRUNI(0,YMIN**2,YMAX**2))
1418 YJAC = HALF * (YMAX**2-YMIN**2) / Y
1420 ELSEIF (IFL.EQ.164) THEN
1422 C---in J/psi photoproduction Y and Q2 are given by the Equivalent Photon
1430 IF (NTRY.GT.NETRY) CALL HWWARN('HWHBRN',50,*10)
1432 Y = (YMIN/YMAX)**HWR()*YMAX
1434 IF (ONE+(ONE-Y)**2.LT.TWO*HWR()) GOTO 20
1436 YJAC=(TWO*LOG(YMAX/YMIN)-TWO*(YMAX-YMIN)
1438 & +HALF*(YMAX**2-YMIN**2))
1446 Y = EXP(HWRUNI(0,LOG(YMIN),LOG(YMAX)))
1448 YJAC = Y * LOG(YMAX/YMIN)
1452 Y = HWRUNI(0,YMIN,YMAX)
1462 Q2INF = ME**2*Y**2 / (ONE-Y)
1464 Q2SUP = MP**2 + SMA*Y - WMIN**2
1466 IF (IFL.EQ.164) THEN
1468 Q2INF = MAX(Q2INF,Q2WWMN)
1470 Q2SUP = MIN(Q2SUP,Q2WWMX)
1474 Q2INF = MAX(Q2INF,Q2MIN)
1476 Q2SUP = MIN(Q2SUP,Q2MAX)
1482 IF (Q2INF .GT. Q2SUP) GOTO 888
1486 IF (.NOT.CHARGD) THEN
1488 IF (IFL.EQ.164) THEN
1490 Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP)))
1492 Q2JAC = LOG(Q2SUP/Q2INF)
1494 ELSEIF (Q2INF.LT.RMASS(4)**2) THEN
1496 Q2 = EXP(HWRUNI(0,LOG(Q2INF),LOG(Q2SUP)))
1498 Q2JAC = Q2 * LOG(Q2SUP/Q2INF)
1502 Q2 = Q2INF*Q2SUP/HWRUNI(0,Q2INF,Q2SUP)
1504 Q2JAC = Q2**2 * (Q2SUP-Q2INF)/(Q2SUP*Q2INF)
1512 Q2=(Q2INF+EMW2)*(Q2SUP+EMW2)/(HWRUNI(0,Q2INF,Q2SUP)+EMW2)-EMW2
1514 Q2JAC=(Q2+EMW2)**2*(Q2SUP-Q2INF)/((Q2SUP+EMW2)*(Q2INF+EMW2))
1518 W2 = MP**2 + SMA*Y - Q2
1520 C---s_hat generation
1524 SHSUP = (MIN(SQRT(W2)-MREMIN,EMLMAX))**2
1528 IF (SHINF .GT. SHSUP) GOTO 888
1532 IF (IPRO.EQ.91) THEN
1534 IF (.NOT.CHARGD) THEN
1536 SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP)
1538 SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF)
1542 SHAT = EXP(HWRUNI(0,LOG(SHINF),LOG(SHSUP)))
1544 SHJAC = SHAT*(LOG(SHSUP/SHINF))
1552 IF (SHINF.GT.EMW2+10*GAMW*EMW) THEN
1554 SHAT = SHINF*SHSUP/HWRUNI(0,SHINF,SHSUP)
1556 SHJAC = SHAT**2 * (SHSUP-SHINF)/(SHSUP*SHINF)
1558 ELSEIF (SHSUP.LT.EMW2-10*EMW*GAMW) THEN
1560 SHAT = HWRUNI(0,SHINF,SHSUP)
1566 TMIN=ATAN((SHINF-EMW2)/(GAMW*EMW))
1568 TMAX=ATAN((SHSUP-EMW2)/(GAMW*EMW))
1570 SHAT = GAMW*EMW*TAN(HWRUNI(0,TMIN,TMAX))+EMW2
1572 SHJAC=((SHAT-EMW2)**2+(GAMW*EMW)**2)/(GAMW*EMW)*(TMAX-TMIN)
1590 IF (.NOT.CHARGD) THEN
1606 Q1CM(I) = HWUPCM( RSHAT, MFIN1(I), MFIN2(I) )
1608 IF (Q1CM(I) .LT. PTMIN) THEN
1616 CTHLIM = SQRT(ONE - (PTMIN / Q1CM(I))**2)
1618 GAMMA2 = SHAT + MFIN1(I)**2 - MFIN2(I)**2
1620 LAMBDA = (SHAT-MFIN1(I)**2-MFIN2(I)**2)**2 -
1622 + 4.D0*MFIN1(I)**2*MFIN2(I)**2
1624 ZMIF(I) = (GAMMA2 - SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1626 ZMIF(I) = MAX(ZMIF(I),ZERO)
1628 ZMAF(I) = (GAMMA2 + SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1630 ZMAF(I) = MIN(ZMAF(I),ONE)
1632 ZMIN = MIN( ZMIN, ZMIF(I) )
1634 ZMAX = MAX( ZMAX, ZMAF(I) )
1640 IF (IFL.EQ.164) ZMAX=MIN(ZMAX,ZJMAX)
1644 Q1 = HWUPCM(RSHAT,MF1,MF2)
1648 IF (Q1.LT.PTMIN) GOTO 888
1650 CTHLIM = SQRT(ONE-(PTMIN/Q1)**2)
1652 GAMMA2 = SHAT+MF1**2-MF2**2
1654 LAMBDA = (SHAT-MF1**2-MF2**2)**2-4.D0*MF1**2*MF2**2
1656 ZMIN = (GAMMA2-SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1658 ZMIN = MAX(ZMIN,1D-6)
1660 ZMAX = (GAMMA2+SQRT(LAMBDA)*CTHLIM)/TWO/SHAT
1662 ZMAX = MIN(ZMAX,ONE-1D-6)
1668 IF (ZMIN .GT. ZMAX) GOTO 888
1670 ZLMIN = LOG(ZMIN/(ONE-ZMIN))
1672 ZINT = LOG(ZMAX/(ONE-ZMAX)) - LOG(ZMIN/(ONE-ZMIN))
1674 ZL = ZLMIN+HWR()*ZINT
1676 Z = EXP(ZL)/(ONE+EXP(ZL))
1678 ZJAC = Z*(ONE-Z)*ZINT
1684 IF ((Y.LT.YMIN.OR.Y.GT.YMAX).OR.(Q2.LT.Q2INF.OR.Q2.GT.Q2SUP).OR.
1686 + (SHAT.LT.SHINF.OR.SHAT.GT.SHSUP).OR.(Z.LT.ZMIN.OR.Z.GT.ZMAX))
1692 PHI = HWRUNI(0,ZERO,2*PIFAC)
1696 IF (IFL.EQ.164) PHIJAC=ONE
1700 AJACOB = YJAC * Q2JAC * SHJAC * ZJAC * PHIJAC
1704 IF (IQK.NE.0.OR.IPRO.EQ.5) GOTO 999
1706 C---contributing subprocesses: filling of logical vector INSIDE
1722 IF ( Y.LT.YMIF(I) ) GOTO 200
1726 Q2MAF(I) = MP**2 + SMA*Y - WMIF(I)**2
1728 Q2MAF(I) = MIN( Q2MAF(I), Q2MAX)
1730 IF (Q2INF .GT. Q2MAF(I)) GOTO 200
1732 IF (Q2.LT.Q2INF .OR. Q2.GT.Q2MAF(I)) GOTO 200
1736 EMMAWF(I) = SQRT(W2) - MREMIF(I)
1738 EMMAWF(I) = MIN( EMMAWF(I), EMLMAX )
1742 IF (EMMIF(I) .GT. EMMAWF(I)) GOTO 200
1744 IF (SHAT.LT.EMMIF(I)**2.OR.SHAT.GT.EMMAWF(I)**2) GOTO 200
1748 IF (ZMIF(I) .GT. ZMAF(I)) GOTO 200
1750 IF (Z.LT.ZMIF(I) .OR. Z.GT.ZMAF(I)) GOTO 200
1762 C---UNCOMMENT THIS LINE TO GET A DEBUGGING WARNING FOR NO PHASE-SPACE
1764 C CALL HWWARN('HWHBRN',DEBUG,*777)