1 gaskelld 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 gaskelld 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 gaskelld 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 gaskelld 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 gaskelld 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 gaskelld 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 gaskelld 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 type(event_main):: main
134 type(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 gaskelld 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 gaskelld 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 gaskelld 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 gaskelld 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 cdg if (vertex.e.E .gt. edge.e.E.max) then
223 cdg force2=.true.
224 cdg ntail=2
225 cdg 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 gaskelld 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 gaskelld 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 gaskelld 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 .or. doing_rho .or.
288 > doing_semi) then
289 if (debug(5)) write(6,*)'gen_rad: at 2d'
290 Egamma_min(1) = 0.
291 Egamma_max(1) = gen%sumEgen%max - vertex%e%E
292 CDJG if (ntail.ne.0) Egamma_min(1) = gen.sumEgen.min - vertex.e.E
293 endif
294
295 gaskelld 1.1 ! ... Constrain Egamma<Egamma1_max(from radc_init), and add energy 'slop'
296 Egamma_max(1) = min(Egamma_max(1),Egamma1_max)
297 Egamma_min(1) = Egamma_min(1) - dE_edge_test
298 Egamma_max(1) = Egamma_max(1) + dE_edge_test
299
300 ! ... override Egamma_max with limit from dbase (if entered)
301 lim1=egamma_max(1)
302 if (hardwired_rad) then
303 Egamma_max(1) = Egamma_gen_max
304 if (lim1.gt.Egamma_gen_max) write(6,*)
305 > 'SIMC found an Egamma(1) limit above your hardwired value'
306 endif
307
308 ! ... radiate
309 if (debug(5)) write(6,*)'gen_rad: init Egamma_min,max =',
310 > Egamma_min(1),Egamma_max(1)
311 call basicrad(1*peaked_basis_flag,Egamma_min(1),Egamma_max(1),
312 > Egamma_used(1),basicrad_weight,basicrad_val_reciprocal)
313 if (basicrad_weight.le.0) then
314 if (debug(4)) write(6,*)'gen_rad: failed at 4'
315 return
316 gaskelld 1.1 endif
317 if (debug(5)) write(6,*)'Egamma_used',Egamma_used(1)
318
319 ! ... adjust kinematics and call complete_ev. If complete_ev fails, return.
320
321 vertex%Ein = vertex%Ein - Egamma_used(1)
322 if (debug(5)) write(6,*) 'calling complete_ev from radc'
323 call complete_ev(main,vertex,evsuccess)
324 if (.not.evsuccess) then !remove previous radiation, try again.
325 !! if (ntried.le.5000) write(6,'(1x,''COMPLETE_EV1 failed in GENERATE_RAD at event'',i7,''!'')') nevent
326 return
327 endif
328 if (debug(4)) write(6,*)'gen_rad: at 5'
329 rad_weight = rad_weight * basicrad_weight
330 if (rad_flag.le.1) then
331 rad_weight = peaked_rad_weight(vertex,Egamma_used(1),Egamma_min(1),
332 > Egamma_max(1),basicrad_val_reciprocal,basicrad_weight)
333 else
334 rad_weight=rad_weight*extrad_phi(1,vertex%Ein,vertex%e%E,Egamma_used(1))
335 endif
336 endif !end of tail 1 (incoming electron)
337 gaskelld 1.1
338 ! The VERTEX kinematics are now set, so check that we are inside the
339 ! VERTEXedges (SF limits).
340
341 if (doing_heavy) then
342 if (debug(5)) write(6,*)'gen_rad: Em, min, max =',vertex%Em,
343 > VERTEXedge%Em%min,VERTEXedge%Em%max
344 if (debug(5)) write(6,*)'gen_rad: Pm, min, max =',vertex%Pm,
345 > VERTEXedge%Pm%min,VERTEXedge%Pm%max
346 if (vertex%Em .lt. VERTEXedge%Em%min .or.
347 > vertex%Em .gt. VERTEXedge%Em%max .or.
348 > vertex%Pm .lt. VERTEXedge%Pm%min .or.
349 > vertex%Pm .gt. VERTEXedge%Pm%max) then
350 if (debug(4)) write(6,*)'gen_rad: failed at 6'
351 return
352 endif
353 endif
354
355 ! RADIATE TAIL #2: the scattered electron tail
356
357 if (doing_tail(2) .and. (ntail.eq.0.or.ntail.eq.2)) then
358 gaskelld 1.1 if (debug(5)) write(6,*) 'we are at 2'
359
360 ! ... Find min/max Egammas that will make in into electron acceptance
361
362 Egamma_min(2) = vertex%e%E - edge%e%E%max
363 Egamma_max(2) = vertex%e%E - edge%e%E%min
364
365 ! ... Measured Em must be within acceptance after radiation. Gives upper
366 ! ... limit on sum of photon energies, and lower limit if hadron arm does
367 ! ... not radiate. Em_MEAS = Em_VERT+(Egamma1+Egamma2+Egamma3)+/-delta_trec.
368
369 if (doing_eep) then
370 Egamma_max(2) = min(Egamma_max(2),
371 > (edge%Em%max - vertex%Em) - Egamma_used(1) + max_delta_Trec )
372 if (ntail.ne.0 .or. .not.rad_proton_this_ev)
373 > Egamma_min(2) = max(Egamma_min(2),
374 > (edge%Em%min - vertex%Em) - Egamma_used(1) - max_delta_Trec )
375 endif
376
377 Egamma_max(2) = min(Egamma_max(2),Egamma_tot_max-Egamma_used(1))
378 Egamma_min(2) = Egamma_min(2) - dE_edge_test
379 gaskelld 1.1 Egamma_max(2) = Egamma_max(2) + dE_edge_test
380
381 ! ... override Egamma_max with limit from dbase (if entered)
382 lim2=egamma_max(2)
383 if (hardwired_rad) then
384 Egamma_max(2) = Egamma_gen_max
385 if (lim2.gt.Egamma_gen_max) write(6,*)
386 > 'SIMC found an Egamma(2) limit above your hardwired value'
387 endif
388
389 ! ... radiate
390 call basicrad(2*peaked_basis_flag,Egamma_min(2),Egamma_max(2),
391 > Egamma_used(2),basicrad_weight,basicrad_val_reciprocal)
392 if (basicrad_weight.le.0) then
393 ! write(6,*)'gen_rad: Maybe you need to change Em limits!'
394 ! write(6,*)'min/max(2)=',egamma_min(2),egamma_max(2)
395 if (debug(4)) write(6,*)'gen_rad: failed at 7'
396 return
397 endif
398 if(debug(5))write(6,*)'gen_rad: peaked_basis_flag=',peaked_basis_flag
399 if(debug(5))write(6,*)'gen_rad: Egamma min,max,used=',
400 gaskelld 1.1 > Egamma_min(2),Egamma_max(2),Egamma_used(2)
401 if(debug(5))write(6,*)'gen_rad: basicrad_weight=',basicrad_weight
402 rad_weight = rad_weight * basicrad_weight
403 if (rad_flag.le.1) then
404 rad_weight = peaked_rad_weight(vertex,Egamma_used(2),Egamma_min(2),
405 > Egamma_max(2),basicrad_val_reciprocal,basicrad_weight)
406 else
407 rad_weight=rad_weight*extrad_phi(2,vertex%Ein,vertex%e%E,Egamma_used(2))
408 endif
409 endif !end of tail 2 (outgoing electron)
410
411 ! RADIATE TAIL #3: the scattered proton tail
412
413 if (rad_proton_this_ev .and. (ntail.eq.0.or.ntail.eq.3)) then
414 if (debug(5)) write(6,*) 'we are at 3'
415
416 ! ... Find min/max Egammas that will make in into hadron acceptance
417
418 Egamma_min(3) = vertex%p%E - edge%p%E%max
419 Egamma_max(3) = vertex%p%E - edge%p%E%min
420
421 gaskelld 1.1 ! ... Measured Em must be within acceptance after radiation. Gives upper
422 ! ... and lower limit on sum of photon energies (and thus on hadron radiation,
423 ! ... since other arms are done).
424 ! ... Em_MEAS = Em_VERT+(Egamma1+Egamma2+Egamma3)+/-delta_trec.
425
426 if (doing_eep) then
427 Egamma_max(3) = min(Egamma_max(3), (edge%Em%max - vertex%Em) -
428 > Egamma_used(1) - Egamma_used(2) + max_delta_Trec)
429 Egamma_min(3) = max(Egamma_min(3), (edge%Em%min - vertex%Em) -
430 > Egamma_used(1) - Egamma_used(2) - max_delta_Trec )
431 endif
432
433 Egamma_max(3) = min(Egamma_max(3),
434 > Egamma_tot_max-Egamma_used(1)-Egamma_used(2) )
435 Egamma_min(3) = Egamma_min(3) - dE_edge_test
436 Egamma_max(3) = Egamma_max(3) + dE_edge_test
437
438 ! ... override Egamma_max with limit from dbase (if entered)
439 lim3=egamma_max(3)
440 if (hardwired_rad) then
441 Egamma_max(3) = Egamma_gen_max
442 gaskelld 1.1 if (lim3.gt.Egamma_gen_max) write(6,*)
443 > 'SIMC found an Egamma(3 ) limit above your hardwired value'
444 endif
445
446 ! ... radiate
447 call basicrad(3*peaked_basis_flag,Egamma_min(3),Egamma_max(3),
448 > Egamma_used(3),basicrad_weight,basicrad_val_reciprocal)
449 if (basicrad_weight.le.0) then
450 if (debug(4)) write(6,*)'gen_rad: failed at 8'
451 return
452 endif
453 rad_weight = rad_weight * basicrad_weight
454 if (rad_flag.le.1) then
455 rad_weight = peaked_rad_weight(vertex,Egamma_used(3),Egamma_min(3),
456 > Egamma_max(3),basicrad_val_reciprocal,basicrad_weight)
457 else
458 rad_weight=rad_weight*extrad_phi(3,vertex%Ein,vertex%e%E,Egamma_used(3))
459 endif
460 endif !end of tail 3 (outgoing hadron)
461
462 ! Now obtain description of ORIG (orig = vertex + radiation) event.
463 gaskelld 1.1 ! ... Note that no kinematics have been changed yet, EXCEPT vertex.Ein, which
464 ! ... has been reduced by Egamma_used(1). orig.Ein is the initial (unradiated)
465 ! ... beam energy, so add back Egamma_used(1). Reduce orig.e.E and orig.p.E
466 ! ... here, and they will be the values passed to the monte carlo (for mult.
467 ! ... scattering and monte carlo model).
468
469 ! ... Remember that we have taken Ein(GEN)-Ebeam_vertex_ave into
470 ! ... account by shifting edge.Em)
471
472 c do i = 1, neventfields
473 c orig.all(i) = vertex.all(i)
474 c enddo
475 orig = vertex
476
477 if (debug(4)) write(6,*)'gen_rad: at 10: orig.p.E =',orig%p%E
478
479 ! ... remove radiation from incoming electron.
480
481 orig%Ein = vertex%Ein + Egamma_used(1)
482
483 ! ... adjust the e'/hadron momenta (if they radiated significantly)
484 gaskelld 1.1
485 orig%e%E = vertex%e%E - Egamma_used(2)
486 if(orig%e%E.le.0d0) then
487 if (debug(4)) write(6,*)'gen_rad: Negative electron energy -failed'
488 return
489 endif
490 orig%e%P = orig%e%E
491 orig%e%delta = (orig%e%P-spec%e%P)/spec%e%P*100.
492
493 orig%p%E = vertex%p%E - Egamma_used(3)
494 if(orig%p%E.le.Mh) then
495 if (debug(4)) write(6,*)'gen_rad: failed at 11'
496 return
497 endif
498 orig%p%P = sqrt(orig%p%E**2-Mh2)
499 orig%p%delta = (orig%p%P-spec%p%P)/spec%p%P*100.
500
501 ! ... Record some information for other routines.
502 phot1=Egamma_used(1)
503 phot2=Egamma_used(2)
504 phot3=Egamma_used(3)
505 gaskelld 1.1 ntup%radphot = phot1 + phot2 + phot3
506 ntup%radarm=ntail
507 if (force2) ntup%radarm=12
508
509 ! Complete determination of the event weight due to radiation.
510
511 main%gen_weight = main%gen_weight*rad_weight/hardcorfac
512 if (force2) main%gen_weight=main%gen_weight*frac(2)
513 if (debug(5)) write(6,*) 'mgw',main%gen_weight,rad_weight,hardcorfac
514 success = .true.
515
516 if (debug(2)) write(6,*)'gen_rad: exiting...'
517 return
518 end
519
520 !--------------------------------------------------------------------------
521
522 real*8 function peaked_rad_weight(vertex,Egamma,
523 > emin,emax,basicrad_val_reciprocal,basicrad_weight)
524
525 implicit none
526 gaskelld 1.1 include 'radc.inc'
527 include 'structures.inc'
528
529 real*8 Egamma, basicrad_val_reciprocal, basicrad_weight
530 real*8 t1, t2, r, phi_ext, ein, eout, emin, emax
531 real*8 brem, bremos, extrad_phi, gamma
532 real*8 dhard, dsoft_intmin, dsoft_intmax
533 real*8 dsoft_ext1, dsoft_ext2, dsoft_extmin
534 real*8 dsoft_int_primemin, dsoft_int_primemax
535 real*8 dsoft_ext1_prime, dsoft_ext2_prime
536 real*8 dsoft_ext_prime, dsoft_extmax, eul
537 real*8 dsoft_int,dsoft_ext,dsoft_int_prime
538 type(event):: vertex
539
540 real*8 zero
541 parameter (zero=0.0d0) !double precision zero for subroutine calls
542
543 basicrad_val_reciprocal=basicrad_val_reciprocal+0. !avoid unused variable error
544 ! Compute a more precise value for the radiative probability at the
545 ! selected Egamma, more precise than the BASICRAD distribution used to
546 ! generate Egamma that is.
547 gaskelld 1.1 !
548 ! NB: This subroutine only deals with calculations done in the PEAKING
549 ! APPROXIMATION basis -
550 ! i.e. we only get here if RAD_FLAG = 0 or 1
551 !
552 ! ... (sort of) KLUGE -- if Egamma < res limit, these things might blow up,
553 ! ... so just screw the correction and leave
554 ein=vertex%Ein
555 eout=vertex%e%E
556 eul=0.577215665
557 c if (Egamma.lt.Egamma_res_limit) then
558 c peaked_rad_weight = basicrad_weight
559 c return
560 c endif
561
562 ! ... External
563
564 phi_ext = 1.0
565 if (extrad_flag.le.2) then
566 phi_ext = extrad_phi(0,ein,eout,Egamma)
567 if (rad_flag.eq.1) then
568 gaskelld 1.1 peaked_rad_weight = basicrad_weight * phi_ext
569 return
570 endif
571 if(emin.gt.0)then
572 dsoft_extmin = log(g_ext/c_ext(0)/emin**g_ext)
573 endif
574 dsoft_extmax = log(g_ext/c_ext(0)/emax**g_ext)
575 dsoft_ext_prime = -g_ext/Egamma
576 else
577 t1 = bt(1)/etatzai
578 t2 = bt(2)/etatzai
579 call extrad_friedrich(ein,Egamma,t1,dsoft_ext1,dsoft_ext1_prime)
580 call extrad_friedrich(eout,Egamma,t2,dsoft_ext2,dsoft_ext2_prime)
581 dsoft_ext = dsoft_ext1 + dsoft_ext2
582 dsoft_ext_prime = dsoft_ext1_prime + dsoft_ext2_prime
583 endif
584
585 ! ... Internal
586 ! ........ use full BREM calculation of deltas
587
588 if (rad_flag.eq.0) then
589 gaskelld 1.1 if (.not.use_offshell_rad) then
590 if(emin.gt.0)then
591 r = brem(ein,eout,emin,rad_proton_this_ev,
592 > dsoft_intmin,dhard,dsoft_int_primemin)
593 else
594 dsoft_intmin=1.0
595 endif
596 r = brem(ein,eout,emax,rad_proton_this_ev,
597 > dsoft_intmax,dhard,dsoft_int_primemax)
598 else
599 if(emin.gt.0)then
600 r = bremos(emin, zero, zero, ein,
601 > vertex%e%E*vertex%ue%x, vertex%e%E*vertex%ue%y,
602 > vertex%e%E*vertex%ue%z,
603 c > vertex%Pmx, vertex%Pmy, vertex%Pmz,
604 > zero,zero,zero,
605 > vertex%p%P*vertex%up%x, vertex%p%P*vertex%up%y,
606 > vertex%p%P*vertex%up%z, vertex%p%E,
607 > rad_proton_this_ev,
608 > dsoft_intmin, dhard, dsoft_int_primemin)
609 else
610 gaskelld 1.1 dsoft_intmin=1.0
611 endif
612 r = bremos(emax, zero, zero, ein,
613 > vertex%e%E*vertex%ue%x, vertex%e%E*vertex%ue%y,
614 > vertex%e%E*vertex%ue%z,
615 c > vertex%Pmx, vertex%Pmy, vertex%Pmz
616 > zero,zero,zero,
617 > vertex%p%P*vertex%up%x, vertex%p%P*vertex%up%y,
618 > vertex%p%P*vertex%up%z, vertex%p%E,
619 > rad_proton_this_ev,
620 > dsoft_intmax, dhard, dsoft_int_primemax)
621 endif
622
623 ! ........ use basic calculation of internal deltas
624
625 else
626 dsoft_int = log(g_int/c_int(0)/Egamma**g_int)
627 dsoft_int_prime = -g_int/Egamma
628 endif
629
630 ! ... All together now
631 gaskelld 1.1
632 if(emin.gt.0)then
633 peaked_rad_weight = c_ext(0)/g_ext*(exp(-dsoft_intmax)*
634 > emax**g_ext-exp(-dsoft_intmin)*emin**g_ext)
635 else
636 peaked_rad_weight = c_ext(0)/g_ext*(exp(-dsoft_intmax)*emax**g_ext)
637 endif
638 peaked_rad_weight = peaked_rad_weight * exp(-eul*g(4))/gamma(1.+g(4))
639 > * gamma(1.+g(4)-bt(1)-bt(2))*gamma(1.+bt(1))
640 > * gamma(1.+bt(2))/gamma(1.+g(4))
641 if (peaked_rad_weight.lt.0) peaked_rad_weight = 0
642
643 return
644 end
645
646 !---------------------------------------------------------------------------
647
648 subroutine extrad_friedrich(Ei,Ecutoff,trad,dbrem,dbrem_prime)
649
650 implicit none
651 include 'simulate.inc'
652 gaskelld 1.1 include 'radc.inc'
653
654 real*8 Ei, Ecutoff, trad, x, dbrem, dbrem_prime
655
656 x = Ecutoff/Ei
657 dbrem = trad*(-(etatzai-0.5) - etatzai*log(x) + etatzai*x - 0.5*x**2)
658 dbrem_prime = -trad/Ei * (etatzai/x - etatzai + x)
659
660 return
661 end
662
663 !--------------------------------------------------------------------------
664
665 real*8 function extrad_phi(itail,E1,E2,Egamma)
666
667 implicit none
668 include 'radc.inc'
669
670 integer itail
671 real*8 E1, E2, Egamma, E(2), x, t, gamma
672
673 gaskelld 1.1 E(1) = E1
674 E(2) = E2
675
676 ! Compute the multiplicative external correction function PHI
677
678 ! ... CASE 1: phi = 1
679
680 extrad_phi = 1.0
681
682 ! ... CASE 2: phi correctly computed according to Dave's formulas
683
684 if (extrad_flag.eq.2) then
685 if (itail.eq.0) then
686 extrad_phi = 1. - (bt(1)/E(1)+bt(2)/E(2)) / (g(1)+g(2)) * Egamma
687 else if (itail.eq.1.or.itail.eq.2) then
688 extrad_phi = 1. - bt(itail)/E(itail) / g(itail) * Egamma
689 endif
690
691 ! ... CASE 3: phi computed in Friedrich prescription
692
693 else if (extrad_flag.eq.3) then
694 gaskelld 1.1 if (itail.eq.0) stop 'Idiot! a multiplicative factor EXTRAD_PHI is not defined for peaking approx and EXTRAD_FLAG>2!'
695 if (itail.eq.1.or.itail.eq.2) then
696 x = Egamma/E(itail)
697 t = bt(itail)/etatzai
698 extrad_phi = extrad_phi * (1. - x + x**2/etatzai) *
699 > exp(t*((etatzai-0.5)-etatzai*x+x**2/2.)) * gamma(1.+bt(itail))
700 endif
701 endif
702
703 return
704 end
705
706 !------------------------------------------------------------------------
707
708 real*8 function schwinger(Ecutoff,vertex,include_hard,dsoft,dhard)
709
710 implicit none
711 include 'simulate.inc'
712 include 'radc.inc'
713
714 real*8 Ecutoff,lq,s2,b,spence,dsoft,dhard,spen
715 gaskelld 1.1 logical include_hard
716 type(event):: vertex
717
718
719 lq = log(vertex%Q2/Me**2)-1.0
720 s2 = sin(vertex%e%theta/2.)**2
721 b = 1.+2.*vertex%nu*s2/(targ%A*amu)
722 spence = spen(1.-s2)-2.5893784
723 dsoft = alpi*lq*log( vertex%Ein/(etta**2)*vertex%e%E*b/(Ecutoff**2) )
724 dhard = -alpi*(2.166666*lq + spence - log(vertex%Ein/vertex%e%E)**2/2.0)
725
726 if (use_expon.eq.0) then
727 schwinger = exp(dsoft)
728 else
729 schwinger = 1.+dsoft
730 endif
731
732 if (include_hard) schwinger = schwinger/(1.-dhard)
733
734 return
735 end
736 gaskelld 1.1
737 !--------------------------------------------------------------------------
738
739 real*8 function spen(x)
740
741 implicit none
742
743 integer i
744 real*8 x,s,y
745
746 ! Approximation formula for the Spence-function or -integral according to
747 ! abramowitz and stegun : formula 27.7.2
748
749 y = 1.0
750 s = 0.0
751 i = 0
752 do while (i.le.100 .and. abs(y).gt.abs(s)*1.d-4)
753 i = i+1
754 y = x*y
755 s = s+y/float(i**2)
756 enddo
757 gaskelld 1.1 spen = s
758
759 return
760 end
761
762 !--------------------------------------------------------------------------
763
764 real*8 function lambda_dave(itail,plus_flag,doing_proton,e1,e2,e3,p3,th)
765
766 implicit none
767 include 'constants.inc'
768
769 integer itail,plus_flag
770 real*8 e1,e2,e3,p3,th
771 real*8 plus_term
772 logical doing_proton, warned/.false./
773
774 ! initialize to silence gfortran error
775 lambda_dave=0.0
776
777 ! The extended peaking approximation
778 gaskelld 1.1
779 plus_term = 0.0
780 if (plus_flag.eq.1 .and. itail.lt.3) then
781 plus_term = log((1.-cos(th))/2.)
782
783 ! ... only add in term due to ep interference if we're using proton
784 ! ... radiation
785
786 if (doing_proton) plus_term = plus_term + 2.*log(e1/e2)
787 endif
788
789 ! Compute lambdas
790
791 if (itail.eq.1) then
792 lambda_dave = alpi*(2.*log(2.*e1/Me) -1. + plus_term)
793 else if (itail.eq.2) then
794 lambda_dave = alpi*(2.*log(2.*e2/Me) -1. + plus_term)
795 else if (itail.eq.3) then
796 if (doing_proton) then
|
Return to
radc.f
CVS log
Up to [HallC] / simc_gfortran