subroutine target_init(using_Eloss,using_Coulomb,the_cent,thp_cent, > Pp_cent,Mh2,ebeam,Pe_cent) implicit none include 'constants.inc' include 'target.inc' integer i real*8 the_cent, thp_cent, Pp_cent, Mh2, ebeam, energy real*8 Pe_cent real*8 za2, fc logical using_Eloss, using_Coulomb real*8 zero parameter (zero=0.0d0) !double precision zero for subroutines calls. ! The radiation length of the target targ.L1 = log(184.15) - log(targ.Z)/3.0 targ.L2 = log(1194.) - 2.*log(targ.Z)/3.0 if(targ.Z.eq.1)then targ.L1=5.31 targ.L2=6.144 endif za2 = (targ.Z*alpha)**2 fc = za2*(1.202+za2*(-1.0369+za2*1.008/(za2+1))) ! ... The radiation length, in both g/cm2 and cm; For H,D,4He, using ! ... PDB values (1999), other calculated directly from Tsai formula ! ... (For 4He, the Tsai formula is 10% lower). if (nint(targ.A).eq.1) then targ.X0 = 61.28 else if (nint(targ.A).eq.2) then targ.X0 = 122.4 else if (nint(targ.A).eq.4) then targ.X0 = 94.32 else if (nint(targ.A).eq.3) then write(6,*) 'Using Tsai formula for He-3 radiation length!!! (is there a better value?)' targ.X0 = 716.405*targ.A/targ.Z/(targ.Z*(targ.L1-fc)+targ.L2) else targ.X0 = 716.405*targ.A/targ.Z/(targ.Z*(targ.L1-fc)+targ.L2) endif targ.X0_cm = targ.X0/targ.rho ! ... 'Average' ionization losses (MOST PROBABLE, REALLY). Printed out by simc ! ... (so ave or m.p. are both OK), and Ebeam_vertex_ave (which should be m.p.) energy = ebeam call trip_thru_target (1,zero,energy,zero,targ.Eloss(1).ave, > targ.teff(1).ave,Me,4) energy = sqrt(Pe_cent**2+Me2) call trip_thru_target (2,zero,energy,the_cent,targ.Eloss(2).ave, > targ.teff(2).ave,Me,4) energy = sqrt(Pp_cent**2 + Mh2) call trip_thru_target (3,zero,energy,thp_cent,targ.Eloss(3).ave, > targ.teff(3).ave,sqrt(Mh2),4) if (.not.using_Eloss) then do i = 1, 3 targ.Eloss(i).ave = 0.0 enddo endif ! ... Coulomb potential energy (NB All of the following are positive) if (using_Coulomb) then targ.Coulomb.ave=6./5.*(targ.Z-1.)*alpha*hbarc/(1.18*targ.A**(1./3.)) targ.Coulomb_constant = 5./12. * targ.Coulomb.ave targ.Coulomb.min = targ.Coulomb_constant * 2.0 targ.Coulomb.max = targ.Coulomb_constant * 3.0 else targ.Coulomb.ave = 0.0 targ.Coulomb_constant = 0.0 targ.Coulomb.min = 0.0 targ.Coulomb.max = 0.0 endif return end !---------------------------------------------------------------------- subroutine limits_init(H) implicit none include 'simulate.inc' include 'radc.inc' include 'histograms.inc' integer i real*8 r real*8 Ebeam_min, Ebeam_max real*8 t1,t2 !temp. variables. real*8 slop_Coulomb, slop_Ebeam, slop_Ee, slop_Ep record /cuts/ the_phys, thp_phys, z, pp record /histograms/ H !----------------------------------------------------------------------- ! RECORDS store information relating to various quantities at several ! 'STATES' in event description: ! (of course, not all qties are _recorded_ in each of these stages) ! ! orig.* qties are TRUE qties but before ANY radiation ! vertex.* qties are those used to determine spec fn and sigcc ! weighting, so TRUE qties before tails 2,3 radiated ! recon.* qties are defined AT the spectrometers, AFTER being ! run through the single arm montecarlos ! ! main.SP.* qties are defined AT the spectrometers, so TRUE ! qties AFTER any modifications due to non-radiative ! interaction with the target (this includes Coulomb ! accel/decel of initial/final electron) ! main.FP.* spect. focal plane values-if using spect. model ! main.RECON.* qties are defined AT the spectrometers, AFTER being ! run through the single arm montecarlos ! ! Note that in the absence of radiation, GEN=VERTEX ! Note that if a spectrometer is not used, RECON=SP, FP is empty. ! Note that the generated qties we have to compare with gen limits later ! are unaffected by radiation in tail1, so equal to vertex qties. ! ! ! EDGES on event qties at any stage can be derived from these cuts, by taking ! the transformations (and associated uncertainties) between each stage into ! account. The general usefulness of defining these edges is to improve ! efficiency --> provide checks at intermediate stages to allow us to abort ! an event early, and provide some foreknowledge in generation routines of ! what will make it to the end. We must be certain that these derived EDGES ! are wide enough that the events making it through the cuts don't come ! uncomfortably close to them. ! ! SPedge.* Limits at spectrometer = Acceptance + slop ! edge.* Limits after interaction = SPedge + musc or eloss ! = VERTEXedge + coulomb + rad. ! VERTEXedge.* Limits at vertex (from energy conservation). ! gen.* Generation limits. ! ! SLOPS ! slop.MC.* Due to spectrometer resolution (spectrometer MC). ! slop_* Slop in calculated physics quantities (from spect. ! resolution + range of eloss/coulomb/... ! !---------------------------------------------------------------------- ! ... get slop values for the difference between orig and recon ! ... SPECTROMETER quantities (recon inaccuracies due to the montecarlos) if (using_E_arm_montecarlo) then if (electron_arm.eq.1) then slop.MC.e.delta.used = slop_param_d_HMS slop.MC.e.yptar.used = slop_param_t_HMS slop.MC.e.xptar.used = slop_param_p_HMS else if (electron_arm.eq.2) then slop.MC.e.delta.used = slop_param_d_SOS slop.MC.e.yptar.used = slop_param_t_SOS slop.MC.e.xptar.used = slop_param_p_SOS else if (electron_arm.eq.3) then slop.MC.e.delta.used = slop_param_d_HRSR slop.MC.e.yptar.used = slop_param_t_HRSR slop.MC.e.xptar.used = slop_param_p_HRSR else if (electron_arm.eq.4) then slop.MC.e.delta.used = slop_param_d_HRSL slop.MC.e.yptar.used = slop_param_t_HRSL slop.MC.e.xptar.used = slop_param_p_HRSL endif endif if (using_P_arm_montecarlo) then if (hadron_arm.eq.1) then slop.MC.p.delta.used = slop_param_d_HMS slop.MC.p.yptar.used = slop_param_t_HMS slop.MC.p.xptar.used = slop_param_p_HMS else if (hadron_arm.eq.2) then slop.MC.p.delta.used = slop_param_d_SOS slop.MC.p.yptar.used = slop_param_t_SOS slop.MC.p.xptar.used = slop_param_p_SOS else if (hadron_arm.eq.3) then slop.MC.p.delta.used = slop_param_d_HRSR slop.MC.p.yptar.used = slop_param_t_HRSR slop.MC.p.xptar.used = slop_param_p_HRSR else if (hadron_arm.eq.4) then slop.MC.p.delta.used = slop_param_d_HRSL slop.MC.p.yptar.used = slop_param_t_HRSL slop.MC.p.xptar.used = slop_param_p_HRSL endif endif ! ... Add slop to SPedges. Used as input to final edge.*.*.* and to ! ... calculate minimum/maximum theta_phys for extreme_trip_thru_target. SPedge.e.delta.min = SPedge.e.delta.min - slop.MC.e.delta.used SPedge.e.delta.max = SPedge.e.delta.max + slop.MC.e.delta.used SPedge.e.yptar.min = SPedge.e.yptar.min - slop.MC.e.yptar.used SPedge.e.yptar.max = SPedge.e.yptar.max + slop.MC.e.yptar.used SPedge.e.xptar.min = SPedge.e.xptar.min - slop.MC.e.xptar.used SPedge.e.xptar.max = SPedge.e.xptar.max + slop.MC.e.xptar.used SPedge.p.delta.min = SPedge.p.delta.min - slop.MC.p.delta.used SPedge.p.delta.max = SPedge.p.delta.max + slop.MC.p.delta.used SPedge.p.yptar.min = SPedge.p.yptar.min - slop.MC.p.yptar.used SPedge.p.yptar.max = SPedge.p.yptar.max + slop.MC.p.yptar.used SPedge.p.xptar.min = SPedge.p.xptar.min - slop.MC.p.xptar.used SPedge.p.xptar.max = SPedge.p.xptar.max + slop.MC.p.xptar.used ! Compute TRUE edges -- distortions in the target come into play edge.e.E.min = (1.+SPedge.e.delta.min/100.)*spec.e.P + > targ.Coulomb.min - dE_edge_test edge.e.E.max = (1.+SPedge.e.delta.max/100.)*spec.e.P + > targ.Coulomb.max + dE_edge_test pp.min = (1.+SPedge.p.delta.min/100.)*spec.p.P - dE_edge_test pp.max = (1.+SPedge.p.delta.max/100.)*spec.p.P + dE_edge_test pp.min = max(0.001d0,pp.min) !avoid p=0 (which can lead to div by zero)( edge.p.E.min = sqrt(pp.min**2 + Mh2) edge.p.E.max = sqrt(pp.max**2 + Mh2) ! ... extreme theta_e,theta_p,z values for extreme_trip_thru_target. the_phys.max = acos( (spec.e.cos_th-spec.e.sin_th*SPedge.e.yptar.max)/ > sqrt(1.+SPedge.e.yptar.max**2+SPedge.e.xptar.max**2) ) the_phys.min = acos( (spec.e.cos_th-spec.e.sin_th*SPedge.e.yptar.min)/ > sqrt(1.+SPedge.e.yptar.min**2) ) thp_phys.max = acos( (spec.p.cos_th-spec.p.sin_th*SPedge.p.yptar.max)/ > sqrt(1.+SPedge.p.yptar.max**2+SPedge.p.xptar.max**2) ) thp_phys.min = acos( (spec.p.cos_th-spec.p.sin_th*SPedge.p.yptar.min)/ > sqrt(1.+SPedge.p.yptar.min**2) ) z.min = -0.5*targ.length z.max = 0.5*targ.length call extreme_trip_thru_target(Ebeam, the_phys, thp_phys, edge.e.E, > pp, z, Mh) if (.not.using_Eloss) then do i = 1, 3 targ.Eloss(i).min = 0.0 targ.Eloss(i).max = 0.0 enddo endif if (.not.mc_smear) then targ.musc_max(1)=0. targ.musc_max(2)=0. targ.musc_max(3)=0. endif edge.e.E.min = edge.e.E.min + targ.Eloss(2).min edge.e.E.max = edge.e.E.max + targ.Eloss(2).max edge.p.E.min = edge.p.E.min + targ.Eloss(3).min edge.p.E.max = edge.p.E.max + targ.Eloss(3).max edge.e.yptar.min = SPedge.e.yptar.min - targ.musc_max(2) edge.e.yptar.max = SPedge.e.yptar.max + targ.musc_max(2) edge.e.xptar.min = SPedge.e.xptar.min - targ.musc_max(2) edge.e.xptar.max = SPedge.e.xptar.max + targ.musc_max(2) edge.p.yptar.min = SPedge.p.yptar.min - targ.musc_max(3) edge.p.yptar.max = SPedge.p.yptar.max + targ.musc_max(3) edge.p.xptar.min = SPedge.p.xptar.min - targ.musc_max(3) edge.p.xptar.max = SPedge.p.xptar.max + targ.musc_max(3) ! Edges on values of Em and Pm BEFORE reconstruction. Need to apply slop to ! take into account all transformations from ORIGINAL TRUE values to ! RECONSTRUCTED TRUE values --> that includes: (a) reconstruction slop on ! spectrometer values (b) variation in the beam energy (c) difference between ! actual losses and the corrections made for them ! ... Save the 'measured' beam energy, which will be used in recon. Ebeam_vertex_ave = Ebeam + targ.Coulomb.ave - targ.Eloss(1).ave ! ... Calculate slop in energies and momenta. Ebeam_max = Ebeam + dEbeam/2. - targ.Eloss(1).min + targ.Coulomb.max Ebeam_min = Ebeam - dEbeam/2. - targ.Eloss(1).max + targ.Coulomb.min slop_Coulomb = targ.Coulomb.max - targ.Coulomb.ave slop_Ebeam = dEbeam/2. + slop_Coulomb slop_Ee = slop.MC.e.delta.used/100.*spec.e.P + slop_Coulomb r = sqrt(edge.p.E.max**2 - Mh2) slop_Ep = sqrt( (r + slop.MC.p.delta.used/100.*spec.p.P)**2 + Mh2 ) - > edge.p.E.max slop_Ebeam = slop_Ebeam + (targ.Eloss(1).max-targ.Eloss(1).min) slop_Ee = slop_Ee + (targ.Eloss(2).max-targ.Eloss(2).min) slop_Ep = slop_Ep + (targ.Eloss(3).max-targ.Eloss(3).min) if (doing_heavy) then ! 'reconstructed' Em cuts. slop.total.Em.used = slop_Ebeam + slop_Ee + slop_Ep + dE_edge_test edge.Em.min = cuts.Em.min - slop.total.Em.used edge.Em.max = cuts.Em.max + slop.total.Em.used edge.Em.min = max(0.d0,edge.Em.min) endif ! Edges on Em, Pm, etc... VERTEXedge.* values are vertex limits. edge.* values ! are limits at spectrometer (after eloss and e'/hadron radiation). Remember ! that min/max values are initialized to wide open values in STRUCTURES. ! Need edge.Em to limit radiated photon energies. ! Need VERTEXedge.Em for calulating limits on radiatied photons. ! Need VERTEXedge.Pm for A(e,e'p) only (for generation limits on Ee', Ep'). ! Note that Pm_theory(*).min/max might allow for positive and negative Pm, ! or it could have positive only. We want *.Pm.min to be zero if zero is ! allowed, so we have to check for Pm_theory.min being negative. ! For (e,e'p), need Trec (for A-1 system) in order to set radiation limits. ! Calculate Trec max from limits on Pm (since P of recoiling A-1 system is Pm). ! For pions, want Trec to be T of recoiling nucleon (the struck nucleon). ! Can only get limits from general energy conservation, and can only get ! upper limit, since the lower limit is determined by the allowed radiation, ! which is not calculated yet (and needs Trec to be calculated). ! For (e,e'pi) and (e,e'K), there are two 'recoil' systems. Trec is the ! spectator (A-1) recoil system. Trec_struck is the recoiling baryon ! from the 'struck' nucleon (i.e. the neutron from pi+ production from a ! proton, and the hyperon from the kaon production). ! Get radiation limits. In all cases, total energy conservation gives ! the primary limit. The doing_heavy case has a second condition. Radiated ! photons change Em from the vertex value to the measured value. They also ! change Pm, which can modify Trec. So the max. change in Em is for a ! minimum generated Em (VERTEXedge.Em.min) that ends up as a maximum measured Em ! (edge.Em.max), with slop to take into account the modification to Trec. if (doing_hyd_elast) then VERTEXedge.Em.min = 0.0 VERTEXedge.Em.max = 0.0 VERTEXedge.Pm.min = 0.0 VERTEXedge.Pm.max = 0.0 VERTEXedge.Trec.min = 0.0 VERTEXedge.Trec.max = 0.0 Egamma_tot_max = Ebeam_max + targ.M - edge.e.E.min - edge.p.E.min else if (doing_deuterium) then VERTEXedge.Em.min = Mp + Mn - targ.M !2.2249 MeV, I hope. VERTEXedge.Em.max = Mp + Mn - targ.M if (Pm_theory(1).min .le. 0.0 .and. Pm_theory(1).max .ge. 0.0) then VERTEXedge.Pm.min = 0.0 else VERTEXedge.Pm.min = Pm_theory(1).min endif VERTEXedge.Pm.max = max(abs(Pm_theory(1).min),abs(Pm_theory(1).max)) VERTEXedge.Trec.min=sqrt(targ.Mrec**2+VERTEXedge.Pm.min**2)-targ.Mrec VERTEXedge.Trec.max=sqrt(targ.Mrec**2+VERTEXedge.Pm.max**2)-targ.Mrec Egamma_tot_max = Ebeam_max + targ.Mtar_struck - VERTEXedge.Em.min - > edge.e.E.min - edge.p.E.min - VERTEXedge.Trec.min else if (doing_heavy) then VERTEXedge.Pm.min=1.d10 VERTEXedge.Pm.max=-1.d10 do i = 1, nrhoPm if (Pm_theory(i).min .le. 0.0 .and. Pm_theory(i).max .ge. 0.0) then t1 = 0.0 else t1 = Pm_theory(1).min endif t2=max(abs(Pm_theory(i).min),abs(Pm_theory(i).max)) VERTEXedge.Pm.min = min(VERTEXedge.Pm.min,t1) VERTEXedge.Pm.max = max(VERTEXedge.Pm.max,t2) enddo VERTEXedge.Em.min = E_Fermi VERTEXedge.Em.max = 1.d10 VERTEXedge.Trec.min=sqrt(targ.Mrec**2+VERTEXedge.Pm.min**2)-targ.Mrec VERTEXedge.Trec.max=sqrt(targ.Mrec**2+VERTEXedge.Pm.max**2)-targ.Mrec ! Radiation gives reconstructed Em larger than the true Em, so apply ! limits on measured Em to vertex Em + Egamma. t1 = Ebeam_max + targ.Mtar_struck - VERTEXedge.Em.min - > edge.e.E.min - edge.p.E.min - VERTEXedge.Trec.min t2 = (edge.Em.max - VERTEXedge.Em.min) + > (VERTEXedge.Trec.max - VERTEXedge.Trec.min) Egamma_tot_max = min(t1,t2) ! ... Note that we check the hardwired_rad and using_rad flags here, instead of ! ... waiting until the end as in the other cases. This is because doing_heavy ! ... uses Egamma_tot_max in calculation of VERTEXedge.Em.min/max limits. if (hardwired_rad) Egamma_tot_max = Egamma_gen_max if (.not.using_rad) Egamma_tot_max = 0.0 VERTEXedge.Em.min = max(VERTEXedge.Em.min,edge.Em.min-Egamma_tot_max) VERTEXedge.Em.max = min(VERTEXedge.Em.max,edge.Em.max) else if (doing_hydpi .or. doing_hydkaon .or. doing_hyddvcs) then VERTEXedge.Em.min = 0.0 VERTEXedge.Em.max = 0.0 VERTEXedge.Pm.min = 0.0 VERTEXedge.Pm.max = 0.0 VERTEXedge.Trec.min = 0. !T_recoil for A-1 system. VERTEXedge.Trec.max = 0. VERTEXedge.Trec_struck.min = 0. VERTEXedge.Trec_struck.max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck - > edge.e.E.min - edge.p.E.min Egamma_tot_max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck - > edge.e.E.min - edge.p.E.min - VERTEXedge.Trec_struck.min else if (doing_deutpi .or. doing_deutkaon .or. doing_deutdvcs) then VERTEXedge.Em.min = Mp + Mn - targ.M !2.2249 MeV, I hope. VERTEXedge.Em.max = Mp + Mn - targ.M VERTEXedge.Pm.min = 0.0 VERTEXedge.Pm.max = pval(nump) VERTEXedge.Trec.min = sqrt(targ.Mrec**2+VERTEXedge.Pm.min**2)-targ.Mrec VERTEXedge.Trec.max = sqrt(targ.Mrec**2+VERTEXedge.Pm.max**2)-targ.Mrec VERTEXedge.Trec_struck.min = 0.0 VERTEXedge.Trec_struck.max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck - > edge.e.E.min - edge.p.E.min - VERTEXedge.Em.min - > VERTEXedge.Trec.min Egamma_tot_max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck - > edge.e.E.min - edge.p.E.min - VERTEXedge.Em.min - > VERTEXedge.Trec.min - VERTEXedge.Trec_struck.min else if (doing_hepi .or. doing_hekaon .or. doing_hedvcs) then VERTEXedge.Em.min = targ.Mtar_struck + targ.Mrec - targ.M !*IF CORRECT Mrec!!! ! VERTEXedge.Em.max = eval(nume) VERTEXedge.Em.max = 100.0 !temp. for generate_he_em VERTEXedge.Pm.min = 0.0 VERTEXedge.Pm.max = pval(nump) VERTEXedge.Trec.min = sqrt(targ.Mrec**2+VERTEXedge.Pm.min**2)-targ.Mrec VERTEXedge.Trec.max = sqrt(targ.Mrec**2+VERTEXedge.Pm.max**2)-targ.Mrec VERTEXedge.Trec_struck.min = 0.0 VERTEXedge.Trec_struck.max = Ebeam_max + targ.Mtar_struck - > targ.Mrec_struck - edge.e.E.min - edge.p.E.min Egamma_tot_max = Ebeam_max + targ.Mtar_struck - > targ.Mrec_struck - edge.e.E.min - edge.p.E.min - > VERTEXedge.Trec_struck.min - VERTEXedge.Em.min endif if (doing_kaon .or. doing_pion .or. doing_dvcs) then write(6,*) 'E_bind =',VERTEXedge.Em.min,'MeV in limits_init (QF only)' endif ! ... override calculated limits with hardwired value if desired. if (hardwired_rad) Egamma_tot_max = Egamma_gen_max if (.not.using_rad) Egamma_tot_max = 0.0 if (doing_tail(1)) Egamma1_max = Egamma_tot_max if (doing_tail(2)) Egamma2_max = Egamma_tot_max if (doing_tail(3)) Egamma3_max = Egamma_tot_max ! ... compute edge on summed _generated_ energies based on predicted VERTEX ! ... and TRUE values of Em (Ee'+Ep' for doing_heavy, Ee' for pion/kaon, ! ... Ee' for D(e,e'p), not used for hydrogen elastic.) if (doing_hyd_elast) then !NO generated energies. gen.sumEgen.min = 0.0 gen.sumEgen.max = 0.0 else if (doing_heavy) then !generated TOTAL (e+p) energy limits. gen.sumEgen.max = Ebeam_max + targ.Mtar_struck - VERTEXedge.Trec.min - VERTEXedge.Em.min gen.sumEgen.min = Ebeam_min + targ.Mtar_struck - VERTEXedge.Trec.max - VERTEXedge.Em.max - Egamma1_max ! ... Secondary limits for sumEgen when generating electron + proton. gen.sumEgen.max = min(gen.sumEgen.max,edge.e.E.max+edge.p.E.max+Egamma_tot_max) gen.sumEgen.min = max(gen.sumEgen.min,edge.e.E.min+edge.p.E.min) else !generated ELECTRON energy limits. if (doing_deuterium) then gen.sumEgen.max = Ebeam_max + targ.Mtar_struck - VERTEXedge.Em.min - > edge.p.E.min - VERTEXedge.Trec.min gen.sumEgen.min = Ebeam_min + targ.Mtar_struck - VERTEXedge.Em.max - > edge.p.E.max - VERTEXedge.Trec.max - Egamma_tot_max else if (doing_hydpi .or. doing_hydkaon .or. doing_hyddvcs) then gen.sumEgen.max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck - > edge.p.E.min - VERTEXedge.Trec_struck.min gen.sumEgen.min = Ebeam_min + targ.Mtar_struck - targ.Mrec_struck - > edge.p.E.max - VERTEXedge.Trec_struck.max - Egamma_tot_max else if (doing_deutpi.or.doing_hepi.or.doing_deutkaon.or. > doing_hekaon.or.doing_deutdvcs.or.doing_hedvcs) then gen.sumEgen.max = Ebeam_max + targ.Mtar_struck - targ.Mrec_struck - > edge.p.E.min - VERTEXedge.Trec_struck.min - VERTEXedge.Em.min - > VERTEXedge.Trec.min gen.sumEgen.min = Ebeam_min + targ.Mtar_struck - targ.Mrec_struck - > edge.p.E.max - VERTEXedge.Trec_struck.max - VERTEXedge.em.max - > VERTEXedge.Trec.max - Egamma_tot_max endif ! ... Secondary limits for sumEgen when only generating electron. gen.sumEgen.max = min(gen.sumEgen.max, edge.e.E.max+Egamma2_max) gen.sumEgen.min = max(gen.sumEgen.min, edge.e.E.min) endif gen.sumEgen.min = gen.sumEgen.min - dE_edge_test gen.sumEgen.max = gen.sumEgen.max + dE_edge_test gen.sumEgen.min = max(0.d0,gen.sumEgen.min) ! ... E arm GENERATION limits from sumEgen. ! ... Not used for doing_hyd_elast, but define for the hardwired histograms. if (doing_hyd_elast) then gen.e.E.min = edge.e.E.min gen.e.E.max = edge.e.E.max + Egamma2_max else if(doing_deuterium.or.doing_pion.or.doing_kaon.or.doing_dvcs) then gen.e.E.min = gen.sumEgen.min gen.e.E.max = gen.sumEgen.max else if (doing_heavy) then gen.e.E.min = gen.sumEgen.min - edge.p.E.max - Egamma3_max gen.e.E.max = gen.sumEgen.max - edge.p.E.min endif ! ... Apply limits from direct comparison to acceptance. gen.e.E.min = max(gen.e.E.min, edge.e.E.min) gen.e.E.max = min(gen.e.E.max, edge.e.E.max + Egamma2_max) gen.e.delta.min = (gen.e.E.min/spec.e.P-1.)*100. gen.e.delta.max = (gen.e.E.max/spec.e.P-1.)*100. gen.e.yptar.min = edge.e.yptar.min gen.e.yptar.max = edge.e.yptar.max gen.e.xptar.min = edge.e.xptar.min gen.e.xptar.max = edge.e.xptar.max ! ... P arm GENERATION limits from sumEgen. Not used for any case ! ... except doing_heavy, but need to define for code that writes out limits. if (doing_hyd_elast.or.doing_deuterium.or.doing_pion.or. > doing_kaon.or.doing_dvcs) then gen.p.E.min = edge.p.E.min gen.p.E.max = edge.p.E.max + Egamma3_max else if (doing_heavy)then gen.p.E.min = gen.sumEgen.min - edge.e.E.max - Egamma2_max gen.p.E.max = gen.sumEgen.max - edge.e.E.min endif ! ... Apply limits from direct comparison to acceptance. gen.p.E.min = max(gen.p.E.min, edge.p.E.min) gen.p.E.max = min(gen.p.E.max, edge.p.E.max + Egamma3_max) gen.p.delta.min = (sqrt(gen.p.E.min**2-Mh2)/spec.p.P-1.)*100. gen.p.delta.max = (sqrt(gen.p.E.max**2-Mh2)/spec.p.P-1.)*100. gen.p.yptar.min = edge.p.yptar.min gen.p.yptar.max = edge.p.yptar.max gen.p.xptar.min = edge.p.xptar.min gen.p.xptar.max = edge.p.xptar.max ! Axis specs for the diagnostic histograms H.gen.e.delta.min = gen.e.delta.min H.gen.e.yptar.min = gen.e.yptar.min H.gen.e.xptar.min = -gen.e.xptar.max H.gen.e.delta.bin = (gen.e.delta.max-gen.e.delta.min)/float(nHbins) H.gen.e.yptar.bin = (gen.e.yptar.max-gen.e.yptar.min)/float(nHbins) H.gen.e.xptar.bin = (gen.e.xptar.max-gen.e.xptar.min)/float(nHbins) H.gen.p.delta.min = gen.p.delta.min H.gen.p.yptar.min = gen.p.yptar.min H.gen.p.xptar.min = -gen.p.xptar.max H.gen.p.delta.bin = (gen.p.delta.max-gen.p.delta.min)/float(nHbins) H.gen.p.yptar.bin = (gen.p.yptar.max-gen.p.yptar.min)/float(nHbins) H.gen.p.xptar.bin = (gen.p.xptar.max-gen.p.xptar.min)/float(nHbins) H.gen.Em.min = VERTEXedge.Em.min H.gen.Em.bin = (max(100.d0,VERTEXedge.Em.max) - VERTEXedge.Em.min)/float(nHbins) H.gen.Pm.min = VERTEXedge.Pm.min H.gen.Pm.bin = (max(100.d0,VERTEXedge.Pm.max) - VERTEXedge.Pm.min)/float(nHbins) H.geni.e.delta.min = H.gen.e.delta.min H.geni.e.yptar.min = H.gen.e.yptar.min H.geni.e.xptar.min = H.gen.e.xptar.min H.geni.e.delta.bin = H.gen.e.delta.bin H.geni.e.yptar.bin = H.gen.e.yptar.bin H.geni.e.xptar.bin = H.gen.e.xptar.bin H.geni.p.delta.min = H.gen.p.delta.min H.geni.p.yptar.min = H.gen.p.yptar.min H.geni.p.xptar.min = H.gen.p.xptar.min H.geni.p.delta.bin = H.gen.p.delta.bin H.geni.p.yptar.bin = H.gen.p.yptar.bin H.geni.p.xptar.bin = H.gen.p.xptar.bin H.geni.Em.min = H.gen.Em.min H.geni.Em.bin = H.gen.Em.bin H.geni.Pm.min = H.gen.Pm.min H.geni.Pm.bin = H.gen.Pm.bin H.RECON.e.delta.min = H.gen.e.delta.min H.RECON.e.yptar.min = H.gen.e.yptar.min H.RECON.e.xptar.min = H.gen.e.xptar.min H.RECON.e.delta.bin = H.gen.e.delta.bin H.RECON.e.yptar.bin = H.gen.e.yptar.bin H.RECON.e.xptar.bin = H.gen.e.xptar.bin H.RECON.p.delta.min = H.gen.p.delta.min H.RECON.p.yptar.min = H.gen.p.yptar.min H.RECON.p.xptar.min = H.gen.p.xptar.min H.RECON.p.delta.bin = H.gen.p.delta.bin H.RECON.p.yptar.bin = H.gen.p.yptar.bin H.RECON.p.xptar.bin = H.gen.p.xptar.bin H.RECON.Em.min = H.gen.Em.min H.RECON.Em.bin = H.gen.Em.bin H.RECON.Pm.min = H.gen.Pm.min H.RECON.Pm.bin = H.gen.Pm.bin return end !------------------------------------------------------------------------ subroutine radc_init implicit none include 'simulate.inc' include 'radc.inc' include 'brem.inc' !-------------------------------------------------------------- ! ! 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 !-------------------------------------------------------------- ! Check setting of EXTRAD_FLAG if (debug(2)) write(6,*)'radc_init: entering...' if (extrad_flag.eq.0) then if (rad_flag.eq.0) then extrad_flag = 3 else if (rad_flag.eq.1 .or. rad_flag.eq.2 .or. rad_flag.eq.3) then extrad_flag = 1 endif else if (extrad_flag.lt.0) then stop 'Imbecile! check your stupid setting of EXTRAD_FLAG' endif ! 'etatzai' parameter if (debug(4)) write(6,*)'radc_init: at 1' etatzai = (12.0+(targ.Z+1.)/(targ.Z*targ.L1+targ.L2))/9.0 if (debug(4)) write(6,*)'radc_init: etatzai = ',etatzai ! Initialize brem flags (brem doesn't include the normal common blocks) produce_output = debug(1) exponentiate = use_expon include_hard = .true. calculate_spence = .true. if (debug(2)) write(6,*)'radc_init: ending...' return end !--------------------------------------------------------------------- subroutine radc_init_ev (main,vertex) implicit none include 'structures.inc' include 'radc.inc' integer i real*8 r, Ecutoff, dsoft, dhard, dsoft_prime real*8 lambda_dave, schwinger, brem, bremos record /event_main/ main record /event/ vertex ! Note that calculate_central calls this with (main,ev) rather than ! (main,vertex). Since these are just local variables, calling it vertex ! here does not cause any problem, and makes it easier to follow ! modifications to vertex.* variables in later calls. real*8 zero parameter (zero=0.0d0) !double precision zero for subroutine calls. ! Compute some quantities that will be needed for rad corr on this event ! ... factor for limiting energy of external radiation along incident electron ! etta = 1.0 + 2*vertex.ein*sin(vertex.e.theta/2.)**2/(targ.A*amu) ! ... moron move! let's can that etta factor ... etta = 1.0 ! ... the bt's do i=1,2 bt(i) = etatzai*main.target.teff(i) enddo ! ... the lambda's (effective bt's for internal radiation) do i=1,3 lambda(i) = lambda_dave(i,1,doing_tail(3),vertex.Ein,vertex.e.E,vertex.p.E, > vertex.p.P,vertex.e.theta) enddo rad_proton_this_ev = lambda(3).gt.0 ! ... get the hard correction factor. don't care about Ecutoff! Just want dhard here Ecutoff = 450. if (intcor_mode.eq.0) then r = schwinger(Ecutoff,vertex,.true.,dsoft,dhard) else if (.not.use_offshell_rad) then r = brem(vertex.Ein,vertex.e.E,Ecutoff,rad_proton_this_ev,dsoft,dhard, > dsoft_prime) else r = bremos(Ecutoff, zero, zero, vertex.Ein, vertex.e.P*vertex.ue.x, > vertex.e.P*vertex.ue.y, vertex.e.P*vertex.ue.z, 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, dhard, dsoft_prime) endif endif hardcorfac = 1./(1.-dhard) g(4)=-dsoft_prime*Ecutoff+bt(1)+bt(2) ! ... initialize the parameters needed for our "basic" calculation call basicrad_init_ev (vertex.Ein,vertex.e.E,vertex.p.E) ! ... the relative magnitudes of the three tails (we may not need them) do i=1,3 frac(i) = g(i)/g(0) enddo return end !----------------------------------------------------------------------- subroutine basicrad_init_ev (e1,e2,e3) implicit none include 'simulate.inc' include 'radc.inc' real*8 one parameter (one=1.) integer i real*8 e1,e2,e3,e(3),gamma if (debug(2)) write(6,*)'basicrad_init_ev: entering...' e(1) = e1 e(2) = e2 e(3) = e3 ! bt's for internal + external ! ??? One possibility for shutting off tails 1 or ! 2 is to set g(1/2) = 0 here ... note that lambda(3) is set to 0 in ! lambda_dave at the moment, AND proton terms are removed from brem if proton ! radiation off. something analogous and similarly consistent would have to be ! done for the other tails, right now they're just nixed in generate_rad. also ! check ALL ntail.eq.0 checks in kinema constraint lines of generate_rad g(1) = lambda(1) + bt(1) g(2) = lambda(2) + bt(2) g(3) = lambda(3) g(0) = g(1)+g(2)+g(3) ! Internal constants c_int(1) = lambda(1)/(e(1)*e(2))**(lambda(1)/2.) c_int(2) = lambda(2)/(e(1)*e(2))**(lambda(2)/2.) * csa NOTE: GN says these masses s/b Mp! c_int(3) = lambda(3)/(Mp*e(3))**(lambda(3)/2.) do i = 1, 3 c_int(i) = c_int(i) * exp(-euler*lambda(i)) / gamma(one+lambda(i)) enddo g_int = lambda(1) + lambda(2) + lambda(3) c_int(0) = c_int(1)*c_int(2) * g_int / lambda(1)/lambda(2) ! ... proton radiation could be off if (lambda(3).gt.0) c_int(0) = c_int(0) * c_int(3)/lambda(3) c_int(0) = c_int(0) * gamma(one+lambda(1)) * gamma(one+lambda(2)) > * gamma(one+lambda(3)) / gamma(one+g_int) ! External constants do i = 1, 2 c_ext(i) = bt(i)/e(i)**bt(i)/gamma(one+bt(i)) enddo c_ext(3) = 0.0 g_ext = bt(1) + bt(2) c_ext(0) = c_ext(1)*c_ext(2) * g_ext / bt(1)/bt(2) c_ext(0) = c_ext(0)*gamma(one+bt(1))*gamma(one+bt(2))/gamma(one+g_ext) ! Internal + external constants do i = 1, 2 c(i) = c_int(i) * c_ext(i) * g(i)/lambda(i)/bt(i) > * gamma(one+lambda(i))*gamma(one+bt(i))/gamma(one+g(i)) enddo c(3) = c_int(3) ! Finally, constant for combined tails c(0) = c(1)*c(2) * g(0)/g(1)/g(2) ! ... proton radiation could be off if (g(3).gt.0) c(0) = c(0) * c(3)/g(3) c(0)=c(0)*gamma(one+g(1))*gamma(one+g(2))*gamma(one+g(3))/gamma(one+g(0)) c(4)=g(4)/(e1*e2)**g(4)/gamma(one+g(4)) if(g(3).gt.0) c(4)=c(4)/e3**g(4) if (debug(2)) write(6,*)'basicrad_init_ev: ending...' return end !---------------------------------------------------------------------- ! theory_init: ! ! Input spectral function(NOT called for H(e,e'p) or pion/kaon production!) ! ! Note: Some flags that existed in old A(e,e'p) code (e.g. DD's version) ! have been removed. Code has been changed to correspond to the following ! settings for the old flags: ! norm_Ji_tail = 0 ! doing_Ji_tail = 0 ! Em_start_Ji_tail = 0. ! fixed_gamma = .true. subroutine theory_init(success) implicit none include 'simulate.inc' real*8 Pm_values(ntheorymax), absorption integer m,n,iok logical success ! ... open the file if ( nint(targ.A) .eq. 2) then theory_file='h2.theory' else if ( nint(targ.A) .eq. 12) then theory_file='c12.theory' else if ( nint(targ.A) .eq. 56) then theory_file='fe56.theory' else if ( nint(targ.A) .eq. 197) then theory_file='au197.theory' else write(6,*) 'No Theory File (spectral function) for A = ',targ.A write(6,*) 'Defaulting to c12.theory' theory_file='c12.theory' endif open(unit=1,file=theory_file,status='old',readonly,shared,iostat=iok) ! ... read in the theory file read(1,*,err=40) nrhoPm, absorption, E_Fermi do m=1, nrhoPm read(1,*,err=40) nprot_theory(m), Em_theory(m), Emsig_theory(m), > bs_norm_theory(m) nprot_theory(m) = nprot_theory(m) * absorption enddo read(1,*,err=40) Pm_values(1),theory(1,1) do m=1, nrhoPm n=2 read(1,*,err=40,end=50) Pm_values(2),theory(m,2) do while (Pm_values(n).gt.Pm_values(n-1)) n=n+1 read(1,*,err=40,end=50) Pm_values(n),theory(m,n) enddo ! ........ figure out details of the Pm axes 50 Pm_theory(m).n=n-1 Pm_theory(m).bin=Pm_values(2)-Pm_values(1) Pm_theory(m).min=Pm_values(1)-Pm_theory(m).bin/2. Pm_theory(m).max=Pm_values(Pm_theory(m).n)+Pm_theory(m).bin/2. if (abs(Pm_theory(m).min+Pm_theory(m).bin*Pm_theory(m).n - 1 Pm_theory(m).max) .gt. 0.1) then write(6,'(1x,''ERROR: theory_init found unequal Pm bins in distribution number '',i2,''!'')') m close(1) return endif ! ........ prepare for the next loop Pm_values(1) = Pm_values(n) theory(m+1,1) = theory(m,n) enddo ! ... are we doing deuterium? (i.e. only using a 1D spectral function) doing_deuterium = nrhoPm.eq.1 .and. E_Fermi.lt.1.0 ! ... renormalize the momentum distributions do m=1, nrhoPm do n=1, Pm_theory(m).n theory(m,n) = theory(m,n)/bs_norm_theory(m) enddo enddo ! ... and calculate the integral of the Em distribution above E_Fermi to ! ... renormalize if (doing_heavy) then do m=1, nrhoPm Em_int_theory(m) = 1. Em_int_theory(m) = (pi/2. + atan((Em_theory(m)-E_Fermi)/ > (0.5*Emsig_theory(m))))/pi enddo endif ! ... we made it success=.true. close(1) return ! ... oops 40 continue write(6,*) 'ERROR: theory_init failed to read in the theory file' return end