---

ReactorDetector.cpp

#include <math.h>
#include <iostream.h>
#include <string.h>
#include "ReactorEvent.hh"
#include "ReactorFortran.hh"
#include "ReactorConstants.hh"
#include "ReactorXC.hh"
#include "ReactorDetector.hh"
//
// constructor
//
ReactorDetector::ReactorDetector(){
  //
  // initialize PMT arrays
  //
  Xpmt = new double[0];
  Ypmt = new double[0];
  Zpmt = new double[0];
  PMTdata = new PMT[0];
  PMTpositron = new PMT[0];
  PMTneutron = new PMT[0];
  //
  // Set parameters
  //
  SetNominalParameters();
  SetPMT();
  //
  // initialize gamma arrays
  //
  Ngamma = 0;
  Rgamma = new double[Ngamma];
  Egamma = new double[Ngamma];
  Tgamma = new double[Ngamma];
  Ogamma = new char[Ngamma];

}
//
// destructor
//
ReactorDetector::~ReactorDetector(){
  delete[] Xpmt;
  delete[] Ypmt;
  delete[] Zpmt;
  delete[] PMTdata;
  delete[] PMTpositron;
  delete[] PMTneutron;
  delete[] Rgamma;
  delete[] Egamma;
  delete[] Tgamma;
  delete[] Ogamma;
}
//
// copy constructor
//
ReactorDetector::ReactorDetector(const ReactorDetector& Detector){

  R0=Detector.R0;
  R1=Detector.R1;
  R2=Detector.R2;
  GdConcentration=Detector.GdConcentration;
  refracGd=Detector.refracGd;
  refracSc=Detector.refracSc;
  mfpGd=Detector.mfpGd;
  mfpSc=Detector.mfpSc;
  MeanNeuDispl = Detector.MeanNeuDispl;
  FastNeutronOption = Detector.FastNeutronOption;
  Temperature = Detector.Temperature;
  tGd=Detector.tGd;
  tSc=Detector.tSc;
  GdCaptureFraction=Detector.GdCaptureFraction;
  GammasPerGd=Detector.GammasPerGd;
  attenlGd=Detector.attenlGd;
  attenlSc=Detector.attenlSc;
  PhotonsPerMeV=Detector.PhotonsPerMeV;
  ShortPosFrac0=Detector.ShortPosFrac0;
  ShortPosFrac1=Detector.ShortPosFrac1;
  LongPosFrac=Detector.LongPosFrac;
  ShortPosTime0=Detector.ShortPosTime0;
  ShortPosTime1=Detector.ShortPosTime1;
  LongPosTime=Detector.LongPosTime;
  PMTcoverage=Detector.PMTcoverage;
  PMTdiameter=Detector.PMTdiameter;
  PMTqe=Detector.PMTqe;
  ScintDecayTime=Detector.ScintDecayTime;

  Npmt = Detector.Npmt;
  Ncosbins = Detector.Ncosbins;
  Nphibins = Detector.Nphibins;

  Xpmt = new double[Npmt];
  Ypmt = new double[Npmt];  
  Zpmt = new double[Npmt];

  double* pold = Detector.Xpmt+Npmt;
  double* pnew = Xpmt+Npmt;
  while(pnew>Ypmt)*--pnew=*--pold;

  pold = Detector.Ypmt+Npmt;
  pnew = Ypmt+Npmt;
  while(pnew>Ypmt)*--pnew=*--pold;

  pold = Detector.Zpmt+Npmt;
  pnew = Zpmt+Npmt;
  while(pnew>Zpmt)*--pnew=*--pold;

  PHpositron = Detector.PHpositron;
  PHneutron = Detector.PHneutron;

  PHposReco = Detector.PHposReco;
  PHneuReco = Detector.PHneuReco;

  TimePositron = Detector.TimePositron;
  TimeNeutron = Detector.TimeNeutron;

  for(int i=0;i<3;i++){
    DipolePositron[i] = Detector.DipolePositron[i];
    DipoleNeutron[i] = Detector.DipoleNeutron[i];
  }

  for(int i=0;i<2;i++){
    MeanNpe[i]=Detector.MeanNpe[i];
    MeanAtten[i] = Detector.MeanAtten[i];
    MeanSfactor[i] = Detector.MeanSfactor[i];
  }

  Absorber = Detector.Absorber;
  NeutronSteps = Detector.NeutronSteps;

  Ngamma = Detector.Ngamma;
  Rgamma = new double[Ngamma];
  pold = Detector.Rgamma+3*Ngamma;
  pnew = Rgamma+3*Ngamma;
  while(pnew>Rgamma)*--pnew=*--pold;

  Egamma = new double[Ngamma];
  pold = Detector.Egamma+Ngamma;
  pnew = Egamma+Ngamma;
  while(pnew>Egamma)*--pnew=*--pold;

  Tgamma = new double[Ngamma];
  pold = Detector.Tgamma+Ngamma;
  pnew = Tgamma+Ngamma;
  while(pnew>Tgamma)*--pnew=*--pold;

  Ogamma = new char[Ngamma];
  char* rold = Detector.Ogamma+Ngamma;
  char* rnew = Ogamma+Ngamma;
  while(rnew>Ogamma)*--rnew=*--rold;

  PMTsize = Detector.PMTsize;
  PMTlen = Detector.PMTlen;
  PMTdata = new PMT[PMTsize];
  PMT* qold = Detector.PMTdata+PMTsize;
  PMT* qnew = PMTdata+PMTsize;
  while(qnew>PMTdata)*--qnew=*--qold;

  PMTpositron = new PMT[Npmt];
  qold = Detector.PMTpositron+Npmt;
  qnew = PMTpositron+Npmt;
  while(qnew>PMTpositron)*--qnew=*--qold;

  PMTneutron = new PMT[Npmt];
  qold = Detector.PMTneutron+Npmt;
  qnew = PMTneutron+Npmt;
  while(qnew>PMTneutron)*--qnew=*--qold;

}    
//
// = operator
//
ReactorDetector& ReactorDetector::operator=(const ReactorDetector& rhs){

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

  R0=rhs.R0;
  R1=rhs.R1;
  R2=rhs.R2;
  GdConcentration=rhs.GdConcentration;
  refracGd=rhs.refracGd;
  refracSc=rhs.refracSc;
  mfpGd=rhs.mfpGd;
  mfpSc=rhs.mfpSc;
  MeanNeuDispl = rhs.MeanNeuDispl;
  FastNeutronOption = rhs.FastNeutronOption;
  Temperature = rhs.Temperature;
  tGd=rhs.tGd;
  tSc=rhs.tSc;
  GdCaptureFraction=rhs.GdCaptureFraction;
  GammasPerGd=rhs.GammasPerGd;
  attenlGd=rhs.attenlGd;
  attenlSc=rhs.attenlSc;
  PhotonsPerMeV=rhs.PhotonsPerMeV;
  ShortPosFrac0=rhs.ShortPosFrac0;
  ShortPosFrac1=rhs.ShortPosFrac1;
  LongPosFrac=rhs.LongPosFrac;
  ShortPosTime0=rhs.ShortPosTime0;
  ShortPosTime1=rhs.ShortPosTime1;
  LongPosTime=rhs.LongPosTime;
  PMTcoverage=rhs.PMTcoverage;
  PMTdiameter=rhs.PMTdiameter;
  PMTqe=rhs.PMTqe;
  ScintDecayTime=rhs.ScintDecayTime;

  Npmt = rhs.Npmt;
  Ncosbins = rhs.Ncosbins;
  Nphibins = rhs.Nphibins;

  delete[] Xpmt;
  delete[] Ypmt;
  delete[] Zpmt;

  Xpmt = new double[Npmt];
  Ypmt = new double[Npmt];  
  Zpmt = new double[Npmt];

  double* pold = rhs.Xpmt+Npmt;
  double* pnew = Xpmt+Npmt;
  while(pnew>Ypmt)*--pnew=*--pold;

  pold = rhs.Ypmt+Npmt;
  pnew = Ypmt+Npmt;
  while(pnew>Ypmt)*--pnew=*--pold;

  pold = rhs.Zpmt+Npmt;
  pnew = Zpmt+Npmt;
  while(pnew>Zpmt)*--pnew=*--pold;

  PHpositron = rhs.PHpositron;
  PHneutron = rhs.PHneutron;

  PHposReco = rhs.PHposReco;
  PHneuReco = rhs.PHneuReco;

  TimePositron = rhs.TimePositron;
  TimeNeutron = rhs.TimeNeutron;

  for(int i=0;i<3;i++){
    DipolePositron[i] = rhs.DipolePositron[i];
    DipoleNeutron[i] = rhs.DipoleNeutron[i];
  }

  for (int i=0;i<2;i++){
    MeanNpe[i] = rhs.MeanNpe[i];
    MeanAtten[i] = rhs.MeanAtten[i];
    MeanSfactor[i] = rhs.MeanSfactor[i];
  }

  Absorber = rhs.Absorber;
  NeutronSteps = rhs.NeutronSteps;

  Ngamma = rhs.Ngamma;
  delete[] Rgamma;
  Rgamma = new double[Ngamma];
  pold = rhs.Rgamma+3*Ngamma;
  pnew = Rgamma+3*Ngamma;
  while(pnew>Rgamma)*--pnew=*--pold;

  delete[] Egamma;
  Egamma = new double[Ngamma];
  pold = rhs.Egamma+Ngamma;
  pnew = Egamma+Ngamma;
  while(pnew>Egamma)*--pnew=*--pold;

  delete[] Tgamma;
  Tgamma = new double[Ngamma];
  pold = rhs.Tgamma+Ngamma;
  pnew = Tgamma+Ngamma;
  while(pnew>Tgamma)*--pnew=*--pold;

  delete[] Ogamma;
  Ogamma = new char[Ngamma];
  char* rold = rhs.Ogamma+Ngamma;
  char* rnew = Ogamma+Ngamma;
  while(rnew>Ogamma)*--rnew=*--rold;

  PMTsize = rhs.PMTsize;
  PMTlen = rhs.PMTlen;
  delete[] PMTdata;
  PMTdata = new PMT[PMTsize];
  PMT* qold = rhs.PMTdata+PMTsize;
  PMT* qnew = PMTdata+PMTsize;
  while(qnew>PMTdata)*--qnew=*--qold;

  delete[] PMTpositron;
  PMTpositron = new PMT[Npmt];
  qold = rhs.PMTpositron+Npmt;
  qnew = PMTpositron+Npmt;
  while(qnew>PMTpositron)*--qnew=*--qold;

  delete[] PMTneutron;
  PMTneutron = new PMT[Npmt];
  qold = rhs.PMTneutron+Npmt;
  qnew = PMTneutron+Npmt;
  while(qnew>PMTneutron)*--qnew=*--qold;

  return *this;

}
//
// PMT sub-class  constructors,destructor,operator
//    
ReactorDetector::PMT::PMT(){
  npmt=0;
  for (int i=0;i<3;i++)xyz[i]=0; 
  TimeGen=0;
  TimeSmr=0;
  PHgen=0;
  PHsmr=0;
  PathGd=0;
  PathSc=0;
  Sfactor=0;
  Atten=0;
  Npe=0;
  Xpe=0;
  Parent='U';

  SubHits=0;
  SubHitsPtr = new PMT*[0];
}
  
