#include "PrimaryGeneratorAction.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4Proton.hh"
#include "G4SystemOfUnits.hh"
#include "G4Event.hh"
#include "G4PhysicalConstants.hh"
#include "G4RandomTools.hh"
#include "G4UnitsTable.hh"

#include "ReactionGenerator.hh"

#include <cmath>

PrimaryGeneratorAction::PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction(),
  fParticleGun(nullptr),
  fReaction(nullptr),
  fTargetRadius(20.0 * mm),
  fTargetHalfLength(175.0 * mm)
{
  G4int nParticle = 1;
  fParticleGun = new G4ParticleGun(nParticle);
  fParticleGun->SetParticleDefinition(G4Proton::ProtonDefinition());
  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
  fParticleGun->SetParticleEnergy(1.0 * MeV);

  fReaction = new ReactionGenerator();
}

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  delete fParticleGun;
  delete fReaction;
}

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
{
  ReactionOutput reaction = fReaction->SampleEvent();

  G4double kineticEnergy = reaction.kineticEnergy;
  if (kineticEnergy < 0) kineticEnergy = 0;

  G4ThreeVector direction(reaction.px, reaction.py, reaction.pz);
  if (direction.mag() > 0) {
    direction = direction.unit();
  } else {
    direction = G4ThreeVector(0., 0., 1.);
  }

  G4double r = fTargetRadius * std::sqrt(G4UniformRand());
  G4double phi = twopi * G4UniformRand();
  G4double x = r * std::cos(phi);
  G4double y = r * std::sin(phi);
  G4double z = (2.0 * G4UniformRand() - 1.0) * fTargetHalfLength;

  fParticleGun->SetParticlePosition(G4ThreeVector(x, y, z));
  fParticleGun->SetParticleMomentumDirection(direction);
  fParticleGun->SetParticleEnergy(kineticEnergy);
  fParticleGun->GeneratePrimaryVertex(event);
}