// VertexGen_AlphaN.cc
// Contact person: Aleksandra Bialek <abialek@ualberta.ca>  and Ian Coulter (icoulter@hep.upenn.edu)
// See VertexGen_AlphaN.hh for more details
//———————————————————————//
#include <CLHEP/Units/PhysicalConstants.h>
using namespace CLHEP;

#include <RAT/VertexGen_AlphaN.hh>
#include <RAT/Log.hh>

#include <G4ParticleTable.hh>
#include <Randomize.hh>

#include <TMath.h>

#include <math.h>
#include <cmath>
#include <sstream>
#include <iostream>

using std::string;
#define G4std  std

namespace RAT
{

 VertexGen_AlphaN::VertexGen_AlphaN(const char *arg_dbname):GLG4VertexGen(arg_dbname)
 {
     fDBName=0;
     fMaxLevels=0;
 }

 ///-------------------------------------------------------------------------
 VertexGen_AlphaN::~VertexGen_AlphaN()
 {

 }
 ///-------------------------------------------------------------------------
 void VertexGen_AlphaN::BeginOfRun() {
   //now check if DB file have information for given material and target
   try {
      fDBName=DB::Get()->GetLink(fTarget,fMaterial);
   } catch (DBNotFoundError &e) {
      Log::Die("AlphaN: Cannot find"+fMaterial);
   }

   try {
      fDeexLevels = fDBName->GetIArray("DeexLevels");
   } catch (DBNotFoundError &e) {
      Log::Die("no DeexLevels");
   }
   fMaxLevels = fDeexLevels.size();

   try {
      fProb = fDBName->GetDArray("Prob");
   } catch (DBNotFoundError &e) {
      Log::Die("no Prob");
   }
   fDeexEn.resize(fMaxLevels-1);
   fDeexT.resize(fMaxLevels-1);
   fDeexPart.resize(fMaxLevels-1);

   string nameDB;
   fNotime=false;
   for(int i=0; i<fMaxLevels-1; i++)//note: there is no ground state here.
   {
      std::ostringstream levl;
      levl<<(fDeexLevels[i]+1);

      nameDB="DeexEn_"+levl.str();// energy of the deexcitation particle
      fDeexEn[i]  = fDBName->GetDArray(nameDB);

      nameDB="DeexPart_"+levl.str();// deexcitation particle name
      fDeexPart[i] = fDBName->GetSArray(nameDB);

      nameDB="DeexT_"+levl.str();   // halflife time of the deexcitation level (if given)
      try {
         fDeexT[i] = fDBName->GetDArray(nameDB);
      } catch (DBNotFoundError &e) {
         fNotime=true;
         info<<"time 0 taken";
      }
   }


   char histName[90];
   string histLevel;

   fNeutronEn.resize(fMaxLevels);
   fNeutronInten.resize(fMaxLevels);
   for(int i=0; i<fMaxLevels; i++)
   {
      std::ostringstream lev;
      lev<<fDeexLevels[i];
      histLevel="EnNeutron_"+lev.str();
      try {
         fNeutronEn[i]=fDBName->GetDArray(histLevel);
      } catch (DBNotFoundError &e) {
         Log::Die("There is no spectrum of neutrons (Energy)");
      }
      histLevel="Intens_"+lev.str();
      try{
         fNeutronInten[i]=fDBName->GetDArray(histLevel);
      } catch (DBNotFoundError &e) {
         Log::Die("There is no spectrum of neutrons (Prob)");
      }

      //load the neutrons spectra to histograms
      int maxbin= fNeutronInten[i].size();

      sprintf(histName,"nSpectra_%d",i);
      fSpectra[i] = new TH1F (histName,histName,maxbin,0,10);

      for(int j=0; j<maxbin; j++)
      {
         fSpectra[i]->Fill(fNeutronEn[i][j],fNeutronInten[i][j]);
      }
   }

   fLoadHistos= true;

   for(int i=0; i<fMaxLevels; i++)
   {
      std::ostringstream nameLow;
      nameLow << "EAlphaLow" << i;
      try{
         fAlphaLimitLowArray.push_back(fDBName->GetDArray(nameLow.str()));
      } catch (DBNotFoundError &e) {
         Log::Die("Alpha limits not found.");
      }

      std::ostringstream nameHigh;
      nameHigh << "EAlphaHigh" << i;
      try{
         fAlphaLimitHighArray.push_back(fDBName->GetDArray(nameHigh.str()));
      } catch (DBNotFoundError &e) {
         Log::Die("Alpha limits not found.");
      }
   }

 }
 ///-------------------------------------------------------------------------
 void VertexGen_AlphaN::GeneratePrimaryVertex(G4Event *event, G4ThreeVector &dx, G4double dt)
 {
    G4int deexPart  = 0, prob=0;
    G4int deexlevel=-1;
    G4String   particleName;
    double nEnergy=0.0, ke=0.0, time=0.0;
    bool notime2=false;

    if(!fLoadHistos)
    {
       info<<"Check the data base file";
       Log::Die("No neutron spectrum found");
       return;
    }

    prob = 100*G4UniformRand(); // random generator over the energy levels
    for(int i=0; i<fMaxLevels; i++)
    {
       if(prob<=fProb[i])
       {
          nEnergy = fSpectra[i]->GetRandom(); //random energy from the neutron spectrum
          if(i!=0)//for all levels but ground state
          {
             deexPart = fDeexEn[i-1].size();
             deexlevel= i-1;
          }
          else     deexPart=0; //for the ground state there is no deexcitation particles

          break;
       }
    }
    G4ThreeVector mom;
    //setting information for particles
    for(int j=0; j<deexPart+1; j++)
    {
       if(j==0)// information for neutron
       {
          particleName = "neutron";
          ke = nEnergy;
          time = 0.0;
       }
       else //information for each deexcitation particle if exists
       {
          particleName = fDeexPart[deexlevel][j-1];
          ke = fDeexEn[deexlevel][j-1];

          if (fDeexEn[deexlevel].size()>fDeexT[deexlevel].size()) notime2=true;

          if(fNotime) time = 0.0;
          else if(!(notime2 && particleName=="e-")) time = -(fDeexT[deexlevel][j-1]/log(2.))*log(G4UniformRand());

       }

       mom = GetMomentum(ke, particleName);
       G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
       // Create primary particle
       G4PrimaryParticle* particle = new G4PrimaryParticle(theParticleTable->FindParticle(particleName),mom.x()*MeV,mom.y()*MeV,mom.z()*MeV);
       particle->SetProperTime(time);
       G4PrimaryVertex* vertex = new G4PrimaryVertex (dx,time*ns+dt);
       vertex->SetPrimary(particle);
       event->AddPrimaryVertex(vertex);
     }

    // Now create alpha based on neutron energy, get the energy of alpha
    // (or at least, the energy it deposits in the scintillator)
    double aEnergy = GetAlphaEnergy(nEnergy, deexlevel);

    // If the alpha energy is 0, don't bother generating an alpha
    if(aEnergy!=0){
      // Create the alpha vertex at the start time of the event and same position as the neutron
      G4ThreeVector mom= GetMomentum(aEnergy, "alpha");
      G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
      G4PrimaryParticle* particle = new G4PrimaryParticle(theParticleTable->FindParticle("alpha"),mom.x()*MeV,mom.y()*MeV,mom.z()*MeV);
      particle->SetProperTime(0.0);
      G4PrimaryVertex* vertex = new G4PrimaryVertex(dx,dt);
      vertex->SetPrimary(particle);
      event->AddPrimaryVertex(vertex);
    }
 }