ReactorDetector::PMT::~PMT(){
  delete[] SubHitsPtr;
}

ReactorDetector::PMT::PMT(const ReactorDetector::PMT& p){
  
  npmt=p.npmt; 
  for (int i=0;i<3;i++)xyz[i]=p.xyz[i]; 
  TimeGen=p.TimeGen;
  TimeSmr=p.TimeSmr;
  PHgen=p.PHgen;
  PHsmr=p.PHsmr;
  PathGd=p.PathGd;
  PathSc=p.PathSc;
  Sfactor=p.Sfactor;
  Atten=p.Atten;
  Npe=p.Npe;
  Xpe=p.Xpe;
  Parent=p.Parent;
  //
  SubHits=p.SubHits;
  SubHitsPtr = new PMT*[SubHits];
  PMT** pold = p.SubHitsPtr+SubHits;
  PMT** pnew = SubHitsPtr+SubHits;
  while(pnew>SubHitsPtr)*--pnew=*--pold;
}
//
// PMT sub-class equal operator
//    
ReactorDetector::PMT& 
ReactorDetector::PMT::operator=(const ReactorDetector::PMT& rhs){

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

  npmt=rhs.npmt;
  for (int i=0;i<3;i++)xyz[i]=rhs.xyz[i]; 
  TimeGen=rhs.TimeGen;
  TimeSmr=rhs.TimeSmr;
  PHgen=rhs.PHgen;
  PHsmr=rhs.PHsmr;
  PathGd=rhs.PathGd;
  PathSc=rhs.PathSc;
  Sfactor=rhs.Sfactor;
  Atten=rhs.Atten;
  Npe=rhs.Npe;
  Xpe=rhs.Xpe;
  Parent=rhs.Parent;
  //
  SubHits=rhs.SubHits;
  delete[] SubHitsPtr;
  SubHitsPtr = new PMT*[SubHits];
  PMT** pold = rhs.SubHitsPtr+SubHits;
  PMT** pnew = SubHitsPtr+SubHits;
  while(pnew>SubHitsPtr)*--pnew=*--pold;

  return *this;
}
ReactorDetector::PMT 
ReactorDetector::PMT::operator+(const ReactorDetector::PMT& rhs){

  PMT sum(*this);
  if(npmt!=rhs.npmt)return sum;
  sum.npmt = npmt;

  double w1=Xpe/(Xpe+rhs.Xpe);
  double w2=1-w1;

  for (int i=0;i<3;i++)sum.xyz[i]=w1*xyz[i]+w2*rhs.xyz[i]; 
  sum.TimeGen=w1*TimeGen+w2*rhs.TimeGen;
  sum.TimeSmr=w1*TimeSmr+w2*rhs.TimeSmr;
  sum.PHgen=w1*PHgen+w2*rhs.PHgen;
  sum.PHsmr=w1*PHsmr+w2*rhs.PHsmr;
  sum.PathGd=w1*PathGd+w2*rhs.PathGd;
  sum.PathSc=w1*PathSc+w2*rhs.PathSc;
  sum.Sfactor=w1*Sfactor+w1*rhs.Sfactor;
  sum.Atten=w1*Atten+w2*rhs.Atten;
  sum.Npe=Npe+rhs.Npe;
  sum.Xpe=Xpe+rhs.Xpe;
  if(Parent==rhs.Parent)sum.Parent=rhs.Parent;
  else sum.Parent='U';

  sum.SubHits = SubHits+rhs.SubHits;
  delete[] sum.SubHitsPtr;
  sum.SubHitsPtr = new PMT*[sum.SubHits];
  PMT** p = sum.SubHitsPtr+sum.SubHits;
  PMT** q = SubHitsPtr+SubHits;
  PMT** r = rhs.SubHitsPtr+rhs.SubHits;
  while(p>sum.SubHitsPtr+SubHits)*--p=*--r;
  while(p>sum.SubHitsPtr)*--p=*--q;

  return sum;
}
ReactorDetector::PMT& 
ReactorDetector::PMT::operator+=(const ReactorDetector::PMT& rhs){

  if(npmt!=rhs.npmt)return *this;

  double w1=Xpe/(Xpe+rhs.Xpe);
  double w2=1-w1;

  for (int i=0;i<3;i++)xyz[i]=w1*xyz[i]+w2*rhs.xyz[i]; 
  TimeGen=w1*TimeGen+w2*rhs.TimeGen;
  TimeSmr=w1*TimeSmr+w2*rhs.TimeSmr;
  PHgen=w1*PHgen+w2*rhs.PHgen;
  PHsmr=w1*PHsmr+w2*rhs.PHsmr;
  PathGd=w1*PathGd+w2*rhs.PathGd;
  PathSc=w1*PathSc+w2*rhs.PathSc;
  Sfactor=w1*Sfactor+w1*rhs.Sfactor;
  Atten=w1*Atten+w2*rhs.Atten;
  Npe+=rhs.Npe;
  Xpe+=rhs.Xpe;
  if(Parent!=rhs.Parent)Parent='U';
  //
  PMT** SubHitsTmp = new PMT*[SubHits];
  PMT** p = SubHitsTmp+SubHits;
  PMT** q = SubHitsPtr+SubHits;
  while(p>SubHitsTmp)*--p=*--q;
  delete[] SubHitsPtr;
  SubHitsPtr = new PMT*[SubHits+rhs.SubHits];
  p = SubHitsPtr+SubHits+rhs.SubHits;
  q = SubHitsTmp+SubHits;
  PMT** r = rhs.SubHitsPtr+rhs.SubHits;
  while(p>SubHitsPtr+SubHits)*--p=*--r;
  while(p>SubHitsPtr)*--p=*--q;
  SubHits += rhs.SubHits;
  delete[] SubHitsTmp;
  //
  return *this;
}
void ReactorDetector::PMT::SetSubHitsPtr(int subhits){
  delete[] SubHitsPtr;
  SubHits = subhits;
  SubHitsPtr = new PMT*[subhits];
  return;
}


