Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyxtot.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyxtot.f
1 
2 C*********************************************************************
3 
4 C...PYXTOT
5 C...Parametrizes total, elastic and diffractive cross-sections
6 C...for different energies and beams. Donnachie-Landshoff for
7 C...total and Schuler-Sjostrand for elastic and diffractive.
8 C...Process code IPROC:
9 C...= 1 : p + p;
10 C...= 2 : pbar + p;
11 C...= 3 : pi+ + p;
12 C...= 4 : pi- + p;
13 C...= 5 : pi0 + p;
14 C...= 6 : phi + p;
15 C...= 7 : J/psi + p;
16 C...= 11 : rho + rho;
17 C...= 12 : rho + phi;
18 C...= 13 : rho + J/psi;
19 C...= 14 : phi + phi;
20 C...= 15 : phi + J/psi;
21 C...= 16 : J/psi + J/psi;
22 C...= 21 : gamma + p (DL);
23 C...= 22 : gamma + p (VDM).
24 C...= 23 : gamma + pi (DL);
25 C...= 24 : gamma + pi (VDM);
26 C...= 25 : gamma + gamma (DL);
27 C...= 26 : gamma + gamma (VDM).
28 
29  SUBROUTINE pyxtot
30 
31 C...Double precision and integer declarations.
32  IMPLICIT DOUBLE PRECISION(a-h, o-z)
33  IMPLICIT INTEGER(i-n)
34  INTEGER pyk,pychge,pycomp
35 C...Commonblocks.
36  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
37  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
38  common/pypars/mstp(200),parp(200),msti(200),pari(200)
39  common/pyint1/mint(400),vint(400)
40  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
41  common/pyint7/sigt(0:6,0:6,0:5)
42  SAVE /pydat1/,/pydat2/,/pypars/,/pyint1/,/pyint5/,/pyint7/
43 C...Local arrays.
44  dimension nproc(30),xpar(30),ypar(30),ihada(20),ihadb(20),
45  &pmhad(4),bhad(4),betp(4),ifitsd(20),ifitdd(20),ceffs(10,8),
46  &ceffd(10,9),sigtmp(6,0:5),expo(4)
47 
48 C...Common constants.
49  DATA eps/0.0808d0/, eta/-0.4525d0/, alp/0.25d0/, cres/2d0/,
50  &pmrc/1.062d0/, smp/0.880d0/, facel/0.0511d0/, facsd/0.0336d0/,
51  &facdd/0.0084d0/
52 
53 C...Number of multiple processes to be evaluated (= 0 : undefined).
54  DATA nproc/7*1,3*0,6*1,4*0,4*3,2*6,4*0/
55 C...X and Y parameters of sigmatot = X * s**epsilon + Y * s**(-eta).
56  DATA xpar/2*21.70d0,3*13.63d0,10.01d0,0.970d0,3*0d0,
57  &8.56d0,6.29d0,0.609d0,4.62d0,0.447d0,0.0434d0,4*0d0,
58  &0.0677d0,0.0534d0,0.0425d0,0.0335d0,2.11d-4,1.31d-4,4*0d0/
59  DATA ypar/
60  &56.08d0,98.39d0,27.56d0,36.02d0,31.79d0,-1.51d0,-0.146d0,3*0d0,
61  &13.08d0,-0.62d0,-0.060d0,0.030d0,-0.0028d0,0.00028d0,4*0d0,
62  &0.129d0,0.115d0,0.081d0,0.072d0,2.15d-4,1.70d-4,4*0d0/
63 
64 C...Beam and target hadron class:
65 C...= 1 : p/n ; = 2 : pi/rho/omega; = 3 : phi; = 4 : J/psi.
66  DATA ihada/2*1,3*2,3,4,3*0,3*2,2*3,4,4*0/
67  DATA ihadb/7*1,3*0,2,3,4,3,2*4,4*0/
68 C...Characteristic class masses, slope parameters, beta = sqrt(X).
69  DATA pmhad/0.938d0,0.770d0,1.020d0,3.097d0/
70  DATA bhad/2.3d0,1.4d0,1.4d0,0.23d0/
71  DATA betp/4.658d0,2.926d0,2.149d0,0.208d0/
72  DATA expo/2.575d0,2.575d0,2.4d0,2.55d0/
73 
74 C...Fitting constants used in parametrizations of diffractive results.
75  DATA ifitsd/2*1,3*2,3,4,3*0,5,6,7,8,9,10,4*0/
76  DATA ifitdd/2*1,3*2,3,4,3*0,5,6,7,8,9,10,4*0/
77  DATA ((ceffs(j1,j2),j2=1,8),j1=1,10)/
78  &0.213d0, 0.0d0, -0.47d0, 150d0, 0.213d0, 0.0d0, -0.47d0, 150d0,
79  &0.213d0, 0.0d0, -0.47d0, 150d0, 0.267d0, 0.0d0, -0.47d0, 100d0,
80  &0.213d0, 0.0d0, -0.47d0, 150d0, 0.232d0, 0.0d0, -0.47d0, 110d0,
81  &0.213d0, 7.0d0, -0.55d0, 800d0, 0.115d0, 0.0d0, -0.47d0, 110d0,
82  &0.267d0, 0.0d0, -0.46d0, 75d0, 0.267d0, 0.0d0, -0.46d0, 75d0,
83  &0.232d0, 0.0d0, -0.46d0, 85d0, 0.267d0, 0.0d0, -0.48d0, 100d0,
84  &0.115d0, 0.0d0, -0.50d0, 90d0, 0.267d0, 6.0d0, -0.56d0, 420d0,
85  &0.232d0, 0.0d0, -0.48d0, 110d0, 0.232d0, 0.0d0, -0.48d0, 110d0,
86  &0.115d0, 0.0d0, -0.52d0, 120d0, 0.232d0, 6.0d0, -0.56d0, 470d0,
87  &0.115d0, 5.5d0, -0.58d0, 570d0, 0.115d0, 5.5d0, -0.58d0, 570d0/
88  DATA ((ceffd(j1,j2),j2=1,9),j1=1,10)/
89  &3.11d0, -7.34d0, 9.71d0, 0.068d0, -0.42d0, 1.31d0,
90  &-1.37d0, 35.0d0, 118d0, 3.11d0, -7.10d0, 10.6d0,
91  &0.073d0, -0.41d0, 1.17d0, -1.41d0, 31.6d0, 95d0,
92  &3.12d0, -7.43d0, 9.21d0, 0.067d0, -0.44d0, 1.41d0,
93  &-1.35d0, 36.5d0, 132d0, 3.13d0, -8.18d0, -4.20d0,
94  &0.056d0, -0.71d0, 3.12d0, -1.12d0, 55.2d0, 1298d0,
95  &3.11d0, -6.90d0, 11.4d0, 0.078d0, -0.40d0, 1.05d0,
96  &-1.40d0, 28.4d0, 78d0, 3.11d0, -7.13d0, 10.0d0,
97  &0.071d0, -0.41d0, 1.23d0, -1.34d0, 33.1d0, 105d0,
98  &3.12d0, -7.90d0, -1.49d0, 0.054d0, -0.64d0, 2.72d0,
99  &-1.13d0, 53.1d0, 995d0, 3.11d0, -7.39d0, 8.22d0,
100  &0.065d0, -0.44d0, 1.45d0, -1.36d0, 38.1d0, 148d0,
101  &3.18d0, -8.95d0, -3.37d0, 0.057d0, -0.76d0, 3.32d0,
102  &-1.12d0, 55.6d0, 1472d0, 4.18d0, -29.2d0, 56.2d0,
103  &0.074d0, -1.36d0, 6.67d0, -1.14d0, 116.2d0, 6532d0/
104 
105 C...Parameters. Combinations of the energy.
106  aem=paru(101)
107  pmth=parp(102)
108  s=vint(2)
109  srt=vint(1)
110  seps=s**eps
111  seta=s**eta
112  slog=log(s)
113 
114 C...Ratio of gamma/pi (for rescaling in parton distributions).
115  vint(281)=(xpar(22)*seps+ypar(22)*seta)/
116  &(xpar(5)*seps+ypar(5)*seta)
117  vint(317)=1d0
118  IF(mint(50).NE.1) RETURN
119 
120 C...Order flavours of incoming particles: KF1 < KF2.
121  IF(iabs(mint(11)).LE.iabs(mint(12))) THEN
122  kf1=iabs(mint(11))
123  kf2=iabs(mint(12))
124  iord=1
125  ELSE
126  kf1=iabs(mint(12))
127  kf2=iabs(mint(11))
128  iord=2
129  ENDIF
130  isgn12=isign(1,mint(11)*mint(12))
131 
132 C...Find process number (for lookup tables).
133  IF(kf1.GT.1000) THEN
134  iproc=1
135  IF(isgn12.LT.0) iproc=2
136  ELSEIF(kf1.GT.100.AND.kf2.GT.1000) THEN
137  iproc=3
138  IF(isgn12.LT.0) iproc=4
139  IF(kf1.EQ.111) iproc=5
140  ELSEIF(kf1.GT.100) THEN
141  iproc=11
142  ELSEIF(kf2.GT.1000) THEN
143  iproc=21
144  IF(mint(123).EQ.2.OR.mint(123).EQ.3) iproc=22
145  ELSEIF(kf2.GT.100) THEN
146  iproc=23
147  IF(mint(123).EQ.2.OR.mint(123).EQ.3) iproc=24
148  ELSE
149  iproc=25
150  IF(mint(123).EQ.2.OR.mint(123).EQ.3.OR.mint(123).EQ.7) iproc=26
151  ENDIF
152 
153 C... Number of multiple processes to be stored; beam/target side.
154  npr=nproc(iproc)
155  mint(101)=1
156  mint(102)=1
157  IF(npr.EQ.3) THEN
158  mint(100+iord)=4
159  ELSEIF(npr.EQ.6) THEN
160  mint(101)=4
161  mint(102)=4
162  ENDIF
163  n1=0
164  IF(mint(101).EQ.4) n1=4
165  n2=0
166  IF(mint(102).EQ.4) n2=4
167 
168 C...Do not do any more for user-set or undefined cross-sections.
169  IF(mstp(31).LE.0) RETURN
170  IF(npr.EQ.0) CALL pyerrm(26,
171  &'(PYXTOT:) cross section for this process not yet implemented')
172 
173 C...Parameters. Combinations of the energy.
174  aem=paru(101)
175  pmth=parp(102)
176  s=vint(2)
177  srt=vint(1)
178  seps=s**eps
179  seta=s**eta
180  slog=log(s)
181 
182 C...Loop over multiple processes (for VDM).
183  DO 110 i=1,npr
184  IF(npr.EQ.1) THEN
185  ipr=iproc
186  ELSEIF(npr.EQ.3) THEN
187  ipr=i+4
188  IF(kf2.LT.1000) ipr=i+10
189  ELSEIF(npr.EQ.6) THEN
190  ipr=i+10
191  ENDIF
192 
193 C...Evaluate hadron species, mass, slope contribution and fit number.
194  iha=ihada(ipr)
195  ihb=ihadb(ipr)
196  pma=pmhad(iha)
197  pmb=pmhad(ihb)
198  bha=bhad(iha)
199  bhb=bhad(ihb)
200  isd=ifitsd(ipr)
201  idd=ifitdd(ipr)
202 
203 C...Skip if energy too low relative to masses.
204  DO 100 j=0,5
205  sigtmp(i,j)=0d0
206  100 CONTINUE
207  IF(srt.LT.pma+pmb+parp(104)) goto 110
208 
209 C...Total cross-section. Elastic slope parameter and cross-section.
210 C...change of elastic slope parameter for rho, phi and J/psi (08/2010)
211 C SIGTMP(I,0)=XPAR(IPR)*SEPS+YPAR(IPR)*SETA
212 C If we include the q2 dependence of the deltafor rho, we need also
213 C to change xpar
214 C seps=s**(0.173+0.068*dlog(VINT(307)/(PMVIRT**2)+1.D0))
215  IF(iha.eq.2) then
216  pmvirt=pmas(pycomp(113),1)
217  sigtmp(i,0)=xpar(ipr)*seps+ypar(ipr)*seta
218  bel=5.84/(1+(parp(165))*(vint(307)/(pmvirt**2))**parp(166))+4.5
219  ELSEIF(iha.eq.3) then
220  pmvirt=pmas(pycomp(333),1)
221  sigtmp(i,0)=xpar(ipr)*seps+ypar(ipr)*seta
222  bel=5.84/(1+(parp(165))*(vint(307)/(pmvirt**2))**parp(166))+4.5
223  elseif(iha.eq.4) then
224  pmvirt=pmas(pycomp(443),1)
225  seps=s**0.2d0
226  bel=4.5
227  sigtmp(i,0)=xpar(ipr)*seps*0.550d0
228  ELSE
229  sigtmp(i,0)=xpar(ipr)*seps+ypar(ipr)*seta
230  bel=2d0*bha+2d0*bhb+4d0*seps-4.2d0
231  write(*,*), iha, ipr, sigtmp(i,0), i
232  ENDIF
233  sigtmp(i,1)=facel*sigtmp(i,0)**2/bel
234 
235 C...Diffractive scattering A + B -> X + B.
236  bsd=2d0*bhb
237  sqml=(pma+pmth)**2
238  sqmu=s*ceffs(isd,1)+ceffs(isd,2)
239  sum1=log((bsd+2d0*alp*log(s/sqml))/
240  & (bsd+2d0*alp*log(s/sqmu)))/(2d0*alp)
241  bxb=ceffs(isd,3)+ceffs(isd,4)/s
242  sum2=cres*log(1d0+((pma+pmrc)/(pma+pmth))**2)/
243  & (bsd+2d0*alp*log(s/((pma+pmth)*(pma+pmrc)))+bxb)
244  sigtmp(i,2)=facsd*xpar(ipr)*betp(ihb)*max(0d0,sum1+sum2)
245 
246 C...Diffractive scattering A + B -> A + X.
247  bsd=2d0*bha
248  sqml=(pmb+pmth)**2
249  sqmu=s*ceffs(isd,5)+ceffs(isd,6)
250  sum1=log((bsd+2d0*alp*log(s/sqml))/
251  & (bsd+2d0*alp*log(s/sqmu)))/(2d0*alp)
252  bax=ceffs(isd,7)+ceffs(isd,8)/s
253  sum2=cres*log(1d0+((pmb+pmrc)/(pmb+pmth))**2)/
254  & (bsd+2d0*alp*log(s/((pmb+pmth)*(pmb+pmrc)))+bax)
255  sigtmp(i,3)=facsd*xpar(ipr)*betp(iha)*max(0d0,sum1+sum2)
256 
257 C...Order single diffractive correctly.
258  IF(iord.EQ.2) THEN
259  sigsav=sigtmp(i,2)
260  sigtmp(i,2)=sigtmp(i,3)
261  sigtmp(i,3)=sigsav
262  ENDIF
263 
264 C...Double diffractive scattering A + B -> X1 + X2.
265  yeff=log(s*smp/((pma+pmth)*(pmb+pmth))**2)
266  deff=ceffd(idd,1)+ceffd(idd,2)/slog+ceffd(idd,3)/slog**2
267  sum1=(deff+yeff*(log(max(1d-10,yeff/deff))-1d0))/(2d0*alp)
268  IF(yeff.LE.0) sum1=0d0
269  sqmu=s*(ceffd(idd,4)+ceffd(idd,5)/slog+ceffd(idd,6)/slog**2)
270  slup=log(max(1.1d0,s/(alp*(pma+pmth)**2*(pmb+pmth)*(pmb+pmrc))))
271  sldn=log(max(1.1d0,s/(alp*sqmu*(pmb+pmth)*(pmb+pmrc))))
272  sum2=cres*log(1d0+((pmb+pmrc)/(pmb+pmth))**2)*log(slup/sldn)/
273  & (2d0*alp)
274  slup=log(max(1.1d0,s/(alp*(pmb+pmth)**2*(pma+pmth)*(pma+pmrc))))
275  sldn=log(max(1.1d0,s/(alp*sqmu*(pma+pmth)*(pma+pmrc))))
276  sum3=cres*log(1d0+((pma+pmrc)/(pma+pmth))**2)*log(slup/sldn)/
277  & (2d0*alp)
278  bxx=ceffd(idd,7)+ceffd(idd,8)/srt+ceffd(idd,9)/s
279  slrr=log(s/(alp*(pma+pmth)*(pma+pmrc)*(pmb+pmth)*(pmb+pmrc)))
280  sum4=cres**2*log(1d0+((pma+pmrc)/(pma+pmth))**2)*
281  & log(1d0+((pmb+pmrc)/(pmb+pmth))**2)/max(0.1d0,2d0*alp*slrr+bxx)
282  sigtmp(i,4)=facdd*xpar(ipr)*max(0d0,sum1+sum2+sum3+sum4)
283 
284 C...Non-diffractive by unitarity.
285  sigtmp(i,5)=sigtmp(i,0)-sigtmp(i,1)-sigtmp(i,2)-sigtmp(i,3)-
286  & sigtmp(i,4)
287  110 CONTINUE
288 
289 C...Put temporary results in output array: only one process.
290  IF(mint(101).EQ.1.AND.mint(102).EQ.1) THEN
291  DO 120 j=0,5
292  sigt(0,0,j)=sigtmp(1,j)
293  120 CONTINUE
294 
295 C...Beam multiple processes.
296 C In principle the power is rho:2.575, phi:2.4 and J/psi:2.55
297  ELSEIF(mint(101).EQ.4.AND.mint(102).EQ.1) THEN
298  DO 140 i=1,4
299  IF(mint(107).EQ.2) THEN
300  vint(317)=(pmhad(i)**2/(pmhad(i)**2+vint(307)))**expo(i)
301  ELSE
302  vint(317)=16d0*parp(15)**2*vint(154)**2/
303  & ((4d0*parp(15)**2+vint(307))*(4d0*vint(154)**2+vint(307)))
304  ENDIF
305  IF(mstp(20).GT.0) THEN
306  vint(317)=vint(317)*(vint(2)/(vint(2)+vint(307)))**mstp(20)
307  ENDIF
308  IF(mint(107).EQ.2) THEN
309  conv=(aem/parp(160+i))*vint(317)
310  ELSEIF(vint(154).GT.parp(15)) THEN
311  conv=(aem/paru(1))*(kchg(i,1)/3d0)**2*parp(18)**2*
312  & (1d0/parp(15)**2-1d0/vint(154)**2)*vint(317)
313  ELSE
314  conv=0d0
315  ENDIF
316  i1=max(1,i-1)
317  DO 130 j=0,5
318  sigt(i,0,j)=conv*sigtmp(i1,j)
319  130 CONTINUE
320  140 CONTINUE
321  DO 150 j=0,5
322  sigt(0,0,j)=sigt(1,0,j)+sigt(2,0,j)+sigt(3,0,j)+sigt(4,0,j)
323  150 CONTINUE
324 
325 C...Target multiple processes.
326  ELSEIF(mint(101).EQ.1.AND.mint(102).EQ.4) THEN
327  DO 170 i=1,4
328  IF(mint(108).EQ.2) THEN
329  vint(317)=(pmhad(i)**2/(pmhad(i)**2+vint(308)))**expo(i)
330  ELSE
331  vint(317)=16d0*parp(15)**2*vint(154)**2/
332  & ((4d0*parp(15)**2+vint(308))*(4d0*vint(154)**2+vint(308)))
333  ENDIF
334  IF(mstp(20).GT.0) THEN
335  vint(317)=vint(317)*(vint(2)/(vint(2)+vint(308)))**mstp(20)
336  ENDIF
337  IF(mint(108).EQ.2) THEN
338  conv=(aem/parp(160+i))*vint(317)
339  ELSEIF(vint(154).GT.parp(15)) THEN
340  conv=(aem/paru(1))*(kchg(i,1)/3d0)**2*parp(18)**2*
341  & (1d0/parp(15)**2-1d0/vint(154)**2)*vint(317)
342  ELSE
343  conv=0d0
344  ENDIF
345  iv=max(1,i-1)
346  DO 160 j=0,5
347  sigt(0,i,j)=conv*sigtmp(iv,j)
348  160 CONTINUE
349  170 CONTINUE
350  DO 180 j=0,5
351  sigt(0,0,j)=sigt(0,1,j)+sigt(0,2,j)+sigt(0,3,j)+sigt(0,4,j)
352  180 CONTINUE
353 
354 C...Both beam and target multiple processes.
355  ELSE
356  IF(mint(107).EQ.2) THEN
357  vint(317)=(pmhad(2)**2/(pmhad(2)**2+vint(307)))**2
358  ELSE
359  vint(317)=16d0*parp(15)**2*vint(154)**2/
360  & ((4d0*parp(15)**2+vint(307))*(4d0*vint(154)**2+vint(307)))
361  ENDIF
362  IF(mint(108).EQ.2) THEN
363  vint(317)=vint(317)*(pmhad(2)**2/(pmhad(2)**2+vint(308)))**2
364  ELSE
365  vint(317)=vint(317)*16d0*parp(15)**2*vint(154)**2/
366  & ((4d0*parp(15)**2+vint(308))*(4d0*vint(154)**2+vint(308)))
367  ENDIF
368  IF(mstp(20).GT.0) THEN
369  vint(317)=vint(317)*(vint(2)/(vint(2)+vint(307)+
370  & vint(308)))**mstp(20)
371  ENDIF
372  DO 210 i1=1,4
373  DO 200 i2=1,4
374  IF(mint(107).EQ.2) THEN
375  conv=(aem/parp(160+i1))*vint(317)
376  ELSEIF(vint(154).GT.parp(15)) THEN
377  conv=(aem/paru(1))*(kchg(i1,1)/3d0)**2*parp(18)**2*
378  & (1d0/parp(15)**2-1d0/vint(154)**2)*vint(317)
379  ELSE
380  conv=0d0
381  ENDIF
382  IF(mint(108).EQ.2) THEN
383  conv=conv*(aem/parp(160+i2))
384  ELSEIF(vint(154).GT.parp(15)) THEN
385  conv=conv*(aem/paru(1))*(kchg(i2,1)/3d0)**2*parp(18)**2*
386  & (1d0/parp(15)**2-1d0/vint(154)**2)
387  ELSE
388  conv=0d0
389  ENDIF
390  IF(i1.LE.2) THEN
391  iv=max(1,i2-1)
392  ELSEIF(i2.LE.2) THEN
393  iv=max(1,i1-1)
394  ELSEIF(i1.EQ.i2) THEN
395  iv=2*i1-2
396  ELSE
397  iv=5
398  ENDIF
399  DO 190 j=0,5
400  jv=j
401  IF(i2.GT.i1.AND.(j.EQ.2.OR.j.EQ.3)) jv=5-j
402  sigt(i1,i2,j)=conv*sigtmp(iv,jv)
403  190 CONTINUE
404  200 CONTINUE
405  210 CONTINUE
406  DO 230 j=0,5
407  DO 220 i=1,4
408  sigt(i,0,j)=sigt(i,1,j)+sigt(i,2,j)+sigt(i,3,j)+sigt(i,4,j)
409  sigt(0,i,j)=sigt(1,i,j)+sigt(2,i,j)+sigt(3,i,j)+sigt(4,i,j)
410  220 CONTINUE
411  sigt(0,0,j)=sigt(1,0,j)+sigt(2,0,j)+sigt(3,0,j)+sigt(4,0,j)
412  230 CONTINUE
413  ENDIF
414 
415 C...Scale up uniformly for Donnachie-Landshoff parametrization.
416  IF(iproc.EQ.21.OR.iproc.EQ.23.OR.iproc.EQ.25) THEN
417  rfac=(xpar(iproc)*seps+ypar(iproc)*seta)/sigt(0,0,0)
418  DO 260 i1=0,n1
419  DO 250 i2=0,n2
420  DO 240 j=0,5
421  sigt(i1,i2,j)=rfac*sigt(i1,i2,j)
422  240 CONTINUE
423  250 CONTINUE
424  260 CONTINUE
425  ENDIF
426 
427  RETURN
428  END