SingleGen_module.cc

Go to the documentation of this file.

///
/// \brief   Generator for specific kinematics - particularly single particles
/// \author  brebel@fnal.gov
/// \date
///

// C++ includes.
#include <iostream>
#include <string>
#include <cmath>
#include <memory>
#include <vector>

// ROOT includes
#include "TDatabasePDG.h"
#include "TFile.h"
#include "TH1F.h"

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/RandomNumberGenerator.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvASoft includes
#include "Geometry/Geometry.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nugen/EventGeneratorBase/evgenbase.h"
#include "SummaryData/RunData.h"
#include "SummaryData/SubRunData.h"

#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandGaussQ.h"
#include "CLHEP/Random/RandExponential.h"


namespace evgen {
  class SingleParticle;

  /// A module to check the results from the Monte Carlo generator
  class SingleGen : public art::EDProducer {
  public:
    explicit SingleGen(fhicl::ParameterSet const& pset);
    virtual ~SingleGen();                        

    void produce(art::Event& evt);  
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginRun(art::Run& run);
    void endSubRun(art::SubRun& sr);

  private:

    void Sample(simb::MCTruth &truth);

    /// Depending on configuration of fPMeaning
    /// calculate the momentum of the particle
    double getMomentum(const double& pkf, const double& m) const;
    
    static const int    kUNIF = 0;
    static const int    kGAUS = 1;
    static const int    kHIST = 2;

    std::vector<int>    fPDG;        ///< PDG code of particles to generate
    int                 fPMeaning;   ///< Meaning of P0, fSigmaP and fPDist. 
                                     ///< By default (fP0Meaning=0), P0 and sigmaP is momentum
                                     ///< If fPMeaning=1, then P0 and sigmaP is kinetic energy in GeV
                                     ///< If fPMeaning=2, then P0and sigmaP is total energy in GeV
    int                 fCycle;      ///< cycle number of MC production
    std::vector<double> fP0;         ///< Central momentum (GeV/c) to generate
    std::vector<double> fSigmaP;     ///< Variation in momenta (GeV/c)
    std::vector<double> fSigmaP2;     ///< Variation in momenta (GeV/c)
    std::vector<int>    fPDist;      ///< How to distribute momenta (0=uniform, 1=gaussian)

    std::vector<double> fX0;         ///< Central x position (cm)
    std::vector<double> fY0;         ///< Central y position (cm)
    std::vector<double> fZ0;         ///< Central z position (cm)
    std::vector<double> fT0;         ///< Central t position (ns)
    std::vector<double> fSigmaX;     ///< Variation in x position (cm)
    std::vector<double> fSigmaY;     ///< Variation in y position (cm)
    std::vector<double> fSigmaZ;     ///< Variation in z position (cm)
    std::vector<double> fSigmaT;     ///< Variation in t position (ns)
    std::vector<double> fSigmaX2;     ///< Variation in x position (cm)
    std::vector<double> fSigmaY2;     ///< Variation in y position (cm)
    std::vector<double> fSigmaZ2;     ///< Variation in z position (cm)
    std::vector<double> fSigmaT2;     ///< Variation in t position (ns)
    std::vector<int>    fPosDist;    ///< How to distribute xyz (0=uniform, 1=gaussian)
    std::vector<double> fCosZ0;      ///< Cosine of central angle wrt z-axis
    std::vector<double> fSigmaCosZ;  ///< Size of variation of cosz
    std::vector<double> fSigmaCosZ2;  ///< Size of variation of cosz
    std::vector<int>    fCosZDist;    ///< How to distribute cosz (0=uniform, 1=gaussian)
    std::vector<double> fPhiXY0;     ///< Central angle in the x-y plane (degrees)
    std::vector<double> fSigmaPhiXY; ///< Size of variation in phixy (degrees)
    std::vector<double> fSigmaPhiXY2; ///< Size of variation in phixy (degrees)
    std::vector<int>    fPhiDist;     ///< How to distribute phixy (0=uniform, 1=gaussian)
    std::vector<bool>   fCoverage4pi; ///< Boolean to override other angle args and achieve 4pi coverage. 

    std::string fROOTPath;

    std::vector<float> sums;
    float sum;
    //int bins;
    int random_index;
    float random_nums_mag;
    float random_nums_theta;
    float random_nums_phi;
    float random_nums_posx;
    float random_nums_posy;
    float random_nums_posz;

    TH1F *pmag_mu;
    TH1F *theta_mu;
    TH1F *phi_mu;
    TH1F *pmag_pi;
    TH1F *theta_pi;
    TH1F *phi_pi;
    TH1F *pmag_pr;
    TH1F *theta_pr;
    TH1F *phi_pr;
    TH1F *pmag_ph;
    TH1F *theta_ph;
    TH1F *phi_ph;
    TH1F *pmag_el;
    TH1F *theta_el;
    TH1F *phi_el;