//
// Detector simulation function
//
void ReactorDetector::LightsOut(ReactorEvent& Event){

  PMTlen=0;
  delete[] PMTdata;
  PMTdata = new PMT[PMTsize];
  
  GenerateSecondaries(Event);
  GenerateGammas(Event);
  GenerateResponse(Event);

  ConsolidatePMT();

  CalculateDipoleMoment('P');
  CalculateDipoleMoment('N');

  ImproveDipoleMoment('P');
  ImproveDipoleMoment('N');
}
//
// Access to PMT data
//
void ReactorDetector::GetPMTdata(int npmt,double* XYZpmt,
                                 int& npos,double* PHpos,double* Tpos,
                                 int& nneu,double* PHneu,double* Tneu){
  
  *(XYZpmt+0) = *(Xpmt+npmt);
  *(XYZpmt+1) = *(Ypmt+npmt);
  *(XYZpmt+2) = *(Zpmt+npmt);

  npos=0;
  nneu=0;
  for(PMT* PMTptr=PMTdata;PMTptr<PMTdata+PMTlen;PMTptr++){
    if(npmt==(*PMTptr).npmt){
      if((*PMTptr).Parent=='P'){
        *(PHpos+npos) = (*PMTptr).PHsmr;
        *(Tpos+npos) = (*PMTptr).TimeSmr;
        npos++;
      }
      if((*PMTptr).Parent=='N'){
        *(PHneu+nneu) = (*PMTptr).PHsmr;
        *(Tneu+nneu) = (*PMTptr).TimeSmr;
        nneu++;
      }
    } 
  }
  return;
}
void ReactorDetector::GenerateSecondaries(ReactorEvent& Event){

  //
  // Displace the positron vertex by its range along the positron direction.
  //
  double kepos = Event.GetPositronKE();
  double rangepos;
  double tpos;
  PositronRange(kepos,rangepos,tpos);
  //
  double cospos = Event.GetPositronCosTheta();
  double phipos = Event.GetPositronPhi();
  double rint[3];
  for(int i=0;i<3;i++)rint[i] = *(Event.GetNubarVertexXYZ()+i);
  //
  double rpos[3];
  rpos[0]=rint[0]+rangepos*sqrt(1-cospos*cospos)*cos(phipos);
  rpos[1]=rint[1]+rangepos*sqrt(1-cospos*cospos)*sin(phipos);
  rpos[2]=rint[2]+rangepos*cospos;
  double rrpos = sqrt(rpos[0]*rpos[0]+rpos[1]*rpos[1]+rpos[2]*rpos[2]);
  //
  //
  //  
  double rneu[3];
  double tneu;
  //
  if(FastNeutronOption){
    //
    // Throw the neutron vertex via a simple Gaussian
    //
    float rn[3];
    int len=3;
    rnormx_(rn,len,ranlux_);
    //
    // Random for throwing neutron time
    //
    len=1;
    float rv[1];
    ranlux_(rv,len);
    //
    //  Simple scheme that is not right at boundaries
    //  Neglects discontinuity in mean free path;
    //  Time will be taken as protportional to sqrt(path)
    //
    double rrint = sqrt(rint[0]*rint[0]+rint[1]*rint[1]+rint[2]*rint[2]);
    double path;
    if(rrint<=R0){
      for (int i=0;i<3;i++)rneu[i]=rint[i]+mfpGd*rn[i]/sqrt(3.);
      path = mfpGd*sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
      tneu = tGd*sqrt(path/mfpGd)*log(1/rv[0]);
    }
    else if(rrint<R2){
      for (int i=0;i<3;i++)rneu[i]=rint[i]+mfpSc*rn[i]/sqrt(3.);
      path = mfpSc*sqrt(rn[0]*rn[0]+rn[1]*rn[1]+rn[2]*rn[2]);
      tneu = tSc*sqrt(path/mfpSc)*log(1/rv[0]);
    }
    else{
      for (int i=0;i<3;i++)rneu[i]=999999;
      tneu = 999999.;
    }
    //
    // Shift z displacement by mean Chooz value
    //
    rneu[2]+=MeanNeuDispl;
  }
  else{
    double pneu[3];
    for(int i=0;i<3;i++)pneu[i]=*(Event.GetNeutronPxyz()+i);
    for(int i=0;i<3;i++)rneu[i]=rint[i];
    double ke = Event.GetNeutronKE();
    NeutronSteps=0;
    Absorber=0;
    NeutronPropagator(R0,R2,ke,Absorber,NeutronSteps,pneu,rneu,tneu);
  }
  //
  //  subevent times
 //
  double T[2] = {tpos,tneu};
  TimePositron = T[0];
  TimeNeutron = T[1];
  //
  // set positron,neutron vertex in Event class
  //
  double Rpos[7];
  Rpos[0] = tpos;
  for(int i=1;i<4;i++)Rpos[i]=rpos[i-1];
  Rpos[4] = rrpos;
  Rpos[5] = rpos[2]/rrpos;
  Rpos[6] = atan2(rpos[1],rpos[0]);
  if(Rpos[6]<0)Rpos[6]+=2*RC.pi;
  Event.SetPositronVertex(Rpos);
  //
  double Rneu[7];
  Rneu[0] = tneu;
  for(int i=1;i<4;i++)Rneu[i]=rneu[i-1];
  double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
  Rneu[4] = rrneu;
  Rneu[5] = rneu[2]/rrneu;
  Rneu[6] = atan2(rneu[1],rneu[0]);
  if(Rneu[6]<0)Rneu[6]+=2*RC.pi;
  Event.SetNeutronVertex(Rneu);
}

