---

ReactorDetector.cpp

Go to the documentation of this file.

#include <math.h>
#include <iostream.h>
#include <string.h>
#include "ReactorEvent.hh"
#include "ReactorFortran.hh"
#include "ReactorConstants.hh"
#include "ReactorXC.hh"
#include "ReactorDetector.hh"
#include "MyTlSource.hh"
#include "MyBiSource.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();
  SetNeutronPropStuff();
  //
  // initialize gamma arrays
  //
  Ngamma = 0;
  Rgamma = new double[Ngamma];
  Egamma = new double[Ngamma];
  Tgamma = new double[Ngamma];
  Ogamma = new std::string[Ngamma];

  Nneutron = 0;
  Absorber = new double[Nneutron];
  NeutronSteps = new double[Nneutron];

  pmtRadCount = 0;
}
//
// destructor
//
ReactorDetector::~ReactorDetector(){
  delete[] Xpmt;
  delete[] Ypmt;
  delete[] Zpmt;
  delete[] PMTdata;
  delete[] PMTpositron;
  delete[] PMTneutron;
  delete[] Rgamma;
  delete[] Egamma;
  delete[] Tgamma;
  delete[] Ogamma;
  delete[] Absorber;
  delete[] NeutronSteps;
}
//
// 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;
  for(int i=0;i<4;i++){
    a[i] = Detector.a[i];
    z[i] = Detector.z[i];
    nn[i] = Detector.nn[i];
  }

  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;

  PHelectron = Detector.PHelectron;
  PHpositron = Detector.PHpositron;
  PHneutron = Detector.PHneutron;
  PHthallium = Detector.PHthallium;
  PHbismuth = Detector.PHbismuth;
  PHcalifornium = Detector.PHcalifornium;
  PHmuon = Detector.PHmuon;

  PHeleReco = Detector.PHeleReco;
  PHposReco = Detector.PHposReco;
  PHneuReco = Detector.PHneuReco;
  PHthaReco = Detector.PHthaReco;
  PHbisReco = Detector.PHbisReco;
  PHcalReco = Detector.PHcalReco;
  PHmuoReco = Detector.PHmuoReco;

  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<7;i++){
    MeanNpe[i]=Detector.MeanNpe[i];
    MeanAtten[i] = Detector.MeanAtten[i];
    MeanSfactor[i] = Detector.MeanSfactor[i];
  }

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

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

  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 std::string[Ngamma];
  std::string* rold = Detector.Ogamma+Ngamma;
  std::string* 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;
  for(int i=0;i<4;i++){
    a[i] = rhs.a[i];
    z[i] = rhs.z[i];
    nn[i] = rhs.nn[i];
  }

  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;

  PHelectron = rhs.PHelectron;
  PHpositron = rhs.PHpositron;
  PHneutron = rhs.PHneutron;
  PHthallium = rhs.PHthallium;
  PHbismuth = rhs.PHbismuth;
  PHcalifornium = rhs.PHcalifornium;
  PHmuon = rhs.PHmuon;

  PHeleReco = rhs.PHeleReco;
  PHposReco = rhs.PHposReco;
  PHneuReco = rhs.PHneuReco;
  PHthaReco = rhs.PHthaReco;
  PHbisReco = rhs.PHbisReco;
  PHcalReco = rhs.PHcalReco;
  PHmuoReco = rhs.PHmuoReco;

  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<7;i++){
    MeanNpe[i] = rhs.MeanNpe[i];
    MeanAtten[i] = rhs.MeanAtten[i];
    MeanSfactor[i] = rhs.MeanSfactor[i];
  }

  Nneutron = rhs.Nneutron;

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

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

  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 std::string[Ngamma];
  std::string* rold = rhs.Ogamma+Ngamma;
  std::string* 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& nele,double* PHele,double* Tele,
                                 int& npos,double* PHpos,double* Tpos,
                                 int& nneu,double* PHneu,double* Tneu,
                                 int& ntha,double* PHtha,double* Ttha,
                                 int& nbis,double* PHbis,double* Tbis,
                                 int& ncal,double* PHcal,double* Tcal,
                                 int& nmuo,double* PHmuo,double* Tmuo){
  
  *(XYZpmt+0) = *(Xpmt+npmt);
  *(XYZpmt+1) = *(Ypmt+npmt);
  *(XYZpmt+2) = *(Zpmt+npmt);

  nele=0;
  npos=0;
  nneu=0;
  ntha=0;
  nbis=0;
  ncal=0;
  nmuo=0;

  for(PMT* PMTptr=PMTdata;PMTptr<PMTdata+PMTlen;PMTptr++){
    if(npmt==(*PMTptr).npmt){
      if((*PMTptr).Parent=="E"){
        *(PHele+nele) = (*PMTptr).PHsmr;
        *(Tele+nele) = (*PMTptr).TimeSmr;
        nele++;
      }
      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++;
      }
      if((*PMTptr).Parent=="Tl"){
        *(PHtha+ntha) = (*PMTptr).PHsmr;
        *(Ttha+ntha) = (*PMTptr).TimeSmr;
        ntha++;
      }
      if((*PMTptr).Parent=="Bi"){
        *(PHbis+nbis) = (*PMTptr).PHsmr;
        *(Tbis+nbis) = (*PMTptr).TimeSmr;
        nbis++;
      }
      if((*PMTptr).Parent=="Cf"){
        *(PHcal+ncal) = (*PMTptr).PHsmr;
        *(Tcal+ncal) = (*PMTptr).TimeSmr;
        ncal++;
      }
      if((*PMTptr).Parent=="muon"){
        *(PHmuo+nmuo) = (*PMTptr).PHsmr;
        *(Tmuo+nmuo) = (*PMTptr).TimeSmr;
        nmuo++;
      }
    } 
  }
  return;
}
void ReactorDetector::GenerateSecondaries(ReactorEvent& Event){

  //cout << "GenerateSecondaries(Event)" << endl;
  //
  // Displace the positron vertex by its range along the positron direction.
  //
  int Npositron = Event.GetNumPositron();
  //cout << Npositron << " positron(s):" << endl;

  double Rpos[Event.Rdim*Npositron];

  // loop over the positrons
  for(int n=0;n<Npositron;n++){

    double kepos = Event.GetPositronKE(n);
    double tpos = 0.;
    double rangepos = 0.;
    PositronRange(kepos,rangepos,tpos);
    tpos+=Event.GetPositronTime(n);
    //
    double cospos = Event.GetPositronCosTheta(n);
    double phipos = Event.GetPositronPhi(n);

    double rint[3];
    for(int i=0;i<3;i++){
      rint[i] = *(Event.GetPositronVertexXYZ()+(i+n*Event.Rdim));
    }
    //
    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]);
    //
    // set positron vertex in Event class
    //
    Rpos[0+n*Event.Rdim] = tpos;
    for(int i=1;i<4;i++)Rpos[i+n*Event.Rdim]=rpos[i-1];
    Rpos[4+n*Event.Rdim] = rrpos;
    Rpos[5+n*Event.Rdim] = rpos[2]/rrpos;
    Rpos[6+n*Event.Rdim] = atan2(rpos[1],rpos[0]);
    if(Rpos[6+n*Event.Rdim]<0)Rpos[6+n*Event.Rdim]+=2*RC.pi;
    Event.SetPositronVertex(n,Rpos);
  }
  //  for(int n=0;n<Event.Rdim*Npositron;n++){
  //    cout << " Rpos[" << n << "] " << Rpos[n] << endl;
  //  }
  // done with positrons

  int Nneutron = Event.GetNumNeutron();
  //cout << Nneutron << " neutron(s):" << endl;

  delete[] Absorber;
  delete[] NeutronSteps;

  Absorber = new double[Nneutron];
  NeutronSteps = new double[Nneutron];

  double Rneu[Event.Rdim*Nneutron];

  // loop over the neutrons
  for(int n=0;n<Nneutron;n++){
    //
    double rneu[3];
    double tneu = 0.;
    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 proportional to sqrt(path)
      //
      double rint[3];
      for(int i=0;i<3;i++){
        rint[i] = *(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
      }
      //
      double rrint = sqrt(rint[0]*rint[0]+rint[1]*rint[1]+rint[2]*rint[2]);
      double path = 0.;
      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+n*Event.Pdim));
        //cout << "pneu[" << i << "] " << pneu[i] << endl;
      }
      for(int i=0;i<3;i++){
        rneu[i]=*(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
        //cout << "rneu[" << i << "] " << rneu[i] << endl;
      }
      double keneu = Event.GetNeutronKE(n);
      NeutronSteps[n]=0.;
      Absorber[n]=0.;
      NeutronPropagator(R0,R2,keneu,Absorber[n],NeutronSteps[n],
                        pneu,rneu,tneu,Temperature,a,z,nn);
      //cout << " Absorber[" << n << "] " << Absorber[n] << endl;
      //cout << " NeutronSteps[" << n << "] " << NeutronSteps[n] << endl;
      tneu+=Event.GetNeutronTime(n);
      //cout << " tneu " << tneu << endl;
    }
    //
    // set neutron vertex in Event class
    //
    Rneu[0+n*Event.Rdim] = tneu;
    for(int i=1;i<4;i++)Rneu[i+n*Event.Rdim]=rneu[i-1];
    double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
    Rneu[4+n*Event.Rdim] = rrneu;
    Rneu[5+n*Event.Rdim] = rneu[2]/rrneu;
    Rneu[6+n*Event.Rdim] = atan2(rneu[1],rneu[0]);
    if(Rneu[6+n*Event.Rdim]<0)Rneu[6+n*Event.Rdim]+=2*RC.pi;
    Event.SetNeutronVertex(n,Rneu);
  }
  //  for(int n=0;n<Event.Rdim*Nneutron;n++){
  //    cout << " Rneu[" << n << "] " << Rneu[n] << endl;
  //  }
  // done with neutrons
}

