EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pepsi_radgen_extras.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pepsi_radgen_extras.f
1 ! 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 "leptou.inc"
10  include "mc_set.inc"
11 
12  double precision dq2, dx, df2, df1
13  DOUBLE PRECISION f2allm, dnp(5), df2nf2p, dratio, dr
14  DOUBLE PRECISION disoc, dhe3d(6), dhe4d(6)
15  DOUBLE PRECISION gamma2, dnu, dw2, pmass2
16  integer a, z, iflavour
17  real xpq(-6:6), xdpq(-6:6)
18 ! ... charge of quark flavours; 2nd index is 1:proton, 2:neutron
19 ! ... (effective isospin rotation for neutron)
20  real qflavour(6,2)
21  save qflavour
22  data qflavour/1.,2.,1.,2.,1.,2.,
23  + 2.,1.,1.,2.,1.,2./
24 c parameters for ratio f2(n)/f2(p)
25 c measured at nmc, (nmc amaudruz et al. cern-ppe/91-167))
26  data dnp / 0.976d0, -1.34d0, 1.319d0,
27  & -2.133d0, 1.533d0/
28 c fit to latest result from taeksu(he3d2_shin280498.fitpar)
29 c to be corrected for non isoscalarity !
30  DATA dhe3d / 0.98733d0, 0.50409d0, -0.22521d0,
31  & -0.66976d2, -0.62318d0, 0.12104d1/
32 
33 c parameters for ratio f2(4he)/f2(d), not corrected for isoscalarity
34 c fit to nmc+slac data (hed2_03.fitpar) - abr
35  DATA dhe4d / 0.1051d1, -0.1896d0, -0.1026d0,
36  & -0.1704d2, 0.3086d0 , 0.8509d1/
37 
38  pmass2=massp**2
39 
40  IF(genset_fstruct(1:4).EQ.'ALLM') THEN
41  if ((a.eq.1).and.(z.eq.1)) then
42  df2=f2allm(dx,dq2)
43  endif
44 
45  if ((a.eq.2).and.(z.eq.1)) then
46  df2=f2allm(dx,dq2)
47  df2nf2p=dnp(1)+dx*(dnp(2)+dx*(dnp(3)+dx*(dnp(4)+dx*dnp(5))))
48  df2=df2*0.5*(df2nf2p+1.)
49  endif
50 
51  if ((a.eq.3).and.(z.eq.2)) then
52 c...f2-deut(from f2-proton(allm) and f2n/f2p)*(he-3/f2-deut.)-ratio:
53  df2=f2allm(dx,dq2)
54  df2nf2p=dnp(1)+dx*(dnp(2)+dx*(dnp(3)+dx*(dnp(4)+dx*dnp(5))))
55  df2=df2*0.5*(df2nf2p+1.)
56  dratio = dhe3d(1) + dhe3d(2)*dx + dhe3d(3)*exp(dhe3d(4)*dx)
57  & + dhe3d(5)*dx**dhe3d(6)
58 c undo isoscalarity correction
59  disoc=((a-z)*df2nf2p+z)/(0.5*a*(1+df2nf2p))
60  df2=df2*dratio*disoc
61  endif
62 
63  if ((a.eq.4).and.(z.eq.2)) then
64 c...f2-deut(from f2-proton(allm) and f2n/f2p)*(f2-4he/f2-deut)-ratio
65  df2=f2allm(dx,dq2)
66  df2nf2p=dnp(1)+dx*(dnp(2)+dx*(dnp(3)+dx*(dnp(4)+dx*dnp(5))))
67  df2=df2*0.5*(df2nf2p+1.)
68  dratio = dhe4d(1) + dhe4d(2)*dx + dhe4d(3)*exp(dhe4d(4)*dx)
69  & + dhe4d(5)*dx**dhe4d(6)
70  df2=df2*dratio
71  endif
72  ENDIF
73 
74  IF (genset_fstruct(1:3).EQ.'PDF') THEN
75  call parton(real(x),real(q2),xpq,xdpq)
76  df2 = 0.
77 c! ... for the moment assume sum over pdfs is F1 this neglects R
78  do iflavour = 1, 6
79  df2 = df2 + (z*(qflavour(iflavour,1)**2)/9.
80  + + (a-z)*(qflavour(iflavour,2)**2)/9. )
81  + *((xpq(iflavour)/dx)+(xpq(-iflavour)/dx))
82  enddo
83  ENDIF
84 
85  Call mkr(dq2,dx,dr)
86  dw2=pmass2+dq2*(1-dx)/dx
87  dnu=(dw2-pmass2+dq2)/(2.*massp)
88  gamma2=dq2/(dnu**2)
89  df1=(1.d0+gamma2)/(2.d0*dx)/(1.d0+dr)*df2
90 
91  return
92  end
93 
94 !---------------------------------------------------------------
95 ! Calculate R = sigma_L/sigma_T
96 
97  SUBROUTINE mkr(DQ2,DX,DR)
98  IMPLICIT NONE
99  include "mc_set.inc"
100 
101  DOUBLE PRECISION dq2, dx
102  DOUBLE PRECISION dr,deltar
103 
104  IF ( genset_r .EQ. '1990' ) THEN
105 * whitlow et al., phys.lett.b 250(1990),193
106  CALL r1990(dq2,dx,dr,deltar)
107  ELSE IF ( genset_r .EQ. '1998' ) THEN
108 * e143, hep-ex/9808028
109  CALL r1998(dq2,dx,dr,deltar)
110  ELSE IF ( genset_r .eq. '0' ) THEN
111 * pure transverse(sigma_l=0)
112  dr = 0.d0
113  ELSE
114  write(*,*)( 'MKR: invalid choice for R parametrization' )
115  ENDIF
116 
117  RETURN
118  END
119 
120 c------------------------------------------------------------------
121 
122  SUBROUTINE r1990(DQ2,DX,DR,DELTAR)
123 
124  IMPLICIT NONE
125 
126  DOUBLE PRECISION dq2, dx
127  DOUBLE PRECISION dr, deltar
128 
129  REAL r
130  REAL qq35, xx
131  REAL fac, rlog, q2thr
132  REAL r_a, r_b, r_c
133 c
134 c data-definition of r-calculation, see
135 c l.w.whitlow, slac-report-357,
136 c ph.d. thesis, stanford university,
137 c march 1990.
138  REAL ar1990(3), br1990(3), cr1990(3)
139  DATA ar1990 / .06723, .46714, 1.89794 /
140  DATA br1990 / .06347, .57468, -.35342 /
141  DATA cr1990 / .05992, .50885, 2.10807 /
142 
143  deltar = 0.
144 
145  xx=real(dx)
146  IF (dq2.LT.0.35) THEN
147  qq35=0.35
148  ELSE
149  qq35=real(dq2)
150  ENDIF
151 c
152 c *** If q2 < 0.35 then variable "R" is calculated at the fixed q2 of 0.35
153 c
154  fac = 1+12.*(qq35/(1.+qq35))*(.125**2/(xx**2+.125**2))
155  rlog = fac/log(qq35/.04)
156  q2thr = 5.*(1.-xx)**5
157 
158  r_a = ar1990(1)*rlog +
159  & ar1990(2)/sqrt(sqrt(qq35**4+ar1990(3)**4))
160  r_b = br1990(1)*rlog +
161  & br1990(2)/qq35 + br1990(3)/(qq35**2+.3**2)
162  r_c = cr1990(1)*rlog +
163  & cr1990(2)/sqrt((qq35-q2thr)**2+cr1990(3)**2)
164  r = (r_a+r_b+r_c)/3.
165 
166  IF (dq2.GE.0.35) THEN
167  dr=dble(r)
168  ELSE
169  dr=dble(r)*dq2/0.35
170  ENDIF
171 
172 c print*,'R:',r
173 
174  END
175 
176 c-----------------------------------------------------------------------
177  SUBROUTINE r1998(DQ2,DX,DR,DELTAR)
178 
179 c new fit to r hep-ex/9808028 e143 collab.
180 c it is based on the old 3 paramter forms
181 c 0.005<x<0.86, 0.5<q2<130 gev2
182 c e143 x-section measurement 0.03<x<0.4
183 c with overall norm uncertainty 2.5%
184 c
185 c u. stoesslein, october 1998
186 c
187 
188  IMPLICIT NONE
189 
190  DOUBLE PRECISION dq2,dx,dr,deltar
191  DOUBLE PRECISION q2,q2max,fac,rlog,q2thr
192  DOUBLE PRECISION r_a_new,r_a,r_b_new,r_b,r_c
193 
194  DOUBLE PRECISION a(6),b(6),c(6)
195 
196  DATA a/ .0485, 0.5470, 2.0621, -.3804, 0.5090, -.0285/
197  DATA b/ .0481, 0.6114, -.3509, -.4611, 0.7172, -.0317/
198  DATA c/ .0577, 0.4644, 1.8288,12.3708,-43.1043,41.7415/
199 
200  DOUBLE PRECISION dr1998
201  EXTERNAL dr1998
202 
203 * use r(0.35) if q2 is below 0.35
204  q2=dq2
205  q2max=0.35
206  IF(q2.LT.q2max) q2=q2max
207 
208  fac = 1+12.*(q2/(1.+q2))*(.125**2/(dx**2+.125**2))
209  rlog = fac/log(q2/.04)
210  q2thr = c(4)*dx+c(5)*dx*dx+c(6)*dx*dx*dx
211 
212 * new additional terms
213  r_a_new = (1.+a(4)*dx+a(5)*dx*dx)*dx**(a(6))
214  r_a = a(1)*rlog + a(2)/sqrt(sqrt(q2**4+a(3)**4))*r_a_new
215  r_b_new = (1.+b(4)*dx+b(5)*dx*dx)*dx**(b(6))
216  r_b = b(1)*rlog + (b(2)/q2 + b(3)/(q2**2+0.3**2))*r_b_new
217  r_c = c(1)*rlog + c(2)/sqrt((q2-q2thr)**2+c(3)**2)
218  dr = (r_a+r_b+r_c)/3.
219 
220 * straight line fit extrapolation to r(q2=0)=0
221  if (dq2.lt.q2max) dr = dr*dq2/q2max
222 
223 * i assume the fit uncertainty only for measured q2 range
224  if (q2 .GT. 0.5) then
225  deltar = dr1998(dx,q2)
226  else
227  deltar=dr
228  endif
229 
230  RETURN
231  END
232 
233 c--------------------------------------------------------------------
234  DOUBLE PRECISION FUNCTION dr1998(DX,DQ2)
235 
236 * parameterize uncertainty in r1998
237 * associated to the fitting procedure only
238 
239  IMPLICIT NONE
240  DOUBLE PRECISION dx,dq2
241 
242  dr1998 = 0.0078-0.013*dx+(0.070-0.39*dx+0.7*dx*dx)/(1.7+dq2)
243 
244  RETURN
245  END
246 
247 !---------------------------------------------------------------
248 ! Calculate the asymmetries A1 and A2
249 
250  subroutine mkasym (dQ2, dX, A, Z, dA1, dA2)
251 
252  implicit none
253  include "leptou.inc"
254 
255  double precision dq2, dx, da1, da2, df1, df2
256  double precision dppar, ddpar, massp, pmass2
257  integer a, z, iflavour
258  real xpq(-6:6), xdpq(-6:6), g1
259 ! ... charge of quark flavours; 2nd index is 1:proton, 2:neutron
260 ! ... (effective isospin rotation for neutron)
261  real qflavour(6,2)
262  save qflavour
263  data qflavour/1.,2.,1.,2.,1.,2.,
264  + 2.,1.,1.,2.,1.,2./
265 
266  massp=0.938272013
267  pmass2=massp**2
268  dppar=0.53
269  ddpar=0.22
270 
271  if (lst(15).le.100) then
272  da1 = 0.d0
273  da2 = 0.d0
274  return
275  endif
276 
277  call mkf2(dq2, dx, a, z, df2, df1)
278 
279  call parton(real(x),real(q2),xpq,xdpq)
280  g1 = 0.
281  do iflavour = 1, 6
282  g1 = g1 + (z*(qflavour(iflavour,1)**2)/9.
283  + + (a-z)*(qflavour(iflavour,2)**2)/9. )
284  + *((xdpq(iflavour)/dx)+(xdpq(-iflavour)/dx))
285  enddo
286 
287  da1 = dble(g1)/df1
288 
289 c if (1.le.dq2.and.dq2.le.1.5) then
290 c write(*,*)'In mkasym', dq2, dx, da1, g1, f1
291 c endif
292 
293  IF(a.EQ.1.and.z.eq.1) THEN
294  da2=dppar*massp*dx/(sqrt(dq2))
295  ELSEIF(a.EQ.2.and.z.eq.1) THEN
296  da2=ddpar*massp*dx/(sqrt(dq2))
297  ELSEIF(a.EQ.1.and.z.eq.0) THEN
298  da2=(ddpar-dppar)*massp*dx/(sqrt(dq2))
299  ELSE
300  da2 = 0.d0
301  ENDIF
302 
303  return
304  end
305 
306 !---------------------------------------------------------------
307 ! Calculate the dilution factor
308 
309  subroutine fdilut (dQ2, dx, A, DF)
310 
311  implicit none
312 
313  DOUBLE PRECISION dq2, dx, df
314  DOUBLE PRECISION dnp, dfn, dfp
315  DOUBLE PRECISION dz, df2nf2p
316  INTEGER a
317  dimension dnp(7)
318 *
319 c ... fit to nmc f2n/f2p data (86/87+89 t1,t14)
320  data dnp/
321  + 0.67225d+00,
322  + 0.16254d+01,
323  + -0.15436d+00,
324  + 0.48301d-01,
325  + 0.41979d+00,
326  + 0.47331d-01,
327  + -0.17816d+00/
328 
329 ! Definitions of 'intrinsic' dilutions for neutron and proton (GNOME confused)
330 
331  dfn=1.
332  dfp=1.
333 
334 ! Only 3He has a dilution, and we determine it as F2n/(2*F2p + F2n)
335 
336  if (a.ne.3) then
337  df = 1.
338  else
339  dz = 1./2.*dlog(1.+dexp(2.0-1000.*dx))
340  df2nf2p = dnp(1)*(1.0-dx)**dnp(2)+dnp(3)*dx**dnp(4)
341  1 +(dnp(5)+dnp(6)*dz+dnp(7)*dz**2)
342  df = dfn*(1./((2./df2nf2p)+1))
343  endif
344 
345  return
346 
347  END
348 
349 !---------------------------------------------------------------
350 ! Function to calculate F2 from ALLM parametrisation
351 
352  double precision FUNCTION f2allm(x,q2)
353 
354  implicit double precision (a-h,o-z)
355 
356  REAL m02,m12,lam2,m22
357  common/allm/sp,ap,bp,sr,ar,br,s,xp,xr,f2p,f2r
358 c pomeron
359  parameter(
360  , s11 = 0.28067, s12 = 0.22291, s13 = 2.1979,
361  , a11 = -0.0808 , a12 = -0.44812, a13 = 1.1709,
362  , b11 = 0.60243**2, b12 = 1.3754**2, b13 = 1.8439,
363  , m12 = 49.457 )
364 
365 c reggeon
366  parameter(
367  , s21 = 0.80107, s22 = 0.97307, s23 = 3.4942,
368  , a21 = 0.58400, a22 = 0.37888, a23 = 2.6063,
369  , b21 = 0.10711**2, b22 = 1.9386**2, b23 = 0.49338,
370  , m22 = 0.15052 )
371 c
372  parameter( m02=0.31985, lam2=0.065270, q02=0.46017 +lam2 )
373  parameter( alfa=112.2, xmp2=0.8802)
374 c
375  w2=q2*(1./x -1.)+xmp2
376  w=dsqrt(w2)
377 c
378  IF(q2.EQ.0.) THEN
379  s=0.
380  z=1.
381 c
382 c pomeron
383 c
384  xp=1./(1.+(w2-xmp2)/(q2+m12))
385  ap=a11
386  bp=b11
387  sp=s11
388  f2p=sp*xp**ap
389 c
390 c reggeon
391 c
392  xr=1./(1.+(w2-xmp2)/(q2+m22))
393  ar=a21
394  br=b21
395  sr=s21
396  f2r=sr*xr**ar
397 c
398  ELSE
399  s=dlog(dlog((q2+q02)/lam2)/dlog(q02/lam2))
400  z=1.-x
401 c
402 c pomeron
403 c
404  xp=1./(1.+(w2-xmp2)/(q2+m12))
405  ap=a11+(a11-a12)*(1./(1.+s**a13)-1.)
406  bp=b11+b12*s**b13
407  sp=s11+(s11-s12)*(1./(1.+s**s13)-1.)
408  f2p=sp*xp**ap*z**bp
409 c
410 c reggeon
411 c
412  xr=1./(1.+(w2-xmp2)/(q2+m22))
413  ar=a21+a22*s**a23
414  br=b21+b22*s**b23
415  sr=s21+s22*s**s23
416  f2r=sr*xr**ar*z**br
417 
418 c
419  ENDIF
420 
421 c cin=alfa/(q2+m02)*(1.+4.*xmp2*q2/(q2+w2-xmp2)**2)/z
422 c sigal=cin*(f2p+f2r)
423 c f2allm=sigal/alfa*(q2**2*(1.-x))/(q2+4.*xmp2*x**2)
424  f2allm = q2/(q2+m02)*(f2p+f2r)
425 
426  RETURN
427  END