pythiaeRHIC/pythia/pyinit.f

0001
0002 C*********************************************************************
0003
0004 C...PYINIT
0005 C...Initializes the generation procedure; finds maxima of the
0006 C...differential cross-sections to be used for weighting.
0007
0008       SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0016       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0017       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0018       COMMON/PYDAT4/CHAF(500,2)
0019       CHARACTER CHAF*16
0020       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0021       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0022       COMMON/PYINT1/MINT(400),VINT(400)
0023       COMMON/PYINT2/ISET(500),KFPR(500,2),COEF(500,20),ICOL(40,4,2)
0024       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0025       SAVE /PYDAT1/,/PYDAT2/,/PYDAT3/,/PYDAT4/,/PYSUBS/,/PYPARS/,
0026      &/PYINT1/,/PYINT2/,/PYINT5/
0027 C...Local arrays and character variables.
0028       DIMENSION ALAMIN(20),NFIN(20)
0029       CHARACTER*(*) FRAME,BEAM,TARGET
0030       CHARACTER CHFRAM*12,CHBEAM*12,CHTARG*12,CHLH(2)*6
0031
0032 C...Interface to PDFLIB.
0033       COMMON/W50511/NPTYPE,NGROUP,NSET,MODE,NFL,LO,TMAS
0034       COMMON/W50512/QCDL4,QCDL5
0035       SAVE /W50511/,/W50512/
0036       DOUBLE PRECISION VALUE(20),TMAS,QCDL4,QCDL5
0037       CHARACTER*20 PARM(20)
0038       DATA VALUE/20*0D0/,PARM/20*' '/
0039
0040 C...Data:Lambda and n_f values for parton distributions..
0041       DATA ALAMIN/0.177D0,0.239D0,0.247D0,0.2322D0,0.248D0,0.248D0,
0042      &0.192D0,0.326D0,2*0.2D0,0.2D0,0.2D0,0.29D0,0.2D0,0.4D0,5*0.2D0/,
0043      &NFIN/20*4/
0044       DATA CHLH/'lepton','hadron'/
0045
0046 C...Check that BLOCK DATA PYDATA has been loaded.
0047       CALL PYCKBD
0048
0049 C...Reset MINT and VINT arrays. Write headers.
0050       MSTI(53)=0
0051       DO 100 J=1,400
0052         MINT(J)=0
0053         VINT(J)=0D0
0054   100 CONTINUE
0055       IF(MSTU(12).NE.12345) CALL PYLIST(0)
0056       IF(MSTP(122).GE.1) WRITE(MSTU(11),5100)
0057
0058 C...Reset error counters.
0059       MSTU(23)=0
0060       MSTU(27)=0
0061       MSTU(30)=0
0062
0063 C...Reset processes that should not be on.
0064       MSUB(96)=0
0065       MSUB(97)=0
0066
0067 C...Select global FSR/ISR/UE parameter set = 'tune'
0068 C...See routine PYTUNE for details
0069       IF (MSTP(5).NE.0) THEN
0070         MSTP5=MSTP(5)
0071         CALL PYTUNE(MSTP5)
0072       ENDIF
0073
0074 C...Call user process initialization routine.
0075       IF(FRAME(1:1).EQ.'u'.OR.FRAME(1:1).EQ.'U') THEN
0076         MSEL=0
0077         CALL UPINIT
0078         MSEL=0
0079       ENDIF
0080
0081 C...Maximum 4 generations; set maximum number of allowed flavours.
0082       MSTP(1)=MIN(4,MSTP(1))
0083       MSTU(114)=MIN(MSTU(114),2*MSTP(1))
0084       MSTP(58)=MIN(MSTP(58),2*MSTP(1))
0085
0086 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton.
0087       DO 120 I=-20,20
0088         VINT(180+I)=0D0
0089         IA=IABS(I)
0090         IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN
0091           DO 110 J=1,MSTP(1)
0092             IB=2*J-1+MOD(IA,2)
0093             IF(IB.GE.6.AND.MSTP(9).EQ.0) GOTO 110
0094             IPM=(5-ISIGN(1,I))/2
0095             IDC=J+MDCY(IA,2)+2
0096             IF(MDME(IDC,1).EQ.1.OR.MDME(IDC,1).EQ.IPM) VINT(180+I)=
0097      &      VINT(180+I)+VCKM((IA+1)/2,(IB+1)/2)
0098   110     CONTINUE
0099         ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN
0100           VINT(180+I)=1D0
0101         ENDIF
0102   120 CONTINUE
0103
0104 C...Initialize parton distributions: PDFLIB.
0105       IF(MSTP(52).EQ.2) THEN
0106         PARM(1)='NPTYPE'
0107         VALUE(1)=1
0108         PARM(2)='NGROUP'
0109         VALUE(2)=MSTP(51)/1000
0110         PARM(3)='NSET'
0111         VALUE(3)=MOD(MSTP(51),1000)
0112         PARM(4)='TMAS'
0113         VALUE(4)=PMAS(6,1)
0114         CALL PDFSET(PARM,VALUE)
0115         MINT(93)=1000000+MSTP(51)
0116       ENDIF
0117
0118 C...Choose Lambda value to use in alpha-strong.
0119       MSTU(111)=MSTP(2)
0120       IF(MSTP(3).GE.2) THEN
0121         ALAM=0.2D0
0122         NF=4
0123         IF(MSTP(52).EQ.1.AND.MSTP(51).GE.1.AND.MSTP(51).LE.20) THEN
0124           ALAM=ALAMIN(MSTP(51))
0125           NF=NFIN(MSTP(51))
0126         ELSEIF(MSTP(52).EQ.2.AND.NFL.EQ.5) THEN
0127           ALAM=QCDL5
0128           NF=5
0129         ELSEIF(MSTP(52).EQ.2) THEN
0130           ALAM=QCDL4
0131           NF=4
0132         ENDIF
0133         PARP(1)=ALAM
0134         PARP(61)=ALAM
0135         PARP(72)=ALAM
0136         PARU(112)=ALAM
0137         MSTU(112)=NF
0138         IF(MSTP(3).EQ.3) PARJ(81)=ALAM
0139       ENDIF
0140
0141 C...Initialize the SUSY generation: couplings, masses,
0142 C...decay modes, branching ratios, and so on.
0143       CALL PYMSIN
0144 C...Initialize widths and partial widths for resonances.
0145       CALL PYINRE
0146 C...Set Z0 mass and width for e+e- routines.
0147       PARJ(123)=PMAS(23,1)
0148       PARJ(124)=PMAS(23,2)
0149
0150 C...Identify beam and target particles and frame of process.
0151       CHFRAM=FRAME//' '
0152       CHBEAM=BEAM//' '
0153       CHTARG=TARGET//' '
0154       CALL PYINBM(CHFRAM,CHBEAM,CHTARG,WIN)
0155       IF(MINT(65).EQ.1) GOTO 170
0156
0157 C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives.
0158 C...For e-gamma allow 2 alternatives.
0159       MINT(121)=1
0160       IF(MSTP(14).EQ.10.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
0161         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
0162      &  (IABS(MINT(11)).GT.100.OR.IABS(MINT(12)).GT.100)) MINT(121)=3
0163         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=6
0164         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
0165      &  (IABS(MINT(11)).EQ.11.OR.IABS(MINT(12)).EQ.11)) MINT(121)=2
0166       ELSEIF(MSTP(14).EQ.20.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
0167         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
0168      &  (IABS(MINT(11)).GT.100.OR.IABS(MINT(12)).GT.100)) MINT(121)=3
0169         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=9
0170       ELSEIF(MSTP(14).EQ.25.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
0171         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
0172      &  (IABS(MINT(11)).GT.100.OR.IABS(MINT(12)).GT.100)) MINT(121)=2
0173         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=4
0174       ELSEIF(MSTP(14).EQ.30.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
0175         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
0176      &  (IABS(MINT(11)).GT.100.OR.IABS(MINT(12)).GT.100)) MINT(121)=4
0177         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=13
0178       ENDIF
0179       MINT(123)=MSTP(14)
0180       IF((MSTP(14).EQ.10.OR.MSTP(14).EQ.20.OR.MSTP(14).EQ.25.OR.
0181      &MSTP(14).EQ.30).AND.MSEL.NE.1.AND.MSEL.NE.2) MINT(123)=0
0182       IF(MSTP(14).GE.11.AND.MSTP(14).LE.19) THEN
0183         IF(MSTP(14).EQ.11) MINT(123)=0
0184         IF(MSTP(14).EQ.12.OR.MSTP(14).EQ.14) MINT(123)=5
0185         IF(MSTP(14).EQ.13.OR.MSTP(14).EQ.17) MINT(123)=6
0186         IF(MSTP(14).EQ.15) MINT(123)=2
0187         IF(MSTP(14).EQ.16.OR.MSTP(14).EQ.18) MINT(123)=7
0188         IF(MSTP(14).EQ.19) MINT(123)=3
0189       ELSEIF(MSTP(14).GE.21.AND.MSTP(14).LE.24) THEN
0190         IF(MSTP(14).EQ.21) MINT(123)=0
0191         IF(MSTP(14).EQ.22.OR.MSTP(14).EQ.23) MINT(123)=4
0192         IF(MSTP(14).EQ.24) MINT(123)=1
0193       ELSEIF(MSTP(14).GE.26.AND.MSTP(14).LE.29) THEN
0194         IF(MSTP(14).EQ.26.OR.MSTP(14).EQ.28) MINT(123)=8
0195         IF(MSTP(14).EQ.27.OR.MSTP(14).EQ.29) MINT(123)=9
0196       ENDIF
0197
0198 C...Set up kinematics of process.
0199       CALL PYINKI(0)
0200
0201 C...Set up kinematics for photons inside leptons.
0202       IF(MINT(141).NE.0.OR.MINT(142).NE.0) CALL PYGAGA(1,WTGAGA)
0203
0204 C...Precalculate flavour selection weights.
0205       CALL PYKFIN
0206
0207 C...Loop over gamma-p or gamma-gamma alternatives.
0208       CKIN3=CKIN(3)
0209       MSAV48=0
0210       DO 160 IGA=1,MINT(121)
0211         CKIN(3)=CKIN3
0212         MINT(122)=IGA
0213
0214 C...Select partonic subprocesses to be included in the simulation.
0215         CALL PYINPR
0216         MINT(101)=1
0217         MINT(102)=1
0218         MINT(103)=MINT(11)
0219         MINT(104)=MINT(12)
0220
0221 C...Count number of subprocesses on.
0222         MINT(48)=0
0223         DO 130 ISUB=1,500
0224           IF(MINT(50).EQ.0.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.
0225      &    MSUB(ISUB).EQ.1.AND.MINT(121).GT.1) THEN
0226             MSUB(ISUB)=0
0227           ELSEIF(MINT(50).EQ.0.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.
0228      &    MSUB(ISUB).EQ.1) THEN
0229             WRITE(MSTU(11),5200) ISUB,CHLH(MINT(41)),CHLH(MINT(42))
0230             CALL PYSTOP(1)
0231           ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN
0232             WRITE(MSTU(11),5300) ISUB
0233             CALL PYSTOP(1)
0234           ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN
0235             WRITE(MSTU(11),5400) ISUB
0236             CALL PYSTOP(1)
0237           ELSEIF(MSUB(ISUB).EQ.1) THEN
0238             MINT(48)=MINT(48)+1
0239           ENDIF
0240   130   CONTINUE
0241
0242 C...Stop or raise warning flag if no subprocesses on.
0243         IF(MINT(121).EQ.1.AND.MINT(48).EQ.0) THEN
0244           IF(MSTP(127).NE.1) THEN
0245             WRITE(MSTU(11),5500)
0246             CALL PYSTOP(1)
0247           ELSE
0248             WRITE(MSTU(11),5700)
0249             MSTI(53)=1
0250           ENDIF
0251         ENDIF
0252         MINT(49)=MINT(48)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94)
0253         MSAV48=MSAV48+MINT(48)
0254
0255 C...Reset variables for cross-section calculation.
0256         DO 150 I=0,500
0257           DO 140 J=1,3
0258             NGEN(I,J)=0
0259             XSEC(I,J)=0D0
0260   140     CONTINUE
0261   150   CONTINUE
0262
0263 C...Find parametrized total cross-sections.
0264         CALL PYXTOT
0265         VINT(318)=VINT(317)
0266
0267 C...Maxima of differential cross-sections.
0268         IF(MSTP(121).LE.1) CALL PYMAXI
0269
0270 C...Initialize possibility of pileup events.
0271         IF(MINT(121).GT.1) MSTP(131)=0
0272         IF(MSTP(131).NE.0) CALL PYPILE(1)
0273
0274 C...Initialize multiple interactions with variable impact parameter.
0275         IF(MINT(50).EQ.1) THEN
0276           PTMN=PARP(82)*(VINT(1)/PARP(89))**PARP(90)
0277           IF(MOD(MSTP(81),10).EQ.0.AND.(CKIN(3).GT.PTMN.OR.
0278      &    ((MSEL.NE.1.AND.MSEL.NE.2)))) MSTP(82)=MIN(1,MSTP(82))
0279           IF((MINT(49).NE.0.OR.MSTP(131).NE.0).AND.MSTP(82).GE.2) THEN
0280             MINT(35)=1
0281             CALL PYMULT(1)
0282             MINT(35)=3
0283             CALL PYMIGN(1)
0284           ENDIF
0285         ENDIF
0286
0287 C...Save results for gamma-p and gamma-gamma alternatives.
0288         IF(MINT(121).GT.1) CALL PYSAVE(1,IGA)
0289   160 CONTINUE
0290
0291 C...Initialization finished.
0292       IF(MSAV48.EQ.0) THEN
0293         IF(MSTP(127).NE.1) THEN
0294           WRITE(MSTU(11),5500)
0295           CALL PYSTOP(1)
0296         ELSE
0297           WRITE(MSTU(11),5700)
0298           MSTI(53)=1
0299         ENDIF
0300       ENDIF
0301   170 IF(MSTP(122).GE.1) WRITE(MSTU(11),5600)
0302
0303 C...Formats for initialization information.
0304  5100 FORMAT('1',18('*'),1X,'PYINIT: initialization of PYTHIA ',
0305      &'routines',1X,17('*'))
0306  5200 FORMAT(1X,'Error: process number ',I3,' not meaningful for ',A6,
0307      &'-',A6,' interactions.'/1X,'Execution stopped!')
0308  5300 FORMAT(1X,'Error: requested subprocess',I4,' not implemented.'/
0309      &1X,'Execution stopped!')
0310  5400 FORMAT(1X,'Error: requested subprocess',I4,' not existing.'/
0311      &1X,'Execution stopped!')
0312  5500 FORMAT(1X,'Error: no subprocess switched on.'/
0313      &1X,'Execution stopped.')
0314  5600 FORMAT(/1X,22('*'),1X,'PYINIT: initialization completed',1X,
0315      &22('*'))
0316  5700 FORMAT(1X,'Error: no subprocess switched on.'/
0317      &1X,'Execution will stop if you try to generate events.')
0318
0319       RETURN
0320       END