Notes on changes to pol_hms_single for radiative correction and
how elastic events are treated
----------------------------------------------------------------

A. In radc.f subroutine generate_rad 

1.  Add argument "elastic_event" to generate_rad argument list.
    When photon is radiated from beam electron then if 'elastic event"
    the outgoing electron energy is recalculated using the
    vertex beam energy = initial beam energy - radiated photon energy
    
2.  Add include 'constants.inc' to have mass of proton.

3.  Add variables ebeam_max,ebeam_min . 

4.  Modify calculation of Egamma_min, Egamma_max when ntail=1 .
    In mc_hms_single.f modified the E_max and E_min that are sent
    to be E_min = p_spec*(1.-0.305)and E_max = p_spec*(1.+0.305)
    instead of E_min = p_spec*(1.+0.01*dpp_lo) and
    E_max = p_spec*(1.+0.01*dpp_hi). 
    Before E_gamma_min(1) = 0 now = Ein - E_max
           
5.  When ntail=1 and elastic ep event then calculate ebeam_max and
    ebeam_min and Egamma_min(1) = Ein - ebeam_max and
    egamma_max(1)  = Ein - ebeam_min 

6.  When ntail = 1 after Ein = Ein - Egamma_used(1) statement then
    if elastic event recalculate Ep 

7.  if ntail = 1 Fixed mistake in call the peaked_rad_weight which left out
    the argument Ein needed in the subroutine.

8.  If ntail = 2 change Egamma_min(2) = 0.0 to Egamma_min(2) = Ep - E_max

9.   if ntail = 2 Fixed mistake in call the peaked_rad_weight which left out
    the argument Ein needed in the subroutine.
  

B. in radc.f function peaked_rad_wieght

1. Add Ein to arguments of the function. It was needed by function and
   previously Ein was set to 0 ( since it wasn't specifically set)
   and led to variables calculated as NaN sinc 1/Ein was used in a formula.

C. In mc_hms_single.f

1. Decide to treat elastic ep and inelastic ep events separately.
   So if scattering on proton then added the following code

	  proton_wgt = 1.
	  elastic_event=0.
          if ( a(imat) .eq. 1 ) then
          if ( grnd() .le. 0.5 ) then 
	    dpp  = grnd()*(dpp_hi - dpp_lo) + dpp_lo
	    proton_wgt = 2.
	  else
            theta_pol = acos( (cos_ts + dydz*sin_ts)
     +                        / sqrt( 1. + dxdz**2 + dydz**2 ) )
	    escat = 938.27*e/(e*(1-cos(theta_pol))+938.27)
	    dpp = (escat-p_spec)/p_spec*100.
	    proton_wgt = 2.  
	    elastic_event=1.0
          endif
          endif

   So random number is thrown if scattering on proton ( a(imat) .eq. 1).
   If the random number larger than 0.5 treat event as an inelastic
   scattering and randomly pick the scattered electron momentum. If 
   random number less than 0.5 then treat event as elastic ep scattering
   and calculate the scattered electron momentum. Now add the variable
   proton_weight which is part of the overall weight for the event.
   Also add variable elastic_event which flags the event for subroutines
   qfs and generate_rad .

2. The radiative correction codes want the unit vector of the scattered electron
   in the beam system with +ue(1) DOWN, +ue(2) beam left and +ue(3) along the beam.
   Previously the code calculated the unit vector in the HMS optics system.

3. Previously divided  the arguements of subroutine radc_init_ev and generate_rad
   that dealt with momentum by 1000 to convert to GeV. But the subroutines really 
   want momentum in MeV.

4. Previously the ep_max and ep_min ( the maximum and minimum scattered electron
   energies) sent to generate_rad where set to 
	ep_min = p_spec*(1.+0.01*dpp_lo)
	ep_max = p_spec*(1.+0.01*dpp_hi)

   To match SIMC decided to have  separate variables and 
   increase range to +- 30% delta 
	ep_rad_min = p_spec*(1.-0.305)
	ep_rad_max = p_spec*(1.+0.305)

5. Add variable "ebeam" . Set it equal to "e" . This is just because I found
   it difficult to track the variable "e" throughout the program. ebeam is
   the measured beam energy in the lab.

6. The radiative correction code determines that the beam electron radiated
   or that the scattered electron radiated. Both can't happen in one event.
   It retursn the radiated photon energy as erad ( photon with beam electron)
   or eprad ( photon with scattered electron) . 
   Calculate the ebeam_vertex = ebeam - erad  and eprime_vertex = ep
   If the scattered electron radiates then ep = ep - erad and delta 
   is also modified dpps = dpps - eprad/p_spec*100. 
   If the event is an elastic ep event and the beam electron radiates
   then the scatter electron energy is recalculated.

7. The argument to qfs subroutine explicitly use the vertex beam energy and 
   transferred energy, enu_vertex = ebeam_vertex - eprime_vertex
   If addition the flag "elastic_event" is sent so that qfs can treat
   elastic and inelastic ep events differently.

8. For proton target and inelastic event calculated W_vertex . If 
    W_vertex < 1.073 then skip event.

D. In qfs_new12_sub.f  

1. Set so that standard for A=1 target the variable dorg is set to 'y'.
   This means that the ep elastic cross section is calculated. If dorg='n'
   is used then the cross section is smeared over W by a gaussian.

2. Pass the flag 'elflag'. When elflag=1 then ep elastic event.

3. Modify so that  when A = 1 , if elflag = 1 then return the
   elastic cross section. If A = 1 and elflag = 0 then
   return the inelastic ep cross section. Previously qfs
   would return the sum of the elastic and inelastic. For
   inelastic ep events with W below the pion threshold the
   inelastic cross section would be zero but an elastic 
   cross section would be calculated.

4. Previously when dorg='y' there was a cut that W was near
   W = 938. Removed these cuts, though with new version
   elastic events should always have W_vertex = 938 so now
   elastic events would have always passed this cut.