// VertexGen_ProtonES.cc
// Contact person:Christopher Jones <christopher.jones@physics.ox.ac.uk>
// See VertexGen_ProtonES.hh for more details
//———————————————————————//
// - Local include
#include "VertexGen_ProtonES.hh"
#include "ProtonESgen.hh"

using namespace CLHEP;

// - RAT includes
#include <RAT/GLG4StringUtil.hh>
#include <RAT/Log.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>

// - Geant4 includes
#include <G4ParticleTable.hh>
#include <G4PrimaryVertex.hh>
#include <G4PrimaryParticle.hh>
#include <G4Event.hh>
#include <CLHEP/Units/PhysicalConstants.h>
#include <Randomize.hh>


namespace RAT {

VertexGen_ProtonES::VertexGen_ProtonES(const char *arg_dbname): GLG4VertexGen(arg_dbname), fNuDir(0.,0.,0.)
{
  // Particles involved : Not sure if I need to include all flavours of neutrino explicitly
  fProton = G4ParticleTable::GetParticleTable()->FindParticle("proton");

  // Other flavours of neutrino (and anti_neutrino) are simulated with this generator. They are all labeled as nu_e as no label for a "general" neutrino exists in Geant4.
  // The flavour of neutrino (or anti-neutrino) is specified in the RAT macro with only one type of neutrino specified per RAT macro. 
  fNu = G4ParticleTable::GetParticleTable()->FindParticle("nu_e");

}

VertexGen_ProtonES::~VertexGen_ProtonES()
{
  //Destructor
}

void VertexGen_ProtonES::GeneratePrimaryVertex(G4Event* argEvent,G4ThreeVector& nudir, G4double nTime)
{
  G4PrimaryVertex* vertex= new G4PrimaryVertex(nudir,nTime);
  G4ThreeVector ev_nu_dir(fNuDir); // By default use specified direction

  if (ev_nu_dir.mag2() == 0.0)
  {
    // Pick isotropic direction
    double theta = acos(2.0 * G4UniformRand() - 1.0);
    double phi = G4UniformRand() * twopi;
    ev_nu_dir.setRThetaPhi(1.0, theta, phi);
  }

  G4LorentzVector mom_proton,mom_nu;

  fProtonESgen.GenerateEvent(ev_nu_dir , mom_nu, mom_proton);

  G4PrimaryParticle* protonParticle =
    new G4PrimaryParticle(fProton,             			// particle code
			  mom_proton.px(),     			// x component of momentum
			  mom_proton.py(),     			// y component of momentum
			  mom_proton.pz());    			// z component of momentum

  protonParticle->SetMass(fProton->GetPDGMass());
  vertex->SetPrimary(protonParticle);


  G4PrimaryParticle* neutrinoParent =
    new G4PrimaryParticle(fNu,     // particle code
              mom_nu.px(),         // x component of momentum
              mom_nu.py(),         // y component of momentum
              mom_nu.pz());        // z component of momentum
  neutrinoParent->SetMass(fNu->GetPDGMass());
  vertex->SetPrimary(neutrinoParent);

  ParentPrimaryVertexInformation *vertinfo = new ParentPrimaryVertexInformation();
  vertinfo->SetVertexParentParticle(neutrinoParent);
  vertex->SetUserInformation(vertinfo);
  argEvent->AddPrimaryVertex(vertex);

}//end of Generate Primary Vertex

void VertexGen_ProtonES::SetState(G4String newValues)
{
  newValues = util_strip_default(newValues); // from GLG4StringUtil
  if (newValues.length() == 0)
  {
    info << "Current state of this VertexGen_ProtonES:\n" << " \"" << GetState() << "\"\n";
    info <<  "Format of argument to VertexGen_ProtonES::SetState: \n"
      " \"fNuDir_x fNuDir_y fNuDir_z\"\n"
      " where fNuDir is the initial direction of the reactor antineutrino.\n"
      " Does not have to be normalized.  Set to \"0. 0. 0.\" for isotropic\n"
      " neutrino direction.\n";
    return;
  }

  std::istringstream is(newValues.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();
    }
}//end of SetState

G4String VertexGen_ProtonES::GetState()
{
  std::ostringstream os;
  os << fNuDir.x() << "\t" << fNuDir.y() << "\t" << fNuDir.z() << std::ends;
  G4String rv(os.str());
  return rv;
}



} // namespace RAT