---

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"
//
// 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];

  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;

  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];
  }

  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 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];
  }

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

  NeutronSteps = new double[Nneutron];
  delete[] NeutronSteps;
  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 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){

  //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);
    //
    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 protportional 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);
      //cout << " Absorber[" << n << "] " << Absorber[n] << endl;
      //cout << " NeutronSteps[" << n << "] " << NeutronSteps[n] << 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;
  //
  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;

  //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 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];
        //cout << i << " " << j << " " << k << "  rdotr " << rdotr << endl;
        //
        // nothing to do if r>R1: cycle through useless info
        //
        if(sqrt(rdotr)>R1){
          double dummy = *egptr++;
          dummy = *tgptr++;
          char 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;
        //
        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++;
        }
      }
      //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;
      }
      //
      // increment PMT index
      //
      pmtindex++;
    }
  }
}

void ReactorDetector::GenerateGammas(ReactorEvent& Event){

  //cout << "GenerateGammas(Event)" << endl;
  //
  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;
  //
  // radioactive source gammas
  //
  if(Event.GetPmtRad() == "on"){

    int Isotope;
    //
    // thallium 208
    Isotope = 208;
    MyTlSource Tl208(Isotope);

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

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

      // get the energy for each Tl gamma
      float e = Tl208.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 " << pmtRadCount << " radius " << sqrt(rdotr) << endl;
      }

      // set the position
      for(int i=0;i<3;i++) *rgptr++ = r[i];

      // pick a time relative to the neutrino capture, up to 10000 ns
      time = rv[3]*10000.;
      //cout << "time " << time << endl;
      *tgptr++ = time;

      *ogptr++ = 'T';

      PMTWeight[n] = 1./float(attempts);
      //cout << "weight " << PMTWeight[n] << endl;

      ngamma++;
    }
    //
    // thallium 210
    /*
      Isotope = 210;
      MyTlSource Tl210(Isotope);

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

      // get the energy for each Tl gamma
      for(int n=0;n<Tl210.GetNumTlGammas();n++){

      *rgptr++ = Xpmt[pmtId];
      *rgptr++ = Ypmt[pmtId];
      *rgptr++ = Zpmt[pmtId];
      *tgptr++ = 0; 
      *egptr++ = Tl210.GetTlGammaEnergy(n); 
      *ogptr++ = 'Tl';

      ngamma++;
      }
    */
  }
  //
  // 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 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;

  /*
  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;
  }
  */

  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;
}    

---