void ReactorDetector::GenerateResponse(ReactorEvent& Event){

  PHpositron = 0;
  PHneutron = 0;
  PHposReco = 0;
  PHneuReco = 0;
  for(int i=0;i<2;i++){
    MeanSfactor[i] = 0;
    MeanAtten[i] = 0;
    MeanNpe[i] = 0;
  }
  //
  //
  //
  int pmtindex=0;
  double* xpmt=Xpmt;
  double* ypmt=Ypmt;
  double* zpmt=Zpmt;
  
  PMT* PMTptr = PMTdata;

  for(int i=0;i<Ncosbins;i++){
    for(int j=0;j<Nphibins;j++){
      double rpmt[3];
      rpmt[0] = *xpmt++;
      rpmt[1] = *ypmt++;
      rpmt[2] = *zpmt++;
      double r[3];
      double dr[3];
      double tr[3];
      //      
      double sumpos[6]={0,0,0,0,0,0};
      double sumneu[6]={0,0,0,0,0,0};
      //
      double* rgptr=Rgamma;
      double* egptr=Egamma;
      double* tgptr=Tgamma;
      char* ogptr=Ogamma;
      //
      int nposgt0=0;
      int nneugt0=0;
      //
      for(int k=0;k<Ngamma;k++){
        //
        for (int l=0;l<3;l++)r[l]=*rgptr++;
        double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
        //
        // nothing to do if r>R1
        //
        if(sqrt(rdotr)>R1)continue;
        //
        for (int l=0;l<3;l++)dr[l]=rpmt[l]-r[l];
        double ddr = sqrt(dr[0]*dr[0]+dr[1]*dr[1]+dr[2]*dr[2]);
        for (int l=0;l<3;l++)tr[l]=dr[l]/ddr;
        double tdotr = r[0]*tr[0]+r[1]*tr[1]+r[2]*tr[2];
        //
        // solve for intersection with inner radius
        //
        double sGd;
        double sSc;
        //
        // check for no intersections
        //
        double test = R0*R0+tdotr*tdotr-rdotr;
        bool inter = (test>=0);
        //
        // A straight line from the interaction vertex to the pmt
        // goes through the inner radius at least once.  In this case,
        // the light will see Gd-scint. and unloaded scint.
        //
        if(inter){
          double sint1 = -tdotr+sqrt(test);
          double sint2 = -tdotr-sqrt(test);
          if(sint1>0&&sint2<0){     // 1 intersection at soln 1
            sGd = sint1;
          }
          else if(sint2>0&&sint1<0){     // 1 intersection at soln 2
            sGd = sint2;
          }
          else if(sint2>0&&sint1>0){     // 2 intersections 
            sGd = sint1-sint2;
            if(sGd<0)sGd*=-1;
          }
          else{                          // no intersections towards PMT
            sGd = 0;
          }
          sSc = ddr-sGd;
        }
        //
        //  otherwise the light only goes through unloaded scintillator
        //
        else{
          sGd = 0;
          sSc = ddr;
        }
        //
        //  increment PH-- make atten. relative to center
        //
        double dsGd = sGd-R0;
        double dsSc = sSc-R2+R0;
        double atten0 = exp(-sGd/attenlGd-sSc/attenlSc);     
        double atten = exp(-dsGd/attenlGd-dsSc/attenlSc);     
        //
        // solid angle factor
        //
        double sfactor=0;
        for (int l=0;l<3;l++){sfactor+=tr[l]*rpmt[l]/R2;}
        sfactor *= R2*R2/ddr/ddr;
        //
        // photostatistics
        //
        double phgen = *egptr++;
        //
        float xpe = (phgen*PhotonsPerMeV);
        xpe *= (atten0*sfactor);
        xpe *= (PMTcoverage/Npmt*PMTqe);
        //
        int npe;
        int error;
        rnpssn_(xpe,npe,error);
        //
        double phsmr = npe/PMTqe/PMTcoverage/PhotonsPerMeV;     
        //        
        double tgen = *tgptr++;
        double tsmr = refracGd*sGd/RC.c+refracSc*sSc/RC.c + tgen; 
        //
        if(*ogptr++ == 'P'){
          nposgt0++;
          sumpos[0]+=phsmr;
          sumpos[1]+=phgen;
          sumpos[2]+=phgen*atten0*sfactor;
          sumpos[3]+=npe;
          sumpos[4]+=sfactor;
          sumpos[5]+=atten;
          //
          PHpositron+=phsmr;
        }
        else{
          nneugt0++;
          sumneu[0]+=phsmr;
          sumneu[1]+=phgen;
          sumneu[2]+=phgen*atten0*sfactor;
          sumneu[3]+=npe;
          sumneu[4]+=sfactor;
          sumneu[5]+=atten;
          //      
          PHneutron+=phsmr;
        }
        //        
        PMT PMTtmp;
        PMTtmp.npmt = pmtindex;
        for(int l=0;l<3;l++)PMTtmp.xyz[l] = rpmt[l];
        PMTtmp.PHgen = phgen;
        PMTtmp.PHsmr = phsmr;
        PMTtmp.PathGd = sGd;
        PMTtmp.PathSc = sSc;
        PMTtmp.Atten = atten;
        PMTtmp.Sfactor = sfactor;
        PMTtmp.TimeGen = tgen;
        PMTtmp.TimeSmr = tsmr;
        PMTtmp.Xpe = xpe;
        PMTtmp.Npe = npe;
        PMTtmp.Parent = *ogptr;
        
        if(PMTptr<PMTdata+PMTsize){
          *PMTptr++=PMTtmp;
          PMTlen++;
        }
        else{
          //
          // expand the data array first
          //
          PMTptr = ResizePMT();
          *PMTptr++=PMTtmp;
          PMTlen++;
        }
      }
      //
      //  Crude estimate of reconstructed PH
      //  Includes poisson npe smearing and fluctuation in atten length
      //
      if(sumpos[2]>0){
        PHposReco+=sumpos[0]*sumpos[1]/sumpos[2];
        MeanNpe[0]+=sumpos[3]/Npmt;
        MeanSfactor[0]+=sumpos[4]/Npmt/nposgt0;
        MeanAtten[0]+=sumpos[5]/Npmt/nposgt0;
      }
      if(sumneu[2]>0){
        PHneuReco+=sumneu[0]*sumneu[1]/sumneu[2];
        MeanNpe[1]+=sumneu[3]/Npmt;
        MeanSfactor[1]+=sumneu[4]/Npmt/nneugt0;
        MeanAtten[1]+=sumneu[5]/Npmt/nneugt0;
      }
      //
      // increment PMT index
      //
      pmtindex++;
    }
  }
}
void ReactorDetector::GenerateGammas(ReactorEvent& Event){
  //
  int maxgamma = 25;
  int ngamma=0;
  double* egamma = new double[maxgamma];
  double* tgamma = new double[maxgamma];
  double* rgamma = new double[3*maxgamma];
  char* ogamma = new char[maxgamma];
  //
  double* egptr=egamma;
  double* rgptr=rgamma;
  double* tgptr=tgamma;
  char* ogptr=ogamma;
  //
  // get positron, neutron vertices
  //
  double rpos[3];
  double rneu[3];
  for(int i=0;i<3;i++)rpos[i] = *(Event.GetPositronVertexXYZ()+i);
  for(int i=0;i<3;i++)rneu[i] = *(Event.GetNeutronVertexXYZ()+i);
  double tpos = Event.GetPositronTime();
  double tneu = Event.GetNeutronTime();
  double kepos = Event.GetPositronKE();
  //
  // positron gammas first
  // the first is a "fake" gamma with E=positron KE generated at positron
  // vertex.
  //
  for(int i=0;i<3;i++)*rgptr++=rpos[i];
  *tgptr++ = tpos;
  *egptr++ = kepos;
  *ogptr++ = 'P';
  ngamma++;
  //
  // Now the two annihilation gammas
  //
  int len=4;
  float rv[4];
  ranlux_(rv,len);
  double cosg = -1+2*rv[0];
  double phig = 2*RC.pi*rv[1];
  double t[3];
  t[0] = cos(phig)*sqrt(1-cosg*cosg);
  t[1] = sin(phig)*sqrt(1-cosg*cosg);
  t[2] = cosg;
  double range;
  for (int i=0;i<2;i++){
    range = GammaComptonRange(RC.Melectron);
    range*=log(1/rv[i+2]);
    if(i==0){
      for(int j=0;j<3;j++)*rgptr++ = rpos[j]+range*t[j];
    }
    if(i==1){
      for(int j=0;j<3;j++)*rgptr++ = rpos[j]-range*t[j];
    }
    *egptr++=RC.Melectron;
    *tgptr++=tpos+range/RC.c*refracSc;; 
    *ogptr++ = 'P';
    ngamma++;
  }
  //
  // Generate neutron gammas
  //
  //
  // pick H or Gd, then which Gd peak
  //
  double emax=RC.Hpeak;
  if(FastNeutronOption){
    Absorber=1;
    double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
    if(rrneu<R0){//in Gd region
      len=2;
      ranlux_(rv,len);
      if(rv[0]<=GdCaptureFraction){
        //
        // model Gd spectrum by Gd155 + Gd157
        //
        if(rv[1]<=RC.Gd155frac){
          emax = RC.Gd155peak;
          Absorber=155;
        }
        else{
          emax = RC.Gd157peak;
          Absorber=157;
        }
      }
    }
  }
  else{
    if(Absorber==1)emax=RC.Hpeak;
    if(Absorber==12)emax=RC.Cpeak;
    if(Absorber==155)emax=RC.Gd155peak;
    if(Absorber==157)emax=RC.Gd157peak;
  }
  //
  // One photon from an H or C capture
  //
  if(emax==RC.Hpeak || emax==RC.Cpeak){
    len=3;
    ranlux_(rv,len);
    range = GammaComptonRange(emax);
    range*=log(1/rv[0]);
    cosg = -1+2*rv[1];
    phig = 2*RC.pi*rv[2];
    t[0] = cos(phig)*sqrt(1-cosg*cosg);
    t[1] = sin(phig)*sqrt(1-cosg*cosg);
    t[2] = cosg;
    for(int j=0;j<3;j++)*rgptr++ = rneu[j]+range*t[j];
    *egptr++ = emax;
    *tgptr++ = tneu+range/RC.c*refracSc;
    *ogptr++ = 'N';
    ngamma++;
  }
  //
  // otherwise model Gd photons
  //
  else{
    int ngd;
    double pgd[3*maxgamma];
    MyGdDecayModel(emax,ngd,pgd);
    float rv2[ngd];
    ranlux_(rv2,ngd);
    double* qgd=pgd;
    double check[4]={0,0,0,0};
    for(int i=0;i<ngd;i++){
      double p[3];
      for(int j=0;j<3;j++)p[j]=*qgd++;
      double e = sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
      range = GammaComptonRange(e);
      range*=log(1/rv[i]);
      for(int j=0;j<3;j++)*rgptr++ = rneu[j]+range*p[j]/e;
      //
      for(int j=0;j<3;j++)check[j]+=p[j];
      check[3]+=e;
      //
      *egptr++ = e;
      *tgptr++ = tneu+range/RC.c*refracGd;
      *ogptr++ = 'N';
      ngamma++;
    }
  }
  //
  // Transfer temp data to storage
  //
  delete[] Rgamma;
  delete[] Egamma;
  delete[] Tgamma;
  delete[] Ogamma;
  
  Ngamma = ngamma;
  Rgamma = new double[3*ngamma];
  Egamma = new double[ngamma];
  Tgamma = new double[ngamma];
  Ogamma = new char[ngamma];

  rgptr=rgamma+3*ngamma;
  double* ptr = Rgamma+3*ngamma;
  while(ptr>Rgamma)*--ptr=*--rgptr;

  egptr=egamma+ngamma;
  ptr = Egamma+ngamma;
  while(ptr>Egamma)*--ptr=*--egptr;

  tgptr=tgamma+ngamma;
  ptr = Tgamma+ngamma;
  while(ptr>Tgamma)*--ptr=*--tgptr;

  ogptr=ogamma+ngamma;
  char* qtr = Ogamma+ngamma;
  while(qtr>Ogamma)*--qtr=*--ogptr;

  delete[] rgamma;
  delete[] egamma;
  delete[] tgamma;
  delete[] ogamma;
    
}


