// InHepEvt.cc
// Contact person: Matthew Strait <strait@hep.uchicago.edu>
// See InHepEvt.hh for more details
//———————————————————————//

#include <RAT/InHepEvt.hh>
#include <RAT/Log.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>
#include <RAT/Factory.hh>
#include <RAT/GLG4StringUtil.hh>
#include <RAT/GLG4TimeGen.hh>
#include <RAT/GLG4PosGen.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/MCParticle.hh>

#include <G4Event.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4ParticleTable.hh>
#include <G4ThreeVector.hh>
#include <G4RunManager.hh>

#include <TTree.h>
#include <TFile.h>
#include <TROOT.h>

#include <stdexcept>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>
using namespace std;

namespace RAT {

// The HEPEVT standard doesn't say anything about whether you can have
// blank lines or not, nor does it give any method of having comments.
// Nevertheless, we will tolerate empty lines, lines with only spaces
// and tabs, and lines that look like they are trying to be comments.
// Return true if this line is one of those cases and false otherwise.
static bool IsCommentLine(const string& line)
{
  if(line.find_first_not_of(" \t") == string::npos) return true;

  switch(line[0]){
    case '#': case 'C': case '!': case '\'': case '"': case '/':
      return true;
    default:
      return false;
  }
}

void InHepEvt::ReadHepEvtFile(std::ifstream & hep)
{
  // The global absolute time, incremented each time we read in a
  // particle. If the whole HEPEVT file is supposed to start at a
  // particular global time, that should be given as the delta-t
  // of the first particle, relative to RAT's default time set in
  // data/DATE.ratdb
  double gtime = 0;

  // These must persist across iterations because the HEPEVT standard
  // specifies that if a line doesn't have a entry for mass, delta-t,
  // position (x,y,z) or polarization (x,y,z), you use the value from
  // the most recent line that does. (Is this really supposed to apply
  // across events? It's hard to tell.)
  double deltat = 0;
  HepEvt::Particle part;
  memset(&part, 0, sizeof(part));

  while(true){
    string line;
    HepEvt::Event event;

    if(!std::getline(hep, line)) break;
    if(IsCommentLine(line)) continue;

    int nhep = -1;
    {
      stringstream ss(line);
      if(!(ss >> nhep))
        Log::Die(Form("Bad HEPEVT line:\n%s", line.c_str()));
    }
    if(nhep < 1)
      Log::Die(Form("Bad HEPEVT entry with %d particles", nhep));

    part.time = 0;
    for(int i = 0; i < nhep; i++){
      if(!std::getline(hep, line))
        Log::Die(Form("HEPEVT file ended unexpectedly on particle %d "
          "of %d in event %u", i, nhep, (unsigned)fHepEvts.size()+1));

      if(IsCommentLine(line)) continue;

      stringstream ss(line);

      if(!(ss >> part.stat >> part.pdg >> part.daughter[0]
              >> part.daughter[1] >> part.px >> part.py >> part.pz))
        Log::Die(Form("HEPEVT info missing in \"%s\"\n", line.c_str()));

      // Apparently until C++11, >> will leave the value alone if
      // conversion fails, but starting with C++11, it writes zero into
      // the value. Holy compatibility breaking, Batman! So be paranoid.
      double tmpflt;
      if(ss >> tmpflt) part.mass = tmpflt;
      if(ss >> tmpflt) deltat    = tmpflt; // note special case
      if(ss >> tmpflt) part.x    = tmpflt;
      if(ss >> tmpflt) part.y    = tmpflt;
      if(ss >> tmpflt) part.z    = tmpflt;
      if(ss >> tmpflt) part.plx  = tmpflt;
      if(ss >> tmpflt) part.ply  = tmpflt;
      if(ss >> tmpflt) part.plz  = tmpflt;

      gtime += deltat;

      // Always set the first particle's time in an event to zero and
      // use the delta-t of the first particle to set the event time.
      if(i == 0) event.time = gtime;
      else       part.time += deltat;

      event.part.push_back(part);
    }
    fHepEvts.push_back(event);
    info << fHepEvts.size() << " HEPEVT events read so far" << newline;
  }
  if(fHepEvts.empty()) Log::Die("No events in HEPEVT file!");
}

void InHepEvt::GenerateEvent(G4Event *event)
{
  const HepEvt::Event & evt = fHepEvts[fCurrentEvent];

  std::vector<G4PrimaryParticle*> g4parts;

  // Make all the G4PrimaryParticles out of fHepEvts
  for(unsigned int i = 0; i < evt.part.size(); i++){
    const HepEvt::Particle & part = evt.part[i];

    // HEPEVT files use GeV, but RAT uses MeV.
    G4PrimaryParticle * const particle = new G4PrimaryParticle(
      part.pdg, part.px*1000, part.py*1000, part.pz*1000);
    particle->SetMass(part.mass*1000);

    // To the best of my knowledge, polarization has no effect anywhere
    // in RAT, but we may as well support it here just in case.
    particle->SetPolarization(part.plx, part.ply, part.plz);

    g4parts.push_back(particle);
  }

  // Set the particle daughters, position and time, and push them into
  // the G4Event
  for(unsigned int i = 0; i < evt.part.size(); i++){
    // I'm not sure if anything in RAT cares about particle daughters,
    // but let's set them just in case.
    const HepEvt::Particle & part = evt.part[i];
    for(unsigned int d = 0; d < 2; d++){
      const unsigned int dindex = part.daughter[d];
      if(dindex == 0) continue;
      if(dindex > evt.part.size()){
        warn << "HepEvt: Ignoring out-of-bounds daughter in event "
             << fCurrentEvent << newline;
        continue;
      }

      if(part.stat != 2)
        warn << "HepEvt: Non-decayed particle (status " << part.stat
             << ") with a daughter in event " << fCurrentEvent
             << "\nContinuing, but please check that." << newline;

      g4parts[i]->SetDaughter(g4parts[part.daughter[d]-1]);
    }

    // If the status is not 1, this particle exists only so it can have
    // daughters (or it is some other special case) and so we don't make
    // a G4PrimaryVertex with it.
    if(part.stat != 1) continue;

    G4PrimaryVertex * const vertex = new G4PrimaryVertex(
      G4ThreeVector(part.x, part.y, part.z), part.time);

    vertex->SetPrimary(g4parts[i]);

    event->AddPrimaryVertex(vertex);
  }

  fCurrentEvent++;
  if(fCurrentEvent >= fHepEvts.size()) return BecomeDone();
  fNextTime = fHepEvts[fCurrentEvent].time;
}

void InHepEvt::BecomeDone()
{
  G4RunManager::GetRunManager()->AbortRun(/* soft abort */ true);
  fDone = true;
}

void InHepEvt::ResetTime(double offset)
{
  if(fCurrentEvent < fHepEvts.size()){
    if(fTimeGen){
      fNextTime = fTimeGen->GenerateEventTime() + offset;
    }
    else{
      fNextTime = fHepEvts[fCurrentEvent].time + offset;
      if(fCurrentEvent > 0) fNextTime -= fHepEvts[fCurrentEvent-1].time;
    }
  }
  else{
    fNextTime = 1e300; // Should never be consulted
  }
}

void InHepEvt::SetState(G4String filename)
{
  ifstream infile(filename.data());
  if(!infile.is_open())
    Log::Die(Form("Could not open %s\n", filename.data()));

  ReadHepEvtFile(infile);
  infile.close();
}

void InHepEvt::SetTimeState(G4String state)
{
  if (fTimeGen)
    fTimeGen->SetState(state);
  else
    warn << "InHepEvt error: Cannot set time state, no time generator "
            "selected\n";
}

G4String InHepEvt::GetTimeState() const
{
  if (fTimeGen)
    return fTimeGen->GetState();
  else
    return G4String("InHepEvt error: no time generator selected");
}

void InHepEvt::SetPosState(G4String state)
{
  if (fPosGen)
    fPosGen->SetState(state);
  else
    warn << "InHepEvt error: Cannot set position state, no position "
            "generator selected\n";
}

G4String InHepEvt::GetPosState() const
{
  if (fPosGen)
    return fPosGen->GetState();
  else
    return G4String("InHepEvt error: no pos generator selected");
}


} // namespace RAT