Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhitest.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhitest.f
1 
2 ***********************************************************************
3 
4  SUBROUTINE pyhitest(MTEST)
5 
6 C...Purpose: to provide a simple program (disguised as a 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  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
15  SAVE /ludat3/
16  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
17  SAVE /pyhisubs/
18  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
19  SAVE /pyhipars/
20 
21 C...Common initial values. Loop over initiating conditions.
22  mstp(122)=1
23  IF(mtest.LE.0) mstp(122)=0
24  mdcy(lucomp(111),1)=0
25  nerr=0
26  DO 130 iproc=1,7
27 
28 C...Reset process type, kinematics cuts, and the flags used.
29  msel=0
30  DO 100 isub=1,200
31  100 msub(isub)=0
32  ckin(1)=2.
33  ckin(3)=0.
34  mstp(2)=1
35  mstp(33)=0
36  mstp(81)=1
37  mstp(82)=1
38  mstp(111)=1
39  mstp(131)=0
40  mstp(133)=0
41  parp(131)=0.01
42 
43 C...Prompt photon production at fixed target.
44  IF(iproc.EQ.1) THEN
45  pzsum=300.
46  pesum=sqrt(pzsum**2+ulmass(211)**2)+ulmass(2212)
47  pqsum=2.
48  msel=10
49  ckin(3)=5.
50  CALL pyhiinit('FIXT','pi+','p',pzsum)
51 
52 C...QCD processes at ISR energies.
53  ELSEIF(iproc.EQ.2) THEN
54  pesum=63.
55  pzsum=0.
56  pqsum=2.
57  msel=1
58  ckin(3)=5.
59  CALL pyhiinit('CMS','p','p',pesum)
60 
61 C...W production + multiple interactions at CERN Collider.
62  ELSEIF(iproc.EQ.3) THEN
63  pesum=630.
64  pzsum=0.
65  pqsum=0.
66  msel=12
67  ckin(1)=20.
68  mstp(82)=4
69  mstp(2)=2
70  mstp(33)=3
71  CALL pyhiinit('CMS','p','pbar',pesum)
72 
73 C...W/Z gauge boson pairs + overlayed events at the Tevatron.
74  ELSEIF(iproc.EQ.4) THEN
75  pesum=1800.
76  pzsum=0.
77  pqsum=0.
78  msub(22)=1
79  msub(23)=1
80  msub(25)=1
81  ckin(1)=200.
82  mstp(111)=0
83  mstp(131)=1
84  mstp(133)=2
85  parp(131)=0.04
86  CALL pyhiinit('CMS','p','pbar',pesum)
87 
88 C...Higgs production at LHC.
89  ELSEIF(iproc.EQ.5) THEN
90  pesum=17000.
91  pzsum=0.
92  pqsum=0.
93  msel=16
94  pmas(25,1)=300.
95  ckin(1)=200.
96  mstp(81)=0
97  mstp(111)=0
98  CALL pyhiinit('CMS','p','pbar',pesum)
99 
100 C...Z' production at SSC.
101  ELSEIF(iproc.EQ.6) THEN
102  pesum=40000.
103  pzsum=0.
104  pqsum=0.
105  msel=21
106  pmas(32,1)=600.
107  ckin(1)=400.
108  mstp(81)=0
109  mstp(111)=0
110  CALL pyhiinit('CMS','p','pbar',pesum)
111 
112 C...W pair production at 1 TeV e+e- collider.
113  ELSEIF(iproc.EQ.7) THEN
114  pesum=1000.
115  pzsum=0.
116  pqsum=0.
117  msub(25)=1
118  CALL pyhiinit('CMS','e+','e-',pesum)
119  ENDIF
120 
121 C...Generate 20 events of each required type.
122  DO 120 iev=1,20
123  CALL pyhithia
124  pesumm=pesum
125  IF(iproc.EQ.4) pesumm=msti(41)*pesum
126 
127 C...Check conservation of energy/momentum/flavour.
128  merr=0
129  deve=abs(plu(0,4)-pesumm)+abs(plu(0,3)-pzsum)
130  devt=abs(plu(0,1))+abs(plu(0,2))
131  devq=abs(plu(0,6)-pqsum)
132  IF(deve.GT.1e-3*pesum.OR.devt.GT.max(0.01,1e-5*pesum).OR.
133  &devq.GT.0.1) merr=1
134  IF(merr.NE.0) WRITE(mstu(11),1000) iproc,iev
135 
136 C...Check that all KF codes are known ones, and that partons/particles
137 C...satisfy energy-momentum-mass relation.
138  DO 110 i=1,n
139  IF(k(i,1).GT.20) goto 110
140  IF(lucomp(k(i,2)).EQ.0) THEN
141  WRITE(mstu(11),1100) i
142  merr=merr+1
143  ENDIF
144  pd=p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2-p(i,5)**2*
145  &sign(1.,p(i,5))
146  IF(abs(pd).GT.max(0.1,0.002*p(i,4)**2,0.002*p(i,5)**2).OR.
147  &(p(i,5).GE.0..AND.p(i,4).LT.0.)) THEN
148  WRITE(mstu(11),1200) i
149  merr=merr+1
150  ENDIF
151  110 CONTINUE
152 
153 C...Listing of erronoeus events, and first event of each type.
154  IF(merr.GE.1) nerr=nerr+1
155  IF(nerr.GE.10) THEN
156  WRITE(mstu(11),1300)
157  CALL lulist(1)
158  stop
159  ENDIF
160  IF(mtest.GE.1.AND.(merr.GE.1.OR.iev.EQ.1)) THEN
161  IF(merr.GE.1) WRITE(mstu(11),1400)
162  CALL lulist(1)
163  ENDIF
164  120 CONTINUE
165 
166 C...List statistics for each process type.
167  IF(mtest.GE.1) CALL pyhistat(1)
168  130 CONTINUE
169 
170 C...Summarize result of run.
171  IF(nerr.EQ.0) WRITE(mstu(11),1500)
172  IF(nerr.GT.0) WRITE(mstu(11),1600) nerr
173  RETURN
174 
175 C...Formats for information.
176  1000 FORMAT(/5x,'Energy/momentum/flavour nonconservation for process',
177  &i2,', event',i4)
178  1100 FORMAT(/5x,'Entry no.',i4,' in following event not known code')
179  1200 FORMAT(/5x,'Entry no.',i4,' in following event has faulty ',
180  &'kinematics')
181  1300 FORMAT(/5x,'This is the tenth error experienced! Something is ',
182  &'wrong.'/5x,'Execution will be stopped after listing of event.')
183  1400 FORMAT(5x,'Faulty event follows:')
184  1500 FORMAT(//5x,'End result of run: no errors detected.')
185  1600 FORMAT(//5x,'End result of run:',i2,' errors detected.'/
186  &5x,'This should not have happened!')
187  END