!-------------------------------------------------------------------------- subroutine basicrad(itail,Egamma_lo,Egamma_hi,Egamma,weight,val_reciprocal) implicit none include 'radc.inc' integer itail real*8 Egamma_lo, Egamma_hi, Egamma, weight, val_reciprocal real*8 power_lo, power_hi real*8 x, y, ymin real*8 grnd !-------------------------------------------------------------------------- ! ! Generate Egamma inside the requested region, using the function ! g * Egamma**(g-1) ! ------------------------------- ! (Egamma_max**g - Egamma_min**g) ! which has the BASICRAD shape, is "invertible" as a probability ! distribution, ! and integrates to 1 over the requested region { Egamma_min to Egamma_max }. ! Each radiative event we generate with this then has to be weighted with ! the _actual_ radiative tail form (which doesn't necessarily integrate to 1) ! _divided_by_ this generating function, so we return the RECIPROCAL of the ! function's value (at the selected Egamma) as VAL_RECIPROCAL. If we're not ! using anything fancier than BASICRAD, we can immediately determine what ! that ! quotient is: ! C ! weight = --- * (Egamma_max**g - Egamma_min**g) ! g ! We return that as (you guessed it!) WEIGHT. ! If we want to use BASICRAD x the phi function (correction to external tail ! shape), all we need to do is multiply WEIGHT by PHI evaluated at the ! selected photon energy. That's dont outside this routine, in GENERATE_RAD. ! Note that if we're not using the BASICRAD prescription at all (except as ! a generation tool), WEIGHT is not used again. ! ! Note this routine handles the possibilities Egamma_lo < 0, Egamma_hi < ! Egamma_lo, and Egamma_hi < 0 !--------------------------------------------------------------------------- ! Initialize Egamma = 0.0 weight = 0.0 val_reciprocal = 0.0 ! ... is radiation in this direction turned off? if(itail.eq.0)itail=4 if (g(itail).le.0) then weight = 1.0 return endif ! ... do we have a chance? if (Egamma_hi.le.Egamma_lo .or. Egamma_hi.le.0) return ! ... want to evaluate powers as INfrequently as possible! power_hi = Egamma_hi**g(itail) power_lo = 0.0 if (Egamma_lo.gt.0) power_lo = Egamma_lo**g(itail) ! ... obtain y region: from (Egamma_lo/Egamma_hi)**g to 1, generate flat y. ymin = power_lo/power_hi y = ymin + grnd()*(1.-ymin) x = y**(1./g(itail)) Egamma = x*Egamma_hi ! The value of our generating function at the selected point if (Egamma.gt.0) val_reciprocal = Egamma**(1.-g(itail)) * > (power_hi-power_lo) / g(itail) ! Event weight = probability radiation falls inside requested range, ! in full BASICRAD prescription weight = c(itail)/g(itail) * (power_hi-power_lo) if(itail.eq.4)itail=0 return end !------------------------------------------------------------------------- real*8 function gamma(x) implicit none integer i, n, s real*8 x, y ! Compute gamma function for xin ! The series computes gamma(1+y), and must be fed y between 0 and 1 ! This can be accomplished using the relation ! gamma(y+n+1) = (y+n)*gamma(y+n) = ... = (y+n)*...*(y+1)*gamma(y+1) gamma = 1.0 n = nint((x-1)-0.5) y = x-1 - n if (n.ne.0) then s = sign(1,n) do i = s, n, s gamma = gamma*(y+1+i)**s enddo endif gamma = gamma * (1. - 0.5748646*y + 0.9512363*y**2 - > 0.6998588*y**3 + 0.4245549*y**4 - 0.1010678*y**5) return end !-------------------------------------------------------------------------- subroutine generate_rad(main,vertex,orig,success) implicit none include 'simulate.inc' include 'radc.inc' integer i, peaked_basis_flag real*8 x, rad_weight real*8 peaked_rad_weight, basicrad_val_reciprocal real*8 basicrad_weight, extrad_phi real*8 max_delta_Trec, ebeam_min, ebeam_max logical force2 !force tail #2 to radiate. logical evsuccess,success record /event_main/ main record /event/ vertex, orig real*8 grnd !---Generate Radiation----------------------------------------------------- ! This routine radiates the incoming event ! Here are comments from radc_init that describe rad_flag and extrad_flag: ! ! First, about those (mysterious) 2 main 'radiative option' flags ... ! ! The significance of RAD_FLAG: ! RAD_FLAG = 0 .. use best available formulas, generate in ! .. (ntail,Egamma) basis ! = 1 .. use BASICRAD only, generate in (ntail,Egamma) ! .. basis ! = 2 .. use BASICRAD only, generate in (Egamma(1,2,3)) ! .. basis but prevent overlap of tails (bogus, note) ! = 3 .. use BASICRAD only, generate in (Egamma(1,2,3)) ! .. allowing radiation in all 3 directions ! .. simultaneously ! The (ntail,Egamma) basis can be called the PEAKED basis since it allows ! only 3 photon directions. PEAKED_BASIS_FLAG is set to zero below when ! the peaked basis is being used, in this way we can conveniently tell ! the BASICRAD routine to use the full Egamma formula to generate the gamma ! energy whenever it's called. ! ! (See N. Makins' thesis, section 4.5.6, for more help on this point) ! ! The significance of EXTRAD_FLAG: ! EXTRAD_FLAG = 1 .. use only BASIC external radiation formulas ! .. (phi = 1) ! = 2 .. use BASIC ext rad formulas x phi ! = 3 .. use Friedrich approximation the way we always ! .. have ! = 0 .. use DEFAULTS: 3 for RAD_FLAG = 0, 1 otherwise; note ! that the defaults mimic the hardwired 'settings' ! in SIMULATE, which doesnt read EXTRAD_FLAG ! but determines what to do based on RAD_FLAG !--------------------------------------------------------------------------- ! Initialize. ! ... ntup.radphot=sum energy of all radiated photons. ! ... ntup.radarm=ntail, or 12 if tail 2 is forced. lim1=0 lim2=0 lim3=0 phot1=0 phot2=0 phot3=0 if (debug(2)) write(6,*)'gen_rad: entering...' success = .false. rad_weight = 1 do i = 1,3 Egamma_used(i) = 0.0 enddo ntup.radphot=0. ntup.radarm=0. force2=.false. ! Which tail has to be radiated? Set ntail=0 if all tails allowed to radiate. peaked_basis_flag = 1 if (rad_flag.le.1) then peaked_basis_flag = 0 x = grnd() if (x.ge.frac(1)+frac(2)) then ntail = 3 else if (x.ge.frac(1)) then ntail = 2 else ntail = 1 endif else if (rad_flag.eq.2) then ntail = int(grnd()*3.) + 1 if (ntail.eq.4) ntail=3 else if (rad_flag.eq.3) then ntail = 0 else stop 'Idiot! rad_flag is set stupidly' endif if (debug(5)) write(6,*) 'vertexradc',vertex.e.E,edge.e.E.max ! If the scattered electron energy is too high to be dectected, only events ! where electron radiates will be accepted. Force ntail=2, and apply extra ! weight to event equal to frac(2) (the normal probability that ntail=2) if (vertex.e.E .gt. edge.e.E.max) then force2=.true. ntail=2 endif ! Need maximum change in recoil momentum for some of the later limits (e,e'p). max_delta_Trec = max((vertex.Trec - VERTEXedge.Trec.min), > (VERTEXedge.Trec.max - vertex.Trec) ) ! RADIATE TAIL #1: the incident electron tail if (doing_tail(1) .and. (ntail.eq.0.or.ntail.eq.1)) then if (debug(5)) write(6,*) 'gen_rad: at 1' ! ... CASE 1: A(e,e'p) for A>2 - the effect of radiation on this tail alters ! ... the VERTEX (actual) value of Em, so we want to make sure we're inside ! ... any source spectral function cuts (like the Fermi surface). ! ... Set Egamma_max to avoid getting below Em.min, and if Em is above Em.max, ! ... set minimum to get below Em.max. Radiating a photon decreases Em ! ... (relative to pre-radiated value), but also shifts Pm, and thus Trec. ! ... Allow Trec shift as well as Em. Note that there are NO limits based on ! ... spectrometer Em limits (cuts) because measured Em used Ebeam before rad. ! ... If only the incoming electron arm radiates, then we can compare Em ! ... to edge.Em.min (the SPECTROMETER Em limit) as well. if (doing_heavy) then if (debug(5)) write(6,*)'gen_rad: at 2a' if (debug(5)) write(6,*)'vertex.Em=',vertex.Em if (debug(5)) write(6,*)'VERTEXedge.Em.max=',VERTEXedge.Em.max if (debug(5)) write(6,*)'max_delta_Trec=',max_delta_Trec Egamma_min(1)=vertex.Em - VERTEXedge.Em.max - max_delta_Trec Egamma_max(1)=vertex.Em - VERTEXedge.Em.min + max_delta_Trec if (ntail.ne.0) Egamma_max(1)=min(Egamma_max(1), > vertex.Em - edge.Em.min + max_delta_Trec) ! ... CASE 2: H(e,e'p) - Require that Ebeam after radiation can generate ! ... electron inside of electron arm acceptance (Ein_max for Ee'_max, etc...). ! ... Also require MEASURED Em < edge.Em.max (Em_meas=egamma1+egamma2+egamma3) ! ... JRA: Note that we get ebeam max/min from E=mE'/(m+E'*(1-cos(theta)), ! ... assuming that E'_max(min) gives E_max(min). However, for very large ! ... angles, or large E'_max, you can get a negative value for E_max. When ! ... this happens, just open up ebeam_max. In the long run, we need to get ! ... a true E_max, not just E for E'_max. else if (doing_hyd_elast) then if (debug(5)) write(6,*)'gen_rad: at 2b' ebeam_max = Mp*edge.e.E.max / (Mp-edge.e.E.max*(1.-vertex.ue.z)) if (ebeam_max.lt.0) ebeam_max=1.d10 ebeam_min = Mp*edge.e.E.min / (Mp-edge.e.E.min*(1.-vertex.ue.z)) Egamma_min(1) = vertex.Ein - ebeam_max Egamma_max(1) = vertex.Ein - ebeam_min Egamma_max(1) = min(Egamma_max(1),edge.Em.max) ! ... CASE 3: D(e,e'p) - limit radiation so that E_beam > E' at the vertex. else if (doing_deuterium) then Egamma_max(1) = min(Egamma1_max, gen.sumEgen.max - vertex.e.E) if (ntail.ne.0) Egamma_min(1) = gen.sumEgen.min - vertex.e.E ! ... CASE 4: Pion production. Too complicated to calculate min/max ebeam ! ... from scattered particles. Use sumEgen-vertex.e.E as limit for remaining ! ... generated energy available. else if (doing_pion .or. doing_kaon .or. doing_rho .or. 1 doing_semi) then if (debug(5)) write(6,*)'gen_rad: at 2d' Egamma_min(1) = 0. Egamma_max(1) = gen.sumEgen.max - vertex.e.E CDJG if (ntail.ne.0) Egamma_min(1) = gen.sumEgen.min - vertex.e.E endif ! ... Constrain Egamma 'SIMC found an Egamma(1) limit above your hardwired value' endif ! ... radiate if (debug(5)) write(6,*)'gen_rad: init Egamma_min,max =', > Egamma_min(1),Egamma_max(1) call basicrad(1*peaked_basis_flag,Egamma_min(1),Egamma_max(1), > Egamma_used(1),basicrad_weight,basicrad_val_reciprocal) if (basicrad_weight.le.0) then if (debug(4)) write(6,*)'gen_rad: failed at 4' return endif if (debug(5)) write(6,*)'Egamma_used',Egamma_used(1) ! ... adjust kinematics and call complete_ev. If complete_ev fails, return. vertex.Ein = vertex.Ein - Egamma_used(1) if (debug(5)) write(6,*) 'calling complete_ev from radc' call complete_ev(main,vertex,evsuccess) if (.not.evsuccess) then !remove previous radiation, try again. !! if (ntried.le.5000) write(6,'(1x,''COMPLETE_EV1 failed in GENERATE_RAD at event'',i7,''!'')') nevent return endif if (debug(4)) write(6,*)'gen_rad: at 5' rad_weight = rad_weight * basicrad_weight if (rad_flag.le.1) then rad_weight = peaked_rad_weight(vertex,Egamma_used(1),Egamma_min(1), > Egamma_max(1),basicrad_val_reciprocal,basicrad_weight) else rad_weight=rad_weight*extrad_phi(1,vertex.Ein,vertex.e.E,Egamma_used(1)) endif endif !end of tail 1 (incoming electron) ! The VERTEX kinematics are now set, so check that we are inside the ! VERTEXedges (SF limits). if (doing_heavy) then if (debug(5)) write(6,*)'gen_rad: Em, min, max =',vertex.Em, > VERTEXedge.Em.min,VERTEXedge.Em.max if (debug(5)) write(6,*)'gen_rad: Pm, min, max =',vertex.Pm, > VERTEXedge.Pm.min,VERTEXedge.Pm.max if (vertex.Em .lt. VERTEXedge.Em.min .or. > vertex.Em .gt. VERTEXedge.Em.max .or. > vertex.Pm .lt. VERTEXedge.Pm.min .or. > vertex.Pm .gt. VERTEXedge.Pm.max) then if (debug(4)) write(6,*)'gen_rad: failed at 6' return endif endif ! RADIATE TAIL #2: the scattered electron tail if (doing_tail(2) .and. (ntail.eq.0.or.ntail.eq.2)) then if (debug(5)) write(6,*) 'we are at 2' ! ... Find min/max Egammas that will make in into electron acceptance Egamma_min(2) = vertex.e.E - edge.e.E.max Egamma_max(2) = vertex.e.E - edge.e.E.min ! ... Measured Em must be within acceptance after radiation. Gives upper ! ... limit on sum of photon energies, and lower limit if hadron arm does ! ... not radiate. Em_MEAS = Em_VERT+(Egamma1+Egamma2+Egamma3)+/-delta_trec. if (doing_eep) then Egamma_max(2) = min(Egamma_max(2), > (edge.Em.max - vertex.Em) - Egamma_used(1) + max_delta_Trec ) if (ntail.ne.0 .or. .not.rad_proton_this_ev) > Egamma_min(2) = max(Egamma_min(2), > (edge.Em.min - vertex.Em) - Egamma_used(1) - max_delta_Trec ) endif Egamma_max(2) = min(Egamma_max(2),Egamma_tot_max-Egamma_used(1)) Egamma_min(2) = Egamma_min(2) - dE_edge_test Egamma_max(2) = Egamma_max(2) + dE_edge_test ! ... override Egamma_max with limit from dbase (if entered) lim2=egamma_max(2) if (hardwired_rad) then Egamma_max(2) = Egamma_gen_max if (lim2.gt.Egamma_gen_max) write(6,*) > 'SIMC found an Egamma(2) limit above your hardwired value' endif ! ... radiate call basicrad(2*peaked_basis_flag,Egamma_min(2),Egamma_max(2), > Egamma_used(2),basicrad_weight,basicrad_val_reciprocal) if (basicrad_weight.le.0) then ! write(6,*)'gen_rad: Maybe you need to change Em limits!' ! write(6,*)'min/max(2)=',egamma_min(2),egamma_max(2) if (debug(4)) write(6,*)'gen_rad: failed at 7' return endif if(debug(5))write(6,*)'gen_rad: peaked_basis_flag=',peaked_basis_flag if(debug(5))write(6,*)'gen_rad: Egamma min,max,used=', > Egamma_min(2),Egamma_max(2),Egamma_used(2) if(debug(5))write(6,*)'gen_rad: basicrad_weight=',basicrad_weight rad_weight = rad_weight * basicrad_weight if (rad_flag.le.1) then rad_weight = peaked_rad_weight(vertex,Egamma_used(2),Egamma_min(2), > Egamma_max(2),basicrad_val_reciprocal,basicrad_weight) else rad_weight=rad_weight*extrad_phi(2,vertex.Ein,vertex.e.E,Egamma_used(2)) endif endif !end of tail 2 (outgoing electron) ! RADIATE TAIL #3: the scattered proton tail if (rad_proton_this_ev .and. (ntail.eq.0.or.ntail.eq.3)) then if (debug(5)) write(6,*) 'we are at 3' ! ... Find min/max Egammas that will make in into hadron acceptance Egamma_min(3) = vertex.p.E - edge.p.E.max Egamma_max(3) = vertex.p.E - edge.p.E.min ! ... Measured Em must be within acceptance after radiation. Gives upper ! ... and lower limit on sum of photon energies (and thus on hadron radiation, ! ... since other arms are done). ! ... Em_MEAS = Em_VERT+(Egamma1+Egamma2+Egamma3)+/-delta_trec. if (doing_eep) then Egamma_max(3) = min(Egamma_max(3), (edge.Em.max - vertex.Em) - > Egamma_used(1) - Egamma_used(2) + max_delta_Trec) Egamma_min(3) = max(Egamma_min(3), (edge.Em.min - vertex.Em) - > Egamma_used(1) - Egamma_used(2) - max_delta_Trec ) endif Egamma_max(3) = min(Egamma_max(3), > Egamma_tot_max-Egamma_used(1)-Egamma_used(2) ) Egamma_min(3) = Egamma_min(3) - dE_edge_test Egamma_max(3) = Egamma_max(3) + dE_edge_test ! ... override Egamma_max with limit from dbase (if entered) lim3=egamma_max(3) if (hardwired_rad) then Egamma_max(3) = Egamma_gen_max if (lim3.gt.Egamma_gen_max) write(6,*) > 'SIMC found an Egamma(3 ) limit above your hardwired value' endif ! ... radiate call basicrad(3*peaked_basis_flag,Egamma_min(3),Egamma_max(3), > Egamma_used(3),basicrad_weight,basicrad_val_reciprocal) if (basicrad_weight.le.0) then if (debug(4)) write(6,*)'gen_rad: failed at 8' return endif rad_weight = rad_weight * basicrad_weight if (rad_flag.le.1) then rad_weight = peaked_rad_weight(vertex,Egamma_used(3),Egamma_min(3), > Egamma_max(3),basicrad_val_reciprocal,basicrad_weight) else rad_weight=rad_weight*extrad_phi(3,vertex.Ein,vertex.e.E,Egamma_used(3)) endif endif !end of tail 3 (outgoing hadron) ! Now obtain description of ORIG (orig = vertex + radiation) event. ! ... Note that no kinematics have been changed yet, EXCEPT vertex.Ein, which ! ... has been reduced by Egamma_used(1). orig.Ein is the initial (unradiated) ! ... beam energy, so add back Egamma_used(1). Reduce orig.e.E and orig.p.E ! ... here, and they will be the values passed to the monte carlo (for mult. ! ... scattering and monte carlo model). ! ... Remember that we have taken Ein(GEN)-Ebeam_vertex_ave into ! ... account by shifting edge.Em) do i = 1, neventfields orig.all(i) = vertex.all(i) enddo if (debug(4)) write(6,*)'gen_rad: at 10: orig.p.E =',orig.p.E ! ... remove radiation from incoming electron. orig.Ein = vertex.Ein + Egamma_used(1) ! ... adjust the e'/hadron momenta (if they radiated significantly) orig.e.E = vertex.e.E - Egamma_used(2) orig.e.P = orig.e.E orig.e.delta = (orig.e.P-spec.e.P)/spec.e.P*100. orig.p.E = vertex.p.E - Egamma_used(3) if(orig.p.E.le.Mh) then if (debug(4)) write(6,*)'gen_rad: failed at 11' return endif orig.p.P = sqrt(orig.p.E**2-Mh2) orig.p.delta = (orig.p.P-spec.p.P)/spec.p.P*100. ! ... Record some information for other routines. phot1=Egamma_used(1) phot2=Egamma_used(2) phot3=Egamma_used(3) ntup.radphot = phot1 + phot2 + phot3 ntup.radarm=ntail if (force2) ntup.radarm=12 ! Complete determination of the event weight due to radiation. main.gen_weight = main.gen_weight*rad_weight/hardcorfac if (force2) main.gen_weight=main.gen_weight*frac(2) if (debug(5)) write(6,*) 'mgw',main.gen_weight,rad_weight,hardcorfac success = .true. if (debug(2)) write(6,*)'gen_rad: exiting...' return end !-------------------------------------------------------------------------- real*8 function peaked_rad_weight(vertex,Egamma, > emin,emax,basicrad_val_reciprocal,basicrad_weight) implicit none include 'radc.inc' include 'structures.inc' real*8 Egamma, basicrad_val_reciprocal, basicrad_weight real*8 t1, t2, r, phi_ext, ein, eout, emin, emax real*8 brem, bremos, extrad_phi, gamma real*8 dhard, dsoft_intmin, dsoft_intmax real*8 dsoft_ext1, dsoft_ext2, dsoft_extmin real*8 dsoft_int_primemin, dsoft_int_primemax real*8 dsoft_ext1_prime, dsoft_ext2_prime real*8 dsoft_ext_prime, dsoft_extmax, eul real*8 dsoft_int,dsoft_ext,dsoft_int_prime record /event/ vertex real*8 zero parameter (zero=0.0d0) !double precision zero for subroutine calls basicrad_val_reciprocal=basicrad_val_reciprocal+0. !avoid unused variable error ! Compute a more precise value for the radiative probability at the ! selected Egamma, more precise than the BASICRAD distribution used to ! generate Egamma that is. ! ! NB: This subroutine only deals with calculations done in the PEAKING ! APPROXIMATION basis - ! i.e. we only get here if RAD_FLAG = 0 or 1 ! ! ... (sort of) KLUGE -- if Egamma < res limit, these things might blow up, ! ... so just screw the correction and leave ein=vertex.Ein eout=vertex.e.E eul=0.577215665 c if (Egamma.lt.Egamma_res_limit) then c peaked_rad_weight = basicrad_weight c return c endif ! ... External phi_ext = 1.0 if (extrad_flag.le.2) then phi_ext = extrad_phi(0,ein,eout,Egamma) if (rad_flag.eq.1) then peaked_rad_weight = basicrad_weight * phi_ext return endif if(emin.gt.0)then dsoft_extmin = log(g_ext/c_ext(0)/emin**g_ext) endif dsoft_extmax = log(g_ext/c_ext(0)/emax**g_ext) dsoft_ext_prime = -g_ext/Egamma else t1 = bt(1)/etatzai t2 = bt(2)/etatzai call extrad_friedrich(ein,Egamma,t1,dsoft_ext1,dsoft_ext1_prime) call extrad_friedrich(eout,Egamma,t2,dsoft_ext2,dsoft_ext2_prime) dsoft_ext = dsoft_ext1 + dsoft_ext2 dsoft_ext_prime = dsoft_ext1_prime + dsoft_ext2_prime endif ! ... Internal ! ........ use full BREM calculation of deltas if (rad_flag.eq.0) then if (.not.use_offshell_rad) then if(emin.gt.0)then r = brem(ein,eout,emin,rad_proton_this_ev, > dsoft_intmin,dhard,dsoft_int_primemin) else dsoft_intmin=1.0 endif r = brem(ein,eout,emax,rad_proton_this_ev, > dsoft_intmax,dhard,dsoft_int_primemax) else if(emin.gt.0)then r = bremos(emin, zero, zero, ein, > vertex.e.E*vertex.ue.x, vertex.e.E*vertex.ue.y, > vertex.e.E*vertex.ue.z, c > vertex.Pmx, vertex.Pmy, vertex.Pmz, > zero,zero,zero, > vertex.p.P*vertex.up.x, vertex.p.P*vertex.up.y, > vertex.p.P*vertex.up.z, vertex.p.E, > rad_proton_this_ev, > dsoft_intmin, dhard, dsoft_int_primemin) else dsoft_intmin=1.0 endif r = bremos(emax, zero, zero, ein, > vertex.e.E*vertex.ue.x, vertex.e.E*vertex.ue.y, > vertex.e.E*vertex.ue.z, c > vertex.Pmx, vertex.Pmy, vertex.Pmz > zero,zero,zero, > vertex.p.P*vertex.up.x, vertex.p.P*vertex.up.y, > vertex.p.P*vertex.up.z, vertex.p.E, > rad_proton_this_ev, > dsoft_intmax, dhard, dsoft_int_primemax) endif ! ........ use basic calculation of internal deltas else dsoft_int = log(g_int/c_int(0)/Egamma**g_int) dsoft_int_prime = -g_int/Egamma endif ! ... All together now if(emin.gt.0)then peaked_rad_weight = c_ext(0)/g_ext*(exp(-dsoft_intmax)* > emax**g_ext-exp(-dsoft_intmin)*emin**g_ext) else peaked_rad_weight = c_ext(0)/g_ext*(exp(-dsoft_intmax)*emax**g_ext) endif peaked_rad_weight = peaked_rad_weight * exp(-eul*g(4))/gamma(1.+g(4)) > * gamma(1.+g(4)-bt(1)-bt(2))*gamma(1.+bt(1)) > * gamma(1.+bt(2))/gamma(1.+g(4)) if (peaked_rad_weight.lt.0) peaked_rad_weight = 0 return end !--------------------------------------------------------------------------- subroutine extrad_friedrich(Ei,Ecutoff,trad,dbrem,dbrem_prime) implicit none include 'simulate.inc' include 'radc.inc' real*8 Ei, Ecutoff, trad, x, dbrem, dbrem_prime x = Ecutoff/Ei dbrem = trad*(-(etatzai-0.5) - etatzai*log(x) + etatzai*x - 0.5*x**2) dbrem_prime = -trad/Ei * (etatzai/x - etatzai + x) return end !-------------------------------------------------------------------------- real*8 function extrad_phi(itail,E1,E2,Egamma) implicit none include 'radc.inc' integer itail real*8 E1, E2, Egamma, E(2), x, t, gamma E(1) = E1 E(2) = E2 ! Compute the multiplicative external correction function PHI ! ... CASE 1: phi = 1 extrad_phi = 1.0 ! ... CASE 2: phi correctly computed according to Dave's formulas if (extrad_flag.eq.2) then if (itail.eq.0) then extrad_phi = 1. - (bt(1)/E(1)+bt(2)/E(2)) / (g(1)+g(2)) * Egamma else if (itail.eq.1.or.itail.eq.2) then extrad_phi = 1. - bt(itail)/E(itail) / g(itail) * Egamma endif ! ... CASE 3: phi computed in Friedrich prescription else if (extrad_flag.eq.3) then if (itail.eq.0) stop 'Idiot! a multiplicative factor EXTRAD_PHI is not defined for peaking approx and EXTRAD_FLAG>2!' if (itail.eq.1.or.itail.eq.2) then x = Egamma/E(itail) t = bt(itail)/etatzai extrad_phi = extrad_phi * (1. - x + x**2/etatzai) * > exp(t*((etatzai-0.5)-etatzai*x+x**2/2.)) * gamma(1.+bt(itail)) endif endif return end !------------------------------------------------------------------------ real*8 function schwinger(Ecutoff,vertex,include_hard,dsoft,dhard) implicit none include 'simulate.inc' include 'radc.inc' real*8 Ecutoff,lq,s2,b,spence,dsoft,dhard,spen logical include_hard record /event/ vertex lq = log(vertex.Q2/Me**2)-1.0 s2 = sin(vertex.e.theta/2.)**2 b = 1.+2.*vertex.nu*s2/(targ.A*amu) spence = spen(1.-s2)-2.5893784 dsoft = alpi*lq*log( vertex.Ein/(etta**2)*vertex.e.E*b/(Ecutoff**2) ) dhard = -alpi*(2.166666*lq + spence - log(vertex.Ein/vertex.e.E)**2/2.0) if (use_expon.eq.0) then schwinger = exp(dsoft) else schwinger = 1.+dsoft endif if (include_hard) schwinger = schwinger/(1.-dhard) return end !-------------------------------------------------------------------------- real*8 function spen(x) implicit none integer i real*8 x,s,y ! Approximation formula for the Spence-function or -integral according to ! abramowitz and stegun : formula 27.7.2 y = 1.0 s = 0.0 i = 0 do while (i.le.100 .and. abs(y).gt.abs(s)*1.d-4) i = i+1 y = x*y s = s+y/float(i**2) enddo spen = s return end !-------------------------------------------------------------------------- real*8 function lambda_dave(itail,plus_flag,doing_proton,e1,e2,e3,p3,th) implicit none include 'constants.inc' integer itail,plus_flag real*8 e1,e2,e3,p3,th real*8 plus_term logical doing_proton, warned/.false./ ! The extended peaking approximation plus_term = 0.0 if (plus_flag.eq.1 .and. itail.lt.3) then plus_term = log((1.-cos(th))/2.) ! ... only add in term due to ep interference if we're using proton ! ... radiation if (doing_proton) plus_term = plus_term + 2.*log(e1/e2) endif ! Compute lambdas if (itail.eq.1) then lambda_dave = alpi*(2.*log(2.*e1/Me) -1. + plus_term) else if (itail.eq.2) then lambda_dave = alpi*(2.*log(2.*e2/Me) -1. + plus_term) else if (itail.eq.3) then if (doing_proton) then lambda_dave = alpi*(log((e3+p3)/(e3-p3)) - 2.) ! ... at low energies, this may fritz, prevent human from causing itself ! ... too much damage if (lambda_dave.lt.0) then lambda_dave = 0.0 if (.not.warned) then write(6,'(1x,a)') '+-----------------------------------+' write(6,'(1x,a)') '| WARNING: lambda_dave(3)<0, think |' write(6,'(1x,a)') '| about running with proton rad OFF |' write(6,'(1x,a)') '+-----------------------------------+' warned = .true. endif endif else lambda_dave = 0.0 endif endif return end