void ReactorDetector::GenerateResponse(ReactorEvent& Event){

  //cout << "GenerateResponse(Event)" << endl;
  //
  PHelectron = 0;
  PHpositron = 0;
  PHneutron = 0;
  PHthallium = 0;
  PHbismuth = 0;
  PHcalifornium = 0;
  PHmuon = 0;

  PHeleReco = 0;
  PHposReco = 0;
  PHneuReco = 0;
  PHthaReco = 0;
  PHbisReco = 0;
  PHcalReco = 0;
  PHmuoReco = 0;

  for(int i=0;i<7;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;

  //cout << Ngamma << " gammas" << endl;
  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++;
      //for (int l=0;l<3;l++) cout << "  rpmt[" << l << "] " << rpmt[l] << endl;
      double r[3];
      double dr[3];
      double tr[3];
      //      
      double sumele[6]={0,0,0,0,0,0};
      double sumpos[6]={0,0,0,0,0,0};
      double sumneu[6]={0,0,0,0,0,0};
      double sumtha[6]={0,0,0,0,0,0};
      double sumbis[6]={0,0,0,0,0,0};
      double sumcal[6]={0,0,0,0,0,0};
      double summuo[6]={0,0,0,0,0,0};
      //
      double* rgptr=Rgamma;
      double* egptr=Egamma;
      double* tgptr=Tgamma;
      std::string* ogptr=Ogamma;
      //
      int nelegt0=0;
      int nposgt0=0;
      int nneugt0=0;
      int nthagt0=0;
      int nbisgt0=0;
      int ncalgt0=0;
      int nmuogt0=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];
        //cout << i << " " << j << " " << k << "  r " << sqrt(rdotr) << endl;
        //
        // nothing to do if r>R1: cycle through useless info
        //
        if(sqrt(rdotr)>R1){
          double dummy = *egptr++;
          dummy = *tgptr++;
          std::string sdummy = *ogptr++;
          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]);
        //cout << "  ddr " << ddr;
        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];
        //cout << "  tdotr " << tdotr;
        //
        // 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.
        //
        //cout << "  " << inter << endl;
        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;
        }
        //cout << "       sGd " << sGd << " sSc " << sSc << endl;
        //
        //  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++;
        //cout << "phgen " << phgen << endl;
        //
        float xpe = (phgen*PhotonsPerMeV);
        xpe *= (atten0*sfactor);
        xpe *= (PMTcoverage/Npmt*PMTqe);
        //cout << "xpe " << xpe << endl;
        //
        int npe;
        int error;
        rnpssn_(xpe,npe,error);
        //cout << "npe " << npe << endl;
        //
        double phsmr = npe/PMTqe/PMTcoverage/PhotonsPerMeV;   
        //cout << "phsmr " << phsmr << endl;
        //        
        double tgen = *tgptr++;
        //cout << "tgen " << tgen << endl;
        double tsmr = refracGd*sGd/RC.c+refracSc*sSc/RC.c + tgen; 
        //cout << "tsmr " << tsmr << endl;
        //
        std::string type = *ogptr++;
        //cout << "type " << type << endl;
        if(type == "E"){
          nelegt0++;
          sumele[0]+=phsmr;
          sumele[1]+=phgen;
          sumele[2]+=phgen*atten0*sfactor;
          sumele[3]+=npe;
          sumele[4]+=sfactor;
          sumele[5]+=atten;
          //
          PHelectron+=phsmr;
        }
        else if(type == "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 if(type == "N"){
          nneugt0++;
          sumneu[0]+=phsmr;
          sumneu[1]+=phgen;
          sumneu[2]+=phgen*atten0*sfactor;
          sumneu[3]+=npe;
          sumneu[4]+=sfactor;
          sumneu[5]+=atten;
          //      
          PHneutron+=phsmr;
        }
        else if(type == "Tl"){
          nthagt0++;
          sumtha[0]+=phsmr;
          sumtha[1]+=phgen;
          sumtha[2]+=phgen*atten0*sfactor;
          sumtha[3]+=npe;
          sumtha[4]+=sfactor;
          sumtha[5]+=atten;
          //      
          PHthallium+=phsmr;
        }
        else if(type == "Bi"){
          nbisgt0++;
          sumbis[0]+=phsmr;
          sumbis[1]+=phgen;
          sumbis[2]+=phgen*atten0*sfactor;
          sumbis[3]+=npe;
          sumbis[4]+=sfactor;
          sumbis[5]+=atten;
          //      
          PHbismuth+=phsmr;
        }
        else if(type == "Cf"){
          ncalgt0++;
          sumcal[0]+=phsmr;
          sumcal[1]+=phgen;
          sumcal[2]+=phgen*atten0*sfactor;
          sumcal[3]+=npe;
          sumcal[4]+=sfactor;
          sumcal[5]+=atten;
          //      
          PHcalifornium+=phsmr;
        }
        else if(type == "muon"){
          nmuogt0++;
          summuo[0]+=phsmr;
          summuo[1]+=phgen;
          summuo[2]+=phgen*atten0*sfactor;
          summuo[3]+=npe;
          summuo[4]+=sfactor;
          summuo[5]+=atten;
          //      
          PHmuon+=phsmr;
        }
        else{
          cout << "  ERROR: GenerateResponse has unrecognized type " 
               << type << endl;
        }
        //        
        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 = type;
        
        if(PMTptr<PMTdata+PMTsize){
          *PMTptr++=PMTtmp;
          PMTlen++;
        }
        else{
          //
          // expand the data array first
          //
          PMTptr = ResizePMT();
          *PMTptr++=PMTtmp;
          PMTlen++;
        }
      }
      //cout << "PHpositron " << PHpositron << endl;
      //cout << "PHneutron  " << PHneutron << endl;
      //
      //  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;
      }
      if(sumtha[2]>0){
        PHthaReco+=sumtha[0]*sumtha[1]/sumtha[2];
        MeanNpe[2]+=sumtha[3]/Npmt;
        MeanSfactor[2]+=sumtha[4]/Npmt/nthagt0;
        MeanAtten[2]+=sumtha[5]/Npmt/nthagt0;
      }
      if(sumbis[2]>0){
        PHbisReco+=sumbis[0]*sumbis[1]/sumbis[2];
        MeanNpe[3]+=sumbis[3]/Npmt;
        MeanSfactor[3]+=sumbis[4]/Npmt/nbisgt0;
        MeanAtten[3]+=sumbis[5]/Npmt/nbisgt0;
      }
      if(sumcal[2]>0){
        PHcalReco+=sumcal[0]*sumcal[1]/sumcal[2];
        MeanNpe[4]+=sumcal[3]/Npmt;
        MeanSfactor[4]+=sumcal[4]/Npmt/ncalgt0;
        MeanAtten[4]+=sumcal[5]/Npmt/ncalgt0;
      }
      if(summuo[2]>0){
        PHmuoReco+=summuo[0]*summuo[1]/summuo[2];
        MeanNpe[5]+=summuo[3]/Npmt;
        MeanSfactor[5]+=summuo[4]/Npmt/nmuogt0;
        MeanAtten[5]+=summuo[5]/Npmt/nmuogt0;
      }
      if(sumele[2]>0){
        PHeleReco+=sumele[0]*sumele[1]/sumele[2];
        MeanNpe[6]+=sumele[3]/Npmt;
        MeanSfactor[6]+=sumele[4]/Npmt/nelegt0;
        MeanAtten[6]+=sumele[5]/Npmt/nelegt0;
      }
      //
      // increment PMT index
      //
      pmtindex++;
    }
  }
}

