////////////////////////////////////////////////////////////////////////
// VertexGen_Flasher.cc
// Contact person: John Walker <john.walker@liverpool.ac.uk>
// See VertexGen_Flasher.hh for more details
//
// This file contains code used by Gen_Flasher.cc to construct event vertices
////////////////////////////////////////////////////////////////////////

#include <RAT/VertexGen_Flasher.hh>
#include <RAT/GLG4VertexGen.hh>
#include <RAT/Log.hh>
#include <RAT/PhotonThinning.hh>
#include <RAT/DB.hh>
#include <RAT/EventInfo.hh>
#include <RAT/DS/Calib.hh>
#include <RAT/ParentPrimaryVertexInformation.hh>

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

#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>

#include <Randomize.hh>

#include <vector>
#include <cmath>
using namespace std;

namespace RAT {


VertexGen_Flasher::VertexGen_Flasher(const char *arg_dbname) : GLG4VertexGen(arg_dbname)
{
  // Not yet initialised - call the Init function to setup on first event
  fInitDB = false;
}

void VertexGen_Flasher::Init()
{
  // This function should be run once on the first event (once DB all set up)
  warn << "VertexGen_Flasher::Init: Initialising simulation settings \n";

  // Photon definition (assume /run/initialize already done)
  fPhotonDef = G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
  // Is there a photon thinning factor - this doesn't necessarily make sense for optical simulations
  G4double photon_thin =  RAT::PhotonThinning::GetFactor();
  if(photon_thin > 1.0){
    // OUTPUT WARNING
    warn << "VertexGen_Flasher::Init: !!!WARNING!!! - You are using photon thinning with factor " << photon_thin << "\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";
  }

  // Load database values
  DBLinkPtr lFlasher = DB::Get()->GetLink("FLASHER");

  // Optional setting
  try {
    G4String temp = lFlasher->GetS("pulse_mode");
    if(temp == "poisson"){
      fPoisson = true;
    }else if(temp == "fixed"){
      fPoisson = false;
    }else{
      Log::Die(dformat("VertexGen_Flasher::Init: Flasher.pulse_mode = \"%s\" is invalid.", temp.c_str()));
    }
  } catch (DBNotFoundError &e) { };

  // Wavelength info
  try{
    fFlasherWavelength = lFlasher->GetS("wl_dist");
    info << "VertexGen_Flasher::Init WARNING: Setting wavelength distribution to " + fFlasherWavelength + "\n";
  } catch (DBNotFoundError &e) {
    info << "VertexGen_Flasher::Init: Using single wavelength distribution as set in Flasher.ratdb\n";
    fFlasherWavelength = "mono";
  }
  fMono_wl = 0.0;
  if(fFlasherWavelength=="mono"){
    fMono_wl = lFlasher->GetD("mono_wl");
    info << "VertexGen_Flasher::Init WARNING: Setting single wavelength to " << fMono_wl << "\n";
  }
  if(fFlasherWavelength=="LED"||fFlasherWavelength=="flat"){
    if(fFlasherWavelength=="LED"){
      // Get wavelengths
      fwave = lFlasher->GetDArray("led_dist_wl");
      // Now get the amplitudes
      fwaveInt = lFlasher->GetDArray("led_dist_wl_intensity");
    }
    if(fFlasherWavelength=="flat"){
      // Get wavelengths
      fwave = lFlasher->GetDArray("flat_dist_wl");
      // Now get the amplitudes
      fwaveInt = lFlasher->GetDArray("flat_dist_wl_intensity");
    }
    fMaxWl = fwave.back();
    fMinWl = fwave.front();
    // make a check that the values in the dist_wl make sense and increase monotonically
    for(unsigned int iv=1;iv<fwave.size(); iv++){
      if(fwave[iv]<=fwave[iv-1]){
        warn << "VertexGen_Flasher::Init: "+fFlasherWavelength+"_dist_wl[" << iv << "] = " << fwave[iv] << " [" << iv-1 << "] = "  << fwave[iv-1] << "!\n";
        Log::Die("VertexGen_Flasher::Init: "+fFlasherWavelength+"_dist_wl must increase monotonically. Please check values.");
      }
    }
    // now create a cumulative magnitude vector
    G4double wavemagsum = 0;
    fwaveCumMag.clear();
    fwaveCumMag.push_back(0.0);
    detail << " Cumulative wavelength ";
    for(unsigned int i=1;i<fwave.size();i++){
      wavemagsum +=(fwaveInt[i]+fwaveInt[i-1])/2.0 * (fwave[i]-fwave[i-1]);
      fwaveCumMag.push_back(wavemagsum);
      detail << fwaveCumMag[i] << " ";
    }
    detail << "\n";
  }
  if(fFlasherWavelength=="helium"){
    // Get wavelengths
    fwave = lFlasher->GetDArray("helium_dist_wl");
    // Now get the probabilities
    fwaveProb = lFlasher->GetDArray("helium_dist_wl_prob");
    fMaxWl = fwave.back();
    fMinWl = fwave.front();
    // make a check that the wavelength values make sense and increase monotonically
    for(unsigned int iv=1;iv<fwave.size(); iv++){
      if(fwave[iv]<=fwave[iv-1]){
        warn << "VertexGen_Flasher::Init: helium_dist_wl[" << iv << "] = " << fwave[iv] << " [" << iv-1 << "] = "  << fwave[iv-1] << "!\n";
        Log::Die("VertexGen_Flasher::Init: helium_dist_wl must increase monotonically. Please check values.");
      }
    }
    // now create a cumulative magnitude vector
    G4double waveprobsum = 0;
    fwaveCumProb.clear();
    fwaveCumProb.push_back(0.0);
    detail << " Cumulative probability ";
    for(unsigned int i=0;i<fwave.size();i++){
      waveprobsum +=fwaveProb[i];
      fwaveCumProb.push_back(waveprobsum);
      detail << fwaveCumProb[i] << " ";
    }
    detail << "\n";
  }
  if(fFlasherWavelength=="blackbody"){
    // Get temperature
    fBlackbody_T = lFlasher->GetI("blackbody_temp");
    info << "VertexGen_Flasher::Init WARNING: Setting blackbody temperature to " << fBlackbody_T << "\n";
    fwave.clear();
    for(unsigned int i=200; i<=800; i=i+10){
      fwave.push_back(i);
    }
    // now create intensity vector
    G4double intensitysum = 0;
    fwaveInt.clear();
    for(unsigned int i=0;i<fwave.size();i++){
      G4double S = ((2.0*CLHEP::pi*CLHEP::c_squared*CLHEP::h_Planck)/(pow(fwave[i]*CLHEP::nm,5)))*(1/(exp((CLHEP::h_Planck*CLHEP::c_light)/(fwave[i]*CLHEP::nm*CLHEP::k_Boltzmann*fBlackbody_T))-1));
      intensitysum += S;
      fwaveInt.push_back(S);
    }
    detail << " Blackbody intensity ";
    for(unsigned int i=0;i<fwaveInt.size();i++){
      fwaveInt[i] = fwaveInt[i]/intensitysum;
      detail << fwaveInt[i] << " ";
    }
    detail << "\n";
    // now create a cumulative magnitude vector
    G4double wavemagsum = 0;
    fwaveCumMag.clear();
    fwaveCumMag.push_back(0.0);
    detail << " Cumulative wavelength ";
    for(unsigned int i=1;i<fwave.size();i++){
      wavemagsum +=(fwaveInt[i]+fwaveInt[i-1])/2.0 * (fwave[i]-fwave[i-1]);
      fwaveCumMag.push_back(wavemagsum);
      detail << fwaveCumMag[i] << " ";
    }
    detail << "\n";
  }

  // Intensity settings
  try{
    fFlasherDist = lFlasher->GetS("distribution_type");
    info << "VertexGen_Flasher::Init WARNING: Setting intensity distribution type to " + fFlasherDist + "\n";
  }catch (DBNotFoundError &e) {
    info << "VertexGen_Flasher::Init: Using uniform intensity distribution as set in Flasher.ratdb\n";
    fFlasherDist = "uniform";
  }
  fMinPhotonsPerEvent = static_cast<G4int>(static_cast<G4double>(lFlasher->GetI("intensity_min"))/photon_thin);
  debug << "VertexGen_Flasher::Init:  min intensity = " << fMinPhotonsPerEvent << "\n";
  fMaxPhotonsPerEvent = static_cast<G4int>(static_cast<G4double>(lFlasher->GetI("intensity_max"))/photon_thin);
  debug << "VertexGen_Flasher::Init:  max intensity = " << fMaxPhotonsPerEvent << "\n";
  if(fFlasherWavelength=="blackbody"){
    fNphotonNhitParam0 = lFlasher->GetD("blackbody_nphot_nhit_param_0");
    fNphotonNhitParam1 = lFlasher->GetD("blackbody_nphot_nhit_param_1");
    fNphotonNhitParam2 = lFlasher->GetD("blackbody_nphot_nhit_param_2");
  }
  else if(fFlasherWavelength=="flat"){
    fNphotonNhitParam0 = lFlasher->GetD("flat_nphot_nhit_param_0");
    fNphotonNhitParam1 = lFlasher->GetD("flat_nphot_nhit_param_1");
    fNphotonNhitParam2 = lFlasher->GetD("flat_nphot_nhit_param_2");
  }
  else if(fFlasherWavelength=="helium"){
    fNphotonNhitParam0 = lFlasher->GetD("helium_nphot_nhit_param_0");
    fNphotonNhitParam1 = lFlasher->GetD("helium_nphot_nhit_param_1");
    fNphotonNhitParam2 = lFlasher->GetD("helium_nphot_nhit_param_2");
  }
  else if(fFlasherWavelength=="mono" && abs(fMono_wl-500.0)<0.1){
    info << "VertexGen_Flasher::Init: Using 500nm photon distribution\n";
    fNphotonNhitParam0 = lFlasher->GetD("mono500_nphot_nhit_param_0");
    fNphotonNhitParam1 = lFlasher->GetD("mono500_nphot_nhit_param_1");
    fNphotonNhitParam2 = lFlasher->GetD("mono500_nphot_nhit_param_2");
  }
  else if(fFlasherWavelength=="mono" && abs(fMono_wl-350.0)<0.1){
    info << "VertexGen_Flasher::Init: Using 350nm photon distribution\n";
    fNphotonNhitParam0 = lFlasher->GetD("mono350_nphot_nhit_param_0");
    fNphotonNhitParam1 = lFlasher->GetD("mono350_nphot_nhit_param_1");
    fNphotonNhitParam2 = lFlasher->GetD("mono350_nphot_nhit_param_2");
  }
  else{ // for other mono wavelengths or "LED" use mono500 parameters.
    fNphotonNhitParam0 = lFlasher->GetD("mono500_nphot_nhit_param_0");
    fNphotonNhitParam1 = lFlasher->GetD("mono500_nphot_nhit_param_1");
    fNphotonNhitParam2 = lFlasher->GetD("mono500_nphot_nhit_param_2");
  }
  fSNODistParam0 = lFlasher->GetD("sno_dist_param_0");
  fSNODistParam1 = lFlasher->GetD("sno_dist_param_1");
  fSNODistParam2 = lFlasher->GetD("sno_dist_param_2");
  fSNODistParam3 = lFlasher->GetD("sno_dist_param_3");
  fSNODistParam4 = lFlasher->GetD("sno_dist_param_4");
  fSNODistParam5 = lFlasher->GetD("sno_dist_param_5");
  fSNODistParam6 = lFlasher->GetD("sno_dist_param_6");
  fSNODistParam7 = lFlasher->GetD("sno_dist_param_7");
  fSNODistParam8 = lFlasher->GetD("sno_dist_param_8");
  fSNODistParam9 = lFlasher->GetD("sno_dist_param_9");
  fSNODistBound0 = lFlasher->GetD("sno_dist_bound_0");
  fSNODistBound1 = lFlasher->GetD("sno_dist_bound_1");
  fSNODistBound2 = lFlasher->GetD("sno_dist_bound_2");
  fSNODistBound3 = lFlasher->GetD("sno_dist_bound_3");
  fSNOMaxNPhoton = lFlasher->GetI("sno_max_nphot");

  // timing info
  fMono_time = 0;
  try{
    fFlasherTime = lFlasher->GetS("time_dist");
    info << "VertexGen_Flasher::Init WARNING: Setting time distribution to " + fFlasherTime + "\n";
  }catch (DBNotFoundError &e) {
    info << "VertexGen_Flasher::Init: Using fixed time distribution as set in Flasher.ratdb\n";
    fFlasherTime = "DEFAULT";
  }
  if(fFlasherTime=="mono"){
    fMono_time = 1;
    info << "VertexGen_Flasher::Init WARNING: Setting time distribution to mono \n";
  }
  else {
    // firstly get the times and store endpoints
    ftime = lFlasher->GetDArray("dist_time");
    fMaxTime = ftime.back();
    fMinTime = ftime.front();
    // make a check that the values in the dist_time make sense and increase monotonically
    for(unsigned int iv=1;iv<ftime.size(); iv++){
      if(ftime[iv]<=ftime[iv-1]){
        warn << "VertexGen_Flasher::Init: dist_time[" << iv << "] = " << ftime[iv] << " [" << iv-1 << "] = "  << ftime[iv-1] << "!\n";
        Log::Die("VertexGen_Flasher::Init: dist_time must increase monotonically. Please check values.");
      }
    }
    // now get the amplitudes
    ftimeInt = lFlasher->GetDArray("dist_time_intensity");
    // now create a cumulative magnitude vector
    G4double timemagsum = 0;
    ftimeCumMag.clear();
    ftimeCumMag.push_back(0.0);
    detail << " Cumulative time ";
    for(unsigned int i=1;i<ftime.size();i++){
      timemagsum +=(ftimeInt[i]+ftimeInt[i-1])/2.0 * (ftime[i]-ftime[i-1]);
      ftimeCumMag.push_back(timemagsum);
      detail << ftimeCumMag[i] << " ";
    }
    detail << "\n";
  }


  // angular info
  try{
    fFlasherAngle = lFlasher->GetS("angle_dist");
    info << "VertexGen_Flasher::Init: Setting angular distribution to " + fFlasherAngle + "\n";
  }catch (DBNotFoundError &e) {
    info << "VertexGen_Flasher::Init: Using fixed angular distribution as set in Flasher.ratdb\n";
    fFlasherAngle="DEFAULT";
  }
  fIso_angle = 0;
  fZero_angle = 0;
  if(fFlasherAngle=="iso"){
    fIso_angle = 1;
    info << "VertexGen_Flasher::Init: Setting angular distribution to isotropic \n";
  }
  else if(fFlasherAngle=="zero"){
    fZero_angle = 1;
    info << "VertexGen_Flasher::Init: Setting angular distribution to zero (ie no dispersion) \n";
    if(fIso_angle && fZero_angle) Log::Die("VertexGen_Flasher: Beam is set to be both isotropic and Zero dispersion!");
  }
  else{
    string angle_units = lFlasher->GetS("dist_angle_option");
    if(angle_units == "degrees" ){
      fConvertAngle = CLHEP::twopi/360.0;
    }else if(angle_units == "radians"){
      fConvertAngle = 1;
    }else{
      fConvertAngle = 1;
      warn << "VertexGen_Flasher::Init: Unknown angular unit, assuming radians\n";
    }
    // firstly get the angles and store endpoints
    fang = lFlasher->GetDArray("dist_angle");
    fMaxAngle = fang.back();
    fMinAngle = fang.front();
    // make a check that the values in the dist_angle make sense and increase monotonically
    for(unsigned int iv=1;iv<fang.size(); iv++){
      if(fang[iv]<=fang[iv-1]){
        warn << "VertexGen_Flasher::Init: dist_angle[" << iv << "] = " << fang[iv] << " [" << iv-1 << "] = "  << fang[iv-1] << "!\n";
        Log::Die("VertexGen_Flasher::Init: dist_angle must increase monotonically. Please check values.");
      }
    }
    if(fMaxAngle-fMinAngle<-0 || fMaxAngle<0 || fMinAngle<0){
      // problem with angular range
      Log::Die("VertexGen_Flasher: Nonsensical angular range for " + fFlasherAngle);
    }
    // now get the amplitudes
    fangInt = lFlasher->GetDArray("dist_angle_intensity");
    // now create a cumulative magnitude vector and account for solid angle weighting
    G4double angmagsum = 0;
    fangCumMag.clear();
    fangCumMag.push_back(0.0);
    detail << " Cumulative angle ";
    for(unsigned int i=1;i<fang.size();i++){
      angmagsum +=(fangInt[i]+fangInt[i-1])/2.0 * (fang[i]-fang[i-1]);
      fangCumMag.push_back(angmagsum);
      detail << fangCumMag[i] << " ";
    }
    detail << "\n";
  }

  // Should now be all initialised so set flag
  fInitDB = true;
  debug << "VertexGen_Flasher::Init: Init completed." << newline;
}


void VertexGen_Flasher::BeginOfRun() {
  Init();
}

void VertexGen_Flasher::GeneratePrimaryVertex( G4Event *event,
                                               G4ThreeVector& dx,
                                               G4double dt )
{
  if( !fInitDB ) Init();

  G4ThreeVector dir(fDirX,fDirY,fDirZ);
  G4ThreeVector normal = dir.unit();
  G4ThreeVector perp = normal.orthogonal().unit();

  // How many photons?
  // Number of photons will be randomly sampled from a uniform distribution or the SNO flasher spectrum
  G4double fPhotonsPerEvent = 0;
  G4int nphot = 0;
  if(fFlasherDist=="uniform"){
    fPhotonsPerEvent = fMinPhotonsPerEvent+G4UniformRand()*(fMaxPhotonsPerEvent-fMinPhotonsPerEvent);
    if(fPoisson){
      nphot = static_cast<G4int>( CLHEP::RandPoisson::shoot(fPhotonsPerEvent) );
      if(nphot<=0){
        detail << "VertexGen_Flasher::GeneratePrimaryVertex:  Poisson throw of intensity results in zero photons, simulating 1 \n";
        nphot = 1;
      }
    }else{
      nphot = static_cast<G4int>(fPhotonsPerEvent);
    }
  }
  else if(fFlasherDist=="SNO"){
    while(nphot<=0 || nphot>fSNOMaxNPhoton){
      G4double rand_num = CLHEP::RandFlat::shoot(0.0,1.0);
      G4double nhits = 0.0;
      if(rand_num<fSNODistBound0){
        nhits = log((rand_num*fSNODistParam1/exp(fSNODistParam0))+exp(fSNODistParam1*0.0))/fSNODistParam1;
      }
      else if(rand_num<fSNODistBound1){
        nhits = log(((rand_num-fSNODistBound0)*fSNODistParam3/exp(fSNODistParam2))+exp(fSNODistParam3*120.0))/fSNODistParam3;
      }
      else if(rand_num<fSNODistBound2){
        nhits = log(((rand_num-fSNODistBound1)*fSNODistParam5/exp(fSNODistParam4))+exp(fSNODistParam5*200.0))/fSNODistParam5;
      }
      else if(rand_num<fSNODistBound3){
        nhits = log(((rand_num-fSNODistBound2)*fSNODistParam7/exp(fSNODistParam6))+exp(fSNODistParam7*350.0))/fSNODistParam7;
      }
      else{
        nhits = log(((rand_num-fSNODistBound3)*fSNODistParam9/exp(fSNODistParam8))+exp(fSNODistParam9*525.0))/fSNODistParam9;
      }
      fPhotonsPerEvent = exp((nhits+fNphotonNhitParam2)/fNphotonNhitParam0) - fNphotonNhitParam1;
      nphot = static_cast<G4int>(fPhotonsPerEvent);
    }
  }
  else{
    Log::Die("VertexGen_Flasher::GeneratePrimaryVertex: Intensity distribution type not valid!");
  }

  // sample the wavelength up front for writing out to calib branch
  G4double wavelength = 0;
  if(fFlasherWavelength=="mono"){
    wavelength = fMono_wl; // nm
    debug << "VertexGen_Flasher::GeneratePrimaryVertex:  Flasher wavelength fixed = " << wavelength <<"\n";
  }
  else if(fFlasherWavelength=="LED"||fFlasherWavelength=="flat"||fFlasherWavelength=="blackbody"){
    debug << "VertexGen_Flasher::GeneratePrimaryVertex:  Flasher wavelength randomly chosen = ";
    wavelength = SampleRandom(fwave, fwaveInt, fwaveCumMag);
    debug << wavelength <<"\n";
  }
  else{
    debug << "VertexGen_Flasher::GeneratePrimaryVertex:  Flasher wavelength randomly chosen = ";
    G4double rand = G4UniformRand();
    for(unsigned int i=0; i<fwaveCumProb.size()-1; i++){
      if(rand>=fwaveCumProb[i] && rand<=fwaveCumProb[i+1]){
        wavelength = fwave[i];
        break;
      }
    }
    debug << wavelength <<"\n";
  }
  if(wavelength==0){
    Log::Die("VertexGen_Flasher::GeneratePrimaryVertex: Wavelength not set.");
  }
  // and sample the time
  G4double time = 0.;
  if(!fMono_time){
    time = SampleRandom(ftime, ftimeInt, ftimeCumMag);
    //time = (fRandTime->shoot() * (fMaxTime - fMinTime) + fMinTime);
  }
  // Record in extra G4Event information
  EventInfo *exinfo = dynamic_cast<EventInfo*>(event->GetUserInformation());
  DS::Calib *calib = new DS::Calib();
  calib->SetSourceName("Flasher"); // set name as fibre ID now
  calib->SetID(1); // This number is non zero if there is a source, but what ID to assign???
  calib->SetMode(static_cast<G4int>(wavelength)); // sampled wavelength distribution
  calib->SetIntensity(nphot);
  calib->SetDir( TVector3(normal.x(), normal.y(), normal.z()) );
  calib->SetPos( TVector3(dx.x(), dx.y(), dx.z()) );
  calib->SetTime(time);
  exinfo->SetCalib(calib);

  // Add vertices for each photon
  for (G4int photon_i=0; photon_i < nphot; photon_i++) {
    // photon time
    G4double t1 = 0;  // stays as 0 for Mono_time
    if(!fMono_time){
      t1 = SampleRandom(ftime, ftimeInt, ftimeCumMag);
    }
    if(!fMono_wl){ // set above for Mono_wls
      if(fFlasherWavelength=="LED"||fFlasherWavelength=="flat"||fFlasherWavelength=="blackbody"){
        wavelength = SampleRandom(fwave, fwaveInt, fwaveCumMag);
      }
      else if(fFlasherWavelength=="helium"){
        G4double rand = G4UniformRand();
        for(unsigned int j=0; j<fwaveCumProb.size()-1; j++){
          if(rand>=fwaveCumProb[j] && rand<=fwaveCumProb[j+1]){
            wavelength = fwave[j];
            break;
          }
        }
      }
    }

    G4double energy = CLHEP::hbarc * CLHEP::twopi / (wavelength * CLHEP::nm);
    G4double momentum = energy; // GEANT uses momentum in same units as energy

    // photon direction
    G4double theta;
    if (fIso_angle){
      theta = acos(0.9999 * (2.0*G4UniformRand()-1.0));
    }else if(fZero_angle){
      theta = 0;
    }else{
      theta = fConvertAngle*SampleRandom(fang, fangInt, fangCumMag);
    }

    // Set angle back to 0 for first photon (which is parent)
    if ( photon_i == 0 )
      theta = 0.0;

    G4double phi = CLHEP::twopi * G4UniformRand();

    // Calculate momentum
    G4ThreeVector mom(normal);
    mom.rotate(theta, perp); // Rotate away from Flasher direction (polar angle)
    mom.rotate(phi, normal); // Rotate around Flasher direction (phi)
    mom.setMag(momentum);    // Scale to right magnitude

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

    vertex->SetPrimary(particle);

    // Set parent information for first photon
    if( photon_i == 0 ){
      ParentPrimaryVertexInformation *vertinfo = new ParentPrimaryVertexInformation();
      vertinfo->SetVertexParentParticle( particle );
      vertex->SetUserInformation( vertinfo );
    }

    event->AddPrimaryVertex(vertex);
  }

}


G4double VertexGen_Flasher::SampleRandom( std::vector<G4double> mydistx,
                                          std::vector<G4double> mydistraw,
                                          std::vector<G4double> mydistcum )
{
  // Return value sampled from distribution with cumulative distribution mydistcum, raw distribution mydistraw for x points mydistx
  // (assume linear interpolation between specified points in spectrum)
  G4double start = mydistcum.front();
  G4double stop = mydistcum.back();
  // now throw a random number between these limits
  G4double random = start+G4UniformRand()*(stop-start);
  // now loop through cumulative distribution to choose energy
  G4int istep = 0;
  while(mydistcum[istep+1]<random){
    istep++;
  }

  // The desired value is bracketed by istep and istep+1 points.
  // Since the cumulative function is an integral of a linearly
  // interpolated spectrum, it is quadratic between data points.
  // Find x assuming quadratic shape.

  G4double spec_slope = (mydistraw[istep+1] - mydistraw[istep]) / (mydistx[istep+1] - mydistx[istep]);

  G4double x;
  G4double delta_y = random - mydistcum[istep];
  if (fabs(spec_slope) < 1e-6) {
    // ?? not sure I get this step?
    x = mydistx[istep] + delta_y/mydistraw[istep];

  }
  else {
    G4double a = 0.5 * spec_slope;
    G4double b = mydistraw[istep] - spec_slope * mydistx[istep];
    G4double c = -a * mydistx[istep] * mydistx[istep] - b * mydistx[istep] - delta_y;
    x = (-b + sqrt(b*b - 4 * a * c))/(2*a);
  }

  return x;

}


void VertexGen_Flasher::SetDirection( G4ThreeVector &dir )
{
  fDirX = dir.x();
  fDirY = dir.y();
  fDirZ = dir.z();
}


void VertexGen_Flasher::SetState( G4String /*newValues*/ )
{
  warn << "VertexGen_Flasher::SetState error: Cannot set vertex state" << newline;
}


G4String VertexGen_Flasher::GetState()
{
  return G4String("VertexGen_Flasher::GetState: vertex generator fixed");
}

} // namespace RAT