PYTHIA/jetset/lutest.F

   1
   2 C*********************************************************************
   3
   4       SUBROUTINE LUTEST(MTEST)
   5
   6 C...Purpose: to provide a simple program (disguised as subroutine) to
   7 C...run at installation as a check that the program works as intended.
   8       COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
   9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  10       SAVE /LUJETS/,/LUDAT1/
  11       DIMENSION PSUM(5),PINI(6),PFIN(6)
  12
  13 C...Loop over events to be generated.
  14       IF(MTEST.GE.1) CALL LUTABU(20)
  15       NERR=0
  16       DO 180 IEV=1,600
  17
  18 C...Reset parameter values. Switch on some nonstandard features.
  19       MSTJ(1)=1
  20       MSTJ(3)=0
  21       MSTJ(11)=1
  22       MSTJ(42)=2
  23       MSTJ(43)=4
  24       MSTJ(44)=2
  25       PARJ(17)=0.1
  26       PARJ(22)=1.5
  27       PARJ(43)=1.
  28       PARJ(54)=-0.05
  29       MSTJ(101)=5
  30       MSTJ(104)=5
  31       MSTJ(105)=0
  32       MSTJ(107)=1
  33       IF(IEV.EQ.301.OR.IEV.EQ.351.OR.IEV.EQ.401) MSTJ(116)=3
  34
  35 C...Ten events each for some single jets configurations.
  36       IF(IEV.LE.50) THEN
  37         ITY=(IEV+9)/10
  38         MSTJ(3)=-1
  39         IF(ITY.EQ.3.OR.ITY.EQ.4) MSTJ(11)=2
  40         IF(ITY.EQ.1) CALL LU1ENT(1,1,15.,0.,0.)
  41         IF(ITY.EQ.2) CALL LU1ENT(1,3101,15.,0.,0.)
  42         IF(ITY.EQ.3) CALL LU1ENT(1,-2203,15.,0.,0.)
  43         IF(ITY.EQ.4) CALL LU1ENT(1,-4,30.,0.,0.)
  44         IF(ITY.EQ.5) CALL LU1ENT(1,21,15.,0.,0.)
  45
  46 C...Ten events each for some simple jet systems; string fragmentation.
  47       ELSEIF(IEV.LE.130) THEN
  48         ITY=(IEV-41)/10
  49         IF(ITY.EQ.1) CALL LU2ENT(1,1,-1,40.)
  50         IF(ITY.EQ.2) CALL LU2ENT(1,4,-4,30.)
  51         IF(ITY.EQ.3) CALL LU2ENT(1,2,2103,100.)
  52         IF(ITY.EQ.4) CALL LU2ENT(1,21,21,40.)
  53         IF(ITY.EQ.5) CALL LU3ENT(1,2101,21,-3203,30.,0.6,0.8)
  54         IF(ITY.EQ.6) CALL LU3ENT(1,5,21,-5,40.,0.9,0.8)
  55         IF(ITY.EQ.7) CALL LU3ENT(1,21,21,21,60.,0.7,0.5)
  56         IF(ITY.EQ.8) CALL LU4ENT(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)
  57
  58 C...Seventy events with independent fragmentation and momentum cons.
  59       ELSEIF(IEV.LE.200) THEN
  60         ITY=1+(IEV-131)/16
  61         MSTJ(2)=1+MOD(IEV-131,4)
  62         MSTJ(3)=1+MOD((IEV-131)/4,4)
  63         IF(ITY.EQ.1) CALL LU2ENT(1,4,-5,40.)
  64         IF(ITY.EQ.2) CALL LU3ENT(1,3,21,-3,40.,0.9,0.4)
  65         IF(ITY.EQ.3) CALL LU4ENT(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)
  66         IF(ITY.GE.4) CALL LU4ENT(1,2,-3,3,-2,40.,0.4,0.64,0.6,0.12,0.2)
  67
  68 C...A hundred events with random jets (check invariant mass).
  69       ELSEIF(IEV.LE.300) THEN
  70   100   DO 110 J=1,5
  71         PSUM(J)=0.
  72   110   CONTINUE
  73         NJET=2.+6.*RLU(0)
  74         DO 130 I=1,NJET
  75         KFL=21
  76         IF(I.EQ.1) KFL=INT(1.+4.*RLU(0))
  77         IF(I.EQ.NJET) KFL=-INT(1.+4.*RLU(0))
  78         EJET=5.+20.*RLU(0)
  79         THETA=ACOS(2.*RLU(0)-1.)
  80         PHI=6.2832*RLU(0)
  81         IF(I.LT.NJET) CALL LU1ENT(-I,KFL,EJET,THETA,PHI)
  82         IF(I.EQ.NJET) CALL LU1ENT(I,KFL,EJET,THETA,PHI)
  83         IF(I.EQ.1.OR.I.EQ.NJET) MSTJ(93)=1
  84         IF(I.EQ.1.OR.I.EQ.NJET) PSUM(5)=PSUM(5)+ULMASS(KFL)
  85         DO 120 J=1,4
  86         PSUM(J)=PSUM(J)+P(I,J)
  87   120   CONTINUE
  88   130   CONTINUE
  89         IF(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2.LT.
  90      &  (PSUM(5)+PARJ(32))**2) GOTO 100
  91
  92 C...Fifty e+e- continuum events with matrix elements.
  93       ELSEIF(IEV.LE.350) THEN
  94         MSTJ(101)=2
  95         CALL LUEEVT(0,40.)
  96
  97 C...Fifty e+e- continuum event with varying shower options.
  98       ELSEIF(IEV.LE.400) THEN
  99         MSTJ(42)=1+MOD(IEV,2)
 100         MSTJ(43)=1+MOD(IEV/2,4)
 101         MSTJ(44)=MOD(IEV/8,3)
 102         CALL LUEEVT(0,90.)
 103
 104 C...Fifty e+e- continuum events with coherent shower, including top.
 105       ELSEIF(IEV.LE.450) THEN
 106         MSTJ(104)=6
 107         CALL LUEEVT(0,500.)
 108
 109 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
 110       ELSEIF(IEV.LE.500) THEN
 111         CALL LUONIA(5,9.46)
 112
 113 C...One decay each for some heavy mesons.
 114       ELSEIF(IEV.LE.560) THEN
 115         ITY=IEV-501
 116         KFLS=2*(ITY/20)+1
 117         KFLB=8-MOD(ITY/5,4)
 118         KFLC=KFLB-MOD(ITY,5)
 119         CALL LU1ENT(1,100*KFLB+10*KFLC+KFLS,0.,0.,0.)
 120
 121 C...One decay each for some heavy baryons.
 122       ELSEIF(IEV.LE.600) THEN
 123         ITY=IEV-561
 124         KFLS=2*(ITY/20)+2
 125         KFLA=8-MOD(ITY/5,4)
 126         KFLB=KFLA-MOD(ITY,5)
 127         KFLC=MAX(1,KFLB-1)
 128         CALL LU1ENT(1,1000*KFLA+100*KFLB+10*KFLC+KFLS,0.,0.,0.)
 129       ENDIF
 130
 131 C...Generate event. Find total momentum, energy and charge.
 132       DO 140 J=1,4
 133       PINI(J)=PLU(0,J)
 134   140 CONTINUE
 135       PINI(6)=PLU(0,6)
 136       CALL LUEXEC
 137       DO 150 J=1,4
 138       PFIN(J)=PLU(0,J)
 139   150 CONTINUE
 140       PFIN(6)=PLU(0,6)
 141
 142 C...Check conservation of energy, momentum and charge;
 143 C...usually exact, but only approximate for single jets.
 144       MERR=0
 145       IF(IEV.LE.50) THEN
 146         IF((PFIN(1)-PINI(1))**2+(PFIN(2)-PINI(2))**2.GE.4.) MERR=MERR+1
 147         EPZREM=PINI(4)+PINI(3)-PFIN(4)-PFIN(3)
 148         IF(EPZREM.LT.0..OR.EPZREM.GT.2.*PARJ(31)) MERR=MERR+1
 149         IF(ABS(PFIN(6)-PINI(6)).GT.2.1) MERR=MERR+1
 150       ELSE
 151         DO 160 J=1,4
 152         IF(ABS(PFIN(J)-PINI(J)).GT.0.0001*PINI(4)) MERR=MERR+1
 153   160   CONTINUE
 154         IF(ABS(PFIN(6)-PINI(6)).GT.0.1) MERR=MERR+1
 155       ENDIF
 156       IF(MERR.NE.0) WRITE(MSTU(11),5000) (PINI(J),J=1,4),PINI(6),
 157      &(PFIN(J),J=1,4),PFIN(6)
 158
 159 C...Check that all KF codes are known ones, and that partons/particles
 160 C...satisfy energy-momentum-mass relation. Store particle statistics.
 161       DO 170 I=1,N
 162       IF(K(I,1).GT.20) GOTO 170
 163       IF(LUCOMP(K(I,2)).EQ.0) THEN
 164         WRITE(MSTU(11),5100) I
 165         MERR=MERR+1
 166       ENDIF
 167       PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
 168       IF(ABS(PD).GT.MAX(0.1,0.001*P(I,4)**2).OR.P(I,4).LT.0.) THEN
 169         WRITE(MSTU(11),5200) I
 170         MERR=MERR+1
 171       ENDIF
 172   170 CONTINUE
 173       IF(MTEST.GE.1) CALL LUTABU(21)
 174
 175 C...List all erroneous events and some normal ones.
 176       IF(MERR.NE.0.OR.MSTU(24).NE.0.OR.MSTU(28).NE.0) THEN
 177         CALL LULIST(2)
 178       ELSEIF(MTEST.GE.1.AND.MOD(IEV-5,100).EQ.0) THEN
 179         CALL LULIST(1)
 180       ENDIF
 181
 182 C...Stop execution if too many errors.
 183       IF(MERR.NE.0) NERR=NERR+1
 184       IF(NERR.GE.10) THEN
 185         WRITE(MSTU(11),5300) IEV
 186         STOP
 187       ENDIF
 188   180 CONTINUE
 189
 190 C...Summarize result of run.
 191       IF(MTEST.GE.1) CALL LUTABU(22)
 192       IF(NERR.EQ.0) WRITE(MSTU(11),5400)
 193       IF(NERR.GT.0) WRITE(MSTU(11),5500) NERR
 194
 195 C...Reset commonblock variables changed during run.
 196       MSTJ(2)=3
 197       PARJ(17)=0.
 198       PARJ(22)=1.
 199       PARJ(43)=0.5
 200       PARJ(54)=0.
 201       MSTJ(105)=1
 202       MSTJ(107)=0
 203
 204 C...Format statements for output.
 205  5000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
 206      &'in following event'/' sum of',9X,'px',11X,'py',11X,'pz',11X,
 207      &'E',8X,'charge'/' before',2X,4(1X,F12.5),1X,F8.2/' after',3X,
 208      &4(1X,F12.5),1X,F8.2)
 209  5100 FORMAT(/5X,'Entry no.',I4,' in following event not known code')
 210  5200 FORMAT(/5X,'Entry no.',I4,' in following event has faulty ',
 211      &'kinematics')
 212  5300 FORMAT(/5X,'Ten errors experienced by event ',I3/
 213      &5X,'Something is seriously wrong! Execution stopped now!')
 214  5400 FORMAT(//5X,'End result of LUTEST: no errors detected.')
 215  5500 FORMAT(//5X,'End result of LUTEST:',I2,' errors detected.'/
 216      &5X,'This should not have happened!')
 217
 218       RETURN
 219       END