|
version 1.3, 2005/01/26 22:11:51
|
version 1.10, 2007/08/31 17:39:02
|
|
|
|
| ! gen.xwid and gen.ywid are the intrinsic beam widths (one sigma value). | ! 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). | ! Use a gaussian beam distribution with +/- 3.0 sigma (add raster afterwards). |
| | |
| |
C DJG As best I can figger, main.target.y is positive when the beam is high in the lab |
| |
C DJG and main.target.x is positive when the beam is right when looking downstream. |
| |
C DJG Don't ask me why, but it seems to be this way. |
| |
C DJG Note that this means that +fry points down. I will make frx point left. |
| |
|
| main.target.x = gauss1(nsig_max)*gen.xwid+targ.xoffset | main.target.x = gauss1(nsig_max)*gen.xwid+targ.xoffset |
| main.target.y = gauss1(nsig_max)*gen.ywid+targ.yoffset | 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=1 - old bedpost raster rectangle - targ.fr1/fr2 are the x/y raster half-widths. |
| ! fr_pattern=2 - circle - targ.fr1/fr2 are the inner and outer radii. | ! fr_pattern=2 - circle - targ.fr1/fr2 are the inner and outer radii. |
| |
! fr_pattern=3 - new flat raster rectangle - targ.fr1/fr2 are the x/y raster half-widths. |
| | |
| if (targ.fr_pattern .eq. 1) then !square raster |
if (targ.fr_pattern .eq. 1) then !old bedpost, square raster |
| t3=grnd()*pi | t3=grnd()*pi |
| t4=grnd()*pi | t4=grnd()*pi |
| t5=cos(t3)*targ.fr1 | t5=cos(t3)*targ.fr1 |
|
|
|
| t4=sqrt(grnd())*(targ.fr2-targ.fr1)+targ.fr1 | t4=sqrt(grnd())*(targ.fr2-targ.fr1)+targ.fr1 |
| t5=cos(t3)*t4 | t5=cos(t3)*t4 |
| t6=sin(t3)*t4 | t6=sin(t3)*t4 |
| |
elseif (targ.fr_pattern .eq. 3) then !new, flat square raster |
| |
t3=2.*grnd()-1.0 |
| |
t4=2.*grnd()-1.0 |
| |
t5=targ.fr1*t3 |
| |
t6=targ.fr2*t4 |
| else !no raster | else !no raster |
| t5=0.0 | t5=0.0 |
| t6=0.0 | t6=0.0 |
|
|
|
| main.target.y = main.target.y+t6 | main.target.y = main.target.y+t6 |
| main.target.z = (0.5-grnd())*targ.length+targ.zoffset | main.target.z = (0.5-grnd())*targ.length+targ.zoffset |
| main.target.rastery = t6 !'raster' contribution to vert. pos. | main.target.rastery = t6 !'raster' contribution to vert. pos. |
| |
main.target.rasterx = t5 ! points right as you look downstream - need to flip sign later. |
| | |
| ! Take fluctuations of the beam energy into account, and remember to correct | ! Take fluctuations of the beam energy into account, and remember to correct |
| ! for ionization loss in the target and Coulomb acceleration of incoming | ! for ionization loss in the target and Coulomb acceleration of incoming |
|
|
|
| 1 .or.doing_deutrho .or.doing_deutsemi) then !Em = binding energy | 1 .or.doing_deutrho .or.doing_deutsemi) then !Em = binding energy |
| vertex.Em = Mp + Mn - targ.M | vertex.Em = Mp + Mn - targ.M |
| m_spec = targ.M - targ.Mtar_struck + vertex.Em != Mn(Mp) for pi+(-) | m_spec = targ.M - targ.Mtar_struck + vertex.Em != Mn(Mp) for pi+(-) |
| c write(6,*) 'in event.f',targ.M,m_spec,pfer |
|
| efer = targ.M - sqrt(m_spec**2+pfer**2) | efer = targ.M - sqrt(m_spec**2+pfer**2) |
| endif | endif |
| if (doing_hepi .or. doing_hekaon .or. doing_hedelta .or. doing_herho) then | if (doing_hepi .or. doing_hekaon .or. doing_hedelta .or. doing_herho) then |
|
|
|
| | |
| real*8 a, b, c, r, t, QA, QB, QC, radical | 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 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 new_y_x,new_y_y,new_y_z,dummy |
| |
real*8 targx,targy,targz |
| real*8 px,py,pz,qx,qy,qz | real*8 px,py,pz,qx,qy,qz |
| real*8 oop_x,oop_y | real*8 oop_x,oop_y |
| real*8 krel,krelx,krely,krelz | real*8 krel,krelx,krely,krelz |
| real*8 MM | real*8 MM |
| real*8 dot |
|
| | |
| real*8 Ehad2,E_rec | real*8 Ehad2,E_rec |
| real*8 W2 | real*8 W2 |
|
|
|
| | |
| radical = QB**2 - 4.*QA*QC | radical = QB**2 - 4.*QA*QC |
| if (radical.lt.0) return | if (radical.lt.0) return |
| |
|
| vertex.p.E = (-QB - sqrt(radical))/2./QA | vertex.p.E = (-QB - sqrt(radical))/2./QA |
| Ehad2 = (-QB + sqrt(radical))/2./QA | Ehad2 = (-QB + sqrt(radical))/2./QA |
| | |
|
|
|
| | |
| | |
| if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then | if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then |
| W2 = targ.M**2 + 2.*targ.M*vertex.nu - vertex.Q2 |
W2 = targ.Mtar_struck**2 + 2.*targ.Mtar_struck*vertex.nu - vertex.Q2 |
| main.W = sqrt(abs(W2)) * W2/abs(W2) | main.W = sqrt(abs(W2)) * W2/abs(W2) |
| main.epsilon=1./(1. + 2.*(1+vertex.nu**2/vertex.Q2)*tan(vertex.e.theta/2.)**2) | main.epsilon=1./(1. + 2.*(1+vertex.nu**2/vertex.Q2)*tan(vertex.e.theta/2.)**2) |
| dot = vertex.up.x*vertex.uq.x+vertex.up.y*vertex.uq.y+vertex.up.z*vertex.uq.z |
main.theta_pq=acos(vertex.up.x*vertex.uq.x+vertex.up.y*vertex.uq.y+vertex.up.z*vertex.uq.z) |
| main.theta_pq=acos(dot) |
|
| main.t = vertex.Q2 - Mh2 + 2*vertex.nu*vertex.p.E - | main.t = vertex.Q2 - Mh2 + 2*vertex.nu*vertex.p.E - |
| > 2*vertex.p.P*vertex.q*cos(main.theta_pq) | > 2*vertex.p.P*vertex.q*cos(main.theta_pq) |
| main.tmin = vertex.Q2 - Mh2 + 2*vertex.p.E*vertex.nu - | main.tmin = vertex.Q2 - Mh2 + 2*vertex.p.E*vertex.nu - |
|
|
|
| write(6,*)'comp_ev: phi_pq =',main.phi_pq | write(6,*)'comp_ev: phi_pq =',main.phi_pq |
| write(6,*)'comp_ev: E_hadron =',vertex.p.E | write(6,*)'comp_ev: E_hadron =',vertex.p.E |
| endif | endif |
| |
|
| |
|
| |
if(using_tgt_field) then !calculate some azimuthal angles that only make |
| |
!sense for polarized target |
| |
! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND TARGET POLARIZATION. |
| |
! 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 |
| |
|
| |
c Target in-plane, so targy=0 |
| |
targx = -targ_pol*sin(abs(targ_Bangle)) ! replay coordinates |
| |
targy = 0.0 |
| |
targz = targ_pol*cos(abs(targ_Bangle)) |
| |
|
| |
c Target out of plane, so targx=0 |
| |
c targx = 0.0 ! replay coordinates |
| |
c targy = -targ_pol*sin(abs(targ_Bangle)) |
| |
c targz = targ_pol*cos(abs(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. |
| |
! 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 |
| |
|
| |
C Taret in plane |
| |
targx = -targ_pol*sin(abs(targ_Bangle)) ! 'replay' coordinates |
| |
targy = 0.0 |
| |
targz = targ_pol*cos(abs(targ_Bangle)) |
| |
|
| |
C Target out of plane |
| |
|
| |
c targx = 0.0 ! 'replay' coordinates |
| |
c targy = -targ_pol*sin(abs(targ_Bangle)) |
| |
c targz = targ_pol*cos(abs(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 |
| |
|
| |
|
| |
CDJG Calculate the "Collins" (phi_pq+phi_targ) and "Sivers"(phi_pq-phi_targ) angles |
| |
vertex.phi_s = main.phi_pq-main.phi_targ |
| |
if(vertex.phi_s .lt. 0.) vertex.phi_s = 2*pi+vertex.phi_s |
| |
|
| |
vertex.phi_c = main.phi_pq+main.phi_targ |
| |
if(vertex.phi_c .gt. 2.*pi) vertex.phi_c = vertex.phi_c-2*pi |
| |
if(vertex.phi_c .lt. 0.0) vertex.phi_c = 2*pi+vertex.phi_c |
| |
|
| |
|
| |
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) |
| |
|
| |
endif !polarized-target specific azimuthal angles |
| |
|
| |
|
| endif !end of pion/kaon specific stuff. | endif !end of pion/kaon specific stuff. |
| | |
| if (debug(4)) write(6,*)'comp_ev: at 8' | if (debug(4)) write(6,*)'comp_ev: at 8' |
|
|
|
| CDJG if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) then | CDJG if ((vertex.Emiss**2-vertex.Pmiss**2).lt.0) then |
| CDJG I should be testing that the missing mass is above two pion | CDJG I should be testing that the missing mass is above two pion |
| CDJG threshold! Otherwise, it's just exclusive | CDJG threshold! Otherwise, it's just exclusive |
| if((vertex.Emiss**2-vertex.Pmiss**2).lt.(targ.M+Mpi0)**2)then |
c if ((vertex.Emiss**2-vertex.Pmiss**2).lt.(Mp+Mpi0)**2) then |
| c if((vertex.Emiss**2-vertex.Pmiss**2).lt.(targ.Mtar_struck+Mpi0)**2)then |
if (((targ.Mtar_struck+vertex.nu-vertex.p.E)**2-vertex.Pmiss**2).lt.(Mp+Mpi0)**2) then |
| success=.false. | success=.false. |
| return | return |
| endif | endif |
|
|
|
| | |
| if(doing_semi) then | if(doing_semi) then |
| vertex.zhad = vertex.p.E/vertex.nu | vertex.zhad = vertex.p.E/vertex.nu |
| vertex.pt2 = vertex.p.P**2*(1.0-dot**2) |
vertex.pt2 = vertex.p.P**2*(1.0-cos(main.theta_pq)**2) |
| if(vertex.zhad.gt.1.0) then | if(vertex.zhad.gt.1.0) then |
| success=.false. | success=.false. |
| return | return |
|
|
|
| | |
| real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z | real*8 p_new_x,p_new_y,new_x_x,new_x_y,new_x_z |
| real*8 new_y_x,new_y_y,new_y_z,dummy | real*8 new_y_x,new_y_y,new_y_z,dummy |
| |
real*8 targ_new_x,targ_new_y |
| real*8 px,py,pz,qx,qy,qz | real*8 px,py,pz,qx,qy,qz |
| |
real*8 targx,targy,targz |
| real*8 W2 | real*8 W2 |
| real*8 oop_x,oop_y | real*8 oop_x,oop_y |
| real*8 dot |
|
| real*8 mm,mmA,mm2,mmA2,t | real*8 mm,mmA,mm2,mmA2,t |
| | |
| logical success | logical success |
|
|
|
| recon.uq.x = - recon.e.P*recon.ue.x / recon.q | recon.uq.x = - recon.e.P*recon.ue.x / recon.q |
| recon.uq.y = - recon.e.P*recon.ue.y / 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 | recon.uq.z =(recon.Ein - recon.e.P*recon.ue.z)/ recon.q |
| |
if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then |
| |
W2 = targ.mtar_struck**2 + 2.*targ.mtar_struck*recon.nu - recon.Q2 |
| |
else |
| W2 = targ.M**2 + 2.*targ.M*recon.nu - recon.Q2 | W2 = targ.M**2 + 2.*targ.M*recon.nu - recon.Q2 |
| |
endif |
| |
|
| |
W2 = mp**2 + 2.*mp*recon.nu - recon.Q2 |
| recon.W = sqrt(abs(W2)) * W2/abs(W2) | recon.W = sqrt(abs(W2)) * W2/abs(W2) |
| recon.xbj = recon.Q2/2./Mp/recon.nu | recon.xbj = recon.Q2/2./Mp/recon.nu |
| if (debug(4)) write(6,*)'comp_rec_ev: at 5' | if (debug(4)) write(6,*)'comp_rec_ev: at 5' |
|
|
|
| ! Compute some pion/kaon stuff. | ! Compute some pion/kaon stuff. |
| | |
| recon.epsilon=1./(1. + 2.*(1+recon.nu**2/recon.Q2)*tan(recon.e.theta/2.)**2) | recon.epsilon=1./(1. + 2.*(1+recon.nu**2/recon.Q2)*tan(recon.e.theta/2.)**2) |
| dot = recon.up.x*recon.uq.x+recon.up.y*recon.uq.y+recon.up.z*recon.uq.z |
recon.theta_pq=acos(min(1.0,recon.up.x*recon.uq.x+recon.up.y*recon.uq.y+recon.up.z*recon.uq.z)) |
| recon.theta_pq=acos(min(1.0,dot)) |
|
| | |
| ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE. | ! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND REACTION PLANE. |
| ! Therefore, define a new system with the z axis parallel to q, and | ! Therefore, define a new system with the z axis parallel to q, and |
|
|
|
| recon.phi_pq = 2*pi - recon.phi_pq | recon.phi_pq = 2*pi - recon.phi_pq |
| endif | endif |
| | |
| |
|
| |
if(using_tgt_field) then !calculate polarized-target specific azimuthal angles |
| |
|
| |
! CALCULATE ANGLE PHI BETWEEN SCATTERING PLANE AND TARGET POL. |
| |
! 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 |
| |
|
| |
C In-plane target |
| |
targx = -targ_pol*sin(abs(targ_Bangle)) ! 'replay' coordinates |
| |
targy = 0.0 |
| |
targz = targ_pol*cos(abs(targ_Bangle)) |
| |
|
| |
C out of plane target |
| |
c targx = 0.0 ! 'replay' coordinates |
| |
c targy = -targ_pol*sin(abs(targ_Bangle)) |
| |
c targz = targ_pol*cos(abs(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 |
| |
|
| |
|
| |
! CALCULATE ANGLE BETA BETWEEN REACTION PLANE AND (TRANSVERSE) TARGET |
| |
! POLARIZATION. |
| |
! 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 |
| |
|
| |
targx = -targ_pol*sin(abs(targ_Bangle)) ! 'replay' coordinates |
| |
targy = 0.0 |
| |
targz = targ_pol*cos(abs(targ_Bangle)) |
| |
|
| |
c targx = 0.0 ! 'replay' coordinates |
| |
c targy = -targ_pol*sin(abs(targ_Bangle)) |
| |
c targz = targ_pol*cos(abs(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 |
| |
|
| |
CDJG Calculate the "Collins" (phi_pq+phi_targ) and "Sivers"(phi_pq-phi_targ) angles |
| |
recon.phi_s = recon.phi_pq-recon.phi_targ |
| |
if(recon.phi_s .lt. 0.) recon.phi_s = 2*pi+recon.phi_s |
| |
|
| |
recon.phi_c = recon.phi_pq+recon.phi_targ |
| |
if(recon.phi_c .gt. 2.*pi) recon.phi_c = recon.phi_c-2*pi |
| |
if(recon.phi_c .lt. 0.) recon.phi_c = 2*pi+recon.phi_c |
| |
|
| |
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) |
| |
endif !polarized-target specific stuff |
| |
|
| if (debug(2)) then | if (debug(2)) then |
| write(6,*)'comp_rec_ev: nu =',recon.nu | write(6,*)'comp_rec_ev: nu =',recon.nu |
| write(6,*)'comp_rec_ev: Q2 =',recon.Q2 | write(6,*)'comp_rec_ev: Q2 =',recon.Q2 |
|
|
|
| | |
| oop_x = -recon.uq.y ! oop_x = q_y *(z_hat x y_hat) = q_y *-(x_hat) | 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) | oop_y = recon.uq.x ! oop_y = q_x *(z_hat x x_hat) = q_x * (y_hat) |
| c recon.PmPar = (recon.Pmx*recon.uq.x + recon.Pmy*recon.uq.y + recon.Pmz*recon.uq.z) |
recon.PmPar = (recon.Pmx*recon.uq.x + recon.Pmy*recon.uq.y + recon.Pmz*recon.uq.z) |
| recon.PmPar = (recon.Pmy*recon.uq.y + recon.Pmz*recon.uq.z)/sqrt(recon.uq.y**2+recon.uq.z**2) |
|
| recon.PmOop = (recon.Pmx*oop_x + recon.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2) | recon.PmOop = (recon.Pmx*oop_x + recon.Pmy*oop_y) / sqrt(oop_x**2+oop_y**2) |
| c recon.PmPer = sqrt( max(0.d0, recon.Pm**2 - recon.PmPar**2 - recon.PmOop**2 ) ) |
recon.PmPer = sqrt( max(0.d0, recon.Pm**2 - recon.PmPar**2 - recon.PmOop**2 ) ) |
| recon.PmPer = (-recon.Pmz*recon.uq.y + recon.Pmy*recon.uq.z)/sqrt(recon.uq.z**2+recon.uq.y**2) |
|
| | |
| if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then | if (doing_pion .or. doing_kaon .or. doing_delta .or. doing_rho .or. doing_semi) then |
| recon.Em = recon.nu + targ.Mtar_struck - recon.p.E | recon.Em = recon.nu + targ.Mtar_struck - recon.p.E |
|
|
|
| | |
| if(doing_semi.or.doing_rho) then | if(doing_semi.or.doing_rho) then |
| recon.zhad = recon.p.E/recon.nu | recon.zhad = recon.p.E/recon.nu |
| recon.pt2 = recon.p.P**2*(1.0-dot**2) |
recon.pt2 = recon.p.P**2*(1.0-cos(recon.theta_pq)**2) |
| endif | endif |
| | |
| ! Calculate Trec, Em. Trec for (A-1) system (eep), or for struck nucleon (pi/K) | ! Calculate Trec, Em. Trec for (A-1) system (eep), or for struck nucleon (pi/K) |
|
|
|
| | |
| integer i, iPm1 | integer i, iPm1 |
| real*8 a, b, r, frac, peepi, peeK, peedelta, peerho, peepiX | real*8 a, b, r, frac, peepi, peeK, peedelta, peerho, peepiX |
| real*8 survivalprob |
real*8 survivalprob, semi_dilution |
| real*8 weight, width, sigep, deForest, tgtweight | real*8 weight, width, sigep, deForest, tgtweight |
| logical force_sigcc, success | logical force_sigcc, success |
| record /event_main/ main | record /event_main/ main |
|
|
|
| main.sigcc = peepiX(vertex,vertex0,main,survivalprob,.FALSE.) | main.sigcc = peepiX(vertex,vertex0,main,survivalprob,.FALSE.) |
| main.sigcc_recon = 1.0 | main.sigcc_recon = 1.0 |
| main.sigcent = peepiX(vertex,vertex0,main,survivalprob,.TRUE.) | main.sigcent = peepiX(vertex,vertex0,main,survivalprob,.TRUE.) |
| |
ntup.dilu = semi_dilution(vertex,main) |
| | |
| else | else |
| main.sigcc = 1.0 | main.sigcc = 1.0 |
Return to
event.f
CVS log
Up to [HallC] / simc_semi