---

ReactorScint.cpp

#include "ReactorConstants.hh"
#include "ReactorDetector.hh"
#include "ReactorScint.hh"
#include "ReactorFortran.hh"
#include "ReactorXC.hh"
#include "DetectorDefaults.hh"
#include <cmath>
#include <iostream>
using namespace std;
ReactorScint::ReactorScint(){

  Points = 0;
  RxyzPoints = new double[4*Points];
  PxyzPoints = new double[4*Points];
  NextTime = RxyzPoints;
  NextXYZ = RxyzPoints+1;
  NextEnergy = PxyzPoints;
  TooSoftToTrack = 0.1;

}
//
//  This constructor puts all the light at the point and time specified by
//  rxyz.    
//
ReactorScint::ReactorScint(double e,double t,double* rxyz){

  RxyzPoints = new double[4];
  PxyzPoints = new double[4];
  //
  for(int i=0;i<3;i++){
    *(RxyzPoints+i+1) = *(rxyz+i);
    *(PxyzPoints+i+1) = 0;
  }
  *RxyzPoints=t;
  *PxyzPoints=e;
  //
  Points = 1;
  NextEnergy = PxyzPoints;
  NextTime = RxyzPoints;
  NextXYZ = RxyzPoints+1;
  //
}
//
// This constructor takes a gamma and generates an array of light outputs via
// Compton scatters
//
ReactorScint::ReactorScint(double e,double t,double* rxyz,double* pxyz,
                           double rzero, double rone, double rtwo){
  //
  // make a guess at max number of hits
  //
  Points = 25;
  RxyzPoints = new double[4*Points];
  PxyzPoints = new double[4*Points];
  NextTime = RxyzPoints;
  NextXYZ = RxyzPoints+1;
  NextEnergy = PxyzPoints;
  TooSoftToTrack = 0.1;

  int scatters=0;
  double k[4];
  double r[4];
  for(int i=0;i<3;i++){
    k[i+1]=pxyz[i];
    r[i+1]=rxyz[i];
  }
  r[0] = t;
  k[0] = e;
  //
  double egamma = e;
  double tgamma = t;
  //cout << "   initial egamma " << egamma << " tgamma " << tgamma << endl;
  double range = 0.;
  bool inside = true;
  double r1 = rone;
  double r2 = rtwo;
  double* Rptr=RxyzPoints;
  double* Pptr=PxyzPoints;
  //
  while(egamma>TooSoftToTrack && scatters<Points && inside==true){

    //cout << " scatter " << scatters << endl;

    double rsq = sqrt(r[1]*r[1] + r[2]*r[2] + r[3]*r[3]);
    //cout << "   radius " << rsq << endl;
    //
    // check that we're inside the detector
    //
    if(rsq > r2){
      inside = false;
      //cout << " radius " << rsq << " > " << r2 
      //   << ": outside steel vessel" << endl;
    }
    else{
      //
      // count the scatter
      //
      scatters++;
      //
      // set the material
      //
      std::string material = "scintillator";
      if(rsq > r1) material = "oil";
      //cout << "   material " << material << endl;
      //
      GammaCompton(k);
      //
      range = GammaComptonRange(material,egamma);
      float rv[1];
      int len=1;
      ranlux_(rv,len);
      range*=log(1/rv[0]);
      //cout << "   range " << range << endl;
      //
      for(int i=1;i<4;i++){
        *(Pptr+i) = k[i]; 
        r[i]+=range*k[i]/k[0];
        *(Rptr+i) = r[i];
      }
      //
      double egammalast = egamma;
      egamma = k[0];
      tgamma += range/c*REFRACSC;
      //cout << "   egamma " << egamma << " tgamma " << tgamma << endl;
      //
      //cout << "   egammalast-egamma " << egammalast-egamma << endl;
      *Pptr = egammalast-egamma;
      *Rptr = tgamma;
      //
      Pptr+=4;
      Rptr+=4;
    }
  }
  //
  // if still inside put remaining energy in the last point
  //
  if(inside){
    scatters++;
    //
    for(int i=1;i<4;i++){
      *(Pptr+i) = k[i]; 
      r[i]+=range*k[i]/k[0];
      *(Rptr+i) = r[i];
    }
    *Pptr++=egamma;
    *Rptr++=tgamma;
  }
  //
  // Trim the size
  //
  double* Ptemp = new double[4*scatters];
  double* Rtemp = new double[4*scatters];
  for(int i=0;i<4*scatters;i++){
    *(Ptemp+i) = *(PxyzPoints+i);
    *(Rtemp+i) = *(RxyzPoints+i);
  }
  Points=scatters;
  delete[] PxyzPoints;
  delete[] RxyzPoints;
  RxyzPoints = new double[4*Points];
  PxyzPoints = new double[4*Points];
  for(int i=0;i<4*scatters;i++){
    *(PxyzPoints+i)=*(Ptemp+i);
    *(RxyzPoints+i)=*(Rtemp+i);
  }
  NextEnergy = PxyzPoints;
  NextTime = RxyzPoints;
  NextXYZ = RxyzPoints+1;
  delete[] Ptemp;
  delete[] Rtemp;

  /*
  cout << " " << scatters << " scatters:" << endl;
  for(int i = 0;i<scatters*4;i++)
    cout << "  PxyzPoints[" << i << "] " << PxyzPoints[i] << endl;
  for(int i = 0;i<scatters*4;i++)
    cout << "   RxyzPoints[" << i << "] " << RxyzPoints[i] << endl;
  */

}

