1 gaskelld 1.1 real*8 function semi_dilution(vertex,main)
2
3
4 * Sept. 2003 D. Gaskell
5 * Purpose:
6 * This routine calculates p(e,e'pi+)X semi-iinclusive cross sections from
7 * the CTEQ5M parton distribution functions and a simple parameterization
8 * of the favored and unfavored fragmentation functions.
9 * output:
10 * sigma_eepiX !d3sigma/dEe'dOmegae'Omegapi (microbarn/MeV/sr^2)
11 *
12 * For now, just does PI+ !!!!!!!!!!!!!!!!!!
13 *
14 * October 3 2003 D. Gaskell
15 * Replace simple fragmentation function paramterization with a slightly
16 * more sophisticated treatment from Binnewies et. al.,PRD 52, p.4947 (1995).
17 * This parameterization gives: D_u^(pi+ + pi-) = D_d^(pi+ + pi-) = D+ + D-
18 * The separate favored and unfavored fragmentation functions are given by
19 * the ratio D-/D+ from HERMES data (my fit to P. Geiger's results).
20 * The D-/D+ fit is valid from z=0.25 to 0.9 or so, but behaves pretty well
21 * at lower z (approaches simple form of Field and Feynman (1-z)/(1+z)).
22 gaskelld 1.1 *
23 *
24 * October 4 2003 D. Gaskell
25 * Add PI- functionality
26 *
27 * October 9 2003 D. Gaskell
28 * Add deuterium functionality. Assume just an incoherent sum of parton
29 * distributions from proton and neutron. Assume isospin symmetry so
30 * d_neutron(x) = u_proton(x)
31 * dbar_neutron(x) = ubar_proton(x)
32 * u_neutron(x) = d_proton(x)
33 * ubar_neutron(x) = dbar_proton(x)
34 *
35 *
36 * October 14 - Now adapted to get the dilution factor.
37
38
39
40 implicit none
41 include 'simulate.inc'
42
43 gaskelld 1.1 type(event_main):: main
44 type(event):: vertex
45
46 * PDFs
47 integer iset !which set (1=cteq5m)
48 integer ipart !particle u=1, ubar=-1, d=2, dbar=-2
49 real*8 u,d,ubar,dbar
50 real*8 qu,qd ! u and d quark charges
51 real*8 D_fav, D_unfav, D_sum, R_D !favored,unfavored,sum,ratio of FFs
52 real*8 lambda, Q2zero ! scales for FF param
53 real*8 sv !scaling variable for FF param
54 real*8 N,a1,a2 !parameters for FF param
55
56 real*8 xbj,y,s, Q2gev, Qgev, pt2gev
57 real*8 b ! pt2 parameter for FFs
58
59 real*8 sum_sq, dsigdz, sigsemi, jacobian, sigma_eepiX
60 real*8 sighad, sige
61
62 real*8 sum_sq_prot,sum_sq_neut,dsigdz_prot,dsigdz_neut
63 real*8 sighad_prot,sighad_neut, sige_prot, sige_neut
64 gaskelld 1.1 real*8 sigsemi_prot,sigsemi_neut,dilution
65
66 real*8 N15,He4,Al,Ni,Cu
67
68 real*8 Ctq5Pdf
69 c external Ctq5Pdf
70
71 logical first
72
73
74 parameter (iset=1)
75 parameter (qu=2./3.)
76 parameter (qd=-1./3.)
77
78 parameter (b=3.76) !GeV^-2
79
80 parameter (lambda=0.227) !0.227 GeV for NLO
81 parameter (Q2zero=2.0) !Gev^2 for u,d,s,g
82
83 data first /.TRUE./
84
85 gaskelld 1.1 if(first) then
86 call SetCtq5(iset) ! initialize Cteq5 (we're using cteq5m)
87 first=.FALSE.
88 endif
89
90 s = (2.*vertex%Ein*mp + mp**2)/1.d6 !convert to GeV2
91 xbj = vertex%Q2/2./mp/vertex%nu
92 if(xbj.gt.1.0) then
93 write(6,*) 'XBj is too large!', xbj
94 xbj=1.0
95 endif
96 y = vertex%nu/vertex%Ein
97
98 C DJG convert some stuff to GeV
99
100 Q2gev = vertex%q2/1.d6
101 Qgev = sqrt(Q2gev)
102 pt2gev = vertex%pt2/1.d6
103
104 C Get the PDFs
105 ipart=1
106 gaskelld 1.1 u = Ctq5pdf (ipart , xbj, Qgev)
107
108 ipart=-1
109 ubar = Ctq5pdf (ipart , xbj, Qgev)
110
111 ipart=2
112 d = Ctq5pdf (ipart , xbj, Qgev)
113
114 ipart=-2
115 dbar = Ctq5pdf (ipart , xbj, Qgev)
116
117 C Simple paramaterization from Kretzer et al (EPJC 22 p. 269)
118 C for Q2=2.5.
119
120 c D_fav = 0.689*vertex.zhad**(-1.039)*(1.0-vertex.zhad)**1.241
121 c D_unfav = 0.217*vertex.zhad**(-1.805)*(1.0-vertex.zhad)**2.037
122
123 C
124
125 C Paramaterization from Binneweis et al
126
127 gaskelld 1.1 sv = log( log(Q2gev/lambda**2)/log(Q2zero/lambda**2) )
128
129 C Form of parameterization is D = N z^a1 (1-z)^a2
130
131 if(doing_semipi) then
132
133 N = 1.150 - 1.522*sv + 1.378*sv**2 - 0.527*sv**3
134 a1 = -0.740 - 1.680*sv + 1.546*sv**2 - 0.596*sv**3
135 a2 = 1.430 + 0.543*sv - 0.023*sv**2
136
137 elseif(doing_semika) then
138
139 STOP 'Kaons fragmentation functions not implemented yet'
140
141 endif
142
143 D_sum = N*vertex%zhad**a1*(1.0-vertex%zhad)**a2
144
145 C Ratio of D-/D+ from P. Geiger's thesis (HERMES)
146
147 R_D = (1.0-vertex%zhad)**0.083583/(1.0+vertex%zhad)**1.9838
148 gaskelld 1.1
149 D_fav = D_sum/(1.0+R_D)
150 D_unfav = D_sum/(1.0+1.0/R_D)
151
152 sum_sq = qu**2*(u+ubar) + qd**2*(d+dbar)
153
154 sum_sq_prot = sum_sq
155 sum_sq_neut = qu**2*(d+dbar) + qd**2*(u+ubar)
156
157 if (doing_deutsemi) then
158 sum_sq = sum_sq + qu**2*(d+dbar) + qd**2*(u+ubar)
159 endif
160
161 if(sum_sq.gt.0.) then
162 if(doing_hplus) then
163 dsigdz = (qu**2*u*D_fav + qu**2*ubar*D_unfav +
164 > qd**2*d*D_unfav + qd**2*dbar*D_fav)/sum_sq
165 else
166 dsigdz = (qu**2*u*D_unfav + qu**2*ubar*D_fav +
167 > qd**2*d*D_fav + qd**2*dbar*D_unfav)/sum_sq
168 endif
169 gaskelld 1.1 if(doing_deutsemi) then
170 if(doing_hplus) then
171 dsigdz = dsigdz + (qu**2*d*D_fav + qu**2*dbar*D_unfav +
172 > qd**2*u*D_unfav + qd**2*ubar*D_fav)/sum_sq
173 else
174 dsigdz = dsigdz + (qu**2*d*D_unfav + qu**2*dbar*D_fav +
175 > qd**2*u*D_fav + qd**2*ubar*D_unfav)/sum_sq
176 endif
177 endif
178 else
179 dsigdz = 0.0
180 endif
181
182 if(sum_sq_prot.gt.0.) then
183 if(doing_hplus) then
184 dsigdz_prot = (qu**2*u*D_fav + qu**2*ubar*D_unfav +
185 > qd**2*d*D_unfav + qd**2*dbar*D_fav)/sum_sq_prot
186
187 dsigdz_neut = (qu**2*d*D_fav + qu**2*dbar*D_unfav +
188 > qd**2*u*D_unfav + qd**2*ubar*D_fav)/sum_sq_neut
189 else
190 gaskelld 1.1 dsigdz_prot = (qu**2*u*D_unfav + qu**2*ubar*D_fav +
191 > qd**2*d*D_fav + qd**2*dbar*D_unfav)/sum_sq_prot
192
193 dsigdz_neut = (qu**2*d*D_unfav + qu**2*dbar*D_fav +
194 > qd**2*u*D_fav + qd**2*ubar*D_unfav)/sum_sq_neut
195 endif
196 else
197 dsigdz_prot = 0.0
198 dsigdz_neut = 0.0
199 endif
200
201 sighad = dsigdz*b*exp(-b*pt2gev)/2./pi
202
203 sige = alpha**2*(1.+(1.-y)**2)*(vertex%e%E/1000.)/
204 > (s*xbj*y**2*(mp/1000.)*(vertex%nu/1000.))
205 > * sum_sq
206
207
208 C DJG This dsig/(dOmega_e dE_e dz dpt**2 dPhi_had) in microbarn/GeV**3/sr
209 sigsemi = sige*sighad*(hbarc/1000.)**2*10000.0
210
211 gaskelld 1.1 sighad_prot = dsigdz_prot*b*exp(-b*pt2gev)/2./pi
212
213 sige_prot = alpha**2*(1.+(1.-y)**2)*(vertex%e%E/1000.)/
214 > (s*xbj*y**2*(mp/1000.)*(vertex%nu/1000.))
215 > * sum_sq_prot
216
217
218 C DJG This dsig/(dOmega_e dE_e dz dpt**2 dPhi_had) in microbarn/GeV**3/sr
219 sigsemi_prot = sige_prot*sighad_prot*(hbarc/1000.)**2*10000.0
220
221
222 sighad_neut = dsigdz_neut*b*exp(-b*pt2gev)/2./pi
223
224 sige_neut = alpha**2*(1.+(1.-y)**2)*(vertex%e%E/1000.)/
225 > (s*xbj*y**2*(mp/1000.)*(vertex%nu/1000.))
226 > * sum_sq_neut
227
228
229 C DJG This dsig/(dOmega_e dE_e dz dpt**2 dPhi_had) in microbarn/GeV**3/sr
230 sigsemi_neut = sige_neut*sighad_neut*(hbarc/1000.)**2*10000.0
231
232 gaskelld 1.1
233
234 N15 = 8.0*sigsemi_neut + 7.0*sigsemi_prot
235 He4 = 2.0*sigsemi_neut + 2.0*sigsemi_prot
236 Al = 14.0*sigsemi_neut + 13.0*sigsemi_prot
237 Cu = 35.*sigsemi_neut + 29.0*sigsemi_prot
238 Ni = 31.0*sigsemi_neut + 28.0*sigsemi_prot
239
240
241 if(sigsemi_prot.gt.0.0.and.sigsemi_neut.gt.0.0) then
242 dilution = sigsemi/(sigsemi+0.3333*N15+0.3165*He4+
243 > 0.0144*Al + 0.0031*Cu + 0.0013*Ni)
244 else
245 dilution = 999
246 endif
247
248 C Need to convert to dsig/ (dOmega_e dE_e dE_h dCos(theta) dPhi_had
249 C This is just given by 1/omega * 2*p_h**2*cos(theta)
250
251 c jacobian = 1./(vertex.nu/1000.)*2.*(vertex.p.P/1000.)**2
252 c 1 *cos(vertex.theta_pq)
253 gaskelld 1.1
254
255 c sigma_eepiX = sigsemi*jacobian/1.d6
256
257 c main.davejac = jacobian
258 c ntup.sigcm = sighad
259
260 c peepiX = sigma_eepiX
261 semi_dilution = dilution
262
263 return
264
265 end
|
Return to
semi_dilution.f
CVS log
Up to [HallC] / simc_gfortran