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