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
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 type(cutstype):: the_phys, thp_phys, z, pp
97 type(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 > doing_hyddelta .or. doing_hydrho .or. doing_semi) 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 .or. doing_semi) 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 c gen%sumEgen%max = Ebeam_max - VERTEXedge%Trec%min - VERTEXedge%Trec_struck%min
438 c gen%sumEgen%min = Ebeam_min - VERTEXedge%Trec%max - VERTEXedge%Trec_struck%max
439 gen%sumEgen%max = Ebeam_max + targ%Mtar_struck - targ%Mrec_struck
440 gen%sumEgen%min = edge%e%E%min + edge%p%E%min
441
442 gaskelld 1.1 else !generated ELECTRON energy limits.
443 gen%sumEgen%max = Ebeam_max + targ%Mtar_struck - targ%Mrec_struck -
444 > edge%p%E%min - VERTEXedge%Em%min - VERTEXedge%Trec%min -
445 > VERTEXedge%Trec_struck%min
446 gen%sumEgen%min = Ebeam_min + targ%Mtar_struck - targ%Mrec_struck -
447 > edge%p%E%max - VERTEXedge%Em%max - VERTEXedge%Trec%max -
448 > VERTEXedge%Trec_struck%max - Egamma_tot_max
449 gen%sumEgen%max = min(gen%sumEgen%max, edge%e%E%max+Egamma2_max)
450 gen%sumEgen%min = max(gen%sumEgen%min, edge%e%E%min)
451 endif
452
453 gen%sumEgen%min = gen%sumEgen%min - dE_edge_test
454 gen%sumEgen%max = gen%sumEgen%max + dE_edge_test
455 gen%sumEgen%min = max(0.d0,gen%sumEgen%min)
456
457 ! ... E arm GENERATION limits from sumEgen.
458 ! ... Not used for doing_hyd_elast, but define for the hardwired histograms.
459
460 if (doing_hyd_elast) then
461 gen%e%E%min = edge%e%E%min
462 gen%e%E%max = edge%e%E%max + Egamma2_max
463 gaskelld 1.1 else if (doing_deuterium .or. doing_pion .or. doing_kaon
464 > .or. doing_rho .or. doing_delta) then
465 gen%e%E%min = gen%sumEgen%min
466 gen%e%E%max = gen%sumEgen%max
467 else if (doing_heavy .or. doing_semi) then
468 gen%e%E%min = gen%sumEgen%min - edge%p%E%max - Egamma3_max
469 gen%e%E%max = gen%sumEgen%max - edge%p%E%min
470 endif
471
472 ! ... Apply limits from direct comparison to acceptance.
473 gen%e%E%min = max(gen%e%E%min, edge%e%E%min)
474 gen%e%E%max = min(gen%e%E%max, edge%e%E%max + Egamma2_max)
475
476 gen%e%delta%min = (gen%e%E%min/spec%e%P-1.)*100.
477 gen%e%delta%max = (gen%e%E%max/spec%e%P-1.)*100.
478 gen%e%yptar%min = edge%e%yptar%min
479 gen%e%yptar%max = edge%e%yptar%max
480 gen%e%xptar%min = edge%e%xptar%min
481 gen%e%xptar%max = edge%e%xptar%max
482
483 ! ... P arm GENERATION limits from sumEgen. Not used for any case
484 gaskelld 1.1 ! ... except doing_heavy, but need to define for code that writes out limits.
485
486 if (doing_hyd_elast.or.doing_deuterium.or.doing_pion.or.doing_kaon .or.
487 > doing_rho .or. doing_delta) then
488 gen%p%E%min = edge%p%E%min
489 gen%p%E%max = edge%p%E%max + Egamma3_max
490 else if (doing_heavy .or. doing_semi)then
491 gen%p%E%min = gen%sumEgen%min - edge%e%E%max - Egamma2_max
492 gen%p%E%max = gen%sumEgen%max - edge%e%E%min
493 endif
494 ! ... Apply limits from direct comparison to acceptance.
495 gen%p%E%min = max(gen%p%E%min, edge%p%E%min)
496 gen%p%E%max = min(gen%p%E%max, edge%p%E%max + Egamma3_max)
497
498 gen%p%delta%min = (sqrt(gen%p%E%min**2-Mh2)/spec%p%P-1.)*100.
499 gen%p%delta%max = (sqrt(gen%p%E%max**2-Mh2)/spec%p%P-1.)*100.
500 gen%p%yptar%min = edge%p%yptar%min
501 gen%p%yptar%max = edge%p%yptar%max
502 gen%p%xptar%min = edge%p%xptar%min
503 gen%p%xptar%max = edge%p%xptar%max
504
505 gaskelld 1.1 ! Axis specs for the diagnostic histograms
506 H%gen%e%delta%min = gen%e%delta%min
507 H%gen%e%yptar%min = gen%e%yptar%min
508 H%gen%e%xptar%min = -gen%e%xptar%max
509 H%gen%e%delta%bin = (gen%e%delta%max-gen%e%delta%min)/float(nHbins)
510 H%gen%e%yptar%bin = (gen%e%yptar%max-gen%e%yptar%min)/float(nHbins)
511 H%gen%e%xptar%bin = (gen%e%xptar%max-gen%e%xptar%min)/float(nHbins)
512 H%gen%p%delta%min = gen%p%delta%min
513 H%gen%p%yptar%min = gen%p%yptar%min
514 H%gen%p%xptar%min = -gen%p%xptar%max
515 H%gen%p%delta%bin = (gen%p%delta%max-gen%p%delta%min)/float(nHbins)
516 H%gen%p%yptar%bin = (gen%p%yptar%max-gen%p%yptar%min)/float(nHbins)
517 H%gen%p%xptar%bin = (gen%p%xptar%max-gen%p%xptar%min)/float(nHbins)
518 H%gen%Em%min = VERTEXedge%Em%min
519 H%gen%Em%bin = (max(100.d0,VERTEXedge%Em%max) - VERTEXedge%Em%min)/float(nHbins)
520 H%gen%Pm%min = VERTEXedge%Pm%min
521 H%gen%Pm%bin = (max(100.d0,VERTEXedge%Pm%max) - VERTEXedge%Pm%min)/float(nHbins)
522
523 H%geni%e%delta%min = H%gen%e%delta%min
524 H%geni%e%yptar%min = H%gen%e%yptar%min
525 H%geni%e%xptar%min = H%gen%e%xptar%min
526 gaskelld 1.1 H%geni%e%delta%bin = H%gen%e%delta%bin
527 H%geni%e%yptar%bin = H%gen%e%yptar%bin
528 H%geni%e%xptar%bin = H%gen%e%xptar%bin
529 H%geni%p%delta%min = H%gen%p%delta%min
530 H%geni%p%yptar%min = H%gen%p%yptar%min
531 H%geni%p%xptar%min = H%gen%p%xptar%min
532 H%geni%p%delta%bin = H%gen%p%delta%bin
533 H%geni%p%yptar%bin = H%gen%p%yptar%bin
534 H%geni%p%xptar%bin = H%gen%p%xptar%bin
535 H%geni%Em%min = H%gen%Em%min
536 H%geni%Em%bin = H%gen%Em%bin
537 H%geni%Pm%min = H%gen%Pm%min
538 H%geni%Pm%bin = H%gen%Pm%bin
539
540 H%RECON%e%delta%min = H%gen%e%delta%min
541 H%RECON%e%yptar%min = H%gen%e%yptar%min
542 H%RECON%e%xptar%min = H%gen%e%xptar%min
543 H%RECON%e%delta%bin = H%gen%e%delta%bin
544 H%RECON%e%yptar%bin = H%gen%e%yptar%bin
545 H%RECON%e%xptar%bin = H%gen%e%xptar%bin
546 H%RECON%p%delta%min = H%gen%p%delta%min
547 gaskelld 1.1 H%RECON%p%yptar%min = H%gen%p%yptar%min
548 H%RECON%p%xptar%min = H%gen%p%xptar%min
549 H%RECON%p%delta%bin = H%gen%p%delta%bin
550 H%RECON%p%yptar%bin = H%gen%p%yptar%bin
551 H%RECON%p%xptar%bin = H%gen%p%xptar%bin
552 H%RECON%Em%min = H%gen%Em%min
553 H%RECON%Em%bin = H%gen%Em%bin
554 H%RECON%Pm%min = H%gen%Pm%min
555 H%RECON%Pm%bin = H%gen%Pm%bin
556
557 return
558 end
559
560 !------------------------------------------------------------------------
561
562 subroutine radc_init
563
564 implicit none
565 include 'simulate.inc'
566 include 'radc.inc'
567 include 'brem.inc'
568 gaskelld 1.1
569 !--------------------------------------------------------------
570 !
571 ! First, about those (mysterious) 2 main 'radiative option' flags ...
572 !
573 ! The significance of RAD_FLAG:
574 ! RAD_FLAG = 0 .. use best available formulas, generate in
575 ! .. (ntail,Egamma) basis
576 ! = 1 .. use BASICRAD only, generate in (ntail,Egamma)
577 ! .. basis
578 ! = 2 .. use BASICRAD only, generate in (Egamma(1,2,3))
579 ! .. basis but prevent overlap of tails (bogus, note)
580 ! = 3 .. use BASICRAD only, generate in (Egamma(1,2,3))
581 ! .. allowing radiation in all 3 directions
582 ! .. simultaneously
583 ! The (ntail,Egamma) basis can be called the PEAKED basis since it allows
584 ! only 3 photon directions. PEAKED_BASIS_FLAG is set to zero below when
585 ! the peaked basis is being used, in this way we can conveniently tell
586 ! the BASICRAD routine to use the full Egamma formula to generate the gamma
587 ! energy whenever it's called.
588 !
589 gaskelld 1.1 ! (See N. Makins' thesis, section 4.5.6, for more help on this point)
590 !
591 ! The significance of EXTRAD_FLAG:
592 ! EXTRAD_FLAG = 1 .. use only BASIC external radiation formulas
593 ! .. (phi = 1)
594 ! = 2 .. use BASIC ext rad formulas x phi
595 ! = 3 .. use Friedrich approximation the way we always
596 ! .. have
597 ! = 0 .. use DEFAULTS: 3 for RAD_FLAG = 0, 1 otherwise; note
598 ! that the defaults mimic the hardwired 'settings'
599 ! in SIMULATE, which doesnt read EXTRAD_FLAG
600 ! but determines what to do based on RAD_FLAG
601 !--------------------------------------------------------------
602
603 ! Check setting of EXTRAD_FLAG
604
605 if (debug(2)) write(6,*)'radc_init: entering...'
606 if (extrad_flag.eq.0) then
607 if (rad_flag.eq.0) then
608 extrad_flag = 3
609 else if (rad_flag.eq.1 .or. rad_flag.eq.2 .or. rad_flag.eq.3) then
610 gaskelld 1.1 extrad_flag = 1
611 endif
612 else if (extrad_flag.lt.0) then
613 stop 'Imbecile! check your stupid setting of EXTRAD_FLAG'
614 endif
615
616 ! 'etatzai' parameter
617
618 if (debug(4)) write(6,*)'radc_init: at 1'
619 etatzai = (12.0+(targ%Z+1.)/(targ%Z*targ%L1+targ%L2))/9.0
620 if (debug(4)) write(6,*)'radc_init: etatzai = ',etatzai
621
622 ! Initialize brem flags (brem doesn't include the normal common blocks)
623 produce_output = debug(1)
624 c exponentiate = use_expon
625 if(use_expon.eq.1) then
626 exponentiate=.true.
627 else
628 exponentiate=.false.
629 endif
630
631 gaskelld 1.1 include_hard = .true.
632 calculate_spence = .true.
633
634 if (debug(2)) write(6,*)'radc_init: ending...'
635 return
636 end
637
638 !---------------------------------------------------------------------
639
640 subroutine radc_init_ev (main,vertex)
641
642 implicit none
643 include 'structures.inc'
644 include 'radc.inc'
645
646 integer i
647 real*8 r, Ecutoff, dsoft, dhard, dsoft_prime
648 real*8 lambda_dave, schwinger, brem, bremos
649 type(event_main):: main
650 type(event):: vertex
651
652 gaskelld 1.1 ! Note that calculate_central calls this with (main,ev) rather than
653 ! (main,vertex). Since these are just local variables, calling it vertex
654 ! here does not cause any problem, and makes it easier to follow
655 ! modifications to vertex.* variables in later calls.
656
657 real*8 zero
658 parameter (zero=0.0d0) !double precision zero for subroutine calls.
659
660 ! Compute some quantities that will be needed for rad corr on this event
661
662 ! ... factor for limiting energy of external radiation along incident electron
663 ! etta = 1.0 + 2*vertex.ein*sin(vertex.e.theta/2.)**2/(targ%A*amu)
664 ! ... moron move! let's can that etta factor ...
665
666 etta = 1.0
667
668 ! ... the bt's
669
670 do i=1,2
671 bt(i) = etatzai*main%target%teff(i)
672 enddo
673 gaskelld 1.1
674 ! ... the lambda's (effective bt's for internal radiation)
675
676 do i=1,3
677 lambda(i) = lambda_dave(i,1,doing_tail(3),vertex%Ein,vertex%e%E,vertex%p%E,
678 > vertex%p%P,vertex%e%theta)
679 enddo
680 rad_proton_this_ev = lambda(3).gt.0
681
682 ! ... get the hard correction factor. don't care about Ecutoff! Just want dhard here
683
684 Ecutoff = 450.
685 if (intcor_mode.eq.0) then
686 r = schwinger(Ecutoff,vertex,.true.,dsoft,dhard)
687 else
688 if (.not.use_offshell_rad) then
689 r = brem(vertex%Ein,vertex%e%E,Ecutoff,rad_proton_this_ev,dsoft,dhard,
690 > dsoft_prime)
691 else
692 r = bremos(Ecutoff, zero, zero, vertex%Ein, vertex%e%P*vertex%ue%x,
693 > vertex%e%P*vertex%ue%y, vertex%e%P*vertex%ue%z, zero, zero, zero,
694 gaskelld 1.1 > vertex%p%P*vertex%up%x, vertex%p%P*vertex%up%y, vertex%p%P*vertex%up%z,
695 > vertex%p%E, rad_proton_this_ev, dsoft, dhard, dsoft_prime)
696 endif
697 endif
698 hardcorfac = 1./(1.-dhard)
699 g(4)=-dsoft_prime*Ecutoff+bt(1)+bt(2)
700
701 ! ... initialize the parameters needed for our "basic" calculation
702
703 call basicrad_init_ev (vertex%Ein,vertex%e%E,vertex%p%E)
704
705 ! ... the relative magnitudes of the three tails (we may not need them)
706
707 do i=1,3
708 frac(i) = g(i)/g(0)
709 enddo
710
711 return
712 end
713
714 !-----------------------------------------------------------------------
715 gaskelld 1.1
716 subroutine basicrad_init_ev (e1,e2,e3)
717
718 implicit none
719 include 'simulate.inc'
720 include 'radc.inc'
721
722 real*8 one
723 parameter (one=1.)
724
725 integer i
726 real*8 e1,e2,e3,e(3),gamma
727
728 if (debug(2)) write(6,*)'basicrad_init_ev: entering...'
729 e(1) = e1
730 e(2) = e2
731 e(3) = e3
732
733 ! bt's for internal + external
734 ! ??? One possibility for shutting off tails 1 or
735 ! 2 is to set g(1/2) = 0 here ... note that lambda(3) is set to 0 in
736 gaskelld 1.1 ! lambda_dave at the moment, AND proton terms are removed from brem if proton
737 ! radiation off. something analogous and similarly consistent would have to be
738 ! done for the other tails, right now they're just nixed in generate_rad. also
739 ! check ALL ntail.eq.0 checks in kinema constraint lines of generate_rad
740
741 g(1) = lambda(1) + bt(1)
742 g(2) = lambda(2) + bt(2)
743 g(3) = lambda(3)
744 g(0) = g(1)+g(2)+g(3)
745
746 ! Internal constants
747
748 c_int(1) = lambda(1)/(e(1)*e(2))**(lambda(1)/2.)
749 c_int(2) = lambda(2)/(e(1)*e(2))**(lambda(2)/2.)
750
751 * csa NOTE: GN says these masses s/b Mp!
752
753 c_int(3) = lambda(3)/(Mp*e(3))**(lambda(3)/2.)
754 do i = 1, 3
755 c_int(i) = c_int(i) * exp(-euler*lambda(i)) / gamma(one+lambda(i))
756 enddo
757 gaskelld 1.1 g_int = lambda(1) + lambda(2) + lambda(3)
758 c_int(0) = c_int(1)*c_int(2) * g_int / lambda(1)/lambda(2)
759
760 ! ... proton radiation could be off
761
762 if (lambda(3).gt.0) c_int(0) = c_int(0) * c_int(3)/lambda(3)
763 c_int(0) = c_int(0) * gamma(one+lambda(1)) * gamma(one+lambda(2))
764 > * gamma(one+lambda(3)) / gamma(one+g_int)
765
766 ! External constants
767
768 do i = 1, 2
769 c_ext(i) = bt(i)/e(i)**bt(i)/gamma(one+bt(i))
770 enddo
771 c_ext(3) = 0.0
772 g_ext = bt(1) + bt(2)
773 c_ext(0) = c_ext(1)*c_ext(2) * g_ext / bt(1)/bt(2)
774 c_ext(0) = c_ext(0)*gamma(one+bt(1))*gamma(one+bt(2))/gamma(one+g_ext)
775
776 ! Internal + external constants
777
778 gaskelld 1.1 do i = 1, 2
779 c(i) = c_int(i) * c_ext(i) * g(i)/lambda(i)/bt(i)
780 > * gamma(one+lambda(i))*gamma(one+bt(i))/gamma(one+g(i))
781 enddo
782 c(3) = c_int(3)
783
784 ! Finally, constant for combined tails
785
786 c(0) = c(1)*c(2) * g(0)/g(1)/g(2)
787
788 ! ... proton radiation could be off
789
790 if (g(3).gt.0) c(0) = c(0) * c(3)/g(3)
791 c(0)=c(0)*gamma(one+g(1))*gamma(one+g(2))*gamma(one+g(3))/gamma(one+g(0))
792 c(4)=g(4)/(e1*e2)**g(4)/gamma(one+g(4))
793 if(g(3).gt.0) c(4)=c(4)/e3**g(4)
794 if (debug(2)) write(6,*)'basicrad_init_ev: ending...'
795 return
796 end
797
798 !----------------------------------------------------------------------
799 gaskelld 1.1 ! theory_init:
800 !
801 ! Input spectral function(NOT called for H(e,e'p) or pion/kaon production!)
802 !
803 ! Note: Some flags that existed in old A(e,e'p) code (e.g. DD's version)
804 ! have been removed. Code has been changed to correspond to the following
805 ! settings for the old flags:
806 ! norm_Ji_tail = 0
807 ! doing_Ji_tail = 0
808 ! Em_start_Ji_tail = 0.
809 ! fixed_gamma = .true.
810
811 subroutine theory_init(success)
812
813 implicit none
814 include 'simulate.inc'
815
816 real*8 Pm_values(ntheorymax), absorption
817 integer m,n,iok
818 logical success
819
820 gaskelld 1.1 ! ... open the file
821 if ( nint(targ%A) .eq. 2) then
822 theory_file='h2.theory'
823 else if ( nint(targ%A) .eq. 12) then
824 theory_file='c12.theory'
825 else if ( nint(targ%A) .eq. 56) then
826 theory_file='fe56.theory'
827 else if ( nint(targ%A) .eq. 197) then
828 theory_file='au197.theory'
829 else
830 write(6,*) 'No Theory File (spectral function) for A = ',targ%A
831 write(6,*) 'Defaulting to c12.theory'
832 theory_file='c12.theory'
833 endif
834 c open(unit=1,file=theory_file,status='old',readonly,shared,iostat=iok)
835 open(unit=1,file=theory_file,status='old',iostat=iok)
836
837 ! ... read in the theory file
838 read(1,*,err=40) nrhoPm, absorption, E_Fermi
839 do m=1, nrhoPm
840 read(1,*,err=40) nprot_theory(m), Em_theory(m), Emsig_theory(m),
841 gaskelld 1.1 > bs_norm_theory(m)
842 nprot_theory(m) = nprot_theory(m) * absorption
843 enddo
844 read(1,*,err=40) Pm_values(1),theory(1,1)
845 do m=1, nrhoPm
846 n=2
847 read(1,*,err=40,end=50) Pm_values(2),theory(m,2)
848 do while (Pm_values(n).gt.Pm_values(n-1))
849 n=n+1
850 read(1,*,err=40,end=50) Pm_values(n),theory(m,n)
851 enddo
852
853 ! ........ figure out details of the Pm axes
854 50 Pm_theory(m)%n=n-1
855 Pm_theory(m)%bin=Pm_values(2)-Pm_values(1)
856 Pm_theory(m)%min=Pm_values(1)-Pm_theory(m)%bin/2.
857 Pm_theory(m)%max=Pm_values(Pm_theory(m)%n)+Pm_theory(m)%bin/2.
858 if (abs(Pm_theory(m)%min+Pm_theory(m)%bin*Pm_theory(m)%n -
859 > Pm_theory(m)%max) .gt. 0.1) then
860 write(6,'(1x,''ERROR: theory_init found unequal Pm bins in distribution number '',i2,''!'')') m
861 close(1)
862 gaskelld 1.1 return
863 endif
864
865 ! ........ prepare for the next loop
866 Pm_values(1) = Pm_values(n)
867 theory(m+1,1) = theory(m,n)
868 enddo
869
870 ! ... are we doing deuterium? (i.e. only using a 1D spectral function)
871 doing_deuterium = nrhoPm.eq.1 .and. E_Fermi.lt.1.0
872
873 ! ... renormalize the momentum distributions
874 do m=1, nrhoPm
875 do n=1, Pm_theory(m)%n
876 theory(m,n) = theory(m,n)/bs_norm_theory(m)
877 enddo
878 enddo
879
880 ! ... and calculate the integral of the Em distribution above E_Fermi to
881 ! ... renormalize
882 if (doing_heavy) then
883 gaskelld 1.1 do m=1, nrhoPm
884 Em_int_theory(m) = 1.
885 Em_int_theory(m) = (pi/2. + atan((Em_theory(m)-E_Fermi)/
886 > (0.5*Emsig_theory(m))))/pi
887 enddo
888 endif
889
890 ! ... we made it
891 success=.true.
892 close(1)
893 return
894
895 ! ... oops
896 40 continue
897 write(6,*) 'ERROR: theory_init failed to read in the theory file'
898
899 return
900 end
|
Return to
init.f
CVS log
Up to [HallC] / simc_gfortran