hijing/src/lutest.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUTEST(MTEST)
0005
0006 C...Purpose: to provide a simple program (disguised as subroutine) to
0007 C...run at installation as a check that the program works as intended.
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       DIMENSION PSUM(5),PINI(6),PFIN(6)
0013
0014 C...Loop over events to be generated.
0015       IF(MTEST.GE.1) CALL LUTABU(20)
0016       NERR=0
0017       DO 170 IEV=1,600
0018
0019 C...Reset parameter values. Switch on some nonstandard features.
0020       MSTJ(1)=1
0021       MSTJ(3)=0
0022       MSTJ(11)=1
0023       MSTJ(42)=2
0024       MSTJ(43)=4
0025       MSTJ(44)=2
0026       PARJ(17)=0.1
0027       PARJ(22)=1.5
0028       PARJ(43)=1.
0029       PARJ(54)=-0.05
0030       MSTJ(101)=5
0031       MSTJ(104)=5
0032       MSTJ(105)=0
0033       MSTJ(107)=1
0034       IF(IEV.EQ.301.OR.IEV.EQ.351.OR.IEV.EQ.401) MSTJ(116)=3
0035
0036 C...Ten events each for some single jets configurations.
0037       IF(IEV.LE.50) THEN
0038         ITY=(IEV+9)/10
0039         MSTJ(3)=-1
0040         IF(ITY.EQ.3.OR.ITY.EQ.4) MSTJ(11)=2
0041         IF(ITY.EQ.1) CALL LU1ENT(1,1,15.,0.,0.)
0042         IF(ITY.EQ.2) CALL LU1ENT(1,3101,15.,0.,0.)
0043         IF(ITY.EQ.3) CALL LU1ENT(1,-2203,15.,0.,0.)
0044         IF(ITY.EQ.4) CALL LU1ENT(1,-4,30.,0.,0.)
0045         IF(ITY.EQ.5) CALL LU1ENT(1,21,15.,0.,0.)
0046
0047 C...Ten events each for some simple jet systems; string fragmentation.
0048       ELSEIF(IEV.LE.130) THEN
0049         ITY=(IEV-41)/10
0050         IF(ITY.EQ.1) CALL LU2ENT(1,1,-1,40.)
0051         IF(ITY.EQ.2) CALL LU2ENT(1,4,-4,30.)
0052         IF(ITY.EQ.3) CALL LU2ENT(1,2,2103,100.)
0053         IF(ITY.EQ.4) CALL LU2ENT(1,21,21,40.)
0054         IF(ITY.EQ.5) CALL LU3ENT(1,2101,21,-3203,30.,0.6,0.8)
0055         IF(ITY.EQ.6) CALL LU3ENT(1,5,21,-5,40.,0.9,0.8)
0056         IF(ITY.EQ.7) CALL LU3ENT(1,21,21,21,60.,0.7,0.5)
0057         IF(ITY.EQ.8) CALL LU4ENT(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)
0058
0059 C...Seventy events with independent fragmentation and momentum cons.
0060       ELSEIF(IEV.LE.200) THEN
0061         ITY=1+(IEV-131)/16
0062         MSTJ(2)=1+MOD(IEV-131,4)
0063         MSTJ(3)=1+MOD((IEV-131)/4,4)
0064         IF(ITY.EQ.1) CALL LU2ENT(1,4,-5,40.)
0065         IF(ITY.EQ.2) CALL LU3ENT(1,3,21,-3,40.,0.9,0.4)
0066         IF(ITY.EQ.3) CALL LU4ENT(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)
0067         IF(ITY.GE.4) CALL LU4ENT(1,2,-3,3,-2,40.,0.4,0.64,0.6,0.12,0.2)
0068
0069 C...A hundred events with random jets (check invariant mass).
0070       ELSEIF(IEV.LE.300) THEN
0071   100   DO 110 J=1,5
0072   110   PSUM(J)=0.
0073         NJET=2.+6.*RLU(0)
0074         DO 120 I=1,NJET
0075         KFL=21
0076         IF(I.EQ.1) KFL=INT(1.+4.*RLU(0))
0077         IF(I.EQ.NJET) KFL=-INT(1.+4.*RLU(0))
0078         EJET=5.+20.*RLU(0)
0079         THETA=ACOS(2.*RLU(0)-1.)
0080         PHI=6.2832*RLU(0)
0081         IF(I.LT.NJET) CALL LU1ENT(-I,KFL,EJET,THETA,PHI)
0082         IF(I.EQ.NJET) CALL LU1ENT(I,KFL,EJET,THETA,PHI)
0083         IF(I.EQ.1.OR.I.EQ.NJET) PSUM(5)=PSUM(5)+ULMASS(KFL)
0084         DO 120 J=1,4
0085   120   PSUM(J)=PSUM(J)+P(I,J)
0086         IF(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2.LT.
0087      &  (PSUM(5)+PARJ(32))**2) GOTO 100
0088
0089 C...Fifty e+e- continuum events with matrix elements.
0090       ELSEIF(IEV.LE.350) THEN
0091         MSTJ(101)=2
0092         CALL LUEEVT(0,40.)
0093
0094 C...Fifty e+e- continuum event with varying shower options.
0095       ELSEIF(IEV.LE.400) THEN
0096         MSTJ(42)=1+MOD(IEV,2)
0097         MSTJ(43)=1+MOD(IEV/2,4)
0098         MSTJ(44)=MOD(IEV/8,3)
0099         CALL LUEEVT(0,90.)
0100
0101 C...Fifty e+e- continuum events with coherent shower, including top.
0102       ELSEIF(IEV.LE.450) THEN
0103         MSTJ(104)=6
0104         CALL LUEEVT(0,500.)
0105
0106 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
0107       ELSEIF(IEV.LE.500) THEN
0108         CALL LUONIA(5,9.46)
0109
0110 C...One decay each for some heavy mesons.
0111       ELSEIF(IEV.LE.560) THEN
0112         ITY=IEV-501
0113         KFLS=2*(ITY/20)+1
0114         KFLB=8-MOD(ITY/5,4)
0115         KFLC=KFLB-MOD(ITY,5)
0116         CALL LU1ENT(1,100*KFLB+10*KFLC+KFLS,0.,0.,0.)
0117
0118 C...One decay each for some heavy baryons.
0119       ELSEIF(IEV.LE.600) THEN
0120         ITY=IEV-561
0121         KFLS=2*(ITY/20)+2
0122         KFLA=8-MOD(ITY/5,4)
0123         KFLB=KFLA-MOD(ITY,5)
0124         KFLC=MAX(1,KFLB-1)
0125         CALL LU1ENT(1,1000*KFLA+100*KFLB+10*KFLC+KFLS,0.,0.,0.)
0126       ENDIF
0127
0128 C...Generate event. Find total momentum, energy and charge.
0129       DO 130 J=1,4
0130   130 PINI(J)=PLU(0,J)
0131       PINI(6)=PLU(0,6)
0132       CALL LUEXEC
0133       DO 140 J=1,4
0134   140 PFIN(J)=PLU(0,J)
0135       PFIN(6)=PLU(0,6)
0136
0137 C...Check conservation of energy, momentum and charge;
0138 C...usually exact, but only approximate for single jets.
0139       MERR=0
0140       IF(IEV.LE.50) THEN
0141         IF((PFIN(1)-PINI(1))**2+(PFIN(2)-PINI(2))**2.GE.4.) MERR=MERR+1
0142         EPZREM=PINI(4)+PINI(3)-PFIN(4)-PFIN(3)
0143         IF(EPZREM.LT.0..OR.EPZREM.GT.2.*PARJ(31)) MERR=MERR+1
0144         IF(ABS(PFIN(6)-PINI(6)).GT.2.1) MERR=MERR+1
0145       ELSE
0146         DO 150 J=1,4
0147   150   IF(ABS(PFIN(J)-PINI(J)).GT.0001*PINI(4)) MERR=MERR+1
0148         IF(ABS(PFIN(6)-PINI(6)).GT.0.1) MERR=MERR+1
0149       ENDIF
0150       IF(MERR.NE.0) WRITE(MSTU(11),1000) (PINI(J),J=1,4),PINI(6),
0151      &(PFIN(J),J=1,4),PFIN(6)
0152
0153 C...Check that all KF codes are known ones, and that partons/particles
0154 C...satisfy energy-momentum-mass relation. Store particle statistics.
0155       DO 160 I=1,N
0156       IF(K(I,1).GT.20) GOTO 160
0157       IF(LUCOMP(K(I,2)).EQ.0) THEN
0158         WRITE(MSTU(11),1100) I
0159         MERR=MERR+1
0160       ENDIF
0161       PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
0162       IF(ABS(PD).GT.MAX(0.1,0.001*P(I,4)**2).OR.P(I,4).LT.0.) THEN
0163         WRITE(MSTU(11),1200) I
0164         MERR=MERR+1
0165       ENDIF
0166   160 CONTINUE
0167       IF(MTEST.GE.1) CALL LUTABU(21)
0168
0169 C...List all erroneous events and some normal ones.
0170       IF(MERR.NE.0.OR.MSTU(24).NE.0.OR.MSTU(28).NE.0) THEN
0171         CALL LULIST(2)
0172       ELSEIF(MTEST.GE.1.AND.MOD(IEV-5,100).EQ.0) THEN
0173         CALL LULIST(1)
0174       ENDIF
0175
0176 C...Stop execution if too many errors. Endresult of run.
0177       IF(MERR.NE.0) NERR=NERR+1
0178       IF(NERR.GE.10) THEN
0179         WRITE(MSTU(11),1300) IEV
0180         STOP
0181       ENDIF
0182   170 CONTINUE
0183       IF(MTEST.GE.1) CALL LUTABU(22)
0184       WRITE(MSTU(11),1400) NERR
0185
0186 C...Reset commonblock variables changed during run.
0187       MSTJ(2)=3
0188       PARJ(17)=0.
0189       PARJ(22)=1.
0190       PARJ(43)=0.5
0191       PARJ(54)=0.
0192       MSTJ(105)=1
0193       MSTJ(107)=0
0194
0195 C...Format statements for output.
0196  1000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
0197      &'in following event'/' sum of',9X,'px',11X,'py',11X,'pz',11X,
0198      &'E',8X,'charge'/' before',2X,4(1X,F12.5),1X,F8.2/' after',3X,
0199      &4(1X,F12.5),1X,F8.2)
0200  1100 FORMAT(/5X,'Entry no.',I4,' in following event not known code')
0201  1200 FORMAT(/5X,'Entry no.',I4,' in following event has faulty ',
0202      &'kinematics')
0203  1300 FORMAT(/5X,'Ten errors experienced by event ',I3/
0204      &5X,'Something is seriously wrong! Execution stopped now!')
0205  1400 FORMAT(/5X,'Number of erroneous or suspect events in run:',I3/
0206      &5X,'(0 fine, 1 acceptable if a single jet, ',
0207      &'>=2 something is wrong)')
0208
0209       RETURN
0210       END