1 gaskelld 1.1 ! BREM.FOR
2
3 ! Exact soft photon bremstrahlung calculation of
4 ! delta, delta^\prime, and exponentiations thereof
5
6 real*8 function brem(ein,eout,egamma,radiate_proton,bsoft,bhard,dbsoft)
7
8 implicit none
9 include 'brem.inc'
10
11 real*8 pi, am, ame, e2
12 parameter (pi=3.141592653589793)
13 parameter (am= .93827231)
14 parameter (ame= .00051099906)
15 parameter (e2= 1./137.0359895)
16
17 real*8 ein ! electron energy
18 real*8 eout ! final electron energy
19 real*8 egamma ! energy of bremsstrahling photon
20 real*8 eang,ak,akp,ap,pang,eta,ape
21 real*8 q2,de
22 gaskelld 1.1 real*8 aprod,adot,ar1,ar2,alpha
23 real*8 bpi,bpf,bpp,bei,bef,bee,bepii,bepif,bepfi,bepff
24 real*8 dbpi,dbpf,dbpp,dbei,dbef,dbee,dbepii,dbepif,dbepfi,
25 > dbepff !derivatives of b's
26 real*8 b,bz,bzz,db,dbz,dbzz,bhard,bsch,spence,bsoft,dbsoft
27 real*8 inter
28 logical radiate_proton
29
30 ! Convert to GeV
31
32 ak=ein/1000.
33 akp=eout/1000.
34 de=egamma/1000.
35
36 ! electron scattering angle
37
38 eang= 2.*asin((am/(2.*ak)*(ak/akp-1.))**0.5)
39 eta= 1.+2.*ak*sin(eang/2.)**2/am
40 q2= 4.*ak*akp*sin(eang/2.)**2
41
42 ! final proton energy/momentum/scattering angle.
43 gaskelld 1.1
44 ape= am+ak-akp
45 ap= sqrt(ape**2-am**2)
46 pang= acos((ak-akp*cos(eang))/ap)
47
48 if(produce_output) then
49 write(6,*)' q2 ',q2
50 write(6,*)' k ',ak
51 write(6,*)' kp ',akp
52 write(6,*)' eang ',eang*180./pi
53 write(6,*)' ap ',ap
54 write(6,*)' pang ',pang*180./pi
55 write(6,*)' eta ',eta
56 endif
57
58 ! Calculate components of delta soft
59
60 ! ... electron terms
61
62 ! direct initial electron
63
64 gaskelld 1.1 aprod= 1.d0
65 bei= aprod*(-1./(2.*pi))*log(ak/de)
66 dbei= aprod*(1./(2.*pi*de))
67
68 ! direct final electron
69
70 aprod= 1.d0
71 bef= aprod*(-1./(2.*pi))*log(akp/de)
72 dbef= aprod*(1./(2.*pi*de))
73
74 ! e-e interference
75
76 aprod= -1.d0
77 adot= ak*akp*(1.-cos(eang))
78 alpha= 2.*ame**2-2.*adot
79 ar1= 0.5+sqrt(adot**2-ame**4)/alpha
80 ar2= 0.5-sqrt(adot**2-ame**4)/alpha
81 bee= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,ak,akp,de)
82 dbee= aprod*adot/(pi*alpha*(ar1-ar2)*de)*
83 > (log((ar1-1)/ar1)-log((ar2-1)/ar2))
84 if(produce_output) write(6,*) ar1,ar2
85 gaskelld 1.1
86 ! ... proton terms
87
88 if (radiate_proton) then
89
90 ! initial p direct
91
92 aprod= 1.d0
93 bpi= aprod*(-1./(2.*pi))*log(am/de)
94 dbpi= aprod*(1./(2.*pi*de))
95
96 ! final p direct
97
98 aprod= 1.d0
99 bpf= aprod*(-1./(2.*pi))*log(ape/de)
100 dbpf= aprod*(1/(2.*pi*de))
101
102 ! p-p interference
103
104 aprod= -1.d0
105 adot= am*ape
106 gaskelld 1.1 alpha= 2.*am**2-2.*adot
107 ar1= 0.5+sqrt(adot**2-am**4)/alpha
108 ar2= 0.5-sqrt(adot**2-am**4)/alpha
109 bpp= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,am,ape,de)
110 dbpp= aprod*adot/(pi*alpha*(ar1-ar2)*de)*
111 > (log((ar1-1)/ar1)-log((ar2-1)/ar2))
112 if(produce_output) write(6,*) ar1,ar2
113
114 ! ei-pi interference
115
116 aprod= -1.d0
117 adot= ak*am
118 alpha= am**2+ame**2-2.*adot
119 ar1= (am**2-adot+sqrt(adot**2-(ame*am)**2))/alpha
120 ar2= (am**2-adot-sqrt(adot**2-(ame*am)**2))/alpha
121 bepii= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,ak,am,de)
122 dbepii= aprod*adot/(pi*alpha*(ar1-ar2)*de)*
123 > (log((ar1-1)/ar1)-log((ar2-1)/ar2))
124 if(produce_output) write(6,*) ar1,ar2
125
126 ! ef-pf interference
127 gaskelld 1.1
128 aprod= -1.d0
129 adot= akp*ape-akp*ap*cos(eang+pang)
130 alpha= am**2+ame**2-2.*adot
131 ar1= (am**2-adot+sqrt(adot**2-(ame*am)**2))/alpha
132 ar2= (am**2-adot-sqrt(adot**2-(ame*am)**2))/alpha
133 bepff= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,akp,ape,de)
134 dbepff= aprod*adot/(pi*alpha*(ar1-ar2)*de)*
135 > (log((ar1-1)/ar1)-log((ar2-1)/ar2))
136 if(produce_output) write(6,*) ar1,ar2
137
138 ! ei-pf interference
139
140 aprod= 1.d0
141 adot= ak*ape-ak*ap*cos(pang)
142 alpha= am**2+ame**2-2.*adot
143 ar1= (am**2-adot+sqrt(adot**2-(ame*am)**2))/alpha
144 ar2= (am**2-adot-sqrt(adot**2-(ame*am)**2))/alpha
145 bepif= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,ak,ape,de)
146 dbepif= aprod*adot/(pi*alpha*(ar1-ar2)*de)*
147 > (log((ar1-1)/ar1)-log((ar2-1)/ar2))
148 gaskelld 1.1 if(produce_output) write(6,*) ar1,ar2
149
150 ! ef-pi interference
151
152 aprod= 1.d0
153 adot= akp*am
154 alpha= am**2+ame**2-2.*adot
155 ar1= (am**2-adot+sqrt(adot**2-(ame*am)**2))/alpha
156 ar2= (am**2-adot-sqrt(adot**2-(ame*am)**2))/alpha
157 bepfi= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,akp,am,de)
158 dbepfi= aprod*adot/(pi*alpha*(ar1-ar2)*de)*
159 > (log((ar1-1)/ar1)-log((ar2-1)/ar2))
160 if(produce_output) write(6,*) ar1,ar2
161
162 endif ! <radiate_proton>
163
164 ! All together now!
165 b= 2.*e2*(bei+bef+bee)
166 if (radiate_proton) then
167 bzz= 2.*e2*(bpi+bpf+bpp)
168 bz= 2.*e2*(bepii+bepff+bepif+bepfi)
169 gaskelld 1.1 else
170 bzz = 0.0
171 bz = 0.0
172 endif
173 bsoft=b+bz+bzz
174 bhard= -1.*(e2/pi)*(-28/9.+13./6.*log(q2/ame**2))
175 db= 2.*e2*(dbei+dbef+dbee)
176 if (radiate_proton) then
177 dbzz= 2.*e2*(dbpi+dbpf+dbpp)
178 dbz= 2.*e2*(dbepii+dbepff+dbepif+dbepfi)
179 else
180 dbzz = 0.0
181 dbz = 0.0
182 endif
183 dbsoft= db+dbz+dbzz
184
185 if(produce_output) then
186 write(6,*)' ----- results ----- '
187 write(6,*)' b= ',b
188 write(6,*)' bz= ',bz
189 write(6,*)' bzz= ',bzz
190 gaskelld 1.1 write(6,*)' bhard= ',bhard
191 write(6,*)' total= ',bsoft+bhard
192 write(6,*)' exp= ',1.-dexp(-1.*(bsoft))*(1.-bhard)
193 write(6,*)' '
194 write(6,*)' ultra-relativistic limit'
195 call srad(ak,akp,eang,q2,ame,am,ap,de,produce_output)
196 write(6,*)' Schwinger Result'
197 bsch= 2.*e2/pi*((log(ak/de)-13./12.)*(log(q2/ame**2)-1.)+17./36.+
198 > 0.5*(pi**2/6.-spence((cos(eang/2.))**2)))
199 write(6,*)' b= ',bsch
200 endif
201
202 ! ... and the result --> the value of the radiative cross-section,
203 ! dsigma/dEgamma = -dbsoft * exp(-bsoft) * (1-bhard)
204 ! ......... the derivative has dimension 1/[energy] --> convert back to MeV
205 dbsoft = dbsoft/1000.
206 if (exponentiate) then
207 brem = -dbsoft/dexp(bsoft)
208 else
209 brem = 1.-dbsoft
210 endif
211 gaskelld 1.1 if (include_hard) brem = brem*(1.-bhard)
212
213 return
214 end
215
216 !-----------------------------------
217
218 real*8 function inter(calculate_spence,alpha,ar1,ar2,e1,e2,de)
219
220 ! explicit evaluation of integral ... may or may not ignore spence functions
221
222 implicit none
223
224 real*8 pi
225 parameter (pi=3.141592653589793)
226
227 real*8 alpha,ar1,ar2,e1,e2,de
228 real*8 de2,amult,arg1,arg2,arg3,arg4,spence
229 logical calculate_spence
230
231 de2 = e1-e2
232 gaskelld 1.1 amult = -1./(alpha*(ar1-ar2))
233 inter = log(abs((e2/de)+ar1*(de2/de)))*log(abs((ar1-1.)/ar1)) -
234 > log(abs((e2/de)+ar2*(de2/de)))*log(abs((ar2-1.)/ar2))
235
236 if (calculate_spence) then
237 arg1= (de2/(e2+ar1*de2))*(ar1-1.)
238 arg2= (de2/(e2+ar1*de2))*(ar1)
239 arg3= (de2/(e2+ar2*de2))*(ar2-1.)
240 arg4= (de2/(e2+ar2*de2))*(ar2)
241 inter = inter - spence(arg1)+spence(arg2)+spence(arg3)-spence(arg4)
242 endif
243
244 inter= inter*amult/(pi)
245
246 return
247 end
248
249 !-----------------------------------
250
251 subroutine srad(ak,akp,eang,q2,ame,am,ap,de,produce_output)
252
253 gaskelld 1.1 ! calculates ultra-relativistic limit
254
255 implicit none
256
257 real*8 pi, alpha
258 parameter (pi = 3.141592653589793)
259 parameter (alpha = 1./137.0359895)
260
261 real*8 ak,akp,eang,q2,ame,am,ap
262 real*8 eta,dsoft,dhard,de
263 real*8 dsoftz,dsoftzz,ep
264 logical produce_output
265
266 ep= sqrt(ap**2+am**2)
267 dhard= -13./12.*(log(q2/ame**2)-1.)+17./36.
268 eta= 1.+2.*ak*(sin(eang/2.)**2)/am
269
270 dsoft= alpha/pi*log(ak*akp/de**2)*(log(q2/ame**2)-1.)
271 dsoftz= alpha/pi*log(eta)*(log(ak*akp/de**2)+log(am*ep/de**2))
272 dsoftzz= alpha/pi*log(am*ep/de**2)*(log(q2/am**2)-1.)
273 if(dsoftzz.le.0) dsoftzz=0.0
274 gaskelld 1.1
275 if (produce_output) then
276 write(6,*)' '
277 write(6,*)' b= ',dsoft
278 write(6,*)' bz= ',dsoftz
279 write(6,*)' bzz= ',dsoftzz
280 endif
281
282 return
283 end
284
285 !-----------------------------------
286
287 real*8 function spence(ax)
288
289 implicit none
290
291 real*8 pi
292 parameter (pi=3.141592653589793)
293
294 real*8 ax,bx
295 gaskelld 1.1
296 bx= dabs(ax)
297
298 ! ... N.B. Have replaced the former calculation (commented out) with an
299 ! ... approximate expression -- saves a WHALE of CPU!
300 !
301 ! if(bx.lt.1) then
302 ! spence= ssum(ax,100)
303 ! else if (ax.gt.1) then
304 ! spence= -0.5*(log(bx))**2+pi**2/3.-ssum(1./ax,100)
305 ! else if (ax.le.-1) then
306 ! spence= -0.5*(log(bx))**2-pi**2/6.-ssum(1./ax,100)
307 ! else if (ax.eq.1) then
308 ! spence= pi**2/6.
309 ! else if (ax.eq.-1) then
310 ! spence= -pi**2/12.
311 ! endif
312
313 if (bx.le.1) then
314 spence = 0.0
315 else
316 gaskelld 1.1 spence = -0.5*(log(bx))**2
317 endif
318
319 return
320 end
321
322 !-----------------------------------
323
324 real*8 function ssum(as,n)
325
326 implicit none
327
328 integer n,i
329 real*8 as
330
331 ssum= 0.0
332 if(as.ne.0) then
333 if(n.ge.10000) write(6,*)' large n in function ssum (brem.f)'
334 do i= 1,n
335 ssum= ssum+as**i/(1.*i)**2
336 enddo
337 gaskelld 1.1 endif
338
339 return
340 end
341
342 !------------------------------------------------------------------------------
343
344 real*8 function bremos(egamma, ! photon energy
345 > k_ix, k_iy, k_iz, ! incoming electron 3-momentum
346 > k_fx, k_fy, k_fz, ! scattered electron 3-momentum
347 > p_ix, p_iy, p_iz, ! incoming proton 3-momentum (pm)
348 > p_fx, p_fy, p_fz, p_fe, ! scattered proton 4-momentum
349 > radiate_proton, ! proton radiation flag
350 > bsoft, bhard, dbsoft)
351
352 ! Calculation of soft photon radiative correction factor and its derivative
353 ! allowing for both initial and final protons to be offshell.
354 ! conventions identical to on-shell calculation
355
356 implicit none
357 include 'brem.inc'
358 gaskelld 1.1
359 real*8 pi, twopi, ame, e2, mp
360 parameter (pi=3.141592653589793)
361 parameter (twopi=2.*pi)
362 parameter (ame= .00051099906)
363 parameter (e2= 1./137.0359895)
364 parameter (mp= .93827231)
365
366 type four_vector
367 real*8 e, x, y, z
368 end type
369
370 real*8 egamma, de
371 real*8 k_ix,k_iy,k_iz
372 real*8 k_fx,k_fy,k_fz
373 real*8 p_ix,p_iy,p_iz
374 real*8 p_fx,p_fy,p_fz,p_fe
375 real*8 ami,amf,q2
376 real*8 aprod,adot,ar1,ar2,alpha
377 real*8 bpi,bpf,bpp,bei,bef,bee,bepii,bepif,bepfi,bepff
378 real*8 dbpi,dbpf,dbpp,dbei,dbef,dbee,dbepii,dbepif,
379 gaskelld 1.1 > dbepfi,dbepff !derivatives of b's
380 real*8 b,bz,bzz,db,dbz,dbzz
381 real*8 bsoft,bhard,dbsoft,bsch
382 real*8 inter,inter_prime
383 logical radiate_proton
384 type(four_vector):: k_i, k_f, p_i, p_f
385
386 ! Initialize
387
388 ! ... put input into local variables, while converting energies/momenta to GeV.
389 de = egamma/1000.
390 k_i%x = k_ix/1000.
391 k_i%y = k_iy/1000.
392 k_i%z = k_iz/1000.
393 k_f%x = k_fx/1000.
394 k_f%y = k_fy/1000.
395 k_f%z = k_fz/1000.
396 p_i%x = p_ix/1000.
397 p_i%y = p_iy/1000.
398 p_i%z = p_iz/1000.
399 p_f%e = p_fe/1000.
400 gaskelld 1.1 p_f%x = p_fx/1000.
401 p_f%y = p_fy/1000.
402 p_f%z = p_fz/1000.
403
404 ! ... compute electron energies
405 k_i%e= (k_i%x**2+k_i%y**2+k_i%z**2+ame**2)**0.5
406 k_f%e= (k_f%x**2+k_f%y**2+k_f%z**2+ame**2)**0.5
407 ! ........ if he insists on e-m conservation like below, just compute p_i.e
408 c p_i%e = k_f%e+p_f%e-k_i%e
409 p_i%e = mp
410
411 ! ... check energy-momentum conservation
412 ! e_check= dabs(k_f%e+p_f%e-k_i%e-p_i%e)
413 ! x_check= dabs(k_f%x+p_f%x-k_i%x-p_i%x)
414 ! y_check= dabs(k_f%y+p_f%y-k_i%y-p_i%y)
415 ! z_check= dabs(k_f%z+p_f%z-k_i%z-p_i%z)
416 !
417 ! if((e_check.gt.0.0001).or.(x_check.gt.0.0001).or.
418 ! + (y_check.gt.0.0001).or.(z_check.gt.0.0001)) then
419 ! write(6,*)' bad kinematics'
420 ! return
421 gaskelld 1.1 ! endif
422
423 ! ... compute Q2 and masses of initial and final protons
424 q2=-1.*( (k_f%e-k_i%e)**2-(k_f%x-k_i%x)**2-
425 > (k_f%y-k_i%y)**2-(k_f%z-k_i%z)**2)
426 if (produce_output) write(6,*)' q2= ',q2
427 c ami= ((p_i%e)**2-(p_i%x)**2-(p_i%y)**2-(p_i%z)**2)**0.5
428 ami= mp
429 amf= ((p_f%e)**2-(p_f%x)**2-(p_f%y)**2-(p_f%z)**2)**0.5
430
431 ! Calculate components of delta soft
432
433 ! ... electron terms
434
435 ! ........ direct initial electron
436 aprod= 1.d0
437 bei= aprod*(-1./twopi)*log(k_i%e/de)
438 dbei= aprod*(-1./twopi)*(-1./de)
439
440 ! ........ direct final electron
441 aprod= 1.d0
442 gaskelld 1.1 bef= aprod*(-1./twopi)*log(k_f%e/de)
443 dbef= aprod*(-1./twopi)*(-1./de)
444
445 ! ........ e-e interference
446 aprod= -1.d0
447 adot= k_i%e*k_f%e-k_i%x*k_f%x-k_i%y*k_f%y-k_i%z*k_f%z
448 alpha= 2.*ame**2-2.*adot
449 ar1= 0.5+sqrt(4.*adot**2-4.*ame**4)/(2.*alpha)
450 ar2= 0.5-sqrt(4.*adot**2-4.*ame**4)/(2.*alpha)
451 bee= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,k_i%e,k_f%e,de)
452 dbee= aprod*adot*inter_prime(alpha,ar1,ar2,de)
453 if (produce_output) write(6,*) ar1,ar2
454
455 ! ... proton terms
456
457 if (radiate_proton) then
458
459 ! ........ initial p direct
460 aprod= 1.d0
461 bpi= aprod*(-1./twopi)*log(p_i%e/de)
462 dbpi= aprod*(-1./twopi)*(-1./de)
463 gaskelld 1.1
464 ! ........ final p direct
465 aprod= 1.d0
466 bpf= aprod*(-1./twopi)*log(p_f%e/de)
467 dbpf= aprod*(-1./twopi)*(-1./de)
468
469 ! ........ p-p interference
470 aprod= -1.d0
471 adot= p_i%e*p_f%e-p_i%x*p_f%x-p_i%y*p_f%y-p_i%z*p_f%z
472 alpha= ami**2+amf**2-2.*adot
473 ar1= (2.*amf**2-2.*adot+sqrt(4.*adot**2-4.*(ami*amf)**2))/(2.*alpha)
474 ar2= (2.*amf**2-2.*adot-sqrt(4.*adot**2-4.*(ami*amf)**2))/(2.*alpha)
475 bpp= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,p_i%e,p_f%e,de)
476 dbpp= aprod*adot*inter_prime(alpha,ar1,ar2,de)
477 if (produce_output) write(6,*) ar1,ar2
478
479 ! ........ ei-pi interference
480 aprod= -1.d0
481 adot= k_i%e*p_i%e-k_i%x*p_i%x-k_i%y*p_i%y-k_i%z*p_i%z
482 alpha= ami**2+ame**2-2.*adot
483 ar1= (2.*ami**2-2.*adot+sqrt(4.*adot**2-4.*(ame*ami)**2))/(2.*alpha)
484 gaskelld 1.1 ar2= (2.*ami**2-2.*adot-sqrt(4.*adot**2-4.*(ame*ami)**2))/(2.*alpha)
485 bepii= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,k_i%e,p_i%e,de)
486 dbepii= aprod*adot*inter_prime(alpha,ar1,ar2,de)
487 if (produce_output) write(6,*) ar1,ar2
488
489 ! ........ ef-pf interference
490 aprod= -1.d0
491 adot= k_f%e*p_f%e-k_f%x*p_f%x-k_f%y*p_f%y-k_f%z*p_f%z
492 alpha= amf**2+ame**2-2.*adot
493 ar1= (2.*amf**2-2.*adot+sqrt(4.*adot**2-4.*(ame*amf)**2))/(2.*alpha)
494 ar2= (2.*amf**2-2.*adot-sqrt(4.*adot**2-4.*(ame*amf)**2))/(2.*alpha)
495 bepff= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,k_f%e,p_f%e,de)
496 dbepff= aprod*adot*inter_prime(alpha,ar1,ar2,de)
497 if (produce_output) write(6,*) ar1,ar2
498
499 ! ........ ei-pf interference
500 aprod= 1.d0
501 adot= k_i%e*p_f%e-k_i%x*p_f%x-k_i%y*p_f%y-k_i%z*p_f%z
502 alpha= amf**2+ame**2-2.*adot
503 ar1=(2.*amf**2-2.*adot+sqrt(4.*adot**2-4.*(ame*amf)**2))/(2.*alpha)
504 ar2=(2.*amf**2-2.*adot-sqrt(4.*adot**2-4.*(ame*amf)**2))/(2.*alpha)
505 gaskelld 1.1 bepif= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,k_i%e,p_f%e,de)
506 dbepif= aprod*adot*inter_prime(alpha,ar1,ar2,de)
507 if (produce_output) write(6,*) ar1,ar2
508
509 ! ........ ef-pi interference
510 aprod= 1.d0
511 adot= k_f%e*p_i%e-k_f%x*p_i%x-k_f%y*p_i%y-k_f%z*p_i%z
512 alpha= ami**2+ame**2-2.*adot
513 ar1=(2.*ami**2-2.*adot+sqrt(4.*adot**2-4.*(ame*ami)**2))/(2.*alpha)
514 ar2=(2.*ami**2-2.*adot-sqrt(4.*adot**2-4.*(ame*ami)**2))/(2.*alpha)
515 bepfi= aprod*adot*inter(calculate_spence,alpha,ar1,ar2,k_f%e,p_i%e,de)
516 dbepfi= aprod*adot*inter_prime(alpha,ar1,ar2,de)
517 if (produce_output) write(6,*) ar1,ar2
518
519 endif ! <radiate_proton>
520
521 ! All together now!
522 b= 2.*e2*(bei+bef+bee)
523 if (radiate_proton) then
524 bzz= 2.*e2*(bpi+bpf+bpp)
525 bz= 2.*e2*(bepii+bepff+bepif+bepfi)
526 gaskelld 1.1 else
527 bzz = 0.0
528 bz = 0.0
529 endif
530 bsoft = b + bz + bzz
531 bhard= -1.*(e2/pi)*(-28/9.+13./6.*log(q2/ame**2))
532 * bhard= (e2/pi)*(2-1.5*log(q2/ame**2))
533 * bhard= bhard+vac(q2)
534 db= 2.*e2*(dbei+dbef+dbee)
535 if (radiate_proton) then
536 dbzz= 2.*e2*(dbpi+dbpf+dbpp)
537 dbz= 2.*e2*(dbepii+dbepff+dbepif+dbepfi)
538 else
539 dbzz = 0.0
540 dbz = 0.0
541 endif
542 dbsoft= db+dbz+dbzz
543
544 if (produce_output) then
545 write(6,*)' ----- results ----- '
546 write(6,*)' b+bhard= ',b+bhard
547 gaskelld 1.1 write(6,*)' bz= ',bz
548 write(6,*)' bzz= ',bzz
549 write(6,*)' bhard= ',bhard
550 write(6,*)' total= ',1-bsoft-bhard
551 write(6,*)' exp= ',dexp(-1.*(bsoft))*(1.-bhard)
552 write(6,*)' 1-bhard= ',1-bhard
553 write(6,*)' exps= ',dexp(-1.*(bsoft))
554 write(6,*)' expse= ',dexp(-1.*b)
555 write(6,*)' expsp= ',dexp(-1.*(bz+bzz))
556 write(6,*)' '
557 write(6,*)' Schwinger Result'
558 bsch= 2.*e2/pi*((log(k_i%e/de)-13./12.)*(log(q2/ame**2)-1.)+17./36.)
559 write(6,*)' b= ',bsch
560 endif
561
562 ! ... and the result --> the value of the radiative cross-section,
563 ! dsigma/dEgamma = -dbsoft * exp(-bsoft) * (1-bhard)
564 ! ......... the derivative has dimension 1/[energy] --> convert back to MeV
565 dbsoft = dbsoft/1000. !convert back to MeV
566 if (exponentiate) then
567 bremos = -dbsoft/dexp(bsoft)
568 gaskelld 1.1 else
569 bremos = 1.-dbsoft
570 endif
571 if (include_hard) bremos = bremos*(1.-bhard)
572
573 return
574 end
575
576 !-----------------------------------
577 c
578 c explicit evaluation of DERIVATIVE of integral (wrt de)
579 c
580
581 real*8 function inter_prime(alpha,ar1,ar2,de)
582
583 implicit none
584
585 real*8 pi
586 parameter (pi=3.141592653589793)
587
588 real*8 alpha,ar1,ar2,de
589 gaskelld 1.1 real*8 amult
590
591 amult = -1./(alpha*(ar1-ar2))
592 inter_prime = (-1./de)*amult/pi*
593 > (log(abs((ar1-1.)/ar1))-log(abs((ar2-1.)/ar2)))
594
595 return
596 end
597
598 !-----------------------------------
599
600 real*8 function vac(q2)
601
602 implicit none
603
604 real*8 q2
605 real*8 am(10),ae(10)
606 real*8 dele,dell
607 integer*4 j
608
609 real*8 del
610 gaskelld 1.1
611 am(1)= 0.00051099906
612 ae(1) = 1.0
613 am(2)= .1057
614 ae(2)= 1.0
615 am(3)= 1.782
616 ae(3)= 1.0
617 am(4)= 0.3
618 ae(4)= 2./3.
619 am(5)= .3
620 ae(5)= 1./3.
621 am(6)= .430
622 ae(6)= 1./3.
623 am(7)= 1.5
624 ae(7)= 2./3.
625 am(8)= 5.
626 ae(8)= 1./3.
627
628 dele= del(q2,am(1),ae(1))
629 dell= 0.0
630 do 20 j= 1,3
631 gaskelld 1.1 dell= dell+del(q2,am(j),ae(j))
632 20 continue
633
634 vac= dell
635
636 return
637 end
638
639 !-----------------------------------
640
641 real*8 function del(q2,bm,be)
642
643 implicit none
644 real*8 q2,bm,be
645 real*8 x,alpha
646
647 x= 4.*bm**2/q2
648 alpha= be**2/137.0359895
649 del= -2.*alpha/(3.*3.14159265)*(-5./3.+x+(1-x/2)*(1+x)**0.5*
650 > log(((1+x)**0.5+1)/((1+x)**0.5-1)))
651
652 gaskelld 1.1 return
653 end
|
Return to
brem.f
CVS log
Up to [HallC] / simc_gfortran