---

ReactorEvent.cpp

#include <math.h>
#include <iostream.h>
#include "ReactorConstants.hh"
#include "ReactorFortran.hh"
#include "ReactorXC.hh"
#include "ReactorEvent.hh"

ReactorEvent::ReactorEvent(double newRmaxgen,
                           double newEmin,
                           double newEmax,
                           int newXCorder){
  Rmaxgen = newRmaxgen;
  Emin = newEmin;
  Emax = newEmax;
  XCorder = newXCorder;

  XCWeight = 1;
  
  double XCmax = InverseBeta(Emax,-1);
  double FluxMax = RMPflux(Emin);

  double Enu,z,phi;
  double R,zR,phiR;
  bool passed=0;


  while(!passed){
    float r[7];
    int len = 7;
    ranlux_(r,len);
    Enu = Emin+(Emax-Emin)*r[0];
    z = -1+2*r[1];
    phi = 2*k.pi*r[2];

    R = Rmaxgen*exp(log(r[3])/3);
    zR = -1+2*r[4];
    phiR = 2*k.pi*r[5];

    float XCtest = XCmax*FluxMax*r[6];
    XCWeight = InverseBeta(Enu,z);
    FluxWeight=RMPflux(Enu);
    passed= XCWeight*FluxWeight>XCtest;
  }
  

  double E0 = Enu - k.Delta;
  double p0 = sqrt(E0*E0-k.Melectron*k.Melectron);
  double v0 = p0/E0;
  double Ysquared = (k.Delta*k.Delta-k.Melectron*k.Melectron)/2;
  double E1 = E0*(1-Enu/k.Mproton*(1-v0*z))-Ysquared/k.Mproton;
  double p1 = sqrt(E1*E1-k.Melectron*k.Melectron);

  Pnubar[0] = Enu;
  Pnubar[1] = 0;
  Pnubar[2] = 0;
  Pnubar[3] = Enu;
  Pnubar[4] = Enu;
  Pnubar[5] = Enu;

  Ppositron[0] = E1;
  Ppositron[1] =  p1*sqrt(1-z*z)*cos(phi);
  Ppositron[2] =  p1*sqrt(1-z*z)*sin(phi);
  Ppositron[3] =  p1*z;
  Ppositron[4] =  E1-k.Melectron;
  Ppositron[5] =  p1;

  Pneutron[0] = Enu+k.Mproton-E1;
  Pneutron[1] = -Ppositron[1];
  Pneutron[2] = -Ppositron[2];
  Pneutron[3] = Pnubar[3]-Ppositron[3];
  Pneutron[4] = Pneutron[0]-k.Mproton-k.Delta;
  Pneutron[5] = sqrt(Pneutron[0]*Pneutron[0]-
                     (k.Mproton+k.Delta)*(k.Mproton+k.Delta));

  Rnubar[0] = 0;
  Rnubar[1] = R*sqrt(1-zR*zR)*cos(phiR);
  Rnubar[2] = R*sqrt(1-zR*zR)*sin(phiR);
  Rnubar[3] = R*zR;
  Rnubar[4] = R;
  Rnubar[5] = zR;
  Rnubar[6] = phiR;

  for(int n=0;n<Rdim;n++){
    Rneutron[n] = Rnubar[n];
    Rpositron[n] = Rnubar[n];
  }
}

ReactorEvent::~ReactorEvent(){
  ;
}

ReactorEvent::ReactorEvent(const ReactorEvent& Event){

  XCorder = Event.XCorder;  
  Emin = Event.Emin;
  Emax = Event.Emax;
  Rmaxgen = Event.Rmaxgen;

  XCWeight = Event.XCWeight;
  FluxWeight = Event.FluxWeight;

  for(int n=0;n<Pdim;n++){
    Pnubar[n]=Event.Pnubar[n];
    Pneutron[n]=Event.Pneutron[n];
    Ppositron[n]=Event.Ppositron[n];
  }
  for(int n=0;n<Rdim;n++){
    Rnubar[n]=Event.Rnubar[n];
    Rneutron[n]=Event.Rneutron[n];
    Rpositron[n]=Event.Rpositron[n];
  }
}    

ReactorEvent& ReactorEvent::operator=(const ReactorEvent& rhs){

  if(this == &rhs)return *this;

  XCorder = rhs.XCorder;  
  Emin = rhs.Emin;
  Emax = rhs.Emax;
  Rmaxgen = rhs.Rmaxgen;

  XCWeight = rhs.XCWeight;
  FluxWeight = rhs.FluxWeight;
  
  for(int n=0;n<Pdim;n++){
    Pnubar[n]=rhs.Pnubar[n];
    Pneutron[n]=rhs.Pneutron[n];
    Ppositron[n]=rhs.Ppositron[n];
  }
  for(int n=0;n<Rdim;n++){
    Rnubar[n]=rhs.Rnubar[n];
    Rneutron[n]=rhs.Rneutron[n];
    Rpositron[n]=rhs.Rpositron[n];
  }
  return *this;

}    

---