void ReactorDetector::GenerateGammas(ReactorEvent& Event){

  //cout << "GenerateGammas(Event)" << endl;
  //
  // first check for any gammas produced in the event itself
  int ngamma = Event.GetNumGamma();
  //cout << " " << ngamma << " gammas from the Event" << endl;
  //
  int maxgamma = 100+ngamma;
  double* egamma = new double[maxgamma];
  double* tgamma = new double[maxgamma];
  double* rgamma = new double[3*maxgamma];
  std::string* ogamma = new std::string[maxgamma];
  //
  double* egptr=egamma;
  double* rgptr=rgamma;
  double* tgptr=tgamma;
  std::string* ogptr=ogamma;
  //
  // event gammas
  //
  for(int n=0;n<ngamma;n++){

    *egptr++ = Event.GetGammaKE(n);
    //
    for(int i=0;i<3;i++) 
      *rgptr++ = *(Event.GetGammaVertexXYZ()+(i+n*Event.Rdim));
    //
    *tgptr++ = Event.GetGammaTime(n);
    //
    if(Event.GetGammaParentProcess(n) == 1){
      cout << "ERROR: gammas are not produced by inverse-beta decay" << endl;
    }
    else if(Event.GetGammaParentProcess(n) == 2){
      *ogptr++ = "Cf";
    }
    else if(Event.GetGammaParentProcess(n) == 3){
      *ogptr++ = "muon";
    }
    else{
      cout << "Process Id " << Event.GetGammaParentProcess(n) << " is unknown" 
           << endl;
    }
  }
  //
  // radioactive source gammas
  //
  if(Event.GetPmtRad() == "on"){

    // pick a PMT location for the decay
    int len = 2;
    float r[len];
    ranlux_(r,len);
    int pmtId = r[0]*Npmt;

    // pick decay material and isotope 
    std::string Stuff; int Isotope;

    // this assumes even fractions of decays from all isotopes
    //
    if(r[1] < 0.25){
      Stuff = "thallium"; Isotope = 208; // thallium 208
    }
    else if(0.25 <= r[1] && r[1] < 0.5){
      Stuff = "thallium"; Isotope = 210; // thallium 210
    }
    else if(0.5 <= r[1] && r[1] < 0.75){
      Stuff = "bismuth"; Isotope = 212; // bismuth 212
    }
    else{
      Stuff = "bismuth"; Isotope = 214; // bismuth 214
    }
    //cout << Stuff << "(" << Isotope << ")" << endl;
    //
    // Thallium
    //
    if(Stuff == "thallium"){
      MyTlSource Thallium(Isotope);

      // create the weight array
      float PMTWeight[Thallium.GetNumTlGammas()];

      // loop over the decay gammas
      for(int n=0;n<Thallium.GetNumTlGammas();n++){

        // get the energy for each Tl gamma
        float e = Thallium.GetTlGammaEnergy(n); 
        *egptr++ = e;

        // propogate the gamma to the scintillating region
        PMTWeight[n] = 1.;
        int attempts = 1;
        double r[3];
        double time;

        // pick a direction for each gamma
        int len=4;
        float rv[len];
        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;
        //for(int i=0;i<3;i++) cout << "t["<<i<<"] " << t[i] << " ";
        //cout << endl;
        /*
         * propogate the gamma through the oil buffer: this assumes 
         * that the gamma propogates entirely through oil, but the 
         * gammas which reach the scint will travel at the end of 
         * their path through scint; however, the mineral oil and 
         * scint have very similar H and C densities
         */
        double range;
        std::string material = "oil";
        range = GammaComptonRange(material,e);
        range*=log(1/rv[2]);
        //cout << "range " << range << endl;
        //
        r[0] = Xpmt[pmtId]+range*t[0];
        r[1] = Ypmt[pmtId]+range*t[1];
        r[2] = Zpmt[pmtId]+range*t[2];
        //
        double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
        if(sqrt(rdotr)<R1){
          pmtRadCount++;
          cout << "PMTrad Tl " << pmtRadCount << " radius " << sqrt(rdotr) 
               << ", energy " << e << endl;
        }
        //
        // set the position
        for(int i=0;i<3;i++) *rgptr++ = r[i];
        //
        // pick a time relative to the neutrino capture, up to 1000000 ns
        time = rv[3]*1000000.;
        //cout << "time " << time << endl;
        *tgptr++ = time;
        *ogptr++ = "Tl";
        //
        PMTWeight[n] = 1./float(attempts);
        //cout << "weight " << PMTWeight[n] << endl;

        ngamma++;
      }
    }
    //
    // Bismuth
    //
    else if(Stuff == "bismuth"){
      MyBiSource Bismuth(Isotope);

      // create the weight array
      float PMTWeight[Bismuth.GetNumBiGammas()];

      // loop over the decay gammas
      for(int n=0;n<Bismuth.GetNumBiGammas();n++){

        // get the energy for each Bi gamma
        float e = Bismuth.GetBiGammaEnergy(n); 
        *egptr++ = e;

        // propogate the gamma to the scintillating region
        PMTWeight[n] = 1.;
        int attempts = 1;
        double r[3];
        double time;

        // pick a direction for each gamma
        int len=4;
        float rv[len];
        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;
        //for(int i=0;i<3;i++) cout << "t["<<i<<"] " << t[i] << " ";
        //cout << endl;
        /*
         * propogate the gamma through the oil buffer: this assumes 
         * that the gamma propogates entirely through oil, but the 
         * gammas which reach the scint will travel at the end of 
         * their path through scint; however, the mineral oil and 
         * scint have very similar H and C densities
         */
        double range;
        std::string material = "oil";
        range = GammaComptonRange(material,e);
        range*=log(1/rv[2]);
        //cout << "range " << range << endl;
        //
        r[0] = Xpmt[pmtId]+range*t[0];
        r[1] = Ypmt[pmtId]+range*t[1];
        r[2] = Zpmt[pmtId]+range*t[2];
        //
        double rdotr = r[0]*r[0]+r[1]*r[1]+r[2]*r[2];
        if(sqrt(rdotr)<R1){
          pmtRadCount++;
          cout << "PMTrad Bi " << pmtRadCount << " radius " << sqrt(rdotr) 
               << ", energy " << e << endl;
        }
        //
        // set the position
        for(int i=0;i<3;i++) *rgptr++ = r[i];
        //
        // pick a time relative to the neutrino capture, up to 1000000 ns
        time = rv[3]*1000000.;
        //cout << "time " << time << endl;
        *tgptr++ = time;
        *ogptr++ = "Bi";
        //
        PMTWeight[n] = 1./float(attempts);
        //cout << "weight " << PMTWeight[n] << endl;

        ngamma++;
      }
    } 
  }
  //
  // electron gammas, needs improvement
  //
  int Nelectron = Event.GetNumElectron();
  //cout << Nelectron << " electron(s):" << endl;

  // loop over the electrons
  for(int n=0;n<Nelectron;n++){

    double rele[3];
    for(int i=0;i<3;i++){
      rele[i] = *(Event.GetElectronVertexXYZ()+(i+n*Event.Rdim));
      //cout << " rele[" << i << "] " << rele[i] << endl;
    }
    double tele = Event.GetElectronTime(n);
    //cout << " tele " << tele << endl;
    double keele = Event.GetElectronKE(n);
    //cout << " keele " << keele << endl;
    //
    // generate a "fake" gamma with E=electron KE generated at vertex.
    //
    for(int i=0;i<3;i++)*rgptr++=rele[i];
    *tgptr++ = tele;
    *egptr++ = keele;
    *ogptr++ = "E";
    ngamma++;
  } // done with electrons
  //cout << "  " << ngamma << " gammas, so far" << endl;
  //
  // positron gammas
  //
  int Npositron = Event.GetNumPositron();
  //cout << Npositron << " positron(s):" << endl;

  // loop over the positrons
  for(int n=0;n<Npositron;n++){

    double rpos[3];
    for(int i=0;i<3;i++){
      rpos[i] = *(Event.GetPositronVertexXYZ()+(i+n*Event.Rdim));
      //cout << " rpos[" << i << "] " << rpos[i] << endl;
    }
    double tpos = Event.GetPositronTime(n);
    //cout << " tpos " << tpos << endl;
    double kepos = Event.GetPositronKE(n);
    //cout << " kepos " << kepos << endl;
    //
    // 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++){
      std::string material = "scintillator";
      range = GammaComptonRange(material,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++;
    }
  } // done with positrons
  //cout << "  " << ngamma << " gammas, so far" << endl;
  //
  // Generate neutron gammas
  //
  Nneutron = Event.GetNumNeutron();
  //cout << Nneutron << " neutron(s):" << endl;

  // loop over the neutrons
  for(int n=0;n<Nneutron;n++){

    double rneu[3];
    for(int i=0;i<3;i++){
      rneu[i] = *(Event.GetNeutronVertexXYZ()+(i+n*Event.Rdim));
      //cout << " rneu[" << i << "] " << rneu[i] << endl;
    }
    double tneu = Event.GetNeutronTime(n);
    //cout << " tneu " << tneu << endl;
    //
    // pick H or Gd, then which Gd peak
    //
    double emax=RC.Hpeak;
    if(FastNeutronOption){
      //cout << "in FastNeutronOption" << endl;
      Absorber[n]=1;
      double rrneu = sqrt(rneu[0]*rneu[0]+rneu[1]*rneu[1]+rneu[2]*rneu[2]);
      if(rrneu<R0){//in Gd region
        int len=2;
        float rv[2];
        ranlux_(rv,len);
        if(rv[0]<=GdCaptureFraction){
          //
          // model Gd spectrum by Gd155 + Gd157
          //
          if(rv[1]<=RC.Gd155frac){
            emax = RC.Gd155peak;
            Absorber[n]=155;
          }
          else{
            emax = RC.Gd157peak;
            Absorber[n]=157;
          }
        }
      }
    }
    else{
      //cout << "Absorber[" << n << "] " << Absorber[n] << endl;
      if(Absorber[n]==1)emax=RC.Hpeak;
      if(Absorber[n]==12)emax=RC.Cpeak;
      if(Absorber[n]==155)emax=RC.Gd155peak;
      if(Absorber[n]==157)emax=RC.Gd157peak;
    }
    //
    // One photon from an H or C capture
    //
    if(emax==RC.Hpeak || emax==RC.Cpeak){
      //cout << " H or C gammas" << endl;
      int len=3;
      float rv[3];
      ranlux_(rv,len);
      std::string material = "scintillator";
      double range = GammaComptonRange(material,emax);
      range*=log(1/rv[0]);
      double cosg = -1+2*rv[1];
      double phig = 2*RC.pi*rv[2];
      double t[3];
      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++;
      //cout << "  " << ngamma << " gammas" << endl;
    }
    //
    // otherwise model Gd photons
    //
    else{
      //cout << " Gd gammas" << endl;
      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]);
        std::string material = "scintillator";
        double range = GammaComptonRange(material,e);
        range*=log(1/rv2[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++;
        //cout << "  " << ngamma << " gammas" << endl;
      }
    }
  } // done with neutrons
  //cout << "  " << ngamma << " gammas, total" << endl;
  //
  // 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 std::string[ngamma];

  double* ptr;

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

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

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

  std::string* qtr = Ogamma+ngamma;
  while(qtr>Ogamma)*--qtr=*--ogptr;

  /*
  cout << "  " << Ngamma << " gammas" << endl;
  for(int n=0;n<Ngamma;n++){
    for(int i=0;i<3;i++) cout << "  Rgamma[" << 3*n+i << "] " 
                              << Rgamma[3*n+i];
    cout << endl;
    cout << "  Egamma[" << n << "] " << Egamma[n] 
         << "  Tgamma[" << n << "] " << Tgamma[n]
         << "  Ogamma[" << n << "] " << Ogamma[n] << endl;
  }
  */

  // clean up
  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(std::string parent){
  int gamma=0;
  for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
    if (*optr==parent)gamma++;
  }
  return gamma;
}
//
// Fetch array of positron gammas
//
void ReactorDetector::GetGamma(std::string parent,double* xyz){

  double* Rptr=Rgamma;
  double* rptr=xyz;
  //
  for(std::string* optr=Ogamma;optr<Ogamma+Ngamma;optr++){
    if (*optr==parent){
      for(int i=0;i<3;i++)*rptr++ = *Rptr++;
    }
  }
}  
void ReactorDetector::GetGamma(std::string parent,double* e,
                               double* t,double* r,
                               double* c,double* f){

  //cout << "GetGamma(" << parent << ")" << endl;
  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(std::string* 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/sqrt(x*x+y*y+z*z);
      *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;
  //cout << "  Set Npmt " << Npmt << endl;

  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);
  SetNeutronPropStuff();
}
void ReactorDetector::SetNeutronPropStuff(){
  //
  a[0] = RC.A0;
  a[1] = RC.A1;
  a[2] = RC.A2;
  a[3] = RC.A3;
  //
  z[0] = RC.Z0;
  z[1] = RC.Z1;
  z[2] = RC.Z2;
  z[3] = RC.Z3;
  //
  nn[0]=RC.Navogadro*RC.f0*RC.density/RC.A0;
  nn[1]=RC.Navogadro*RC.f1*RC.density/RC.A1;
  nn[2]=RC.Navogadro*GdConcentration/100*RC.f2*RC.density/RC.A2;
  nn[3]=RC.Navogadro*GdConcentration/100*RC.f3*RC.density/RC.A3;
  //
}
double ReactorDetector::GetHitFrac(std::string 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(std::string 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(std::string 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;
}    

---