subroutine limits_update(main,vertex,orig,recon,doing_deuterium, > doing_pion,doing_kaon,doing_dvcs,contrib,slop) include 'structures.inc' include 'radc.inc' record /event_main/ main record /event/ vertex, orig, recon record /contrib/ contrib record /slop/ slop integer i logical doing_deuterium, doing_pion, doing_kaon, doing_dvcs ! Update the "contribution limits" records ! ... GENERATION values call update_range(vertex.e.delta, contrib.gen.e.delta) call update_range(vertex.e.yptar, contrib.gen.e.yptar) call update_range(vertex.e.xptar, contrib.gen.e.xptar) call update_range(vertex.p.delta, contrib.gen.p.delta) call update_range(vertex.p.yptar, contrib.gen.p.yptar) call update_range(vertex.p.xptar, contrib.gen.p.xptar) call update_range(main.Trec, contrib.gen.Trec) ! ........ another tricky shift if (doing_deuterium.or.doing_pion.or.doing_kaon.or.doing_dvcs) then call update_range(vertex.e.E-main.Ein_shift,contrib.gen.sumEgen) else call update_range(vertex.e.E+vertex.p.E-main.Ein_shift,contrib.gen.sumEgen) endif ! ... TRUE values ! ........ tricky shift here, remember this'll get compared with edge.e.E, ! ........ compensate for main.target.Coulomb to ! ........ copy the shift made to the edge in generate call update_range(orig.e.E-main.Ee_shift, contrib.tru.e.E) call update_range(orig.e.xptar, contrib.tru.e.xptar) call update_range(orig.e.yptar, contrib.tru.e.yptar) call update_range(orig.e.xptar, contrib.tru.e.xptar) call update_range(orig.p.E, contrib.tru.p.E) call update_range(orig.p.yptar, contrib.tru.p.yptar) call update_range(orig.p.xptar, contrib.tru.p.xptar) ! ........ another tricky shift call update_range(orig.Em-main.Ein_shift+main.Ee_shift,contrib.tru.Em) call update_range(orig.Pm, contrib.tru.Pm) call update_range(orig.Trec, contrib.tru.Trec) ! ... SPECTROMETER values call update_range(main.SP.e.delta, contrib.SP.e.delta) call update_range(main.SP.e.yptar, contrib.SP.e.yptar) call update_range(main.SP.e.xptar, contrib.SP.e.xptar) call update_range(main.SP.p.delta, contrib.SP.p.delta) call update_range(main.SP.p.yptar, contrib.SP.p.yptar) call update_range(main.SP.p.xptar, contrib.SP.p.xptar) ! ... VERTEX values call update_range(vertex.Trec, contrib.vertex.Trec) call update_range(vertex.Em, contrib.vertex.Em) call update_range(vertex.Pm, contrib.vertex.Pm) ! ... RADIATION stuff ! ??? should be looking at Egamma2+3 cause we do use limits on that, indirectly do i = 1, 3 call update_range(Egamma_used(i), contrib.rad.Egamma(i)) enddo call update_range(Egamma_used(1)+Egamma_used(2)+Egamma_used(3), > contrib.rad.Egamma_total) ! Update the "slop limits" records ! ... MC slops call update_range(main.RECON.e.delta-main.SP.e.delta,slop.MC.e.delta) call update_range(main.RECON.e.yptar-main.SP.e.yptar,slop.MC.e.yptar) call update_range(main.RECON.e.xptar-main.SP.e.xptar,slop.MC.e.xptar) call update_range(main.RECON.p.delta-main.SP.p.delta,slop.MC.p.delta) call update_range(main.RECON.p.yptar-main.SP.p.yptar,slop.MC.p.yptar) call update_range(main.RECON.p.xptar-main.SP.p.xptar,slop.MC.p.xptar) ! ... total slops ! ........ that tricky shift again, slops accounted for by the shift not ! ........ included in slop.total.Em. call update_range(recon.Em-(orig.Em-main.Ein_shift+main.Ee_shift), > slop.total.Em) call update_range(abs(recon.Pm)-abs(orig.Pm), slop.total.Pm) return end !------------------------------------------------------------------- subroutine update_range(val,range) include 'structures.inc' record /range/ range real*8 val range.lo = min(range.lo, val) range.hi = max(range.hi, val) return end !---------------------------------------------------------------------- ! THE routine to GENERATE the (max of 7) random quantities we need to get ! a full description of our event. subroutine generate(main,vertex,orig,success) implicit none include 'simulate.inc' integer i, ii real*8 Emin, Emax real*8 ranprob,ranth1,ranth,ranph,t3,t4,t5,t6 real*8 pferlo,pferhi real*8 gauss1 logical success real*8 grnd !random # generator. record /event_main/ main record /event/ vertex, orig real*8 nsig_max parameter(nsig_max=3.0d0) !max #/sigma for gaussian ran #s. ! Randomize the position of the interaction inside the available region. ! gen.xwid and gen.ywid are the intrinsic beam widths (one sigma value). ! Use a gaussian beam distribution with +/- 3.0 sigma (add raster afterwards). if (debug(2)) write(6,*)'gen: entering' main.target.x = gauss1(nsig_max)*gen.xwid+targ.xoffset main.target.y = gauss1(nsig_max)*gen.ywid+targ.yoffset ! fr_pattern=1 - rectangle - targ.fr1/fr2 are the x/y raster half-widths. ! fr_pattern=2 - circle - targ.fr1/fr2 are the inner and outer radii. if (targ.fr_pattern .eq. 1) then !square raster t3=grnd()*pi t4=grnd()*pi t5=cos(t3)*targ.fr1 t6=cos(t4)*targ.fr2 elseif (targ.fr_pattern .eq. 2) then !circular raster t3=grnd()*2.*pi t4=sqrt(grnd())*(targ.fr2-targ.fr1)+targ.fr1 t5=cos(t3)*t4 t6=sin(t3)*t4 else !no raster t5=0.0 t6=0.0 endif main.target.x = main.target.x+t5 main.target.y = main.target.y+t6 main.target.z = (0.5-grnd())*targ.length+targ.zoffset main.target.rastery = t6 !'raster' contribution to vert. pos. main.target.rasterx = t5 !points left as look downstream. GAW ! Take fluctuations of the beam energy into account, and remember to correct ! for ionization loss in the target and Coulomb acceleration of incoming ! electron. Remove targ.zoffset from the z position of the scattering ! in order to get the position relative to the center of the target. call trip_thru_target (1, main.target.z-targ.zoffset, Ebeam, 0.0d0, > main.target.Eloss(1), main.target.teff(1),Me,1) if (.not.using_Eloss) main.target.Eloss(1) = 0.0 if (using_Coulomb) then main.target.Coulomb=targ.Coulomb_constant*(3.-(grnd())**(2./3.)) else main.target.Coulomb=0.0 endif vertex.Ein = Ebeam + (grnd()-0.5)*dEbeam + > main.target.Coulomb - main.target.Eloss(1) ! ... deterimine known variation in Ein from Ebeam_vertex_ave and in ! ... Ee due to Coulomb energy to compare limits to generated event by event. ! ... (USED TO SHIFT 'LIMIT' VALUES IN UPDATE_RANGE CALLS. ARE THEY NEEDED???) main.Ein_shift = vertex.Ein - Ebeam_vertex_ave main.Ee_shift = main.target.Coulomb - targ.Coulomb.ave ! Initialize success code and fractional weight due to radiation and ! generation tricks 5 success = .false. main.gen_weight = 1.0 ! Generated quantities: (phase_space NOT YET IMPLEMENTED). ! ! phase_space: Generate electron E,yptar,xptar and hadron yptar,xptar?? ! doing_hyd_elast: fixed Em, generate electron angles ! doing_deuterium: fixed Em, generate electron fully and proton angles, calc Ep ! doing_eep, A>2: generate electron and hadron energy and angles (calc Em/Pm). ! doing_pion: generate electron energy/angles, hadron angles, p_fermi, Em. ! doing_kaon: as doing_pion. ! ! The above is summarized in the following table: ! ! ELECTRON HADRON ! ------------------ ------------------ ! E yptar xptar E yptar xptar p_fermi Em ! !H(e,e'p) X X !D(e,e'p) X X X X X !A(e,e'p) X X X X X X !---------------------------------------------------------------------- !H(e,e'pi) X X X X X !A(e,e'pi) X X X X X X X !H(e,e'K) X X X X X !A(e,e'K) X X X X X X X !---------------------------------------------------------------------- !phase_space X X X ? X X ! ! So our procedure is the following: ! 1) Always generate electron yptar and xptar ! 2) generate hadron yptar and xptar for all cases except H(e,e'p), D(e,e'p) ! 3) Generate hadron E for A(e,e'p) ! 4) generate p_fermi, Em for A(e,e'pi) and A(e,e'K) - calculate vertex.Mrec ! (invariant mass M* of the A-1 recoil system) ! 5) generate electron E for all but hydrogen elastic. ! ! After we generate xptar/yptar/energy, we calculate physics angles (theta/phi), ! momenta, unit vectors, etc... here and/or in complete_ev. ! ! Note that there are also jacobians associated with some and/or all of ! the above. ! 1: We generate uniformly in xptar/yptar, not theta/phi. We define the ! phase space volume (genvol contribution) as the product of the xptar/yptar ! range, and have a jacobian for each event taking into account the mapping ! between the solid angle on the unit sphere, and the dxptar/dyptar volume ! (the jacobian is 1/cos**3(dtheta), where dtheta is the angle between the ! event and the central spectrometer vector ! 2: For the D(e,e'p), we take Em as fixed in order to calculate the proton ! energy. There is a jacobian ( |dEp'/dEm| ). It comes from integrating ! over the energy conservation delta function: delta(E_D - E_p - E_n - Em). ! Generate Electron Angles (all cases): vertex.e.yptar=gen.e.yptar.min+grnd()*(gen.e.yptar.max-gen.e.yptar.min) vertex.e.xptar=gen.e.xptar.min+grnd()*(gen.e.xptar.max-gen.e.xptar.min) ! Generate Hadron Angles (all but H(e,e'p)): if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon.or.doing_dvcs) then vertex.p.yptar=gen.p.yptar.min+grnd()*(gen.p.yptar.max-gen.p.yptar.min) vertex.p.xptar=gen.p.xptar.min+grnd()*(gen.p.xptar.max-gen.p.xptar.min) endif ! Generate Hadron Momentum (A(e,e'p) only). if (doing_heavy) then Emin = max(gen.p.E.min, gen.sumEgen.min - gen.e.E.max) Emax = min(gen.p.E.max, gen.sumEgen.max - gen.e.E.min) if (Emin.gt.Emax) goto 100 main.gen_weight=main.gen_weight*(Emax-Emin)/(gen.p.E.max-gen.p.E.min) vertex.p.E = Emin + grnd()*(Emax-Emin) vertex.p.P = sqrt(vertex.p.E**2 - Mh2) vertex.p.delta = 100.*(vertex.p.P-spec.p.P)/spec.p.P endif ! Generate Fermi Momentum and Em for A(e,e'pi) and A(e,e'K). pfer=0.0 pferx=0.0 pfery=0.0 pferz=0.0 vertex.Em=0.0 if(doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon.or.doing_hedvcs)then ranprob=grnd() ii=1 do while (ranprob.gt.mprob(ii) .and. ii.lt.nump) ii=ii+1 enddo if (ii.eq.1) then pferlo=0 else pferlo=(pval(ii-1)+pval(ii))/2 endif if (ii.eq.nump) then pferhi=pval(nump) else pferhi=(pval(ii)+pval(ii+1))/2 endif pfer=pferlo+(pferhi-pferlo)*grnd() ranth1=grnd()*2.-1.0 ranth=acos(ranth1) ranph=grnd()*2.*pi pferx=sin(ranth)*cos(ranph) pfery=sin(ranth)*sin(ranph) pferz=cos(ranth) if(doing_deutpi.or.doing_deutkaon.or.doing_hedvcs)then !Em = binding energy vertex.Em = Mp + Mn - targ.M endif if (doing_hepi .or. doing_hekaon.or.doing_hedvcs) then!Generate Em call generate_em(pfer,vertex.Em) !Generate Em endif endif ! Generate Electron Energy (all but hydrogen elastic) if (doing_deuterium.or.doing_heavy.or.doing_pion.or.doing_kaon.or. 1 doing_dvcs) then Emin=gen.e.E.min Emax=gen.e.E.max if (doing_deuterium .or. doing_pion .or. doing_kaon .or. 1 doing_dvcs) then Emin = max(Emin,gen.sumEgen.min) Emax = min(Emax,gen.sumEgen.max) else if (doing_heavy) then ! A(e,e'p) Emin = max(Emin, gen.sumEgen.min - vertex.p.E) Emax = min(Emax, gen.sumEgen.max - vertex.p.E) endif if (Emin.gt.Emax) goto 100 main.gen_weight=main.gen_weight*(Emax-Emin)/(gen.e.E.max-gen.e.E.min) vertex.e.E = Emin + grnd()*(Emax-Emin) vertex.e.P = vertex.e.E vertex.e.delta = 100.*(vertex.e.P-spec.e.P)/spec.e.P endif !not (doing_hyd_elast) ! Calculate the electron and proton PHYSICS angles from the spectrometer angles. call physics_angles(spec.e.theta,spec.e.phi, & vertex.e.xptar,vertex.e.yptar,vertex.e.theta,vertex.e.phi) call physics_angles(spec.p.theta,spec.p.phi, & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi) ! Compute all non-generated quantities if (debug(5)) write(6,*)'gen: calling comp_ev with false, main, vertex' if (debug(3)) write(6,*)'gen: calling comp_ev with false, main, vertex' if (debug(3)) write(6,*)'gen: Ein, E =',vertex.Ein,vertex.e.E call complete_ev(main,vertex,success) main.sigcc = 1.0 if (debug(2)) write(6,*)'gen: initial success =',success if (.not.success) goto 100 ! ........ temporary storage of Trec for generated event main.Trec = vertex.Trec ! Radiate the event, if requested. If we get an event weight of zero (i.e. ! if we CAN'T radiate our event into the acceptance) then success is ! false. generate_rad also set orig kinematics, and cuts on Em/Pm. ! If not using_rad, must do these here. if (using_rad) then call generate_rad(main,vertex,orig,success) if (debug(2)) write(6,*)'gen: after gen_rad, success =',success else success = .true. if (doing_heavy) success = > (vertex.Em .ge. VERTEXedge.Em.min .and. > vertex.Em .le. VERTEXedge.Em.max .and. > vertex.Pm .ge. VERTEXedge.Pm.min .and. > vertex.Pm .le. VERTEXedge.Pm.max) if (success) then do i = 1, neventfields orig.all(i) = vertex.all(i) enddo endif endif 100 if (debug(2)) write(6,*)'gen: final success =',success if (debug(2)) write(6,*)'gen: ending' return end !--------------------------------------------------------------------- subroutine complete_ev(main,vertex,success) implicit none include 'simulate.inc' real*8 a, b, c, r, t, QA, QB, QC, radical real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z real*8 targ_new_x,targ_new_y real*8 new_y_x,new_y_y,new_y_z,dummy real*8 px,py,pz,qx,qy,qz real*8 targx,targy,targz,pol real*8 cos_dtheta,y_spec,y_event real*8 oop_x,oop_y real*8 krel,krelx,krely,krelz real*8 MM logical success record /event_main/ main record /event/ vertex !----------------------------------------------------------------------- ! Calculate everything left in the /event/ structure, given all necessary ! GENERATION values (some set of xptar,yptar,delta for both arms and p_fermi, ! and p_fermi, depending on the scattering process: see table in generate.f ! ! The SINGLE element of /event/ NOT computed here is sigcc. ! ! Another small anomaly is that main.jacobian IS computed here (if not ! doing_recon) ... this is because all the necessary terms have to be ! computed here anyway to calculate /event/ qties. ! !----------------------------------------------------------------------- ! Initialize success = .false. main.jacobian = 1.0 ! ... unit vector components of outgoing e,p ! ... z is DOWNSTREAM, x is DOWN and y is LEFT looking downstream. if (debug(4)) write(6,*)'comp_ev: at 1' vertex.ue.x = sin(vertex.e.theta)*cos(vertex.e.phi) vertex.ue.y = sin(vertex.e.theta)*sin(vertex.e.phi) vertex.ue.z = cos(vertex.e.theta) if (.not.doing_hyd_elast) then vertex.up.x = sin(vertex.p.theta)*cos(vertex.p.phi) vertex.up.y = sin(vertex.p.theta)*sin(vertex.p.phi) vertex.up.z = cos(vertex.p.theta) endif if (debug(4)) write(6,*)'comp_ev: at 2' ! First finish off the e side ! Calculate scattered electron energy for hydrogen/deuterium (e,e'p) if (doing_hyd_elast) then vertex.e.E = vertex.Ein*Mh/(Mh+vertex.Ein*(1.-vertex.ue.z)) if (vertex.e.E.gt.vertex.Ein) return vertex.e.P = vertex.e.E vertex.e.delta = (vertex.e.P - spec.e.P)*100./spec.e.P if (debug(4)) write(6,*)'comp_ev: at 3' endif ! The q vector if (debug(5)) write(6,*)'comp_ev: Ein,E,uez=',vertex.Ein,vertex.e.E,vertex.ue.z vertex.omega = vertex.Ein - vertex.e.E vertex.Q2 = 2*vertex.Ein*vertex.e.E*(1.-vertex.ue.z) vertex.q = sqrt(vertex.Q2 + vertex.omega**2) vertex.uq.x = - vertex.e.P*vertex.ue.x / vertex.q vertex.uq.y = - vertex.e.P*vertex.ue.y / vertex.q vertex.uq.z =(vertex.Ein - vertex.e.P*vertex.ue.z)/ vertex.q if (abs(vertex.uq.x**2+vertex.uq.y**2+vertex.uq.z**2-1).gt.0.01) > stop 'Error in q vector normalization' if (debug(4)) write(6,*)'comp_ev: at 5' ! Determine vertex.Pm, vertex.Em, vertex.Mrec (needed to complete hadron). ! For D(e,e'p) and A(e,e'p), vertex.Pm has to wait until after getting hadron. ! For A(e,e'p) case, Em and Mrec have to wait until after calculating Pm. ! For electroproduction, Em was generated at same time as pfermi (=0 for A=1). if (doing_hyd_elast) then vertex.Em = 0.0 vertex.Pm = 0.0 vertex.Mrec = 0.0 else if (doing_deuterium) then vertex.Em = Mp + Mn - targ.M !2.2 MeV vertex.Mrec = targ.M - targ.Mtar_struck + vertex.Em !neutron mass. else if (doing_pion .or. doing_kaon .or. doing_dvcs) then vertex.Pm = pfer vertex.Mrec = targ.M - targ.Mtar_struck + vertex.Em endif ! Now complete the p side. if (doing_hyd_elast) then !p = q vertex.up.x = vertex.uq.x vertex.up.y = vertex.uq.y vertex.up.z = vertex.uq.z vertex.p.P = vertex.q vertex.p.theta = acos(vertex.up.z) if (sin(vertex.p.theta).eq.0) write(6,*) 'sin(vertex.p.theta)=',sin(vertex.p.theta) if (abs(vertex.up.x/sin(vertex.p.theta)).gt.1) > write(6,*) 'cos(phi)=',vertex.up.x/sin(vertex.p.theta) vertex.p.phi = acos(min(1.d0,max(-1.d0,vertex.up.x/sin(vertex.p.theta)))) !...Get the cosine of the angle between the central spectrometer setting !...and the real event (dot product of p and p0 unit vectors). For in-plane !...spectrometers, xptar (=dx/dz) is just the x component of p (=dx) !...over the component of p along the spectrometer direction (=dz=cos(dtheta)). !...Knowing xptar, get yptar from cos(dtheta)=1/sqrt(1+xptar**2+yptar**2) !...Choose sign of yptar by comparing y for the spectrometer to y of the !...event (projected to the xptar-yptar plane). cos_dtheta = vertex.up.z * cos(spec.p.theta) > + vertex.up.y * sin(spec.p.theta)*sin(spec.p.phi) > + vertex.up.x * sin(spec.p.theta)*cos(spec.p.phi) vertex.p.xptar = vertex.up.x / cos_dtheta if ((1/cos_dtheta**2 - (1 + vertex.p.xptar**2)).lt.0) write (6,*) > 'complete_ev: yptar**2=',(1/cos_dtheta**2 - (1 + vertex.p.xptar**2)) vertex.p.yptar = sqrt(max(0.d0,1/cos_dtheta**2-(1+vertex.p.xptar**2))) y_spec = sin(spec.p.theta)*sin(spec.p.phi) y_event = vertex.up.y/cos_dtheta !projected to plane perp. to spectrometer. if (y_event .lt. y_spec) vertex.p.yptar = -vertex.p.yptar vertex.p.E = sqrt(vertex.p.P**2+Mh2) vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P if (debug(4)) write(6,*)'comp_ev: at 6' elseif (doing_deuterium) then !need Ep, and a jacobian. a = -1.*vertex.q*(vertex.uq.x*vertex.up.x+vertex.uq.y*vertex.up.y+vertex.uq.z*vertex.up.z) b = vertex.q**2 c = vertex.omega + targ.M t = c**2 - b + Mh2 - vertex.Mrec**2 QA = 4.*(a**2 - c**2) QB = 4.*c*t QC = -4.*a**2*Mh2 - t**2 radical = QB**2 - 4.*QA*QC if (radical.lt.0) return vertex.p.E = (-QB - sqrt(radical))/2./QA if (vertex.p.E.le.Mh) return vertex.p.P = sqrt(vertex.p.E**2 - Mh2) vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P ! ........ the Jacobian here is |dEp'/dEm| main.jacobian = (t*(c-vertex.p.E) + 2*c*vertex.p.E*(vertex.p.E-c)) / 1 (2*(a**2-c**2)*vertex.p.E + c*t) main.jacobian = abs(main.jacobian) elseif (doing_pion .or. doing_kaon .or. doing_dvcs) then ! Hydrogen case. Coherent production (bound final state) is treated as ! hydrogen, but with targ.Mtar_struck=targ.M, targ.Mtar_rec=bound final state. a = -1.*vertex.q*(vertex.uq.x*vertex.up.x+vertex.uq.y*vertex.up.y+vertex.uq.z*vertex.up.z) b = vertex.q**2 c = vertex.omega + targ.M ! For nuclei, correct for fermi motion and missing energy. if (doing_deutpi.or.doing_hepi.or.doing_deutkaon.or.doing_hekaon 1 .or.doing_deutdvcs.or.doing_hedvcs) then a = a - abs(pfer)*(pferx*vertex.up.x+pfery*vertex.up.y+pferz*vertex.up.z) b = b + pfer**2 + 2*vertex.q*abs(pfer)* > (pferx*vertex.uq.x+pfery*vertex.uq.y+pferz*vertex.uq.z) c = c - sqrt(vertex.Mrec**2+pfer**2) endif t = c**2 - b + Mh2 - targ.Mrec_struck**2 QA = 4.*(a**2 - c**2) QB = 4.*c*t QC = -4.*a**2*Mh2 - t**2 radical = QB**2 - 4.*QA*QC if (radical.lt.0) return vertex.p.E = (-QB - sqrt(radical))/2./QA if (vertex.p.E.le.Mh) return vertex.p.P = sqrt(vertex.p.E**2 - Mh2) vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P elseif (doing_phsp) then vertex.p.P = spec.p.P !????? single arm phsp?? vertex.p.E= sqrt(Mh2+vertex.p.P**2) vertex.p.delta = (vertex.p.P - spec.p.P)*100./spec.p.P if (debug(4)) write(6,*)'comp_ev: at 7.5',Mh2,vertex.p.E endif if (debug(4)) write(6,*)'comp_ev: at 7' ! Compute some pion and kaon stuff if (doing_pion .or. doing_kaon.or.doing_dvcs) then main.epsilon=1./(1. + 2.*(1+vertex.omega**2/vertex.Q2)*tan(vertex.e.theta/2.)**2) main.theta_pq=acos(vertex.up.x*vertex.uq.x+vertex.up.y*vertex.uq.y+vertex.up.z*vertex.uq.z) ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE. ! Therefore, define a new system with the z axis parallel to q, and ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane ! particles closer to the downstream beamline than q. ! phi=90 is above the horizontal plane when q points to the right, and ! below the horizontal plane when q points to the left. ! Also take into account the different definitions of x, y and z in ! replay and SIMC: ! As seen looking downstream: replay SIMC (old_simc) ! x right down (left) ! y down left (up) ! z all have z pointing downstream ! ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc qx = -vertex.uq.y !convert to 'replay' coord. system qy = vertex.uq.x qz = vertex.uq.z px = -vertex.up.y py = vertex.up.x pz = vertex.up.z dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2)) new_x_x = -qx*qz/dummy new_x_y = -qy*qz/dummy new_x_z = (qx**2 + qy**2)/dummy dummy = sqrt(qx**2 + qy**2) new_y_x = qy/dummy new_y_y = -qx/dummy new_y_z = 0.0 p_new_x = px*new_x_x + py*new_x_y + pz*new_x_z p_new_y = px*new_y_x + py*new_y_y + pz*new_y_z main.phi_pq = atan2(p_new_y,p_new_x) !atan2(y,x)=atan(y/x) if(main.phi_pq .lt. 0.) main.phi_pq = main.phi_pq + 2.*pi ! if ((p_new_x**2+p_new_y**2).eq.0.) then ! main.phi_pq = 0.0 ! else ! main.phi_pq = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2)) ! endif ! if (p_new_y.lt.0.) then ! main.phi_pq = 2*pi - main.phi_pq ! endif ! if (debug(2)) then write(6,*)'comp_ev: omega =',vertex.omega write(6,*)'comp_ev: Q2 =',vertex.Q2 write(6,*)'comp_ev: theta_e =',vertex.e.theta write(6,*)'comp_ev: epsilon =',main.epsilon write(6,*)'comp_ev: theta_pq =',main.theta_pq write(6,*)'comp_ev: phi_pq =',main.phi_pq endif endif !end of pion/kaon specific stuff. ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND TARGET POLARIZATION. ! Therefore, define a new system with the z axis parallel to q, and ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane ! particles closer to the downstream beamline than q. ! phi=90 is above the horizontal plane when q points to the right, and ! below the horizontal plane when q points to the left. ! Also take into account the different definitions of x, y and z in ! replay and SIMC: ! As seen looking downstream: replay SIMC (old_simc) ! x right down (left) ! y down left (up) ! z all have z pointing downstream ! ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc qx = -vertex.uq.y !convert to 'replay' coord. system qy = vertex.uq.x qz = vertex.uq.z targx = -targ_pol*sin(targ_Bangle) ! replay coordinates targy = 0.0 targz = targ_pol*cos(targ_Bangle) dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2)) new_x_x = -qx*qz/dummy new_x_y = -qy*qz/dummy new_x_z = (qx**2 + qy**2)/dummy dummy = sqrt(qx**2 + qy**2) new_y_x = qy/dummy new_y_y = -qx/dummy new_y_z = 0.0 p_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z p_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z main.phi_targ = atan2(p_new_y,p_new_x) !atan2(y,x)=atan(y/x) if(main.phi_targ.lt.0.) main.phi_targ = 2.*pi+main.phi_targ ! CALCULATE ANGLE BETA BETWEEN REACTION PLANE AND TRANSVERSE TARGET ! POLARIZATION. ! Therefore, define a new system with the z axis parallel to q, and ! with the y-direction defined by q x p_meson and x-axis given ! by y x z ! Also take into account the different definitions of x, y and z in ! replay and SIMC: ! As seen looking downstream: replay SIMC (old_simc) ! x right down (left) ! y down left (up) ! z all have z pointing downstream ! ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc qx = -vertex.uq.y !convert to 'replay' coord. system qy = vertex.uq.x qz = vertex.uq.z pol = 1.0 targx = -targ_pol*sin(targ_Bangle) ! 'replay' coordinates targy = 0.0 targz = targ_pol*cos(targ_Bangle) px = -vertex.up.y py = vertex.up.x pz = vertex.up.z dummy = sqrt((qy*pz-qz*py)**2 + (qz*px-qx*pz)**2 + (qx*py-qy*px)**2) new_y_x = (qy*pz-qz*py)/dummy new_y_y = (qz*px-qx*pz)/dummy new_y_z = (qx*py-qy*px)/dummy dummy = sqrt((new_y_y*qz-new_y_z*qy)**2 + (new_y_z*qx-new_y_x*qz)**2 1 + (new_y_x*qy-new_y_y*qx)**2) new_x_x = (new_y_y*qz - new_y_z*qy)/dummy new_x_y = (new_y_z*qx - new_y_x*qz)/dummy new_x_z = (new_y_x*qy - new_y_y*qx)/dummy targ_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z targ_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z main.beta = atan2(targ_new_y,targ_new_x) if(main.beta .lt. 0.) main.beta = 2*pi+main.beta c if ((targ_new_x**2+targ_new_y**2).eq.0.) then c main.beta = 0.0 c else c main.beta = acos(targ_new_x/sqrt(targ_new_x**2+targ_new_y**2)) c endif c if (targ_new_y.lt.0.) then c main.beta = 2*3.1415926536 - main.beta c endif dummy = sqrt((qx**2+qy**2+qz**2))*sqrt((targx**2+targy**2+targz**2)) main.theta_tarq = acos((qx*targx+qy*targy+qz*targz)/dummy) if (debug(4)) write(6,*)'comp_ev: at 8' ! Compute the Pm vector in in SIMC LAB system, with x down, and y to the left. ! Computer Parallel, Perpendicular, and Out of Plane componenants. ! Parallel is component along q_hat. Out/plane is component along ! (e_hat) x (q_hat) (oop_x,oop_y are components of out/plane unit vector) ! Perp. component is what's left: along (q_hat) x (oop_hat). ! So looking along q, out of plane is down, perp. is left. vertex.Pmx = vertex.p.P*vertex.up.x - vertex.q*vertex.uq.x vertex.Pmy = vertex.p.P*vertex.up.y - vertex.q*vertex.uq.y vertex.Pmz = vertex.p.P*vertex.up.z - vertex.q*vertex.uq.z vertex.Pmiss = sqrt(vertex.Pmx**2+vertex.Pmy**2+vertex.Pmz**2) vertex.Emiss = vertex.omega + targ.M - vertex.p.E !STILL NEED SIGN FOR PmPer!!!!!! oop_x = -vertex.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat) oop_y = vertex.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat) vertex.PmPar = (vertex.Pmx*vertex.uq.x + vertex.Pmy*vertex.uq.y + vertex.Pmz*vertex.uq.z) vertex.PmOop = (vertex.Pmx*oop_x + vertex.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2) vertex.PmPer = sqrt( max(0.d0, vertex.Pm**2 - vertex.PmPar**2 - vertex.PmOop**2 ) ) c c if (debug(4)) write(6,*)'comp_ev: at 9',vertex.Pmx,vertex.Pmy,vertex.Pmz ! Calculate Em, Pm, Mrec, Trec for all cases not already done. ! For doing_heavy, get Mrec from nu+M=Ep+Erec, and Erec**2=Mrec**2+Pm**2 ! For (e,e'pi/K), could go back and determine momentum of recoil struck ! particle, and get Mrec and Trec seperately for struck nucleon(hyperon) ! and A-1 system. For now, just take Trec for the A-1 system, and ignore ! the recoiling struck nucleon (hyperon), so Trec=0 for hydrogen target. if (doing_hyd_elast) then vertex.Trec = 0.0 else if (doing_deuterium) then vertex.Pm = vertex.Pmiss vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec else if (doing_heavy) then vertex.Pm = vertex.Pmiss vertex.Mrec = sqrt(vertex.Emiss**2-vertex.Pmiss**2) vertex.Em = targ.Mtar_struck + vertex.Mrec - targ.M vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec else if (doing_hydpi .or. doing_hydkaon .or. doing_dvcs) then vertex.Trec = 0.0 else if (doing_deutpi.or.doing_deutkaon.or.doing_hepi.or.doing_hekaon) then vertex.Trec = sqrt(vertex.Mrec**2 + vertex.Pm**2) - vertex.Mrec endif if (debug(5)) write(6,*) 'vertex.Pm,vertex.Trec,vertex.Em',vertex.Pm,vertex.Trec,vertex.Em if (debug(4)) write(6,*)'comp_ev: at 10' ! calculate krel for deuteron/heavy pion(kaon). Deuteron is straightforward. ! A>2 case is some approximation for 3He (DJG). if (doing_deutpi .or. doing_deutkaon .or. doing_deutdvcs) then if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) write(6,*) 'BAD MM!!!!! Emiss,Pmiss=',vertex.Emiss, vertex.Pmiss MM = sqrt(max(0.d0,vertex.Emiss**2-vertex.Pmiss**2)) krel = sqrt( max(0.d0,MM**2-4.*targ.Mrec_struck**2) ) else if (doing_hepi .or. doing_hekaon .or. doing_hedvcs) then if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) write(6,*) 'BAD MM!!!!! Emiss,Pmiss=',vertex.Emiss, vertex.Pmiss MM = sqrt(max(0.d0,vertex.Emiss**2-vertex.Pmiss**2)) krelx = vertex.Pmx + 1.5*pferx*pfer krely = vertex.Pmy + 1.5*pfery*pfer krelz = vertex.Pmz + 1.5*pferz*pfer krel = sqrt(krelx**2+krely**2+krelz**2) if(vertex.Em.lt.6.0) krel = -krel !bound state test??? endif ntup.krel = krel ! Determine PHYSICS scattering angles theta/phi for the two spectrometer ! vectors, and the Jacobian which we must include in our xsec computation ! to take into account the fact that we're not generating in the physics ! solid angle d(omega) but in 'spectrometer angles' th, ph. call physics_angles(spec.e.theta,spec.e.phi, & vertex.e.xptar,vertex.e.yptar,vertex.e.theta,vertex.e.phi) call physics_angles(spec.p.theta,spec.p.phi, & vertex.p.xptar,vertex.p.yptar,vertex.p.theta,vertex.p.phi) r = sqrt(1.+vertex.e.yptar**2+vertex.e.xptar**2) !r=cos(theta-theta0) main.jacobian = main.jacobian / r**3 !1/cos**3(theta-theta0) if (doing_heavy .or. doing_pion .or. doing_kaon .or. doing_dvcs) then r = sqrt(1.+vertex.p.yptar**2+vertex.p.xptar**2) main.jacobian = main.jacobian / r**3 !1/cos**3(theta-theta0) endif if (debug(4)) write(6,*)'comp_ev: at 11' ! The effective target thickness that the scattered particles see, and the ! resulting energy losses call trip_thru_target (2, main.target.z-targ.zoffset, vertex.e.E, > vertex.e.theta, main.target.Eloss(2), main.target.teff(2),Me,1) call trip_thru_target (3, main.target.z-targ.zoffset, vertex.p.E, > vertex.p.theta, main.target.Eloss(3), main.target.teff(3),Mh,1) if (.not.using_Eloss) then main.target.Eloss(2) = 0.0 main.target.Eloss(3) = 0.0 endif if (debug(4)) write(6,*)'comp_ev: at 12' ! Initialize parameters necessary for radiative corrections call radc_init_ev(main,vertex) ! Success code success = .true. if(debug(2)) write(6,*)'comp_ev: at end...' return end !--------------------------------------------------------------------- subroutine complete_recon_ev(recon,success) implicit none include 'simulate.inc' real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z real*8 targ_new_x,targ_new_y real*8 new_y_x,new_y_y,new_y_z,dummy real*8 px,py,pz,qx,qy,qz real*8 targx,targy,targz,pol real*8 W2 real*8 oop_x,oop_y real*8 mm,mmA,mm2,mmA2,t logical success record /event/ recon !----------------------------------------------------------------------- ! Calculate everything left in the /event/ structure, given ! recon.Ein, and the following for both recon.e AND recon.p: ! delta,xptar,yptar, z,P,E,theta,phi (with E,P corrected for eloss, if desired). ! ! The SINGLE element of /event/ NOT computed here is sigcc (for (e,e'p)). ! !----------------------------------------------------------------------- ! Initialize success = .false. recon.Ein = Ebeam_vertex_ave !lowered by most probable eloss (init.f) ! ... unit vector components of outgoing e,p ! ... z is DOWNSTREAM, x is DOWN and y is LEFT looking downstream. if (debug(4)) write(6,*)'comp_rec_ev: at 1' recon.ue.x = sin(recon.e.theta)*cos(recon.e.phi) recon.ue.y = sin(recon.e.theta)*sin(recon.e.phi) recon.ue.z = cos(recon.e.theta) recon.up.x = sin(recon.p.theta)*cos(recon.p.phi) recon.up.y = sin(recon.p.theta)*sin(recon.p.phi) recon.up.z = cos(recon.p.theta) if (debug(4)) write(6,*)'comp_rec_ev: at 2' ! The q vector if (debug(5)) write(6,*)'comp_rec_ev: Ein,E,uez=',recon.Ein,recon.e.E,recon.ue.z recon.omega = recon.Ein - recon.e.E recon.Q2 = 2*recon.Ein*recon.e.E*(1-recon.ue.z) recon.q = sqrt(recon.Q2 + recon.omega**2) recon.uq.x = - recon.e.P*recon.ue.x / recon.q recon.uq.y = - recon.e.P*recon.ue.y / recon.q recon.uq.z =(recon.Ein - recon.e.P*recon.ue.z)/ recon.q W2 = targ.M**2 + 2.*targ.M*recon.omega - recon.Q2 recon.W = sqrt(abs(W2)) * W2/abs(W2) if (debug(4)) write(6,*)'comp_rec_ev: at 5' ! Now complete the p side if(doing_phsp)then recon.p.P=spec.p.P recon.p.E=sqrt(Mh2+recon.p.P**2) if (debug(4)) write(6,*)'comp_rec_ev: at 7.5',Mh2,recon.p.E endif ! Compute some pion/kaon stuff. recon.epsilon=1./(1. + 2.*(1+recon.omega**2/recon.Q2)*tan(recon.e.theta/2.)**2) recon.theta_pq=acos(recon.up.x*recon.uq.x+recon.up.y*recon.uq.y+recon.up.z*recon.uq.z) ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE. ! Therefore, define a new system with the z axis parallel to q, and ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane ! particles closer to the downstream beamline than q. ! phi=90 is above the horizontal plane when q points to the right, and ! below the horizontal plane when q points to the left. ! Also take into account the different definitions of x, y and z in ! replay and SIMC: ! As seen looking downstream: replay SIMC (old_simc) ! x right down (left) ! y down left (up) ! z all have z pointing downstream ! ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc qx = -recon.uq.y !convert to 'replay' coord. system qy = recon.uq.x qz = recon.uq.z px = -recon.up.y py = recon.up.x pz = recon.up.z dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2)) new_x_x = -qx*qz/dummy new_x_y = -qy*qz/dummy new_x_z = (qx**2 + qy**2)/dummy dummy = sqrt(qx**2 + qy**2) new_y_x = qy/dummy new_y_y = -qx/dummy new_y_z = 0.0 p_new_x = px*new_x_x + py*new_x_y + pz*new_x_z p_new_y = px*new_y_x + py*new_y_y + pz*new_y_z recon.phi_pq = atan2(p_new_y,p_new_x) if(recon.phi_pq.lt.0.) recon.phi_pq = 2.*pi + recon.phi_pq c if ((p_new_x**2+p_new_y**2).eq.0.) then c recon.phi_pq = 0.0 c else c recon.phi_pq = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2)) c endif c if (p_new_y.lt.0.) then c recon.phi_pq = 2*3.1415926536 - recon.phi_pq c endif ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND TARGET POL. ! Therefore, define a new system with the z axis parallel to q, and ! the x axis inside the q-z-plane: z' = q, y' = q X z, x' = y' X q ! this gives phi the way it is usually defined, i.e. phi=0 for in-plane ! particles closer to the downstream beamline than q. ! phi=90 is above the horizontal plane when q points to the right, and ! below the horizontal plane when q points to the left. ! Also take into account the different definitions of x, y and z in ! replay and SIMC: ! As seen looking downstream: replay SIMC (old_simc) ! x right down (left) ! y down left (up) ! z all have z pointing downstream ! ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc qx = -recon.uq.y !convert to 'replay' coord. system qy = recon.uq.x qz = recon.uq.z pol = 1.0 targx = -targ_pol*sin(targ_Bangle) ! 'replay' coordinates targy = 0.0 targz = targ_pol*cos(targ_Bangle) dummy=sqrt((qx**2+qy**2)*(qx**2+qy**2+qz**2)) new_x_x = -qx*qz/dummy new_x_y = -qy*qz/dummy new_x_z = (qx**2 + qy**2)/dummy dummy = sqrt(qx**2 + qy**2) new_y_x = qy/dummy new_y_y = -qx/dummy new_y_z = 0.0 p_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z p_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z recon.phi_targ = atan2(p_new_y,p_new_x) if(recon.phi_targ.lt.0.) recon.phi_targ = 2.*pi + recon.phi_targ c if ((p_new_x**2+p_new_y**2).eq.0.) then c recon.phi_targ = 0.0 c else c recon.phi_targ = acos(p_new_x/sqrt(p_new_x**2+p_new_y**2)) c endif c if (p_new_y.lt.0.) then c recon.phi_targ = 2*3.1415926536 - recon.phi_targ c endif c recon.phi_targ = atan2(p_new_y,p_new_x) ! CALCULATE ANGLE BETA BETWEEN REACTION PLANE AND TRANSVERSE TARGET ! POLARIZATION. ! Therefore, define a new system with the z axis parallel to q, and ! with the y-direction defined by q x p_meson and x-axis given ! by y x z ! Also take into account the different definitions of x, y and z in ! replay and SIMC: ! As seen looking downstream: replay SIMC (old_simc) ! x right down (left) ! y down left (up) ! z all have z pointing downstream ! ! SO: x_replay=-y_simc, y_replay=x_simc, z_replay= z_simc qx = -recon.uq.y !convert to 'replay' coord. system qy = recon.uq.x qz = recon.uq.z pol = 1.0 targx = -targ_pol*sin(targ_Bangle) ! 'replay' coordinates targy = 0.0 targz = targ_pol*cos(targ_Bangle) px = -recon.up.y py = recon.up.x pz = recon.up.z dummy = sqrt((qy*pz-qz*py)**2 + (qz*px-qx*pz)**2 + (qx*py-qy*px)**2) new_y_x = (qy*pz-qz*py)/dummy new_y_y = (qz*px-qx*pz)/dummy new_y_z = (qx*py-qy*px)/dummy dummy = sqrt((new_y_y*qz-new_y_z*qy)**2 + (new_y_z*qx-new_y_x*qz)**2 1 + (new_y_x*qy-new_y_y*qx)**2) new_x_x = (new_y_y*qz - new_y_z*qy)/dummy new_x_y = (new_y_z*qx - new_y_x*qz)/dummy new_x_z = (new_y_x*qy - new_y_y*qx)/dummy targ_new_x = targx*new_x_x + targy*new_x_y + targz*new_x_z targ_new_y = targx*new_y_x + targy*new_y_y + targz*new_y_z recon.beta = atan2(targ_new_y,targ_new_x) if(recon.beta.lt.0.) recon.beta = 2.*pi + recon.beta c if ((targ_new_x**2+targ_new_y**2).eq.0.) then c recon.beta = 0.0 c else c recon.beta = acos(targ_new_x/sqrt(targ_new_x**2+targ_new_y**2)) c endif c if (targ_new_y.lt.0.) then c recon.beta = 2*3.1415926536 - recon.beta c endif dummy = sqrt((qx**2+qy**2+qz**2))*sqrt((targx**2+targy**2+targz**2)) recon.theta_tarq = acos((qx*targx+qy*targy+qz*targz)/dummy) c write(6,*) 'cheesy poofs', recon.beta, p_new_x, p_new_y,qx,targy c if (debug(2)) then c write(6,*)'comp_rec_ev: omega =',recon.omega c write(6,*)'comp_rec_ev: Q2 =',recon.Q2 c write(6,*)'comp_rec_ev: theta_e =',recon.e.theta c write(6,*)'comp_rec_ev: epsilon =',recon.epsilon c write(6,*)'comp_rec_ev: theta_pq =',recon.theta_pq c write(6,*)'comp_rec_ev: phi_pq =',recon.phi_pq c endif if (debug(4)) write(6,*)'comp_rec_ev: at 8' ! Compute the Pm vector in in SIMC LAB system, with x down, and y to the left. ! Computer Parallel, Perpendicular, and Out of Plane componenants. ! Parallel is component along q_hat. Out/plane is component along ! (e_hat) x (q_hat) (oop_x,oop_y are components of out/plane unit vector) ! Perp. component is what's left: along (q_hat) x (oop_hat). ! So looking along q, out of plane is down, perp. is left. recon.Pmx = recon.p.P*recon.up.x - recon.q*recon.uq.x recon.Pmy = recon.p.P*recon.up.y - recon.q*recon.uq.y recon.Pmz = recon.p.P*recon.up.z - recon.q*recon.uq.z recon.Pm = sqrt(recon.Pmx**2+recon.Pmy**2+recon.Pmz**2) !STILL NEED SIGN FOR PmPer!!!!!! oop_x = -recon.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat) oop_y = recon.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat) recon.PmPar = (recon.Pmx*recon.uq.x + recon.Pmy*recon.uq.y + recon.Pmz*recon.uq.z) recon.PmOop = (recon.Pmx*oop_x + recon.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2) recon.PmPer = sqrt( max(0.d0, recon.Pm**2 - recon.PmPar**2 - recon.PmOop**2 ) ) if (doing_pion .or. doing_kaon .or. doing_dvcs) then recon.Em = recon.omega + targ.Mtar_struck - recon.p.E mm2 = recon.Em**2 - recon.Pm**2 mm = sqrt(abs(mm2)) * abs(mm2)/mm2 mmA2= (recon.omega + targ.M - recon.p.E)**2 - recon.Pm**2 mmA = sqrt(abs(mmA2)) * abs(mmA2)/mmA2 t = recon.Q2 - Mh2 > + 2*(recon.omega*recon.p.E - recon.p.P*recon.q*cos(recon.theta_pq)) ntup.mm = mm ntup.mmA = mmA ntup.t = t endif ! Calculate Trec, Em. Trec for (A-1) system (eep), or for struck nucleon (pi/K) ! Note that there are other ways to calculate 'Em' for the pion/kaon case. if (doing_hyd_elast) then recon.Trec = 0 recon.Em = recon.Ein + Mp - recon.e.E - recon.p.E - recon.Trec else if (doing_deuterium .or. doing_heavy) then recon.Trec = sqrt(recon.Pm**2+targ.Mrec**2) - targ.Mrec recon.Em = recon.Ein + Mp - recon.e.E - recon.p.E - recon.Trec else if (doing_pion) then recon.Em = recon.omega + targ.Mtar_struck - recon.p.E else if (doing_kaon) then recon.Em = recon.omega + targ.Mtar_struck - recon.p.E else if (doing_dvcs) then recon.Em = recon.omega + targ.Mtar_struck - recon.p.E endif if (debug(5)) write(6,*) 'recon.Pm,recon.Trec,recon.Em',recon.Pm,recon.Trec,recon.Em if (debug(4)) write(6,*)'comp_rec_ev: at 10' success=.true. return end !------------------------------------------------------------------------ subroutine complete_main(force_sigcc,main,vertex,recon,success) implicit none include 'simulate.inc' integer i, iPm1 real*8 a, b, r, frac, peepi, peeK real*8 survivalprob real*8 weight, width, sigep, deForest, tgtweight logical force_sigcc, success record /event_main/ main record /event/ vertex, recon !----------------------------------------------------------------------- ! Calculate everything left in the /main/ structure that hasn't been ! required up til now. This routine is called ONLY after we ! know an event IS going to contribute, don't want to be doing these ! calculations if they're not needed! ! ! A small anomaly is that sigcc from the /event/ structure is computed ! here since it's not needed til now, and was not calculated by ! COMPLETE_ev. !----------------------------------------------------------------------- ! Initialize success code if (debug(2)) write(6,*)'comp_main: entering...' success = .false. ! The spectral function weighting if (doing_hyd_elast.or.doing_pion.or.doing_kaon.or.doing_dvcs.or.doing_phsp) then !no SF. main.SF_weight=1.0 else if (doing_deuterium .or. doing_heavy) then main.SF_weight = 0.0 do i=1,nrhoPm weight = 0.0 ! ... use linear interpolation to determine rho(pm) r = (vertex.Pm-Pm_theory(i).min)/Pm_theory(i).bin if (r.ge.0 .and. r.le.Pm_theory(i).n) then iPm1 = nint(r) if (iPm1.eq.0) iPm1 = 1 if (iPm1.eq.Pm_theory(i).n) iPm1 = Pm_theory(i).n-1 frac = r+0.5-float(iPm1) b = theory(i,iPm1) a = theory(i,iPm1+1)-b weight = a*frac + b endif if (doing_heavy) then width=Emsig_theory(i)/2.0 if (vertex.Em.lt.E_Fermi) weight = 0.0 weight=weight/pi/Em_int_theory(i) * width/ > ((vertex.Em-Em_theory(i))**2+width**2) endif main.SF_weight = main.SF_weight+weight*nprot_theory(i) enddo ! endif ! ... if we have come up with weight<=, quit now (avoid weight<0 in ntuple) ! ... unless FORCE_SIGCC is on (to get central sigcc). if (debug(5))write(6,*)'comp_main: calculating cross section' if (main.SF_weight.le.0 .and. .not.force_sigcc) return ! ... Cross section, for BOTH vertex and recon (where possible) ! ... Note that all of these are cross section per nucleon. For A(e,e'p) ! ... this becomes cross section per nucleus because of the weighting of ! ... the spectral function. For pion/kaon production, we explicitly ! ... apply a weight for the number of nucleons involved (protons for pi+ or Y0, ! ... neutrons for pi- or Y-) (Y=hyperon) tgtweight = 1.0 if (doing_phsp) then main.sigcc = 1.0 main.sigcc_recon = 1.0 elseif (doing_hyd_elast) then main.sigcc = sigep(vertex) main.sigcc_recon = sigep(recon) elseif (doing_deuterium.or.doing_heavy) then main.sigcc = deForest(vertex) main.sigcc_recon = deForest(recon) elseif (doing_pion) then main.sigcc = peepi(vertex,main) main.sigcc_recon = 1.0 if (which_pion.eq.1 .or. which_pion.eq.11) then !OK for coherent??? tgtweight = targ.N else tgtweight = targ.Z endif elseif (doing_kaon) then main.sigcc = peeK(vertex,main,survivalprob) main.sigcc_recon = 1.0 if (which_kaon.eq.2 .or. which_kaon.eq.12) then !OK for coherent??? tgtweight = targ.N else tgtweight = targ.Z endif elseif (doing_dvcs) then main.sigcc = 1.0 main.sigcc_recon = 1.0 if (which_dvcs.eq.0) then !OK for coherent??? tgtweight = targ.Z else if(which_dvcs.eq.10) then tgtweight = targ.Z else tgtweight = targ.N endif else main.sigcc = 1.0 main.sigcc_recon = 1.0 endif if (debug(3)) then write(6,*)'======================================' write(6,*)'complete_main:' write(6,*)'theta_e =',vertex.e.theta write(6,*)'q mag =',vertex.q write(6,*)'omega =',vertex.omega write(6,*)'Q2 =',vertex.Q2/1000000. write(6,*)'Ein =',vertex.Ein write(6,*)'hadron mom =',vertex.p.P write(6,*)'e mom =',vertex.e.P write(6,*)'mass =',Mh write(6,*)'epsilon =',main.epsilon write(6,*)'phi_pq =',main.phi_pq write(6,*)'theta_pq =',main.theta_pq write(6,*)'======================================' endif ! The total contributing weight from this event -- this weight is ! proportional to # experimental counts represented by the event. ! Apply survival probability to kaons if we're not modeling decay. main.weight = main.SF_weight*main.jacobian*main.gen_weight*main.sigcc main.weight = main.weight * tgtweight !correct for #/nucleons involved if (doing_kaon .and. .not.doing_decay) > main.weight = main.weight*survivalprob if (debug(5))write(6,*) 'gen_weight = ',main.gen_weight, > main.jacobian,main.sigcc success = .true. if (debug(2)) write(6,*)'comp_main: ending, success =',success return end subroutine physics_angles(theta0,phi0,dx,dy,theta,phi) !Generate physics angles in lab frame. Theta is angle from beamline. !phi is projection on x-y plane (so phi=0 is down, sos=pi/2, hms=3*pi/2. ! !theta=acos( (cos(theta0)-dy*sin(theta0)*sin(phi0))/ sqrt(1+dx**2+dy**2) ) !phi=atan( (dy*cos(theta0)+sin(theta0)*sin(phi0)) / (sin(theta0)*cos(phi0)+dx) ) ! ! Note that these formulae assume phi0=pi/2 or 3*pi/2 (thus the sin(phi0) ! gives a -/+ sign for the HMS/SOS). Thus, we set the cos(phi0) term to zero. real*8 dx,dy !dx/dy (xptar/yptar) for event. real*8 theta0,phi0 !central physics angles of spectrometer. real*8 theta,phi !physics angles for event. real*8 r,sinth,costh,sinph,cosph !intermediate variables. real*8 tmp include 'simulate.inc' costh = cos(theta0) sinth = sin(theta0) sinph = sin(phi0) cosph = cos(phi0) r = sqrt( 1. + dx**2 + dy**2 ) if (abs(cosph).gt.0.0001) then !phi not at +/- pi/2 write(6,*) 'theta,phi will be incorrect if phi0 <> pi/2 or 3*pi/2' write(6,*) 'phi0=',phi0,'=',phi0*180/pi,'degrees' endif tmp=(costh - dy*sinth*sinph) / r if (abs(tmp).gt.1) write(6,*) 'tmp=',tmp theta = acos( (costh - dy*sinth*sinph) / r ) if (dx.ne.0.0) then phi = atan( (dy*costh + sinth*sinph) / dx ) !gives -90 to 90 deg. if (phi.le.0) phi=phi+pi !make 0 to 180 deg. if (sinph.lt.0.) phi=phi+pi !add pi to phi for HMS else phi = phi0 endif return end