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