HIJING/hipyset1_35/lutest_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUTEST_HIJING(MTEST)
   6
   7 C...Purpose: to provide a simple program (disguised as subroutine) to
   8 C...run at installation as a check that the program works as intended.
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11       DIMENSION PSUM(5),PINI(6),PFIN(6)
  12
  13 C...Loop over events to be generated.
  14       IF(MTEST.GE.1) CALL LUTABU_HIJING(20)
  15       NERR=0
  16       DO 170 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_HIJING(1,1,15.,0.,0.)
  41         IF(ITY.EQ.2) CALL LU1ENT_HIJING(1,3101,15.,0.,0.)
  42         IF(ITY.EQ.3) CALL LU1ENT_HIJING(1,-2203,15.,0.,0.)
  43         IF(ITY.EQ.4) CALL LU1ENT_HIJING(1,-4,30.,0.,0.)
  44         IF(ITY.EQ.5) CALL LU1ENT_HIJING(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_HIJING(1,1,-1,40.)
  50         IF(ITY.EQ.2) CALL LU2ENT_HIJING(1,4,-4,30.)
  51         IF(ITY.EQ.3) CALL LU2ENT_HIJING(1,2,2103,100.)
  52         IF(ITY.EQ.4) CALL LU2ENT_HIJING(1,21,21,40.)
  53         IF(ITY.EQ.5) CALL LU3ENT_HIJING(1,2101,21,-3203,30.,0.6,0.8)
  54         IF(ITY.EQ.6) CALL LU3ENT_HIJING(1,5,21,-5,40.,0.9,0.8)
  55         IF(ITY.EQ.7) CALL LU3ENT_HIJING(1,21,21,21,60.,0.7,0.5)
  56         IF(ITY.EQ.8) CALL LU4ENT_HIJING(1,2,21,21,-2,40.,0.4,0.64,0.6,0
  57      $       .12,0.2)
  58
  59 C...Seventy events with independent fragmentation and momentum cons.
  60       ELSEIF(IEV.LE.200) THEN
  61         ITY=1+(IEV-131)/16
  62         MSTJ(2)=1+MOD(IEV-131,4)
  63         MSTJ(3)=1+MOD((IEV-131)/4,4)
  64         IF(ITY.EQ.1) CALL LU2ENT_HIJING(1,4,-5,40.)
  65         IF(ITY.EQ.2) CALL LU3ENT_HIJING(1,3,21,-3,40.,0.9,0.4)
  66         IF(ITY.EQ.3) CALL LU4ENT_HIJING(1,2,21,21,-2,40.,0.4,0.64,0.6,0
  67      $       .12,0.2)
  68         IF(ITY.GE.4) CALL LU4ENT_HIJING(1,2,-3,3,-2,40.,0.4,0.64,0.6,0
  69      $       .12,0.2)
  70
  71 C...A hundred events with random jets (check invariant mass).
  72       ELSEIF(IEV.LE.300) THEN
  73   100   DO 110 J=1,5
  74   110   PSUM(J)=0.
  75         NJET=2.+6.*RLU_HIJING(0)
  76         DO 120 I=1,NJET
  77         KFL=21
  78         IF(I.EQ.1) KFL=INT(1.+4.*RLU_HIJING(0))
  79         IF(I.EQ.NJET) KFL=-INT(1.+4.*RLU_HIJING(0))
  80         EJET=5.+20.*RLU_HIJING(0)
  81         THETA=ACOS(2.*RLU_HIJING(0)-1.)
  82         PHI=6.2832*RLU_HIJING(0)
  83         IF(I.LT.NJET) CALL LU1ENT_HIJING(-I,KFL,EJET,THETA,PHI)
  84         IF(I.EQ.NJET) CALL LU1ENT_HIJING(I,KFL,EJET,THETA,PHI)
  85         IF(I.EQ.1.OR.I.EQ.NJET) PSUM(5)=PSUM(5)+ULMASS_HIJING(KFL)
  86         DO 120 J=1,4
  87   120   PSUM(J)=PSUM(J)+P(I,J)
  88         IF(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2.LT.
  89      &  (PSUM(5)+PARJ(32))**2) GOTO 100
  90
  91 C...Fifty e+e- continuum events with matrix elements.
  92       ELSEIF(IEV.LE.350) THEN
  93         MSTJ(101)=2
  94         CALL LUEEVT_HIJING(0,40.)
  95
  96 C...Fifty e+e- continuum event with varying shower options.
  97       ELSEIF(IEV.LE.400) THEN
  98         MSTJ(42)=1+MOD(IEV,2)
  99         MSTJ(43)=1+MOD(IEV/2,4)
 100         MSTJ(44)=MOD(IEV/8,3)
 101         CALL LUEEVT_HIJING(0,90.)
 102
 103 C...Fifty e+e- continuum events with coherent shower, including top.
 104       ELSEIF(IEV.LE.450) THEN
 105         MSTJ(104)=6
 106         CALL LUEEVT_HIJING(0,500.)
 107
 108 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
 109       ELSEIF(IEV.LE.500) THEN
 110         CALL LUONIA_HIJING(5,9.46)
 111
 112 C...One decay each for some heavy mesons.
 113       ELSEIF(IEV.LE.560) THEN
 114         ITY=IEV-501
 115         KFLS=2*(ITY/20)+1
 116         KFLB=8-MOD(ITY/5,4)
 117         KFLC=KFLB-MOD(ITY,5)
 118         CALL LU1ENT_HIJING(1,100*KFLB+10*KFLC+KFLS,0.,0.,0.)
 119
 120 C...One decay each for some heavy baryons.
 121       ELSEIF(IEV.LE.600) THEN
 122         ITY=IEV-561
 123         KFLS=2*(ITY/20)+2
 124         KFLA=8-MOD(ITY/5,4)
 125         KFLB=KFLA-MOD(ITY,5)
 126         KFLC=MAX(1,KFLB-1)
 127         CALL LU1ENT_HIJING(1,1000*KFLA+100*KFLB+10*KFLC+KFLS,0.,0.,0.)
 128       ENDIF
 129
 130 C...Generate event. Find total momentum, energy and charge.
 131       DO 130 J=1,4
 132   130 PINI(J)=PLU_HIJING(0,J)
 133       PINI(6)=PLU_HIJING(0,6)
 134       CALL LUEXEC_HIJING
 135       DO 140 J=1,4
 136   140 PFIN(J)=PLU_HIJING(0,J)
 137       PFIN(6)=PLU_HIJING(0,6)
 138
 139 C...Check conservation of energy, momentum and charge;
 140 C...usually exact, but only approximate for single jets.
 141       MERR=0
 142       IF(IEV.LE.50) THEN
 143         IF((PFIN(1)-PINI(1))**2+(PFIN(2)-PINI(2))**2.GE.4.) MERR=MERR+1
 144         EPZREM=PINI(4)+PINI(3)-PFIN(4)-PFIN(3)
 145         IF(EPZREM.LT.0..OR.EPZREM.GT.2.*PARJ(31)) MERR=MERR+1
 146         IF(ABS(PFIN(6)-PINI(6)).GT.2.1) MERR=MERR+1
 147       ELSE
 148         DO 150 J=1,4
 149   150   IF(ABS(PFIN(J)-PINI(J)).GT.0001*PINI(4)) MERR=MERR+1
 150         IF(ABS(PFIN(6)-PINI(6)).GT.0.1) MERR=MERR+1
 151       ENDIF
 152       IF(MERR.NE.0) WRITE(MSTU(11),1000) (PINI(J),J=1,4),PINI(6),
 153      &(PFIN(J),J=1,4),PFIN(6)
 154
 155 C...Check that all KF codes are known ones, and that partons/particles
 156 C...satisfy energy-momentum-mass relation. Store particle statistics.
 157       DO 160 I=1,N
 158       IF(K(I,1).GT.20) GOTO 160
 159       IF(LUCOMP_HIJING(K(I,2)).EQ.0) THEN
 160         WRITE(MSTU(11),1100) I
 161         MERR=MERR+1
 162       ENDIF
 163       PD=P(I,4)**2-P(I,1)**2-P(I,2)**2-P(I,3)**2-P(I,5)**2
 164       IF(ABS(PD).GT.MAX(0.1,0.001*P(I,4)**2).OR.P(I,4).LT.0.) THEN
 165         WRITE(MSTU(11),1200) I
 166         MERR=MERR+1
 167       ENDIF
 168   160 CONTINUE
 169       IF(MTEST.GE.1) CALL LUTABU_HIJING(21)
 170
 171 C...List all erroneous events and some normal ones.
 172       IF(MERR.NE.0.OR.MSTU(24).NE.0.OR.MSTU(28).NE.0) THEN
 173         CALL LULIST_HIJING(2)
 174       ELSEIF(MTEST.GE.1.AND.MOD(IEV-5,100).EQ.0) THEN
 175         CALL LULIST_HIJING(1)
 176       ENDIF
 177
 178 C...Stop execution if too many errors. Endresult of run.
 179       IF(MERR.NE.0) NERR=NERR+1
 180       IF(NERR.GE.10) THEN
 181         WRITE(MSTU(11),1300) IEV
 182         STOP
 183       ENDIF
 184   170 CONTINUE
 185       IF(MTEST.GE.1) CALL LUTABU_HIJING(22)
 186       WRITE(MSTU(11),1400) NERR
 187
 188 C...Reset commonblock variables changed during run.
 189       MSTJ(2)=3
 190       PARJ(17)=0.
 191       PARJ(22)=1.
 192       PARJ(43)=0.5
 193       PARJ(54)=0.
 194       MSTJ(105)=1
 195       MSTJ(107)=0
 196
 197 C...Format statements for output.
 198  1000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
 199      &'in following event'/' sum of',9X,'px',11X,'py',11X,'pz',11X,
 200      &'E',8X,'charge'/' before',2X,4(1X,F12.5),1X,F8.2/' after',3X,
 201      &4(1X,F12.5),1X,F8.2)
 202  1100 FORMAT(/5X,'Entry no.',I4,' in following event not known code')
 203  1200 FORMAT(/5X,'Entry no.',I4,' in following event has faulty ',
 204      &'kinematics')
 205  1300 FORMAT(/5X,'Ten errors experienced by event ',I3/
 206      &5X,'Something is seriously wrong! Execution stopped now!')
 207  1400 FORMAT(/5X,'Number of erroneous or suspect events in run:',I3/
 208      &5X,'(0 fine, 1 acceptable if a single jet, ',
 209      &'>=2 something is wrong)')
 210
 211       RETURN
 212       END