Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pysigh.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pysigh.f
1 
2 C***********************************************************************
3 
4 C...PYSIGH
5 C...Differential matrix elements for all included subprocesses
6 C...Note that what is coded is (disregarding the COMFAC factor)
7 C...1) for 2 -> 1 processes: s-hat/pi*d(sigma-hat), where,
8 C...when d(sigma-hat) is given in the zero-width limit, the delta
9 C...function in tau is replaced by a (modified) Breit-Wigner:
10 C...1/pi*s*H_res/((s*tau-m_res^2)^2+H_res^2),
11 C...where H_res = s-hat/m_res*Gamma_res(s-hat);
12 C...2) for 2 -> 2 processes: (s-hat)**2/pi*d(sigma-hat)/d(t-hat);
13 C...i.e., dimensionless quantities
14 C...3) for 2 -> 3 processes: abs(M)^2, where the total cross-section is
15 C...Integral abs(M)^2/(2shat') * (prod_(i=1)^3 d^3p_i/((2pi)^3*2E_i)) *
16 C...(2pi)^4 delta^4(P - sum p_i)
17 C...COMFAC contains the factor pi/s (or equivalent) and
18 C...the conversion factor from GeV^-2 to mb
19 
20  SUBROUTINE pysigh(NCHN,SIGS)
21 
22 C...Double precision and integer declarations
23  IMPLICIT DOUBLE PRECISION(a-h, o-z)
24  IMPLICIT INTEGER(i-n)
25  INTEGER pyk,pychge,pycomp
26 C...Parameter statement to help give large particle numbers.
27  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
28  &kexcit=4000000,kdimen=5000000)
29 C...Commonblocks
30  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
31  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
32  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
33  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
34  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
35  common/pypars/mstp(200),parp(200),msti(200),pari(200)
36  common/pyint1/mint(400),vint(400)
37  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
38  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
39  common/pyint4/mwid(500),wids(500,5)
40  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
41  common/pyint7/sigt(0:6,0:6,0:5)
42  common/pymssm/imss(0:99),rmss(0:99)
43  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
44  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
45  common/pytcsm/itcm(0:99),rtcm(0:99)
46  common/pysgcm/isub,isubsv,mmin1,mmax1,mmin2,mmax2,mmina,mmaxa,
47  &kfac(2,-40:40),comfac,fack,faca,sh,th,uh,sh2,th2,uh2,sqm3,sqm4,
48  &shr,sqpth,taup,be34,cth,x(2),sqmz,sqmw,gmmz,gmmw,
49  &aem,as,xw,xw1,xwc,xwv,poll,polr,polll,polrr
50  SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,
51  &/pyint1/,/pyint2/,/pyint3/,/pyint4/,/pyint5/,/pyint7/,
52  &/pymssm/,/pyssmt/,/pytcsm/,/pysgcm/
53  include "mc_set.inc"
54 C...Local arrays and complex variables
55  dimension xpq(-25:25)
56 
57 C...Map of processes onto which routine to call
58 C...in order to evaluate cross section:
59 C...0 = not implemented;
60 C...1 = standard QCD (including photons);
61 C...2 = heavy flavours;
62 C...3 = W/Z;
63 C...4 = Higgs (2 doublets; including longitudinal W/Z scattering);
64 C...5 = SUSY;
65 C...6 = Technicolor;
66 C...7 = exotics (Z'/W'/LQ/R/f*/H++/Z_R/W_R/G*).
67  dimension mappr(500)
68  DATA (mappr(i),i=1,180)/
69  & 3, 3, 4, 0, 4, 0, 0, 4, 0, 1,
70  1 1, 1, 1, 1, 3, 3, 0, 1, 3, 3,
71  2 0, 3, 3, 4, 3, 4, 0, 1, 1, 3,
72  3 3, 4, 1, 1, 3, 3, 0, 0, 0, 0,
73  4 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
74  5 0, 0, 1, 1, 0, 0, 0, 1, 0, 0,
75  6 0, 0, 0, 0, 0, 0, 0, 1, 3, 3,
76  7 4, 4, 4, 0, 0, 4, 4, 0, 0, 1,
77  8 2, 2, 2, 2, 2, 2, 2, 2, 2, 0,
78  9 1, 1, 1, 1, 1, 1, 0, 0, 1, 0,
79  & 0, 4, 4, 2, 2, 2, 2, 2, 0, 4,
80  1 4, 4, 4, 1, 1, 0, 0, 0, 0, 0,
81  2 4, 4, 4, 4, 0, 0, 0, 0, 0, 0,
82  3 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
83  4 7, 7, 4, 7, 7, 7, 7, 7, 6, 0,
84  5 4, 4, 4, 0, 0, 4, 4, 4, 0, 0,
85  6 4, 7, 7, 7, 6, 6, 7, 7, 7, 0,
86  7 4, 4, 4, 4, 0, 4, 4, 4, 4, 0/
87  DATA (mappr(i),i=181,500)/
88  8 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
89  9 6, 6, 6, 6, 6, 0, 0, 0, 0, 0,
90  & 100*5,
91  & 5, 0, 0, 0, 0, 0, 0, 0, 0, 0,
92  1 30*0,
93  4 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
94  5 7, 7, 7, 7, 0, 0, 0, 0, 0, 0,
95  6 6, 6, 6, 6, 6, 6, 6, 6, 0, 6,
96  7 6, 6, 6, 6, 6, 6, 6, 0, 0, 0,
97  8 6, 6, 6, 6, 6, 6, 6, 6, 0, 0,
98  9 7, 7, 7, 7, 7, 0, 0, 0, 0, 0,
99  & 4, 4, 18*0,
100  2 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
101  3 2, 2, 2, 2, 2, 2, 2, 2, 2, 0,
102  4 20*0,
103  6 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
104  7 2, 2, 2, 2, 2, 2, 2, 2, 2, 0,
105  8 20*0/
106 
107 C...Reset number of channels and cross-section
108  nchn=0
109  sigs=0d0
110 
111 C...Read process to consider.
112  isub=mint(1)
113  isubsv=isub
114  map=mappr(isub)
115 
116 C...Read kinematical variables and limits
117  istsb=iset(isubsv)
118  taumin=vint(11)
119  ystmin=vint(12)
120  ctnmin=vint(13)
121  ctpmin=vint(14)
122  taupmn=vint(16)
123  tau=vint(21)
124  yst=vint(22)
125  cth=vint(23)
126  xt2=vint(25)
127  taup=vint(26)
128  taumax=vint(31)
129  ystmax=vint(32)
130  ctnmax=vint(33)
131  ctpmax=vint(34)
132  taupmx=vint(36)
133 
134 C...Derive kinematical quantities
135  taue=tau
136  IF(istsb.GE.3.AND.istsb.LE.5) taue=taup
137  x(1)=sqrt(taue)*exp(yst)
138  x(2)=sqrt(taue)*exp(-yst)
139  IF(mint(45).EQ.2.AND.istsb.GE.1) THEN
140  IF(x(1).GT.1d0-1d-7) RETURN
141  ELSEIF(mint(45).EQ.3) THEN
142  x(1)=min(1d0-1.1d-10,x(1))
143  ENDIF
144  IF(mint(46).EQ.2.AND.istsb.GE.1) THEN
145  IF(x(2).GT.1d0-1d-7) RETURN
146  ELSEIF(mint(46).EQ.3) THEN
147  x(2)=min(1d0-1.1d-10,x(2))
148  ENDIF
149  sh=max(1d0,tau*vint(2))
150  sqm3=vint(63)
151  sqm4=vint(64)
152  rm3=sqm3/sh
153  rm4=sqm4/sh
154  be34=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
155  rpts=4d0*vint(71)**2/sh
156  be34l=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4-rpts))
157  rm34=max(1d-20,2d0*rm3*rm4)
158  rsqm=1d0+rm34
159  IF(2d0*vint(71)**2/max(1d0,vint(21)*vint(2)).LT.0.0001d0)
160  &rm34=max(rm34,2d0*vint(71)**2/max(1d0,vint(21)*vint(2)))
161  rthm=(4d0*rm3*rm4+rpts)/(1d0-rm3-rm4+be34l)
162  IF(istsb.EQ.0) THEN
163  th=vint(45)
164  uh=-0.5d0*sh*max(rthm,1d0-rm3-rm4+be34*cth)
165  sqpth=max(vint(71)**2,0.25d0*sh*be34**2*vint(59)**2)
166  ELSE
167 C...Kinematics with incoming masses tricky: now depends on how
168 C...subprocess has been set up w.r.t. order of incoming partons.
169  rm1=0d0
170  IF(mint(15).EQ.22.AND.vint(3).LT.0d0) rm1=-vint(3)**2/sh
171  rm2=0d0
172  IF(mint(16).EQ.22.AND.vint(4).LT.0d0) rm2=-vint(4)**2/sh
173  IF(isub.EQ.35) THEN
174  rm2=min(rm1,rm2)
175  rm1=0d0
176  ENDIF
177  be12=sqrt(max(0d0,(1d0-rm1-rm2)**2-4d0*rm1*rm2))
178  tucom=(1d0-rm1-rm2)*(1d0-rm3-rm4)
179  th=-0.5d0*sh*max(rthm,tucom-2d0*rm1*rm4-2d0*rm2*rm3-
180  & be12*be34*cth)
181  uh=-0.5d0*sh*max(rthm,tucom-2d0*rm1*rm3-2d0*rm2*rm4+
182  & be12*be34*cth)
183  sqpth=max(vint(71)**2,0.25d0*sh*be34**2*(1d0-cth**2))
184  ENDIF
185  shr=sqrt(sh)
186  sh2=sh**2
187  th2=th**2
188  uh2=uh**2
189 
190 C...Choice of Q2 scale for hard process (e.g. alpha_s).
191  IF(istsb.EQ.1.OR.istsb.EQ.3.OR.istsb.EQ.5) THEN
192  q2=sh
193  ELSEIF(istsb.EQ.8) THEN
194  IF(mint(107).EQ.4) q2=vint(307)
195  IF(mint(108).EQ.4) q2=vint(308)
196  ELSEIF(mod(istsb,2).EQ.0.OR.istsb.EQ.9) THEN
197  q2in1=0d0
198  IF(mint(11).EQ.22.AND.vint(3).LT.0d0) q2in1=vint(3)**2
199  q2in2=0d0
200  IF(mint(12).EQ.22.AND.vint(4).LT.0d0) q2in2=vint(4)**2
201  IF(mstp(32).EQ.1) THEN
202  q2=2d0*sh*th*uh/(sh**2+th**2+uh**2)
203  ELSEIF(mstp(32).EQ.2) THEN
204  q2=sqpth+0.5d0*(sqm3+sqm4)
205  ELSEIF(mstp(32).EQ.3) THEN
206  q2=min(-th,-uh)
207  ELSEIF(mstp(32).EQ.4) THEN
208  q2=sh
209  ELSEIF(mstp(32).EQ.5) THEN
210  q2=-th
211  ELSEIF(mstp(32).EQ.6) THEN
212  xsf1=x(1)
213  IF(istsb.EQ.9) xsf1=x(1)/vint(143)
214  xsf2=x(2)
215  IF(istsb.EQ.9) xsf2=x(2)/vint(144)
216  q2=(1d0+xsf1*q2in1/sh+xsf2*q2in2/sh)*
217  & (sqpth+0.5d0*(sqm3+sqm4))
218  ELSEIF(mstp(32).EQ.7) THEN
219  q2=(1d0+q2in1/sh+q2in2/sh)*(sqpth+0.5d0*(sqm3+sqm4))
220  ELSEIF(mstp(32).EQ.8) THEN
221  q2=sqpth+0.5d0*(q2in1+q2in2+sqm3+sqm4)
222  ELSEIF(mstp(32).EQ.9) THEN
223  q2=sqpth+q2in1+q2in2+sqm3+sqm4
224  ELSEIF(mstp(32).EQ.10) THEN
225  q2=vint(2)
226 C..Begin JA 040914
227  ELSEIF(mstp(32).EQ.11) THEN
228  q2=0.25*(sqm3+sqm4+2*sqrt(sqm3*sqm4))
229  ELSEIF(mstp(32).EQ.12) THEN
230  q2=parp(193)
231 C..End JA
232  ELSEIF(mstp(32).EQ.13) THEN
233  q2=sqpth
234  ENDIF
235  IF(mint(35).LE.2.AND.istsb.EQ.9) q2=sqpth
236  IF(istsb.EQ.9.AND.mstp(82).GE.2) q2=q2+
237  & (parp(82)*(vint(1)/parp(89))**parp(90))**2
238  ENDIF
239 
240 C...Choice of Q2 scale for parton densities.
241  q2sf=q2
242 C..Begin JA 040914
243  IF(mstp(32).EQ.12.AND.(mod(istsb,2).EQ.0.OR.istsb.EQ.9)
244  & .OR.mstp(39).EQ.8.AND.(istsb.GE.3.AND.istsb.LE.5))
245  & q2=parp(194)
246 C..End JA
247  IF(istsb.GE.3.AND.istsb.LE.5) THEN
248  q2sf=pmas(23,1)**2
249  IF(isub.EQ.8.OR.isub.EQ.76.OR.isub.EQ.77.OR.isub.EQ.124.OR.
250  & isub.EQ.174.OR.isub.EQ.179.OR.isub.EQ.351) q2sf=pmas(24,1)**2
251  IF(isub.EQ.352) q2sf=pmas(pycomp(9900024),1)**2
252  IF(isub.EQ.121.OR.isub.EQ.122.OR.isub.EQ.181.OR.isub.EQ.182.OR.
253  & isub.EQ.186.OR.isub.EQ.187.OR.isub.EQ.401.OR.isub.EQ.402) THEN
254  q2sf=pmas(pycomp(kfpr(isubsv,2)),1)**2
255  IF(mstp(39).EQ.2) q2sf=
256  & max(vint(201)**2+vint(202),vint(206)**2+vint(207))
257  IF(mstp(39).EQ.3) q2sf=sh
258  IF(mstp(39).EQ.4) q2sf=vint(26)*vint(2)
259  IF(mstp(39).EQ.5) q2sf=pmas(pycomp(kfpr(isubsv,1)),1)**2
260 C..Begin JA 040914
261  IF(mstp(39).EQ.6) q2sf=0.25*(vint(201)+sqrt(sh))**2
262  IF(mstp(39).EQ.7) q2sf=
263  & (vint(201)**2+vint(202)+vint(206)**2+vint(207))/2d0
264  IF(mstp(39).EQ.8) q2sf=parp(193)
265 C..End JA
266  ENDIF
267  ENDIF
268  IF(mint(35).GE.3.AND.istsb.EQ.9) q2sf=sqpth
269 
270  q2ps=q2sf
271  q2sf=q2sf*parp(34)
272  IF(mstp(69).GE.1.AND.mint(47).EQ.5) q2sf=vint(2)
273  IF(mstp(69).GE.2) q2sf=vint(2)
274 
275 C...Identify to which class(es) subprocess belongs
276  ismecr=0
277  isqcd=0
278  isjets=0
279  IF (isubsv.EQ.1.OR.isubsv.EQ.2.OR.isubsv.EQ.3.OR.
280  & isubsv.EQ.102.OR.isubsv.EQ.141.OR.isubsv.EQ.142.OR.
281  & isubsv.EQ.144.OR.isubsv.EQ.151.OR.isubsv.EQ.152.OR.
282  & isubsv.EQ.156.OR.isubsv.EQ.157) ismecr=1
283  IF (isubsv.EQ.11.OR.isubsv.EQ.12.OR.isubsv.EQ.13.OR.
284  & isubsv.EQ.28.OR.isubsv.EQ.53.OR.isubsv.EQ.68) isqcd=1
285  IF ((isubsv.EQ.81.OR.isubsv.EQ.82).AND.mint(55).LE.5) isqcd=1
286  IF (isubsv.GE.381.AND.isubsv.LE.386) isqcd=1
287  IF ((isubsv.EQ.387.OR.isubsv.EQ.388).AND.mint(55).LE.5) isqcd=1
288  IF (istsb.EQ.9) isqcd=1
289  IF ((isubsv.GE.86.AND.isubsv.LE.89).OR.isubsv.EQ.107.OR.
290  & (isubsv.GE.14.AND.isubsv.LE.16).OR.(isubsv.GE.29.AND.
291  & isubsv.LE.32).OR.(isubsv.GE.111.AND.isubsv.LE.113).OR.
292  & isubsv.EQ.115.OR.(isubsv.GE.183.AND.isubsv.LE.185).OR.
293  & (isubsv.GE.188.AND.isubsv.LE.190).OR.isubsv.EQ.161.OR.
294  & isubsv.EQ.167.OR.isubsv.EQ.168.OR.(isubsv.GE.393.AND.
295  & isubsv.LE.395).OR.(isubsv.GE.421.AND.isubsv.LE.439).OR.
296  & (isubsv.GE.461.AND.isubsv.LE.479)) isjets=1
297 C...WBF is special case of ISJETS
298  IF (isubsv.EQ.5.OR.isubsv.EQ.8.OR.
299  & (isubsv.GE.71.AND.isubsv.LE.73).OR.
300  & isubsv.EQ.76.OR.isubsv.EQ.77.OR.
301  & (isubsv.GE.121.AND.isubsv.LE.124).OR.
302  & isubsv.EQ.173.OR.isubsv.EQ.174.OR.
303  & isubsv.EQ.178.OR.isubsv.EQ.179.OR.
304  & isubsv.EQ.181.OR.isubsv.EQ.182.OR.
305  & isubsv.EQ.186.OR.isubsv.EQ.187.OR.
306  & isubsv.EQ.351.OR.isubsv.EQ.352) isjets=2
307 C...Some processes with photons also belong here.
308  IF (isubsv.EQ.10.OR.(isubsv.GE.18.AND.isubsv.LE.20).OR.
309  & (isubsv.GE.33.AND.isubsv.LE.36).OR.isubsv.EQ.54.OR.
310  & isubsv.EQ.58.OR.isubsv.EQ.69.OR.isubsv.EQ.70.OR.
311  & isubsv.EQ.80.OR.(isubsv.GE.83.AND.isubsv.LE.85).OR.
312  & (isubsv.GE.106.AND.isubsv.LE.110).OR.isubsv.EQ.114.OR.
313  & (isubsv.GE.131.AND.isubsv.LE.140)) isjets=3
314 
315 C...Choice of Q2 scale for parton-shower activity.
316  IF(mstp(22).GE.1.AND.(isub.EQ.10.OR.isub.EQ.83).AND.
317  &(mint(43).EQ.2.OR.mint(43).EQ.3)) THEN
318  xbj=x(2)
319  IF(mint(43).EQ.3) xbj=x(1)
320  IF(mstp(22).EQ.1) THEN
321  q2ps=-th
322  ELSEIF(mstp(22).EQ.2) THEN
323  q2ps=((1d0-xbj)/xbj)*(-th)
324  ELSEIF(mstp(22).EQ.3) THEN
325  q2ps=sqrt((1d0-xbj)/xbj)*(-th)
326  ELSE
327  q2ps=(1d0-xbj)*max(1d0,-log(xbj))*(-th)
328  ENDIF
329  ENDIF
330 C...For multiple interactions, start from scale defined above
331 C...For all other QCD or "+jets"-type events, start shower from pThard.
332  IF (isjets.EQ.1.OR.isqcd.EQ.1.AND.istsb.NE.9) q2ps=sqpth
333  IF((mstp(68).EQ.1.OR.mstp(68).EQ.3).AND.ismecr.EQ.1) THEN
334 C...Max shower scale = s for ME corrected processes.
335 C...(pT-ordering: max pT2 is s/4)
336  q2ps=vint(2)
337  IF (mint(35).GE.3) q2ps=q2ps*0.25d0
338  ELSEIF(mstp(68).GE.2.AND.isqcd.EQ.0.AND.isjets.EQ.0) THEN
339 C...Max shower scale = s for all non-QCD, non-"+ jet" type processes.
340 C...(pT-ordering: max pT2 is s/4)
341  q2ps=vint(2)
342  IF (mint(35).GE.3) q2ps=q2ps*0.25d0
343  ENDIF
344  IF(mint(35).EQ.2.AND.istsb.EQ.9) q2ps=sqpth
345 
346 C...Elastic and diffractive events not associated with scales so set 0.
347  IF(isubsv.GE.91.AND.isubsv.LE.94) THEN
348  q2sf=0d0
349  q2ps=0d0
350  ENDIF
351 
352 C...Store derived kinematical quantities
353  vint(41)=x(1)
354  vint(42)=x(2)
355  vint(44)=sh
356  vint(43)=sqrt(sh)
357  vint(45)=th
358  vint(46)=uh
359  IF(istsb.NE.8) vint(48)=sqpth
360  IF(istsb.NE.8) vint(47)=sqrt(sqpth)
361  vint(50)=taup*vint(2)
362  vint(49)=sqrt(max(0d0,vint(50)))
363  vint(52)=q2
364  vint(51)=sqrt(q2)
365  vint(54)=q2sf
366  vint(53)=sqrt(q2sf)
367  vint(56)=q2ps
368  vint(55)=sqrt(q2ps)
369 
370 C...Set starting scale for multiple interactions
371  IF (isubsv.EQ.95) THEN
372  xt2gmx=0d0
373  ELSEIF(mstp(86).EQ.3.OR.(mstp(86).EQ.2.AND.isubsv.NE.11.AND.
374  & isubsv.NE.12.AND.isubsv.NE.13.AND.isubsv.NE.28.AND.
375  & isubsv.NE.53.AND.isubsv.NE.68.AND.isubsv.NE.95.AND.
376  & isubsv.NE.96)) THEN
377 C...All accessible phase space allowed.
378  xt2gmx=(1d0-vint(41))*(1d0-vint(42))
379  ELSE
380 C...Scale of hard process sets limit.
381 C...2 -> 1. Limit is tau = x1*x2.
382 C...2 -> 2. Limit is XT2 for hard process + FS masses.
383 C...2 -> n > 2. Limit is tau' = tau of outer process.
384  xt2gmx=vint(25)
385  IF(istsb.EQ.1) xt2gmx=vint(21)
386  IF(istsb.EQ.2)
387  & xt2gmx=(4d0*vint(48)+2d0*vint(63)+2d0*vint(64))/vint(2)
388  IF(istsb.GE.3.AND.istsb.LE.5) xt2gmx=vint(26)
389  ENDIF
390  vint(62)=0.25d0*xt2gmx*vint(2)
391  vint(61)=sqrt(max(0d0,vint(62)))
392 
393 C...Calculate parton distributions
394  IF(istsb.LE.0) goto 160
395  IF(mint(47).GE.2) THEN
396  DO 110 i=3-min(2,mint(45)),min(2,mint(46))
397  xsf=x(i)
398  IF(istsb.EQ.9) xsf=x(i)/vint(142+i)
399  IF(isub.EQ.99) THEN
400  IF(mint(140+i).EQ.0) THEN
401  xsf=vint(309-i)/(vint(2)+vint(309-i)-vint(i+2)**2)
402  ELSE
403  xsf=vint(309-i)/(vint(2)+vint(307)+vint(308))
404  ENDIF
405  vint(40+i)=xsf
406  q2sf=vint(309-i)
407  ENDIF
408  mint(105)=mint(102+i)
409  mint(109)=mint(106+i)
410  vint(120)=vint(2+i)
411  IF(mstp(57).LE.1) THEN
412  CALL pypdfu(mint(10+i),xsf,q2sf,xpq)
413  ELSE
414  CALL pypdfl(mint(10+i),xsf,q2sf,xpq)
415  ENDIF
416 C...Safety margin against heavy flavour very close to threshold,
417 C...e.g. caused by mismatch in c and b masses.
418  IF(q2sf.LT.1.1*pmas(4,1)**2) THEN
419  xpq(4)=0d0
420  xpq(-4)=0d0
421  ENDIF
422  IF(q2sf.LT.1.1*pmas(5,1)**2) THEN
423  xpq(5)=0d0
424  xpq(-5)=0d0
425  ENDIF
426  DO 100 kfl=-25,25
427  xsfx(i,kfl)=xpq(kfl)
428  100 CONTINUE
429  110 CONTINUE
430  ENDIF
431 
432 C...Calculate alpha_em, alpha_strong and K-factor
433  xw=paru(102)
434  xwv=xw
435  IF(mstp(8).GE.2.OR.(isub.GE.71.AND.isub.LE.77)) xw=
436  &1d0-(pmas(24,1)/pmas(23,1))**2
437  xw1=1d0-xw
438  xwc=1d0/(16d0*xw*xw1)
439  aem=pyalem(q2)
440  IF(mstp(8).GE.1) aem=sqrt(2d0)*paru(105)*pmas(24,1)**2*xw/paru(1)
441  IF(mstp(33).NE.3) as=pyalps(parp(34)*q2)
442  fack=1d0
443  faca=1d0
444  IF(mstp(33).EQ.1) THEN
445  fack=parp(31)
446  ELSEIF(mstp(33).EQ.2) THEN
447  fack=parp(31)
448  faca=parp(32)/parp(31)
449  ELSEIF(mstp(33).EQ.3) THEN
450  q2as=parp(33)*q2
451  IF(istsb.EQ.9.AND.mstp(82).GE.2) q2as=q2as+
452  & paru(112)*parp(82)*(vint(1)/parp(89))**parp(90)
453  as=pyalps(q2as)
454  ENDIF
455  vint(138)=1d0
456  vint(57)=aem
457  vint(58)=as
458 
459 C...Set flags for allowed reacting partons/leptons
460  DO 140 i=1,2
461  DO 120 j=-25,25
462  kfac(i,j)=0
463  120 CONTINUE
464  IF(mint(44+i).EQ.1) THEN
465  kfac(i,mint(10+i))=1
466  ELSEIF(mint(40+i).EQ.1.AND.mstp(12).EQ.0) THEN
467  kfac(i,mint(10+i))=1
468  kfac(i,22)=1
469  kfac(i,24)=1
470  kfac(i,-24)=1
471  ELSE
472  DO 130 j=-25,25
473  kfac(i,j)=kfin(i,j)
474  IF(iabs(j).GT.mstp(58).AND.iabs(j).LE.10) kfac(i,j)=0
475  IF(xsfx(i,j).LT.1d-10) kfac(i,j)=0
476  130 CONTINUE
477  ENDIF
478  140 CONTINUE
479 
480 C...Lower and upper limit for fermion flavour loops
481  mmin1=0
482  mmax1=0
483  mmin2=0
484  mmax2=0
485  DO 150 j=-20,20
486  IF(kfac(1,-j).EQ.1) mmin1=-j
487  IF(kfac(1,j).EQ.1) mmax1=j
488  IF(kfac(2,-j).EQ.1) mmin2=-j
489  IF(kfac(2,j).EQ.1) mmax2=j
490  150 CONTINUE
491  mmina=min(mmin1,mmin2)
492  mmaxa=max(mmax1,mmax2)
493 
494 C...Common resonance mass and width combinations
495  sqmz=pmas(23,1)**2
496  sqmw=pmas(24,1)**2
497  gmmz=pmas(23,1)*pmas(23,2)
498  gmmw=pmas(24,1)*pmas(24,2)
499 
500 C...Polarization factors...implemented so far for W+W-(25)
501  polr=(1d0+parj(132))*(1d0-parj(131))
502  poll=(1d0-parj(132))*(1d0+parj(131))
503  polrr=(1d0+parj(132))*(1d0+parj(131))
504  polll=(1d0-parj(132))*(1d0-parj(131))
505 
506 C...Phase space integral in tau
507  comfac=paru(1)*paru(5)/vint(2)
508  IF(mint(41).EQ.2.AND.mint(42).EQ.2) comfac=comfac*fack
509  IF((mint(47).GE.2.OR.(istsb.GE.3.AND.istsb.LE.5)).AND.
510  &istsb.NE.8.AND.istsb.NE.9) THEN
511  atau1=log(taumax/taumin)
512  atau2=(taumax-taumin)/(taumax*taumin)
513  h1=coef(isubsv,1)+(atau1/atau2)*coef(isubsv,2)/tau
514  IF(mint(72).GE.1) THEN
515  taur1=vint(73)
516  gamr1=vint(74)
517  ataud=log(taumax/taumin*(taumin+taur1)/(taumax+taur1))
518  atau3=ataud/taur1
519  IF(ataud.GT.1d-10) h1=h1+
520  & (atau1/atau3)*coef(isubsv,3)/(tau+taur1)
521  ataud=atan((taumax-taur1)/gamr1)-atan((taumin-taur1)/gamr1)
522  atau4=ataud/gamr1
523  IF(ataud.GT.1d-10) h1=h1+
524  & (atau1/atau4)*coef(isubsv,4)*tau/((tau-taur1)**2+gamr1**2)
525  ENDIF
526  IF(mint(72).EQ.2) THEN
527  taur2=vint(75)
528  gamr2=vint(76)
529  ataud=log(taumax/taumin*(taumin+taur2)/(taumax+taur2))
530  atau5=ataud/taur2
531  IF(ataud.GT.1d-10) h1=h1+
532  & (atau1/atau5)*coef(isubsv,5)/(tau+taur2)
533  ataud=atan((taumax-taur2)/gamr2)-atan((taumin-taur2)/gamr2)
534  atau6=ataud/gamr2
535  IF(ataud.GT.1d-10) h1=h1+
536  & (atau1/atau6)*coef(isubsv,6)*tau/((tau-taur2)**2+gamr2**2)
537  ENDIF
538  IF(mint(47).EQ.5.AND.(istsb.LE.2.OR.istsb.GE.5)) THEN
539  atau7=log(max(2d-10,1d0-taumin)/max(2d-10,1d0-taumax))
540  IF(atau7.GT.1d-10) h1=h1+(atau1/atau7)*coef(isubsv,7)*tau/
541  & max(2d-10,1d0-tau)
542  ELSEIF(mint(47).GE.6.AND.(istsb.LE.2.OR.istsb.GE.5)) THEN
543  atau7=log(max(1d-10,1d0-taumin)/max(1d-10,1d0-taumax))
544  IF(atau7.GT.1d-10) h1=h1+(atau1/atau7)*coef(isubsv,7)*tau/
545  & max(1d-10,1d0-tau)
546  ENDIF
547  comfac=comfac*atau1/(tau*h1)
548  ENDIF
549 
550 C...Phase space integral in y*
551  IF((mint(47).EQ.4.OR.mint(47).EQ.5).AND.istsb.NE.8.AND.istsb.NE.9)
552  &THEN
553  ayst0=ystmax-ystmin
554  IF(ayst0.LT.1d-10) THEN
555  comfac=0d0
556  ELSE
557  ayst1=0.5d0*(ystmax-ystmin)**2
558  ayst2=ayst1
559  ayst3=2d0*(atan(exp(ystmax))-atan(exp(ystmin)))
560  h2=(ayst0/ayst1)*coef(isubsv,8)*(yst-ystmin)+
561  & (ayst0/ayst2)*coef(isubsv,9)*(ystmax-yst)+
562  & (ayst0/ayst3)*coef(isubsv,10)/cosh(yst)
563  IF(mint(45).EQ.3) THEN
564  yst0=-0.5d0*log(taue)
565  ayst4=log(max(1d-10,exp(yst0-ystmin)-1d0)/
566  & max(1d-10,exp(yst0-ystmax)-1d0))
567  IF(ayst4.GT.1d-10) h2=h2+(ayst0/ayst4)*coef(isubsv,11)/
568  & max(1d-10,1d0-exp(yst-yst0))
569  ENDIF
570  IF(mint(46).EQ.3) THEN
571  yst0=-0.5d0*log(taue)
572  ayst5=log(max(1d-10,exp(yst0+ystmax)-1d0)/
573  & max(1d-10,exp(yst0+ystmin)-1d0))
574  IF(ayst5.GT.1d-10) h2=h2+(ayst0/ayst5)*coef(isubsv,12)/
575  & max(1d-10,1d0-exp(-yst-yst0))
576  ENDIF
577  comfac=comfac*ayst0/h2
578  ENDIF
579  ENDIF
580 
581 C...2 -> 1 processes: reduction in angular part of phase space integral
582 C...for case of decaying resonance
583  acth0=ctnmax-ctnmin+ctpmax-ctpmin
584  IF((istsb.EQ.1.OR.istsb.EQ.3.OR.istsb.EQ.5)) THEN
585  IF(mdcy(pycomp(kfpr(isubsv,1)),1).EQ.1) THEN
586  IF(kfpr(isub,1).EQ.25.OR.kfpr(isub,1).EQ.37.OR.
587  & kfpr(isub,1).EQ.39) THEN
588  comfac=comfac*0.5d0*acth0
589  ELSE
590  comfac=comfac*0.125d0*(3d0*acth0+ctnmax**3-ctnmin**3+
591  & ctpmax**3-ctpmin**3)
592  ENDIF
593  ENDIF
594 
595 C...2 -> 2 processes: angular part of phase space integral
596  ELSEIF(istsb.EQ.2.OR.istsb.EQ.4) THEN
597  acth1=log((max(rm34,rsqm-ctnmin)*max(rm34,rsqm-ctpmin))/
598  & (max(rm34,rsqm-ctnmax)*max(rm34,rsqm-ctpmax)))
599  acth2=log((max(rm34,rsqm+ctnmax)*max(rm34,rsqm+ctpmax))/
600  & (max(rm34,rsqm+ctnmin)*max(rm34,rsqm+ctpmin)))
601  acth3=1d0/max(rm34,rsqm-ctnmax)-1d0/max(rm34,rsqm-ctnmin)+
602  & 1d0/max(rm34,rsqm-ctpmax)-1d0/max(rm34,rsqm-ctpmin)
603  acth4=1d0/max(rm34,rsqm+ctnmin)-1d0/max(rm34,rsqm+ctnmax)+
604  & 1d0/max(rm34,rsqm+ctpmin)-1d0/max(rm34,rsqm+ctpmax)
605  h3=coef(isubsv,13)+
606  & (acth0/acth1)*coef(isubsv,14)/max(rm34,rsqm-cth)+
607  & (acth0/acth2)*coef(isubsv,15)/max(rm34,rsqm+cth)+
608  & (acth0/acth3)*coef(isubsv,16)/max(rm34,rsqm-cth)**2+
609  & (acth0/acth4)*coef(isubsv,17)/max(rm34,rsqm+cth)**2
610  comfac=comfac*acth0*0.5d0*be34/h3
611 
612 C...2 -> 2 processes: take into account final state Breit-Wigners
613  comfac=comfac*vint(80)
614  ENDIF
615 
616 C...2 -> 3, 4 processes: phace space integral in tau'
617  IF(mint(47).GE.2.AND.istsb.GE.3.AND.istsb.LE.5) THEN
618  ataup1=log(taupmx/taupmn)
619  ataup2=((1d0-tau/taupmx)**4-(1d0-tau/taupmn)**4)/(4d0*tau)
620  h4=coef(isubsv,18)+
621  & (ataup1/ataup2)*coef(isubsv,19)*(1d0-tau/taup)**3/taup
622  IF(mint(47).EQ.5) THEN
623  ataup3=log(max(2d-10,1d0-taupmn)/max(2d-10,1d0-taupmx))
624  h4=h4+(ataup1/ataup3)*coef(isubsv,20)*taup/max(2d-10,1d0-taup)
625  ELSEIF(mint(47).GE.6) THEN
626  ataup3=log(max(1d-10,1d0-taupmn)/max(1d-10,1d0-taupmx))
627  h4=h4+(ataup1/ataup3)*coef(isubsv,20)*taup/max(1d-10,1d0-taup)
628  ENDIF
629  comfac=comfac*ataup1/h4
630  ENDIF
631 
632 C...2 -> 3, 4 processes: effective W/Z parton distributions
633  IF(istsb.EQ.3.OR.istsb.EQ.4) THEN
634  IF(1d0-tau/taup.GT.1d-4) THEN
635  fzw=(1d0+tau/taup)*log(taup/tau)-2d0*(1d0-tau/taup)
636  ELSE
637  fzw=1d0/6d0*(1d0-tau/taup)**3*tau/taup
638  ENDIF
639  comfac=comfac*fzw
640  ENDIF
641 
642 C...2 -> 3 processes: phase space integrals for pT1, pT2, y3, mirror
643  IF(istsb.EQ.5) THEN
644  comfac=comfac*vint(205)*vint(210)*vint(212)*vint(214)/
645  & (128d0*paru(1)**4*vint(220))*(tau**2/taup)
646  ENDIF
647 
648 C...Phase space integral for low-pT and multiple interactions
649  IF(istsb.EQ.9) THEN
650  comfac=paru(1)*paru(5)*fack*0.5d0*vint(2)/sh2
651  atau1=log(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)
652  atau2=2d0*atan(1d0/xt2-1d0)/sqrt(xt2)
653  h1=coef(isubsv,1)+(atau1/atau2)*coef(isubsv,2)/sqrt(tau)
654  comfac=comfac*atau1/h1
655  ayst0=ystmax-ystmin
656  ayst1=0.5d0*(ystmax-ystmin)**2
657  ayst3=2d0*(atan(exp(ystmax))-atan(exp(ystmin)))
658  h2=(ayst0/ayst1)*coef(isubsv,8)*(yst-ystmin)+
659  & (ayst0/ayst1)*coef(isubsv,9)*(ystmax-yst)+
660  & (ayst0/ayst3)*coef(isubsv,10)/cosh(yst)
661  comfac=comfac*ayst0/h2
662  IF(mstp(82).LE.1) comfac=comfac*xt2**2*(1d0/vint(149)-1d0)
663 C...For MSTP(82)>=2 an additional factor (xT2/(xT2+VINT(149))**2 is
664 C...introduced to make cross-section finite for xT2 -> 0
665  IF(mstp(82).GE.2) comfac=comfac*xt2**2/(vint(149)*
666  & (1d0+vint(149)))
667  ENDIF
668 
669 C...Real gamma + gamma: include factor 2 when different nature
670  160 IF(mint(11).EQ.22.AND.mint(12).EQ.22.AND.mint(123).GE.4.AND.
671  &mstp(14).LE.10) comfac=2d0*comfac
672 
673 C...Extra factors to include the effects of
674 C...longitudinal resolved photons (but not direct or DIS ones).
675  DO 170 isde=1,2
676  IF(mint(10+isde).EQ.22.AND.mint(106+isde).GE.1.AND.
677  & mint(106+isde).LE.3) THEN
678  vint(314+isde)=1d0
679  xy=parp(166+isde)
680  IF(mstp(16).EQ.0) THEN
681  IF(vint(304+isde).GT.0d0.AND.vint(304+isde).LT.1d0)
682  & xy=vint(304+isde)
683  ELSE
684  IF(vint(308+isde).GT.0d0.AND.vint(308+isde).LT.1d0)
685  & xy=vint(308+isde)
686  ENDIF
687  q2ga=vint(306+isde)
688  IF(mstp(17).GT.0.AND.xy.GT.0d0.AND.xy.LT.1d0.AND.
689  & q2ga.GT.0d0) THEN
690  reduce=0d0
691  IF(mstp(17).EQ.1) THEN
692  reduce=4d0*q2*q2ga/(q2+q2ga)**2
693  ELSEIF(mstp(17).EQ.2) THEN
694  reduce=4d0*q2ga/(q2+q2ga)
695  ELSEIF(mstp(17).EQ.3) THEN
696  pmvirt=pmas(pycomp(113),1)
697  reduce=4d0*q2ga/(pmvirt**2+q2ga)
698  ELSEIF(mstp(17).EQ.4.AND.mint(106+isde).EQ.1) THEN
699  pmvirt=pmas(pycomp(113),1)
700  reduce=4d0*pmvirt**2*q2ga/(pmvirt**2+q2ga)**2
701  ELSEIF(mstp(17).EQ.4.AND.mint(106+isde).EQ.2) THEN
702  pmvirt=pmas(pycomp(113),1)
703  reduce=4d0*pmvirt**2*q2ga/(pmvirt**2+q2ga)**2
704  ELSEIF(mstp(17).EQ.4.AND.mint(106+isde).EQ.3) THEN
705  pmvsmn=4d0*parp(15)**2
706  pmvsmx=4d0*vint(154)**2
707  redtra=1d0/(pmvsmn+q2ga)-1d0/(pmvsmx+q2ga)
708  redlon=(3d0*pmvsmn+q2ga)/(pmvsmn+q2ga)**3-
709  & (3d0*pmvsmx+q2ga)/(pmvsmx+q2ga)**3
710  reduce=4d0*(q2ga/6d0)*redlon/redtra
711  ELSEIF(mstp(17).EQ.5.AND.mint(106+isde).EQ.1) THEN
712  pmvirt=pmas(pycomp(113),1)
713  reduce=4d0*q2ga/(pmvirt**2+q2ga)
714  ELSEIF(mstp(17).EQ.5.AND.mint(106+isde).EQ.2) THEN
715  pmvirt=pmas(pycomp(113),1)
716  reduce=4d0*q2ga/(pmvirt**2+q2ga)
717  ELSEIF(mstp(17).EQ.5.AND.mint(106+isde).EQ.3) THEN
718  pmvsmn=4d0*parp(15)**2
719  pmvsmx=4d0*vint(154)**2
720  redtra=1d0/(pmvsmn+q2ga)-1d0/(pmvsmx+q2ga)
721  redlon=1d0/(pmvsmn+q2ga)**2-1d0/(pmvsmx+q2ga)**2
722  reduce=4d0*(q2ga/2d0)*redlon/redtra
723 C ........Hermes version of R_VMD, this still needs work to
724 C only apply it for process 91 and separate rho from the rest.
725  ELSEIF(mstp(17).EQ.6.AND.mint(106+isde).EQ.2) THEN
726  IF ((mint(1).le.94).and.(mint(1).gt.90)) then
727  if (vint(67).gt.0.0) then
728  if (vint(67).eq.pmas(pycomp(113),1)) then
729  reduce=(q2ga/vint(67)**2)**parp(166)
730  else
731  reduce=q2ga/vint(67)**2
732  endif
733  else
734  pmvirt=pmas(pycomp(113),1)
735  reduce=q2ga/pmvirt**2
736  endif
737  ENDIF
738  ENDIF
739  beamas=pymass(11)
740  IF(vint(302+isde).GT.0d0) beamas=vint(302+isde)
741  IF((mint(11).EQ.22).and.
742  & (mint(12).EQ.2212.or.mint(12).EQ.2112)) THEN
743  fraclt=1d0/(1d0+(xy**2*(1d0-2d0*beamas**2/q2ga))/
744  & (2d0/(1d0+q2ga/xy**2/ebeamenucl**2)*(1d0-xy-
745  & (q2ga/4d0/ebeamenucl**2))))
746  ELSE
747  fraclt=1d0/(1d0+xy**2/2d0/(1d0-xy)*
748  & (1d0-2d0*beamas**2/q2ga))
749  ENDIF
750  vint(314+isde)=1d0+parp(165)*reduce*fraclt
751  ENDIF
752  ELSE
753  vint(314+isde)=1d0
754  ENDIF
755  comfac=comfac*vint(314+isde)
756  170 CONTINUE
757 
758 C...Evaluate cross sections - done in separate routines by kind
759 C...of physics, to keep PYSIGH of sensible size.
760  IF(map.EQ.1) THEN
761 C...Standard QCD (including photons).
762  CALL pysgqc(nchn,sigs)
763  ELSEIF(map.EQ.2) THEN
764 C...Heavy flavours.
765  CALL pysghf(nchn,sigs)
766  ELSEIF(map.EQ.3) THEN
767 C...W/Z.
768  CALL pysgwz(nchn,sigs)
769  ELSEIF(map.EQ.4) THEN
770 C...Higgs (2 doublets; including longitudinal W/Z scattering).
771  CALL pysghg(nchn,sigs)
772  ELSEIF(map.EQ.5) THEN
773 C...SUSY.
774  CALL pysgsu(nchn,sigs)
775  ELSEIF(map.EQ.6) THEN
776 C...Technicolor.
777  CALL pysgtc(nchn,sigs)
778  ELSEIF(map.EQ.7) THEN
779 C...Exotics (Z'/W'/LQ/R/f*/H++/Z_R/W_R/G*).
780  CALL pysgex(nchn,sigs)
781  ENDIF
782 
783 C...Multiply with parton distributions
784  IF(isub.LE.90.OR.isub.GE.96) THEN
785  DO 180 ichn=1,nchn
786  IF(mint(45).GE.2) THEN
787  kfl1=isig(ichn,1)
788  sigh(ichn)=sigh(ichn)*xsfx(1,kfl1)
789  ENDIF
790  IF(mint(46).GE.2) THEN
791  kfl2=isig(ichn,2)
792  sigh(ichn)=sigh(ichn)*xsfx(2,kfl2)
793  ENDIF
794  sigs=sigs+sigh(ichn)
795  180 CONTINUE
796  ENDIF
797 
798  RETURN
799  END