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

#include <RAT/SolarGen.hh>
#include <RAT/SolarMessenger.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/SunUtil.hh>

// -- ROOT includes
#include <TGraph.h>
#include <TMath.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 <G4SystemOfUnits.hh>

#include <stdexcept>

namespace RAT {

SolarGen::SolarGen() : GLG4Gen("solar"), fStateStr(""), fTimeGenName("poisson"),
    fPosGenName("fill"),fUseEphemerisData(true),fNuDir(0.0,0,-1.0),fFluxScale(1.0),volumeInfoLoaded(false)

{
  fEvTime.SetStart();
  SetPosGenerator(fPosGenName,true);
  SetTimeGenerator(fTimeGenName,true);

#ifdef RATDEBUG
  fUTtime = 0;

  debug << "[SolarGen]:: In constructor." << newline;
#endif
  fMessenger = new SolarMessenger(this);
  fRatePerTarget = 0.0;
}


SolarGen::~SolarGen()
{

#ifdef RATDEBUG
  debug << "[SolarGen]:: In destructor." << newline;
#endif

  // Clean up everything built in this generator
  if (fMessenger) delete fMessenger;

}

void SolarGen::GenerateEvent(G4Event *event)
{


#ifdef RATDEBUG
  debug << "[SolarGen]:: 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();
  GenerateNuDirection();
#ifdef RATDEBUG
  detail << "[SolarGen]::GenerateEvent : Incoming neutrino direction set to  [ " << fNuDir.x() << " " << fNuDir.y() << " " << fNuDir.z() << " ]." << newline;
  detail << "[SolarGen]::GenerateEvent : Event time (before time generator offset): [" <<  fEvTime.GetStrTime() << "] date : [" << fEvTime.GetStrDate() << "] UTtime " << fUTtime<< "."<< newline;
  detail << "[SolarGen]::GenerateEvent : Event position ( " << pos.x()/cm << ","<< pos.y()/cm << ","<< pos.z()/cm << ") cm " << newline;
#endif

  dynamic_cast<VertexGen_ES*>(fVertexGen)->GeneratePrimaryVertex(event, pos, t0, fNuDir);
}

void SolarGen::ResetTime(double offset)
{
  //TODO: One could add here corrections to the time based on the eccentricity of the Earth
  // The effect is too small to make any difference
  fNextTime = fTimeGen->GenerateEventTime() + offset;
  fEvTime.IncrementUT(fNextTime);
#ifdef RATDEBUG
  fUTtime += fNextTime;

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

void SolarGen::SetState(G4String state)
{
  /// 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 flux and the flavor

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

  state = util_strip_default(state);
#ifdef RATDEBUG
  debug << "[SolarGen]::SetState : Setting state with input string after stripping [" << state << "] " << newline;
#endif

  std::vector<std::string> parts = util_split(state, ":");
#ifdef RATDEBUG
  detail << "[SolarGen]::SetState : Input state split into " << parts.size() << " parts." << newline;
  for (unsigned int i = 0; i < parts.size(); ++i) {
    detail << "[SolarGen]::SetState : [" << parts[i] << "]." << newline;
  }
#endif

  // There can be 2 options
  // <flux>:<flavour:{e:mu}>

  VertexGen_ES *local_vtx = NULL;

  switch (parts.size()) {
  case 2:
    fFlux = parts[0]; fNuFlavor = parts[1];
#ifdef RATDEBUG
    detail << "[SolarGen]::SetState : Found choice for flux=" << fFlux << " and nu flavour=" << fNuFlavor << newline;
#endif

    if (fVertexGen) {
      delete fVertexGen;
      fVertexGen = 0;
    }
    // I shouldn't need to go through the GlobalFactory to get the vertex generator
    fVertexGen = new VertexGen_ES("solar");
    local_vtx = dynamic_cast<VertexGen_ES*>(fVertexGen);
    // disallow modifications through the messenger.
    local_vtx->SetEnableMessenger(false);
    // Now specify the proper options to the generator and only then load it.
    local_vtx->SetFlux(fFlux);
    local_vtx->SetNuFlavor(fNuFlavor);
    SetFluxScale(fFluxScale);
    break;
  default:
    Log::Die( "Solar generator syntax error: ["+state+"]. Should only <flux>:<flavor>" );
    break;
  }

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



void SolarGen::GenerateNuDirection() {
  // If the direction is fixed do nothing.
  // If we are using the ephemeris data pick up the start date and generate a Julian date to get the direction.
  // Then do a simple random generation of the date and that will give us the direction.
  if (fUseEphemerisData) {
    // The -1 scaling factor is to account that the neutrino is coming
    // *towards* the detector
    // Transform fNextTime into a date to calculate the direction
    G4int days = fEvTime.GetUTDays();
    G4int secs = fEvTime.GetUTSecs();
    G4int nsecs = fEvTime.GetUTNsecs();

    TVector3 tmp_dir = -1.0*RAT::SunDirection(days,secs,nsecs);
    fNuDir = G4ThreeVector(tmp_dir.x(),tmp_dir.y(),tmp_dir.z());
    // Sample from the livetime distribution to get the zenith angle

#ifdef RATDEBUG
    detail << "[SolarGen]::GenerateNuDirection : Time set." << newline;
    detail << "[SolarGen]::DATE: " << days << " " << secs << " " << nsecs << newline;
    detail << "[SolarGen]::DIRECTION: " << fNuDir.x() << " " << fNuDir.y() << " " << fNuDir.z() << newline;
#endif
  }
}


G4String SolarGen::GetState() const
{
  return fStateStr;
}

G4String SolarGen::GetTimeState() const
{
#ifdef RATDEBUG
  debug << "[SolarGen]::GetTimeState : time state requested." << newline;
#endif
  if (fTimeGen)
    return fTimeGen->GetState();
  else
    return G4String("SolarGen error: no time generator selected");
}
/** This may be reviewed and re-incorporated later.
 * For the moment keep is disallowed to avoid confusion with messing up with the generator
void SolarGen::SetTimeState(G4String state)
{
  if (fTimeGen)
    fTimeGen->SetState(state);
  else
    G4cerr << "SolarGen error: Cannot set time state, no vertex generator selected" << newline;
}

void SolarGen::SetVertexState(G4String state)
{
  if (fVertexGen)
    fVertexGen->SetState(state);
  else
    error << "SolarGen error: Cannot set vertex state, no vertex generator selected" << newline;
}

G4String SolarGen::GetVertexState() const
{
  if (fVertexGen)
    return fVertexGen->GetState();
  else
    return G4String("SolarGen error: no vertex generator selected");
}


void SolarGen::SetPosState(G4String state)
{
  if (fPosGen)
    fPosGen->SetState(state);
  else
    error << "SolarGen error: Cannot set position state, no position generator selected" << newline;
}
 */

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

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




void SolarGen::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 SolarGen::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 SolarGen::SetNuDirection(G4ThreeVector dir) {
  // If this method is called then the incoming neutrino direction is
  // fixed and therefore we won't use ephemeris data to determine
  // the direction to the Sun
  fNuDir = dir;
  fUseEphemerisData = false;
}

void SolarGen::SetFluxScale(G4double scale) {

  fFluxScale = scale;
  // declare them as static since they are common to all instances of the class
  static G4double theVolume = 0.0;
  static G4double theNelec = 0.0;
#ifdef RATDEBUG
  /* Used only for debug*/
  static G4double theMass = 0.0;
#endif

  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 RATDEBUG
    detail << "[SolrGen]::SetFluxScale : 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) {
      detail << "[SolarGen]::SetFluxScale : Interactions will be generated in volume " << thePhyVolume->GetName() << newline;
    } else {
      Log::Die("[SolarGen]::SetFluxScale : 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(); // in cubic meters
#ifdef RATDEBUG
    theMass = theLogVolume->GetMass(false,false);
#endif


    // Get the electron density
    theNelec = theMaterial->GetElectronDensity();

    // rate per target. Only dependent on flux type and neutrino type
    fRatePerTarget = dynamic_cast<VertexGen_ES*>(fVertexGen)->GetHelper()->GetRatePerTarget();
#ifdef RATDEBUG
    detail << "[SolarGen]::SetFluxScale : Rate per Target "
        << fRatePerTarget  << newline;
    detail << "[SolarGen]::SetFluxScale : Number of targets " << theNelec << newline;
#endif
    volumeInfoLoaded = true;
  }

  fEvtRate = fFluxScale*theVolume*theNelec*fRatePerTarget;
  G4String newstate = util_to_string(fEvtRate);
  // Pass this state also to the time generator
  fTimeGen->SetState(newstate);
#ifdef RATDEBUG
  detail << "[SolarGen]::SetFluxScale: Passed state " << newstate << " to time generator." << newline;

  info << "[SolarGen]::SetFluxScale: SSM scale " << fFluxScale << newline;
  info << "[SolarGen]::SetFluxScale: Volume " << theVolume/m3 << " m^3 ; Mass " << theMass/kg << " kg." << newline;
  info << "[SolarGen]::SetFluxScale: Nelec " << theNelec/cm3 << " cm^-3 rate_target " << fRatePerTarget << newline;
  info << "[SolarGen]::SetFluxScale: Total targets : " << theVolume*theNelec << newline;

#endif
  detail << "[SolarGen]::SetFluxScale: Event rate updated to " << fEvtRate << " Hz." << newline;
  info << "[SolarGen]::SetFluxScale: Estimated number of events per year : [ " << fFlux << ","  << fNuFlavor << " ] : SSM scale " << fFluxScale << " : " << fEvtRate*3600.0*24.0*365.25 << newline;

}

void SolarGen::SetPosState(G4String state)
{
  if (fPosGen) delete fPosGen;
  fPosGen = RAT::GlobalFactory<GLG4PosGen>::New(fPosGenName);
  fPosGen->SetState(state);
  volumeInfoLoaded = false;
#ifdef RATDEBUG
  detail << "[SolarGen]::SetPosState: Change in pos generator detected. Refreshing flux scale recalculation." << newline;
#endif
  SetFluxScale(fFluxScale);
}

} // namespace RAT