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