pythiaeRHIC/pythia/pytest.f

0001
0002 C*********************************************************************
0003
0004 C...PYTEST
0005 C...A simple program (disguised as subroutine) to run at installation
0006 C...as a check that the program works as intended.
0007
0008       SUBROUTINE PYTEST(MTEST)
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/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
0019       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0020       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0021       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYDAT3/,/PYSUBS/,/PYPARS/
0022 C...Local arrays.
0023       DIMENSION PSUM(5),PINI(6),PFIN(6)
0024
0025 C...Save defaults for values that are changed.
0026       MSTJ1=MSTJ(1)
0027       MSTJ3=MSTJ(3)
0028       MSTJ11=MSTJ(11)
0029       MSTJ42=MSTJ(42)
0030       MSTJ43=MSTJ(43)
0031       MSTJ44=MSTJ(44)
0032       PARJ17=PARJ(17)
0033       PARJ22=PARJ(22)
0034       PARJ43=PARJ(43)
0035       PARJ54=PARJ(54)
0036       MST101=MSTJ(101)
0037       MST104=MSTJ(104)
0038       MST105=MSTJ(105)
0039       MST107=MSTJ(107)
0040       MST116=MSTJ(116)
0041
0042 C...First part: loop over simple events to be generated.
0043       IF(MTEST.GE.1) CALL PYTABU(20)
0044       NERR=0
0045       DO 180 IEV=1,500
0046
0047 C...Reset parameter values. Switch on some nonstandard features.
0048         MSTJ(1)=1
0049         MSTJ(3)=0
0050         MSTJ(11)=1
0051         MSTJ(42)=2
0052         MSTJ(43)=4
0053         MSTJ(44)=2
0054         PARJ(17)=0.1D0
0055         PARJ(22)=1.5D0
0056         PARJ(43)=1D0
0057         PARJ(54)=-0.05D0
0058         MSTJ(101)=5
0059         MSTJ(104)=5
0060         MSTJ(105)=0
0061         MSTJ(107)=1
0062         IF(IEV.EQ.301.OR.IEV.EQ.351.OR.IEV.EQ.401) MSTJ(116)=3
0063
0064 C...Ten events each for some single jets configurations.
0065         IF(IEV.LE.50) THEN
0066           ITY=(IEV+9)/10
0067           MSTJ(3)=-1
0068           IF(ITY.EQ.3.OR.ITY.EQ.4) MSTJ(11)=2
0069           IF(ITY.EQ.1) CALL PY1ENT(1,1,15D0,0D0,0D0)
0070           IF(ITY.EQ.2) CALL PY1ENT(1,3101,15D0,0D0,0D0)
0071           IF(ITY.EQ.3) CALL PY1ENT(1,-2203,15D0,0D0,0D0)
0072           IF(ITY.EQ.4) CALL PY1ENT(1,-4,30D0,0D0,0D0)
0073           IF(ITY.EQ.5) CALL PY1ENT(1,21,15D0,0D0,0D0)
0074
0075 C...Ten events each for some simple jet systems; string fragmentation.
0076         ELSEIF(IEV.LE.130) THEN
0077           ITY=(IEV-41)/10
0078           IF(ITY.EQ.1) CALL PY2ENT(1,1,-1,40D0)
0079           IF(ITY.EQ.2) CALL PY2ENT(1,4,-4,30D0)
0080           IF(ITY.EQ.3) CALL PY2ENT(1,2,2103,100D0)
0081           IF(ITY.EQ.4) CALL PY2ENT(1,21,21,40D0)
0082           IF(ITY.EQ.5) CALL PY3ENT(1,2101,21,-3203,30D0,0.6D0,0.8D0)
0083           IF(ITY.EQ.6) CALL PY3ENT(1,5,21,-5,40D0,0.9D0,0.8D0)
0084           IF(ITY.EQ.7) CALL PY3ENT(1,21,21,21,60D0,0.7D0,0.5D0)
0085           IF(ITY.EQ.8) CALL PY4ENT(1,2,21,21,-2,40D0,
0086      &    0.4D0,0.64D0,0.6D0,0.12D0,0.2D0)
0087
0088 C...Seventy events with independent fragmentation and momentum cons.
0089         ELSEIF(IEV.LE.200) THEN
0090           ITY=1+(IEV-131)/16
0091           MSTJ(2)=1+MOD(IEV-131,4)
0092           MSTJ(3)=1+MOD((IEV-131)/4,4)
0093           IF(ITY.EQ.1) CALL PY2ENT(1,4,-5,40D0)
0094           IF(ITY.EQ.2) CALL PY3ENT(1,3,21,-3,40D0,0.9D0,0.4D0)
0095           IF(ITY.EQ.3) CALL PY4ENT(1,2,21,21,-2,40D0,
0096      &    0.4D0,0.64D0,0.6D0,0.12D0,0.2D0)
0097           IF(ITY.GE.4) CALL PY4ENT(1,2,-3,3,-2,40D0,
0098      &    0.4D0,0.64D0,0.6D0,0.12D0,0.2D0)
0099
0100 C...A hundred events with random jets (check invariant mass).
0101         ELSEIF(IEV.LE.300) THEN
0102   100     DO 110 J=1,5
0103             PSUM(J)=0D0
0104   110     CONTINUE
0105           NJET=2D0+6D0*PYR(0)
0106           DO 130 I=1,NJET
0107             KFL=21
0108             IF(I.EQ.1) KFL=INT(1D0+4D0*PYR(0))
0109             IF(I.EQ.NJET) KFL=-INT(1D0+4D0*PYR(0))
0110             EJET=5D0+20D0*PYR(0)
0111             THETA=ACOS(2D0*PYR(0)-1D0)
0112             PHI=6.2832D0*PYR(0)
0113             IF(I.LT.NJET) CALL PY1ENT(-I,KFL,EJET,THETA,PHI)
0114             IF(I.EQ.NJET) CALL PY1ENT(I,KFL,EJET,THETA,PHI)
0115             IF(I.EQ.1.OR.I.EQ.NJET) MSTJ(93)=1
0116             IF(I.EQ.1.OR.I.EQ.NJET) PSUM(5)=PSUM(5)+PYMASS(KFL)
0117             DO 120 J=1,4
0118               PSUM(J)=PSUM(J)+P(I,J)
0119   120       CONTINUE
0120   130     CONTINUE
0121           IF(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2.LT.
0122      &    (PSUM(5)+PARJ(32))**2) GOTO 100
0123
0124 C...Fifty e+e- continuum events with matrix elements.
0125         ELSEIF(IEV.LE.350) THEN
0126           MSTJ(101)=2
0127           CALL PYEEVT(0,40D0)
0128
0129 C...Fifty e+e- continuum event with varying shower options.
0130         ELSEIF(IEV.LE.400) THEN
0131           MSTJ(42)=1+MOD(IEV,2)
0132           MSTJ(43)=1+MOD(IEV/2,4)
0133           MSTJ(44)=MOD(IEV/8,3)
0134           CALL PYEEVT(0,90D0)
0135
0136 C...Fifty e+e- continuum events with coherent shower.
0137         ELSEIF(IEV.LE.450) THEN
0138           CALL PYEEVT(0,500D0)
0139
0140 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
0141         ELSE
0142           CALL PYONIA(5,9.46D0)
0143         ENDIF
0144
0145 C...Generate event. Find total momentum, energy and charge.
0146         DO 140 J=1,4
0147           PINI(J)=PYP(0,J)
0148   140   CONTINUE
0149         PINI(6)=PYP(0,6)
0150         CALL PYEXEC
0151         DO 150 J=1,4
0152           PFIN(J)=PYP(0,J)
0153   150   CONTINUE
0154         PFIN(6)=PYP(0,6)
0155
0156 C...Check conservation of energy, momentum and charge;
0157 C...usually exact, but only approximate for single jets.
0158         MERR=0
0159         IF(IEV.LE.50) THEN
0160           IF((PFIN(1)-PINI(1))**2+(PFIN(2)-PINI(2))**2.GE.10D0)
0161      &    MERR=MERR+1
0162           EPZREM=PINI(4)+PINI(3)-PFIN(4)-PFIN(3)
0163           IF(EPZREM.LT.0D0.OR.EPZREM.GT.2D0*PARJ(31)) MERR=MERR+1
0164           IF(ABS(PFIN(6)-PINI(6)).GT.2.1D0) MERR=MERR+1
0165         ELSE
0166           DO 160 J=1,4
0167             IF(ABS(PFIN(J)-PINI(J)).GT.0.0001D0*PINI(4)) MERR=MERR+1
0168   160     CONTINUE
0169           IF(ABS(PFIN(6)-PINI(6)).GT.0.1D0) MERR=MERR+1
0170         ENDIF
0171         IF(MERR.NE.0) WRITE(MSTU(11),5000) (PINI(J),J=1,4),PINI(6),
0172      &  (PFIN(J),J=1,4),PFIN(6)
0173
0174 C...Check that all KF codes are known ones, and that partons/particles
0175 C...satisfy energy-momentum-mass relation. Store particle statistics.
0176         DO 170 I=1,N
0177           IF(K(I,1).GT.20) GOTO 170
0178           IF(PYCOMP(K(I,2)).EQ.0) THEN
0179             WRITE(MSTU(11),5100) I
0180             MERR=MERR+1
0181           ENDIF
0182           PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
0183           IF(ABS(PD).GT.MAX(0.1D0,0.001D0*P(I,4)**2).OR.P(I,4).LT.0D0)
0184      &    THEN
0185             WRITE(MSTU(11),5200) I
0186             MERR=MERR+1
0187           ENDIF
0188   170   CONTINUE
0189         IF(MTEST.GE.1) CALL PYTABU(21)
0190
0191 C...List all erroneous events and some normal ones.
0192         IF(MERR.NE.0.OR.MSTU(24).NE.0.OR.MSTU(28).NE.0) THEN
0193           IF(MERR.GE.1) WRITE(MSTU(11),6400)
0194           CALL PYLIST(2)
0195         ELSEIF(MTEST.GE.1.AND.MOD(IEV-5,100).EQ.0) THEN
0196           CALL PYLIST(1)
0197         ENDIF
0198
0199 C...Stop execution if too many errors.
0200         IF(MERR.NE.0) NERR=NERR+1
0201         IF(NERR.GE.10) THEN
0202           WRITE(MSTU(11),6300)
0203           CALL PYLIST(1)
0204           CALL PYSTOP(9)
0205         ENDIF
0206   180 CONTINUE
0207
0208 C...Summarize result of run.
0209       IF(MTEST.GE.1) CALL PYTABU(22)
0210
0211 C...Reset commonblock variables changed during run.
0212       MSTJ(1)=MSTJ1
0213       MSTJ(3)=MSTJ3
0214       MSTJ(11)=MSTJ11
0215       MSTJ(42)=MSTJ42
0216       MSTJ(43)=MSTJ43
0217       MSTJ(44)=MSTJ44
0218       PARJ(17)=PARJ17
0219       PARJ(22)=PARJ22
0220       PARJ(43)=PARJ43
0221       PARJ(54)=PARJ54
0222       MSTJ(101)=MST101
0223       MSTJ(104)=MST104
0224       MSTJ(105)=MST105
0225       MSTJ(107)=MST107
0226       MSTJ(116)=MST116
0227
0228 C...Second part: complete events of various kinds.
0229 C...Common initial values. Loop over initiating conditions.
0230       MSTP(122)=MAX(0,MIN(2,MTEST))
0231       MDCY(PYCOMP(111),1)=0
0232       DO 230 IPROC=1,8
0233
0234 C...Reset process type, kinematics cuts, and the flags used.
0235         MSEL=0
0236         DO 190 ISUB=1,500
0237           MSUB(ISUB)=0
0238   190   CONTINUE
0239         CKIN(1)=2D0
0240         CKIN(3)=0D0
0241         MSTP(2)=1
0242         MSTP(11)=0
0243         MSTP(33)=0
0244         MSTP(81)=1
0245         MSTP(82)=1
0246         MSTP(111)=1
0247         MSTP(131)=0
0248         MSTP(133)=0
0249         PARP(131)=0.01D0
0250
0251 C...Prompt photon production at fixed target.
0252         IF(IPROC.EQ.1) THEN
0253           PZSUM=300D0
0254           PESUM=SQRT(PZSUM**2+PYMASS(211)**2)+PYMASS(2212)
0255           PQSUM=2D0
0256           MSEL=10
0257           CKIN(3)=5D0
0258           CALL PYINIT('FIXT','pi+','p',PZSUM)
0259
0260 C...QCD processes at ISR energies.
0261         ELSEIF(IPROC.EQ.2) THEN
0262           PESUM=63D0
0263           PZSUM=0D0
0264           PQSUM=2D0
0265           MSEL=1
0266           CKIN(3)=5D0
0267           CALL PYINIT('CMS','p','p',PESUM)
0268
0269 C...W production + multiple interactions at CERN Collider.
0270         ELSEIF(IPROC.EQ.3) THEN
0271           PESUM=630D0
0272           PZSUM=0D0
0273           PQSUM=0D0
0274           MSEL=12
0275           CKIN(1)=20D0
0276           MSTP(82)=4
0277           MSTP(2)=2
0278           MSTP(33)=3
0279           CALL PYINIT('CMS','p','pbar',PESUM)
0280
0281 C...W/Z gauge boson pairs + pileup events at the Tevatron.
0282         ELSEIF(IPROC.EQ.4) THEN
0283           PESUM=1800D0
0284           PZSUM=0D0
0285           PQSUM=0D0
0286           MSUB(22)=1
0287           MSUB(23)=1
0288           MSUB(25)=1
0289           CKIN(1)=200D0
0290           MSTP(111)=0
0291           MSTP(131)=1
0292           MSTP(133)=2
0293           PARP(131)=0.04D0
0294           CALL PYINIT('CMS','p','pbar',PESUM)
0295
0296 C...Higgs production at LHC.
0297         ELSEIF(IPROC.EQ.5) THEN
0298           PESUM=15400D0
0299           PZSUM=0D0
0300           PQSUM=2D0
0301           MSUB(3)=1
0302           MSUB(102)=1
0303           MSUB(123)=1
0304           MSUB(124)=1
0305           PMAS(25,1)=300D0
0306           CKIN(1)=200D0
0307           MSTP(81)=0
0308           MSTP(111)=0
0309           CALL PYINIT('CMS','p','p',PESUM)
0310
0311 C...Z' production at SSC.
0312         ELSEIF(IPROC.EQ.6) THEN
0313           PESUM=40000D0
0314           PZSUM=0D0
0315           PQSUM=2D0
0316           MSEL=21
0317           PMAS(32,1)=600D0
0318           CKIN(1)=400D0
0319           MSTP(81)=0
0320           MSTP(111)=0
0321           CALL PYINIT('CMS','p','p',PESUM)
0322
0323 C...W pair production at 1 TeV e+e- collider.
0324         ELSEIF(IPROC.EQ.7) THEN
0325           PESUM=1000D0
0326           PZSUM=0D0
0327           PQSUM=0D0
0328           MSUB(25)=1
0329           MSUB(69)=1
0330           MSTP(11)=1
0331           CALL PYINIT('CMS','e+','e-',PESUM)
0332
0333 C...Deep inelastic scattering at a LEP+LHC ep collider.
0334         ELSEIF(IPROC.EQ.8) THEN
0335           P(1,1)=0D0
0336           P(1,2)=0D0
0337           P(1,3)=8000D0
0338           P(2,1)=0D0
0339           P(2,2)=0D0
0340           P(2,3)=-80D0
0341           PESUM=8080D0
0342           PZSUM=7920D0
0343           PQSUM=0D0
0344           MSUB(10)=1
0345           CKIN(3)=50D0
0346           MSTP(111)=0
0347           CALL PYINIT('3MOM','p','e-',PESUM)
0348         ENDIF
0349
0350 C...Generate 20 events of each required type.
0351         DO 220 IEV=1,20
0352           CALL PYEVNT
0353           PESUMM=PESUM
0354           IF(IPROC.EQ.4) PESUMM=MSTI(41)*PESUM
0355
0356 C...Check conservation of energy/momentum/flavour.
0357           PINI(1)=0D0
0358           PINI(2)=0D0
0359           PINI(3)=PZSUM
0360           PINI(4)=PESUMM
0361           PINI(6)=PQSUM
0362           DO 200 J=1,4
0363             PFIN(J)=PYP(0,J)
0364   200     CONTINUE
0365           PFIN(6)=PYP(0,6)
0366           MERR=0
0367           DEVE=ABS(PFIN(4)-PINI(4))+ABS(PFIN(3)-PINI(3))
0368           DEVT=ABS(PFIN(1)-PINI(1))+ABS(PFIN(2)-PINI(2))
0369           DEVQ=ABS(PFIN(6)-PINI(6))
0370           IF(DEVE.GT.2D-3*PESUM.OR.DEVT.GT.MAX(0.01D0,1D-4*PESUM).OR.
0371      &    DEVQ.GT.0.1D0) MERR=1
0372           IF(MERR.NE.0) WRITE(MSTU(11),5000) (PINI(J),J=1,4),PINI(6),
0373      &    (PFIN(J),J=1,4),PFIN(6)
0374
0375 C...Check that all KF codes are known ones, and that partons/particles
0376 C...satisfy energy-momentum-mass relation.
0377           DO 210 I=1,N
0378             IF(K(I,1).GT.20) GOTO 210
0379             IF(PYCOMP(K(I,2)).EQ.0) THEN
0380               WRITE(MSTU(11),5100) I
0381               MERR=MERR+1
0382             ENDIF
0383             PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2*
0384      &      SIGN(1D0,P(I,5))
0385             IF(ABS(PD).GT.MAX(0.1D0,0.002D0*P(I,4)**2,0.002D0*P(I,5)**2)
0386      &      .OR.(P(I,5).GE.0D0.AND.P(I,4).LT.0D0)) THEN
0387               WRITE(MSTU(11),5200) I
0388               MERR=MERR+1
0389             ENDIF
0390   210     CONTINUE
0391
0392 C...Listing of erroneous events, and first event of each type.
0393           IF(MERR.GE.1) NERR=NERR+1
0394           IF(NERR.GE.10) THEN
0395             WRITE(MSTU(11),6300)
0396             CALL PYLIST(1)
0397             CALL PYSTOP(9)
0398           ENDIF
0399           IF(MTEST.GE.1.AND.(MERR.GE.1.OR.IEV.EQ.1)) THEN
0400             IF(MERR.GE.1) WRITE(MSTU(11),6400)
0401             CALL PYLIST(1)
0402           ENDIF
0403   220   CONTINUE
0404
0405 C...List statistics for each process type.
0406         IF(MTEST.GE.1) CALL PYSTAT(1)
0407   230 CONTINUE
0408
0409 C...Summarize result of run.
0410       IF(NERR.EQ.0) WRITE(MSTU(11),6500)
0411       IF(NERR.GT.0) WRITE(MSTU(11),6600) NERR
0412
0413 C...Format statements for output.
0414  5000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
0415      &'in following event'/' sum of',9X,'px',11X,'py',11X,'pz',11X,
0416      &'E',8X,'charge'/' before',2X,4(1X,F12.5),1X,F8.2/' after',3X,
0417      &4(1X,F12.5),1X,F8.2)
0418  5100 FORMAT(/5X,'Entry no.',I4,' in following event not known code')
0419  5200 FORMAT(/5X,'Entry no.',I4,' in following event has faulty ',
0420      &'kinematics')
0421  6300 FORMAT(/5X,'This is the tenth error experienced! Something is ',
0422      &'wrong.'/5X,'Execution will be stopped after listing of event.')
0423  6400 FORMAT(5X,'Faulty event follows:')
0424  6500 FORMAT(//5X,'End result of PYTEST: no errors detected.')
0425  6600 FORMAT(//5X,'End result of PYTEST:',I2,' errors detected.'/
0426      &5X,'This should not have happened!')
0427
0428       RETURN
0429       END