// VertexGen_PhotonBomb.cc
// Contact person:Phil Jones <p.g.jones@qmul.ac.uk>
// See VertexGen_PhotonBomb.hh for more details
//———————————————————————//
#include <RAT/VertexGen_PhotonBomb.hh>
#include <RAT/Log.hh>
#include <Randomize.hh>
#include <G4Event.hh>
#include <G4ParticleTable.hh>
#include <G4PrimaryParticle.hh>
#include <G4PrimaryVertex.hh>
#include <G4ThreeVector.hh>
#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <RAT/PhotonThinning.hh>
using namespace CLHEP;

namespace RAT {
  VertexGen_PhotonBomb::VertexGen_PhotonBomb(const char *arg_dbname)
    : GLG4VertexGen(arg_dbname)
  {
    fOpticalPhoton = G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
    SetState("1 400"); // one photon per event, 400 nm
  }

  VertexGen_PhotonBomb::~VertexGen_PhotonBomb()
  {

  }

  void VertexGen_PhotonBomb::
  GeneratePrimaryVertex(G4Event *event,
                        G4ThreeVector &dx,
                        G4double dt)
  {
    const double photonThinning =  RAT::PhotonThinning::GetFactor();
    if( photonThinning > 1.0 )
      {
        warn << "VertexGen_PhotonBomb::GeneratePrimaryVertex: - You are using photon thinning with factor " << photonThinning
	       << "\n Are you sure you want to do this for an Optical Simulation (could result in loss of accuracy if low number of photons) \n";
      }

    for (int i = 0; i < static_cast<int>( static_cast<double>( fNumPhotons ) / photonThinning ); i++) {
      // Pick direction isotropically
      G4ThreeVector mom;
      double theta = acos(2.0 * G4UniformRand() - 1.0);
      double phi = G4UniformRand() * twopi;
      mom.setRThetaPhi(fEnergy, theta, phi); // Momentum == energy units in GEANT4

      G4PrimaryVertex* vertex= new G4PrimaryVertex(dx, dt);
      G4PrimaryParticle* photon = new G4PrimaryParticle(fOpticalPhoton,
                                                        mom.x(),mom.y(),mom.z());
      // Generate random polarization
      phi = (G4UniformRand()*2.0-1.0)*pi;
      G4ThreeVector e1 = mom.orthogonal().unit();
      G4ThreeVector e2 = mom.unit().cross(e1);
      G4ThreeVector pol = e1*cos(phi)+e2*sin(phi);
      photon->SetPolarization(pol.x(), pol.y(), pol.z());
      photon->SetMass(0.0); // Seems odd, but used in GLG4VertexGen_Gun

      vertex->SetPrimary(photon);
      event->AddPrimaryVertex(vertex);
    } // loop over photons
  }

  void VertexGen_PhotonBomb::SetState(G4String newValues)
  {
    if (newValues.length() == 0) {
      // print help and current state
      G4cout << "Current state of this VertexGen_PhotonBomb:\n"
             << " \"" << GetState() << "\"\n" << G4endl;
      G4cout << "Format of argument to VertexGen_PhotonBomb::SetState: \n"
        " \"num_photons wavelength_nm\"\n" << G4endl;
      return;
    }

    std::istringstream is(newValues.c_str());
    int num, wavelength;
    is >> num >> wavelength;
    if (is.fail())
      Log::Die("VertexGen_PhotonBomb: Incorrect vertex setting " + newValues);

    fNumPhotons = num;
    fEnergy = hbarc * twopi / (wavelength * nm);
  }

  G4String VertexGen_PhotonBomb::GetState()
  {
    return dformat("%d\t%f", fNumPhotons, fEnergy);
  }

} // namespace RAT