 double VertexGen_AlphaN::GetAlphaEnergy(double nEnergy, G4int deexlevel)
 {
   // First find the bin containing the neutron energy
   std::vector<double>::iterator nBin = std::lower_bound(fNeutronEn[0].begin(),fNeutronEn[0].end(),nEnergy);
   int iBin = (int) std::distance(fNeutronEn[0].begin(),nBin);

   // Get the limits for this excitation state
   double energyLo = GetAlphaLimits("EAlphaLow",deexlevel,iBin);
   double energyHi = GetAlphaLimits("EAlphaHigh",deexlevel,iBin);

   if(energyHi==0) return 0.0;

   // Choose an alpha energy from a random uniform distribution between the low and high bounds
   double energy= energyLo + (energyHi-energyLo)*G4UniformRand();
   return energy;
 }

 double VertexGen_AlphaN::GetAlphaLimits(G4String limit, G4int deexlevel, int iBin)
 {
   if(limit="EAlphaLow")
      return fAlphaLimitLowArray[deexlevel+1][iBin];
   else return fAlphaLimitHighArray[deexlevel+1][iBin];
 }

 ///-------------------------------------------------------------------------
 G4ThreeVector VertexGen_AlphaN::GetMomentum(double ke, G4String partName)
 {
    G4ThreeVector dir, mom;
    // Random distribution uniform in solid angle
    G4double cos_theta = 2. * G4UniformRand() - 1.;
    G4double sin_theta = sqrt(1. - pow(cos_theta,2));
    G4double phi = twopi * G4UniformRand();

    dir.setX(sin_theta * cos(phi));
    dir.setY(sin_theta *sin(phi));
    dir.setZ(cos_theta);

    fPartDef=G4ParticleTable::GetParticleTable()->FindParticle(partName);
    double mass=fPartDef->GetPDGMass();

    G4double en = ke + mass;
    G4double pmom = sqrt(pow(en,2) - pow(mass,2));
    mom.setX(pmom * dir.x());
    mom.setY(pmom * dir.y());
    mom.setZ(pmom * dir.z());

    return mom;
 }

 ///-------------------------------------------------------------------------
 void VertexGen_AlphaN::SetState(G4String newValues)
 {
     if(newValues.length() == 0) // no arguments given in the macro
     {
        info << "Format of argument to VertexGen_AlphaN::SetState: "
             <<  GetState()<< " \n"
             << "/generator/vtx/set [MATERIAL] [TARGET] \n"
             << "Neutron spectra given for MATERIALs:\n"
             << "LABPPO, TeLABPPO, AV\n"
             << "TARGETs are: 13C and 16O";

        Log::Die("Please set correctly parameters for the generator "+GetState());
        return;
     }
     // arguments indicating the material of the volume where the event is simulated (LAB, TeLAB or AV) and the target of alpna-n reaction (13C or 18O)


     std::istringstream is(newValues.c_str());
     G4String material, target;

     // argument [MATERIAL]
     is>> material;
     fMaterial = material;
     // argument [TARGET]
     is>>target;

     if(is.fail())
        fTarget = "13C";
     else
        fTarget = target;

     info << "Arguments of the VertexGen_AlphaN: \n"
          <<  GetState();

     // Following this, the histograms are creating in the BeginOfRun
 }

 ///-------------------------------------------------------------------------
 G4String VertexGen_AlphaN::GetState()
 {

    // Settings in the input macro
    G4std::ostringstream os;
    os << fMaterial << " " <<fTarget<<"\n"<<G4std::ends;
    G4String rv(os.str());
    return rv;
 }

}//end of namespace RAT