////////////////////////////////////////////////////////////////////
#ifndef __RAT_VertexGen_IBD__
#define __RAT_VertexGen_IBD__


/**
	\class RAT::VertexGen_IBD
	vertex generator that can generate the products of a inverse
    beta-decay reaction from a reactor anti-neutrino or geoneutrino. The
    direction of the neutrino is supplied, and the energy and angle of the
    produced positron and neutron are drawn from the distribution
    produced by the differential cross-section and a reactor anti-neutrino or
    geoneutrino energy spectrum.

    17-10-2016 - N. Barros - Added a couple of updates to follow with the
               overload of the reactor generator
    14-02-2019 - S. Nae - Add the option to select the spectrum by the IBD table
               index field to accomodate for the introduction of geoneutrino
               spectra
*/

#include <RAT/NeutrinoVertexGen.hh>
#ifdef RATDEBUG
#include <RAT/Profiler.hh>
#endif
#include <G4LorentzVector.hh>
#include <G4ThreeVector.hh>

#include <vector>

class TH1D;
class G4ParticleDefinition;

namespace RAT {

class IBDCrossSec;

class VertexGen_IBD : public NeutrinoVertexGen {
public:
  VertexGen_IBD(const char *arg_dbname="IBD");
  virtual ~VertexGen_IBD();
  virtual void GeneratePrimaryVertex( G4Event *argEvent,
				      G4ThreeVector &dx,
				      G4double dt);

  // This one only makes sense for neutrino sources such as ES and IBD
  virtual void GeneratePrimaryVertex( G4Event *argEvent,
                G4ThreeVector &dx,
                G4double dt,
                G4ThreeVector &nu_dir);

  // generates a primary vertex with given particle type, direction, energy,
  // and consistent polarization.
  virtual void SetState( G4String newValues );
  // format: dir_x dir_y dir_z spec_index
  // If dir_x==dir_y==dir_z==0, the directions are isotropic.
  virtual G4String GetState();
  // returns current state formatted as above

  virtual void BeginOfRun();
  virtual void EndOfRun();

  void SetSpectrum();
  void SetSpectrum(TH1D *flux);

  void GenEvent(const G4ThreeVector &nu_dir,
      G4LorentzVector &neutrino,
      G4LorentzVector &positron,
      G4LorentzVector &neutron);

  void GenInteraction(double &E, double &CosThetaLab);


  virtual G4double GetRatePerTarget();



private:
  const G4ParticleDefinition *fANu, *fNeutron, *fEplus;
  G4ThreeVector         fNuDir;
  double                fCrossSecMax;
  double                fFluxMax;
  IBDCrossSec           *fCrossSec;

  static const double   sDeltaM;

  TH1D *fFlux;
  std::vector<double>   fSpectrumX;
  std::vector<double>   fSpectrumY;
  G4String              fSpectrumIndex;

  bool                  fSpectrumLoaded;
  bool                  fStandaloneMode;
#ifdef RATDEBUG
  Profiler              fProfiler;
#endif
};


} // namespace RAT

#endif