    TH1F *posx_mu;
    TH1F *posx_pi;
    TH1F *posx_pr;
    TH1F *posx_ph;
    TH1F *posx_el;
    TH1F *posy_mu;
    TH1F *posy_pi;
    TH1F *posy_pr;
    TH1F *posy_ph;
    TH1F *posy_el;
    TH1F *posz_mu;
    TH1F *posz_pi;
    TH1F *posz_pr;
    TH1F *posz_ph;
    TH1F *posz_el;

    std::vector<float> sums_pmag_mu;
    std::vector<float> sums_theta_mu;
    std::vector<float> sums_phi_mu;
    std::vector<float> sums_pmag_pi;
    std::vector<float> sums_theta_pi;
    std::vector<float> sums_phi_pi;
    std::vector<float> sums_pmag_pr;
    std::vector<float> sums_theta_pr;
    std::vector<float> sums_phi_pr;
    std::vector<float> sums_pmag_ph;
    std::vector<float> sums_theta_ph;
    std::vector<float> sums_phi_ph;
    std::vector<float> sums_pmag_el;
    std::vector<float> sums_theta_el;
    std::vector<float> sums_phi_el;

    std::vector<float> sums_posx_mu;
    std::vector<float> sums_posx_pi;
    std::vector<float> sums_posx_pr;
    std::vector<float> sums_posx_ph;
    std::vector<float> sums_posx_el;
    std::vector<float> sums_posy_mu;
    std::vector<float> sums_posy_pi;
    std::vector<float> sums_posy_pr;
    std::vector<float> sums_posy_ph;
    std::vector<float> sums_posy_el;
    std::vector<float> sums_posz_mu;
    std::vector<float> sums_posz_pi;
    std::vector<float> sums_posz_pr;
    std::vector<float> sums_posz_ph;
    std::vector<float> sums_posz_el;

    int bins_pmag_mu;
    int bins_theta_mu;
    int bins_phi_mu;
    int bins_pmag_pi;
    int bins_theta_pi;
    int bins_phi_pi;
    int bins_pmag_pr;
    int bins_theta_pr;
    int bins_phi_pr;
    int bins_pmag_ph;
    int bins_theta_ph;
    int bins_phi_ph;
    int bins_pmag_el;
    int bins_theta_el;
    int bins_phi_el;

    int bins_posx_mu;
    int bins_posx_pi;
    int bins_posx_pr;
    int bins_posx_ph;
    int bins_posx_el;
    int bins_posy_mu;
    int bins_posy_pi;
    int bins_posy_pr;
    int bins_posy_ph;
    int bins_posy_el;
    int bins_posz_mu;
    int bins_posz_pi;
    int bins_posz_pr;
    int bins_posz_ph;
    int bins_posz_el;

    // My TFS Hists
    //TH1F *h_flatx, *h_flaty, *h_flatz;
    //TH1F *h_pvecx, *h_pvecy, *h_pvecz;
  };
}

////////////////////////////////////////////////////////////////////////

namespace evgen{

  //____________________________________________________________________________
  SingleGen::SingleGen(fhicl::ParameterSet const& pset)
  : fPMeaning(0)
  , fCycle   (pset.get<int>("Cycle", 0))
  {
    this->reconfigure(pset);

    // get the random number seed, use a random default if not specified    
    // in the configuration file.  
    int seed = pset.get< unsigned int >("Seed", evgb::GetRandomNumberSeed());
  
    createEngine( seed );
    
    produces< std::vector<simb::MCTruth> >();
    produces< sumdata::SubRunData, art::InSubRun >();
    produces< sumdata::RunData,    art::InRun    >();
  }

  //____________________________________________________________________________
  SingleGen::~SingleGen() { }

  //____________________________________________________________________________
  void SingleGen::reconfigure(const fhicl::ParameterSet& pset)
  {
    fPDG        = pset.get< std::vector<int>    >("PDG");

    // Default to momentum
    fPMeaning = pset.get< int >("PMeaning", 0);
    if(fPMeaning == 1) {
      std::cout<<"SingleGen: Using Kinetic energy for the meaning of P0, SigmaP and PDist\n";
    }
    else if(fPMeaning == 2) {
      std::cout<<"SingleGen: Using Total energy for the meaning of P0, SigmaP and PDist\n";
    }

    fP0         = pset.get< std::vector<double> >("P0");
    fSigmaP     = pset.get< std::vector<double> >("SigmaP");
    fSigmaP2    = pset.get< std::vector<double> >("SigmaP2");
    fPDist      = pset.get< std::vector<int>    >("PDist");
    fX0         = pset.get< std::vector<double> >("X0");
    fY0         = pset.get< std::vector<double> >("Y0");
    fZ0         = pset.get< std::vector<double> >("Z0");
    fT0         = pset.get< std::vector<double> >("T0");
    fSigmaX     = pset.get< std::vector<double> >("SigmaX");
    fSigmaY     = pset.get< std::vector<double> >("SigmaY");
    fSigmaZ     = pset.get< std::vector<double> >("SigmaZ");
    fSigmaT     = pset.get< std::vector<double> >("SigmaT");
    fSigmaX2    = pset.get< std::vector<double> >("SigmaX2");
    fSigmaY2    = pset.get< std::vector<double> >("SigmaY2");
    fSigmaZ2    = pset.get< std::vector<double> >("SigmaZ2");
    fSigmaT2    = pset.get< std::vector<double> >("SigmaT2");
    fPosDist    = pset.get< std::vector<int>    >("PosDist");
    fCosZ0      = pset.get< std::vector<double> >("CosZ0");
    fSigmaCosZ  = pset.get< std::vector<double> >("SigmaCosZ");
    fSigmaCosZ2 = pset.get< std::vector<double> >("SigmaCosZ2");
    fCosZDist   = pset.get< std::vector<int>    >("CosZDist");
    fPhiXY0     = pset.get< std::vector<double> >("PhiXY0");
    fSigmaPhiXY = pset.get< std::vector<double> >("SigmaPhiXY");
    fSigmaPhiXY2= pset.get< std::vector<double> >("SigmaPhiXY2");
    fPhiDist    = pset.get< std::vector<int>    >("PhiDist");
    fCoverage4pi= pset.get< std::vector<bool>   >("Coverage4pi");

    fROOTPath   = pset.get< std::string  >("ROOTPath");
  }