ReactorDetector::PMT* ReactorDetector::ResizePMT(){

  int oldsize = PMTsize;
  int newsize = 2*oldsize;
  PMT* PMThold = new PMT[newsize];
  PMT* PMTold = PMTdata+oldsize;
  PMT* PMTnew = PMThold+oldsize;
  while(PMTnew>PMThold)*--PMTnew=*--PMTold;
  delete[] PMTdata;
  PMTdata = new PMT[newsize];
  PMTold = PMThold+newsize;
  PMTnew = PMTdata+newsize;
  while(PMTnew>PMTdata)*--PMTnew=*--PMTold;
  delete[] PMThold;
  PMTsize = newsize;
  return PMTdata+oldsize;
}
//
// Fetch number of gammas from positron
//
int ReactorDetector::GetNgamma(char parent){
  int gamma=0;
  for(char* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
    if (*optr==parent)gamma++;
  }
  return gamma;
}
//
// Fetch array of positron gammas
//
void ReactorDetector::GetGamma(char parent,double* xyz){

  double* Rptr=Rgamma;
  double* rptr=xyz;
  //
  for(char* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
    if (*optr==parent){
      for(int i=0;i<3;i++)*rptr++ = *Rptr++;
    }
  }
}  
void ReactorDetector::GetGamma(char parent,double* e,
                                  double* t,double* r,
                                  double* c,double* f){
  double* Eptr=Egamma;
  double* Rptr=Rgamma;
  double* Tptr=Tgamma;
  double* eptr=e;
  double* rptr=r;
  double* tptr=t;
  double* cptr=c;
  double* fptr=f;
  //
  for(char* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
    if (*optr==parent){
      *e++=*Eptr;
      *t++=*Tptr;
      double x=*(Rptr+0);
      double y=*(Rptr+1);
      double z=*(Rptr+2);
      *r++=sqrt(x*x+y*y+z*z);
      *c++=z/(*r);
      *f++=atan2(y,x);
      if((*f)<0)*f+=2*RC.pi;
    }
    Eptr++;
    Tptr++;
    Rptr+=3;
  }
}
void ReactorDetector::SetNominalParameters(){
  //
  // Set nominal detector parameters from ReactorConstants
  //
  R0 = RC.R0;
  R1 = RC.R1;
  R2 = RC.R2;
  //
  GdConcentration=RC.GdConcentration;
  //
  refracGd = RC.refracGd;
  refracSc = RC.refracSc;
  //
  mfpGd = RC.mfpGd; 
  mfpSc = RC.mfpSc;
  //
  MeanNeuDispl = RC.MeanNeuDispl;
  //
  FastNeutronOption = RC.FastNeutronOption;
  Temperature = RC.Temperature;
  //
  tGd = RC.tGd;
  tSc = RC.tSc;
  //
  GdCaptureFraction = RC.GdCaptureFraction;
  //
  // This is an ad-hoc guess for the average number of gammas/per Gd.
  // The program in practice generates at least two gammas per capture
  // and an additional Poisson fluctuated GammasPerGd-2.
  //
  GammasPerGd = RC.GammasPerGd;
  //
  attenlGd = RC.attenlGd;
  attenlSc = RC.attenlSc;
  //
  PhotonsPerMeV = RC.PhotonsPerMeV;
  //
  ShortPosFrac0 = RC.ShortPosFrac0;
  ShortPosFrac1 = RC.ShortPosFrac1;
  LongPosFrac = RC.LongPosFrac;
  ShortPosTime0 = RC.ShortPosTime0;
  ShortPosTime1 = RC.ShortPosTime1;
  LongPosTime = RC.LongPosTime;
  //
  PMTcoverage = RC.PMTcoverage;
  PMTdiameter = RC.PMTdiameter;
  //
  PMTqe = RC.PMTqe;
}
void ReactorDetector::SetPMT(){
  //
  // determine the number of PMT
  //

  Npmt = 16*R2*R2*PMTcoverage/PMTdiameter/PMTdiameter;
  Ncosbins = sqrt(Npmt/RC.pi);
  Nphibins = sqrt(Npmt*RC.pi);
  Npmt = Ncosbins*Nphibins;

  delete[] Xpmt;
  delete[] Ypmt;
  delete[] Zpmt;

  Xpmt = new double[Npmt];
  Ypmt = new double[Npmt];
  Zpmt = new double[Npmt];
  
  //
  // Make a guess at PMT size to avoid too much re-sizing
  //
  delete[] PMTdata;
  PMTlen = 0;
  PMTsize = 10*Npmt;
  PMTdata = new PMT[PMTsize];
  //
  // loop over cotheta, then phi
  //
  double* xpmt=Xpmt;
  double* ypmt=Ypmt;
  double* zpmt=Zpmt;
  
  double dphi = 2*RC.pi/Nphibins;
  double dz = 2.0/Ncosbins;
  
  for(int i=0;i<Ncosbins;i++){
    double z = -1+i*dz+dz/2;
    for(int j=0;j<Nphibins;j++){
      double phi = j*dphi+dphi/2;
      *xpmt++ = R2*sqrt(1-z*z)*cos(phi);
      *ypmt++ = R2*sqrt(1-z*z)*sin(phi);
      *zpmt++ = R2*z;
      //      cout << "z,phi  " << z << " " << phi << endl;
    }
  }     
}
void ReactorDetector::SetParam_GdConcentration(double x){
  GdConcentration =x;
  double y =(GdConcentration/RC.GdConcentrationRef);
  double z = GdCaptureFraction/(1-GdCaptureFraction);
  //
  GdCaptureFraction = y*z/(1+y*z);
  //
  mfpGd*=(1+z)/(1+y*z);
  tGd*=(1+z)/(1+y*z);
}
double ReactorDetector::GetHitFrac(char parent){
  bool hit[Npmt];
  for(int i=0;i<Npmt;i++)hit[i]=0;
  for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
    if(pPMT->Parent==parent&&pPMT->Npe>0)hit[pPMT->npmt]=1;
  }
  int nhit=0;
  for(int i=0;i<Npmt;i++)nhit+=hit[i];
  return 1.0*nhit/Npmt;
}
void ReactorDetector::CalculateDipoleMoment(char parent){

  double d[3] = {0,0,0};
  double q=0;
  for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
    if(pPMT->Parent==parent&&pPMT->Npe>0){
      q+=pPMT->PHsmr;
      for(int i=0;i<3;i++)d[i]+=(pPMT->PHsmr)*(pPMT->xyz[i]);
    }
  }
  if(q>0)for(int i=0;i<3;i++)d[i]/=q;

  if(parent=='P'){
    for(int i=0;i<3;i++)DipolePositron[i] = d[i];
  }
  if(parent=='N'){
    for(int i=0;i<3;i++)DipoleNeutron[i] = d[i];
  }
}
void ReactorDetector::ImproveDipoleMoment(char parent){

  double d[3];
  double k;
  PMT* qPMT;
  if(parent=='P'){
    for(int i=0;i<3;i++)d[i]=DipolePositron[i];
    k = PHpositron;
    qPMT=PMTpositron;
  }
  if(parent=='N'){
    for(int i=0;i<3;i++)d[i]=DipoleNeutron[i];
    k = PHneutron;
    qPMT=PMTneutron;
  }
  //
  double C= RC.PMTqe*RC.PMTdiameter*RC.PMTdiameter/16*RC.PhotonsPerMeV;
  //
  double sumlast = 999999;
  double test=999999;
  double tol=0.01;
  int sign=1;
  bool first=1;
  double dev=0.1;
  while(test>tol){
  //
    double dlast[3];
    for(int i=0;i<3;i++)dlast[i]=d[i];
    double klast=k;
    double sum[3]={0,0,0};
    for(PMT* pPMT=qPMT;pPMT<qPMT+Npmt;pPMT++){
      if(pPMT->Parent!=parent)continue;
      double r[3];
      double rn[3];
      double t[3];
      for(int i=0;i<3;i++)rn[i]=pPMT->xyz[i];
      for(int i=0;i<3;i++)r[i]=rn[i]-d[i];
      double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
      double rdotrn = r[0]*rn[0]+r[1]*rn[1]+r[2]*rn[2];
      for(int i=0;i<3;i++)t[i]=r[i]/sqrt(rdotr);
      //
      double tdotd=d[0]*t[0]+d[1]*t[1]+d[2]*t[2];
      double ddotd=d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
      double test2=tdotd*tdotd-ddotd+R0*R0;
      double sSc;
      double sGd;
      //
      if(test2<0){//no intersections with inner shell
        sSc=sqrt(rdotr);
        sGd=0;
      }
      else{
        double s1=-tdotd+sqrt(test2);
        double s2=-tdotd-sqrt(test2);
        if(s1>0&&s2>0){//2 interesections
          sSc=sqrt(rdotr)-s1+s2;
          sGd=s1-s2;
        }
        else if(s1>0&&s2<0){//1 interestion
          sSc=sqrt(rdotr)-s1;
          sGd=s1;
        }
        else {//no interesctions
          sSc=sqrt(rdotr);
          sGd=0;
        }
      }
      //
      double mn=pPMT->Npe;
      double w0=C*rdotrn/R2/sqrt(rdotr)/rdotr*
        exp(-sGd/RC.attenlGd-sSc/RC.attenlSc);
      sum[0]+=mn;
      sum[1]+=w0;
      sum[2]+=w0*klast-mn*log(w0*klast);
    }
    if(first){
      k = sum[0]/sum[1];
      for(int i=0;i<3;i++)d[i]*=(1+dev);
      first=0;
    }
    else{
      k = sum[0]/sum[1];
      if(sum[2]<sumlast){
        for(int i=0;i<3;i++)d[i]*=(1+sign*dev);
      }
      else{
        sign*=-1;
        dev/=2;
        for(int i=0;i<3;i++)d[i]*=(1+sign*dev);
      }
    }
    test=1-sum[2]/sumlast;
    if(test<0)test*=-1;
    sumlast=sum[2];
  }
  if(parent=='P')for(int i=0;i<3;i++)DipolePositron[i]=d[i];
  if(parent=='N')for(int i=0;i<3;i++)DipoleNeutron[i]=d[i];
}  
//
//
//
void ReactorDetector::ConsolidatePMT(){
  //
  // first sum the total xpe
  //
  double xpepos[Npmt];
  double xpeneu[Npmt];
  int subhitspos[Npmt];
  int subhitsneu[Npmt];
  //
  for(int i=0;i<Npmt;i++)xpepos[i]=0;
  for(int i=0;i<Npmt;i++)xpeneu[i]=0;
  for(int i=0;i<Npmt;i++)subhitspos[i]=0;
  for(int i=0;i<Npmt;i++)subhitsneu[i]=0;
  //
  for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
    int n=pPMT->npmt;
    if(pPMT->Parent=='P'){
      xpepos[n]+=pPMT->Npe;
      subhitspos[n]++;
    }
    if(pPMT->Parent=='N'){
      xpeneu[n]+=pPMT->Npe;
      subhitsneu[n]++;
    }
  }
  //   
  delete[] PMTpositron;
  PMTpositron = new PMT[Npmt];
  delete[] PMTneutron;
  PMTneutron = new PMT[Npmt];
  //
  for(int i=0;i<Npmt;i++){
    (*(PMTpositron+i)).npmt=i;
    (*(PMTpositron+i)).Parent='P';
    (PMTpositron+i)->SetSubHitsPtr(subhitspos[i]);
    //
    (*(PMTneutron+i)).npmt=i;
    (*(PMTneutron+i)).Parent='N';
    (PMTneutron+i)->SetSubHitsPtr(subhitsneu[i]);
  }
  //
  for(PMT* pPMT=PMTdata;pPMT<PMTdata+PMTlen;pPMT++){
    int n=pPMT->npmt;
    PMT* qPMT;
    if(pPMT->Parent=='P'){
      qPMT = PMTpositron+n;
      *qPMT+=*pPMT;
    }
    if(pPMT->Parent=='N'){
      qPMT = PMTneutron+n;
      *qPMT+=*pPMT;
    }
  }
  cout;
}    

---