Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lutest.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  dimension psum(5),pini(6),pfin(6)
13 
14 C...Loop over events to be generated.
15  IF(mtest.GE.1) CALL lutabu(20)
16  nerr=0
17  DO 170 iev=1,600
18 
19 C...Reset parameter values. Switch on some nonstandard features.
20  mstj(1)=1
21  mstj(3)=0
22  mstj(11)=1
23  mstj(42)=2
24  mstj(43)=4
25  mstj(44)=2
26  parj(17)=0.1
27  parj(22)=1.5
28  parj(43)=1.
29  parj(54)=-0.05
30  mstj(101)=5
31  mstj(104)=5
32  mstj(105)=0
33  mstj(107)=1
34  IF(iev.EQ.301.OR.iev.EQ.351.OR.iev.EQ.401) mstj(116)=3
35 
36 C...Ten events each for some single jets configurations.
37  IF(iev.LE.50) THEN
38  ity=(iev+9)/10
39  mstj(3)=-1
40  IF(ity.EQ.3.OR.ity.EQ.4) mstj(11)=2
41  IF(ity.EQ.1) CALL lu1ent(1,1,15.,0.,0.)
42  IF(ity.EQ.2) CALL lu1ent(1,3101,15.,0.,0.)
43  IF(ity.EQ.3) CALL lu1ent(1,-2203,15.,0.,0.)
44  IF(ity.EQ.4) CALL lu1ent(1,-4,30.,0.,0.)
45  IF(ity.EQ.5) CALL lu1ent(1,21,15.,0.,0.)
46 
47 C...Ten events each for some simple jet systems; string fragmentation.
48  ELSEIF(iev.LE.130) THEN
49  ity=(iev-41)/10
50  IF(ity.EQ.1) CALL lu2ent(1,1,-1,40.)
51  IF(ity.EQ.2) CALL lu2ent(1,4,-4,30.)
52  IF(ity.EQ.3) CALL lu2ent(1,2,2103,100.)
53  IF(ity.EQ.4) CALL lu2ent(1,21,21,40.)
54  IF(ity.EQ.5) CALL lu3ent(1,2101,21,-3203,30.,0.6,0.8)
55  IF(ity.EQ.6) CALL lu3ent(1,5,21,-5,40.,0.9,0.8)
56  IF(ity.EQ.7) CALL lu3ent(1,21,21,21,60.,0.7,0.5)
57  IF(ity.EQ.8) CALL lu4ent(1,2,21,21,-2,40.,0.4,0.64,0.6,0.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(1,4,-5,40.)
65  IF(ity.EQ.2) CALL lu3ent(1,3,21,-3,40.,0.9,0.4)
66  IF(ity.EQ.3) CALL lu4ent(1,2,21,21,-2,40.,0.4,0.64,0.6,0.12,0.2)
67  IF(ity.GE.4) CALL lu4ent(1,2,-3,3,-2,40.,0.4,0.64,0.6,0.12,0.2)
68 
69 C...A hundred events with random jets (check invariant mass).
70  ELSEIF(iev.LE.300) THEN
71  100 DO 110 j=1,5
72  110 psum(j)=0.
73  njet=2.+6.*rlu(0)
74  DO 120 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) psum(5)=psum(5)+ulmass(kfl)
84  DO 120 j=1,4
85  120 psum(j)=psum(j)+p(i,j)
86  IF(psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2.LT.
87  & (psum(5)+parj(32))**2) goto 100
88 
89 C...Fifty e+e- continuum events with matrix elements.
90  ELSEIF(iev.LE.350) THEN
91  mstj(101)=2
92  CALL lueevt(0,40.)
93 
94 C...Fifty e+e- continuum event with varying shower options.
95  ELSEIF(iev.LE.400) THEN
96  mstj(42)=1+mod(iev,2)
97  mstj(43)=1+mod(iev/2,4)
98  mstj(44)=mod(iev/8,3)
99  CALL lueevt(0,90.)
100 
101 C...Fifty e+e- continuum events with coherent shower, including top.
102  ELSEIF(iev.LE.450) THEN
103  mstj(104)=6
104  CALL lueevt(0,500.)
105 
106 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
107  ELSEIF(iev.LE.500) THEN
108  CALL luonia(5,9.46)
109 
110 C...One decay each for some heavy mesons.
111  ELSEIF(iev.LE.560) THEN
112  ity=iev-501
113  kfls=2*(ity/20)+1
114  kflb=8-mod(ity/5,4)
115  kflc=kflb-mod(ity,5)
116  CALL lu1ent(1,100*kflb+10*kflc+kfls,0.,0.,0.)
117 
118 C...One decay each for some heavy baryons.
119  ELSEIF(iev.LE.600) THEN
120  ity=iev-561
121  kfls=2*(ity/20)+2
122  kfla=8-mod(ity/5,4)
123  kflb=kfla-mod(ity,5)
124  kflc=max(1,kflb-1)
125  CALL lu1ent(1,1000*kfla+100*kflb+10*kflc+kfls,0.,0.,0.)
126  ENDIF
127 
128 C...Generate event. Find total momentum, energy and charge.
129  DO 130 j=1,4
130  130 pini(j)=plu(0,j)
131  pini(6)=plu(0,6)
132  CALL luexec
133  DO 140 j=1,4
134  140 pfin(j)=plu(0,j)
135  pfin(6)=plu(0,6)
136 
137 C...Check conservation of energy, momentum and charge;
138 C...usually exact, but only approximate for single jets.
139  merr=0
140  IF(iev.LE.50) THEN
141  IF((pfin(1)-pini(1))**2+(pfin(2)-pini(2))**2.GE.4.) merr=merr+1
142  epzrem=pini(4)+pini(3)-pfin(4)-pfin(3)
143  IF(epzrem.LT.0..OR.epzrem.GT.2.*parj(31)) merr=merr+1
144  IF(abs(pfin(6)-pini(6)).GT.2.1) merr=merr+1
145  ELSE
146  DO 150 j=1,4
147  150 IF(abs(pfin(j)-pini(j)).GT.0001*pini(4)) merr=merr+1
148  IF(abs(pfin(6)-pini(6)).GT.0.1) merr=merr+1
149  ENDIF
150  IF(merr.NE.0) WRITE(mstu(11),1000) (pini(j),j=1,4),pini(6),
151  &(pfin(j),j=1,4),pfin(6)
152 
153 C...Check that all KF codes are known ones, and that partons/particles
154 C...satisfy energy-momentum-mass relation. Store particle statistics.
155  DO 160 i=1,n
156  IF(k(i,1).GT.20) goto 160
157  IF(lucomp(k(i,2)).EQ.0) THEN
158  WRITE(mstu(11),1100) i
159  merr=merr+1
160  ENDIF
161  pd=p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2-p(i,5)**2
162  IF(abs(pd).GT.max(0.1,0.001*p(i,4)**2).OR.p(i,4).LT.0.) THEN
163  WRITE(mstu(11),1200) i
164  merr=merr+1
165  ENDIF
166  160 CONTINUE
167  IF(mtest.GE.1) CALL lutabu(21)
168 
169 C...List all erroneous events and some normal ones.
170  IF(merr.NE.0.OR.mstu(24).NE.0.OR.mstu(28).NE.0) THEN
171  CALL lulist(2)
172  ELSEIF(mtest.GE.1.AND.mod(iev-5,100).EQ.0) THEN
173  CALL lulist(1)
174  ENDIF
175 
176 C...Stop execution if too many errors. Endresult of run.
177  IF(merr.NE.0) nerr=nerr+1
178  IF(nerr.GE.10) THEN
179  WRITE(mstu(11),1300) iev
180  stop
181  ENDIF
182  170 CONTINUE
183  IF(mtest.GE.1) CALL lutabu(22)
184  WRITE(mstu(11),1400) nerr
185 
186 C...Reset commonblock variables changed during run.
187  mstj(2)=3
188  parj(17)=0.
189  parj(22)=1.
190  parj(43)=0.5
191  parj(54)=0.
192  mstj(105)=1
193  mstj(107)=0
194 
195 C...Format statements for output.
196  1000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
197  &'in following event'/' sum of',9x,'px',11x,'py',11x,'pz',11x,
198  &'E',8x,'charge'/' before',2x,4(1x,f12.5),1x,f8.2/' after',3x,
199  &4(1x,f12.5),1x,f8.2)
200  1100 FORMAT(/5x,'Entry no.',i4,' in following event not known code')
201  1200 FORMAT(/5x,'Entry no.',i4,' in following event has faulty ',
202  &'kinematics')
203  1300 FORMAT(/5x,'Ten errors experienced by event ',i3/
204  &5x,'Something is seriously wrong! Execution stopped now!')
205  1400 FORMAT(/5x,'Number of erroneous or suspect events in run:',i3/
206  &5x,'(0 fine, 1 acceptable if a single jet, ',
207  &'>=2 something is wrong)')
208 
209  RETURN
210  END