PYTHIA/pythia/pyinit.F

   1 C*********************************************************************
   2 C*********************************************************************
   3 C*                                                                  **
   4 C*                                                 December 1993    **
   5 C*                                                                  **
   6 C*           The Lund Monte Carlo for Hadronic Processes            **
   7 C*                                                                  **
   8 C*                        PYTHIA version 5.7                        **
   9 C*                                                                  **
  10 C*                        Torbjorn Sjostrand                        **
  11 C*                Department of theoretical physics 2               **
  12 C*                        University of Lund                        **
  13 C*               Solvegatan 14A, S-223 62 Lund, Sweden              **
  14 C*                    E-mail torbjorn@thep.lu.se                    **
  15 C*                    phone +46 - 46 - 222 48 16                    **
  16 C*                                                                  **
  17 C*         Several parts are written by Hans-Uno Bengtsson          **
  18 C*     CTEQ 2 parton distributions are by the CTEQ collaboration    **
  19 C*   SaS photon parton distributions together with Gerhard Schuler  **
  20 C*    g + g -> Z + b + bbar matrix element code by Ronald Kleiss    **
  21 C*     g + g and q + qbar -> t + tbar + H code by Zoltan Kunszt     **
  22 C*                                                                  **
  23 C*   The latest program version and documentation is found on WWW   **
  24 C*         http://thep.lu.se/tf2/staff/torbjorn/Welcome.html        **
  25 C*                                                                  **
  26 C*        Copyright Torbjorn Sjostrand and CERN, Geneva 1993        **
  27 C*                                                                  **
  28 C*********************************************************************
  29 C*********************************************************************
  30 C                                                                    *
  31 C  List of subprograms in order of appearance, with main purpose     *
  32 C  (S = subroutine, F = function, B = block data)                    *
  33 C                                                                    *
  34 C  S   PYINIT   to administer the initialization procedure           *
  35 C  S   PYEVNT   to administer the generation of an event             *
  36 C  S   PYSTAT   to print cross-section and other information         *
  37 C  S   PYINRE   to initialize treatment of resonances                *
  38 C  S   PYINBM   to read in beam, target and frame choices            *
  39 C  S   PYINKI   to initialize kinematics of incoming particles       *
  40 C  S   PYINPR   to set up the selection of included processes        *
  41 C  S   PYXTOT   to give total, elastic and diffractive cross-sect.   *
  42 C  S   PYMAXI   to find differential cross-section maxima            *
  43 C  S   PYPILE   to select multiplicity of pileup events              *
  44 C  S   PYSAVE   to save alternatives for gamma-p and gamma-gamma     *
  45 C  S   PYRAND   to select subprocess and kinematics for event        *
  46 C  S   PYSCAT   to set up kinematics and colour flow of event        *
  47 C  S   PYSSPA   to simulate initial state spacelike showers          *
  48 C  S   PYRESD   to perform resonance decays                          *
  49 C  S   PYMULT   to generate multiple interactions                    *
  50 C  S   PYREMN   to add on target remnants                            *
  51 C  S   PYDIFF   to set up kinematics for diffractive events          *
  52 C  S   PYDOCU   to compute cross-sections and handle documentation   *
  53 C  S   PYFRAM   to perform boosts between different frames           *
  54 C  S   PYWIDT   to calculate full and partial widths of resonances   *
  55 C  S   PYOFSH   to calculate partial width into off-shell channels   *
  56 C  S   PYKLIM   to calculate borders of allowed kinematical region   *
  57 C  S   PYKMAP   to construct value of kinematical variable           *
  58 C  S   PYSIGH   to calculate differential cross-sections             *
  59 C  S   PYSTFU   to evaluate structure functions                      *
  60 C  S   PYSTFL   to evaluate structure functions at low x and Q^2     *
  61 C  S   PYSTEL   to evaluate electron structure function              *
  62 C  S   PYSTGA   to evaluate photon structure function (generic)      *
  63 C  S   PYGGAM   to evaluate photon structure function (SaS sets)     *
  64 C  S   PYGVMD   to evaluate VMD part of photon structure functions   *
  65 C  S   PYGANO   to evaluate anomalous part of photon str. func.      *
  66 C  S   PYGBEH   to evaluate Bethe-Heitler part of photon str. func.  *
  67 C  S   PYGDIR   to evaluate direct contribution to photon str. func. *
  68 C  S   PYSTPI   to evaluate pion structure function                  *
  69 C  S   PYSTPR   to evaluate proton structure function                *
  70 C  F   PYCTQ2   to evaluate the CTEQ 2 proton structure function     *
  71 C  F   PYHFTH   to evaluate threshold factor for heavy flavour       *
  72 C  S   PYSPLI   to find flavours left in hadron when one removed     *
  73 C  F   PYGAMM   to evaluate ordinary Gamma function Gamma(x)         *
  74 C  S   PYWAUX   to evaluate auxiliary functions W1(s) and W2(s)      *
  75 C  S   PYI3AU   to evaluate auxiliary function I3(s,t,u,v)           *
  76 C  F   PYSPEN   to evaluate Spence (dilogarithm) function Sp(x)      *
  77 C  S   PYQQBH   to evaluate matrix element for g + g -> Q + Q~ + H   *
  78 C  S   PYTEST   to test the proper functioning of the package        *
  79 C  B   PYDATA   to contain all default values                        *
  80 C  S   PYKCUT   to provide dummy routine for user kinematical cuts   *
  81 C  S   PYEVWT   to provide dummy routine for weighting events        *
  82 C  S   PYUPIN   to initialize a user process                         *
  83 C  S   PYUPEV   to generate a user process event (dummy routine)     *
  84 C  S   PDFSET   dummy routine to be removed when using PDFLIB        *
  85 C  S   STRUCTM  dummy routine to be removed when using PDFLIB        *
  86 C                                                                    *
  87 C*********************************************************************
  88
  89       SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
  90
  91 C...Initializes the generation procedure; finds maxima of the
  92 C...differential cross-sections to be used for weighting.
  93       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  94       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  95       COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
  96       COMMON/LUDAT4/CHAF(500)
  97       CHARACTER CHAF*8
  98       COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
  99       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
 100       COMMON/PYINT1/MINT(400),VINT(400)
 101       COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
 102       COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
 103       COMMON/PYINT9/DXSEC(0:200)
 104       DOUBLE PRECISION DXSEC
 105       SAVE /LUDAT1/,/LUDAT2/,/LUDAT3/,/LUDAT4/
 106       SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT5/,/PYINT9/
 107       DIMENSION ALAMIN(20),NFIN(20)
 108       CHARACTER*(*) FRAME,BEAM,TARGET
 109       CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHLH(2)*6
 110
 111 C...Interface to PDFLIB.
 112       COMMON/W50512/QCDL4,QCDL5
 113       SAVE /W50512/
 114       DOUBLE PRECISION VALUE(20),QCDL4,QCDL5
 115       CHARACTER*20 PARM(20)
 116       DATA VALUE/20*0D0/,PARM/20*' '/
 117
 118 C...Data:Lambda and n_f values for structure functions; months.
 119       DATA ALAMIN/0.20,0.29,0.20,0.40,0.213,0.208,0.208,0.322,
 120      &0.190,0.235,10*0.2/,NFIN/20*4/
 121       DATA CHLH/'lepton','hadron'/
 122
 123 C...Reset MINT and VINT arrays. Write headers.
 124       DO 100 J=1,400
 125       MINT(J)=0
 126       VINT(J)=0.
 127   100 CONTINUE
 128       IF(MSTU(12).GE.1) CALL LULIST(0)
 129       IF(MSTP(122).GE.1) WRITE(MSTU(11),5100)
 130
 131 C...Maximum 4 generations; set maximum number of allowed flavours.
 132       MSTP(1)=MIN(4,MSTP(1))
 133       MSTU(114)=MIN(MSTU(114),2*MSTP(1))
 134       MSTP(58)=MIN(MSTP(58),2*MSTP(1))
 135
 136 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton.
 137       DO 120 I=-20,20
 138       VINT(180+I)=0.
 139       IA=IABS(I)
 140       IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN
 141         DO 110 J=1,MSTP(1)
 142         IB=2*J-1+MOD(IA,2)
 143         IPM=(5-ISIGN(1,I))/2
 144         IDC=J+MDCY(IA,2)+2
 145         IF(MDME(IDC,1).EQ.1.OR.MDME(IDC,1).EQ.IPM) VINT(180+I)=
 146      &  VINT(180+I)+VCKM((IA+1)/2,(IB+1)/2)
 147   110   CONTINUE
 148       ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN
 149         VINT(180+I)=1.
 150       ENDIF
 151   120 CONTINUE
 152
 153 C...Initialize structure functions: PDFLIB.
 154       IF(MSTP(52).EQ.2) THEN
 155         PARM(1)='NPTYPE'
 156         VALUE(1)=1
 157         PARM(2)='NGROUP'
 158         VALUE(2)=MSTP(51)/1000
 159         PARM(3)='NSET'
 160         VALUE(3)=MOD(MSTP(51),1000)
 161         PARM(4)='TMAS'
 162         VALUE(4)=PMAS(6,1)
 163         CALL PDFSET(PARM,VALUE)
 164         MINT(93)=1000000+MSTP(51)
 165       ENDIF
 166
 167 C...Choose Lambda value to use in alpha-strong.
 168       MSTU(111)=MSTP(2)
 169       IF(MSTP(3).GE.2) THEN
 170         ALAM=0.2
 171         NF=4
 172         IF(MSTP(52).EQ.1.AND.MSTP(51).GE.1.AND.MSTP(51).LE.10) THEN
 173           ALAM=ALAMIN(MSTP(51))
 174           NF=NFIN(MSTP(51))
 175         ELSEIF(MSTP(52).EQ.2) THEN
 176           ALAM=QCDL4
 177           NF=4
 178         ENDIF
 179         PARP(1)=ALAM
 180         PARP(61)=ALAM
 181         PARP(72)=ALAM
 182         PARU(112)=ALAM
 183         MSTU(112)=NF
 184         IF(MSTP(3).EQ.3) PARJ(81)=ALAM
 185       ENDIF
 186
 187 C...Initialize widths and partial widths for resonances.
 188       CALL PYINRE
 189
 190 C...Identify beam and target particles and frame of process.
 191       CHFRAM=FRAME//' '
 192       CHBEAM=BEAM//' '
 193       CHTARG=TARGET//' '
 194       CALL PYINBM(CHFRAM,CHBEAM,CHTARG,WIN)
 195       IF(MINT(65).EQ.1) GOTO 170
 196
 197 C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives.
 198 C...For e-gamma allow 2 alternatives.
 199       MINT(121)=1
 200       MINT(123)=MSTP(14)
 201       IF(MSTP(14).EQ.10.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
 202         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
 203      &  (IABS(MINT(11)).GE.28.OR.IABS(MINT(12)).GE.28)) MINT(121)=3
 204         IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=6
 205         IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
 206      &  (IABS(MINT(11)).EQ.11.OR.IABS(MINT(12)).EQ.11)) MINT(121)=2
 207       ENDIF
 208
 209 C...Set up kinematics of process.
 210       CALL PYINKI(0)
 211
 212 C...Loop over gamma-p or gamma-gamma alternatives.
 213       DO 160 IGA=1,MINT(121)
 214       MINT(122)=IGA
 215
 216 C...Select partonic subprocesses to be included in the simulation.
 217       CALL PYINPR
 218
 219 C...Count number of subprocesses on.
 220       MINT(48)=0
 221       DO 130 ISUB=1,200
 222       IF(MINT(50).EQ.0.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.
 223      &MSUB(ISUB).EQ.1) THEN
 224         WRITE(MSTU(11),5200) ISUB,CHLH(MINT(41)),CHLH(MINT(42))
 225         STOP
 226       ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN
 227         WRITE(MSTU(11),5300) ISUB
 228         STOP
 229       ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN
 230         WRITE(MSTU(11),5400) ISUB
 231         STOP
 232       ELSEIF(MSUB(ISUB).EQ.1) THEN
 233         MINT(48)=MINT(48)+1
 234       ENDIF
 235   130 CONTINUE
 236       IF(MINT(48).EQ.0) THEN
 237         WRITE(MSTU(11),5500)
 238         STOP
 239       ENDIF
 240       MINT(49)=MINT(48)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94)
 241
 242 C...Reset variables for cross-section calculation.
 243       DO 150 I=0,200
 244       DO 140 J=1,3
 245       NGEN(I,J)=0
 246       XSEC(I,J)=0.
 247   140 CONTINUE
 248       DXSEC(I)=0D0
 249   150 CONTINUE
 250
 251 C...Find parametrized total cross-sections.
 252       CALL PYXTOT
 253
 254 C...Maxima of differential cross-sections.
 255       IF(MSTP(121).LE.1) CALL PYMAXI
 256
 257 C...Initialize possibility of pileup events.
 258       IF(MINT(121).GT.1) MSTP(131)=0
 259       IF(MSTP(131).NE.0) CALL PYPILE(1)
 260
 261 C...Initialize multiple interactions with variable impact parameter.
 262       IF(MINT(50).EQ.1.AND.(MINT(49).NE.0.OR.MSTP(131).NE.0).AND.
 263      &MSTP(82).GE.2) CALL PYMULT(1)
 264
 265 C...Save results for gamma-p and gamma-gamma alternatives.
 266       IF(MINT(121).GT.1) CALL PYSAVE(1,IGA)
 267   160 CONTINUE
 268
 269 C...Initialization finished.
 270   170 IF(MSTP(122).GE.1) WRITE(MSTU(11),5600)
 271
 272 C...Formats for initialization information.
 273  5100 FORMAT('1',18('*'),1X,'PYINIT: initialization of PYTHIA ',
 274      &'routines',1X,17('*'))
 275  5200 FORMAT(1X,'Error: process number ',I3,' not meaningful for ',A6,
 276      &'-',A6,' interactions.'/1X,'Execution stopped!')
 277  5300 FORMAT(1X,'Error: requested subprocess',I4,' not implemented.'/
 278      &1X,'Execution stopped!')
 279  5400 FORMAT(1X,'Error: requested subprocess',I4,' not existing.'/
 280      &1X,'Execution stopped!')
 281  5500 FORMAT(1X,'Error: no subprocess switched on.'/
 282      &1X,'Execution stopped.')
 283  5600 FORMAT(/1X,22('*'),1X,'PYINIT: initialization completed',1X,
 284      &22('*'))
 285
 286       RETURN
 287       END