// ReactorGen.cc
// Contact person: Nuno Barros <nfbarros@hep.upenn.edu>
// See ReactorGen.hh for more details
//———————————————————————//

#include <RAT/ReactorGen.hh>

#include <RAT/GLG4PosGen.hh>
#include <RAT/GLG4TimeGen.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/Factory.hh>
#include <RAT/DB.hh>
#include <RAT/EventTime.hh>
#include <RAT/Log.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>
#include <RAT/UnitsUtil.hh>
#include <RAT/VertexGen_IBD.hh>
#include <RAT/StringUtil.hh>
#include <RAT/LLAUtil.hh>

// -- ROOT includes
#include <TMath.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TF1.h>
#include <TGraph.h>

// -- Geant4 includes
#include <G4Event.hh>
#include <G4Track.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4ThreeVector.hh>
#include <G4TransportationManager.hh>
#include <G4Navigator.hh>

#include <stdexcept>
#include <string>
#include <cstdlib>
#include <algorithm>

using std::string;
using std::vector;

// using CLHEP units
using CLHEP::cm;
using CLHEP::mm;
using CLHEP::km;
using CLHEP::keV;
using CLHEP::MeV;
using CLHEP::ns;
using CLHEP::s;
using CLHEP::day;

namespace RAT {


/**
 * SNO Longitude = -1*(81+12./60+5./3600.0); = -81.2014 degrees
 * SNO Latitude  = (46+28./60+31.0/3600.0); = 46.4753 degrees
 * SNO Altitude  = (307.0 - 2073.0); //elevation of SNOLAB(Lively) - SNO depth (in meters)
 **/
const G4ThreeVector ReactorGen::SNO_LLA_coord_ = G4ThreeVector(-81.2014, 46.4753,-1766.0);
const G4ThreeVector ReactorGen::SNO_ECEF_coord_ = G4ThreeVector(672.87,-4347.18,4600.51);

ReactorSpectrum::ReactorSpectrum(G4String name,
    SpectrumType stype,
    double emin,
    double emax,
    std::vector<double> spec_e,
    std::vector<double> spec_f,
    std::vector<std::string> isotopes,
    std::vector<double> composition) {
  _name = name;
  _e_min = emin;
  _e_max = emax;
  _spec_type = stype;
  _spectrum = NULL;

  switch(_spec_type) {
  case parameterized:
    // Spectrum is built from the contributions of several isotopes
    BuildSpectrumFromIsotopes(isotopes,composition);
    break;
  case shape:
    _spectrum = BuildSpectrumFromArrays(spec_e,spec_f);
    _spectrum->SetName(_name.c_str());
    _spectrum_norm = _spectrum->Integral("width");
    _spectrum->Scale(1./_spectrum_norm);
    break;
  case unknown:
  default:
    Log::Die("Antineutrino spectrum type not recognized.",5000);
    break;
  }

}

ReactorSpectrum::~ReactorSpectrum() {
  delete _spectrum;
}

void ReactorSpectrum::BuildSpectrumFromIsotopes(std::vector<std::string> &isotopes, std::vector<double> &composition) {
  // -- Parameterized spectrum. The final spectrum requires loading the individual isotopes and summing the spectra up
  DB* db = DB::Get();
  DBLinkPtr linkdb;
  std::string stype;
  std::vector<double> param1;
  std::vector<double> param2;

  std::vector<double> e_fission;
  std::vector<TH1D*> h_spec;
  TF1 * iso_fcn = NULL;
  TH1D *iso_spec;
  const int nbins=1000;
  std::string iso_formula;
  _spectrum = new TH1D(Form("%s_tot_nu_spec",_name.c_str()),"",nbins,_e_min,_e_max);

  double e_total_k = 0.0;
  for (size_t iso_idx = 0; iso_idx < isotopes.size(); ++iso_idx) {
    // Get the isotope table
    linkdb = db->GetLink("REACTOR_ISOTOPE_INFO",isotopes.at(iso_idx));
    // Load the information for this isotope
    stype =linkdb->GetS("spec_type");
    if ( stype == std::string("polynomial")) {
      // grab the polynomial parameters and construct it
      param1 = linkdb->GetDArray("poly_coeff");
      // Construct a TF1 based on these coefficients
      iso_formula = "exp(";
      for (unsigned int i = 0; i < param1.size();++i) {
        iso_formula += Form("[%d]",i);
        for (unsigned int j = 0 ; j < i; ++j) {
          iso_formula += "*x";
        }
        if (i < param1.size()-1) {
          iso_formula += "+";
        }
      }
      iso_formula += ")";
#ifdef REACTORDEV_DEBUG
      debug << "Formula for isotope " << isotopes.at(iso_idx) << " : [" << iso_formula << "]" << newline;
#endif
      iso_fcn = new TF1(Form("%s_fcn",isotopes.at(iso_idx).c_str()),iso_formula.c_str(),_e_min,_e_max);

      for(unsigned int i=0;i<param1.size();++i)
        iso_fcn->SetParameter(i,param1.at(i));

      // Convert a function into a histogram with 100 bins from 1.8 to 9.0
      // NOTE: This is set fixed so that all the histograms for each isotope can safely be added together
      iso_spec = new TH1D(Form("%s_spec",isotopes.at(iso_idx).c_str()),"",nbins,_e_min,_e_max);

      iso_spec->Eval(iso_fcn);

      // -- push back this spectrum into the vector
      h_spec.push_back(iso_spec);
      e_fission.push_back(linkdb->GetD("e_per_fission"));
      e_total_k += composition.at(iso_idx)*e_fission.at(e_fission.size()-1);
#ifdef REACTORDEV_DEBUG
      detail<<"Isotope "<<isotopes.at(iso_idx).c_str()<<" Integral="<< iso_spec->Integral("width")<<" antinu / Erel="<<e_fission.at(e_fission.size()-1)<<" MeV"<<newline;
#endif

      delete iso_fcn;
    } else if (stype == std::string("shape")) {
      // -- If the isotope spectrum is still in terms of shape, it is assumed that it is
      // a normalized one and each value is the content of a bin
      // Build an interpolator and then fill a histogram with it
      param1 = linkdb->GetDArray("spec_e");
      param2 = linkdb->GetDArray("spec_flux");

      // Build a TGraph from these points
      TGraph *gtmp = new TGraph(param1.size(),&param1[0],&param2[0]);
      iso_spec = BuildSpectrumFromArrays(param1,param2);
      iso_spec->SetName(Form("%s_spec",isotopes.at(iso_idx).c_str()));

      h_spec.push_back(iso_spec);
      e_fission.push_back(linkdb->GetD("e_per_fission"));
      e_total_k += composition.at(iso_idx)*e_fission.at(e_fission.size()-1);
      delete gtmp;
    } else {
      Log::Die("Found unknown spectrum type in isotope" + isotopes.at(iso_idx));
    }
  }

  for (size_t iso_idx = 0; iso_idx < isotopes.size(); ++iso_idx) {
#ifdef REACTORDEV_DEBUG
    detail << "[ReactorSpectrum]::BuildSpectrumFromIsotopes : Adding isotope " << iso_idx <<  " name " << isotopes.at(iso_idx) << newline;
#endif
    // Build up the final spectrum by adding the newly produced histogram -- should be #anti-neutrinos of each energy / MeV release
    _spectrum->Add(h_spec.at(iso_idx),composition.at(iso_idx)/e_total_k);
  }


  // By now the final spectrum should be built.
  _spectrum_norm = _spectrum->Integral("width");
  //scale to one before convolving with cross-section
  _spectrum->Scale(1./_spectrum_norm);

  // Clear the temporary histograms in h_spec
  for (size_t i = 0; i < h_spec.size(); ++i) {
    delete h_spec[i];
  }
  h_spec.clear();

#ifdef REACTORDEV_DEBUG
  detail << "[ReactorSpectrum]::BuildSpectrumFromIsotopes : Spectrum built: " << _spectrum->GetName() <<  " with " <<  _spectrum_norm << " neutrinos/MeV release --- integral is now " << _spectrum->Integral("width")<<newline;
  detail << "[ReactorSpectrum]::BuildSpectrumFromIsotopes : Hist address " << _spectrum << newline;
#endif
}


TH1D * ReactorSpectrum::BuildSpectrumFromArrays(std::vector<double> spec_e, std::vector<double> spec_f) {
  // Build a TGraph from these points
  const int nbins=1000;
  TGraph *gtmp = new TGraph(spec_e.size(),&spec_e[0],&spec_f[0]);
  TH1D *iso_spec = new TH1D("tmp_spec","",nbins,_e_min,_e_max);
  for (int bin = 1; bin < nbins; ++bin) {
    double val = gtmp->Eval(iso_spec->GetXaxis()->GetBinCenter(bin));
    iso_spec->SetBinContent(bin,val);
  }
  return iso_spec;
}


ReactorGen::ReactorGen() : GLG4Gen("reactor"), fStateStr(""),
    fCurrentVertexGen(0),
    fTimeGenName("poisson"),fPosGenName("fill"),fVertexGenName("ibd"),
    fVertexStateString(""),
    fReactorListName(""),fPowerPDF(0),
    fNuFlavor("nue"),fUseReactorPosition(true),
    fNuDir(0.0,-1.0,0.0),fFluxScale(1.0),fEvtRate(0.0),volumeInfoLoaded(false),
    fRatePerTarget(0.0), fSelectedReactor("none"),fSelectedCore(0)
{
#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]:: In constructor." << newline;
#endif
  fEvTime.SetStart();
  SetPosGenerator(fPosGenName,true);
  SetTimeGenerator(fTimeGenName,true);
}


ReactorGen::~ReactorGen()
{

#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]:: In destructor." << newline;
#endif

}

void ReactorGen::GenerateEvent(G4Event *event)
{

#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]:: In GenerateEvent. " << event << newline;
#endif

  // To generate an event we need to generate two sampling distributions.
  // Direction: Built from the ephemeris data. Build a live time distribution for the required time span and then sample the direction from there
  // Energy : Use the neutrino spectrum and sample the energy from there

  G4ThreeVector pos;
  fPosGen->GeneratePosition(pos);
  G4double t0 = NextTime();

  // The pointer to the generator that should be used is set here
  GenerateNuDirection();
#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::GenerateEvent : Incoming neutrino direction set to  [ " << fNuDir.x() << " " << fNuDir.y() << " " << fNuDir.z() << " ]." << newline;
  //  detail << "[ReactorGen]::GenerateEvent : Event time (before time generator offset): [" <<  fEvTime.GetStrTime() << "] date : [" << fEvTime.GetStrDate() << "] UTtime " << fUTtime<< "."<< newline;
  detail << "[ReactorGen]::GenerateEvent : Event position ( " << pos.x()/cm << ","<< pos.y()/cm << ","<< pos.z()/cm << ") cm " << newline;
#endif


  fCurrentVertexGen->GeneratePrimaryVertex(event, pos, t0, fNuDir);
#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::GenerateEvent : Event generated. Updating meta info with reactor information." << newline;
#endif

  G4String meta_info = Form("%s %d",fSelectedReactor.data(),fSelectedCore);
  // as a bonus, add the reactor/core information to the meta info
  ParentPrimaryVertexInformation * parent_vtx = dynamic_cast<ParentPrimaryVertexInformation*>(event->GetPrimaryVertex(0)->GetUserInformation());

  // Need to add the name of the reactor and core that generated this neutrino
  parent_vtx->SetVertexParentMetaInfo(meta_info);
  event->GetPrimaryVertex(0)->SetUserInformation(parent_vtx);

}

void ReactorGen::ResetTime(double offset)
{

  fNextTime = fTimeGen->GenerateEventTime() + offset;
  fEvTime.IncrementUT(fNextTime);
#ifdef REACTORDEV_DEBUG
  //fUTtime += fNextTime;

  debug << "[ReactorGen]::ResetTime : Next event time : " << fNextTime/ns <<" ns. " <<newline;
  debug << "                        Corresponding to :" << newline;
  debug << "                        Event time : [" <<  fEvTime.GetStrTime() << "] date : [" << fEvTime.GetStrDate() << "]."<< newline;
#endif
}


double ReactorGen::NumNeutrinosPSecCore(unsigned int ireactor,unsigned int icore) {

  // Conversion factor 1 MeV = 1.6e-13 Watt.s
  static const double MeVToWatt_s = 1.60217733e-13;

  Reactor::CoreData &core = fReactorData.at(ireactor).GetCore(icore);

  double spec_norm = 0.0;
  double n_nus = 0.0;
  double core_power = 0.0;

  // Grab the power of the reactor in MegaWatts
  spec_norm = fReactorSpectra[core.core_spectrum]->GetSpectrumNorm();
  // in nu/MeV
  // Now use the power to estimate the produced neutrino flux
  core_power = core.core_power*1e6/MeVToWatt_s; // core power in MeV/s
  //We use thermal power, and for now include duty_cycle, so it should be ok!
  //later treatment will have time evolution
  n_nus = spec_norm*core_power; // number of neutrinos produced per second at this reactor

#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::NumNeutrinosPSecCore : core_power [MeV/s] * spec_norm [anti-nu/MeV]? =" << core_power <<" * "<< spec_norm <<newline;
#endif




  return n_nus;
}

void ReactorGen::GenerateNuDirection() {

  // Generating the direction is not a trivial issue.
  // Based on the power of the different reactors one should decide which reactor/core
  // should the neutrino be generated from
  // Then, based on that, one directly has the direction

  // Direction should be calculated based on the position of the reactor that is being currently sampled.
  // If the direction is fixed do nothing.
  // If we are using the reactor position data pick up the reactor from the core PDF

  /// Have to find the latitude and longitude of the SNO+ detector

  if (fUseReactorPosition) {
    /// Let me explain what is being done here:
    ///
    /// TH2D::GetRandom2 samples a 2D coordinate from the existing contents, with the probability being scaled
    /// by the histogram contents themselves. The nice thing about this function is that empty
    /// bins are ignored and therefore a value will never be generated for a non-existing, or non-active core
    ///
    /// The trick is in the choice of the range of the histogram, so that the truncated number has precisely the
    /// correct probability
    /// in a bin has everything in the range of
    /// the truncated number

    double x,y;
    fPowerPDF->GetRandom2(x,y);
    //  The truncated part of the random number gives us the position

#ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::GenerateNuDirection : Random numbers generated: x " << x << " y " << y  << newline;
    detail << "[ReactorGen]::GenerateNuDirection : Using reactor " << static_cast<int>(x) << " core " << static_cast<int>(y)  << newline;

    detail << "[ReactorGen]::GenerateNuDirection : Using reactor " <<  fReactorData.at((int)x).GetName() << " core " << static_cast<int>(y)  << newline;
    detail << "[ReactorGen]::GenerateNuDirection : Latitude " <<  fReactorData.at((int)x).GetCore((int)y).core_LLA_coord.x() << " Longitude " << fReactorData.at((int)x).GetCore((int)y).core_LLA_coord.y()
               << " Altitude " << fReactorData.at((int)x).GetCore((int)y).core_LLA_coord.z()<< newline;

#endif


    fNuDir = fReactorData.at((int)x).GetCore((int)y).core_direction;
    fCurrentVertexGen = fVertexGenMap.at(fReactorData.at((int)x).GetCore((int)y).core_spectrum);
    fSelectedReactor = fReactorData.at((int)x).GetName();
    fSelectedCore = (int) y;

#ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::GenerateNuDirection : Time set." << newline;
    //    detail << "[ReactorGen]::DATE: " << days << " " << secs << " " << nsecs << newline;
    detail << "[ReactorGen]::DIRECTION: " << fNuDir.x() << " " << fNuDir.y() << " " << fNuDir.z() << newline;
#endif
  } else if (fReactorData.size() != 1) {
    // There is more than one vertex generator. Something went very wrong
    Log::Die("ReactorGen::GenerateNuDirection : Neutrino direction fixed but more than one reactor listed.",5000);
  } else {
    // We still want to decide on the core
    double x,y;
    fPowerPDF->GetRandom2(x,y);
    fSelectedCore = (int) y;
    fSelectedReactor = fReactorData.at(0).GetName();
    fCurrentVertexGen = fVertexGenMap.at(fReactorData.at(0).GetCore((int)y).core_spectrum);
  }
}


void ReactorGen::SetTimeGenerator(G4String generator, bool force) {
  // Only take action if we are using a different generator
  if ((generator != fTimeGenName) || force) {
    if (fTimeGen) delete fTimeGen;
    fTimeGen = RAT::GlobalFactory<GLG4TimeGen>::New(generator);
    fTimeGenName = generator;
  }
}

void ReactorGen::SetPosGenerator(G4String generator, bool force) {
  // Only take action if we are using a different generator
  if ((generator != fPosGenName) || force) {
    if (fPosGen) delete fPosGen;
    fPosGen = RAT::GlobalFactory<GLG4PosGen>::New(generator);
    // The state for the position generator is for the moment fixed to the
    // active region (scintillator)
    // Center of the detector is always a good choice
    G4String posState = "0.0 0.0 0.0";
    fPosGen->SetState(posState);
    fPosGenName = generator;
  }
}

void ReactorGen::SetNuDirection(G4String newValue) {
  // If this method is called then the incoming neutrino direction is
  // fixed and therefore we won't use the reactor data to determine
  // the direction of the neutrino

  // This repeats the code used by the vertex generator
  newValue = util_strip_default(newValue); // from GLG4StringUtil
  if (newValue.length() == 0) {
    return;
  }

  std::istringstream is(newValue.c_str());
  double x, y, z;
  is >> x >> y >> z;
  if (is.fail())
    return;

  if ( x == 0. && y == 0. && z == 0. )
    fNuDir.set(0., 0., 0.);
  else
    fNuDir = G4ThreeVector(x, y, z).unit();

  fUseReactorPosition = false;
}

void ReactorGen::CalculateTotalRate() {
  // This method calculates the total neutrino *interaction rate* for the selected generator options
  // This means that it takes into account the following:
  //
  // * Power of each active core: --> Estimate the total emitted flux of neutrinos
  // * Spectrum emitted by each core --> Estimate the neutrino energy spectra arriving to the detector, per core
  // * Distance of each active core --> Estimate total flux of neutrinos at the detector (make a point like estimate => 1/d^2)
  // * Type of interaction being used --> Weight the emitted spectra to convert into detected spectra
  //
  // OTHER NOTES:
  // - The convolution of the *detected* neutrino spectrum with the cross section yields the interaction rate per target.
  // - Product of the above with the number of targets yields the event rate of the detector


  // declare them as static since they are common to all instances of the class
  static G4double theVolume = 0.0;
  static G4double fNTarget = 0.0;

  if (!volumeInfoLoaded) {
    // Update the event rate based on this value
    // Use the total cross section, the total flux and the scale provided
    // Here we need to get the material properties and detector properties
    // The center of the detector is always inside the active volume

    G4ThreeVector sample_position;
    fPosGen->GeneratePosition(sample_position);
#ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::CalculateTotalRate : Sample position obtained in " << sample_position.x() << " "
        << sample_position.y() << " " << sample_position.z() << newline;
#endif

    G4Navigator* theNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
    G4VPhysicalVolume* thePhyVolume = theNavigator->LocateGlobalPointAndSetup(sample_position);

    if (thePhyVolume) {
      debug << "[ReactorGen]::CalculateTotalRate : Interactions will be generated in volume " << thePhyVolume->GetName() << newline;
    } else {
      Log::Die("[ReactorGen]::CalculateTotalRate : Couldn't find the physical volume to generate the interactions.");
    }

    // Get the logical volume associated with the active medium
    G4LogicalVolume* theLogVolume = thePhyVolume->GetLogicalVolume();
    G4Material *theMaterial = theLogVolume->GetMaterial();
    // Get the associated solid (may not be a sphere)
    G4VSolid *theSolid = theLogVolume->GetSolid();

    // Get the volume
    theVolume = theSolid->GetCubicVolume(); // G4 uses cubic millimeters

#ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::CalculateTotalRate : Volume estimated to be " << theVolume/CLHEP::m3 << " m3 " << newline;
#endif

    fNTarget = 0.0;
    // -- Here things get more interesting. Depending on the
    // vertex generator, the calculation comes different
    // Get the electron density
    /** This code remains here in preparation to add support to ES reactor events
     * FIXME: Finish implementing ES support
    if (fVertexGenName == G4String("es")) {
  #ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::CalculateTotalRate : Estimating # targets for ES" << newline;
  #endif

      // This makes sense for ES
      fNTarget = theMaterial->GetElectronDensity();
      warn << "[ReactorGen]::CalculateTotalRate : The vertex generator is set to ES. The reactor spectrum below 1.8 MeV is only approximate." << newline;

       } else if (fVertexGenName == G4String("ibd")) {
     **/
#ifdef REACTORDEV_DEBUG
    info << "[ReactorGen]::CalculateTotalRate : Estimating # targets for IBD" << newline;
#endif
    // Here we need to know how many Hydrogens (free protons) we have per unit volume (G4 uses mm3)
    // This makes sense to IBD
    const G4ElementVector * comp_elem = theMaterial->GetElementVector();
    G4int n_elems = theMaterial->GetNumberOfElements();
#ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::CalculateTotalRate : Material has " << n_elems << " elements." << newline;
#endif
    const G4double *n_elem_vol = theMaterial->GetVecNbOfAtomsPerVolume();

    for (int i = 0; i < n_elems; ++i) {
#ifdef REACTORDEV_DEBUG
      detail << "[ReactorGen]::CalculateTotalRate : Checking element " << comp_elem->at(i)->GetName() << newline;
      detail << "[ReactorGen]::CalculateTotalRate : # : " << n_elem_vol[i] << newline;
#endif
      if (comp_elem->at(i)->GetName() == G4String("Hydrogen")) {
        fNTarget = n_elem_vol[i];
        break;
      }
    }
#ifdef REACTORDEV_DEBUG
    detail << "[ReactorGen]::CalculateTotalRate : Estimated a number of atoms of  " << fNTarget << newline;
#endif

    if (fNTarget == 0.0) {
      Log::Die("[ReactorGen]::CalculateTotalRate:Couldn't estimate number of free protons in the simulation volume.");
    }
    /** FIXME: Finish implementing ES support
    } else {
      Log::Die("[ReactorGen]::CalculateTotalRate type of neutrino vertex generator ["+fVertexGenName + "]");
    }
     **/
    fNTarget *= theVolume;
#ifdef REACTORDEV_DEBUG
    info << "[ReactorGen]::CalculateTotalRate : Estimated total number of targets : " << fNTarget/1.e32 <<" x 10^32 protons"<<newline;
#endif

    volumeInfoLoaded = true;

  } // volumeInfoLoaded
#ifdef REACTORDEV_DEBUG
  else {
    info <<  "[ReactorGen]::CalculateTotalRate : Volume info already loaded." << newline;
  }
#endif
  // -- Now calculate the event rate


  fEvtRate = 0.0;
  // I would expect that the flux is necessary here...
  // flux is anti-neutrinos/s/cm^2
  // rate per target keeps track of number of interactions/s/target (but normalized to 1)
  double rate_per_target = 0.0, flux_detector = 0.0;


  //event_rate = flux_scale*Ntargets_total*rate_per_target
  // rate_per_target = integral of neutrino_spectrum*sigma
  /**
   * Normalize the spectrum to 1. Convolve it with the total cross section
   * rate_per_target = convolved_spectrum*total_flux*any_necessary_norm
   */

  // loop over all reactors
#ifdef REACTORDEV_DEBUG
  detail <<  "[ReactorGen]::CalculateTotalRate : Loading the reactor details." << newline;
#endif

  for (size_t i = 0; i < fReactorData.size(); ++i) {
#ifdef REACTORDEV_DEBUG
    detail <<  "[ReactorGen]::CalculateTotalRate : At reactor " << i << newline;
#endif
    Reactor &data = fReactorData.at(i);
    // loop over all cores
    for (size_t j = 0; j < data.GetNumCores(); ++j) {
#ifdef REACTORDEV_DEBUG
      detail <<  "[ReactorGen]::CalculateTotalRate : At reactor " << i << ", core " << j << newline;
#endif
      Reactor::CoreData &core = data.GetCore(j);
      // an individual core in a reactor can be turned off
      if (!core.is_active) {
#ifdef REACTORDEV_DEBUG
        detail <<  "[ReactorGen]::CalculateTotalRate : At reactor " << i << ", core " << j << " ...not active!" << newline;
#endif
        continue;
      }

#ifdef REACTORDEV_DEBUG
      info <<  "[ReactorGen]::CalculateTotalRate : Calculating rate per target for core type " << core.core_spectrum << newline;
      detail <<  "[ReactorGen]::CalculateTotalRate : beware of units, XSec in mm2, Flux will be in km2!"<<newline;
#endif
      // This gives the number of neutrinos/MeV release convoluted with the cross section in cm2;
      rate_per_target = fVertexGenMap.at(core.core_spectrum)->GetRatePerTarget();
#ifdef REACTORDEV_DEBUG
      detail <<  "[ReactorGen]::CalculateTotalRate : RatePerTarget = spectrum (norm=1) x cross-section = " << rate_per_target/CLHEP::cm2 << " cm2"<<newline;
      detail <<  "[ReactorGen]::CalculateTotalRate : Core Flux [nu/sec / 4piD2] = " << NumNeutrinosPSecCore(i,j) <<" / " <<(4*TMath::Pi()*pow(core.core_distance/CLHEP::cm,2.))<<newline;
#endif
      rate_per_target*=NumNeutrinosPSecCore(i,j);

      // This gives the number of neutrinos that actually pass by the detector //10 passes from km2 -> cm2, but cross-section is in mm2!
      flux_detector = rate_per_target/(4*TMath::Pi()*pow(core.core_distance,2.));

      // Have to account for how many hydrogens exist to scale things
      fEvtRate += flux_detector*fNTarget;
#ifdef REACTORDEV_DEBUG
      detail <<  "[ReactorGen]::CalculateTotalRate : flux at detector * Target = " << flux_detector <<" * "<< fNTarget << newline;
      info << " This core contribution is = "<< flux_detector*fNTarget <<" Hz in "<< fEvtRate<<" Hz total "<< newline;
#endif
    }
  }

  G4String newstate = util_to_string(fEvtRate*fFluxScale);
  // Pass this state also to the time generator
  fTimeGen->SetState(newstate);

  detail << "[ReactorGen]::CalculateTotalRate: Event rate updated to " << fEvtRate*fFluxScale << " Hz." << newline;
  info << "[ReactorGen]::CalculateTotalRate: Estimated number of events per year :" << newline;
  info << "            :: [ " << fReactorListName << ","  << fNuFlavor << " ] : Rate (unscaled) : " << fEvtRate*3600.0*24.0*365.25 << " evts |  flux scale " << fFluxScale << " = " << fFluxScale*fEvtRate*3600.0*24.0*365.25 << "evts/yr ("<< fFluxScale*fEvtRate*3600.0*24.0*365.25 *1.e32/fNTarget <<" TNU w/ no oscillation)"<<newline;

}

void ReactorGen::BeginOfRun() {
  // Load all the stuff from the database and fill into the necessary structures
  // All the DB accesses should be done at this point, in order to speed up the accesses later
#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::BeginOfRun : Loading generator details." << newline;
#endif

  DBLinkPtr linkdb;
  std::vector<string> reactor_list;
  vector<bool> auto_active_cores;

  try {
    linkdb = DB::Get()->GetLink("REACTOR_LIST",fReactorListName.data());

#ifdef REACTORDEV_DEBUG
    debug << "[ReactorGen]::BeginOfRun : Loading list " << fReactorListName << newline;
#endif
    // Found the reactor list. Load the specific reactors that are flagged as active in the list
    if(fReactorListName!="ALL")
      reactor_list = linkdb->GetSArray("active_reactors");
    else
    {
      DBLinkGroup fLReactorTable = DB::Get()->GetLinkGroup("REACTOR");
      DBLinkGroup::iterator iIterator;
      for (iIterator = fLReactorTable.begin(); iIterator != fLReactorTable.end(); ++iIterator) {
        //if 'ALL' is set then load all the reactor information into the vectors
        reactor_list.push_back(iIterator->first);
      }
    }

    // Set the auto_active_cores to false on all entries by default
    // the auto_active_cores accounts for the case that we want all the cores to be on by default.
    auto_active_cores.resize(reactor_list.size(),false);
  }
  catch (RAT::DBNotFoundError &e) {
    std::ostringstream msg;
    msg << "ReactorGen::BeginOfRun : Failed to load specified reactor list [" << fReactorListName << "] : " << e.what();
    warn << msg.str() << newline;
    Log::Die("ReactorGen::BeginOfRun : Exception caught while loading the reactor list");
  }
#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]::BeginOfRun : List of active reactors: " << newline;
  for (size_t i = 0; i < reactor_list.size(); ++i) {
    debug << "\t" << reactor_list.at(i) << newline;
  }
#endif

  size_t n_cores;
  vector<double> core_power;
  vector<double> core_distance;
  vector<double> latitude;
  vector<double> longitude;
  vector<double> altitude;
  vector<string> core_spectrum;
  vector<int> active_cores;
  // Build a PDF of all the cores, based on the power that each has
  fPowerPDF = new TH2D("PowerPDF","PDF of reactor power",reactor_list.size(),0,reactor_list.size(),max_cores,0,max_cores);
  // -- Now populate the PDF with the core power
  VertexGen_IBD * vtxGen = NULL;

  double weight_prob =  0.0;
  for (unsigned int i = 0; i < reactor_list.size(); ++i) {
    try {

#ifdef REACTORDEV_DEBUG
      debug << "[ReactorGen]::BeginOfRun : Loading reactor " << reactor_list.at(i) << newline;
#endif
      // -- Grab the table for that specific reactor
      linkdb = DB::Get()->GetLink("REACTOR",reactor_list.at(i));
#ifdef REACTORDEV_DEBUG
      debug << "[ReactorGen]::BeginOfRun : Loaded reactor " << linkdb->GetIndex() << newline;
#endif


      Reactor reactor(reactor_list.at(i));
      // Collect the static information about the reactor
      n_cores = linkdb->GetI("no_cores");
      latitude = linkdb->GetDArray("latitude");
      longitude = linkdb->GetDArray("longitude");
      altitude = linkdb->GetDArray("altitude");

      // Now grab the transient reactor data
      linkdb = DB::Get()->GetLink("REACTOR_STATUS",reactor_list.at(i));

      core_power = linkdb->GetDArray("core_power");
      core_spectrum = linkdb->GetSArray("core_spectrum");

      //the PRIS database uses reactor types: PWR, PHBR, FBR, LWGR, BWR, GCR. We have spectra for PWR, BWR & PHBR. We set (FBR,LWGR,GCR)=PWR.
      std::replace(core_spectrum.begin(), core_spectrum.end(), std::string("FBR"), std::string("PWR"));
      std::replace(core_spectrum.begin(), core_spectrum.end(), std::string("LWGR"), std::string("PWR"));
      std::replace(core_spectrum.begin(), core_spectrum.end(), std::string("GCR"), std::string("PWR"));

      // active_cores is a virtual parameter, based on the value of core_power
      for (size_t j = 0; j < core_power.size(); ++j) {
        active_cores.push_back((core_power.at(j)>0.0)?1:0);
      }

      if((int)n_cores > max_cores)
        Log::Die(Form("[ReactorGen]:: %s reactor has %d cores but maximum allowed is %d : change max_cores in ReactorGen.hh",(reactor_list.at(i)).c_str(),(int)n_cores,max_cores),5000);


#ifdef REACTORDEV_DEBUG
      detail << "[ReactorGen]::BeginOfRun : Loaded core information : " << newline;
      detail << "no_cores     : " << n_cores << newline;
#endif
      for (unsigned int j = 0; j < n_cores; ++j) {
        Reactor::CoreData core;
        core.reactor_name = reactor_list.at(i);
        core.core_direction = Direction(longitude.at(j),latitude.at(j),altitude.at(j));
        //in km!
        core.core_distance = GetReactorDistanceLLA(longitude.at(j),latitude.at(j),altitude.at(j))*km;

#ifdef REACTORDEV_DEBUG
        detail << ":: core        : " << j << newline;
        detail << "core_power     : " << core_power.at(j) << newline;
        detail << "latitude       : " << latitude.at(j) << newline;
        detail << "longitude      : " << longitude.at(j) << newline;
        detail << "altitude       : " << altitude.at(j) << newline;
        detail << "core_spectrum  : " << core_spectrum.at(j) << newline;
        detail << "distance [km]  : " << core.core_distance/km << newline;
        detail << newline;
#endif

        core.core_LLA_coord = G4ThreeVector(longitude.at(j),latitude.at(j),altitude.at(j));
        core.core_power = core_power.at(j);
        core.core_spectrum = core_spectrum.at(j);
        if(fReactorListName=="ALL")
          core.is_active = true;
        else
          core.is_active = (active_cores.at(j)==1)?true:false;
        reactor.SetCore(core,j);

        // Check if there is a vertex generator for this core type
        // If there is none, add a new one
        if (fVertexGenMap.find(core.core_spectrum) == fVertexGenMap.end()) {
          // There is no vertex generator for this
          detail << "ReactorGen:: Instantiating a new vertex generator for core type [" << core.core_spectrum << "]" << newline;
          fVertexGenMap[core.core_spectrum] = dynamic_cast<VertexGen_IBD *>(GlobalFactory<GLG4VertexGen>::New(fVertexGenName));
#ifdef REACTORDEV_DEBUG
          debug << "[ReactorGen]::BeginOfRun : Loaded vertex gen " << fVertexGenName << newline;
#endif
          vtxGen = fVertexGenMap[core.core_spectrum];
          vtxGen->SetStandaloneMode(false);

          // Pass the state that by now should already be set to the generator
          if (fVertexStateString.length() != 0) {
#ifdef REACTORDEV_DEBUG
            debug << "ReactorGen::BeginOfRun : Passing state string to newly created vertex generator for spectrum [" << core.core_spectrum << "]" << newline;
            debug << "ReactorGen::BeginOfRun : [" << fVertexStateString << "]"<<newline;
#endif
            vtxGen->SetState(fVertexStateString);
            // this also means that the direction is fixed
            SetNuDirection(fVertexStateString);

          }
        }

        // The weight should only be added if the core is active
        if (core.is_active) {
          // The content should be the power weighted probablity
          // Reproduce the calculations of Jeff Lidgard
          weight_prob = core.core_power/(core.core_distance*core.core_distance);
#ifdef REACTORDEV_DEBUG
          debug << "[ReactorGen]::BeginOfRun : Loading spectrum  " << core.core_spectrum << " with weight prob " << weight_prob << newline;
#endif
          fPowerPDF->SetBinContent(i+1,j+1,weight_prob);


          // Load the spectrum information for this core type
          LoadSpectrum(core.core_spectrum);
        } else {
          fPowerPDF->SetBinContent(i+1,j+1,0.0);
        }
      }
      fReactorData.push_back(reactor);
    }
    catch (RAT::DBNotFoundError &e) {
      std::ostringstream msg;
      msg << "ReactorGen::BeginOfRun : Failed to load details of reactor [" << reactor_list.at(i) << "] : " << e.what();
      warn << msg.str() << newline;
      Log::Die("ReactorGen::BeginOfRun : Exception caught while loading the reactor properties");

    }

  }

  // Calculate the total rate for the selected reactors
  // and accordingly update the time generator
  CalculateTotalRate();

}

void ReactorGen::EndOfRun() {
  // Do the reverse of the BeginOfRun
  // Do it in the reverse order in case there's some dependencies
  // that should be met


  // Empty the list of vertex generators
  std::map<string,VertexGen_IBD*>::const_iterator it = fVertexGenMap.begin();

  for (it = fVertexGenMap.begin(); it != fVertexGenMap.end(); ++it) {
    it->second->EndOfRun();
    delete it->second;
  }
  fVertexGenMap.clear();
  // Clear

  // -- Clean up the reactor data
  fReactorData.clear();

  // -- delete the PDF of active cores
  delete fPowerPDF;

}


void ReactorGen::LoadSpectrum(G4String name) {
#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]::LoadSpectrum : Loading spectrum " <<  name<< newline;
#endif


  // First thing to do is check that this spectrum has already been loaded.
  if (fReactorSpectra.find(name) == fReactorSpectra.end()) {
    try {
#ifdef REACTORDEV_DEBUG
    debug << "[ReactorGen]::LoadSpectrum : Spectrum " <<  name << " does not exist yet. Creating." << newline;
#endif

    // It has not yet been loaded
    DBLinkPtr dblink = DB::Get()->GetLink("REACTOR_SPECTRUM",name);
    vector<double> Enu;
    vector<double> Flux;
    vector<double> composition;
    vector<string> isotope_list;
    string spec_type;
    double emin,emax;
    double FNorm = 1.;  // necessary for Satoko's inputs
    // Load all the necessary entries to construct the spectrum
    spec_type = dblink->GetS("spectrum_type");
    if (ReactorSpectrum::SpectrumNameToType(spec_type) == ReactorSpectrum::parameterized) {
      isotope_list = dblink->GetSArray("param_isotope");
      composition = dblink->GetDArray("param_composition");
    } else if (ReactorSpectrum::SpectrumNameToType(spec_type) == ReactorSpectrum::shape) {
      Enu = dblink->GetDArray("spec_e");
      Flux = dblink->GetDArray("spec_flux");
      FNorm = dblink->GetD("flux_norm");
    } else {
      Log::Die(Form("[ReactorGen]::LoadSpectrum :Invalid spectrum type [%s] for entry [%s]",spec_type.c_str(),name.data()),5000);
    }
    emin = dblink->GetD("emin");
    emax = dblink->GetD("emax");

#ifdef REACTORDEV_DEBUG
    debug << "[ReactorGen]::LoadSpectrum : Loaded entry " << name << " with spectrum type " << spec_type << newline;
    debug << "[ReactorGen]::LoadSpectrum : Building spectrum with the following info:" << newline;
    debug << "Name          : " << name << newline;
    debug << "Spectrum Type : " << spec_type << newline;
    debug << "Emin          : " << emin << newline;
    debug << "Emax          : " << emax << newline;
    debug << "Enu (#)       : " << Enu.size() << newline;
    debug << "Flux (#)      : " << Flux.size() << newline;
    debug << "FluxNorm (#)  : " << FNorm << newline;
    debug << "Isotopes      : ";
    for (size_t i = 0; i < isotope_list.size(); ++i) {
      debug << isotope_list.at(i) << " , ";
    }
    debug << newline;
    debug << "Composition   : ";
    for (size_t i = 0; i < composition.size(); ++i) {
      debug << composition.at(i) << " , ";
    }
    debug << newline;
    debug << "multiplying given flux by flux_norm, to be normalized later" << newline;
#endif
    for(size_t iflux=0; iflux<Flux.size(); iflux++)
      Flux.at(iflux)/=FNorm;

    fReactorSpectra[name] = new ReactorSpectrum(name,ReactorSpectrum::SpectrumNameToType(spec_type),
        emin,emax,
        Enu,Flux,isotope_list,composition);
#ifdef REACTORDEV_DEBUG
    debug << "[ReactorGen]::LoadSpectrum : Spectrum created. Assigning to vertex generator." << newline;
    debug << "[ReactorGen]::LoadSpectrum : Spec address " << fReactorSpectra.at(name) << newline;
    debug << "[ReactorGen]::LoadSpectrum : Hist address " << ((fReactorSpectra.at(name))->GetSpectrum()) << newline;
    debug << "[ReactorGen]::LoadSpectrum : Name : " << (fReactorSpectra.at(name)->GetSpectrum())->GetName() << newline;
#endif
    // Assign this spectrum to the vertex generator
    // Warning!!! Ownership is retained. This should then be cleaned out
    fVertexGenMap.at(name)->SetSpectrum(fReactorSpectra.at(name)->GetSpectrum());

    } catch(DBNotFoundError &e) {
      std::ostringstream msg;
      msg << "ReactorGen::LoadSpectrum : Failed to load details of spectrum [" << name << "] : " << e.what();
      warn << msg.str() << newline;

      throw;
    }
  }
#ifdef REACTORDEV_DEBUG
  else {
    debug << "[ReactorGen]::LoadSpectrum : Spectrum " << name << " already loaded." << newline;
  }
#endif
}

void ReactorGen::SetState(G4String state)
{
  // Accept multiple forms of the state setting
  // First look for specific forms by parsing the first entry
  // if everything else fails assume that the state is the initial state setting that comes through
  // GLG4PrimaryGeneratorAction


#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::SetState : Setting state with input string  [" << state << "] " << newline;
#endif

  // strip weird characters
  state = util_strip_default(state);
#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::SetState : Setting state with input string  (stripped) [" << state << "] " << newline;
#endif

  // now strip into command and rest
  G4String command,rest,gen_name,gen_state;
  pop_first_word(state, command, rest);
  // Direction is passed directly to the vertex gen

  if (command == "time_gen") {
#ifdef REACTORDEV_DEBUG
    info << "[ReactorGen]::SetState : Setting time generator with input string  (stripped) [" << rest << "] " << newline;
#endif

    pop_first_word(rest,gen_name,gen_state);
    if (fTimeGen) {
      delete fTimeGen;
      fTimeGenName = gen_name;
      fTimeGen = RAT::GlobalFactory<GLG4TimeGen>::New(fTimeGenName);
    }
    if (gen_state.length() != 0)
      fTimeGen->SetState(gen_state);
  } else if (command == "pos_gen") {
#ifdef REACTORDEV_DEBUG
    info << "[ReactorGen]::SetState : Setting pos generator with input string  (stripped) [" << rest << "] " << newline;
#endif
    // split the generator name
    pop_first_word(rest,gen_name,gen_state);
    // New generator name, or new generator and state
    if (fPosGen) {
      delete fPosGen;
      fPosGenName = gen_name;
      fPosGen = RAT::GlobalFactory<GLG4PosGen>::New(fPosGenName);
    }
    if (gen_state.length() != 0) {
      fPosGen->SetState(gen_state);
    }
  } else if (command == "reactor_list") {
    // This only allows one other entry
    info << "[ReactorGen]::SetState : Setting reactor list to [" << rest << "]" << newline;
    SetReactorList(rest);
  } else if (command == "flux_scale") {
    char *pEnd;
    double scale = strtod(rest,&pEnd);
    info << "[ReactorGen]::SetState : Setting flux scale to " << scale << newline;
    SetFluxScale(scale);
  } else {
    // This is the state string passed on the command /generator/add reactor [state]
    // Re-Tokenize the state string
    vector<string> parts = util_split(state, ":");

#ifdef REACTORDEV_DEBUG
    for (size_t i = 0; i < parts.size(); i++) {
      detail << "[ReactorGen]::SetState : Part " << i << " :  [" << parts.at(i) << "] " << newline;
    }
#endif

    /// This state is called from GSim. We have another, specific messenger to take care of other options
    /// This is where things are a bit different from the standard combo generator.
    ///
    /// The state choices are the reactor list entry and vertex generator
    /// optionally one can also indicate the time generator and the
    /// position generator

    // Now parse the options according to the selection
    switch (parts.size()) {
    case 4:
      fTimeGenName = parts[3];
      SetTimeGenerator(fTimeGenName,true);
    case 3:
      fPosGenName = parts[2];
      SetPosGenerator(fPosGenName,true);
    case 2:
      fVertexGenName = parts[1];
    case 1:
      fReactorListName = parts[0];
      break;
    default:
      Log::Die( "Reactor generator syntax error: ["+state+"].\n Should follow the pattern <reactor_list>[:vertex_gen][:pos_gen][:time_gen]" );
      break;
    }

    fStateStr = state; // Save for later call to GetState()
  }
}

G4String ReactorGen::GetState() const
{
  return fStateStr + PrintStateHelp();
}

G4String ReactorGen::GetTimeState() const
{
#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]::GetTimeState : time state requested." << newline;
#endif
  if (fTimeGen)
    return fTimeGen->GetState();
  else
    return G4String("ReactorGen error: no time generator selected");
}

void ReactorGen::SetTimeState(G4String state)
{
  detail << "[ReactorGen]::SetTimeState : Detected changes in the time generator. New state : [ " << state << " ] " << newline;
  if (fTimeGen) {
    delete fTimeGen;
    fTimeGen = RAT::GlobalFactory<GLG4TimeGen>::New(fTimeGenName);
    fTimeGen->SetState(state);
  }
  else
    Log::Die("[ReactorGen]::SetTimeState error: no time generator selected.");
}

void ReactorGen::SetPosState(G4String state)
{
  if (fPosGen) delete fPosGen;
  fPosGen = RAT::GlobalFactory<GLG4PosGen>::New(fPosGenName);
  fPosGen->SetState(state);
  volumeInfoLoaded = false;
#ifdef REACTORDEV_DEBUG
  detail << "[ReactorGen]::SetPosState: Change in pos generator detected. Refreshing flux scale recalculation." << newline;
#endif
  // Recalculate the flux based on the changes in the position generator
}

G4String ReactorGen::GetPosState() const
{
#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]::GetPosState : position state requested." << newline;
#endif

  if (fPosGen)
    return fPosGen->GetState();
  else
    Log::Die("[ReactorGen]::GetPosState error: no pos generator selected.");
  return "";
}


G4String ReactorGen::GetVertexState() const
{
#ifdef REACTORDEV_DEBUG
  debug << "[ReactorGen]::GetVertexState : vertex state requested." << newline;
#endif

  return fVertexStateString;
}

void ReactorGen::SetVertexState(G4String state)
{
  debug << "[ReactorGen]::SetVertexState : Detected changes in the vertex generator. New state : [ " << state << " ] " << newline;
  // These options are defined before the run is started, but the Vertex generators are all created at BeginOfRun
  // The best way is to store the state here and then pass to all generators
  fVertexStateString = state;
}

string ReactorGen::PrintStateHelp() {
  string help = "\n\nReactor generator options.\n\n";
  help += "Usage: /generator/set <command> [command options]\n\n";
  help += "Possible options:\n\n";
  help += "time_gen     : Sets the time generator associated with this instance of the reactor generator.\n";
  help += "               Command options: <gen_name> [gen_state(optional)]\n";
  help += "               Example 1: /generator/set time_gen poisson\n";
  help += "               Example 2: /generator/set time_gen uniform 1.0\n";
  help += "pos_gen      : Sets the position generator associated with this instance of the reactor generator.\n";
  help += "               Command options: <gen_name> [gen_state(optional)]\n";
  help += "               Example 1: /generator/set pos_gen fill\n";
  help += "               Example 2: /generator/set pos_gen fill inner_av\n";
  help += "reactor_list : Sets the reactor list to be used by the generator.\n";
  help += "               Command options: <reactor_list> <list_name>\n";
  help += "               Example: /generator/set reactor_list ALL\n";
  help += "flux_scale   : Set a global scale for the neutrino flux. This scales the flux from each reactor by this amount.\n";
  help += "               A scale of 1.0 (default) uses the nominal scale calculated from the reactor power and spectrum.\n";
  help += "               Command options: <flux_scale> <scale_factor>\n";
  help += "               Example: /generator/set flux_scale 100.0\n\n";
  help += "After these options can be mixed with commands for the underlying low level generators.\n";

  return help;
}

} // namespace RAT