GFluxGenerator.cxx

Go to the documentation of this file.

// ROOT includes
#include "TH2.h"

// GENIE includes
#include "GENIE/Framework/Messenger/Messenger.h"
#include "GENIE/Framework/Conventions/Units.h"

// NOvASoft includes
#include "GenieSNova/src/GFluxGenerator.h"

using namespace genie::supernova;


//.......................................................................
GFluxGenerator::GFluxGenerator(int pdg, TH2* hist) : fPdg(pdg), fT(0), fE(0)
{
    // Set the limits
    fLimsT.Set(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax());
    fLimsE.Set(hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());

    // Set generator integral
    fIntegral = hist->Integral("width");

    // Init the time sampler
    TH1D* histT = hist->ProjectionX(Form("HistT_%d", pdg));
    fTSampler.Init(histT);
    int NBins = hist->GetNbinsX();

    // Init the energy sampler (per each time bin)
    fESamplers.resize(NBins);
    for (int nbin = 0; nbin < NBins; ++nbin) {
        TH1D* h = hist->ProjectionY(Form("HistE_%d_%d", pdg, nbin), nbin, nbin);
        fESamplers[nbin].Init(h);
    }
}


//.......................................................................
void GFluxGenerator::GenerateNext()
{
    fE = 0;
    do {
        // Generate next event time
        fT = fTSampler.GenerateNext();
        if (this->End()) return;

        // Get current time bin
        int bin = fTSampler.GetBin();
        LOG("GFluxGenerator", pDEBUG) << "[" << fPdg
                                      << "] T=" << fT
                                      << " (bin="  << bin  <<")";

        // Generate event energy using current bin distribution
        fE = fESamplers[bin].GenerateNext();
        LOG("GFluxGenerator", pDEBUG) << "[" << fPdg
                                      << "] E=" << fE
                                      << " vs Emin=" << fEmin;

    } while (fE <= fEmin);

    fE *= genie::units::MeV;

    LOG("GFluxGenerator", pDEBUG) << "[" << fPdg
                                  << "]: generated T=" << fT
                                  << "; E=" << fE;
}


//.......................................................................
void GFluxGenerator::FillParticle(GenParticle& particle)
{
    particle.SetPDG(fPdg);
    particle.SetE(fE);
    particle.SetT(fT);
}


//.......................................................................
void GFluxGenerator::Reload()
{
    fTSampler.Reload();
}


//.......................................................................
void GFluxGenerator::SetProbScale(double scale)
{
    fTSampler.SetProbScale(scale);
}


//.......................................................................
void GFluxGenerator::SetEmin(double Emin)
{
    fEmin = Emin;
}


//.......................................................................
bool GFluxGenerator::End()
{
    return fTSampler.End();
}


//.......................................................................
int GFluxGenerator::PDG() const
{
    return fPdg;
}


//.......................................................................
double GFluxGenerator::T() const
{
    return fT;
}


//.......................................................................
double GFluxGenerator::E() const
{
    return fE;
}


//.......................................................................
const OrderedSampler& GFluxGenerator::GetSamplerT() const
{
    return fTSampler;
}


//.......................................................................
const RandomSampler& GFluxGenerator::GetSamplerE(size_t i) const
{
    return fESamplers[i];
}