Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhithia.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhithia.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhithia
5 
6 C...Administers the generation of a high-pt event via calls to a number
7 C...of subroutines; also computes cross-sections.
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/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
15  SAVE /pyhisubs/
16  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
17  SAVE /pyhipars/
18  common/pyhiint1/mint(400),vint(400)
19  SAVE /pyhiint1/
20  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
21  SAVE /pyhiint2/
22  common/pyhiint5/ngen(0:200,3),xsec(0:200,3)
23  SAVE /pyhiint5/
24 
25 C...Loop over desired number of overlayed events (normally 1).
26  mint(7)=0
27  mint(8)=0
28  novl=1
29  IF(mstp(131).NE.0) CALL pyhiovly(2)
30  IF(mstp(131).NE.0) novl=mint(81)
31  mint(83)=0
32  mint(84)=mstp(126)
33  mstu(70)=0
34  DO 190 iovl=1,novl
35  IF(mint(84)+100.GE.mstu(4)) THEN
36  CALL luerrm(11,
37  & '(PYHITHIA:) no more space in LUJETS for overlayed events')
38  IF(mstu(21).GE.1) goto 200
39  ENDIF
40  mint(82)=iovl
41 
42 C...Generate variables of hard scattering.
43  100 CONTINUE
44  IF(iovl.EQ.1) ngen(0,2)=ngen(0,2)+1
45  mint(31)=0
46  mint(51)=0
47  CALL pyhirand
48  isub=mint(1)
49  IF(iovl.EQ.1) THEN
50  ngen(isub,2)=ngen(isub,2)+1
51 
52 C...Store information on hard interaction.
53  DO 110 j=1,200
54  msti(j)=0
55  110 pari(j)=0.
56  msti(1)=mint(1)
57  msti(2)=mint(2)
58  msti(11)=mint(11)
59  msti(12)=mint(12)
60  msti(15)=mint(15)
61  msti(16)=mint(16)
62  msti(17)=mint(17)
63  msti(18)=mint(18)
64  pari(11)=vint(1)
65  pari(12)=vint(2)
66  IF(isub.NE.95) THEN
67  DO 120 j=13,22
68  120 pari(j)=vint(30+j)
69  pari(33)=vint(41)
70  pari(34)=vint(42)
71  pari(35)=pari(33)-pari(34)
72  pari(36)=vint(21)
73  pari(37)=vint(22)
74  pari(38)=vint(26)
75  pari(41)=vint(23)
76  ENDIF
77  ENDIF
78 
79  IF(mstp(111).EQ.-1) goto 160
80  IF(isub.LE.90.OR.isub.GE.95) THEN
81 C...Hard scattering (including low-pT):
82 C...reconstruct kinematics and colour flow of hard scattering.
83  CALL pyhiscat
84  IF(mint(51).EQ.1) goto 100
85 
86 C...Showering of initial state partons (optional).
87  ipu1=mint(84)+1
88  ipu2=mint(84)+2
89  IF(mstp(61).GE.1.AND.mint(43).NE.1.AND.isub.NE.95)
90  & CALL pyhisspa(ipu1,ipu2)
91  nsav1=n
92 
93 C...Multiple interactions.
94  IF(mstp(81).GE.1.AND.mint(43).EQ.4.AND.isub.NE.95)
95  & CALL pyhimult(6)
96  mint(1)=isub
97  nsav2=n
98 
99 C...Hadron remnants and primordial kT.
100  CALL pyhiremn(ipu1,ipu2)
101  IF(mint(51).EQ.1) goto 100
102  nsav3=n
103 
104 C...Showering of final state partons (optional).
105  ipu3=mint(84)+3
106  ipu4=mint(84)+4
107  IF(mstp(71).GE.1.AND.isub.NE.95.AND.k(ipu3,1).GT.0.AND.
108  & k(ipu3,1).LE.10.AND.k(ipu4,1).GT.0.AND.k(ipu4,1).LE.10) THEN
109  qmax=sqrt(parp(71)*vint(52))
110  IF(isub.EQ.5) qmax=sqrt(pmas(23,1)**2)
111  IF(isub.EQ.8) qmax=sqrt(pmas(24,1)**2)
112  CALL lushow(ipu3,ipu4,qmax)
113  ENDIF
114 
115 C...Sum up transverse and longitudinal momenta.
116  IF(iovl.EQ.1) THEN
117  pari(65)=2.*pari(17)
118  DO 130 i=mstp(126)+1,n
119  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 130
120  pt=sqrt(p(i,1)**2+p(i,2)**2)
121  pari(69)=pari(69)+pt
122  IF(i.LE.nsav1.OR.i.GT.nsav3) pari(66)=pari(66)+pt
123  IF(i.GT.nsav1.AND.i.LE.nsav2) pari(68)=pari(68)+pt
124  130 CONTINUE
125  pari(67)=pari(68)
126  pari(71)=vint(151)
127  pari(72)=vint(152)
128  pari(73)=vint(151)
129  pari(74)=vint(152)
130  ENDIF
131 
132 C...Decay of final state resonances.
133  IF(mstp(41).GE.1.AND.isub.NE.95) CALL pyhiresd
134 
135  ELSE
136 C...Diffractive and elastic scattering.
137  CALL pyhidiff
138  IF(iovl.EQ.1) THEN
139  pari(65)=2.*pari(17)
140  pari(66)=pari(65)
141  pari(69)=pari(65)
142  ENDIF
143  ENDIF
144 
145 C...Recalculate energies from momenta and masses (if desired).
146  IF(mstp(113).GE.1) THEN
147  DO 140 i=mint(83)+1,n
148  140 IF(k(i,1).GT.0.AND.k(i,1).LE.10) p(i,4)=sqrt(p(i,1)**2+
149  & p(i,2)**2+p(i,3)**2+p(i,5)**2)
150  ENDIF
151 
152 C...Rearrange partons along strings, check invariant mass cuts.
153  mstu(28)=0
154  CALL luprep(mint(84)+1)
155  IF(mstp(112).EQ.1.AND.mstu(28).EQ.3) goto 100
156  IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) THEN
157  DO 150 i=mint(84)+1,n
158  IF(k(i,2).NE.94) goto 150
159  k(i+1,3)=mod(k(i+1,4)/mstu(5),mstu(5))
160  k(i+2,3)=mod(k(i+2,4)/mstu(5),mstu(5))
161  150 CONTINUE
162  CALL luedit(12)
163  CALL luedit(14)
164  IF(mstp(125).EQ.0) CALL luedit(15)
165  IF(mstp(125).EQ.0) mint(4)=0
166  ENDIF
167 
168 C...Introduce separators between sections in LULIST event listing.
169  IF(iovl.EQ.1.AND.mstp(125).LE.0) THEN
170  mstu(70)=1
171  mstu(71)=n
172  ELSEIF(iovl.EQ.1) THEN
173  mstu(70)=3
174  mstu(71)=2
175  mstu(72)=mint(4)
176  mstu(73)=n
177  ENDIF
178 
179 C...Perform hadronization (if desired).
180  IF(mstp(111).GE.1) CALL luexec
181  IF(mstp(125).EQ.0.OR.mstp(125).EQ.1) CALL luedit(14)
182 
183 C...Calculate Monte Carlo estimates of cross-sections.
184  160 IF(iovl.EQ.1) THEN
185  IF(mstp(111).NE.-1) ngen(isub,3)=ngen(isub,3)+1
186  ngen(0,3)=ngen(0,3)+1
187  xsec(0,3)=0.
188  DO 170 i=1,200
189  IF(i.EQ.96) THEN
190  xsec(i,3)=0.
191  ELSEIF(msub(95).EQ.1.AND.(i.EQ.11.OR.i.EQ.12.OR.i.EQ.13.OR.
192  & i.EQ.28.OR.i.EQ.53.OR.i.EQ.68)) THEN
193  xsec(i,3)=xsec(96,2)*ngen(i,3)/max(1.,float(ngen(96,1))*
194  & float(ngen(96,2)))
195  ELSEIF(ngen(i,1).EQ.0) THEN
196  xsec(i,3)=0.
197  ELSEIF(ngen(i,2).EQ.0) THEN
198  xsec(i,3)=xsec(i,2)*ngen(0,3)/(float(ngen(i,1))*
199  & float(ngen(0,2)))
200  ELSE
201  xsec(i,3)=xsec(i,2)*ngen(i,3)/(float(ngen(i,1))*
202  & float(ngen(i,2)))
203  ENDIF
204  170 xsec(0,3)=xsec(0,3)+xsec(i,3)
205  IF(msub(95).EQ.1) THEN
206  ngens=ngen(91,3)+ngen(92,3)+ngen(93,3)+ngen(94,3)+ngen(95,3)
207  xsecs=xsec(91,3)+xsec(92,3)+xsec(93,3)+xsec(94,3)+xsec(95,3)
208  xmaxs=xsec(95,1)
209  IF(msub(91).EQ.1) xmaxs=xmaxs+xsec(91,1)
210  IF(msub(92).EQ.1) xmaxs=xmaxs+xsec(92,1)
211  IF(msub(93).EQ.1) xmaxs=xmaxs+xsec(93,1)
212  IF(msub(94).EQ.1) xmaxs=xmaxs+xsec(94,1)
213  fac=1.
214  IF(ngens.LT.ngen(0,3)) fac=(xmaxs-xsecs)/(xsec(0,3)-xsecs)
215  xsec(11,3)=fac*xsec(11,3)
216  xsec(12,3)=fac*xsec(12,3)
217  xsec(13,3)=fac*xsec(13,3)
218  xsec(28,3)=fac*xsec(28,3)
219  xsec(53,3)=fac*xsec(53,3)
220  xsec(68,3)=fac*xsec(68,3)
221  xsec(0,3)=xsec(91,3)+xsec(92,3)+xsec(93,3)+xsec(94,3)+
222  & xsec(95,1)
223  ENDIF
224 
225 C...Store final information.
226  mint(5)=mint(5)+1
227  msti(3)=mint(3)
228  msti(4)=mint(4)
229  msti(5)=mint(5)
230  msti(6)=mint(6)
231  msti(7)=mint(7)
232  msti(8)=mint(8)
233  msti(13)=mint(13)
234  msti(14)=mint(14)
235  msti(21)=mint(21)
236  msti(22)=mint(22)
237  msti(23)=mint(23)
238  msti(24)=mint(24)
239  msti(25)=mint(25)
240  msti(26)=mint(26)
241  msti(31)=mint(31)
242  pari(1)=xsec(0,3)
243  pari(2)=xsec(0,3)/mint(5)
244  pari(31)=vint(141)
245  pari(32)=vint(142)
246  IF(isub.NE.95.AND.mint(7)*mint(8).NE.0) THEN
247  pari(42)=2.*vint(47)/vint(1)
248  DO 180 is=7,8
249  pari(36+is)=p(mint(is),3)/vint(1)
250  pari(38+is)=p(mint(is),4)/vint(1)
251  i=mint(is)
252  pr=max(1e-20,p(i,5)**2+p(i,1)**2+p(i,2)**2)
253  pari(40+is)=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/
254  & sqrt(pr),1e20)),p(i,3))
255  pr=max(1e-20,p(i,1)**2+p(i,2)**2)
256  pari(42+is)=sign(log(min((sqrt(pr+p(i,3)**2)+abs(p(i,3)))/
257  & sqrt(pr),1e20)),p(i,3))
258  pari(44+is)=p(i,3)/sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
259  pari(46+is)=ulangl(p(i,3),sqrt(p(i,1)**2+p(i,2)**2))
260  pari(48+is)=ulangl(p(i,1),p(i,2))
261  180 CONTINUE
262  ENDIF
263  pari(61)=vint(148)
264  IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN
265  mstu(161)=mint(21)
266  mstu(162)=0
267  ELSE
268  mstu(161)=mint(21)
269  mstu(162)=mint(22)
270  ENDIF
271  ENDIF
272 
273 C...Prepare to go to next overlayed event.
274  msti(41)=iovl
275  IF(iovl.GE.2.AND.iovl.LE.10) msti(40+iovl)=isub
276  IF(mstu(70).LT.10) THEN
277  mstu(70)=mstu(70)+1
278  mstu(70+mstu(70))=n
279  ENDIF
280  mint(83)=n
281  mint(84)=n+mstp(126)
282  190 CONTINUE
283 
284 C...Information on overlayed events.
285  IF(mstp(131).EQ.1.AND.mstp(133).GE.1) THEN
286  pari(91)=vint(132)
287  pari(92)=vint(133)
288  pari(93)=vint(134)
289  IF(mstp(133).EQ.2) pari(93)=pari(93)*xsec(0,3)/vint(131)
290  ENDIF
291 
292 C...Transform to the desired coordinate frame.
293  200 CALL pyhifram(mstp(124))
294 
295  RETURN
296  END