Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhirand.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhirand.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhirand
5 
6 C...Generates quantities characterizing the high-pT scattering at the
7 C...parton level according to the matrix elements. Chooses incoming,
8 C...reacting partons, their momentum fractions and one of the possible
9 C...subprocesses.
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/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
23  SAVE /pyhiint3/
24  common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3)
25  SAVE /pyhiint4/
26  common/pyhiint5/ngen(0:200,3),xsec(0:200,3)
27  SAVE /pyhiint5/
28 
29 C...Initial values, specifically for (first) semihard interaction.
30  mint(17)=0
31  mint(18)=0
32  vint(143)=1.
33  vint(144)=1.
34  IF(msub(95).EQ.1.OR.mint(82).GE.2) CALL pyhimult(2)
35  isub=0
36  100 mint(51)=0
37 
38 C...Choice of process type - first event of overlay.
39  IF(mint(82).EQ.1.AND.(isub.LE.90.OR.isub.GT.96)) THEN
40  rsub=xsec(0,1)*rlu(0)
41  DO 110 i=1,200
42  IF(msub(i).NE.1) goto 110
43  isub=i
44  rsub=rsub-xsec(i,1)
45  IF(rsub.LE.0.) goto 120
46  110 CONTINUE
47  120 IF(isub.EQ.95) isub=96
48 
49 C...Choice of inclusive process type - overlayed events.
50  ELSEIF(mint(82).GE.2.AND.isub.EQ.0) THEN
51  rsub=vint(131)*rlu(0)
52  isub=96
53  IF(rsub.GT.vint(106)) isub=93
54  IF(rsub.GT.vint(106)+vint(104)) isub=92
55  IF(rsub.GT.vint(106)+vint(104)+vint(103)) isub=91
56  ENDIF
57  IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)+1
58  IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)+1
59  mint(1)=isub
60 
61 C...Find resonances (explicit or implicit in cross-section).
62  mint(72)=0
63  kfr1=0
64  IF(iset(isub).EQ.1.OR.iset(isub).EQ.3) THEN
65  kfr1=kfpr(isub,1)
66  ELSEIF(isub.GE.71.AND.isub.LE.77) THEN
67  kfr1=25
68  ENDIF
69  IF(kfr1.NE.0) THEN
70  taur1=pmas(kfr1,1)**2/vint(2)
71  gamr1=pmas(kfr1,1)*pmas(kfr1,2)/vint(2)
72  mint(72)=1
73  mint(73)=kfr1
74  vint(73)=taur1
75  vint(74)=gamr1
76  ENDIF
77  IF(isub.EQ.141) THEN
78  kfr2=23
79  taur2=pmas(kfr2,1)**2/vint(2)
80  gamr2=pmas(kfr2,1)*pmas(kfr2,2)/vint(2)
81  mint(72)=2
82  mint(74)=kfr2
83  vint(75)=taur2
84  vint(76)=gamr2
85  ENDIF
86 
87 C...Find product masses and minimum pT of process,
88 C...optionally with broadening according to a truncated Breit-Wigner.
89  vint(63)=0.
90  vint(64)=0.
91  mint(71)=0
92  vint(71)=ckin(3)
93  IF(mint(82).GE.2) vint(71)=0.
94  IF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN
95  DO 130 i=1,2
96  IF(kfpr(isub,i).EQ.0) THEN
97  ELSEIF(mstp(42).LE.0) THEN
98  vint(62+i)=pmas(kfpr(isub,i),1)**2
99  ELSE
100  vint(62+i)=ulmass(kfpr(isub,i))**2
101  ENDIF
102  130 CONTINUE
103  IF(min(vint(63),vint(64)).LT.ckin(6)**2) mint(71)=1
104  IF(mint(71).EQ.1) vint(71)=max(ckin(3),ckin(5))
105  ENDIF
106 
107  IF(iset(isub).EQ.0) THEN
108 C...Double or single diffractive, or elastic scattering:
109 C...choose m^2 according to 1/m^2 (diffractive), constant (elastic)
110  is=int(1.5+rlu(0))
111  vint(63)=vint(3)**2
112  vint(64)=vint(4)**2
113  IF(isub.EQ.92.OR.isub.EQ.93) vint(62+is)=parp(111)**2
114  IF(isub.EQ.93) vint(65-is)=parp(111)**2
115  sh=vint(2)
116  sqm1=vint(3)**2
117  sqm2=vint(4)**2
118  sqm3=vint(63)
119  sqm4=vint(64)
120  sqla12=(sh-sqm1-sqm2)**2-4.*sqm1*sqm2
121  sqla34=(sh-sqm3-sqm4)**2-4.*sqm3*sqm4
122  thter1=sqm1+sqm2+sqm3+sqm4-(sqm1-sqm2)*(sqm3-sqm4)/sh-sh
123  thter2=sqrt(max(0.,sqla12))*sqrt(max(0.,sqla34))/sh
124  thl=0.5*(thter1-thter2)
125  thu=0.5*(thter1+thter2)
126  thm=min(max(thl,parp(101)),thu)
127  jtmax=0
128  IF(isub.EQ.92.OR.isub.EQ.93) jtmax=isub-91
129  DO 140 jt=1,jtmax
130  mint(13+3*jt-is*(2*jt-3))=1
131  sqmmin=vint(59+3*jt-is*(2*jt-3))
132  sqmi=vint(8-3*jt+is*(2*jt-3))**2
133  sqmj=vint(3*jt-1-is*(2*jt-3))**2
134  sqmf=vint(68-3*jt+is*(2*jt-3))
135  squa=0.5*sh/sqmi*((1.+(sqmi-sqmj)/sh)*thm+sqmi-sqmf-
136  & sqmj**2/sh+(sqmi+sqmj)*sqmf/sh+(sqmi-sqmj)**2/sh**2*sqmf)
137  quar=sh/sqmi*(thm*(thm+sh-sqmi-sqmj-sqmf*(1.-(sqmi-sqmj)/sh))+
138  & sqmi*sqmj-sqmj*sqmf*(1.+(sqmi-sqmj-sqmf)/sh))
139  sqmmax=squa+sqrt(max(0.,squa**2-quar))
140  IF(abs(quar/squa**2).LT.1.e-06) sqmmax=0.5*quar/squa
141  sqmmax=min(sqmmax,(vint(1)-sqrt(sqmf))**2)
142  vint(59+3*jt-is*(2*jt-3))=sqmmin*(sqmmax/sqmmin)**rlu(0)
143  140 CONTINUE
144 C...Choose t-hat according to exp(B*t-hat+C*t-hat^2).
145  sqm3=vint(63)
146  sqm4=vint(64)
147  sqla34=(sh-sqm3-sqm4)**2-4.*sqm3*sqm4
148  thter1=sqm1+sqm2+sqm3+sqm4-(sqm1-sqm2)*(sqm3-sqm4)/sh-sh
149  thter2=sqrt(max(0.,sqla12))*sqrt(max(0.,sqla34))/sh
150  thl=0.5*(thter1-thter2)
151  thu=0.5*(thter1+thter2)
152  b=vint(121)
153  c=vint(122)
154  IF(isub.EQ.92.OR.isub.EQ.93) THEN
155  b=0.5*b
156  c=0.5*c
157  ENDIF
158  thm=min(max(thl,parp(101)),thu)
159  expth=0.
160  tharg=b*(thm-thu)
161  IF(tharg.GT.-20.) expth=exp(tharg)
162  150 th=thu+log(expth+(1.-expth)*rlu(0))/b
163  th=max(thm,min(thu,th))
164  ratlog=min((b+c*(th+thm))*(th-thm),(b+c*(th+thu))*(th-thu))
165  IF(ratlog.LT.log(rlu(0))) goto 150
166  vint(21)=1.
167  vint(22)=0.
168  vint(23)=min(1.,max(-1.,(2.*th-thter1)/thter2))
169 
170 C...Note: in the following, by In is meant the integral over the
171 C...quantity multiplying coefficient cn.
172 C...Choose tau according to h1(tau)/tau, where
173 C...h1(tau) = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau+tau_R) +
174 C...I0/I3*c3*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
175 C...I0/I4*c4*1/(tau+tau_R') +
176 C...I0/I5*c5*tau/((s*tau-m'^2)^2+(m'*Gamma')^2), and
177 C...c0 + c1 + c2 + c3 + c4 + c5 = 1
178  ELSEIF(iset(isub).GE.1.AND.iset(isub).LE.4) THEN
179  CALL pyhiklim(1)
180  IF(mint(51).NE.0) goto 100
181  rtau=rlu(0)
182  mtau=1
183  IF(rtau.GT.coef(isub,1)) mtau=2
184  IF(rtau.GT.coef(isub,1)+coef(isub,2)) mtau=3
185  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)) mtau=4
186  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4))
187  & mtau=5
188  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4)+
189  & coef(isub,5)) mtau=6
190  CALL pyhikmap(1,mtau,rlu(0))
191 
192 C...2 -> 3, 4 processes:
193 C...Choose tau' according to h4(tau,tau')/tau', where
194 C...h4(tau,tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau', and
195 C...c0 + c1 = 1.
196  IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) THEN
197  CALL pyhiklim(4)
198  IF(mint(51).NE.0) goto 100
199  rtaup=rlu(0)
200  mtaup=1
201  IF(rtaup.GT.coef(isub,15)) mtaup=2
202  CALL pyhikmap(4,mtaup,rlu(0))
203  ENDIF
204 
205 C...Choose y* according to h2(y*), where
206 C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
207 C...I0/I3*c3*1/cosh(y*), I0 = y*max-y*min, and c1 + c2 + c3 = 1.
208  CALL pyhiklim(2)
209  IF(mint(51).NE.0) goto 100
210  ryst=rlu(0)
211  myst=1
212  IF(ryst.GT.coef(isub,7)) myst=2
213  IF(ryst.GT.coef(isub,7)+coef(isub,8)) myst=3
214  CALL pyhikmap(2,myst,rlu(0))
215 
216 C...2 -> 2 processes:
217 C...Choose cos(theta-hat) (cth) according to h3(cth), where
218 C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
219 C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
220 C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
221 C...and c0 + c1 + c2 + c3 + c4 = 1.
222  CALL pyhiklim(3)
223  IF(mint(51).NE.0) goto 100
224  IF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN
225  rcth=rlu(0)
226  mcth=1
227  IF(rcth.GT.coef(isub,10)) mcth=2
228  IF(rcth.GT.coef(isub,10)+coef(isub,11)) mcth=3
229  IF(rcth.GT.coef(isub,10)+coef(isub,11)+coef(isub,12)) mcth=4
230  IF(rcth.GT.coef(isub,10)+coef(isub,11)+coef(isub,12)+
231  & coef(isub,13)) mcth=5
232  CALL pyhikmap(3,mcth,rlu(0))
233  ENDIF
234 
235 C...Low-pT or multiple interactions (first semihard interaction).
236  ELSEIF(iset(isub).EQ.5) THEN
237  CALL pyhimult(3)
238  isub=mint(1)
239  ENDIF
240 
241 C...Choose azimuthal angle.
242  vint(24)=paru(2)*rlu(0)
243 
244 C...Check against user cuts on kinematics at parton level.
245  mint(51)=0
246  IF(isub.LE.90.OR.isub.GT.100) CALL pyhiklim(0)
247  IF(mint(51).NE.0) goto 100
248  IF(mint(82).EQ.1.AND.mstp(141).GE.1) THEN
249  mcut=0
250  IF(msub(91)+msub(92)+msub(93)+msub(94)+msub(95).EQ.0)
251  & CALL pyhikcut(mcut)
252  IF(mcut.NE.0) goto 100
253  ENDIF
254 
255 C...Calculate differential cross-section for different subprocesses.
256  CALL pyhisigh(nchn,sigs)
257 
258 C...Calculations for Monte Carlo estimate of all cross-sections.
259  IF(mint(82).EQ.1.AND.isub.LE.90.OR.isub.GE.96) THEN
260  xsec(isub,2)=xsec(isub,2)+sigs
261  ELSEIF(mint(82).EQ.1) THEN
262  xsec(isub,2)=xsec(isub,2)+xsec(isub,1)
263  ENDIF
264 
265 C...Multiple interactions: store results of cross-section calculation.
266  IF(mint(43).EQ.4.AND.mstp(82).GE.3) THEN
267  vint(153)=sigs
268  CALL pyhimult(4)
269  ENDIF
270 
271 C...Weighting using estimate of maximum of differential cross-section.
272  viol=sigs/xsec(isub,1)
273  IF(viol.LT.rlu(0)) goto 100
274 
275 C...Check for possible violation of estimated maximum of differential
276 C...cross-section used in weighting.
277  IF(mstp(123).LE.0) THEN
278  IF(viol.GT.1.) THEN
279  WRITE(mstu(11),1000) viol,ngen(0,3)+1
280  WRITE(mstu(11),1100) isub,vint(21),vint(22),vint(23),vint(26)
281  stop
282  ENDIF
283  ELSEIF(mstp(123).EQ.1) THEN
284  IF(viol.GT.vint(108)) THEN
285  vint(108)=viol
286 C IF(VIOL.GT.1.) THEN
287 C WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
288 C WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),
289 C & VINT(26)
290 C ENDIF
291  ENDIF
292  ELSEIF(viol.GT.vint(108)) THEN
293  vint(108)=viol
294  IF(viol.GT.1.) THEN
295  xdif=xsec(isub,1)*(viol-1.)
296  xsec(isub,1)=xsec(isub,1)+xdif
297  IF(msub(isub).EQ.1.AND.(isub.LE.90.OR.isub.GT.96))
298  & xsec(0,1)=xsec(0,1)+xdif
299 C WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
300 C WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
301 C IF(ISUB.LE.9) THEN
302 C WRITE(MSTU(11),1300) ISUB,XSEC(ISUB,1)
303 C ELSEIF(ISUB.LE.99) THEN
304 C WRITE(MSTU(11),1400) ISUB,XSEC(ISUB,1)
305 C ELSE
306 C WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1)
307 C ENDIF
308  vint(108)=1.
309  ENDIF
310  ENDIF
311 
312 C...Multiple interactions: choose impact parameter.
313  vint(148)=1.
314  IF(mint(43).EQ.4.AND.(isub.LE.90.OR.isub.GE.96).AND.mstp(82).GE.3)
315  &THEN
316  CALL pyhimult(5)
317  IF(vint(150).LT.rlu(0)) goto 100
318  ENDIF
319  IF(mint(82).EQ.1.AND.msub(95).EQ.1) THEN
320  IF(isub.LE.90.OR.isub.GE.95) ngen(95,1)=ngen(95,1)+1
321  IF(isub.LE.90.OR.isub.GE.96) ngen(96,2)=ngen(96,2)+1
322  ENDIF
323  IF(isub.LE.90.OR.isub.GE.96) mint(31)=mint(31)+1
324 
325 C...Choose flavour of reacting partons (and subprocess).
326  rsigs=sigs*rlu(0)
327  qt2=vint(48)
328  rqqbar=parp(87)*(1.-(qt2/(qt2+(parp(88)*parp(82))**2))**2)
329  IF(isub.NE.95.AND.(isub.NE.96.OR.mstp(82).LE.1.OR.
330  &rlu(0).GT.rqqbar)) THEN
331  DO 190 ichn=1,nchn
332  kfl1=isig(ichn,1)
333  kfl2=isig(ichn,2)
334  mint(2)=isig(ichn,3)
335  rsigs=rsigs-sigh(ichn)
336  IF(rsigs.LE.0.) goto 210
337  190 CONTINUE
338 
339 C...Multiple interactions: choose qqbar preferentially at small pT.
340  ELSEIF(isub.EQ.96) THEN
341  CALL pyhispli(mint(11),21,kfl1,kfldum)
342  CALL pyhispli(mint(12),21,kfl2,kfldum)
343  mint(1)=11
344  mint(2)=1
345  IF(kfl1.EQ.kfl2.AND.rlu(0).LT.0.5) mint(2)=2
346 
347 C...Low-pT: choose string drawing configuration.
348  ELSE
349  kfl1=21
350  kfl2=21
351  rsigs=6.*rlu(0)
352  mint(2)=1
353  IF(rsigs.GT.1.) mint(2)=2
354  IF(rsigs.GT.2.) mint(2)=3
355  ENDIF
356 
357 C...Reassign QCD process. Partons before initial state radiation.
358  210 IF(mint(2).GT.10) THEN
359  mint(1)=mint(2)/10
360  mint(2)=mod(mint(2),10)
361  ENDIF
362  mint(15)=kfl1
363  mint(16)=kfl2
364  mint(13)=mint(15)
365  mint(14)=mint(16)
366  vint(141)=vint(41)
367  vint(142)=vint(42)
368 
369 C...Format statements for differential cross-section maximum violations.
370  1000 FORMAT(1x,'Error: maximum violated by',1p,e11.3,1x,
371  &'in event',1x,i7,'.'/1x,'Execution stopped!')
372  1100 FORMAT(1x,'ISUB = ',i3,'; Point of violation:'/1x,'tau =',1p,
373  &e11.3,', y* =',e11.3,', cthe = ',0p,f11.7,', tau'' =',1p,e11.3)
374  1200 FORMAT(1x,'Warning: maximum violated by',1p,e11.3,1x,
375  &'in event',1x,i7)
376  1300 FORMAT(1x,'XSEC(',i1,',1) increased to',1p,e11.3)
377  1400 FORMAT(1x,'XSEC(',i2,',1) increased to',1p,e11.3)
378  1500 FORMAT(1x,'XSEC(',i3,',1) increased to',1p,e11.3)
379 
380  RETURN
381  END