Return to event_orig.f CVS log

Up to [HallC] / Poltar

 Import simc poltar

	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 ! <i=1,nrhoPm>
	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