---

ReactorScint.cpp

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

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

}
//
//  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 photon and generates an array of light outputs via
// Compton scatters
//
ReactorScint::ReactorScint(ReactorConstants& RC,
                           double t,double* rxyz,double* pxyz){

  //
  // 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;
  
  int scatters=0;
  double k[4];
  double r[4];
  for(int i=1;i<4;i++){
    k[i]=pxyz[i];
    r[i]=rxyz[i-1];
  }
  r[0] = t;
  k[0] = pxyz[0];
  scatters++;
  //
  double egamma = pxyz[0];
  double* Rptr=RxyzPoints;
  double* Pptr=PxyzPoints;
  while(egamma>TooSoftToTrack&&scatters<Points){
    GammaCompton(RC,k);
    //
    double range = GammaComptonRange("scintillator",egamma);
    float rv[1];
    int len=1;
    ranlux_(rv,len);
    range*=log(1/rv[0]);
    //
    scatters++;
    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];
    *Pptr = egammalast-egamma;
    *Rptr = *(Rptr-4)+range/RC.c*RC.refracSc;
    //
    Pptr+=4;
    Rptr+=4;
  }
  //
  // put remaining energy in the last point
  //
  *Pptr=egamma;
  *Rptr=-1;
  for(int i = 1;i<4;i++)*Pptr++=0;
  for(int i = 1;i<4;i++)*Rptr++=0;
  //
  // 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;
}

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+=4;
  return NextXYZ-4;
}

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

---