1 jones 1.1 !--------------------------------------------------------------------------
2
3 subroutine basicrad(itail,Egamma_lo,Egamma_hi,Egamma,weight,val_reciprocal)
4
5 implicit none
6 include 'radc.inc'
7
8 integer itail
9 real*8 Egamma_lo, Egamma_hi, Egamma, weight, val_reciprocal
10 real*8 power_lo, power_hi
11 real*8 x, y, ymin
12
13 real*8 grnd
14
15 !--------------------------------------------------------------------------
16 !
17 ! Generate Egamma inside the requested region, using the function
18 ! g * Egamma**(g-1)
19 ! -------------------------------
20 ! (Egamma_max**g - Egamma_min**g)
21 ! which has the BASICRAD shape, is "invertible" as a probability
22 jones 1.1 ! distribution,
23 ! and integrates to 1 over the requested region { Egamma_min to Egamma_max }.
24 ! Each radiative event we generate with this then has to be weighted with
25 ! the _actual_ radiative tail form (which doesn't necessarily integrate to 1)
26 ! _divided_by_ this generating function, so we return the RECIPROCAL of the
27 ! function's value (at the selected Egamma) as VAL_RECIPROCAL. If we're not
28 ! using anything fancier than BASICRAD, we can immediately determine what
29 ! that
30 ! quotient is:
31 ! C
32 ! weight = --- * (Egamma_max**g - Egamma_min**g)
33 ! g
34 ! We return that as (you guessed it!) WEIGHT.
35 ! If we want to use BASICRAD x the phi function (correction to external tail
36 ! shape), all we need to do is multiply WEIGHT by PHI evaluated at the
37 ! selected photon energy. That's dont outside this routine, in GENERATE_RAD.
38 ! Note that if we're not using the BASICRAD prescription at all (except as
39 ! a generation tool), WEIGHT is not used again.
40 !
41 ! Note this routine handles the possibilities Egamma_lo < 0, Egamma_hi <
42 ! Egamma_lo, and Egamma_hi < 0
43 jones 1.1 !---------------------------------------------------------------------------
44
45 ! Initialize
46
47 Egamma = 0.0
48 weight = 0.0
49 val_reciprocal = 0.0
50
51 ! ... is radiation in this direction turned off?
52
53 if(itail.eq.0)itail=4
54 if (g(itail).le.0) then
55 weight = 1.0
56 return
57 endif
58
59 ! ... do we have a chance?
60
61 if (Egamma_hi.le.Egamma_lo .or. Egamma_hi.le.0) return
62
63 ! ... want to evaluate powers as INfrequently as possible!
64 jones 1.1
65 power_hi = Egamma_hi**g(itail)
66 power_lo = 0.0
67 if (Egamma_lo.gt.0) power_lo = Egamma_lo**g(itail)
68
69 ! ... obtain y region: from (Egamma_lo/Egamma_hi)**g to 1, generate flat y.
70
71 ymin = power_lo/power_hi
72 y = ymin + grnd()*(1.-ymin)
73 x = y**(1./g(itail))
74 Egamma = x*Egamma_hi
75
76 ! The value of our generating function at the selected point
77
78 if (Egamma.gt.0) val_reciprocal = Egamma**(1.-g(itail)) *
79 > (power_hi-power_lo) / g(itail)
80
81 ! Event weight = probability radiation falls inside requested range,
82 ! in full BASICRAD prescription
83
84 weight = c(itail)/g(itail) * (power_hi-power_lo)
85 jones 1.1
86 if(itail.eq.4)itail=0
87 return
88 end
89
90 !-------------------------------------------------------------------------
91
92 real*8 function gamma(x)
93
94 implicit none
95
96 integer i, n, s
97 real*8 x, y
98
99 ! Compute gamma function for xin
100 ! The series computes gamma(1+y), and must be fed y between 0 and 1
101 ! This can be accomplished using the relation
102 ! gamma(y+n+1) = (y+n)*gamma(y+n) = ... = (y+n)*...*(y+1)*gamma(y+1)
103
104 gamma = 1.0
105 n = nint((x-1)-0.5)
106 jones 1.1 y = x-1 - n
107 if (n.ne.0) then
108 s = sign(1,n)
109 do i = s, n, s
110 gamma = gamma*(y+1+i)**s
111 enddo
112 endif
113 gamma = gamma * (1. - 0.5748646*y + 0.9512363*y**2 -
114 > 0.6998588*y**3 + 0.4245549*y**4 - 0.1010678*y**5)
115 return
116 end
117
118 !--------------------------------------------------------------------------
119
120 subroutine generate_rad(main,vertex,orig,success)
121
122 implicit none
123 include 'simulate.inc'
124 include 'radc.inc'
125
126 integer i, peaked_basis_flag
127 jones 1.1 real*8 x, rad_weight
128 real*8 peaked_rad_weight, basicrad_val_reciprocal
129 real*8 basicrad_weight, extrad_phi
130 real*8 max_delta_Trec, ebeam_min, ebeam_max
131 logical force2 !force tail #2 to radiate.
132 logical evsuccess,success
133 record /event_main/ main
134 record /event/ vertex, orig
135
136 real*8 grnd
137
138 !---Generate Radiation-----------------------------------------------------
139 ! This routine radiates the incoming event
140 ! Here are comments from radc_init that describe rad_flag and extrad_flag:
141 !
142 ! First, about those (mysterious) 2 main 'radiative option' flags ...
143 !
144 ! The significance of RAD_FLAG:
145 ! RAD_FLAG = 0 .. use best available formulas, generate in
146 ! .. (ntail,Egamma) basis
147 ! = 1 .. use BASICRAD only, generate in (ntail,Egamma)
148 jones 1.1 ! .. basis
149 ! = 2 .. use BASICRAD only, generate in (Egamma(1,2,3))
150 ! .. basis but prevent overlap of tails (bogus, note)
151 ! = 3 .. use BASICRAD only, generate in (Egamma(1,2,3))
152 ! .. allowing radiation in all 3 directions
153 ! .. simultaneously
154 ! The (ntail,Egamma) basis can be called the PEAKED basis since it allows
155 ! only 3 photon directions. PEAKED_BASIS_FLAG is set to zero below when
156 ! the peaked basis is being used, in this way we can conveniently tell
157 ! the BASICRAD routine to use the full Egamma formula to generate the gamma
158 ! energy whenever it's called.
159 !
160 ! (See N. Makins' thesis, section 4.5.6, for more help on this point)
161 !
162 ! The significance of EXTRAD_FLAG:
163 ! EXTRAD_FLAG = 1 .. use only BASIC external radiation formulas
164 ! .. (phi = 1)
165 ! = 2 .. use BASIC ext rad formulas x phi
166 ! = 3 .. use Friedrich approximation the way we always
167 ! .. have
168 ! = 0 .. use DEFAULTS: 3 for RAD_FLAG = 0, 1 otherwise; note
169 jones 1.1 ! that the defaults mimic the hardwired 'settings'
170 ! in SIMULATE, which doesnt read EXTRAD_FLAG
171 ! but determines what to do based on RAD_FLAG
172 !---------------------------------------------------------------------------
173
174 ! Initialize.
175 ! ... ntup.radphot=sum energy of all radiated photons.
176 ! ... ntup.radarm=ntail, or 12 if tail 2 is forced.
177
178 lim1=0
179 lim2=0
180 lim3=0
181 phot1=0
182 phot2=0
183 phot3=0
184
185 if (debug(2)) write(6,*)'gen_rad: entering...'
186 success = .false.
187 rad_weight = 1
188 do i = 1,3
189 Egamma_used(i) = 0.0
190 jones 1.1 enddo
191 ntup.radphot=0.
192 ntup.radarm=0.
193 force2=.false.
194
195 ! Which tail has to be radiated? Set ntail=0 if all tails allowed to radiate.
196
197 peaked_basis_flag = 1
198 if (rad_flag.le.1) then
199 peaked_basis_flag = 0
200 x = grnd()
201 if (x.ge.frac(1)+frac(2)) then
202 ntail = 3
203 else if (x.ge.frac(1)) then
204 ntail = 2
205 else
206 ntail = 1
207 endif
208 else if (rad_flag.eq.2) then
209 ntail = int(grnd()*3.) + 1
210 if (ntail.eq.4) ntail=3
211 jones 1.1 else if (rad_flag.eq.3) then
212 ntail = 0
213 else
214 stop 'Idiot! rad_flag is set stupidly'
215 endif
216 if (debug(5)) write(6,*) 'vertexradc',vertex.e.E,edge.e.E.max
217
218 ! If the scattered electron energy is too high to be dectected, only events
219 ! where electron radiates will be accepted. Force ntail=2, and apply extra
220 ! weight to event equal to frac(2) (the normal probability that ntail=2)
221
222 if (vertex.e.E .gt. edge.e.E.max) then
223 force2=.true.
224 ntail=2
225 endif
226 ! Need maximum change in recoil momentum for some of the later limits (e,e'p).
227
228 max_delta_Trec = max((vertex.Trec - VERTEXedge.Trec.min),
229 > (VERTEXedge.Trec.max - vertex.Trec) )
230
231 ! RADIATE TAIL #1: the incident electron tail
232 jones 1.1
233 if (doing_tail(1) .and. (ntail.eq.0.or.ntail.eq.1)) then
234 if (debug(5)) write(6,*) 'gen_rad: at 1'
235
236 ! ... CASE 1: A(e,e'p) for A>2 - the effect of radiation on this tail alters
237 ! ... the VERTEX (actual) value of Em, so we want to make sure we're inside
238 ! ... any source spectral function cuts (like the Fermi surface).
239 ! ... Set Egamma_max to avoid getting below Em.min, and if Em is above Em.max,
240 ! ... set minimum to get below Em.max. Radiating a photon decreases Em
241 ! ... (relative to pre-radiated value), but also shifts Pm, and thus Trec.
242 ! ... Allow Trec shift as well as Em. Note that there are NO limits based on
243 ! ... spectrometer Em limits (cuts) because measured Em used Ebeam before rad.
244
245 ! ... If only the incoming electron arm radiates, then we can compare Em
246 ! ... to edge.Em.min (the SPECTROMETER Em limit) as well.
247
248 if (doing_heavy) then
249 if (debug(5)) write(6,*)'gen_rad: at 2a'
250 if (debug(5)) write(6,*)'vertex.Em=',vertex.Em
251 if (debug(5)) write(6,*)'VERTEXedge.Em.max=',VERTEXedge.Em.max
252 if (debug(5)) write(6,*)'max_delta_Trec=',max_delta_Trec
253 jones 1.1 Egamma_min(1)=vertex.Em - VERTEXedge.Em.max - max_delta_Trec
254 Egamma_max(1)=vertex.Em - VERTEXedge.Em.min + max_delta_Trec
255 if (ntail.ne.0) Egamma_max(1)=min(Egamma_max(1),
256 > vertex.Em - edge.Em.min + max_delta_Trec)
257
258 ! ... CASE 2: H(e,e'p) - Require that Ebeam after radiation can generate
259 ! ... electron inside of electron arm acceptance (Ein_max for Ee'_max, etc...).
260 ! ... Also require MEASURED Em < edge.Em.max (Em_meas=egamma1+egamma2+egamma3)
261
262 ! ... JRA: Note that we get ebeam max/min from E=mE'/(m+E'*(1-cos(theta)),
263 ! ... assuming that E'_max(min) gives E_max(min). However, for very large
264 ! ... angles, or large E'_max, you can get a negative value for E_max. When
265 ! ... this happens, just open up ebeam_max. In the long run, we need to get
266 ! ... a true E_max, not just E for E'_max.
267
268 else if (doing_hyd_elast) then
269 if (debug(5)) write(6,*)'gen_rad: at 2b'
270 ebeam_max = Mp*edge.e.E.max / (Mp-edge.e.E.max*(1.-vertex.ue.z))
271 if (ebeam_max.lt.0) ebeam_max=1.d10
272 ebeam_min = Mp*edge.e.E.min / (Mp-edge.e.E.min*(1.-vertex.ue.z))
273 Egamma_min(1) = vertex.Ein - ebeam_max
274 jones 1.1 Egamma_max(1) = vertex.Ein - ebeam_min
275 Egamma_max(1) = min(Egamma_max(1),edge.Em.max)
276
277 ! ... CASE 3: D(e,e'p) - limit radiation so that E_beam > E' at the vertex.
278
279 else if (doing_deuterium) then
280 Egamma_max(1) = min(Egamma1_max, gen.sumEgen.max - vertex.e.E)
281 if (ntail.ne.0) Egamma_min(1) = gen.sumEgen.min - vertex.e.E
282
283 ! ... CASE 4: Pion production. Too complicated to calculate min/max ebeam
284 ! ... from scattered particles. Use sumEgen-vertex.e.E as limit for remaining
285 ! ... generated energy available.
286
287 else if (doing_pion .or. doing_kaon) then
288 if (debug(5)) write(6,*)'gen_rad: at 2d'
289 Egamma_min(1) = 0.
290 Egamma_max(1) = gen.sumEgen.max - vertex.e.E
291 if (ntail.ne.0) Egamma_min(1) = gen.sumEgen.min - vertex.e.E
292 endif
293
294 ! ... Constrain Egamma<Egamma1_max(from radc_init), and add energy 'slop'
295 jones 1.1 Egamma_max(1) = min(Egamma_max(1),Egamma1_max)
296 Egamma_min(1) = Egamma_min(1) - dE_edge_test
297 Egamma_max(1) = Egamma_max(1) + dE_edge_test
298
299 ! ... override Egamma_max with limit from dbase (if entered)
300 lim1=egamma_max(1)
301 if (hardwired_rad) then
302 Egamma_max(1) = Egamma_gen_max
303 if (lim1.gt.Egamma_gen_max) write(6,*)
304 > 'SIMC found an Egamma(1) limit above your hardwired value'
305 endif
306
307 ! ... radiate
308 if (debug(5)) write(6,*)'gen_rad: init Egamma_min,max =',
309 > Egamma_min(1),Egamma_max(1)
310 call basicrad(1*peaked_basis_flag,Egamma_min(1),Egamma_max(1),
311 > Egamma_used(1),basicrad_weight,basicrad_val_reciprocal)
312 if (basicrad_weight.le.0) then
313 if (debug(4)) write(6,*)'gen_rad: failed at 4'
314 return
315 endif
316 jones 1.1 if (debug(5)) write(6,*)'Egamma_used',Egamma_used(1)
317
318 ! ... adjust kinematics and call complete_ev. If complete_ev fails, return.
319
320 vertex.Ein = vertex.Ein - Egamma_used(1)
321 if (debug(5)) write(6,*) 'calling complete_ev from radc'
322 call complete_ev(main,vertex,evsuccess)
323 if (.not.evsuccess) then !remove previous radiation, try again.
324 !! if (ntried.le.5000) write(6,'(1x,''COMPLETE_EV1 failed in GENERATE_RAD at event'',i7,''!'')') nevent
325 return
326 endif
327 if (debug(4)) write(6,*)'gen_rad: at 5'
328 rad_weight = rad_weight * basicrad_weight
329 if (rad_flag.le.1) then
330 rad_weight = peaked_rad_weight(vertex,Egamma_used(1),Egamma_min(1),
331 > Egamma_max(1),basicrad_val_reciprocal,basicrad_weight)
332 else
333 rad_weight=rad_weight*extrad_phi(1,vertex.Ein,vertex.e.E,Egamma_used(1))
334 endif
335 endif !end of tail 1 (incoming electron)
336
337 jones 1.1 ! The VERTEX kinematics are now set, so check that we are inside the
338 ! VERTEXedges (SF limits).
339
340 if (doing_heavy) then
341 if (debug(5)) write(6,*)'gen_rad: Em, min, max =',vertex.Em,
342 > VERTEXedge.Em.min,VERTEXedge.Em.max
343 if (debug(5)) write(6,*)'gen_rad: Pm, min, max =',vertex.Pm,
344 > VERTEXedge.Pm.min,VERTEXedge.Pm.max
345 if (vertex.Em .lt. VERTEXedge.Em.min .or.
346 > vertex.Em .gt. VERTEXedge.Em.max .or.
347 > vertex.Pm .lt. VERTEXedge.Pm.min .or.
348 > vertex.Pm .gt. VERTEXedge.Pm.max) then
349 if (debug(4)) write(6,*)'gen_rad: failed at 6'
350 return
351 endif
352 endif
353
354 ! RADIATE TAIL #2: the scattered electron tail
355
356 if (doing_tail(2) .and. (ntail.eq.0.or.ntail.eq.2)) then
357 if (debug(5)) write(6,*) 'we are at 2'
358 jones 1.1
359 ! ... Find min/max Egammas that will make in into electron acceptance
360
361 Egamma_min(2) = vertex.e.E - edge.e.E.max
362 Egamma_max(2) = vertex.e.E - edge.e.E.min
363
364 ! ... Measured Em must be within acceptance after radiation. Gives upper
365 ! ... limit on sum of photon energies, and lower limit if hadron arm does
366 ! ... not radiate. Em_MEAS = Em_VERT+(Egamma1+Egamma2+Egamma3)+/-delta_trec.
367
368 if (doing_eep) then
369 Egamma_max(2) = min(Egamma_max(2),
370 > (edge.Em.max - vertex.Em) - Egamma_used(1) + max_delta_Trec )
371 if (ntail.ne.0 .or. .not.rad_proton_this_ev)
372 > Egamma_min(2) = max(Egamma_min(2),
373 > (edge.Em.min - vertex.Em) - Egamma_used(1) - max_delta_Trec )
374 endif
375
376 Egamma_max(2) = min(Egamma_max(2),Egamma_tot_max-Egamma_used(1))
377 Egamma_min(2) = Egamma_min(2) - dE_edge_test
378 Egamma_max(2) = Egamma_max(2) + dE_edge_test
379 jones 1.1
380 ! ... override Egamma_max with limit from dbase (if entered)
381 lim2=egamma_max(2)
382 if (hardwired_rad) then
383 Egamma_max(2) = Egamma_gen_max
384 if (lim2.gt.Egamma_gen_max) write(6,*)
385 > 'SIMC found an Egamma(2) limit above your hardwired value'
386 endif
387
388 ! ... radiate
389 call basicrad(2*peaked_basis_flag,Egamma_min(2),Egamma_max(2),
390 > Egamma_used(2),basicrad_weight,basicrad_val_reciprocal)
391 if (basicrad_weight.le.0) then
392 ! write(6,*)'gen_rad: Maybe you need to change Em limits!'
393 ! write(6,*)'min/max(2)=',egamma_min(2),egamma_max(2)
394 if (debug(4)) write(6,*)'gen_rad: failed at 7'
395 return
396 endif
397 if(debug(5))write(6,*)'gen_rad: peaked_basis_flag=',peaked_basis_flag
398 if(debug(5))write(6,*)'gen_rad: Egamma min,max,used=',
399 > Egamma_min(2),Egamma_max(2),Egamma_used(2)
400 jones 1.1 if(debug(5))write(6,*)'gen_rad: basicrad_weight=',basicrad_weight
401 rad_weight = rad_weight * basicrad_weight
402 if (rad_flag.le.1) then
403 rad_weight = peaked_rad_weight(vertex,Egamma_used(2),Egamma_min(2),
404 > Egamma_max(2),basicrad_val_reciprocal,basicrad_weight)
405 else
406 rad_weight=rad_weight*extrad_phi(2,vertex.Ein,vertex.e.E,Egamma_used(2))
407 endif
408 endif !end of tail 2 (outgoing electron)
409
410 ! RADIATE TAIL #3: the scattered proton tail
411
412 if (rad_proton_this_ev .and. (ntail.eq.0.or.ntail.eq.3)) then
413 if (debug(5)) write(6,*) 'we are at 3'
414
415 ! ... Find min/max Egammas that will make in into hadron acceptance
416
417 Egamma_min(3) = vertex.p.E - edge.p.E.max
418 Egamma_max(3) = vertex.p.E - edge.p.E.min
419
420 ! ... Measured Em must be within acceptance after radiation. Gives upper
421 jones 1.1 ! ... and lower limit on sum of photon energies (and thus on hadron radiation,
422 ! ... since other arms are done).
423 ! ... Em_MEAS = Em_VERT+(Egamma1+Egamma2+Egamma3)+/-delta_trec.
424
425 if (doing_eep) then
426 Egamma_max(3) = min(Egamma_max(3), (edge.Em.max - vertex.Em) -
427 > Egamma_used(1) - Egamma_used(2) + max_delta_Trec)
428 Egamma_min(3) = max(Egamma_min(3), (edge.Em.min - vertex.Em) -
429 > Egamma_used(1) - Egamma_used(2) - max_delta_Trec )
430 endif
431
432 Egamma_max(3) = min(Egamma_max(3),
433 > Egamma_tot_max-Egamma_used(1)-Egamma_used(2) )
434 Egamma_min(3) = Egamma_min(3) - dE_edge_test
435 Egamma_max(3) = Egamma_max(3) + dE_edge_test
436
437 ! ... override Egamma_max with limit from dbase (if entered)
438 lim3=egamma_max(3)
439 if (hardwired_rad) then
440 Egamma_max(3) = Egamma_gen_max
441 if (lim3.gt.Egamma_gen_max) write(6,*)
442 jones 1.1 > 'SIMC found an Egamma(3 ) limit above your hardwired value'
443 endif
444
445 ! ... radiate
446 call basicrad(3*peaked_basis_flag,Egamma_min(3),Egamma_max(3),
447 > Egamma_used(3),basicrad_weight,basicrad_val_reciprocal)
448 if (basicrad_weight.le.0) then
449 if (debug(4)) write(6,*)'gen_rad: failed at 8'
450 return
451 endif
452 rad_weight = rad_weight * basicrad_weight
453 if (rad_flag.le.1) then
454 rad_weight = peaked_rad_weight(vertex,Egamma_used(3),Egamma_min(3),
455 > Egamma_max(3),basicrad_val_reciprocal,basicrad_weight)
456 else
457 rad_weight=rad_weight*extrad_phi(3,vertex.Ein,vertex.e.E,Egamma_used(3))
458 endif
459 endif !end of tail 3 (outgoing hadron)
460
461 ! Now obtain description of ORIG (orig = vertex + radiation) event.
462 ! ... Note that no kinematics have been changed yet, EXCEPT vertex.Ein, which
463 jones 1.1 ! ... has been reduced by Egamma_used(1). orig.Ein is the initial (unradiated)
464 ! ... beam energy, so add back Egamma_used(1). Reduce orig.e.E and orig.p.E
465 ! ... here, and they will be the values passed to the monte carlo (for mult.
466 ! ... scattering and monte carlo model).
467
468 ! ... Remember that we have taken Ein(GEN)-Ebeam_vertex_ave into
469 ! ... account by shifting edge.Em)
470
471 do i = 1, neventfields
472 orig.all(i) = vertex.all(i)
473 enddo
474 if (debug(4)) write(6,*)'gen_rad: at 10: orig.p.E =',orig.p.E
475
476 ! ... remove radiation from incoming electron.
477
478 orig.Ein = vertex.Ein + Egamma_used(1)
479
480 ! ... adjust the e'/hadron momenta (if they radiated significantly)
481
482 orig.e.E = vertex.e.E - Egamma_used(2)
483 orig.e.P = orig.e.E
484 jones 1.1 orig.e.delta = (orig.e.P-spec.e.P)/spec.e.P*100.
485
486 orig.p.E = vertex.p.E - Egamma_used(3)
487 if(orig.p.E.le.Mh) then
488 if (debug(4)) write(6,*)'gen_rad: failed at 11'
489 return
490 endif
491 orig.p.P = sqrt(orig.p.E**2-Mh2)
492 orig.p.delta = (orig.p.P-spec.p.P)/spec.p.P*100.
493
494 ! ... Record some information for other routines.
495 phot1=Egamma_used(1)
496 phot2=Egamma_used(2)
497 phot3=Egamma_used(3)
498 ntup.radphot = phot1 + phot2 + phot3
499 ntup.radarm=ntail
500 if (force2) ntup.radarm=12
501
502 ! Complete determination of the event weight due to radiation.
503
504 main.gen_weight = main.gen_weight*rad_weight/hardcorfac
505 jones 1.1 if (force2) main.gen_weight=main.gen_weight*frac(2)
506 if (debug(5)) write(6,*) 'mgw',main.gen_weight,rad_weight,hardcorfac
507 success = .true.
508
509 if (debug(2)) write(6,*)'gen_rad: exiting...'
510 return
511 end
512
513 !--------------------------------------------------------------------------
514
515 real*8 function peaked_rad_weight(vertex,Egamma,
516 > emin,emax,basicrad_val_reciprocal,basicrad_weight)
517
518 implicit none
519 include 'radc.inc'
520 include 'structures.inc'
521
522 real*8 Egamma, basicrad_val_reciprocal, basicrad_weight
523 real*8 t1, t2, r, phi_ext, ein, eout, emin, emax
524 real*8 brem, bremos, extrad_phi, gamma
525 real*8 dhard, dsoft_intmin, dsoft_intmax
526 jones 1.1 real*8 dsoft_ext1, dsoft_ext2, dsoft_extmin
527 real*8 dsoft_int_primemin, dsoft_int_primemax
528 real*8 dsoft_ext1_prime, dsoft_ext2_prime
529 real*8 dsoft_ext_prime, dsoft_extmax, eul
530 real*8 dsoft_int,dsoft_ext,dsoft_int_prime
531 record /event/ vertex
532
533 real*8 zero
534 parameter (zero=0.0d0) !double precision zero for subroutine calls
535
536 basicrad_val_reciprocal=basicrad_val_reciprocal+0. !avoid unused variable error
537 ! Compute a more precise value for the radiative probability at the
538 ! selected Egamma, more precise than the BASICRAD distribution used to
539 ! generate Egamma that is.
540 !
541 ! NB: This subroutine only deals with calculations done in the PEAKING
542 ! APPROXIMATION basis -
543 ! i.e. we only get here if RAD_FLAG = 0 or 1
544 !
545 ! ... (sort of) KLUGE -- if Egamma < res limit, these things might blow up,
546 ! ... so just screw the correction and leave
547 jones 1.1 ein=vertex.Ein
548 eout=vertex.e.E
549 eul=0.577215665
550 c if (Egamma.lt.Egamma_res_limit) then
551 c peaked_rad_weight = basicrad_weight
552 c return
553 c endif
554
555 ! ... External
556
557 phi_ext = 1.0
558 if (extrad_flag.le.2) then
559 phi_ext = extrad_phi(0,ein,eout,Egamma)
560 if (rad_flag.eq.1) then
561 peaked_rad_weight = basicrad_weight * phi_ext
562 return
563 endif
564 if(emin.gt.0)then
565 dsoft_extmin = log(g_ext/c_ext(0)/emin**g_ext)
566 endif
567 dsoft_extmax = log(g_ext/c_ext(0)/emax**g_ext)
568 jones 1.1 dsoft_ext_prime = -g_ext/Egamma
569 else
570 t1 = bt(1)/etatzai
571 t2 = bt(2)/etatzai
572 call extrad_friedrich(ein,Egamma,t1,dsoft_ext1,dsoft_ext1_prime)
573 call extrad_friedrich(eout,Egamma,t2,dsoft_ext2,dsoft_ext2_prime)
574 dsoft_ext = dsoft_ext1 + dsoft_ext2
575 dsoft_ext_prime = dsoft_ext1_prime + dsoft_ext2_prime
576 endif
577
578 ! ... Internal
579 ! ........ use full BREM calculation of deltas
580
581 if (rad_flag.eq.0) then
582 if (.not.use_offshell_rad) then
583 if(emin.gt.0)then
584 r = brem(ein,eout,emin,rad_proton_this_ev,
585 > dsoft_intmin,dhard,dsoft_int_primemin)
586 else
587 dsoft_intmin=1.0
588 endif
589 jones 1.1 r = brem(ein,eout,emax,rad_proton_this_ev,
590 > dsoft_intmax,dhard,dsoft_int_primemax)
591 else
592 if(emin.gt.0)then
593 r = bremos(emin, zero, zero, ein,
594 > vertex.e.E*vertex.ue.x, vertex.e.E*vertex.ue.y,
595 > vertex.e.E*vertex.ue.z,
596 c > vertex.Pmx, vertex.Pmy, vertex.Pmz,
597 > zero,zero,zero,
598 > vertex.p.P*vertex.up.x, vertex.p.P*vertex.up.y,
599 > vertex.p.P*vertex.up.z, vertex.p.E,
600 > rad_proton_this_ev,
601 > dsoft_intmin, dhard, dsoft_int_primemin)
602 else
603 dsoft_intmin=1.0
604 endif
605 r = bremos(emax, zero, zero, ein,
606 > vertex.e.E*vertex.ue.x, vertex.e.E*vertex.ue.y,
607 > vertex.e.E*vertex.ue.z,
608 c > vertex.Pmx, vertex.Pmy, vertex.Pmz
609 > zero,zero,zero,
610 jones 1.1 > vertex.p.P*vertex.up.x, vertex.p.P*vertex.up.y,
611 > vertex.p.P*vertex.up.z, vertex.p.E,
612 > rad_proton_this_ev,
613 > dsoft_intmax, dhard, dsoft_int_primemax)
614 endif
615
616 ! ........ use basic calculation of internal deltas
617
618 else
619 dsoft_int = log(g_int/c_int(0)/Egamma**g_int)
620 dsoft_int_prime = -g_int/Egamma
621 endif
622
623 ! ... All together now
624
625 if(emin.gt.0)then
626 peaked_rad_weight = c_ext(0)/g_ext*(exp(-dsoft_intmax)*
627 > emax**g_ext-exp(-dsoft_intmin)*emin**g_ext)
628 else
629 peaked_rad_weight = c_ext(0)/g_ext*(exp(-dsoft_intmax)*emax**g_ext)
630 endif
631 jones 1.1 peaked_rad_weight = peaked_rad_weight * exp(-eul*g(4))/gamma(1.+g(4))
632 > * gamma(1.+g(4)-bt(1)-bt(2))*gamma(1.+bt(1))
633 > * gamma(1.+bt(2))/gamma(1.+g(4))
634 if (peaked_rad_weight.lt.0) peaked_rad_weight = 0
635
636 return
637 end
638
639 !---------------------------------------------------------------------------
640
641 subroutine extrad_friedrich(Ei,Ecutoff,trad,dbrem,dbrem_prime)
642
643 implicit none
644 include 'simulate.inc'
645 include 'radc.inc'
646
647 real*8 Ei, Ecutoff, trad, x, dbrem, dbrem_prime
648
649 x = Ecutoff/Ei
650 dbrem = trad*(-(etatzai-0.5) - etatzai*log(x) + etatzai*x - 0.5*x**2)
651 dbrem_prime = -trad/Ei * (etatzai/x - etatzai + x)
652 jones 1.1
653 return
654 end
655
656 !--------------------------------------------------------------------------
657
658 real*8 function extrad_phi(itail,E1,E2,Egamma)
659
660 implicit none
661 include 'radc.inc'
662
663 integer itail
664 real*8 E1, E2, Egamma, E(2), x, t, gamma
665
666 E(1) = E1
667 E(2) = E2
668
669 ! Compute the multiplicative external correction function PHI
670
671 ! ... CASE 1: phi = 1
672
673 jones 1.1 extrad_phi = 1.0
674
675 ! ... CASE 2: phi correctly computed according to Dave's formulas
676
677 if (extrad_flag.eq.2) then
678 if (itail.eq.0) then
679 extrad_phi = 1. - (bt(1)/E(1)+bt(2)/E(2)) / (g(1)+g(2)) * Egamma
680 else if (itail.eq.1.or.itail.eq.2) then
681 extrad_phi = 1. - bt(itail)/E(itail) / g(itail) * Egamma
682 endif
683
684 ! ... CASE 3: phi computed in Friedrich prescription
685
686 else if (extrad_flag.eq.3) then
687 if (itail.eq.0) stop 'Idiot! a multiplicative factor EXTRAD_PHI is not defined for peaking approx and EXTRAD_FLAG>2!'
688 if (itail.eq.1.or.itail.eq.2) then
689 x = Egamma/E(itail)
690 t = bt(itail)/etatzai
691 extrad_phi = extrad_phi * (1. - x + x**2/etatzai) *
692 > exp(t*((etatzai-0.5)-etatzai*x+x**2/2.)) * gamma(1.+bt(itail))
693 endif
694 jones 1.1 endif
695
696 return
697 end
698
699 !------------------------------------------------------------------------
700
701 real*8 function schwinger(Ecutoff,vertex,include_hard,dsoft,dhard)
702
703 implicit none
704 include 'simulate.inc'
705 include 'radc.inc'
706
707 real*8 Ecutoff,lq,s2,b,spence,dsoft,dhard,spen
708 logical include_hard
709 record /event/ vertex
710
711
712 lq = log(vertex.Q2/Me**2)-1.0
713 s2 = sin(vertex.e.theta/2.)**2
714 b = 1.+2.*vertex.omega*s2/(targ.A*amu)
715 jones 1.1 spence = spen(1.-s2)-2.5893784
716 dsoft = alpi*lq*log( vertex.Ein/(etta**2)*vertex.e.E*b/(Ecutoff**2) )
717 dhard = -alpi*(2.166666*lq + spence - log(vertex.Ein/vertex.e.E)**2/2.0)
718
719 if (use_expon.eq.0) then
720 schwinger = exp(dsoft)
721 else
722 schwinger = 1.+dsoft
723 endif
724
725 if (include_hard) schwinger = schwinger/(1.-dhard)
726
727 return
728 end
729
730 !--------------------------------------------------------------------------
731
732 real*8 function spen(x)
733
734 implicit none
735
736 jones 1.1 integer i
737 real*8 x,s,y
738
739 ! Approximation formula for the Spence-function or -integral according to
740 ! abramowitz and stegun : formula 27.7.2
741
742 y = 1.0
743 s = 0.0
744 i = 0
745 do while (i.le.100 .and. abs(y).gt.abs(s)*1.d-4)
746 i = i+1
747 y = x*y
748 s = s+y/float(i**2)
749 enddo
750 spen = s
751
752 return
753 end
754
755 !--------------------------------------------------------------------------
756
757 jones 1.1 real*8 function lambda_dave(itail,plus_flag,doing_proton,e1,e2,e3,p3,th)
758
759 implicit none
760 include 'constants.inc'
761
762 integer itail,plus_flag
763 real*8 e1,e2,e3,p3,th
764 real*8 plus_term
765 logical doing_proton, warned/.false./
766
767 ! The extended peaking approximation
768
769 plus_term = 0.0
770 if (plus_flag.eq.1 .and. itail.lt.3) then
771 plus_term = log((1.-cos(th))/2.)
772
773 ! ... only add in term due to ep interference if we're using proton
774 ! ... radiation
775
776 if (doing_proton) plus_term = plus_term + 2.*log(e1/e2)
777 endif
778 jones 1.1
779 ! Compute lambdas
780
781 if (itail.eq.1) then
782 lambda_dave = alpi*(2.*log(2.*e1/Me) -1. + plus_term)
783 else if (itail.eq.2) then
784 lambda_dave = alpi*(2.*log(2.*e2/Me) -1. + plus_term)
785 else if (itail.eq.3) then
786 if (doing_proton) then
787 lambda_dave = alpi*(log((e3+p3)/(e3-p3)) - 2.)
788 ! ... at low energies, this may fritz, prevent human from causing itself
789 ! ... too much damage
790
791 if (lambda_dave.lt.0) then
792 lambda_dave = 0.0
793 if (.not.warned) then
794 write(6,'(1x,a)') '+-----------------------------------+'
795 write(6,'(1x,a)') '| WARNING: lambda_dave(3)<0, think |'
796 write(6,'(1x,a)') '| about running with proton rad OFF |'
797 write(6,'(1x,a)') '+-----------------------------------+'
798 warned = .true.
799 jones 1.1 endif
800 endif
801 else
802 lambda_dave = 0.0
803 endif
804 endif
805
806 return
807 end
|
Return to
radc.f
CVS log
Up to [HallC] / Poltar