Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pythia_radgen_extras.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pythia_radgen_extras.f
1 C A collection of routines needed by RADGEN
2 
3 !---------------------------------------------------------------
4 ! Calculate the F2 structure function
5 
6  subroutine mkf2 (DQ2,DX,A,Z,DF2,DF1)
7 
8  implicit none
9  include "py6strf.inc"
10  include "mc_set.inc"
11 
12  double precision dq2, dx, df2, df1, df2nf2p
13  DOUBLE PRECISION f2allm, dnp(5), dr, xpq
14  DOUBLE PRECISION gamma2, dnu, w2, pmass2
15  integer a, z
16  dimension xpq(-25:25)
17 
18 C parameters for ratio F2(n)/F2(p)
19 C measured at NMC, (NMC Amaudruz et al. CERN-PPE/91-167))
20  data dnp / 0.976d0, -1.34d0, 1.319d0,
21  & -2.133d0, 1.533d0/
22 
23  pmass2=massp**2
24  w2=pmass2+(dq2*(1/dx-1))
25  dnu=(w2-pmass2+dq2)/(2.*massp)
26  gamma2=dq2/(dnu**2)
27 
28  if ((a.eq.1).and.(z.eq.1)) then
29 *: ALLM:
30  IF(genset_fstruct(1:4).EQ.'ALLM') THEN
31  Call mkr(dq2,dx,dr)
32  df2=f2allm(dx,dq2)
33  df1=(1.d0+gamma2)/(2.d0*dx)/(1.d0+dr)*df2
34  ELSEIF(genset_fstruct(1:4).EQ.'F2PY') THEN
35  call f2pyth(dx,dq2,df1,df2,z)
36  ELSEIF(genset_fstruct(1:4).EQ.'PDF') THEN
37  call pypdfu(2212,dx,dq2,xpq)
38  df2=1d0/9d0*(xpq(1)+xpq(-1)+xpq(3)+xpq(-3))+
39  & 4d0/9d0*(xpq(2)+xpq(-2))
40  df1=(1.d0+gamma2)/(2.d0*dx)*df2
41  ELSE
42 *: error:
43  write(*,*)('invalid parametrisation choice in mkf2')
44  ENDIF
45  ELSEIF(a.EQ.2.and.z.eq.1)THEN
46 *: ALLM:
47  IF(genset_fstruct(1:4).EQ.'ALLM') THEN
48  Call mkr(dq2,dx,dr)
49  df2=f2allm(dx,dq2)
50  df2nf2p=dnp(1)+dx*(dnp(2)+dx*(dnp(3)+dx*(dnp(4)+dx*dnp(5))))
51  df2=df2*0.5*(df2nf2p+1.)
52  df1=(1.d0+gamma2)/(2.d0*dx)/(1.d0+dr)*df2
53  ELSE
54 *: error
55  write(*,*)('MKF2: invalid parametrisation choice FStruct')
56  ENDIF
57 *
58 ! ... neutron = 2*(deuterium_per_nucleon) - proton:
59  ELSEIF(a.EQ.1.and.z.eq.0)THEN
60 *: ALLM:
61  IF(genset_fstruct(1:4).EQ.'ALLM') THEN
62  df2=f2allm(dx,dq2)
63  Call mkr(dq2,dx,dr)
64  df2nf2p=dnp(1)+dx*(dnp(2)+dx*(dnp(3)+dx*(dnp(4)+dx*dnp(5))))
65  df2=df2*df2nf2p
66  df1=(1.d0+gamma2)/(2.d0*dx)/(1.d0+dr)*df2
67  ELSEIF(genset_fstruct(1:4).EQ.'F2PY') THEN
68  call f2pyth(dx,dq2,df1,df2,z)
69  ELSEIF(genset_fstruct(1:4).EQ.'PDF') THEN
70  call pypdfu(2112,dx,dq2,xpq)
71  df2=1d0/9d0*(xpq(1)+xpq(-1)+xpq(3)+xpq(-3))+
72  & 4d0/9d0*(xpq(2)+xpq(-2))
73  df1=(1.d0+gamma2)/(2.d0*dx)*df2
74  ELSE
75 *: error
76  write(*,*)('MKF2: invalid parametrisation choice FStruct')
77  ENDIF
78  ELSEIF((a.gt.1).and.(z.gt.1)) then
79  IF(genset_fstruct(1:4).EQ.'PDF') THEN
80  call pypdfu(2212,dx,dq2,xpq)
81  df2=1d0/9d0*(xpq(1)+xpq(-1)+xpq(3)+xpq(-3))+
82  & 4d0/9d0*(xpq(2)+xpq(-2))
83  df1=(1.d0+gamma2)/(2.d0*dx)*df2
84  ELSE
85 *: error:
86  write(*,*)('invalid parametrisation choice in mkf2')
87  ENDIF
88  ELSE
89 *: error:
90  write(*,*)('MKF2: invalid target type'), a, z
91  ENDIF
92 
93  return
94  end
95 
96 !---------------------------------------------------------------
97 ! Calculate R = sigma_L/sigma_T
98 
99  SUBROUTINE mkr(DQ2,DX,DR)
100  IMPLICIT NONE
101  include "mc_set.inc"
102  include "py6strf.inc"
103 
104  DOUBLE PRECISION dq2, dx
105  DOUBLE PRECISION dr,deltar
106 
107  IF ( genset_r .EQ. '1990' ) THEN
108 * Whitlow et al., Phys.Lett.B 250(1990),193
109  CALL r1990(dq2,dx,dr,deltar)
110  py6r=dr
111  ELSE IF ( genset_r .EQ. '1998' ) THEN
112 * E143, hep-ex/9808028
113  CALL r1998(dq2,dx,dr,deltar)
114  py6r=dr
115  ELSE IF ( genset_r .eq. '0' ) THEN
116 * pure transverse (sigma_L=0)
117  dr = 0.d0
118  py6r=0.d0
119  ELSE
120  write(*,*)( 'MKR: invalid choice for R parametrization' )
121  ENDIF
122 
123  RETURN
124  END
125 
126 C------------------------------------------------------------------
127 
128  SUBROUTINE r1990(DQ2,DX,DR,DELTAR)
129 
130  IMPLICIT NONE
131 
132  DOUBLE PRECISION dq2, dx
133  DOUBLE PRECISION dr, deltar
134 
135  REAL r
136  REAL qq35, xx
137  REAL fac, rlog, q2thr
138  REAL r_a, r_b, r_c
139 C
140 C Data-Definition of R-Calculation, see
141 C L.W.WHITLOW, SLAC-REPORT-357,
142 C PH.D. THESIS, STANFORD UNIVERSITY,
143 C MARCH 1990.
144  REAL ar1990(3), br1990(3), cr1990(3)
145  DATA ar1990 / .06723, .46714, 1.89794 /
146  DATA br1990 / .06347, .57468, -.35342 /
147  DATA cr1990 / .05992, .50885, 2.10807 /
148 
149  deltar = 0.
150 
151  xx=real(dx)
152  IF (dq2.LT.0.35) THEN
153  qq35=0.35
154  ELSE
155  qq35=real(dq2)
156  ENDIF
157 C
158 C *** If Q2 < 0.35 then variable "R" is calculated at the fixed Q2 of 0.35
159 C
160  fac = 1+12.*(qq35/(1.+qq35))*(.125**2/(xx**2+.125**2))
161  rlog = fac/log(qq35/.04)
162  q2thr = 5.*(1.-xx)**5
163 
164  r_a = ar1990(1)*rlog +
165  & ar1990(2)/sqrt(sqrt(qq35**4+ar1990(3)**4))
166  r_b = br1990(1)*rlog +
167  & br1990(2)/qq35 + br1990(3)/(qq35**2+.3**2)
168  r_c = cr1990(1)*rlog +
169  & cr1990(2)/sqrt((qq35-q2thr)**2+cr1990(3)**2)
170  r = (r_a+r_b+r_c)/3.
171 
172  IF (dq2.GE.0.35) THEN
173  dr=dble(r)
174  ELSE
175  dr=dble(r)*dq2/0.35
176  ENDIF
177 
178 c print*,'R:',R
179 
180  END
181 
182 C-----------------------------------------------------------------------
183  SUBROUTINE r1998(DQ2,DX,DR,DELTAR)
184 
185 C new fit to R hep-ex/9808028 E143 Collab.
186 C it is based on the old 3 paramter forms
187 C 0.005<x<0.86, 0.5<Q2<130 GeV2
188 C E143 x-section measurement 0.03<x<0.4
189 C with overall norm uncertainty 2.5%
190 C
191 C U. Stoesslein, October 1998
192 C
193 
194  IMPLICIT NONE
195 
196  DOUBLE PRECISION dq2,dx,dr,deltar
197  DOUBLE PRECISION q2,q2max,fac,rlog,q2thr
198  DOUBLE PRECISION r_a_new,r_a,r_b_new,r_b,r_c
199 
200  DOUBLE PRECISION a(6),b(6),c(6)
201 
202  DATA a/ .0485, 0.5470, 2.0621, -.3804, 0.5090, -.0285/
203  DATA b/ .0481, 0.6114, -.3509, -.4611, 0.7172, -.0317/
204  DATA c/ .0577, 0.4644, 1.8288,12.3708,-43.1043,41.7415/
205 
206  DOUBLE PRECISION dr1998
207  EXTERNAL dr1998
208 
209 * use R(0.35) if q2 is below 0.35
210  q2=dq2
211  q2max=0.35
212  IF(q2.LT.q2max) q2=q2max
213 
214  fac = 1+12.*(q2/(1.+q2))*(.125**2/(dx**2+.125**2))
215  rlog = fac/log(q2/.04)
216  q2thr = c(4)*dx+c(5)*dx*dx+c(6)*dx*dx*dx
217 
218 * new additional terms
219  r_a_new = (1.+a(4)*dx+a(5)*dx*dx)*dx**(a(6))
220  r_a = a(1)*rlog + a(2)/sqrt(sqrt(q2**4+a(3)**4))*r_a_new
221  r_b_new = (1.+b(4)*dx+b(5)*dx*dx)*dx**(b(6))
222  r_b = b(1)*rlog + (b(2)/q2 + b(3)/(q2**2+0.3**2))*r_b_new
223  r_c = c(1)*rlog + c(2)/sqrt((q2-q2thr)**2+c(3)**2)
224  dr = (r_a+r_b+r_c)/3.
225 
226 * straight line fit extrapolation to R(Q2=0)=0
227  if (dq2.lt.q2max) dr = dr*dq2/q2max
228 
229 * I assume the fit uncertainty only for measured Q2 range
230  if (q2 .GT. 0.5) then
231  deltar = dr1998(dx,q2)
232  else
233  deltar=dr
234  endif
235 
236  RETURN
237  END
238 
239 C--------------------------------------------------------------------
240  DOUBLE PRECISION FUNCTION dr1998(DX,DQ2)
241 
242 * Parameterize uncertainty in R1998
243 * associated to the fitting procedure only
244 
245  IMPLICIT NONE
246  DOUBLE PRECISION dx,dq2
247 
248  dr1998 = 0.0078-0.013*dx+(0.070-0.39*dx+0.7*dx*dx)/(1.7+dq2)
249 
250  RETURN
251  END
252 
253 !---------------------------------------------------------------
254 ! Calculate the asymmetries A1 and A2. Routine is empty because Pythia
255 ! is unpolarised, but radgen expects it
256 
257  subroutine mkasym (dQ2, dX, A, Z, dA1, dA2)
258 
259  implicit none
260 
261  double precision dq2, dx, da1, da2
262  integer a, z
263 
264  da1 = 0.d0
265  da2 = 0.d0
266 
267  return
268  end
269 
270 !---------------------------------------------------------------
271 ! Calculate the dilution factor
272 
273  double precision function fdilut (dQ2, dx, A)
274 
275  implicit none
276 
277  DOUBLE PRECISION dq2, dx, df
278  DOUBLE PRECISION dnp, dfn, dfp
279  DOUBLE PRECISION dz, df2nf2p
280  INTEGER a
281  dimension dnp(7)
282 *
283 C ... fit to NMC F2n/F2p data (86/87+89 T1,T14)
284  data dnp/
285  + 0.67225d+00,
286  + 0.16254d+01,
287  + -0.15436d+00,
288  + 0.48301d-01,
289  + 0.41979d+00,
290  + 0.47331d-01,
291  + -0.17816d+00/
292 
293 ! Definitions of 'intrinsic' dilutions for neutron and proton (GNOME confused)
294 
295  dfn=1.
296  dfp=1.
297 
298 ! Only 3He has a dilution, and we determine it as F2n/(2*F2p + F2n)
299 
300  if (a.ne.3) then
301  df = 1.
302  else
303  dz = 1./2.*dlog(1.+dexp(2.0-1000.*dx))
304  df2nf2p = dnp(1)*(1.0-dx)**dnp(2)+dnp(3)*dx**dnp(4)
305  1 +(dnp(5)+dnp(6)*dz+dnp(7)*dz**2)
306  df = dfn*(1./((2./df2nf2p)+1))
307  endif
308 
309  fdilut = df
310  return
311 
312  END
313 
314 !---------------------------------------------------------------
315 ! Function to calculate F2 from ALLM parametrisation
316 
317  double precision FUNCTION f2allm(x,q2)
318 
319  implicit double precision (a-h,o-z)
320 
321  REAL m02,m12,lam2,m22
322  common/allm/sp,ap,bp,sr,ar,br,s,xp,xr,f2p,f2r
323 C POMERON
324  parameter(
325  , s11 = 0.28067, s12 = 0.22291, s13 = 2.1979,
326  , a11 = -0.0808 , a12 = -0.44812, a13 = 1.1709,
327  , b11 = 0.60243**2, b12 = 1.3754**2, b13 = 1.8439,
328  , m12 = 49.457 )
329 
330 C REGGEON
331  parameter(
332  , s21 = 0.80107, s22 = 0.97307, s23 = 3.4942,
333  , a21 = 0.58400, a22 = 0.37888, a23 = 2.6063,
334  , b21 = 0.10711**2, b22 = 1.9386**2, b23 = 0.49338,
335  , m22 = 0.15052 )
336 C
337  parameter( m02=0.31985, lam2=0.065270, q02=0.46017 +lam2 )
338  parameter( alfa=112.2, xmp2=0.8802)
339 C
340  w2=q2*(1./x -1.)+xmp2
341  w=dsqrt(w2)
342 C
343  IF(q2.EQ.0.) THEN
344  s=0.
345  z=1.
346 C
347 C POMERON
348 C
349  xp=1./(1.+(w2-xmp2)/(q2+m12))
350  ap=a11
351  bp=b11
352  sp=s11
353  f2p=sp*xp**ap
354 C
355 C REGGEON
356 C
357  xr=1./(1.+(w2-xmp2)/(q2+m22))
358  ar=a21
359  br=b21
360  sr=s21
361  f2r=sr*xr**ar
362 C
363  ELSE
364  s=dlog(dlog((q2+q02)/lam2)/dlog(q02/lam2))
365  z=1.-x
366 C
367 C POMERON
368 C
369  xp=1./(1.+(w2-xmp2)/(q2+m12))
370  ap=a11+(a11-a12)*(1./(1.+s**a13)-1.)
371  bp=b11+b12*s**b13
372  sp=s11+(s11-s12)*(1./(1.+s**s13)-1.)
373  f2p=sp*xp**ap*z**bp
374 C
375 C REGGEON
376 C
377  xr=1./(1.+(w2-xmp2)/(q2+m22))
378  ar=a21+a22*s**a23
379  br=b21+b22*s**b23
380  sr=s21+s22*s**s23
381  f2r=sr*xr**ar*z**br
382 
383 C
384  ENDIF
385 
386 c CIN=ALFA/(Q2+M02)*(1.+4.*XMP2*Q2/(Q2+W2-XMP2)**2)/Z
387 c SIGal=CIN*(F2P+F2R)
388 c f2allm=sigal/alfa*(q2**2*(1.-x))/(q2+4.*xmp2*x**2)
389  f2allm = q2/(q2+m02)*(f2p+f2r)
390 
391  RETURN
392  END
393 C--------------------------------------------------------------------
394 
395 C.....Total inclusive Pythia cross section according to the model in
396 C.....C.Friberg, T.Sjoestrand, J. High Energy Phys. JHEP 0009, 010 (2000)
397 C.....(neglecting GVMD cross sections)
398 C.....P.Liebing, 08/16/2003
399 
400  SUBROUTINE f2pyth(x,q2,f1,f2,z)
401 
402  IMPLICIT NONE
403  common/pyint1/mint(400),vint(400)
404  INTEGER mint
405  DOUBLE PRECISION vint
406  SAVE/pyint1/
407  common/pypars/mstp(200),parp(200),msti(200),pari(200)
408  INTEGER mstp,msti
409  DOUBLE PRECISION parp,pari
410  SAVE/pypars/
411  DOUBLE PRECISION f2dis,f1dis,sdis,svmd1,svmd2,sigvm,x,q2,dipol
412  DOUBLE PRECISION mrho2,alpha,eps,eta,pmass,pmass2,gevmb,rvmd
413  DOUBLE PRECISION w2,gamma2,nu,convf2,convf1,conv,sigh,pi
414  DOUBLE PRECISION f1,f2,df2allm,df2allml,df1allm,df1allml,r,rl
415  DOUBLE PRECISION w2l,gamma2l,nul,xl,sw2
416  INTEGER z
417 C...Local arrays.
418  DOUBLE PRECISION xpq,xpar,ypar
419  INTEGER i
420  dimension xpq(-25:25),xpar(4),ypar(4)
421 
422 C...X and Y parameters of sigmatot = X * s**epsilon + Y * s**(-eta).
423  DATA xpar/2*13.63d0,10.01d0,0.970d0/
424  DATA ypar/2*31.79d0,-1.51d0,-0.146d0/
425  DATA eps/0.0808d0/,eta/-0.4525d0/
426  DATA mrho2/0.591822d0/,alpha/7.297352533d-3/,pmass/0.93827d0/,
427  & gevmb/0.389379292d0/,pi/3.14159265358979324d0/
428 
429  EXTERNAL pypdfu
430  f2dis=0d0
431  f1dis=0d0
432  sdis=0d0
433  svmd1=0d0
434  svmd2=0d0
435  sigvm=0d0
436  rvmd=0d0
437  dipol=0d0
438  w2=0d0
439  gamma2=0d0
440  nu=0d0
441  convf2=0d0
442  convf1=0d0
443  conv=0d0
444  sigh=0d0
445  f1=0d0
446  f2=0d0
447  df1allm=0d0
448  df2allm=0d0
449  df1allml=0d0
450  df2allml=0d0
451  r=0d0
452  rl=0d0
453  sw2=0d0
454 c....Kinematic factors and constants
455  pmass2=pmass**2
456  nul=-1d0
457  gamma2l=0d0
458  w2l=-1d0
459  xl=0d0
460 
461  if ((x.gt.0d0).and.(x.le.1d0)) then
462  w2=pmass2+(q2*(1./x-1.))
463  else
464  f1=0d0
465  f2=0d0
466  return
467  endif
468 
469  if (w2.lt.4d0) then
470  w2l=w2
471  w2=4d0
472  xl=x
473  nul=(w2l-pmass2+q2)/(2d0*pmass)
474  endif
475 
476  nu=(w2-pmass2+q2)/(2d0*pmass)
477 
478  if (nu.gt.0d0) then
479  gamma2=q2/(nu**2)
480  if (nul.gt.0d0) then
481  gamma2l=q2/(nul**2)
482  endif
483  else
484  f1=0d0
485  f2=0d0
486  return
487  endif
488 
489 c....For W2<4, we don't trust the PYTHIA F2, so we calculate
490 C....F2-ALLM(W2,Q2). The real kinematics have and "l" at the end,
491 C....the kinematics without the "l" are the ones we get by setting W2 to 4
492 c....Output: f2allm=F2-ALLM(W2=4,Q2),f2allml=F2-ALLM(W2=w2l,Q2)
493  if (w2l.gt.0d0) then
494  sw2=((w2l-pmass2)/(4d0-pmass2))**10
495  endif
496 c....This factor is needed to convert the Pythia virtual gamma cross
497 c....section for VMD to the same level as F2
498 c....The kinematic factors making the (ep) cross section out of F2 are
499 c....provided by RADGEN
500  conv=q2*(1d0-x)/(4d0*pi**2*alpha)/gevmb
501 c....Pythia PDF call, sum PDFs to F2
502  if (z.eq.1) then
503  call pypdfu(2212,x,q2,xpq)
504  elseif (z.eq.0) then
505  call pypdfu(2112,x,q2,xpq)
506  endif
507  f2dis=1d0/9d0*(xpq(1)+xpq(-1)+xpq(3)+xpq(-3))+
508  & 4d0/9d0*(xpq(2)+xpq(-2))
509 c....Suppression factor for DIS
510  if (mstp(19).eq.0) then
511  sdis=1.
512  else
513  sdis=q2/(q2+mrho2)
514  if (mstp(19).gt.1) then
515  sdis=sdis**2
516  endif
517  endif
518 C....Sum of Hadronic (Vector Meson) cross sections * Photon couplings
519 C....const.
520  sigh=0.
521  do 10 i=1,4
522  sigh=sigh+alpha/parp(160+i)*(xpar(i)*w2**eps+ypar(i)*w2**eta)
523  10 continue
524 C....W2/Q2 suppression of VMD and (1+epsilon R_VMD)
525  svmd1=(w2/(w2+q2))**mstp(20)
526  if (mstp(20).eq.0) then
527  dipol=2.575d0
528  else
529  dipol=2d0
530  endif
531  if (mstp(17).eq.6) then
532  rvmd=parp(165)*(q2/mrho2)**parp(166)
533  else
534 c ...Attention: This is only good for MSTP(17)=4, i.e., the Pythia
535 c ...default
536  rvmd=parp(165)*(4.*mrho2*q2)/(mrho2+q2)**2
537  endif
538 C .... Dipole factor for VMD
539  svmd2=(mrho2/(mrho2+q2))**dipol
540 C.....virtual photon xsec for VMD
541  sigvm=svmd1*svmd2*sigh
542  convf2=(1d0+rvmd)/(1d0+gamma2)
543  convf1=1d0/(2d0*x)
544 c.....Total "F2"
545  f2=sdis*f2dis+conv*convf2*sigvm
546  f1dis=(1.d0+gamma2)/(2.d0*x)*f2dis
547  f1=sdis*f1dis+conv*convf1*sigvm
548  if (w2l.gt.0d0) then
549 C.....Here we scale F2-ALLM(W2=w2l,Q2) by the factor
550 C.....F2-PYTH(W2=4,Q2)/F2-ALLM(W2=4,Q2) (normalize ALLM to PYTHIA model)
551  f2=sw2*f2
552  f1=sw2*f1
553  endif
554  RETURN
555  END