Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pytest.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pytest.f
1 
2 C*********************************************************************
3 
4 C...PYTEST
5 C...A simple program (disguised as subroutine) to run at installation
6 C...as a check that the program works as intended.
7 
8  SUBROUTINE pytest(MTEST)
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pyk,pychge,pycomp
14 C...Commonblocks.
15  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
16  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
17  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
19  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
20  common/pypars/mstp(200),parp(200),msti(200),pari(200)
21  SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/
22 C...Local arrays.
23  dimension psum(5),pini(6),pfin(6)
24 
25 C...Save defaults for values that are changed.
26  mstj1=mstj(1)
27  mstj3=mstj(3)
28  mstj11=mstj(11)
29  mstj42=mstj(42)
30  mstj43=mstj(43)
31  mstj44=mstj(44)
32  parj17=parj(17)
33  parj22=parj(22)
34  parj43=parj(43)
35  parj54=parj(54)
36  mst101=mstj(101)
37  mst104=mstj(104)
38  mst105=mstj(105)
39  mst107=mstj(107)
40  mst116=mstj(116)
41 
42 C...First part: loop over simple events to be generated.
43  IF(mtest.GE.1) CALL pytabu(20)
44  nerr=0
45  DO 180 iev=1,500
46 
47 C...Reset parameter values. Switch on some nonstandard features.
48  mstj(1)=1
49  mstj(3)=0
50  mstj(11)=1
51  mstj(42)=2
52  mstj(43)=4
53  mstj(44)=2
54  parj(17)=0.1d0
55  parj(22)=1.5d0
56  parj(43)=1d0
57  parj(54)=-0.05d0
58  mstj(101)=5
59  mstj(104)=5
60  mstj(105)=0
61  mstj(107)=1
62  IF(iev.EQ.301.OR.iev.EQ.351.OR.iev.EQ.401) mstj(116)=3
63 
64 C...Ten events each for some single jets configurations.
65  IF(iev.LE.50) THEN
66  ity=(iev+9)/10
67  mstj(3)=-1
68  IF(ity.EQ.3.OR.ity.EQ.4) mstj(11)=2
69  IF(ity.EQ.1) CALL py1ent(1,1,15d0,0d0,0d0)
70  IF(ity.EQ.2) CALL py1ent(1,3101,15d0,0d0,0d0)
71  IF(ity.EQ.3) CALL py1ent(1,-2203,15d0,0d0,0d0)
72  IF(ity.EQ.4) CALL py1ent(1,-4,30d0,0d0,0d0)
73  IF(ity.EQ.5) CALL py1ent(1,21,15d0,0d0,0d0)
74 
75 C...Ten events each for some simple jet systems; string fragmentation.
76  ELSEIF(iev.LE.130) THEN
77  ity=(iev-41)/10
78  IF(ity.EQ.1) CALL py2ent(1,1,-1,40d0)
79  IF(ity.EQ.2) CALL py2ent(1,4,-4,30d0)
80  IF(ity.EQ.3) CALL py2ent(1,2,2103,100d0)
81  IF(ity.EQ.4) CALL py2ent(1,21,21,40d0)
82  IF(ity.EQ.5) CALL py3ent(1,2101,21,-3203,30d0,0.6d0,0.8d0)
83  IF(ity.EQ.6) CALL py3ent(1,5,21,-5,40d0,0.9d0,0.8d0)
84  IF(ity.EQ.7) CALL py3ent(1,21,21,21,60d0,0.7d0,0.5d0)
85  IF(ity.EQ.8) CALL py4ent(1,2,21,21,-2,40d0,
86  & 0.4d0,0.64d0,0.6d0,0.12d0,0.2d0)
87 
88 C...Seventy events with independent fragmentation and momentum cons.
89  ELSEIF(iev.LE.200) THEN
90  ity=1+(iev-131)/16
91  mstj(2)=1+mod(iev-131,4)
92  mstj(3)=1+mod((iev-131)/4,4)
93  IF(ity.EQ.1) CALL py2ent(1,4,-5,40d0)
94  IF(ity.EQ.2) CALL py3ent(1,3,21,-3,40d0,0.9d0,0.4d0)
95  IF(ity.EQ.3) CALL py4ent(1,2,21,21,-2,40d0,
96  & 0.4d0,0.64d0,0.6d0,0.12d0,0.2d0)
97  IF(ity.GE.4) CALL py4ent(1,2,-3,3,-2,40d0,
98  & 0.4d0,0.64d0,0.6d0,0.12d0,0.2d0)
99 
100 C...A hundred events with random jets (check invariant mass).
101  ELSEIF(iev.LE.300) THEN
102  100 DO 110 j=1,5
103  psum(j)=0d0
104  110 CONTINUE
105  njet=2d0+6d0*pyr(0)
106  DO 130 i=1,njet
107  kfl=21
108  IF(i.EQ.1) kfl=int(1d0+4d0*pyr(0))
109  IF(i.EQ.njet) kfl=-int(1d0+4d0*pyr(0))
110  ejet=5d0+20d0*pyr(0)
111  theta=acos(2d0*pyr(0)-1d0)
112  phi=6.2832d0*pyr(0)
113  IF(i.LT.njet) CALL py1ent(-i,kfl,ejet,theta,phi)
114  IF(i.EQ.njet) CALL py1ent(i,kfl,ejet,theta,phi)
115  IF(i.EQ.1.OR.i.EQ.njet) mstj(93)=1
116  IF(i.EQ.1.OR.i.EQ.njet) psum(5)=psum(5)+pymass(kfl)
117  DO 120 j=1,4
118  psum(j)=psum(j)+p(i,j)
119  120 CONTINUE
120  130 CONTINUE
121  IF(psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2.LT.
122  & (psum(5)+parj(32))**2) goto 100
123 
124 C...Fifty e+e- continuum events with matrix elements.
125  ELSEIF(iev.LE.350) THEN
126  mstj(101)=2
127  CALL pyeevt(0,40d0)
128 
129 C...Fifty e+e- continuum event with varying shower options.
130  ELSEIF(iev.LE.400) THEN
131  mstj(42)=1+mod(iev,2)
132  mstj(43)=1+mod(iev/2,4)
133  mstj(44)=mod(iev/8,3)
134  CALL pyeevt(0,90d0)
135 
136 C...Fifty e+e- continuum events with coherent shower.
137  ELSEIF(iev.LE.450) THEN
138  CALL pyeevt(0,500d0)
139 
140 C...Fifty Upsilon decays to ggg or gammagg with coherent shower.
141  ELSE
142  CALL pyonia(5,9.46d0)
143  ENDIF
144 
145 C...Generate event. Find total momentum, energy and charge.
146  DO 140 j=1,4
147  pini(j)=pyp(0,j)
148  140 CONTINUE
149  pini(6)=pyp(0,6)
150  CALL pyexec
151  DO 150 j=1,4
152  pfin(j)=pyp(0,j)
153  150 CONTINUE
154  pfin(6)=pyp(0,6)
155 
156 C...Check conservation of energy, momentum and charge;
157 C...usually exact, but only approximate for single jets.
158  merr=0
159  IF(iev.LE.50) THEN
160  IF((pfin(1)-pini(1))**2+(pfin(2)-pini(2))**2.GE.10d0)
161  & merr=merr+1
162  epzrem=pini(4)+pini(3)-pfin(4)-pfin(3)
163  IF(epzrem.LT.0d0.OR.epzrem.GT.2d0*parj(31)) merr=merr+1
164  IF(abs(pfin(6)-pini(6)).GT.2.1d0) merr=merr+1
165  ELSE
166  DO 160 j=1,4
167  IF(abs(pfin(j)-pini(j)).GT.0.0001d0*pini(4)) merr=merr+1
168  160 CONTINUE
169  IF(abs(pfin(6)-pini(6)).GT.0.1d0) merr=merr+1
170  ENDIF
171  IF(merr.NE.0) WRITE(mstu(11),5000) (pini(j),j=1,4),pini(6),
172  & (pfin(j),j=1,4),pfin(6)
173 
174 C...Check that all KF codes are known ones, and that partons/particles
175 C...satisfy energy-momentum-mass relation. Store particle statistics.
176  DO 170 i=1,n
177  IF(k(i,1).GT.20) goto 170
178  IF(pycomp(k(i,2)).EQ.0) THEN
179  WRITE(mstu(11),5100) i
180  merr=merr+1
181  ENDIF
182  pd=p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2-p(i,5)**2
183  IF(abs(pd).GT.max(0.1d0,0.001d0*p(i,4)**2).OR.p(i,4).LT.0d0)
184  & THEN
185  WRITE(mstu(11),5200) i
186  merr=merr+1
187  ENDIF
188  170 CONTINUE
189  IF(mtest.GE.1) CALL pytabu(21)
190 
191 C...List all erroneous events and some normal ones.
192  IF(merr.NE.0.OR.mstu(24).NE.0.OR.mstu(28).NE.0) THEN
193  IF(merr.GE.1) WRITE(mstu(11),6400)
194  CALL pylist(2)
195  ELSEIF(mtest.GE.1.AND.mod(iev-5,100).EQ.0) THEN
196  CALL pylist(1)
197  ENDIF
198 
199 C...Stop execution if too many errors.
200  IF(merr.NE.0) nerr=nerr+1
201  IF(nerr.GE.10) THEN
202  WRITE(mstu(11),6300)
203  CALL pylist(1)
204  CALL pystop(9)
205  ENDIF
206  180 CONTINUE
207 
208 C...Summarize result of run.
209  IF(mtest.GE.1) CALL pytabu(22)
210 
211 C...Reset commonblock variables changed during run.
212  mstj(1)=mstj1
213  mstj(3)=mstj3
214  mstj(11)=mstj11
215  mstj(42)=mstj42
216  mstj(43)=mstj43
217  mstj(44)=mstj44
218  parj(17)=parj17
219  parj(22)=parj22
220  parj(43)=parj43
221  parj(54)=parj54
222  mstj(101)=mst101
223  mstj(104)=mst104
224  mstj(105)=mst105
225  mstj(107)=mst107
226  mstj(116)=mst116
227 
228 C...Second part: complete events of various kinds.
229 C...Common initial values. Loop over initiating conditions.
230  mstp(122)=max(0,min(2,mtest))
231  mdcy(pycomp(111),1)=0
232  DO 230 iproc=1,8
233 
234 C...Reset process type, kinematics cuts, and the flags used.
235  msel=0
236  DO 190 isub=1,500
237  msub(isub)=0
238  190 CONTINUE
239  ckin(1)=2d0
240  ckin(3)=0d0
241  mstp(2)=1
242  mstp(11)=0
243  mstp(33)=0
244  mstp(81)=1
245  mstp(82)=1
246  mstp(111)=1
247  mstp(131)=0
248  mstp(133)=0
249  parp(131)=0.01d0
250 
251 C...Prompt photon production at fixed target.
252  IF(iproc.EQ.1) THEN
253  pzsum=300d0
254  pesum=sqrt(pzsum**2+pymass(211)**2)+pymass(2212)
255  pqsum=2d0
256  msel=10
257  ckin(3)=5d0
258  CALL pyinit('FIXT','pi+','p',pzsum)
259 
260 C...QCD processes at ISR energies.
261  ELSEIF(iproc.EQ.2) THEN
262  pesum=63d0
263  pzsum=0d0
264  pqsum=2d0
265  msel=1
266  ckin(3)=5d0
267  CALL pyinit('CMS','p','p',pesum)
268 
269 C...W production + multiple interactions at CERN Collider.
270  ELSEIF(iproc.EQ.3) THEN
271  pesum=630d0
272  pzsum=0d0
273  pqsum=0d0
274  msel=12
275  ckin(1)=20d0
276  mstp(82)=4
277  mstp(2)=2
278  mstp(33)=3
279  CALL pyinit('CMS','p','pbar',pesum)
280 
281 C...W/Z gauge boson pairs + pileup events at the Tevatron.
282  ELSEIF(iproc.EQ.4) THEN
283  pesum=1800d0
284  pzsum=0d0
285  pqsum=0d0
286  msub(22)=1
287  msub(23)=1
288  msub(25)=1
289  ckin(1)=200d0
290  mstp(111)=0
291  mstp(131)=1
292  mstp(133)=2
293  parp(131)=0.04d0
294  CALL pyinit('CMS','p','pbar',pesum)
295 
296 C...Higgs production at LHC.
297  ELSEIF(iproc.EQ.5) THEN
298  pesum=15400d0
299  pzsum=0d0
300  pqsum=2d0
301  msub(3)=1
302  msub(102)=1
303  msub(123)=1
304  msub(124)=1
305  pmas(25,1)=300d0
306  ckin(1)=200d0
307  mstp(81)=0
308  mstp(111)=0
309  CALL pyinit('CMS','p','p',pesum)
310 
311 C...Z' production at SSC.
312  ELSEIF(iproc.EQ.6) THEN
313  pesum=40000d0
314  pzsum=0d0
315  pqsum=2d0
316  msel=21
317  pmas(32,1)=600d0
318  ckin(1)=400d0
319  mstp(81)=0
320  mstp(111)=0
321  CALL pyinit('CMS','p','p',pesum)
322 
323 C...W pair production at 1 TeV e+e- collider.
324  ELSEIF(iproc.EQ.7) THEN
325  pesum=1000d0
326  pzsum=0d0
327  pqsum=0d0
328  msub(25)=1
329  msub(69)=1
330  mstp(11)=1
331  CALL pyinit('CMS','e+','e-',pesum)
332 
333 C...Deep inelastic scattering at a LEP+LHC ep collider.
334  ELSEIF(iproc.EQ.8) THEN
335  p(1,1)=0d0
336  p(1,2)=0d0
337  p(1,3)=8000d0
338  p(2,1)=0d0
339  p(2,2)=0d0
340  p(2,3)=-80d0
341  pesum=8080d0
342  pzsum=7920d0
343  pqsum=0d0
344  msub(10)=1
345  ckin(3)=50d0
346  mstp(111)=0
347  CALL pyinit('3MOM','p','e-',pesum)
348  ENDIF
349 
350 C...Generate 20 events of each required type.
351  DO 220 iev=1,20
352  CALL pyevnt
353  pesumm=pesum
354  IF(iproc.EQ.4) pesumm=msti(41)*pesum
355 
356 C...Check conservation of energy/momentum/flavour.
357  pini(1)=0d0
358  pini(2)=0d0
359  pini(3)=pzsum
360  pini(4)=pesumm
361  pini(6)=pqsum
362  DO 200 j=1,4
363  pfin(j)=pyp(0,j)
364  200 CONTINUE
365  pfin(6)=pyp(0,6)
366  merr=0
367  deve=abs(pfin(4)-pini(4))+abs(pfin(3)-pini(3))
368  devt=abs(pfin(1)-pini(1))+abs(pfin(2)-pini(2))
369  devq=abs(pfin(6)-pini(6))
370  IF(deve.GT.2d-3*pesum.OR.devt.GT.max(0.01d0,1d-4*pesum).OR.
371  & devq.GT.0.1d0) merr=1
372  IF(merr.NE.0) WRITE(mstu(11),5000) (pini(j),j=1,4),pini(6),
373  & (pfin(j),j=1,4),pfin(6)
374 
375 C...Check that all KF codes are known ones, and that partons/particles
376 C...satisfy energy-momentum-mass relation.
377  DO 210 i=1,n
378  IF(k(i,1).GT.20) goto 210
379  IF(pycomp(k(i,2)).EQ.0) THEN
380  WRITE(mstu(11),5100) i
381  merr=merr+1
382  ENDIF
383  pd=p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2-p(i,5)**2*
384  & sign(1d0,p(i,5))
385  IF(abs(pd).GT.max(0.1d0,0.002d0*p(i,4)**2,0.002d0*p(i,5)**2)
386  & .OR.(p(i,5).GE.0d0.AND.p(i,4).LT.0d0)) THEN
387  WRITE(mstu(11),5200) i
388  merr=merr+1
389  ENDIF
390  210 CONTINUE
391 
392 C...Listing of erroneous events, and first event of each type.
393  IF(merr.GE.1) nerr=nerr+1
394  IF(nerr.GE.10) THEN
395  WRITE(mstu(11),6300)
396  CALL pylist(1)
397  CALL pystop(9)
398  ENDIF
399  IF(mtest.GE.1.AND.(merr.GE.1.OR.iev.EQ.1)) THEN
400  IF(merr.GE.1) WRITE(mstu(11),6400)
401  CALL pylist(1)
402  ENDIF
403  220 CONTINUE
404 
405 C...List statistics for each process type.
406  IF(mtest.GE.1) CALL pystat(1)
407  230 CONTINUE
408 
409 C...Summarize result of run.
410  IF(nerr.EQ.0) WRITE(mstu(11),6500)
411  IF(nerr.GT.0) WRITE(mstu(11),6600) nerr
412 
413 C...Format statements for output.
414  5000 FORMAT(/' Momentum, energy and/or charge were not conserved ',
415  &'in following event'/' sum of',9x,'px',11x,'py',11x,'pz',11x,
416  &'E',8x,'charge'/' before',2x,4(1x,f12.5),1x,f8.2/' after',3x,
417  &4(1x,f12.5),1x,f8.2)
418  5100 FORMAT(/5x,'Entry no.',i4,' in following event not known code')
419  5200 FORMAT(/5x,'Entry no.',i4,' in following event has faulty ',
420  &'kinematics')
421  6300 FORMAT(/5x,'This is the tenth error experienced! Something is ',
422  &'wrong.'/5x,'Execution will be stopped after listing of event.')
423  6400 FORMAT(5x,'Faulty event follows:')
424  6500 FORMAT(//5x,'End result of PYTEST: no errors detected.')
425  6600 FORMAT(//5x,'End result of PYTEST:',i2,' errors detected.'/
426  &5x,'This should not have happened!')
427 
428  RETURN
429  END