hijing/src/pyhiinit.f

0001       SUBROUTINE PYHIINIT(FRAME,BEAM,TARGET,WIN)
0002
0003 C...Initializes the generation procedure; finds maxima of the
0004 C...differential cross-sections to be used for weighting.
0005       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0006       SAVE /LUDAT1/
0007       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0008       SAVE /LUDAT2/
0009       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
0010       SAVE /LUDAT3/
0011       COMMON/LUDAT4/CHAF(500)
0012       CHARACTER CHAF*8
0013       SAVE /LUDAT4/
0014       COMMON/PYHISUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
0015       SAVE /PYHISUBS/
0016       COMMON/PYHIPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0017       SAVE /PYHIPARS/
0018       COMMON/PYHIINT1/MINT(400),VINT(400)
0019       SAVE /PYHIINT1/
0020       COMMON/PYHIINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
0021       SAVE /PYHIINT2/
0022       COMMON/PYHIINT5/NGEN(0:200,3),XSEC(0:200,3)
0023       SAVE /PYHIINT5/
0024       CHARACTER*(*) FRAME,BEAM,TARGET
0025       CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHMO(12)*3,CHLH(2)*6
0026       DATA CHMO/'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep',
0027      &'Oct','Nov','Dec'/, CHLH/'lepton','hadron'/
0028
0029 C...Write headers.
0030 C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1000) MSTP(181),MSTP(182),
0031 C     &MSTP(185),CHMO(MSTP(184)),MSTP(183)
0032       CALL LULIST(0)
0033 C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1100)
0034
0035 C...Identify beam and target particles and initialize kinematics.
0036       CHFRAM=FRAME//' '
0037       CHBEAM=BEAM//' '
0038       CHTARG=TARGET//' '
0039       CALL PYHIINKI(CHFRAM,CHBEAM,CHTARG,WIN)
0040
0041 C...Select partonic subprocesses to be included in the simulation.
0042       IF(MSEL.NE.0) THEN
0043         DO 100 I=1,200
0044   100   MSUB(I)=0
0045       ENDIF
0046       IF(MINT(43).EQ.1.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
0047 C...Lepton+lepton -> gamma/Z0 or W.
0048         IF(MINT(11)+MINT(12).EQ.0) MSUB(1)=1
0049         IF(MINT(11)+MINT(12).NE.0) MSUB(2)=1
0050       ELSEIF(MSEL.EQ.1) THEN
0051 C...High-pT QCD processes:
0052         MSUB(11)=1
0053         MSUB(12)=1
0054         MSUB(13)=1
0055         MSUB(28)=1
0056         MSUB(53)=1
0057         MSUB(68)=1
0058         IF(MSTP(82).LE.1.AND.CKIN(3).LT.PARP(81)) MSUB(95)=1
0059         IF(MSTP(82).GE.2.AND.CKIN(3).LT.PARP(82)) MSUB(95)=1
0060       ELSEIF(MSEL.EQ.2) THEN
0061 C...All QCD processes:
0062         MSUB(11)=1
0063         MSUB(12)=1
0064         MSUB(13)=1
0065         MSUB(28)=1
0066         MSUB(53)=1
0067         MSUB(68)=1
0068         MSUB(91)=1
0069         MSUB(92)=1
0070         MSUB(93)=1
0071         MSUB(95)=1
0072       ELSEIF(MSEL.GE.4.AND.MSEL.LE.8) THEN
0073 C...Heavy quark production.
0074         MSUB(81)=1
0075         MSUB(82)=1
0076         DO 110 J=1,MIN(8,MDCY(21,3))
0077   110   MDME(MDCY(21,2)+J-1,1)=0
0078         MDME(MDCY(21,2)+MSEL-1,1)=1
0079       ELSEIF(MSEL.EQ.10) THEN
0080 C...Prompt photon production:
0081         MSUB(14)=1
0082         MSUB(18)=1
0083         MSUB(29)=1
0084       ELSEIF(MSEL.EQ.11) THEN
0085 C...Z0/gamma* production:
0086         MSUB(1)=1
0087       ELSEIF(MSEL.EQ.12) THEN
0088 C...W+/- production:
0089         MSUB(2)=1
0090       ELSEIF(MSEL.EQ.13) THEN
0091 C...Z0 + jet:
0092         MSUB(15)=1
0093         MSUB(30)=1
0094       ELSEIF(MSEL.EQ.14) THEN
0095 C...W+/- + jet:
0096         MSUB(16)=1
0097         MSUB(31)=1
0098       ELSEIF(MSEL.EQ.15) THEN
0099 C...Z0 & W+/- pair production:
0100         MSUB(19)=1
0101         MSUB(20)=1
0102         MSUB(22)=1
0103         MSUB(23)=1
0104         MSUB(25)=1
0105       ELSEIF(MSEL.EQ.16) THEN
0106 C...H0 production:
0107         MSUB(3)=1
0108         MSUB(5)=1
0109         MSUB(8)=1
0110         MSUB(102)=1
0111       ELSEIF(MSEL.EQ.17) THEN
0112 C...H0 & Z0 or W+/- pair production:
0113         MSUB(24)=1
0114         MSUB(26)=1
0115       ELSEIF(MSEL.EQ.21) THEN
0116 C...Z'0 production:
0117         MSUB(141)=1
0118       ELSEIF(MSEL.EQ.22) THEN
0119 C...H+/- production:
0120         MSUB(142)=1
0121       ELSEIF(MSEL.EQ.23) THEN
0122 C...R production:
0123         MSUB(143)=1
0124       ENDIF
0125
0126 C...Count number of subprocesses on.
0127       MINT(44)=0
0128       DO 120 ISUB=1,200
0129       IF(MINT(43).LT.4.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.
0130      &MSUB(ISUB).EQ.1) THEN
0131         WRITE(MSTU(11),1200) ISUB,CHLH(MINT(41)),CHLH(MINT(42))
0132         STOP
0133       ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN
0134         WRITE(MSTU(11),1300) ISUB
0135         STOP
0136       ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN
0137         WRITE(MSTU(11),1400) ISUB
0138         STOP
0139       ELSEIF(MSUB(ISUB).EQ.1) THEN
0140         MINT(44)=MINT(44)+1
0141       ENDIF
0142   120 CONTINUE
0143       IF(MINT(44).EQ.0) THEN
0144         WRITE(MSTU(11),1500)
0145         STOP
0146       ENDIF
0147       MINT(45)=MINT(44)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94)
0148
0149 C...Maximum 4 generations; set maximum number of allowed flavours.
0150       MSTP(1)=MIN(4,MSTP(1))
0151       MSTU(114)=MIN(MSTU(114),2*MSTP(1))
0152       MSTP(54)=MIN(MSTP(54),2*MSTP(1))
0153
0154 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton.
0155       DO 140 I=-20,20
0156       VINT(180+I)=0.
0157       IA=IABS(I)
0158       IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN
0159         DO 130 J=1,MSTP(1)
0160         IB=2*J-1+MOD(IA,2)
0161         IPM=(5-ISIGN(1,I))/2
0162         IDC=J+MDCY(IA,2)+2
0163   130   IF(MDME(IDC,1).EQ.1.OR.MDME(IDC,1).EQ.IPM) VINT(180+I)=
0164      &  VINT(180+I)+VCKM((IA+1)/2,(IB+1)/2)
0165       ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN
0166         VINT(180+I)=1.
0167       ENDIF
0168   140 CONTINUE
0169
0170 C...Choose Lambda value to use in alpha-strong.
0171       MSTU(111)=MSTP(2)
0172       IF(MSTP(3).GE.1) THEN
0173         ALAM=PARP(1)
0174         IF(MSTP(51).EQ.1) ALAM=0.2
0175         IF(MSTP(51).EQ.2) ALAM=0.29
0176         IF(MSTP(51).EQ.3) ALAM=0.2
0177         IF(MSTP(51).EQ.4) ALAM=0.4
0178         IF(MSTP(51).EQ.11) ALAM=0.16
0179         IF(MSTP(51).EQ.12) ALAM=0.26
0180         IF(MSTP(51).EQ.13) ALAM=0.36
0181         PARP(1)=ALAM
0182         PARP(61)=ALAM
0183         PARU(112)=ALAM
0184         PARJ(81)=ALAM
0185       ENDIF
0186
0187 C...Initialize widths and partial widths for resonances.
0188       CALL PYHIINRE
0189
0190 C...Reset variables for cross-section calculation.
0191       DO 150 I=0,200
0192       DO 150 J=1,3
0193       NGEN(I,J)=0
0194   150 XSEC(I,J)=0.
0195       VINT(108)=0.
0196
0197 C...Find parametrized total cross-sections.
0198       IF(MINT(43).EQ.4) CALL PYHIXTOT
0199
0200 C...Maxima of differential cross-sections.
0201       IF(MSTP(121).LE.0) CALL PYHIMAXI
0202
0203 C...Initialize possibility of overlayed events.
0204       IF(MSTP(131).NE.0) CALL PYHIOVLY(1)
0205
0206 C...Initialize multiple interactions with variable impact parameter.
0207       IF(MINT(43).EQ.4.AND.(MINT(45).NE.0.OR.MSTP(131).NE.0).AND.
0208      &MSTP(82).GE.2) CALL PYHIMULT(1)
0209 C      IF(MSTP(122).GE.1) WRITE(MSTU(11),1600)
0210
0211 C...Formats for initialization information.
0212  1000 FORMAT(///20X,'The Lund Monte Carlo - PYHITHIA version ',I1,
0213      &'.',I1/
0214      &20X,'**  Last date of change:  ',I2,1X,A3,1X,I4,'  **'/)
0215  1100 FORMAT('1',18('*'),1X,'PYHIINIT: initialization of PYHITHIA ',
0216      &'(hijing pythia) routines',1X,17('*'))
0217  1200 FORMAT(1X,'Error: process number ',I3,' not meaningful for ',A6,
0218      &'-',A6,' interactions.'/1X,'Execution stopped!')
0219  1300 FORMAT(1X,'Error: requested subprocess',I4,' not implemented.'/
0220      &1X,'Execution stopped!')
0221  1400 FORMAT(1X,'Error: requested subprocess',I4,' not existing.'/
0222      &1X,'Execution stopped!')
0223  1500 FORMAT(1X,'Error: no subprocess switched on.'/
0224      &1X,'Execution stopped.')
0225  1600 FORMAT(/1X,22('*'),1X,'PYHIINIT: initialization completed',1X,
0226      &22('*'))
0227
0228       RETURN
0229       END