TEvtGen/EvtGenModels/Pythia.F

   1 C--------------------------------------------------------------------------
   2 C
   3 C Environment:
   4 C      This software is part of the EvtGen package developed jointly
   5 C      for the BaBar and CLEO collaborations.  If you use all or part
   6 C      of it, please give an appropriate acknowledgement.
   7 C
   8 C Copyright Information: See BelEvtGen/COPYRIGHT
   9 C      Copyright (C) 1998      Caltech, UCSB
  10 C
  11 C Module: EvtJetSetInit.F
  12 C
  13 C Description:
  14 C
  15 C Modification history:
  16 C
  17 C    DJL/RYD     August 11, 1998         Module created
  18 C    RS          October 28, 2002        copied and modified from Jetset
  19 C
  20 C------------------------------------------------------------------------
  21       SUBROUTINE EVTPYTHIAINIT(FNAME)
  22       IMPLICIT NONE
  23       EXTERNAL PYDATA
  24
  25       CHARACTER*(*) FNAME
  26
  27       PRINT *,'PYUPDA : : ',FNAME
  28       OPEN(54,STATUS='OLD',FILE=FNAME)
  29       CALL PYUPDA(3,54)
  30       CLOSE(54)
  31       WRITE (*,*) 'PYUPDA DONE'
  32       RETURN
  33       END
  34       SUBROUTINE INITPYTHIA(GAGA)
  35       IMPLICIT NONE
  36       INTEGER PYCOMP
  37       EXTERNAL PYDATA,PYCOMP
  38       LOGICAL FIRST
  39       SAVE FIRST
  40       DATA FIRST / .TRUE. /
  41       DOUBLE PRECISION MAXIMUM
  42       COMMON/CBBEAM/MAXIMUM
  43       INTEGER IDC,GAGA
  44 C...Decay information.
  45       INTEGER MDCY,MDME,KFDP
  46       DOUBLE PRECISION BRAT
  47       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
  48
  49       CHARACTER CH1*16,CH2*16
  50
  51       INTEGER DC,CDC
  52 C     order of DC is d u s c b t b' t' ...
  53 C     ... e nu_e mu nu_mu tau nu_tau tau' nu_tau'
  54       COMMON/DECAYCH/DC(18)
  55
  56       IF(FIRST) THEN
  57          MAXIMUM=10.6
  58          WRITE (*,*) "Decay of Z0/gamma to"
  59 C...  Only allow Z0 decay to quarks (i.e. no leptonic final states).
  60 C     RS veto also b quarks
  61          DO IDC=MDCY(23,2),MDCY(23,2)+MDCY(23,3)-1
  62 C     first switch channel off
  63             MDME(IDC,1)=MIN(0,MDME(IDC,1))
  64             CDC=IABS(KFDP(IDC,1))
  65 C     then switch on again if requested
  66             CALL PYNAME(CDC,CH1)
  67             CALL PYNAME(-CDC,CH2)
  68             IF(DC(CDC).EQ.1) THEN
  69                MDME(IDC,1)=1
  70                WRITE (*,*) "  ",ch1,"+ ",ch2," ALLOWED"
  71             ELSE
  72                WRITE (*,*) "  ",ch1,"+ ",ch2," DISABLED"
  73             ENDIF
  74          ENDDO
  75 C     CALL PYSTAT(2)
  76
  77          IF(GAGA.EQ.0) THEN
  78             CALL PYINIT('CMS','E+','E-',MAXIMUM)
  79          ELSE
  80             CALL PYINIT('CMS','GAMMA/E+','GAMMA/E-',MAXIMUM)
  81          ENDIF
  82 C     CALL PYLIST(0)
  83 C         WRITE (*,*) 'done PYTHIA initialization'//
  84 C     >        ' with varying beam energy'
  85 C         WRITE (*,*) 'maximum allowed energy is ',MAXIMUM,' GeV'
  86          FIRST=.FALSE.
  87       ENDIF
  88       RETURN
  89       END
  90
  91 C--------------------------------------------------------------------------
  92 C
  93 C Environment:
  94 C      This software is part of the EvtGen package developed jointly
  95 C      for the BaBar and CLEO collaborations.  If you use all or part
  96 C      of it, please give an appropriate acknowledgement.
  97 C
  98 C Copyright Information: See BelEvtGen/COPYRIGHT
  99 C      Copyright (C) 1998      Caltech, UCSB
 100 C
 101 C Module: continuum.F
 102 C
 103 C Description:
 104 C
 105 C Modification history:
 106 C
 107 C    DJL/RYD     August 11, 1998         Module created
 108 C     26-Sep-2002 - RS : changed to PYTHIA
 109 C
 110 C------------------------------------------------------------------------
 111       SUBROUTINE PYCONTINUUM(ENERGY,NDAUG,KF,PX,PY,PZ,E)
 112       IMPLICIT NONE
 113       DOUBLE PRECISION MAXIMUM
 114       COMMON/CBBEAM/MAXIMUM
 115
 116       DOUBLE PRECISION P,V
 117       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
 118
 119       INTEGER MSTU,MSTJ
 120       DOUBLE PRECISION PARU,PARJ
 121       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
 122
 123       INTEGER MSTP,MSTI
 124       DOUBLE PRECISION PARP,PARI
 125       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
 126 C...Decay information.
 127       INTEGER MDCY,MDME,KFDP
 128       DOUBLE PRECISION BRAT
 129       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
 130
 131       INTEGER N,NPAD,K,IMSTJ
 132
 133       REAL*8 PXSUM,PYSUM,PZSUM
 134
 135       REAL*8 ENERGY
 136
 137       INTEGER PYCOMP
 138       EXTERNAL PYCOMP
 139
 140       INTEGER KF(100),I,NDAUG
 141       REAL*8  PX(100),PY(100),PZ(100),E(100)
 142
 143 C     RS Particles should not decay at this stage
 144 C     but remember the previous setting
 145       IMSTJ=MSTJ(21)
 146       MSTJ(21)=0
 147
 148 C      IF(FLAVOR.NE.LastProdFLav) THEN
 149 C...  Only allow Z0 decay to decay to certain flavour
 150 C         DO IDC=MDCY(23,2),MDCY(23,2)+MDCY(23,3)-1
 151 C            IF(IABS(KFDP(IDC,1)).NE.Flavor)
 152 C     >           MDME(IDC,1)=MIN(0,MDME(IDC,1))
 153 C         ENDDO
 154 C      ENDIF
 155
 156  4    PARP(171)=Energy/Maximum
 157       CALL PYEVNT
 158
 159 C--- only primary quarks and stable particles remain in the event record
 160       CALL PYEDIT(5)
 161
 162 C     RS and allow decays
 163       MSTJ(21)=IMSTJ
 164
 165       NDAUG=0
 166 C     sum to check that we preserve momentum
 167       PXSUM=0.0
 168       PYSUM=0.0
 169       PZSUM=0.0
 170
 171       DO 1,I=1,N
 172          NDAUG=NDAUG+1
 173          KF(NDAUG)=K(I,2)
 174          PX(NDAUG)=P(I,1)
 175          PY(NDAUG)=P(I,2)
 176          PZ(NDAUG)=P(I,3)
 177          E(NDAUG) =P(I,4)
 178          IF(K(I,1).LT.10)THEN
 179             PXSUM=PXSUM+PX(NDAUG)
 180             PYSUM=PYSUM+PY(NDAUG)
 181             PZSUM=PZSUM+PZ(NDAUG)
 182          ENDIF
 183 1     CONTINUE
 184
 185       IF (ABS(PXSUM).GT.0.001.OR.
 186      +       ABS(PYSUM).GT.0.001.OR.
 187      +       ABS(PZSUM).GT.0.001) THEN
 188
 189
 190          PRINT *, 'Momentum not conserved in jetset fragmentation:'
 191          PRINT *,'dPx:',pxsum,' dPy:',pysum,' dPz:',pzsum
 192
 193          CALL PYLIST(2)
 194
 195          GOTO 4
 196
 197       ENDIF
 198       RETURN
 199       END
 200 C--------------------------------------------------------------------------
 201 C
 202 C Environment:
 203 C      This software is part of the EvtGen package developed jointly
 204 C      for the BaBar and CLEO collaborations.  If you use all or part
 205 C      of it, please give an appropriate acknowledgement.
 206 C
 207 C Copyright Information: See BelEvtGen/COPYRIGHT
 208 C      Copyright (C) 1998      Caltech, UCSB
 209 C
 210 C Module: jetset1.F
 211 C
 212 C Description:
 213 C
 214 C Modification history:
 215 C
 216 C    DJL/RYD     August 11, 1998         Module created
 217 C     26-Sep-2002 - RS : changed to PYTHIA
 218 C
 219 C------------------------------------------------------------------------
 220       SUBROUTINE PYTHIADEC(IP,M,NDAUG,KF,KM,PX,PY,PZ,E)
 221
 222 C
 223 C interface to JETSET 7.4 to have one particle decayed
 224 C including possibly fragmentation, if the decay products include
 225 C partons.
 226 C
 227       IMPLICIT NONE
 228
 229       INTEGER MSTU,MSTJ
 230       DOUBLE PRECISION PARU,PARJ
 231       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
 232       SAVE /PYDAT1/
 233       DOUBLE PRECISION P,V
 234       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
 235
 236       INTEGER N,NPAD,K
 237
 238       INTEGER NDMAX
 239       PARAMETER (NDMAX=100)
 240
 241       INTEGER PYCOMP
 242       EXTERNAL PYCOMP
 243
 244       INTEGER IP,KF(NDMAX),I,NDAUG,KM(NDMAX)
 245       INTEGER KP,KID,IPART1,IPART
 246       REAL*8  M,PX(NDMAX),PY(NDMAX),PZ(NDMAX),E(NDMAX)
 247
 248       INTEGER IP1,IMSTJ
 249       DOUBLE PRECISION QMAX
 250
 251 C      CALL PY1ENT(1,IP,M,0D0,0D0)
 252
 253       K(1,1)=1
 254       K(1,2)=ip
 255       P(1,5)=m
 256       P(1,4)=m
 257       P(1,1)=0.0
 258       P(1,2)=0.0
 259       P(1,3)=0.0
 260       n=1
 261
 262 C     RS Particles should not decay at this stage
 263 C     but remember the previous setting
 264       IMSTJ=MSTJ(21)
 265       MSTJ(21)=0
 266
 267       CALL PYDECY(1)
 268
 269 C     RS and allow decays
 270       MSTJ(21)=IMSTJ
 271
 272
 273
 274
 275 C code copied from LUEXEC to avoid error with shower
 276 C switched on
 277 C...Decay products may develop a shower.
 278         IF(MSTJ(92).GT.0) THEN
 279           IP1=MSTJ(92)
 280           QMAX=SQRT(MAX(0.,(P(IP1,4)+P(IP1+1,4))**2-(P(IP1,1)+P(IP1+1,
 281      &    1))**2-(P(IP1,2)+P(IP1+1,2))**2-(P(IP1,3)+P(IP1+1,3))**2))
 282           CALL PYSHOW(IP1,IP1+1,QMAX)
 283           CALL PYPREP(IP1)
 284           MSTJ(92)=0
 285         ELSEIF(MSTJ(92).LT.0) THEN
 286           IP1=-MSTJ(92)
 287           CALL PYSHOW(IP1,-3,P(IP,5))
 288           CALL PYPREP(IP1)
 289           MSTJ(92)=0
 290         ENDIF
 291
 292 c
 293 c for debugging:
 294 c      call lulist(1)
 295 c
 296       mstj(21)=0
 297       call pyexec
 298       mstj(21)=2
 299
 300 c find partons, delete secondary partons, set mother pointers
 301
 302       ndaug = 0
 303       ipart1 = 1
 304       ipart = 1
 305
 306       do 10 i=2,n
 307          kp = k(i,3)
 308          kid = k(i,2)
 309          if (abs(kid) .ge. 1 .and. abs(kid) .le. 8
 310      1        .or. kid .eq. 21
 311      2        .or. kid .ge. 91 .and. kid .le. 94)  then
 312             if (ipart1 .eq. 1) ipart1 = i
 313             ipart = i
 314             if (kp .ne. 1)  goto 10
 315             kp = 0
 316          else
 317             if (kp .gt. ipart)  then
 318                goto 10
 319             elseif (kp .ge. ipart1)  then
 320                kp = ipart1-1
 321             else
 322                kp = 0
 323             endif
 324          endif
 325          ndaug = ndaug + 1
 326          km(ndaug)=kp
 327          kf(ndaug)=kid
 328          px(ndaug)=p(i,1)
 329          py(ndaug)=p(i,2)
 330          pz(ndaug)=p(i,3)
 331          e(ndaug)=p(i,4)
 332
 333 C        print '( 2I5,I12,4F12.4 )',ndaug,km(ndaug),kf(ndaug),
 334 C     1    px(ndaug),py(ndaug),pz(ndaug),e(ndaug)
 335
 336  10   CONTINUE
 337
 338 C      NDAUG=0
 339 C      DO I=2,N
 340 C        NDAUG = NDAUG + 1
 341 C        KM(NDAUG)=K(I,3)
 342 C        KF(NDAUG)=K(I,2)
 343 C        PX(NDAUG)=P(I,1)
 344 C        PY(NDAUG)=P(I,2)
 345 C        PZ(NDAUG)=P(I,3)
 346 C        E(NDAUG)=P(I,4)
 347 C      ENDDO
 348
 349
 350       RETURN
 351       END