GenieMultiverseSyst.h

Go to the documentation of this file.

#pragma once

#include <boost/multi_index_container.hpp>
#include <boost/multi_index/member.hpp>
#include <boost/multi_index/ordered_index.hpp>
#include <set>
#include <string>

#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/FwdDeclare.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SystShifts.h"
#include "TGraphAsymmErrors.h"

// multi_index abbreviation
using boost::multi_index_container;
using namespace boost::multi_index;

// Forward declare TH2
class TH2;

namespace ana
{
  /// forward declarations
  class SpectrumLoaderBase;
  
  /// Obtain GENIE systematic band by randomly varying the number of sigmas
  /// pulled out from a normal distribution for all the enabled GENIE tunable parameters.
  /// Originally the json format was chosen as the format of knob configuration file.
  /// Due to this bug: https://svn.boost.org/trac/boost/ticket/6785,
  /// turned to pure text file instead.
  class GenieMultiverseSpectra
  {
  public:
    enum BandOptions
    {
      kNoBand,
      kBandFromMean,
      kBandFromNominal,
      kBandFromMPV,
      kBandFromMedian
    };
    GenieMultiverseSpectra(){};
    GenieMultiverseSpectra(unsigned int        nuniverses,
                           SpectrumLoaderBase& loader,
                           const HistAxis&     axis,
                           const Cut&          cut,
                           const SystShifts&   shift,
                           const Var&          wei = kUnweighted,
                           std::string         config_pathname = "knob_config.txt");

    GenieMultiverseSpectra(unsigned int              nuniverses,
                           SpectrumLoaderBase&       loader,
                           const HistAxis&           axis,
                           const Cut&                cut,
                           const SystShifts&         shift,
                           const Var&                wei,
                           std::vector<const ISyst*> systs);

    /// Multi-universe constructor for spill truch variables.
    GenieMultiverseSpectra(unsigned int           nuniverses,
                           SpectrumLoaderBase&    loader,
                           const NuTruthHistAxis& axis,
                           const NuTruthCut&      cut,
                           const SystShifts&      shift,
                           const NuTruthVar&      wei = kNuTruthUnweighted,
                           std::string            config_pathname = "knob_config.txt");

    /// Multi-universe constructor for spill truch variables.
    GenieMultiverseSpectra(unsigned int                     nuniverses,
                           SpectrumLoaderBase&              loader,
                           const NuTruthHistAxis&           axis,
                           const NuTruthCut&                cut,
                           const SystShifts&                shift,
                           const NuTruthVar&                wei,
                           std::vector<const ISyst*> systs);

    GenieMultiverseSpectra(unsigned int nuniverses, 
                      std::vector<std::unique_ptr<ana::Spectrum>>& spectra,
          bool fromfile = false);
    
    GenieMultiverseSpectra(std::vector<std::unique_ptr<ana::Spectrum>>& spectra);


    ~GenieMultiverseSpectra();
    
    /// return the nominal histogram with the RMS of all other universes as
    /// error band
    TH1* ToTH1() const;
    TH1* ToTH1(double pot) const;
    TGraphAsymmErrors* ToAreaNormalizedTH1(double pot,
                                           int col = -1,
                                           int errCol = -1,
                                           float headroom = 1.3,
                                           bool newaxis = true) const;

    TH2* ToTH2() const;
    TH2* ToTH2(double pot) const;

    /// return the nominal universe
    const Spectrum* Nominal() const;
    /// return all universes
    std::vector<const Spectrum*> AllUniverses() const {return {fSpectra.begin(), fSpectra.end()};}
    /// return RMS extremes
    const Spectrum* UpperRMS() const;
    const Spectrum* LowerRMS() const;
    const Spectrum* Mean() const;
    const Spectrum* UpperSigma(BandOptions opt = kBandFromNominal) const;
    const Spectrum* LowerSigma(BandOptions opt = kBandFromNominal) const;
    
    /// Save the Spectrum structure of all universe to a TFile directory.
    void SaveTo(TDirectory* dir, const std::string& name) const;
    
    /// Load the Spectrum structure of all universe from a TFile directory.
    static std::unique_ptr<GenieMultiverseSpectra> LoadFrom(TDirectory* dir, const std::string& name);
    