ReactorScint::~ReactorScint(){

  delete[] RxyzPoints;
  delete[] PxyzPoints;

}

ReactorScint::ReactorScint(const ReactorScint& G){

  Points = G.Points;
  RxyzPoints = new double[4*Points];
  double* pold = G.RxyzPoints+4*Points;
  double* pnew = RxyzPoints+4*Points;
  while(pnew>RxyzPoints)*--pnew = *--pold;

  PxyzPoints = new double[4*Points];
  pold = G.PxyzPoints+4*Points;
  pnew = PxyzPoints+4*Points;
  while(pnew>PxyzPoints)*--pnew = *--pold;

  NextXYZ = RxyzPoints+1;
  NextTime = RxyzPoints;
  NextEnergy = PxyzPoints;

}


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

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

  Points = rhs.Points;
  delete[] RxyzPoints;
  RxyzPoints = new double[4*Points];
  double* pold = rhs.RxyzPoints+4*Points;
  double* pnew = RxyzPoints+4*Points;
  while(pnew>RxyzPoints)*--pnew = *--pold;

  delete[] PxyzPoints;
  PxyzPoints = new double[4*Points];
  pold = rhs.PxyzPoints+4*Points;
  pnew = PxyzPoints+4*Points;
  while(pnew>PxyzPoints)*--pnew = *--pold;

  NextXYZ = RxyzPoints+1;
  NextTime = RxyzPoints;
  NextEnergy = PxyzPoints;

  return *this;
}
double ReactorScint::GetNextEnergy(){

  if(NextEnergy>=PxyzPoints+4*Points){
    return -1.;
  }
  NextEnergy+=4;
  return *(NextEnergy-4);
}
double ReactorScint::GetNextTime(){

  if(NextTime>=RxyzPoints+4*Points){
    return -1.;
  }
  NextTime+=4;
  return *(NextTime-4);
}
double* ReactorScint::GetNextXYZ(){

  if(NextXYZ>=RxyzPoints+4*Points){
    return 0;
  }
  NextXYZ++;
  return NextXYZ-1;
}

void ReactorScint::ResetNext(){
  NextTime = RxyzPoints;
  NextXYZ = NextTime+1;
  NextEnergy = PxyzPoints;
}

---