////////////////////////////////////////////////////////////////////////
/// \class RAT::Gen_Muon
///
/// \brief   Top level event generator of cosmic muons using energy analytical energy spectra from: https://arxiv.org/pdf/astro-ph/0512125
///
/// \author Katharine Dixon ,katharine.dixon@kcl.ac.uk>
///
/// HISTORY:
///   08/05/2024: first version - Katharine Dixon
///
/// \details This generator simulates cosmic muons. The code has been developed from the flux generator in RAT.
/// The difference between this generator and the flux generator is that this generator uses an analytical form for sampling the muon's energy and direction. The sampling method used is an inverse CDF approach.
/// The CDF for the direction sampling is found in MUON.ratdb. This table also defines the constants used in the energy spectrum formula, taken from the paper given above. The constants are:
/// b: 0.4 per km.w.e
/// eps: 693 GeV
/// gam: 3.77
/// h0: 6.011 km.w.e
/// The command to use this generator is: /generator/add muon
/// This generator implements a standard time generator, the default option is the poisson time generator. To change this use the command: /generator/muon/settimegen
/// This generator allows the user to change the radius of the sphere on whose surface the vertices are created. To change this use the command: /generator/muon/setradius double . The default value is the psup h2o outer radius 8500.0 mm.
///
////////////////////////////////////////////////////////////////////////

#ifndef __RAT_Gen_Muon__
#define __RAT_Gen_Muon__

#include <RAT/GLG4Gen.hh>
#include <G4ThreeVector.hh>
#include <vector>
#include <CLHEP/Random/RandGeneral.h>
#include <RAT/DB.hh>
#include <stdexcept>

class G4ParticleDefinition;

namespace RAT {

class Gen_MuonMessenger;

class Gen_Muon : public GLG4Gen {

public:
  class MuonException : public std::runtime_error{
    public:
      MuonException(const std::string& errorStr_) : runtime_error(errorStr_) {}
  };

  Gen_Muon();
  virtual ~Gen_Muon();
  virtual void GenerateEvent(G4Event *event);
  virtual void ResetTime(double offset=0.0);
  virtual bool IsRepeatable() const { return true; };

  virtual void SetState(G4String state);
  virtual G4String GetState() const;

  //! specify time generator, default is poisson
  virtual void SetTimeGen(G4String state);

  //! position setter.
  virtual void SetGeneratorPosition(const double x0_,const double y0_,const double z0_);
  virtual void SetGeneratorPosition(const double theta_,const double phi_);
  virtual void SetGeneratorDirection(const double theta_,const double phi_);

  //! overall time state - time from one event to next
  virtual void SetTimeState(G4String state);
  virtual G4String GetTimeState() const;

  //! specify radius of spherical surface for vertex production, default is 8500.0 mm
  virtual void SetRadius(double radius);
  virtual double GetRadius() const;

  // Sample cos(theta)
  void SampleCostheta(double &CosTheta);

  //! Sampling energy for a given cos(theta)
  void SampleE(double &E, const double CosTheta);

  void GenVertex(double long_radius, G4ThreeVector direction, double &x_vertex, double &y_vertex, double &z_vertex);

  virtual void BeginOfRun();

protected:
  DBLinkPtr fLMuon;

  G4String fStateStr;   // generator initialization, specified by user

  G4ParticleDefinition* fPDef;    // particle definition

  Gen_MuonMessenger* fMessenger;  // allows to change parameters via command line

  double fR1;

  bool fQpositionXYZ;
  bool fQpositionThetaPhi;
  bool fQDirectionThetaPhi;
  double fTheta;
  double fPhi;

  double fTheta_dir, fPhi_dir;

  double fPosx;
  double fPosy;
  double fPosz;

  std::vector<double> fCDFcostheta;
  std::vector<double> fCosThetaBins;

  // constants for energy spectrum
  double fBMuon;
  double fEpsMuon;
  double fGamMuon;
  double fH0Muon;

};

} //namespace RAT


#endif