1 gaskelld 1.1 real*8 function peepiX(vertex,vertex0,main,survivalprob,doing_cent)
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 *
13 * October 3 2003 D. Gaskell
14 * Replace simple fragmentation function paramterization with a slightly
15 * more sophisticated treatment from Binnewies et. al.,PRD 52, p.4947 (1995).
16 * This parameterization gives: D_u^(pi+ + pi-) = D_d^(pi+ + pi-) = D+ + D-
17 * The separate favored and unfavored fragmentation functions are given by
18 * the ratio D-/D+ from HERMES data (my fit to P. Geiger's results).
19 * The D-/D+ fit is valid from z=0.25 to 0.9 or so, but behaves pretty well
20 * at lower z (approaches simple form of Field and Feynman (1-z)/(1+z)).
21 *
22 gaskelld 1.1 *
23 * October 4 2003 D. Gaskell
24 * Add PI- functionality
25 *
26 * October 9 2003 D. Gaskell
27 * Add deuterium functionality. Assume just an incoherent sum of parton
28 * distributions from proton and neutron. Assume isospin symmetry so
29 * d_neutron(x) = u_proton(x)
30 * dbar_neutron(x) = ubar_proton(x)
31 * u_neutron(x) = d_proton(x)
32 * ubar_neutron(x) = dbar_proton(x)
33 *
34 *
35 * March 2004 D. Gaskell
36 * Add calculation of "central" cross section. This is convenient for
37 * bin centering.
38 *
39 * April 8, 2004 D. Gaskell
40 * Add conversion for jacobian so in ntuple: siglab/jacobian = dsig/dE dOe dz dPt2 dPhi
41 * (units are MeV everywhere)
42 * Add strange quark contributions and kaon fragmentation
43 gaskelld 1.1 * April 15, 2004 D. Gaskell
44 * Put in Fermi motion effects. Note that for now, if you use Fermi motion
45 * the "central cross section" calculation stuff probably doesn't make much sense.
46
47
48
49 implicit none
50 include 'simulate.inc'
51
52 record /event_main/ main
53 record /event/ vertex, vertex0
54
55 * PDFs
56 integer iset !which set (1=cteq5m)
57 integer ipart !particle u=1, ubar=-1, d=2, dbar=-2, s=3, sbar=-3
58 real*8 u,d,ubar,dbar,s,sbar
59 real*8 qu,qd,qs ! u, d, and s quark charges
60 real*8 D_fav, D_unfav, D_sum, R_D !favored,unfavored,sum,ratio of FFs
61 real*8 D_sum_s, D_s ! strange frag. functions
62 real*8 lambda, Q2zero ! scales for FF param
63 real*8 sv !scaling variable for FF param
64 gaskelld 1.1 real*8 N,a1,a2 !parameters for FF param
65 real*8 Ns, a1s, a2s !parameters for strange FF param
66
67 C Some local kinematic variables
68 real*8 xbj,sx, Q2gev, Qgev, pt2gev !unitless or GeV
69 real*8 b ! pt2 parameter for FFs
70
71 real*8 nu,qx,qy,qz,mtar,Q2,Eb,Eprime,Epx,Epy,Epz !MeV
72 real*8 pt2,zhad,Ehad,phad,mhad ! all in MeV
73 real*8 cthpq
74
75 real*8 kcent,klo,khi
76 integer i
77
78 real*8 sum_sq, dsigdz, sigsemi, jacobian, fac, sigma_eepiX
79 real*8 sighad, sige
80
81 real*8 F1,F2,W1,W2,sin2th2,cos2th2,W2coeff
82 real*8 Ctq5Pdf
83 c external Ctq5Pdf
84
85 gaskelld 1.1 C Variables for kaon decay stuff
86 real*8 survivalprob ! this will be included in main.weight
87 real*8 zaero,pathlen,p_kaon,betak,gammak
88
89 logical first
90 logical doing_cent,first_cent! flag for "central" cross section calc.
91
92 parameter (iset=1)
93 parameter (qu=2./3.)
94 parameter (qd=-1./3.)
95 parameter (qs=-1./3.)
96
97 c parameter (b=3.76) !GeV^-2
98 c parameter (b=4.68) !GeV^-2
99
100 parameter (lambda=0.227) !0.227 GeV for NLO
101 parameter (Q2zero=2.0) !Gev^2 for u,d,s,g
102
103 data first /.TRUE./
104 data first_cent /.TRUE./
105
106 gaskelld 1.1 b = pt_b_param ! now parameter in input file
107
108 C DJG: Setup stuff for doing "central" cross section calculation. Here, I'm
109 C DJG: assuming we want the cross section at some "point" in Q2 and W space.
110 C DJG: I allow for binning in either z or pt2 - holding the other constant.
111 C DJG: Since the cross section has no phi-dependence, that part is easy -
112 C DJG: I can just ignore it.
113
114 if(doing_cent) then
115 nu = vertex0.nu
116 qx = vertex0.uq.x*vertex0.q
117 qy = vertex0.uq.y*vertex0.q
118 qz = vertex0.uq.z*vertex0.q
119 Q2 = vertex0.Q2
120 Eb = vertex0.Ein
121 Eprime = vertex0.e.E
122 Epx = vertex0.ue.x*Eprime
123 Epy = vertex0.ue.y*Eprime
124 Epz = vertex0.ue.z*Eprime
125 if(sigc_flag.eq.0) then ! binning in z
126 kcent = vertex0.zhad
127 gaskelld 1.1 pt2 = sigc_kin_ind*1.d6
128 zhad = 0.0
129 do i=1,sigc_nbin
130 klo = sigc_kin_min+(i-1)*(sigc_kin_max-sigc_kin_min)/sigc_nbin
131 khi = sigc_kin_min+i*(sigc_kin_max-sigc_kin_min)/sigc_nbin
132 if(vertex.zhad.gt.klo .and. vertex.zhad.le.khi) then
133 zhad = klo + (khi-klo)/2.
134 endif
135 enddo
136 elseif (sigc_flag.eq.1) then ! binning in pt2
137 kcent = vertex0.pt2
138 zhad = sigc_kin_ind
139 pt2 = 0.0
140 do i=1,sigc_nbin
141 klo = sigc_kin_min+(i-1)*(sigc_kin_max-sigc_kin_min)/sigc_nbin
142 klo = klo*1.d6
143 khi = sigc_kin_min+i*(sigc_kin_max-sigc_kin_min)/sigc_nbin
144 khi = khi*1.d6
145 if(vertex.pt2.gt.klo .and. vertex.pt2.le.khi) then
146 pt2 = klo + (khi-klo)/2.
147 endif
148 gaskelld 1.1 c write(6,*) 'chessy poofs',pt2,klo,khi
149 enddo
150 endif
151
152 if(first_cent) then
|
154 gaskelld 1.1 write(6,*) 'Ebeam (GeV):',Eb/1000
155 write(6,*) 'nu (GeV) :',nu/1000
156 write(6,*) 'Q2 (GeV2) :',Q2/1d6
157 if(sigc_flag.eq.0) then
158 write(6,*) 'Pt2 (GeV2) :',pt2/1d6
159 write(6,*) 'Binning in z from',sigc_kin_min,
160 1 'to',sigc_kin_max
161 elseif(sigc_flag.eq.1) then
162 write(6,*) 'z :',zhad
163 write(6,*) 'Binning in pt2 from',sigc_kin_min,
164 1 'to',sigc_kin_max,'GeV2'
165 endif
166 first_cent=.false.
167 endif
168 else
169 nu = vertex.nu
170 qx = vertex.uq.x*vertex.q
171 qy = vertex.uq.y*vertex.q
172 qz = vertex.uq.z*vertex.q
173 Q2 = vertex.Q2
174 Eb = vertex.Ein
175 gaskelld 1.1 Eprime = vertex.e.E
176 Epx = vertex.ue.x*Eprime
177 Epy = vertex.ue.y*Eprime
178 Epz = vertex.ue.z*Eprime
179 pt2 = vertex.pt2
180 zhad = vertex.zhad
181 endif
182
183 if(doing_semipi) then
184 mhad = Mpi
185 else if(doing_semika) then
186 mhad = Mk
187 else
188 write(6,*) 'semi_physics error: unknown hadron type'
189 stop
190 endif
191
192 mtar = targ.Mtar_struck
193
194 Ehad = zhad*nu
195 phad = sqrt(Ehad**2-mhad**2)
196 gaskelld 1.1
197 if(doing_cent) then
198 cthpq = sqrt(1.-pt2/phad**2)
199 else
200 cthpq = cos(vertex.theta_pq)
201 endif
202
203 sx = (2.*Eb*mtar + mtar**2)/1.d6 !convert to GeV2
204 if(do_fermi) then ! xbj = Q2/(2 P.q)
205 xbj = Q2/2./(efer*nu - abs(pfer)*(pferx*qx + pfery*qy + pferz*qz))
206 if(.not.doing_cent) then
207 ntup.xfermi = xbj
208 endif
209 else
210 xbj = Q2/2./mtar/nu
211 endif
212 if(xbj.gt.1.0) then
213 write(6,*) 'XBj is too large!', xbj
214 xbj=1.0
215 endif
216
217 gaskelld 1.1 c if(do_fermi) then ! y = P.q/P.k_beam
218 c y = (efer*nu - pferx*qx - pfery*qy - pferz*qz)/(efer*Eb)
219 c else
220 c y = nu/Eb
221 c endif
222
223 C DJG convert some stuff to GeV
224
225 Q2gev = Q2/1.d6
226 Qgev = sqrt(Q2gev)
227 pt2gev = pt2/1.d6
228
229 C Get the PDFs
230 if(first) then
231 call SetCtq5(iset) ! initialize Cteq5 (we're using cteq5m)
232 first=.FALSE.
233 endif
234
235 ipart=1
236 u = Ctq5pdf (ipart , xbj, Qgev)
237
238 gaskelld 1.1 ipart=-1
239 ubar = Ctq5pdf (ipart , xbj, Qgev)
240
241 ipart=2
242 d = Ctq5pdf (ipart , xbj, Qgev)
243
244 ipart=-2
245 dbar = Ctq5pdf (ipart , xbj, Qgev)
246
247 ipart=3
248 s = Ctq5pdf (ipart , xbj, Qgev)
249
250 ipart=-3
251 sbar = Ctq5pdf (ipart , xbj, Qgev)
252
253 sum_sq = qu**2*(u+ubar) + qd**2*(d+dbar) + qs**2*(s+sbar)
254
255 if (doing_deutsemi) then
256 sum_sq = sum_sq + qu**2*(d+dbar) + qd**2*(u+ubar) + qs**2*(s+sbar)
257 endif
258
259 gaskelld 1.1
260 C Simple paramaterization from Kretzer et al (EPJC 22 p. 269)
261 C for Q2=2.5.
262
263 c D_fav = 0.689*vertex.zhad**(-1.039)*(1.0-vertex.zhad)**1.241
264 c D_unfav = 0.217*vertex.zhad**(-1.805)*(1.0-vertex.zhad)**2.037
265
266 C
267
268 C Paramaterization from Binneweis et al (PRD 52 p.4947)
269
270 C sv is their scaling variable = ln( ln(Q2/Lambda^2)/ln(Q2_0/Lambda^2) )
271 C Q_0 sets the scale - for light quarks (u,d,s) they get Q_0 = sqrt(2) GeV
272 C Lambda is 227 MeV in NLO.
273
274 sv = log( log(Q2gev/lambda**2)/log(Q2zero/lambda**2) )
275
276 C Form of parameterization is D = N z^a1 (1-z)^a2
277
278 if(doing_semipi) then
279
280 gaskelld 1.1 N = 1.150 - 1.522*sv + 1.378*sv**2 - 0.527*sv**3
281 a1 = -0.740 - 1.680*sv + 1.546*sv**2 - 0.596*sv**3
282 a2 = 1.430 + 0.543*sv - 0.023*sv**2
283
284 Ns = 4.250 - 3.147*sv + 0.755*sv**2
285 a1s = -0.770 -0.573*sv + 0.117*sv**2
286 a2s = 4.48 + 0.890*sv - 0.138*sv**2
287
288 C This is Du (pi+ + pi-) = = Dd (pi+ + pi-) = D_favored + D_unfavored
289 D_sum = N*zhad**a1*(1.0-zhad)**a2
290 C This is Ds (pi+ + pi-)
291 D_sum_s = Ns*zhad**a1s*(1.0-zhad)**a2s
292
293
294 C Ratio of D-/D+ from P. Geiger's thesis (HERMES)
295
296 R_D = (1.0-zhad)**0.083583/(1.0+zhad)**1.9838
297
298 D_fav = D_sum/(1.0+R_D)
299 D_unfav = D_sum/(1.0+1.0/R_D)
300
301 gaskelld 1.1 C Assume Ds(pi+) = Ds(pi-) = Dsbar(pi+) = Dsbar(pi-)
302 C Note that this contrdicted by the HERMES data, but shouldn't make much
303 C difference for pions anyway.
304
305 D_s = D_sum_s/2.0
306
307
308 if(sum_sq.gt.0.) then
309 if(doing_hplus) then
310 dsigdz = (qu**2*u*D_fav + qu**2*ubar*D_unfav +
311 1 qd**2*d*D_unfav + qd**2*dbar*D_fav +
|
313 gaskelld 1.1 else
314 dsigdz = (qu**2*u*D_unfav + qu**2*ubar*D_fav +
315 1 qd**2*d*D_fav + qd**2*dbar*D_unfav +
316 2 qs**2*s*D_s + qs**2*sbar*D_s)/sum_sq
317 endif
318 if(doing_deutsemi) then
319 if(doing_hplus) then
320 dsigdz = dsigdz + (qu**2*d*D_fav + qu**2*dbar*D_unfav +
321 1 qd**2*u*D_unfav + qd**2*ubar*D_fav +
322 2 qs**2*s*D_s + qs**2*sbar*D_s)/sum_sq
323 else
324 dsigdz = dsigdz + (qu**2*d*D_unfav + qu**2*dbar*D_fav +
325 1 qd**2*u*D_fav + qd**2*ubar*D_unfav +
326 2 qs**2*s*D_s + qs**2*sbar*D_s)/sum_sq
327 endif
328 endif
329 else
330 dsigdz = 0.0
331 endif
332
333 elseif (doing_semika) then
334 gaskelld 1.1
335 N = 0.310 - 0.038*sv - 0.042*sv**2
336 a1 = -0.980 - 0.260*sv + 0.008*sv**2
337 a2 = 0.970 + 0.978*sv - 0.229*sv**2
338
339 Ns = 1.080 - 0.469*sv + 0.003*sv**2
340 a1s = -0.820 -0.240*sv - 0.035*sv**2
341 a2s = 2.550 + 1.026*sv - 0.246*sv**2
342
343 C This is Du (K+ + K-) = Ds(K+ + K-) = D_favored + D_unfavored (K+ + K-)
344 D_sum = N*zhad**a1*(1.0-zhad)**a2
345 C This is Dd (K+ + K-) (for convenience, I still call it D_sum_s)
346 D_sum_s = Ns*zhad**a1s*(1.0-zhad)**a2s
347
348
349 C Here I make the wild assumption that the ratio of unfavored to favored
350 C fragmentation functions for Kaons is the same as that for pions.
351
352 C Ratio of D-/D+ from P. Geiger's thesis (HERMES)
353
354 R_D = (1.0-zhad)**0.083583/(1.0+zhad)**1.9838
355 gaskelld 1.1
356 D_fav = D_sum/(1.0+R_D)
357 D_unfav = D_sum/(1.0+1.0/R_D)
358
359 C Assume Dd(K+) = Dd(K-) = Ddbar(K+) = Ddbar(K-)
360 C Note that this contrdicted by the HERMES data, but shouldn't make much
361 C difference for pions anyway.
362
363 D_s = D_sum_s/2.0
364
365 if(sum_sq.gt.0.) then
366 if(doing_hplus) then
367 dsigdz = (qu**2*u*D_fav + qu**2*ubar*D_unfav +
368 1 qd**2*d*D_s + qd**2*dbar*D_s +
369 2 qs**2*s*D_unfav + qs**2+sbar*D_fav)/sum_sq
370 else
371 dsigdz = (qu**2*u*D_unfav + qu**2*ubar*D_fav +
372 1 qd**2*d*D_s + qd**2*dbar*D_s +
373 2 qs**2*s*D_fav + qs**2*sbar*D_unfav)/sum_sq
374 endif
375 if(doing_deutsemi) then
376 gaskelld 1.1 if(doing_hplus) then
377 dsigdz = dsigdz + (qu**2*d*D_fav + qu**2*dbar*D_unfav +
378 1 qd**2*u*D_s + qd**2*ubar*D_s +
379 2 qs**2*s*D_unfav + qs**2*sbar*D_fav)/sum_sq
380 else
381 dsigdz = dsigdz + (qu**2*d*D_unfav + qu**2*dbar*D_fav +
382 1 qd**2*u*D_s + qd**2*ubar*D_s +
383 2 qs**2*s*D_fav + qs**2*sbar*D_unfav)/sum_sq
384 endif
385 endif
386 else
387 dsigdz = 0.0
388 endif
389
390
391 endif
392
393 sighad = dsigdz*b*exp(-b*pt2gev)/2./pi
394
395 C DJG: OK - I THINK I've got it right now.
396 C DJG: Should be dsig/dOmega dE
397 gaskelld 1.1
398 c sige = 2.*alpha**2*(y**2/2. + 1. - y - mtar*xbj*y/2./Eb)*
399 c 1 (Eprime/1000.)/(Q2gev*y*(mtar/1000.)*(nu/1000.))
400 c 2 * sum_sq
401
402 F2 = xbj*sum_sq
403 F1 = F2/2./xbj
404
405 W1 = F1/(mtar/1000.)
406 W2 = F2/(nu/1000.)
407
408 sin2th2 = Q2/4./Eb/Eprime
409 cos2th2 = 1.-sin2th2
410 if(do_fermi) then
411 W2coeff = (efer-abs(pfer)*pferz)*(efer-abs(pfer)*(pferx*Epx+pfery*Epy+pferz*Epz)/Eprime)/mtar**2 - sin2th2
412 else
413 W2coeff = cos2th2
414 endif
415
416 sige = 4.*alpha**2*(Eprime/1000)**2/Q2gev**2 * ( W2*W2coeff + 2.*W1*sin2th2)
417
418 gaskelld 1.1 C DJG This dsig/(dOmega_e dE_e dz dpt**2 dPhi_had) in microbarn/GeV**3/sr
419 sigsemi = sige*sighad*(hbarc/1000.)**2*10000.0
420
421
422 C Need to convert to dsig/ (dOmega_e dE_e dE_h dCos(theta) dPhi_had
423 C This is just given by 1/omega * 2*p_h**2*cos(theta)
424
425 jacobian = 1./(nu/1000.)*2.*(phad/1000.)**2*cthpq
426
427 C The 1.d6 converts from microbarn/GeV^2 to microbarn/MeV^2
428 sigma_eepiX = sigsemi*jacobian/1.d6
429
430
431 * Note that there is an additional factor 'fac' included with the fermi-smeared cross
432 * section. This takes into account pieces in the flux factor that are neglected (=1) in
433 * colinear collisions. The flux factor is |v_1-v_2| * 2E_1 * 2E_2.
434 * For a stationary target, v_2=0 and so velocity term is v_1=1 (electron
435 * beam), and E_2=M_2. For collinear boost, the flux factor can be expressed
436 * in a way that is lorenz invariant, and so can be used for lab or C.M.
437 * For a NON-COLLINEAR boost, there are two changes. First, the |v| term
438 * becomes 1 - (z component of pfer)/efer. Second, E_2 isn't just the mass,
439 gaskelld 1.1 * it becomes E_fermi, so we have to remove targ.Mtar_struck (which is used
440 * for E_2 by default) and replace it with efer. Since the flux factor
441 * comes in the denominator, we replace the usual flux factor (gtpr) with
442 * gtpr*fac, where fac = 1/ ( (1-pfer_z/efer)* (efer/mtar_struck) ).
443
444 if(do_fermi) then
445 fac = 1./(1.-pferz*pfer/efer) * mtar/efer
446 else
447 fac = 1.0
448 endif
449
450 sigma_eepiX = sigma_eepiX*fac
451
452 if(doing_cent) then
453 main.johnjac = jacobian*1000.0
454 else
455 main.davejac = jacobian*1000.0
456 ntup.sigcm = sighad
457 endif
458
459
460 gaskelld 1.1
461 peepiX = sigma_eepiX
462
463
464 c If doing_decay=.false., generate survival probability for main.weight.
465 c main.FP.p.path is dist. to back of detector, want decay prob. at aerogel.
466 C NOTE THAT ZAERO IS TAKEN WITH RESPECT TO THE POSITION AT WHICH PATHLEN
467 C IS CALCULATED (USUALLY THE BACK OF THE CALORIMETER). IF THE DRIFTS IN
468 C MC_SOS_HUT OR MC_HMS_HUT ARE CHANGED, THEN THE STARTING POINT MAY BE DIFFERENT AND THESE
469 C NUMBERS MAY BE WRONG. AERO BETWEEN S2Y AND S2X IN SOS AND BETWEEN DC2 S1X in HMS
470
471 C Beta/Gamma for decay need to use momentum after radiation/eloss, not vertex
472 C momentum. Get particle momentum from main.SP.p.delta
473
474 if (.not.doing_decay) then
475 if (hadron_arm.eq.1) then
476 zaero = -331.491 !aero at 40.199 cm, last project=371.69 cm
477 else if (hadron_arm.eq.2) then
478 * zaero = -76. !aero at 270cm,last project=346(cal).
479 zaero = -82.8 !From Rick: aero at 263.2cm,last project=346(cal).
480 else if (hadron_arm.eq.3) then
481 gaskelld 1.1 zaero = -183. !aero at 130cm,last project=313(S2)
482 else if (hadron_arm.eq.4) then
483 zaero = -183.
484 endif
485 pathlen = main.FP.p.path + zaero*(1+main.FP.p.dx**2+main.FP.p.dy**2)
486 p_kaon = spec.p.P*(1.+main.SP.p.delta/100.)
487 betak = spec.p.P/sqrt(spec.p.P**2+Mh2)
488 gammak = 1./sqrt(1.-betak**2)
489 survivalprob = 1./exp(pathlen/(ctau*betak*gammak))
490 decdist = survivalprob !decdist in ntuple
491 endif
492
493 return
494
495 end
|
Return to
semi_physics.f
CVS log
Up to [HallC] / simc_semi