  //____________________________________________________________________________
  void SingleGen::beginRun(art::Run& run)
  {
    /*
    art::ServiceHandle<art::TFileService> tfs;
    h_flatx = tfs->make<TH1F>("UnitX", "Unit V Mom;X Number; Count"  , 42, -1.1, 1.1);
    h_flaty = tfs->make<TH1F>("UnitY", "Unit V Mom;Y Number; Count"  , 42, -1.1, 1.1);
    h_flatz = tfs->make<TH1F>("UnitZ", "Unit V Mom;Z Number; Count"  , 42, -1.1, 1.1);

    h_pvecx = tfs->make<TH1F>("VecX", "Mom Vec;X Mom Vec; Count" , 100, -5, 5);
    h_pvecy = tfs->make<TH1F>("VecY", "Mom Vec;Y Mom Vec; Count" , 100, -5, 5);
    h_pvecz = tfs->make<TH1F>("VecZ", "Mom Vec;Z Mom Vec; Count" , 100, -5, 5);
    */
    // grab the geometry object to see what geometry we are using
    art::ServiceHandle<geo::Geometry> geo;
  
    std::unique_ptr<sumdata::RunData> runcol(new sumdata::RunData(geo->DetId(),
                                                                geo->FileBaseName(),
                                                                geo->ExtractGDML()));
  
    run.put(std::move(runcol));

    bool UseROOT = false;
    for (unsigned int i=0; i<fPDG.size(); ++i)
      if (fPDist[i] == kHIST || fPosDist[i] == kHIST || fCosZDist[i] == kHIST || fPhiDist[i] == kHIST)
        UseROOT = true;

    if (UseROOT){
      // Take histos from file
      TFile f(fROOTPath.c_str());
      assert(!f.IsZombie());
      pmag_mu = (TH1F*)f.Get("pmag_mu");
      theta_mu = (TH1F*)f.Get("theta_mu");
      phi_mu = (TH1F*)f.Get("phi_mu");
      pmag_pi = (TH1F*)f.Get("pmag_pi");
      theta_pi = (TH1F*)f.Get("theta_pi");
      phi_pi = (TH1F*)f.Get("phi_pi");
      pmag_pr = (TH1F*)f.Get("pmag_pr");
      theta_pr = (TH1F*)f.Get("theta_pr");
      phi_pr = (TH1F*)f.Get("phi_pr");
      pmag_ph = (TH1F*)f.Get("pmag_ph");
      theta_ph = (TH1F*)f.Get("theta_ph");
      phi_ph = (TH1F*)f.Get("phi_ph");
      pmag_el = (TH1F*)f.Get("pmag_el");
      theta_el = (TH1F*)f.Get("theta_el");
      phi_el = (TH1F*)f.Get("phi_el");
      
      posx_mu = (TH1F*)f.Get("posx_mu");
      posx_pi = (TH1F*)f.Get("posx_pi");
      posx_pr = (TH1F*)f.Get("posx_pr");
      posx_ph = (TH1F*)f.Get("posx_ph");
      posx_el = (TH1F*)f.Get("posx_el");
      posy_mu = (TH1F*)f.Get("posy_mu");
      posy_pi = (TH1F*)f.Get("posy_pi");
      posy_pr = (TH1F*)f.Get("posy_pr");
      posy_ph = (TH1F*)f.Get("posy_ph");
      posy_el = (TH1F*)f.Get("posy_el");
      posz_mu = (TH1F*)f.Get("posz_mu");
      posz_pi = (TH1F*)f.Get("posz_pi");
      posz_pr = (TH1F*)f.Get("posz_pr");
      posz_ph = (TH1F*)f.Get("posz_ph");
      posz_el = (TH1F*)f.Get("posz_el");
      
      assert(pmag_mu);
      assert(theta_mu);
      assert(phi_mu);
      assert(pmag_pi);
      assert(theta_pi);
      assert(phi_pi);
      assert(pmag_pr);
      assert(theta_pr);
      assert(phi_pr);
      assert(pmag_ph);
      assert(theta_ph);
      assert(phi_ph);
      assert(pmag_el);
      assert(theta_el);
      assert(phi_el);
      
      assert(posx_mu);
      assert(posx_pi);
      assert(posx_pr);
      assert(posx_ph);
      assert(posx_el);
      assert(posy_mu);
      assert(posy_pi);
      assert(posy_pr);
      assert(posy_ph);
      assert(posy_el);
      assert(posz_mu);
      assert(posz_pi);
      assert(posz_pr);
      assert(posz_ph);
      assert(posz_el);
      
      bins_pmag_mu = pmag_mu->GetNbinsX();
      bins_theta_mu = theta_mu->GetNbinsX();
      bins_phi_mu = phi_mu->GetNbinsX();
      bins_pmag_pi = pmag_pi->GetNbinsX();
      bins_theta_pi = theta_pi->GetNbinsX();
      bins_phi_pi = phi_pi->GetNbinsX();
      bins_pmag_pr = pmag_pr->GetNbinsX();
      bins_theta_pr = theta_pr->GetNbinsX();
      bins_phi_pr = phi_pr->GetNbinsX();
      bins_pmag_ph = pmag_ph->GetNbinsX();
      bins_theta_ph = theta_ph->GetNbinsX();
      bins_phi_ph = phi_ph->GetNbinsX();
      bins_pmag_el = pmag_el->GetNbinsX();
      bins_theta_el = theta_el->GetNbinsX();
      bins_phi_el = phi_el->GetNbinsX();
      
      bins_posx_mu = posx_mu->GetNbinsX();
      bins_posx_pi = posx_pi->GetNbinsX();
      bins_posx_pr = posx_pr->GetNbinsX();
      bins_posx_ph = posx_ph->GetNbinsX();
      bins_posx_el = posx_el->GetNbinsX();
      bins_posy_mu = posy_mu->GetNbinsX();
      bins_posy_pi = posy_pi->GetNbinsX();
      bins_posy_pr = posy_pr->GetNbinsX();
      bins_posy_ph = posy_ph->GetNbinsX();
      bins_posy_el = posy_el->GetNbinsX();
      bins_posz_mu = posz_mu->GetNbinsX();
      bins_posz_pi = posz_pi->GetNbinsX();
      bins_posz_pr = posz_pr->GetNbinsX();
      bins_posz_ph = posz_ph->GetNbinsX();
      bins_posz_el = posz_el->GetNbinsX();
      
      // Take the CDFs
      sum = 0;
      for (Int_t i=0;i<bins_pmag_mu;i++){
        sum = pmag_mu->GetBinContent(i) + sum;
        sums_pmag_mu.push_back( sum * ( 1/pmag_mu->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_theta_mu;i++){
        sum = theta_mu->GetBinContent(i) + sum;
        sums_theta_mu.push_back( sum * ( 1/theta_mu->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_phi_mu;i++){
        sum = phi_mu->GetBinContent(i) + sum;
        sums_phi_mu.push_back( sum * ( 1/phi_mu->Integral() ) );
      }
      
      sum = 0;
      for (Int_t i=0;i<bins_pmag_pi;i++){
        sum = pmag_pi->GetBinContent(i) + sum;
        sums_pmag_pi.push_back( sum * ( 1/pmag_pi->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_theta_pi;i++){
        sum = theta_pi->GetBinContent(i) + sum;
        sums_theta_pi.push_back( sum * ( 1/theta_pi->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_phi_pi;i++){
        sum = phi_pi->GetBinContent(i) + sum;
        sums_phi_pi.push_back( sum * ( 1/phi_pi->Integral() ) );
      }
      
      sum = 0;
      for (Int_t i=0;i<bins_pmag_pr;i++){
        sum = pmag_pr->GetBinContent(i) + sum;
        sums_pmag_pr.push_back( sum * ( 1/pmag_pr->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_theta_pr;i++){
        sum = theta_pr->GetBinContent(i) + sum;
        sums_theta_pr.push_back( sum * ( 1/theta_pr->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_phi_pr;i++){
        sum = phi_pr->GetBinContent(i) + sum;
        sums_phi_pr.push_back( sum * ( 1/phi_pr->Integral() ) );
      }
      
      sum = 0;
      for (Int_t i=0;i<bins_pmag_ph;i++){
        sum = pmag_ph->GetBinContent(i) + sum;
        sums_pmag_ph.push_back( sum * ( 1/pmag_ph->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_theta_ph;i++){
        sum = theta_ph->GetBinContent(i) + sum;
        sums_theta_ph.push_back( sum * ( 1/theta_ph->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_phi_ph;i++){
        sum = phi_ph->GetBinContent(i) + sum;
        sums_phi_ph.push_back( sum * ( 1/phi_ph->Integral() ) );
      }
      
      sum = 0;
      for (Int_t i=0;i<bins_pmag_el;i++){
        sum = pmag_el->GetBinContent(i) + sum;
        sums_pmag_el.push_back( sum * ( 1/pmag_el->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_theta_el;i++){
        sum = theta_el->GetBinContent(i) + sum;
        sums_theta_el.push_back( sum * ( 1/theta_el->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_phi_el;i++){
        sum = phi_el->GetBinContent(i) + sum;
        sums_phi_el.push_back( sum * ( 1/phi_el->Integral() ) );
      }
      
      sum = 0;
      for (Int_t i=0;i<bins_posx_mu;i++){
        sum = posx_mu->GetBinContent(i) + sum;
        sums_posx_mu.push_back( sum * ( 1/posx_mu->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posx_pi;i++){
        sum = posx_pi->GetBinContent(i) + sum;
      sums_posx_pi.push_back( sum * ( 1/posx_pi->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posx_pr;i++){
        sum = posx_pr->GetBinContent(i) + sum;
        sums_posx_pr.push_back( sum * ( 1/posx_pr->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posx_ph;i++){
        sum = posx_ph->GetBinContent(i) + sum;
        sums_posx_ph.push_back( sum * ( 1/posx_ph->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posx_el;i++){
        sum = posx_el->GetBinContent(i) + sum;
        sums_posx_el.push_back( sum * ( 1/posx_el->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posy_mu;i++){
        sum = posy_mu->GetBinContent(i) + sum;
        sums_posy_mu.push_back( sum * ( 1/posy_mu->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posy_pi;i++){
        sum = posy_pi->GetBinContent(i) + sum;
        sums_posy_pi.push_back( sum * ( 1/posy_pi->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posy_pr;i++){
        sum = posy_pr->GetBinContent(i) + sum;
        sums_posy_pr.push_back( sum * ( 1/posy_pr->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posy_ph;i++){
        sum = posy_ph->GetBinContent(i) + sum;
        sums_posy_ph.push_back( sum * ( 1/posy_ph->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posy_el;i++){
        sum = posy_el->GetBinContent(i) + sum;
        sums_posy_el.push_back( sum * ( 1/posy_el->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posz_mu;i++){
        sum = posz_mu->GetBinContent(i) + sum;
        sums_posz_mu.push_back( sum * ( 1/posz_mu->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posz_pi;i++){
        sum = posz_pi->GetBinContent(i) + sum;
        sums_posz_pi.push_back( sum * ( 1/posz_pi->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posz_pr;i++){
        sum = posz_pr->GetBinContent(i) + sum;
        sums_posz_pr.push_back( sum * ( 1/posz_pr->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posz_ph;i++){
        sum = posz_ph->GetBinContent(i) + sum;
        sums_posz_ph.push_back( sum * ( 1/posz_ph->Integral() ) );
      }
      sum = 0;
      for (Int_t i=0;i<bins_posz_el;i++){
        sum = posz_el->GetBinContent(i) + sum;
        sums_posz_el.push_back( sum * ( 1/posz_el->Integral() ) );
      }
    }
    
    return;
  }

  //____________________________________________________________________________
  void SingleGen::endSubRun(art::SubRun &sr)
  {
    // store the cycle information
    art::ServiceHandle<geo::Geometry> geo;
    std::unique_ptr< sumdata::SubRunData > sd(new sumdata::SubRunData(fCycle));
    sr.put(std::move(sd));

  }
  
  //____________________________________________________________________________
  void SingleGen::produce(art::Event& evt)
  {
  
    ///unique_ptr allows ownership to be transferred to the art::Event after the ut statement
    std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
  
    simb::MCTruth truth;
    truth.SetOrigin(simb::kSingleParticle);
    Sample(truth);
  
    //     std::cout << "put mctruth into the vector" << std::endl;
    truthcol->push_back(truth);
  
    //     std::cout << "add vector to the event " << truthcol->size() << std::endl;
    evt.put(std::move(truthcol));
  
    return;
  }

  //____________________________________________________________________________
  void SingleGen::Sample(simb::MCTruth& mct) 
  {
    // get the random number generator service and make some CLHEP generators
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine &engine = rng->getEngine();
    CLHEP::RandFlat   flat(engine);
    CLHEP::RandGaussQ gauss(engine);

    /*
    // Searchsorted for indices
    //random_nums = gRandom->Uniform(0, 1);
    random_nums = flat.fire(0, 1);
    for (int j=1; j<bins_pmag_mu; j++){
      if (random_nums<sums_pmag_mu.at(j) && random_nums>sums_pmag_mu.at(j-1))
        random_index = j;
    }
    std::cout << "---------------------- HERE " << pmag_mu->GetBinCenter(random_index) << std::endl;
    */
  
    //
    // every event will have one of each particle species in the fPDG array
    //
    for (unsigned int i=0; i<fPDG.size(); ++i) {    
      // Momentum, kinetic energy or full energy
      double pkf = 0.0;
      if (fPDist[i] == kGAUS) {
        pkf = abs( gauss.fire(fP0[i], fSigmaP[i]) );
      }
      else if (fPDist[i] == kUNIF) {
        pkf = flat.fire(fSigmaP[i], fSigmaP2[i]);
      }
      else if (fPDist[i] == kHIST) {
        random_nums_mag = flat.fire(0, 1);
        switch (fPDG[i]) {
        case 13: { 
          for (int j=1; j<bins_pmag_mu; j++)
            if (random_nums_mag<sums_pmag_mu.at(j) && random_nums_mag>sums_pmag_mu.at(j-1))
              random_index = j;
          pkf = pmag_mu->GetBinCenter(random_index); } break;
        case 211: { 
          for (int j=1; j<bins_pmag_pi; j++)
            if (random_nums_mag<sums_pmag_pi.at(j) && random_nums_mag>sums_pmag_pi.at(j-1))
              random_index = j;
          pkf = pmag_pi->GetBinCenter(random_index); } break;
        case 2212: { 
          for (int j=1; j<bins_pmag_pr; j++)
            if (random_nums_mag<sums_pmag_pr.at(j) && random_nums_mag>sums_pmag_pr.at(j-1))
              random_index = j;
          pkf = pmag_pr->GetBinCenter(random_index); } break;
        case 22: { 
          for (int j=1; j<bins_pmag_ph; j++)
            if (random_nums_mag<sums_pmag_ph.at(j) && random_nums_mag>sums_pmag_ph.at(j-1))
              random_index = j;
          pkf = pmag_ph->GetBinCenter(random_index); } break;
        case 11: { 
          for (int j=1; j<bins_pmag_el; j++)
            if (random_nums_mag<sums_pmag_el.at(j) && random_nums_mag>sums_pmag_el.at(j-1))
              random_index = j;
          pkf = pmag_el->GetBinCenter(random_index); } break;
        default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
        }
      }
      else {
        std::cerr << "Not correct config: Mom" << std::endl;
        exit(-1);
      }
      
      // std::cout << "get the mass" << std::endl;
      const TDatabasePDG* pdgt = TDatabasePDG::Instance();
      const TParticlePDG* pdgp = pdgt->GetParticle(fPDG[i]);
      // Mass in GeV
      double    m = 0.0;
      if (pdgp) m = pdgp->Mass();

      // Momentum in GeV/c
      const double p = getMomentum(pkf, m);

      // Choose position
      double x[4];
      if (fPosDist[i] == kGAUS) {
        x[0] = gauss.fire(fX0[i], fSigmaX[i]);
        x[1] = gauss.fire(fY0[i], fSigmaY[i]);
        x[2] = gauss.fire(fZ0[i], fSigmaZ[i]);
        x[3] = gauss.fire(fT0[i], fSigmaT[i]);
      }
      else if (fPosDist[i] == kUNIF) {
        x[0] = flat.fire(fSigmaX[i], fSigmaX2[i]);
        x[1] = flat.fire(fSigmaY[i], fSigmaY2[i]);
        x[2] = flat.fire(fSigmaZ[i], fSigmaZ2[i]);
        x[3] = flat.fire(fSigmaT[i], fSigmaT2[i]);
      }
      else if (fPosDist[i] == kHIST) {
        x[0] = flat.fire(fSigmaX[i], fSigmaX2[i]);
        x[1] = flat.fire(fSigmaY[i], fSigmaY2[i]);
        x[2] = flat.fire(fSigmaZ[i], fSigmaZ2[i]);
        x[3] = flat.fire(fSigmaT[i], fSigmaT2[i]);

        random_nums_posx = flat.fire(0, 1);
        random_nums_posy = flat.fire(0, 1);
        random_nums_posz = flat.fire(0, 1);
        switch (fPDG[i]) {
        case 13: { 
          for (int j=1; j<bins_posx_mu; j++)
            if (random_nums_posx<sums_posx_mu.at(j) && random_nums_posx>sums_posx_mu.at(j-1))
              random_index = j;
          x[0] = posx_mu->GetBinCenter(random_index);     
          for (int j=1; j<bins_posy_mu; j++)
            if (random_nums_posy<sums_posy_mu.at(j) && random_nums_posy>sums_posy_mu.at(j-1))
              random_index = j;
          x[1] = posy_mu->GetBinCenter(random_index);     
          for (int j=1; j<bins_posz_mu; j++)
            if (random_nums_posz<sums_posz_mu.at(j) && random_nums_posz>sums_posz_mu.at(j-1))
              random_index = j;
          x[2] = posz_mu->GetBinCenter(random_index);     
        } break;
        case 211: { 
          for (int j=1; j<bins_posx_pi; j++)
            if (random_nums_posx<sums_posx_pi.at(j) && random_nums_posx>sums_posx_pi.at(j-1))
              random_index = j;
          x[0] = posx_pi->GetBinCenter(random_index);     
          for (int j=1; j<bins_posy_pi; j++)
            if (random_nums_posy<sums_posy_pi.at(j) && random_nums_posy>sums_posy_pi.at(j-1))
              random_index = j;
          x[1] = posy_pi->GetBinCenter(random_index);     
          for (int j=1; j<bins_posz_pi; j++)
            if (random_nums_posz<sums_posz_pi.at(j) && random_nums_posz>sums_posz_pi.at(j-1))
              random_index = j;
          x[2] = posz_pi->GetBinCenter(random_index);     
        } break;
        case 2212: { 
          for (int j=1; j<bins_posx_pr; j++)
            if (random_nums_posx<sums_posx_pr.at(j) && random_nums_posx>sums_posx_pr.at(j-1))
              random_index = j;
          x[0] = posx_pr->GetBinCenter(random_index);     
          for (int j=1; j<bins_posy_pr; j++)
            if (random_nums_posy<sums_posy_pr.at(j) && random_nums_posy>sums_posy_pr.at(j-1))
              random_index = j;
          x[1] = posy_pr->GetBinCenter(random_index);     
          for (int j=1; j<bins_posz_pr; j++)
            if (random_nums_posz<sums_posz_pr.at(j) && random_nums_posz>sums_posz_pr.at(j-1))
              random_index = j;
          x[2] = posz_pr->GetBinCenter(random_index);     
        } break;
        case 22: { 
          for (int j=1; j<bins_posx_ph; j++)
            if (random_nums_posx<sums_posx_ph.at(j) && random_nums_posx>sums_posx_ph.at(j-1))
              random_index = j;
          x[0] = posx_ph->GetBinCenter(random_index);     
          for (int j=1; j<bins_posy_ph; j++)
            if (random_nums_posy<sums_posy_ph.at(j) && random_nums_posy>sums_posy_ph.at(j-1))
              random_index = j;
          x[1] = posy_ph->GetBinCenter(random_index);     
          for (int j=1; j<bins_posz_ph; j++)
            if (random_nums_posz<sums_posz_ph.at(j) && random_nums_posz>sums_posz_ph.at(j-1))
              random_index = j;
          x[2] = posz_ph->GetBinCenter(random_index);     
        } break;
        case 11: { 
          for (int j=1; j<bins_posx_el; j++)
            if (random_nums_posx<sums_posx_el.at(j) && random_nums_posx>sums_posx_el.at(j-1))
              random_index = j;
          x[0] = posx_el->GetBinCenter(random_index);     
          for (int j=1; j<bins_posy_el; j++)
            if (random_nums_posy<sums_posy_el.at(j) && random_nums_posy>sums_posy_el.at(j-1))
              random_index = j;
          x[1] = posy_el->GetBinCenter(random_index);     
          for (int j=1; j<bins_posz_el; j++)
            if (random_nums_posz<sums_posz_el.at(j) && random_nums_posz>sums_posz_el.at(j-1))
              random_index = j;
          x[2] = posz_el->GetBinCenter(random_index);     
        } break;
        default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
        }         
      }
      else {
        std::cerr << "Not correct config: Pos" << std::endl;
        exit(-1);
      }

      // std::cout << "set the position" << std::endl;
      const TLorentzVector pos(x[0], x[1], x[2], x[3]);
    
      // Choose angles

      double       cosz, phi;
      unsigned int itry;
      for (itry=0; itry<1000000; ++itry) {
        if (fCosZDist[i] == kGAUS) {
          cosz = gauss.fire(fCosZ0[i],  fSigmaCosZ[i]);
        }
        else if (fCosZDist[i] == kUNIF) {
          cosz = flat.fire(fSigmaCosZ[i], fSigmaCosZ2[i]);
        }
        else if (fCosZDist[i] == kHIST) {
          //random_nums_theta = flat.fire(0, 1);
          random_nums_theta = gRandom->Uniform(0, 1);
          switch (fPDG[i]) {
          case 13: { 
            for (int j=1; j<bins_theta_mu; j++)
              if (random_nums_theta<sums_theta_mu.at(j) && random_nums_theta>sums_theta_mu.at(j-1))
                random_index = j;
            cosz = theta_mu->GetBinCenter(random_index); 
          } break;
          case 211: { 
            for (int j=1; j<bins_theta_pi; j++)
              if (random_nums_theta<sums_theta_pi.at(j) && random_nums_theta>sums_theta_pi.at(j-1))
                random_index = j;
            cosz = theta_pi->GetBinCenter(random_index); 
          } break;
          case 2212: { 
            for (int j=1; j<bins_theta_pr; j++)
              if (random_nums_theta<sums_theta_pr.at(j) && random_nums_theta>sums_theta_pr.at(j-1))
                random_index = j;
            cosz = theta_pr->GetBinCenter(random_index); 
          } break;
          case 22: { 
            for (int j=1; j<bins_theta_ph; j++)
              if (random_nums_theta<sums_theta_ph.at(j) && random_nums_theta>sums_theta_ph.at(j-1))
                random_index = j;
            cosz = theta_ph->GetBinCenter(random_index); 
          } break;
          case 11: { 
            for (int j=1; j<bins_theta_el; j++)
              if (random_nums_theta<sums_theta_el.at(j) && random_nums_theta>sums_theta_el.at(j-1))
                random_index = j;
            cosz = theta_el->GetBinCenter(random_index); 
          } break;
          default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
          }       
        }
        else {
          std::cerr << "Not correct config: CosZ" << std::endl;
          exit(-1);
        }
        if (cosz>=-1.0 && cosz<=1.0) break;
      }

      if (cosz<-1.0 || cosz>1.0) {
        mf::LogError("SingleGen") << __FILE__ << ":" << __LINE__
        << " Failed to find an acceptable cos(theta_z)"
        << " after many tries.\n"
        << " Please adjust CosZ0 and SigmaCosZ0 in your"
        << " SingleGen.fcl file and rerun";
        abort();
      }

      if (fPhiDist[i] == kGAUS) {
          phi  = gauss.fire(fPhiXY0[i], fSigmaPhiXY[i]);       
      }
      else if(fPhiDist[i] == kUNIF) {
          phi  = flat.fire(fSigmaPhiXY[i], fSigmaPhiXY2[i]);
      }
      else if(fPhiDist[i] == kHIST) {
        random_nums_phi = flat.fire(0, 1);
        switch (fPDG[i]) {
        case 13: { 
          for (int j=1; j<bins_phi_mu; j++)
            if (random_nums_phi<sums_phi_mu.at(j) && random_nums_phi>sums_phi_mu.at(j-1))
              random_index = j;
          phi = phi_mu->GetBinCenter(random_index); 
        } break;
        case 211: {
          for (int j=1; j<bins_phi_pi; j++)
            if (random_nums_phi<sums_phi_pi.at(j) && random_nums_phi>sums_phi_pi.at(j-1))
              random_index = j;
          phi = phi_pi->GetBinCenter(random_index); 
        } break;
        case 2212: {
          for (int j=1; j<bins_phi_pr; j++)
            if (random_nums_phi<sums_phi_pr.at(j) && random_nums_phi>sums_phi_pr.at(j-1))
              random_index = j;
          phi = phi_pr->GetBinCenter(random_index); 
        } break;
        case 22: {
          for (int j=1; j<bins_phi_ph; j++)
            if (random_nums_phi<sums_phi_ph.at(j) && random_nums_phi>sums_phi_ph.at(j-1))
              random_index = j;
          phi = phi_ph->GetBinCenter(random_index); 
        } break;
        case 11: {
          for (int j=1; j<bins_phi_el; j++)
            if (random_nums_phi<sums_phi_el.at(j) && random_nums_phi>sums_phi_el.at(j-1))
              random_index = j;
          phi = phi_el->GetBinCenter(random_index); 
        } break;
        default : { std::cerr << "Not correct config" << std::endl; exit(-1); } break;
        }         
      }
      else {
        std::cerr << "Not correct config: Phi" << std::endl;
        exit(-1);
      }

    
      const double sinz    = sqrt(1.0-cosz*cosz);
      const double sinphi  = sin(phi*M_PI/180.0);
      const double cosphi  = cos(phi*M_PI/180.0);
      //
      // set track id to -i as these are all primary particles and have
      // id <= 0
      //
      const int    trackid = -1*(i+1);

      const std::string primary("primary");
      simb::MCParticle part(trackid, fPDG[i], primary);

      //std::cout << "!!!!!!!!!!!!!! " << cosz << " " << p << std::endl; 
      TLorentzVector pvec;
      if (fCoverage4pi[i]) {
        // Make a unit vector, of gaussian sampled random numbers.
        TVector3 flat_unit_p = TVector3(gauss.fire(), gauss.fire(), gauss.fire()).Unit();
        pvec.SetPxPyPzE(p*flat_unit_p.Px(),
                        p*flat_unit_p.Py(),
                        p*flat_unit_p.Pz(),
                        sqrt(p*p+m*m));
        /*
        h_flatx -> Fill( flat_unit_p.Px() );
        h_flaty -> Fill( flat_unit_p.Py() );
        h_flatz -> Fill( flat_unit_p.Pz() );
        h_pvecx -> Fill( pvec.Px() );
        h_pvecy -> Fill( pvec.Py() );
        h_pvecz -> Fill( pvec.Pz() );
        */
      } else {
        pvec.SetPxPyPzE(cosphi*sinz*p,
                        sinphi*sinz*p,
                        cosz*p,
                        sqrt(p*p+m*m));
      }
      
      part.AddTrajectoryPoint(pos, pvec);
    
      mct.Add(part);
    } //end loop over particles
  }

  //------------------------------------------------------------------------------
  double SingleGen::getMomentum(const double& pkf, const double& m) const{
    
    if     (fPMeaning == 0) return pkf;
    
    double total_energy = 0.0;
    if     (fPMeaning == 1){
      total_energy = pkf + m;
    }
    else if(fPMeaning == 2){
      total_energy = pkf;
    }
    else{
      total_energy = sqrt(pkf*pkf + m*m);
    }
    
    return sqrt(total_energy*total_energy - m*m);
  }
  
  //------------------------------------------------------------------------------
  DEFINE_ART_MODULE(SingleGen)
  
} //end namespace