    /// Get the number of universes in an object of this class.
    unsigned int GetNUniverses(){return fNUniverses;};

    /// Calculate error band given 2D histograms
    static std::pair<TH2*, TH2*> GetSigmaHistograms(std::vector<TH2*>, BandOptions opt = kBandFromNominal);

    const auto & GetShiftTable() const { return fShiftTable; };
  
  protected:
    ///
    /// spectra of all universes
    ///
    std::vector<Spectrum*> fSpectra;
    
    int fUniverseSeed = 1001; ///< Random seed sequence.
                              ///< Need to have each universe with the same physical parameters.
                              ///< Start from 1001, giving way to the PPFX multi-universe.
    
    mutable BandOptions fCurBandOpt = kNoBand;
    
    ///
    /// function for building knob configuration table
    ///
    void BuildKnobConfigTable(std::string config_pathname);
    
    ///
    /// function to create one universe
    ///
    Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
                             const HistAxis&     axis,
                             const Cut&          cut,
                             const SystShifts&   shift = kNoShift,
                             const Var&          wei = kUnweighted);

    ///
    /// Function to create one universe
    ///
    Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
                             const HistAxis&     axis,
                             const Cut&          cut,
                             std::vector<const ISyst*> systs,
                             const Var&          wei = kUnweighted);

    ///
    /// Function to create one universe for NuTruth variables.
    ///
    Spectrum* CreateUniverse(SpectrumLoaderBase&    loader,
                             const NuTruthHistAxis& axis,
                             const NuTruthCut&      cut,
                             const SystShifts&      shift = kNoShift,
                             const NuTruthVar&      wei = kNuTruthUnweighted);

    ///
    /// Function to create one universe for NuTruth variables.
    ///
    Spectrum* CreateUniverse(SpectrumLoaderBase&   loader,
                             const NuTruthHistAxis& axis,
                             const NuTruthCut&      cut,
                             std::vector<const ISyst*> systs,
                             const NuTruthVar&      wei = kNuTruthUnweighted);

  private:
    /// find the boundaries from RMS/sigma of all universes and store results to
    /// fUpperRMS, fLowerRMS, fUpperSigma, fLowerSigma, etc.
    void FindBandBoundaries(BandOptions opt = kBandFromNominal) const;
    
    /// number of universes to generate
    unsigned int fNUniverses;
    
    bool fLoadFromFile = false;
    /// modes of valid configuration.
    std::set<int> fModes = {0, 1, 2};

    std::vector<std::map<const ISyst*, double>> fShiftTable;

    mutable Spectrum* fUpperRMS       = NULL; ///< spectra for upper rms
    mutable Spectrum* fLowerRMS       = NULL; ///< spectra for lower rms
    mutable Spectrum* fMean           = NULL; ///< spectra for mean
    mutable Spectrum* fUpperSigma     = NULL; ///< spectra for upper sigma
    mutable Spectrum* fLowerSigma     = NULL; ///< spectra for lower sigma
    mutable bool fBoundariesFound     = false;

    /// Find +/- 1 sigma around the mean for a bin of the Spectrum.
    /// A sigma is defined as enclosing 34% of area with the mean.
    static float BinSigma(std::vector<float> events, float nsigmas, float mean = -1);
    
    /// declarations for knob configuration table
    // use boost::multi_index to store a lookup table for knob sampling modes
    struct knob_sampling_mode
    {
      int id;
      std::string name;
      int mode;
      
      knob_sampling_mode(int id_,std::string name_,int mode_):id(id_),name(name_),mode(mode_){}
    };
    
    /* tags for accessing the corresponding indices of employee_set */
    struct id{};
    struct name{};

    typedef multi_index_container<
      knob_sampling_mode,
      indexed_by<
        ordered_unique<
          tag<id>,  BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,int,id)>,
        ordered_unique<
          tag<name>,BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,std::string,name)> >
    > knob_sampling_mode_set;
  protected:
    ///
    /// knob configuration table
    ///
    knob_sampling_mode_set fKnobConfigTable;

  };
  /// typedef for backward compatibility
  typedef GenieMultiverseSpectra MultiverseSpectra;
  
  class GenieMultiverseParameters
  {
  public:
    /// This class generates the parameters of the multi-universes. No more, no less.
    GenieMultiverseParameters(unsigned int nuniverses,
                              std::string  config_pathname = "knob_config.txt");
    GenieMultiverseParameters(unsigned int nuniverses,
                              std::vector<const ISyst*> systs);
    /// Return the tables of knob shift values.
    std::vector<std::map<const ISyst*, double>>        ShiftTables() {return fShiftTables;};
    std::vector<std::map<const ISyst*, double>> NuTruthShiftTables() {return fNuTruthShiftTables;};
    std::vector<SystShifts>                            GetSystShifts() {return fSystShifts;};
  private:
    int          fUniverseSeed; ///< Random seed sequence.
                                ///< Need to have each universe with the same physical parameters.
                                ///< Start from 1001, giving way to the PPFX multi-universe.
    //int          fUniverseSeedTruth;
    bool fIsNuTruth = false;
    
    std::vector<std::map<const ISyst*, double>>        fShiftTables; ///< Container of shift tables for all universes.
    std::vector<std::map<const ISyst*, double>> fNuTruthShiftTables; ///< Container of neutrino truth shift tables for all universes.
    std::vector<SystShifts>                            fSystShifts;
    
    void BuildKnobConfigTable(std::string config_pathname); ///< Function for building knob configuration table.
    
    struct knob_sampling_mode
    {
      int id;
      std::string name;
      int mode;
      
      knob_sampling_mode(int id_,std::string name_,int mode_):id(id_),name(name_),mode(mode_){}
    };
    /* tags for accessing the corresponding indices of employee_set */
    struct id{};
    struct name{};

    typedef multi_index_container<
      knob_sampling_mode,
      indexed_by<
        ordered_unique<
          tag<id>,  BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,int,id)>,
        ordered_unique<
          tag<name>,BOOST_MULTI_INDEX_MEMBER(knob_sampling_mode,std::string,name)> >
    > knob_sampling_mode_set;
    /// knob configuration table
    knob_sampling_mode_set fKnobConfigTable;
  };
  
  ///
  /// Class for shape-only systematic effect.
  ///
  class GenieMultiverseNormalizedSpectra : public GenieMultiverseSpectra
  {
  public:
    GenieMultiverseNormalizedSpectra(unsigned int        nuniverses,
                                     SpectrumLoaderBase& loader,
                                     const HistAxis&     axis,
                                     const Cut&          cut,
                                     const SystShifts&   shift,
                                     const Var&          wei = kUnweighted,
                                     std::string         config_pathname = "knob_config.txt");

    GenieMultiverseNormalizedSpectra(unsigned int              nuniverses,
                                     SpectrumLoaderBase&       loader,
                                     const HistAxis&           axis,
                                     const Cut&                cut,
                                     const SystShifts&         shift,
                                     const Var&                wei,
                                     std::vector<const ISyst*> systs);
    
    GenieMultiverseNormalizedSpectra(
      unsigned int nuniverses, 
      std::vector<std::unique_ptr<ana::Spectrum>> & spectra,
      bool fromfile = false) : GenieMultiverseSpectra(nuniverses, spectra, fromfile){};
    
    const Spectrum* UpperSigma(BandOptions opt = kBandFromNominal);
    const Spectrum* LowerSigma(BandOptions opt = kBandFromNominal);

    /// Save the Spectrum structure of all universe to a TFile directory.
    void SaveTo(TDirectory* dir, const std::string& name);
    
    /// Load the Spectrum structure of all universe from a TFile directory.
    static std::unique_ptr<GenieMultiverseNormalizedSpectra> LoadFrom(TDirectory* dir, const std::string& name);
    
  private:
    double fNominalCount;
    
    /// A flag signaling whether the spectra are normalized already.
    bool fSpectraNormalized = false;
    
    Spectrum* CreateUniverse(SpectrumLoaderBase& loader,
                             const HistAxis&     axis,
                             const Cut&          cut,
                             const SystShifts&   shift = kNoShift,
                             const Var&          wei = kUnweighted);

    /// Normalize spectra. It has to be done AFTER SpectrumLoader::Go() is called.
    void NormalizeSpectra();
  };
}