1 gaskelld 1.1 subroutine target_init(using_Eloss,using_Coulomb,the_cent,thp_cent,
2 > Pp_cent,Mh2,ebeam,Pe_cent)
3
4 implicit none
5 include 'constants.inc'
6 include 'target.inc'
7
8 integer i
9 real*8 the_cent, thp_cent, Pp_cent, Mh2, ebeam, energy
10 real*8 Pe_cent
11 real*8 za2, fc
12 logical using_Eloss, using_Coulomb
13
14 real*8 zero
15 parameter (zero=0.0d0) !double precision zero for subroutines calls.
16
17 ! The radiation length of the target
18
19 targ.L1 = log(184.15) - log(targ.Z)/3.0
20 targ.L2 = log(1194.) - 2.*log(targ.Z)/3.0
21 if(targ.Z.eq.1)then
22 gaskelld 1.1 targ.L1=5.31
23 targ.L2=6.144
24 endif
25 za2 = (targ.Z*alpha)**2
26 fc = za2*(1.202+za2*(-1.0369+za2*1.008/(za2+1)))
27
28 ! ... The radiation length, in both g/cm2 and cm; For H,D,4He, using
29 ! ... PDB values (1999), other calculated directly from Tsai formula
30 ! ... (For 4He, the Tsai formula is 10% lower).
31
32 if (nint(targ.A).eq.1) then
33 targ.X0 = 61.28
34 else if (nint(targ.A).eq.2) then
35 targ.X0 = 122.4
36 else if (nint(targ.A).eq.4) then
37 targ.X0 = 94.32
38 else if (nint(targ.A).eq.3) then
39 write(6,*) 'Using Tsai formula for He-3 radiation length!!! (is there a better value?)'
40 targ.X0 = 716.405*targ.A/targ.Z/(targ.Z*(targ.L1-fc)+targ.L2)
41 else
42 targ.X0 = 716.405*targ.A/targ.Z/(targ.Z*(targ.L1-fc)+targ.L2)
43 gaskelld 1.1 endif
44 targ.X0_cm = targ.X0/targ.rho
45
46 ! ... 'Average' ionization losses (MOST PROBABLE, REALLY). Printed out by simc
47 ! ... (so ave or m.p. are both OK), and Ebeam_vertex_ave (which should be m.p.)
48
49 energy = ebeam
50 call trip_thru_target (1,zero,energy,zero,targ.Eloss(1).ave,
51 > targ.teff(1).ave,Me,4)
52 energy = Pe_cent**2
53 call trip_thru_target (2,zero,energy,the_cent,targ.Eloss(2).ave,
54 > targ.teff(2).ave,Me,4)
55 energy = sqrt(Pp_cent**2 + Mh2)
56 call trip_thru_target (3,zero,energy,thp_cent,targ.Eloss(3).ave,
57 > targ.teff(3).ave,sqrt(Mh2),4)
58 if (.not.using_Eloss) then
59 do i = 1, 3
60 targ.Eloss(i).ave = 0.0
61 enddo
62 endif
63
64 gaskelld 1.1 ! ... Coulomb potential energy (NB All of the following are positive)
65
66 if (using_Coulomb) then
67 targ.Coulomb.ave=6./5.*(targ.Z-1.)*alpha*hbarc/(1.18*targ.A**(1./3.))
68 targ.Coulomb_constant = 5./12. * targ.Coulomb.ave
69 targ.Coulomb.min = targ.Coulomb_constant * 2.0
70 targ.Coulomb.max = targ.Coulomb_constant * 3.0
71 else
72 targ.Coulomb.ave = 0.0
73 targ.Coulomb_constant = 0.0
74 targ.Coulomb.min = 0.0
75 targ.Coulomb.max = 0.0
76 endif
77
78 return
79 end
80
81 !----------------------------------------------------------------------
82
83 subroutine limits_init(H)
84
85 gaskelld 1.1 implicit none
86 include 'simulate.inc'
87 include 'radc.inc'
88 include 'histograms.inc'
89 include 'sf_lookup.inc' !need to know Em range of spec. fcn.
90
91 integer i
92 real*8 r
93 real*8 Ebeam_min, Ebeam_max
94 real*8 t1,t2 !temp. variables.
95 real*8 slop_Coulomb, slop_Ebeam, slop_Ee, slop_Ep
96 record /cuts/ the_phys, thp_phys, z, pp
97 record /histograms/ H
98
99 !-----------------------------------------------------------------------
100 ! RECORDS store information relating to various quantities at several
101 ! 'STATES' in event description:
102 ! (of course, not all qties are _recorded_ in each of these stages)
103 !
104 ! orig.* qties are TRUE qties but before ANY radiation
105 ! vertex.* qties are those used to determine spec fn and sigcc
106 gaskelld 1.1 ! weighting, so TRUE qties before tails 2,3 radiated
107 ! recon.* qties are defined AT the spectrometers, AFTER being
108 ! run through the single arm montecarlos
109 !
110 ! main.SP.* qties are defined AT the spectrometers, so TRUE
111 ! qties AFTER any modifications due to non-radiative
112 ! interaction with the target (this includes Coulomb
113 ! accel/decel of initial/final electron)
114 ! main.FP.* spect. focal plane values-if using spect. model
115 ! main.RECON.* qties are defined AT the spectrometers, AFTER being
116 ! run through the single arm montecarlos
117 !
118 ! Note that in the absence of radiation, GEN=VERTEX
119 ! Note that if a spectrometer is not used, RECON=SP, FP is empty.
120 ! Note that the generated qties we have to compare with gen limits later
121 ! are unaffected by radiation in tail1, so equal to vertex qties.
122 !
123 !
124 ! EDGES on event qties at any stage can be derived from these cuts, by taking
125 ! the transformations (and associated uncertainties) between each stage into
126 ! account. The general usefulness of defining these edges is to improve
127 gaskelld 1.1 ! efficiency --> provide checks at intermediate stages to allow us to abort
128 ! an event early, and provide some foreknowledge in generation routines of
129 ! what will make it to the end. We must be certain that these derived EDGES
130 ! are wide enough that the events making it through the cuts don't come
131 ! uncomfortably close to them.
132 !
133 ! SPedge.* Limits at spectrometer = Acceptance + slop
134 ! edge.* Limits after interaction = SPedge + musc or eloss
135 ! = VERTEXedge + coulomb + rad.
136 ! VERTEXedge.* Limits at vertex (from energy conservation).
137 ! gen.* Generation limits.
138 !
139 ! SLOPS
140 ! slop.MC.* Due to spectrometer resolution (spectrometer MC).
141 ! slop_* Slop in calculated physics quantities (from spect.
142 ! resolution + range of eloss/coulomb/...
143 !
144 !----------------------------------------------------------------------
145
146 ! ... get slop values for the difference between orig and recon
147 ! ... SPECTROMETER quantities (recon inaccuracies due to the montecarlos)
148 gaskelld 1.1
149 if (using_E_arm_montecarlo) then
150 if (electron_arm.eq.1) then
151 slop.MC.e.delta.used = slop_param_d_HMS
152 slop.MC.e.yptar.used = slop_param_t_HMS
153 slop.MC.e.xptar.used = slop_param_p_HMS
154 else if (electron_arm.eq.2) then
155 slop.MC.e.delta.used = slop_param_d_SOS
156 slop.MC.e.yptar.used = slop_param_t_SOS
157 slop.MC.e.xptar.used = slop_param_p_SOS
158 else if (electron_arm.eq.3) then
159 slop.MC.e.delta.used = slop_param_d_HRSR
160 slop.MC.e.yptar.used = slop_param_t_HRSR
161 slop.MC.e.xptar.used = slop_param_p_HRSR
162 else if (electron_arm.eq.4) then
163 slop.MC.e.delta.used = slop_param_d_HRSL
164 slop.MC.e.yptar.used = slop_param_t_HRSL
165 slop.MC.e.xptar.used = slop_param_p_HRSL
166 else if (electron_arm.eq.5 .or. electron_arm.eq.6) then
167 slop.MC.e.delta.used = slop_param_d_SHMS
168 slop.MC.e.yptar.used = slop_param_t_SHMS
169 gaskelld 1.1 slop.MC.e.xptar.used = slop_param_p_SHMS
170 endif
171 endif
172 if (using_P_arm_montecarlo) then
173 if (hadron_arm.eq.1) then
174 slop.MC.p.delta.used = slop_param_d_HMS
175 slop.MC.p.yptar.used = slop_param_t_HMS
176 slop.MC.p.xptar.used = slop_param_p_HMS
177 else if (hadron_arm.eq.2) then
178 slop.MC.p.delta.used = slop_param_d_SOS
179 slop.MC.p.yptar.used = slop_param_t_SOS
180 slop.MC.p.xptar.used = slop_param_p_SOS
181 else if (hadron_arm.eq.3) then
182 slop.MC.p.delta.used = slop_param_d_HRSR
183 slop.MC.p.yptar.used = slop_param_t_HRSR
184 slop.MC.p.xptar.used = slop_param_p_HRSR
185 else if (hadron_arm.eq.4) then
186 slop.MC.p.delta.used = slop_param_d_HRSL
187 slop.MC.p.yptar.used = slop_param_t_HRSL
188 slop.MC.p.xptar.used = slop_param_p_HRSL
189 else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then
190 gaskelld 1.1 slop.MC.p.delta.used = slop_param_d_SHMS
191 slop.MC.p.yptar.used = slop_param_t_SHMS
192 slop.MC.p.xptar.used = slop_param_p_SHMS
193 endif
194 endif
195
196 ! ... Add slop to SPedges. Used as input to final edge.*.*.* and to
197 ! ... calculate minimum/maximum theta_phys for extreme_trip_thru_target.
198
199 SPedge.e.delta.min = SPedge.e.delta.min - slop.MC.e.delta.used
200 SPedge.e.delta.max = SPedge.e.delta.max + slop.MC.e.delta.used
201 SPedge.e.yptar.min = SPedge.e.yptar.min - slop.MC.e.yptar.used
202 SPedge.e.yptar.max = SPedge.e.yptar.max + slop.MC.e.yptar.used
203 SPedge.e.xptar.min = SPedge.e.xptar.min - slop.MC.e.xptar.used
204 SPedge.e.xptar.max = SPedge.e.xptar.max + slop.MC.e.xptar.used
205 SPedge.p.delta.min = SPedge.p.delta.min - slop.MC.p.delta.used
206 SPedge.p.delta.max = SPedge.p.delta.max + slop.MC.p.delta.used
207 SPedge.p.yptar.min = SPedge.p.yptar.min - slop.MC.p.yptar.used
208 SPedge.p.yptar.max = SPedge.p.yptar.max + slop.MC.p.yptar.used
209 SPedge.p.xptar.min = SPedge.p.xptar.min - slop.MC.p.xptar.used
210 SPedge.p.xptar.max = SPedge.p.xptar.max + slop.MC.p.xptar.used
211 gaskelld 1.1
212 ! Compute TRUE edges -- distortions in the target come into play
213
214 edge.e.E.min = (1.+SPedge.e.delta.min/100.)*spec.e.P +
215 > targ.Coulomb.min - dE_edge_test
216 edge.e.E.max = (1.+SPedge.e.delta.max/100.)*spec.e.P +
217 > targ.Coulomb.max + dE_edge_test
218 pp.min = (1.+SPedge.p.delta.min/100.)*spec.p.P - dE_edge_test
219 pp.max = (1.+SPedge.p.delta.max/100.)*spec.p.P + dE_edge_test
220 pp.min = max(0.001d0,pp.min) !avoid p=0 (which can lead to div by zero)(
221 edge.p.E.min = sqrt(pp.min**2 + Mh2)
222 edge.p.E.max = sqrt(pp.max**2 + Mh2)
223
224 ! ... extreme theta_e,theta_p,z values for extreme_trip_thru_target.
225 the_phys.max = acos( (spec.e.cos_th-spec.e.sin_th*SPedge.e.yptar.max)/
226 > sqrt(1.+SPedge.e.yptar.max**2+SPedge.e.xptar.max**2) )
227 the_phys.min = acos( (spec.e.cos_th-spec.e.sin_th*SPedge.e.yptar.min)/
228 > sqrt(1.+SPedge.e.yptar.min**2) )
229 thp_phys.max = acos( (spec.p.cos_th-spec.p.sin_th*SPedge.p.yptar.max)/
230 > sqrt(1.+SPedge.p.yptar.max**2+SPedge.p.xptar.max**2) )
231 thp_phys.min = acos( (spec.p.cos_th-spec.p.sin_th*SPedge.p.yptar.min)/
232 gaskelld 1.1 > sqrt(1.+SPedge.p.yptar.min**2) )
233 z.min = -0.5*targ.length
234 z.max = 0.5*targ.length
235 call extreme_trip_thru_target(Ebeam, the_phys, thp_phys, edge.e.E,
236 > pp, z, Mh)
237 if (.not.using_Eloss) then
238 do i = 1, 3
239 targ.Eloss(i).min = 0.0
240 targ.Eloss(i).max = 0.0
241 enddo
242 endif
243
244 if (.not.mc_smear) then
245 targ.musc_max(1)=0.
246 targ.musc_max(2)=0.
247 targ.musc_max(3)=0.
248 endif
249
250 edge.e.E.min = edge.e.E.min + targ.Eloss(2).min
251 edge.e.E.max = edge.e.E.max + targ.Eloss(2).max
252 edge.p.E.min = edge.p.E.min + targ.Eloss(3).min
253 gaskelld 1.1 edge.p.E.max = edge.p.E.max + targ.Eloss(3).max
254 edge.e.yptar.min = SPedge.e.yptar.min - targ.musc_max(2)
255 edge.e.yptar.max = SPedge.e.yptar.max + targ.musc_max(2)
256 edge.e.xptar.min = SPedge.e.xptar.min - targ.musc_max(2)
257 edge.e.xptar.max = SPedge.e.xptar.max + targ.musc_max(2)
258 edge.p.yptar.min = SPedge.p.yptar.min - targ.musc_max(3)
259 edge.p.yptar.max = SPedge.p.yptar.max + targ.musc_max(3)
260 edge.p.xptar.min = SPedge.p.xptar.min - targ.musc_max(3)
261 edge.p.xptar.max = SPedge.p.xptar.max + targ.musc_max(3)
262
263 ! Edges on values of Em and Pm BEFORE reconstruction. Need to apply slop to
264 ! take into account all transformations from ORIGINAL TRUE values to
265 ! RECONSTRUCTED TRUE values --> that includes: (a) reconstruction slop on
266 ! spectrometer values (b) variation in the beam energy (c) difference between
267 ! actual losses and the corrections made for them
268
269 ! ... Save the 'measured' beam energy, which will be used in recon.
270
271 Ebeam_vertex_ave = Ebeam + targ.Coulomb.ave - targ.Eloss(1).ave
272
273 ! ... Calculate slop in energies and momenta.
274 gaskelld 1.1
275 Ebeam_max = Ebeam + dEbeam/2. - targ.Eloss(1).min + targ.Coulomb.max
276 Ebeam_min = Ebeam - dEbeam/2. - targ.Eloss(1).max + targ.Coulomb.min
277 slop_Coulomb = targ.Coulomb.max - targ.Coulomb.ave
278 slop_Ebeam = dEbeam/2. + slop_Coulomb
279 slop_Ee = slop.MC.e.delta.used/100.*spec.e.P + slop_Coulomb
280 r = sqrt(edge.p.E.max**2 - Mh2)
281 slop_Ep = sqrt( (r + slop.MC.p.delta.used/100.*spec.p.P)**2 + Mh2 ) -
282 > edge.p.E.max
283 slop_Ebeam = slop_Ebeam + (targ.Eloss(1).max-targ.Eloss(1).min)
284 slop_Ee = slop_Ee + (targ.Eloss(2).max-targ.Eloss(2).min)
285 slop_Ep = slop_Ep + (targ.Eloss(3).max-targ.Eloss(3).min)
286
287 if (doing_heavy) then ! 'reconstructed' Em cuts.
288 slop.total.Em.used = slop_Ebeam + slop_Ee + slop_Ep + dE_edge_test
289 edge.Em.min = cuts.Em.min - slop.total.Em.used
290 edge.Em.max = cuts.Em.max + slop.total.Em.used
291 edge.Em.min = max(0.d0,edge.Em.min)
292 endif
293
294 ! Edges on Em, Pm, etc... VERTEXedge.* values are vertex limits. edge.* values
295 gaskelld 1.1 ! are limits at spectrometer (after eloss and e'/hadron radiation). Remember
296 ! that min/max values are initialized to wide open values in STRUCTURES.
297 ! Need edge.Em to limit radiated photon energies.
298 ! Need VERTEXedge.Em for calulating limits on radiatied photons.
299 ! Need VERTEXedge.Pm for A(e,e'p) only (for generation limits on Ee', Ep').
300 ! Note that Pm_theory(*).min/max might allow for positive and negative Pm,
301 ! or it could have positive only. We want *.Pm.min to be zero if zero is
302 ! allowed, so we have to check for Pm_theory.min being negative.
303
304 ! For all cases:
305 ! 1. Start by giving Em, Pm limits (=0 for hydrogen, based on S.F. limits
306 ! for all other cases except A(e,e'p). For this case, Pm limits come
307 ! from spectral function. Em limit (max) has to come from energy
308 ! conservation after everything else is calculated.
309 !
310
311 if (doing_hyd_elast) then
312 VERTEXedge.Em.min = 0.0
313 VERTEXedge.Em.max = 0.0
314 VERTEXedge.Pm.min = 0.0
315 VERTEXedge.Pm.max = 0.0
316 gaskelld 1.1 else if (doing_deuterium) then
317 VERTEXedge.Em.min = Mp + Mn - targ.M !2.2249 MeV, I hope.
318 VERTEXedge.Em.max = Mp + Mn - targ.M
319 VERTEXedge.Pm.min = 0.0
320 VERTEXedge.Pm.max = max(abs(Pm_theory(1).min),abs(Pm_theory(1).max))
321 else if (doing_heavy) then
322 VERTEXedge.Pm.min=0.0
323 VERTEXedge.Pm.max=0.0
324 do i = 1, nrhoPm
325 t1=max(abs(Pm_theory(i).min),abs(Pm_theory(i).max))
326 VERTEXedge.Pm.max = max(VERTEXedge.Pm.max,t1)
327 enddo
328 VERTEXedge.Em.min = E_Fermi
329 VERTEXedge.Em.max = 1000. !Need Egamma_tot_max for good limit.
330 else if (doing_hydpi .or. doing_hydkaon .or. doing_hyddelta .or. doing_hydrho) then
331 VERTEXedge.Em.min = 0.0
332 VERTEXedge.Em.max = 0.0
333 VERTEXedge.Pm.min = 0.0
334 VERTEXedge.Pm.max = 0.0
335 else if (doing_deutpi .or. doing_deutkaon .or. doing_deutdelta .or. doing_deutrho) then
336 VERTEXedge.Em.min = Mp + Mn - targ.M !2.2249 MeV, I hope.
337 gaskelld 1.1 VERTEXedge.Em.max = Mp + Mn - targ.M
338 VERTEXedge.Pm.min = 0.0
339 VERTEXedge.Pm.max = pval(nump)
340 else if (doing_hepi .or. doing_hekaon .or. doing_hedelta .or. doing_herho) then
341 VERTEXedge.Em.min = targ.Mtar_struck + targ.Mrec - targ.M
342 VERTEXedge.Em.max = Emval(numEm)
343 VERTEXedge.Pm.min = 0.0
344 VERTEXedge.Pm.max = pval(nump)
345 else if (doing_semi) then
346 VERTEXedge.Em.min = 0.0
347 VERTEXedge.Em.max = 0.0
348 VERTEXedge.Pm.min = 0.0
349 VERTEXedge.Pm.max = 0.0
350
351 endif
352 write(6,*) 'E_bind =',VERTEXedge.Em.min,'MeV in limits_init (QF only)'
353
354
355 ! Calculate limits for recoiling (A-1) system, if there is one.
356 ! M_{A-1} (Mrec) limits from Em range, since M_{A-1} = M_A - M_N + E_m
357 ! T_{A-1} (Trec) limits from Mrec and Pm limits, since
358 gaskelld 1.1 ! E_rec=sqrt(M_rec**2+P_rec**2), and P_rec = -P_m
359
360 if (doing_hyd_elast .or. doing_hydpi .or. doing_hydkaon .or.
361 1 doing_hyddelta .or. doing_hydrho) then
362 VERTEXedge.Mrec.min = 0.0
363 VERTEXedge.Mrec.max = 0.0
364 VERTEXedge.Trec.min = 0.0
365 VERTEXedge.Trec.max = 0.0
366 else
367 VERTEXedge.Mrec.min = targ.M - targ.Mtar_struck + VERTEXedge.Em.min
368 VERTEXedge.Mrec.max = targ.M - targ.Mtar_struck + VERTEXedge.Em.max
369 VERTEXedge.Trec.min = sqrt(VERTEXedge.Mrec.max**2+VERTEXedge.Pm.min**2) -
370 > VERTEXedge.Mrec.max
371 VERTEXedge.Trec.max = sqrt(VERTEXedge.Mrec.min**2+VERTEXedge.Pm.max**2) -
372 > VERTEXedge.Mrec.min
373 endif
374
375 ! For pion(kaon) production, Trec_struck is T of the recoiling nucleon(hyperon).
376 ! Can only get limits from general energy conservation, and can only get
377 ! upper limit, since the lower limit is determined by the allowed radiation,
378 ! which is not calculated yet (and needs Trec to be calculated).
379 gaskelld 1.1
380 if (doing_eep) then
381 VERTEXedge.Trec_struck.min = 0.
382 VERTEXedge.Trec_struck.max = 0.
383 else
384 VERTEXedge.Trec_struck.min = 0.
385 VERTEXedge.Trec_struck.max = Ebeam_max + targ.Mtar_struck -
386 > targ.Mrec_struck - edge.e.E.min - edge.p.E.min -
387 > VERTEXedge.Em.min - VERTEXedge.Trec.min
388 endif
389
390 ! Get radiation limits. In all cases, total energy conservation gives
391 ! the primary limit. The doing_heavy case has a second condition. Radiated
392 ! photons change Em from the vertex value to the measured value. They also
393 ! change Pm, which can modify Trec. So the max. change in Em is for a
394 ! minimum generated Em (VERTEXedge.Em.min) that ends up as a maximum measured Em
395 ! (edge.Em.max), with slop to take into account the modification to Trec.
396
397 Egamma_tot_max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck -
398 > edge.e.E.min - edge.p.E.min - VERTEXedge.Em.min -
399 > VERTEXedge.Trec.min - VERTEXedge.Trec_struck.min
400 gaskelld 1.1 if (doing_heavy) then
401 t2 = (edge.Em.max - VERTEXedge.Em.min) +
402 > (VERTEXedge.Trec.max - VERTEXedge.Trec.min)
403 Egamma_tot_max = min(Egamma_tot_max,t2)
404 endif
405
406 ! ... override calculated limits with hardwired value if desired.
407 if (hardwired_rad) Egamma_tot_max = Egamma_gen_max
408 if (.not.using_rad) Egamma_tot_max = 0.0
409 if (doing_tail(1)) Egamma1_max = Egamma_tot_max
410 if (doing_tail(2)) Egamma2_max = Egamma_tot_max
411 if (doing_tail(3)) Egamma3_max = Egamma_tot_max
412
413 if (doing_heavy) then !Needed Egamma_tot_max for final limits.
414 VERTEXedge.Em.min = max(VERTEXedge.Em.min,edge.Em.min-Egamma_tot_max)
415 VERTEXedge.Em.max = min(VERTEXedge.Em.max,edge.Em.max)
416 endif
417
418 ! ... compute edge on summed _generated_ energies based on predicted VERTEX
419 ! ... and TRUE values of Em (Ee'+Ep' for doing_heavy, Ee' for pion/kaon,
420 ! ... Ee' for D(e,e'p), not used for hydrogen elastic.)
421 gaskelld 1.1 ! ... Primary limits from energy conservation at vertex, seconday limits from
422 ! ... spectrometer limits modified by radiation.
423
424 if (doing_hyd_elast) then !NO generated energies.
425 gen.sumEgen.min = 0.0
426 gen.sumEgen.max = 0.0
427
428 else if (doing_heavy) then !generated TOTAL (e+p) energy limits.
429 gen.sumEgen.max = Ebeam_max + targ.Mtar_struck -
430 > VERTEXedge.Trec.min - VERTEXedge.Em.min
431 gen.sumEgen.min = Ebeam_min + targ.Mtar_struck -
432 > VERTEXedge.Trec.max - VERTEXedge.Em.max - Egamma1_max
433 gen.sumEgen.max = min(gen.sumEgen.max,edge.e.E.max+edge.p.E.max+Egamma_tot_max)
434 gen.sumEgen.min = max(gen.sumEgen.min,edge.e.E.min+edge.p.E.min)
435
436 else if (doing_semi) then
437 gen.sumEgen.max = Ebeam_max - VERTEXedge.Trec.min - VERTEXedge.Trec_struck.min
438 gen.sumEgen.min = Ebeam_min - VERTEXedge.Trec.max - VERTEXedge.Trec_struck.max
439
440 else !generated ELECTRON energy limits.
441 gen.sumEgen.max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck -
442 gaskelld 1.1 > edge.p.E.min - VERTEXedge.Em.min - VERTEXedge.Trec.min -
443 > VERTEXedge.Trec_struck.min
444 gen.sumEgen.min = Ebeam_min + targ.Mtar_struck - targ.Mrec_struck -
445 > edge.p.E.max - VERTEXedge.Em.max - VERTEXedge.Trec.max -
446 > VERTEXedge.Trec_struck.max - Egamma_tot_max
447 gen.sumEgen.max = min(gen.sumEgen.max, edge.e.E.max+Egamma2_max)
448 gen.sumEgen.min = max(gen.sumEgen.min, edge.e.E.min)
449 endif
450
451 gen.sumEgen.min = gen.sumEgen.min - dE_edge_test
452 gen.sumEgen.max = gen.sumEgen.max + dE_edge_test
453 gen.sumEgen.min = max(0.d0,gen.sumEgen.min)
454
455 ! ... E arm GENERATION limits from sumEgen.
456 ! ... Not used for doing_hyd_elast, but define for the hardwired histograms.
457
458 if (doing_hyd_elast) then
459 gen.e.E.min = edge.e.E.min
460 gen.e.E.max = edge.e.E.max + Egamma2_max
461 else if (doing_deuterium .or. doing_pion .or. doing_kaon
462 1 .or. doing_rho .or. doing_delta) then
463 gaskelld 1.1 gen.e.E.min = gen.sumEgen.min
464 gen.e.E.max = gen.sumEgen.max
465 else if (doing_heavy .or. doing_semi) then
466 gen.e.E.min = gen.sumEgen.min - edge.p.E.max - Egamma3_max
467 gen.e.E.max = gen.sumEgen.max - edge.p.E.min
468 endif
469
470 ! ... Apply limits from direct comparison to acceptance.
471 gen.e.E.min = max(gen.e.E.min, edge.e.E.min)
472 gen.e.E.max = min(gen.e.E.max, edge.e.E.max + Egamma2_max)
473
474 gen.e.delta.min = (gen.e.E.min/spec.e.P-1.)*100.
475 gen.e.delta.max = (gen.e.E.max/spec.e.P-1.)*100.
476 gen.e.yptar.min = edge.e.yptar.min
477 gen.e.yptar.max = edge.e.yptar.max
478 gen.e.xptar.min = edge.e.xptar.min
479 gen.e.xptar.max = edge.e.xptar.max
480
481 ! ... P arm GENERATION limits from sumEgen. Not used for any case
482 ! ... except doing_heavy, but need to define for code that writes out limits.
483
484 gaskelld 1.1 if (doing_hyd_elast.or.doing_deuterium.or.doing_pion.or.doing_kaon .or.
485 1 doing_rho .or. doing_delta) then
486 gen.p.E.min = edge.p.E.min
487 gen.p.E.max = edge.p.E.max + Egamma3_max
488 else if (doing_heavy .or. doing_semi)then
489 gen.p.E.min = gen.sumEgen.min - edge.e.E.max - Egamma2_max
490 gen.p.E.max = gen.sumEgen.max - edge.e.E.min
491 endif
492 ! ... Apply limits from direct comparison to acceptance.
493 gen.p.E.min = max(gen.p.E.min, edge.p.E.min)
494 gen.p.E.max = min(gen.p.E.max, edge.p.E.max + Egamma3_max)
495
496 gen.p.delta.min = (sqrt(gen.p.E.min**2-Mh2)/spec.p.P-1.)*100.
497 gen.p.delta.max = (sqrt(gen.p.E.max**2-Mh2)/spec.p.P-1.)*100.
498 gen.p.yptar.min = edge.p.yptar.min
499 gen.p.yptar.max = edge.p.yptar.max
500 gen.p.xptar.min = edge.p.xptar.min
501 gen.p.xptar.max = edge.p.xptar.max
502
503 ! Axis specs for the diagnostic histograms
504 H.gen.e.delta.min = gen.e.delta.min
505 gaskelld 1.1 H.gen.e.yptar.min = gen.e.yptar.min
506 H.gen.e.xptar.min = -gen.e.xptar.max
507 H.gen.e.delta.bin = (gen.e.delta.max-gen.e.delta.min)/float(nHbins)
508 H.gen.e.yptar.bin = (gen.e.yptar.max-gen.e.yptar.min)/float(nHbins)
509 H.gen.e.xptar.bin = (gen.e.xptar.max-gen.e.xptar.min)/float(nHbins)
510 H.gen.p.delta.min = gen.p.delta.min
511 H.gen.p.yptar.min = gen.p.yptar.min
512 H.gen.p.xptar.min = -gen.p.xptar.max
513 H.gen.p.delta.bin = (gen.p.delta.max-gen.p.delta.min)/float(nHbins)
514 H.gen.p.yptar.bin = (gen.p.yptar.max-gen.p.yptar.min)/float(nHbins)
515 H.gen.p.xptar.bin = (gen.p.xptar.max-gen.p.xptar.min)/float(nHbins)
516 H.gen.Em.min = VERTEXedge.Em.min
517 H.gen.Em.bin = (max(100.d0,VERTEXedge.Em.max) - VERTEXedge.Em.min)/float(nHbins)
518 H.gen.Pm.min = VERTEXedge.Pm.min
519 H.gen.Pm.bin = (max(100.d0,VERTEXedge.Pm.max) - VERTEXedge.Pm.min)/float(nHbins)
520
521 H.geni.e.delta.min = H.gen.e.delta.min
522 H.geni.e.yptar.min = H.gen.e.yptar.min
523 H.geni.e.xptar.min = H.gen.e.xptar.min
524 H.geni.e.delta.bin = H.gen.e.delta.bin
525 H.geni.e.yptar.bin = H.gen.e.yptar.bin
526 gaskelld 1.1 H.geni.e.xptar.bin = H.gen.e.xptar.bin
527 H.geni.p.delta.min = H.gen.p.delta.min
528 H.geni.p.yptar.min = H.gen.p.yptar.min
529 H.geni.p.xptar.min = H.gen.p.xptar.min
530 H.geni.p.delta.bin = H.gen.p.delta.bin
531 H.geni.p.yptar.bin = H.gen.p.yptar.bin
532 H.geni.p.xptar.bin = H.gen.p.xptar.bin
533 H.geni.Em.min = H.gen.Em.min
534 H.geni.Em.bin = H.gen.Em.bin
535 H.geni.Pm.min = H.gen.Pm.min
536 H.geni.Pm.bin = H.gen.Pm.bin
537
538 H.RECON.e.delta.min = H.gen.e.delta.min
539 H.RECON.e.yptar.min = H.gen.e.yptar.min
540 H.RECON.e.xptar.min = H.gen.e.xptar.min
541 H.RECON.e.delta.bin = H.gen.e.delta.bin
542 H.RECON.e.yptar.bin = H.gen.e.yptar.bin
543 H.RECON.e.xptar.bin = H.gen.e.xptar.bin
544 H.RECON.p.delta.min = H.gen.p.delta.min
545 H.RECON.p.yptar.min = H.gen.p.yptar.min
546 H.RECON.p.xptar.min = H.gen.p.xptar.min
547 gaskelld 1.1 H.RECON.p.delta.bin = H.gen.p.delta.bin
548 H.RECON.p.yptar.bin = H.gen.p.yptar.bin
549 H.RECON.p.xptar.bin = H.gen.p.xptar.bin
550 H.RECON.Em.min = H.gen.Em.min
551 H.RECON.Em.bin = H.gen.Em.bin
552 H.RECON.Pm.min = H.gen.Pm.min
553 H.RECON.Pm.bin = H.gen.Pm.bin
554
555 return
556 end
557
558 !------------------------------------------------------------------------
559
560 subroutine radc_init
561
562 implicit none
563 include 'simulate.inc'
564 include 'radc.inc'
565 include 'brem.inc'
566
567 !--------------------------------------------------------------
568 gaskelld 1.1 !
569 ! First, about those (mysterious) 2 main 'radiative option' flags ...
570 !
571 ! The significance of RAD_FLAG:
572 ! RAD_FLAG = 0 .. use best available formulas, generate in
573 ! .. (ntail,Egamma) basis
574 ! = 1 .. use BASICRAD only, generate in (ntail,Egamma)
575 ! .. basis
576 ! = 2 .. use BASICRAD only, generate in (Egamma(1,2,3))
577 ! .. basis but prevent overlap of tails (bogus, note)
578 ! = 3 .. use BASICRAD only, generate in (Egamma(1,2,3))
579 ! .. allowing radiation in all 3 directions
580 ! .. simultaneously
581 ! The (ntail,Egamma) basis can be called the PEAKED basis since it allows
582 ! only 3 photon directions. PEAKED_BASIS_FLAG is set to zero below when
583 ! the peaked basis is being used, in this way we can conveniently tell
584 ! the BASICRAD routine to use the full Egamma formula to generate the gamma
585 ! energy whenever it's called.
586 !
587 ! (See N. Makins' thesis, section 4.5.6, for more help on this point)
588 !
589 gaskelld 1.1 ! The significance of EXTRAD_FLAG:
590 ! EXTRAD_FLAG = 1 .. use only BASIC external radiation formulas
591 ! .. (phi = 1)
592 ! = 2 .. use BASIC ext rad formulas x phi
593 ! = 3 .. use Friedrich approximation the way we always
594 ! .. have
595 ! = 0 .. use DEFAULTS: 3 for RAD_FLAG = 0, 1 otherwise; note
596 ! that the defaults mimic the hardwired 'settings'
597 ! in SIMULATE, which doesnt read EXTRAD_FLAG
598 ! but determines what to do based on RAD_FLAG
599 !--------------------------------------------------------------
600
601 ! Check setting of EXTRAD_FLAG
602
603 if (debug(2)) write(6,*)'radc_init: entering...'
604 if (extrad_flag.eq.0) then
605 if (rad_flag.eq.0) then
606 extrad_flag = 3
607 else if (rad_flag.eq.1 .or. rad_flag.eq.2 .or. rad_flag.eq.3) then
608 extrad_flag = 1
609 endif
610 gaskelld 1.1 else if (extrad_flag.lt.0) then
611 stop 'Imbecile! check your stupid setting of EXTRAD_FLAG'
612 endif
613
614 ! 'etatzai' parameter
615
616 if (debug(4)) write(6,*)'radc_init: at 1'
617 etatzai = (12.0+(targ.Z+1.)/(targ.Z*targ.L1+targ.L2))/9.0
618 if (debug(4)) write(6,*)'radc_init: etatzai = ',etatzai
619
620 ! Initialize brem flags (brem doesn't include the normal common blocks)
621 produce_output = debug(1)
622 exponentiate = use_expon
623 include_hard = .true.
624 calculate_spence = .true.
625
626 if (debug(2)) write(6,*)'radc_init: ending...'
627 return
628 end
629
630 !---------------------------------------------------------------------
631 gaskelld 1.1
632 subroutine radc_init_ev (main,vertex)
633
634 implicit none
635 include 'structures.inc'
636 include 'radc.inc'
637
638 integer i
639 real*8 r, Ecutoff, dsoft, dhard, dsoft_prime
640 real*8 lambda_dave, schwinger, brem, bremos
641 record /event_main/ main
642 record /event/ vertex
643
644 ! Note that calculate_central calls this with (main,ev) rather than
645 ! (main,vertex). Since these are just local variables, calling it vertex
646 ! here does not cause any problem, and makes it easier to follow
647 ! modifications to vertex.* variables in later calls.
648
649 real*8 zero
650 parameter (zero=0.0d0) !double precision zero for subroutine calls.
651
652 gaskelld 1.1 ! Compute some quantities that will be needed for rad corr on this event
653
654 ! ... factor for limiting energy of external radiation along incident electron
655 ! etta = 1.0 + 2*vertex.ein*sin(vertex.e.theta/2.)**2/(targ.A*amu)
656 ! ... moron move! let's can that etta factor ...
657
658 etta = 1.0
659
660 ! ... the bt's
661
662 do i=1,2
663 bt(i) = etatzai*main.target.teff(i)
664 enddo
665
666 ! ... the lambda's (effective bt's for internal radiation)
667
668 do i=1,3
669 lambda(i) = lambda_dave(i,1,doing_tail(3),vertex.Ein,vertex.e.E,vertex.p.E,
670 > vertex.p.P,vertex.e.theta)
671 enddo
672 rad_proton_this_ev = lambda(3).gt.0
673 gaskelld 1.1
674 ! ... get the hard correction factor. don't care about Ecutoff! Just want dhard here
675
676 Ecutoff = 450.
677 if (intcor_mode.eq.0) then
678 r = schwinger(Ecutoff,vertex,.true.,dsoft,dhard)
679 else
680 if (.not.use_offshell_rad) then
681 r = brem(vertex.Ein,vertex.e.E,Ecutoff,rad_proton_this_ev,dsoft,dhard,
682 > dsoft_prime)
683 else
684 r = bremos(Ecutoff, zero, zero, vertex.Ein, vertex.e.P*vertex.ue.x,
685 > vertex.e.P*vertex.ue.y, vertex.e.P*vertex.ue.z, zero, zero, zero,
686 > vertex.p.P*vertex.up.x, vertex.p.P*vertex.up.y, vertex.p.P*vertex.up.z,
687 > vertex.p.E, rad_proton_this_ev, dsoft, dhard, dsoft_prime)
688 endif
689 endif
690 hardcorfac = 1./(1.-dhard)
691 g(4)=-dsoft_prime*Ecutoff+bt(1)+bt(2)
692
693 ! ... initialize the parameters needed for our "basic" calculation
694 gaskelld 1.1
695 call basicrad_init_ev (vertex.Ein,vertex.e.E,vertex.p.E)
696
697 ! ... the relative magnitudes of the three tails (we may not need them)
698
699 do i=1,3
700 frac(i) = g(i)/g(0)
701 enddo
702
703 return
704 end
705
706 !-----------------------------------------------------------------------
707
708 subroutine basicrad_init_ev (e1,e2,e3)
709
710 implicit none
711 include 'simulate.inc'
712 include 'radc.inc'
713
714 real*8 one
715 gaskelld 1.1 parameter (one=1.)
716
717 integer i
718 real*8 e1,e2,e3,e(3),gamma
719
720 if (debug(2)) write(6,*)'basicrad_init_ev: entering...'
721 e(1) = e1
722 e(2) = e2
723 e(3) = e3
724
725 ! bt's for internal + external
726 ! ??? One possibility for shutting off tails 1 or
727 ! 2 is to set g(1/2) = 0 here ... note that lambda(3) is set to 0 in
728 ! lambda_dave at the moment, AND proton terms are removed from brem if proton
729 ! radiation off. something analogous and similarly consistent would have to be
730 ! done for the other tails, right now they're just nixed in generate_rad. also
731 ! check ALL ntail.eq.0 checks in kinema constraint lines of generate_rad
732
733 g(1) = lambda(1) + bt(1)
734 g(2) = lambda(2) + bt(2)
735 g(3) = lambda(3)
736 gaskelld 1.1 g(0) = g(1)+g(2)+g(3)
737
738 ! Internal constants
739
740 c_int(1) = lambda(1)/(e(1)*e(2))**(lambda(1)/2.)
741 c_int(2) = lambda(2)/(e(1)*e(2))**(lambda(2)/2.)
742
743 * csa NOTE: GN says these masses s/b Mp!
744
745 c_int(3) = lambda(3)/(Mp*e(3))**(lambda(3)/2.)
746 do i = 1, 3
747 c_int(i) = c_int(i) * exp(-euler*lambda(i)) / gamma(one+lambda(i))
748 enddo
749 g_int = lambda(1) + lambda(2) + lambda(3)
750 c_int(0) = c_int(1)*c_int(2) * g_int / lambda(1)/lambda(2)
751
752 ! ... proton radiation could be off
753
754 if (lambda(3).gt.0) c_int(0) = c_int(0) * c_int(3)/lambda(3)
755 c_int(0) = c_int(0) * gamma(one+lambda(1)) * gamma(one+lambda(2))
756 > * gamma(one+lambda(3)) / gamma(one+g_int)
757 gaskelld 1.1
758 ! External constants
759
760 do i = 1, 2
761 c_ext(i) = bt(i)/e(i)**bt(i)/gamma(one+bt(i))
762 enddo
763 c_ext(3) = 0.0
764 g_ext = bt(1) + bt(2)
765 c_ext(0) = c_ext(1)*c_ext(2) * g_ext / bt(1)/bt(2)
766 c_ext(0) = c_ext(0)*gamma(one+bt(1))*gamma(one+bt(2))/gamma(one+g_ext)
767
768 ! Internal + external constants
769
770 do i = 1, 2
771 c(i) = c_int(i) * c_ext(i) * g(i)/lambda(i)/bt(i)
772 > * gamma(one+lambda(i))*gamma(one+bt(i))/gamma(one+g(i))
773 enddo
774 c(3) = c_int(3)
775
776 ! Finally, constant for combined tails
777
778 gaskelld 1.1 c(0) = c(1)*c(2) * g(0)/g(1)/g(2)
779
780 ! ... proton radiation could be off
781
782 if (g(3).gt.0) c(0) = c(0) * c(3)/g(3)
783 c(0)=c(0)*gamma(one+g(1))*gamma(one+g(2))*gamma(one+g(3))/gamma(one+g(0))
784 c(4)=g(4)/(e1*e2)**g(4)/gamma(one+g(4))
785 if(g(3).gt.0) c(4)=c(4)/e3**g(4)
786 if (debug(2)) write(6,*)'basicrad_init_ev: ending...'
787 return
788 end
789
790 !----------------------------------------------------------------------
791 ! theory_init:
792 !
793 ! Input spectral function(NOT called for H(e,e'p) or pion/kaon production!)
794 !
795 ! Note: Some flags that existed in old A(e,e'p) code (e.g. DD's version)
796 ! have been removed. Code has been changed to correspond to the following
797 ! settings for the old flags:
798 ! norm_Ji_tail = 0
799 gaskelld 1.1 ! doing_Ji_tail = 0
800 ! Em_start_Ji_tail = 0.
801 ! fixed_gamma = .true.
802
803 subroutine theory_init(success)
804
805 implicit none
806 include 'simulate.inc'
807
808 real*8 Pm_values(ntheorymax), absorption
809 integer m,n,iok
810 logical success
811
812 ! ... open the file
813 if ( nint(targ.A) .eq. 2) then
814 theory_file='h2.theory'
815 else if ( nint(targ.A) .eq. 12) then
816 theory_file='c12.theory'
817 else if ( nint(targ.A) .eq. 56) then
818 theory_file='fe56.theory'
819 else if ( nint(targ.A) .eq. 197) then
820 gaskelld 1.1 theory_file='au197.theory'
821 else
822 write(6,*) 'No Theory File (spectral function) for A = ',targ.A
823 write(6,*) 'Defaulting to c12.theory'
824 theory_file='c12.theory'
825 endif
826 open(unit=1,file=theory_file,status='old',readonly,shared,iostat=iok)
827
828 ! ... read in the theory file
829 read(1,*,err=40) nrhoPm, absorption, E_Fermi
830 do m=1, nrhoPm
831 read(1,*,err=40) nprot_theory(m), Em_theory(m), Emsig_theory(m),
832 > bs_norm_theory(m)
833 nprot_theory(m) = nprot_theory(m) * absorption
834 enddo
835 read(1,*,err=40) Pm_values(1),theory(1,1)
836 do m=1, nrhoPm
837 n=2
838 read(1,*,err=40,end=50) Pm_values(2),theory(m,2)
839 do while (Pm_values(n).gt.Pm_values(n-1))
840 n=n+1
841 gaskelld 1.1 read(1,*,err=40,end=50) Pm_values(n),theory(m,n)
842 enddo
843
844 ! ........ figure out details of the Pm axes
845 50 Pm_theory(m).n=n-1
846 Pm_theory(m).bin=Pm_values(2)-Pm_values(1)
847 Pm_theory(m).min=Pm_values(1)-Pm_theory(m).bin/2.
848 Pm_theory(m).max=Pm_values(Pm_theory(m).n)+Pm_theory(m).bin/2.
849 if (abs(Pm_theory(m).min+Pm_theory(m).bin*Pm_theory(m).n -
850 1 Pm_theory(m).max) .gt. 0.1) then
851 write(6,'(1x,''ERROR: theory_init found unequal Pm bins in distribution number '',i2,''!'')') m
852 close(1)
853 return
854 endif
855
856 ! ........ prepare for the next loop
857 Pm_values(1) = Pm_values(n)
858 theory(m+1,1) = theory(m,n)
859 enddo
860
861 ! ... are we doing deuterium? (i.e. only using a 1D spectral function)
862 gaskelld 1.1 doing_deuterium = nrhoPm.eq.1 .and. E_Fermi.lt.1.0
863
864 ! ... renormalize the momentum distributions
865 do m=1, nrhoPm
866 do n=1, Pm_theory(m).n
867 theory(m,n) = theory(m,n)/bs_norm_theory(m)
868 enddo
869 enddo
870
871 ! ... and calculate the integral of the Em distribution above E_Fermi to
872 ! ... renormalize
873 if (doing_heavy) then
874 do m=1, nrhoPm
875 Em_int_theory(m) = 1.
876 Em_int_theory(m) = (pi/2. + atan((Em_theory(m)-E_Fermi)/
877 > (0.5*Emsig_theory(m))))/pi
878 enddo
879 endif
880
881 ! ... we made it
882 success=.true.
883 gaskelld 1.1 close(1)
884 return
885
886 ! ... oops
887 40 continue
888 write(6,*) 'ERROR: theory_init failed to read in the theory file'
889
890 return
891 end
|
Return to
init.f
CVS log
Up to [HallC] / simc_semi