GSNovaModel.cxx

Go to the documentation of this file.

// C++ includes
#include <algorithm>
#include <exception>

// ROOT includes
#include "TFile.h"
#include "TH2.h"

// GENIE includes
#include "GENIE/Framework/Messenger/Messenger.h"
#include "GenieSNova/src/GSNovaModel.h"

using namespace genie;
using namespace genie::supernova;


//.......................................................................
GSNovaModel::GSNovaModel(std::string filename)
{
    this->ReadModel(filename);
}


//.......................................................................
TH2* ScaleHist(TH2* hist)
{
    // The flux is normalized to a fluence through 1cm2 surface at 1kPc distance
    const double kkpcPerCm         = 3.086e21;        /// kPc to cm
    const double kLumiMag          = 1e50;            /// Luminosity magnitude
    const double kSphereSolidAngle = 4 * TMath::Pi(); /// in str
    const double kWindowArea       = 1;               /// in cm2

    const double norm = (kLumiMag / (kkpcPerCm * kkpcPerCm)) * kWindowArea / kSphereSolidAngle;

    hist->Scale(norm);

    return hist;
}


//.......................................................................
void GSNovaModel::ReadModel(std::string filename)
{
    TFile f(filename.data(), "READ");
    if (!f.IsOpen()) {
        throw std::runtime_error("Error opening SN model file " + filename);
    }

    std::map<int, std::string> histMap = {
                                           { 12, "fluence_e" },
                                           {-12, "fluence_ae"},
                                           { 14, "fluence_x" },
                                           {-14, "fluence_ax"}
                                         };

    int histPDG;
    std::string histName;
    for(auto& histMapEntry: histMap){
        histPDG  = histMapEntry.first;
        histName = histMapEntry.second;

        TH2* hFlux = (TH2*)f.Get(histName.data());
        if (!hFlux) {
            throw std::runtime_error("Failed to read histogram \'"
                                     + histName
                                     + "\' from file \'"
                                     + filename + "\'");
         }

        this->AddGenerator(GFluxGenerator(histPDG, ScaleHist(hFlux)));
    }

    this->SetProbScale(1.0);
    f.Close();
}


//.......................................................................
const genie::PDGCodeList& GSNovaModel::PdgList()
{
    return fPdgList;
}


//.......................................................................
void GSNovaModel::AddGenerator(GFluxGenerator generator)
{
    // Add new generator to the model
    LOG("GSNovaModel", pNOTICE) << "Add generator for pdg=" << generator.PDG()
                                << ", integral="            << generator.GetIntegral();

    //update generators list
    fGeneratorsDead.push_back(generator);

    //update integral
    fIntegral += generator.GetIntegral();

    //update pdg list
    fPdgList.push_back(generator.PDG());

    //update limits
    fLimsE.Update(generator.LimsE());
    fLimsT.Update(generator.LimsT());
}


//.......................................................................
void GSNovaModel::GenerateNext(GenParticle& particle)
{
    // Generate next particle and fill it to the argument
    if (this->End()) return;

    // (Randomly) select from which generator we will take the input
    auto itGenerator = this->SelectGenerator();
    assert (itGenerator != fGeneratorsActive.end());

    GFluxGenerator& gen =* itGenerator;
    LOG("GSNovaModel", pDEBUG) << "Using generator [" << gen.PDG() << "],"
                               << "t="  << gen.T()
                               << ",e=" << gen.E();

    // Fill the generated calues (T,E) to the part
    gen.FillParticle(particle);

    //Generate a new particle in this generator
    gen.GenerateNext();
    if (gen.End()) {
        LOG("GSNovaModel", pNOTICE) << "Finishing generator #" << gen.PDG();

        // Move generator to Dead list
        fGeneratorsDead.splice(fGeneratorsDead.begin(), fGeneratorsActive, itGenerator);
        LOG("GSNovaModel", pNOTICE) << fGeneratorsActive.size() << " left";
    }

    LOG("GSNovaModel", pDEBUG) << "Done";
}


//.......................................................................
bool GSNovaModel::End()
{
    return fGeneratorsActive.empty();
}


//.......................................................................
void GSNovaModel::Reload()
{
    // Move all generators to active state
    fGeneratorsActive.splice(fGeneratorsActive.begin(), fGeneratorsDead);

    // Reload all generators
    for (auto &gen: fGeneratorsActive) {
        LOG("GSNovaModel", pDEBUG) << "Reloading generator " << gen.PDG();
        gen.Reload();
        gen.GenerateNext();
    }
}


//.......................................................................
void GSNovaModel::SetEmin(double Emin)
{
    for (auto& gen: fGeneratorsActive) gen.SetEmin(Emin);
}


//.......................................................................
void GSNovaModel::SetProbScale(double scale)
{
    LOG("GSNovaModel", pNOTICE) << "Set Probability scale:" << scale;

    for (auto& gen: fGeneratorsDead)   gen.SetProbScale(scale);
    for (auto& gen: fGeneratorsActive) gen.SetProbScale(scale);

    this->Reload();
}


//.......................................................................
bool compare(const GFluxGenerator& g1, const GFluxGenerator& g2)
{
    return g1.T() < g2.T();
}


//.......................................................................
const std::list<GFluxGenerator> GSNovaModel::ActiveGenerators() {
    return fGeneratorsActive;
}


//.......................................................................
std::list<GFluxGenerator>::iterator GSNovaModel::SelectGenerator()
{
    // For timed distribution: select the earliest particle in time
    LOG("GSNovaModel", pDEBUG) << "Select from " << fGeneratorsActive.size() << " active generators";

    return std::min_element(fGeneratorsActive.begin(), fGeneratorsActive.end(), compare);
}