JEvtWeightFactorMeffRatio.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTFACTORMEFFRATIO__
#define __JAANET__JEVTWEIGHTFACTORMEFFRATIO__
 
#include <string>
#include <memory>
 
#include "TH1.h"
#include "TH1D.h"
#include "TFile.h"
 
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"
 
#include "JSystem/JStat.hh"
 
#include "JLang/JClonable.hh"
#include "JLang/JException.hh"
#include "JLang/JStringStream.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JProperties.hh"
 
#include "JTools/JRange.hh"
 
#include "JGeometry3D/JCylinder3D.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtWeightFactor.hh"
 
#include "JGizmo/JRootObjectID.hh"
 
 
/**
 * \author bjung
 */
 
namespace JSIRENE {
 
  using JGEOMETRY3D::JCylinder3D;
 
  /** Forward function declarations. **/
  const JCylinder3D getMaximumContainmentVolume();
  double            getVisibleEnergy           (const Trk&, const JCylinder3D&);  
  double            getVisibleEnergy           (const Evt&, const JCylinder3D&);
}
 
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JLANG::JClonable;
  
  using JTOOLS::JRange;
 
  using JGEOMETRY3D::JCylinder3D;
 
  using JSIRENE::getMaximumContainmentVolume;
 
  using JGIZMO::JRootObjectID;
 
  
  /**
   * Implementation of reweighting factor for effective mass ratios.
   */
  struct JEvtWeightFactorMeffRatio final :
    public JClonable<JEvtWeightFactor, JEvtWeightFactorMeffRatio>
  {
    typedef JRange<double>      JRange_t;
    
    /**
     * Indices of options for calculating effective mass ratios.\n
     * The option names correspond to the variable which is used to compute the effective mass ratio.
     */
    enum JMeffRatioOption {
      ENERGY_INITIAL_STATE,      //!< Compute effective mass ratio based on scaling of the initial state energy
      ENERGY_NEUTRINO,           //!< Compute effective mass ratio based on scaling of theprimary neutrino energy
      EVIS,                      //!< Compute effective mass ratio based on scaling of the total visible energy
      EVIS_LEADING_LEPTON_CONTR, //!< Compute effective mass ratio based on scaling of the visible energy contribution of the leading lepton
      EVIS_HADRONIC_CONTR,       //!< Compute effective mass ratio based on scaling of the visible energy contribution of the hadronic shower
      NUMBER_OF_OPTIONS          //!< N.B.: This enum value needs to be specified last!
    };
 
      
    /**
     * Default constructor.
     */
    JEvtWeightFactorMeffRatio() :
      pMeff         (),
      hMeff         ("", "default"),
      option        (ENERGY_NEUTRINO),
      ratio         (0.0),
      range         (0.0, 4.0),
      fiducialVolume(getMaximumContainmentVolume()),
      logE          (true),
      MeffMin       (1.0e-9)
    {
      pMeff = std::make_unique<TH1D>(hMeff.getObjectName(), "", 10, range.getLowerLimit(), range.getUpperLimit());
      
      for (Int_t i = 1; i <= pMeff->GetNbinsX(); ++i) {
        pMeff->SetBinContent(i, 1.0);
      }
 
      pMeff->SetDirectory(0);
 
      check_validity();
      configure();
    }
 
    
    /**
     * Constructor.
     *
     * \param  objectID             effective mass histogram OID
     * \param  option               option
     * \param  ratio                abscissa ratio
     * \param  range                applicable (logarithmic) energy range [GeV]
     * \param  logE                 toggle logarithmic energies
     * \param  fiducialVolume       fiducial volume
     */
    JEvtWeightFactorMeffRatio(const JRootObjectID&   objectID,
                              const JMeffRatioOption option,
                              const double           ratio,
                              const JRange_t&        range,
                              const bool             logE,
                              const JCylinder3D&     fiducialVolume) :
      pMeff         (),
      hMeff         (objectID),
      option        (option),
      ratio         (ratio),
      range         (range),
      fiducialVolume(fiducialVolume),
      logE          (logE),
      MeffMin       (1.0e-9)      
    {
      check_validity();
      configure();
    }
 
 
    /**
     * Copy constructor.
     *
     * \param  factor               effective mass ratio factor
     */
    JEvtWeightFactorMeffRatio(const JEvtWeightFactorMeffRatio& factor) :
      pMeff         (),
      hMeff         (factor.hMeff),
      option        (factor.option),
      ratio         (factor.ratio),
      range         (factor.range),
      fiducialVolume(factor.fiducialVolume),
      logE          (factor.logE),
      MeffMin       (1.0e-9)      
    {
      pMeff.reset((TH1*) factor.pMeff->Clone());
      
      pMeff->SetDirectory(0);
 
      check_validity();
      configure();
    }
 
 
    /**
     * Configure effective mass ratio factor.
     */
    void configure()
    {
      using namespace std;
      using namespace JPP;
 
      if (hMeff.is_valid() && (!pMeff || hMeff.getObjectName() != string(pMeff->GetName()))) {
 
        if (!getFileStatus(hMeff.getFilename().c_str())) {
          THROW(JFileOpenException,    "JEvtWeightFactorMeffRatio::configure(): Could not open file " << hMeff.getFilename());
        }
        
        TFile f(hMeff.getFilename().c_str(), "read");
        
        TH1* h = (TH1*) f.Get(hMeff.getObjectName());
 
        if (h != NULL) {
          pMeff.reset(h);
          pMeff->SetDirectory(0);
        } else {
          THROW(JNullPointerException, "JEvtWeightFactorMeffRatio::configure(): Could not extract histogram " << hMeff);
        }
        
        f.Close();
      }
 
      switch (option) { 
      case (int) ENERGY_INITIAL_STATE :
        mfp = &JEvtWeightFactorMeffRatio::getMeffRatio1;
        break;
      case (int) ENERGY_NEUTRINO:
        mfp = &JEvtWeightFactorMeffRatio::getMeffRatio2;
        break;
      case (int) EVIS:
        mfp = &JEvtWeightFactorMeffRatio::getMeffRatio3;
        break;
      case (int) EVIS_LEADING_LEPTON_CONTR:
        mfp = &JEvtWeightFactorMeffRatio::getMeffRatio4;
        break;
      case (int) EVIS_HADRONIC_CONTR:
        mfp = &JEvtWeightFactorMeffRatio::getMeffRatio5;
        break;
      default:
        THROW(JValueOutOfRange, "JEvtWeightFactorRatio::configure(): Invalid option " << option);
      }
    }
 
 
    /**
     * Retrieve effective mass histogram.
     *
     * \return                      effective mass histogram
     */
    const TH1& getHistogram() const
    {
      if (pMeff) {
        return *pMeff;
      } else {
        THROW(JLANG::JNullPointerException, "JEvtWeightFactorMeffRatio::getHistogram(): Effective mass histogram is not set.");
      }
    }
    
    
    /**
     * Retrieve fiducial volume.
     *
     * \return                      fiducial volume
     */
    JCylinder3D getFiducialVolume() const
    {
      return fiducialVolume;
    }
 
 
    /**
     * Set fiducial volume.
     *
     * \param  fiducialVolume       fiducial volume
     */
    void setFiducialVolume(const JCylinder3D& fiducialVolume)
    {
      this->fiducialVolume = fiducialVolume;
      
      check_validity();
    }
    
 
    /**
     * Check if this effective mass ratio weight factor is valid.
     *
     * \return                      true if valid; else false
     */
    bool is_valid() const override final
    {
      return (pMeff                               &&
              option >= ENERGY_INITIAL_STATE      &&
              option <  NUMBER_OF_OPTIONS         &&
              ratio  >= 0.0                       &&
              range.is_valid()                    &&
              fiducialVolume.getVolume() >= 0.0);
    }
 
 
    /**
     * Perform linear inter- or extrapolation of effective mass for given abscissa value.
     *
     * \param  x                    abscissa value
     * \return                      effective mass
     */
    inline double interpolate(const double x) const
    {
      using namespace std;
 
      const Int_t    Nbins = pMeff->GetNbinsX();
      
      const Int_t    i0 = max(1, min(Nbins, pMeff->FindBin(x)));
      
      const Double_t x0 = pMeff->GetBinCenter (i0);
      const Double_t y0 = pMeff->GetBinContent(i0);
      
      const Int_t    i1 = max(2, min(Nbins-1, (x > x0 ? i0+1 : i0-1)));
      
      const Double_t x1 = pMeff->GetBinCenter(i1);
      const Double_t y1 = pMeff->GetBinCenter(i1);
 
      return y0 + (x-x0) * (y1-y0)/(x1-x0);
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the initial state energy.
     * 
     * \param  x1                   first  abscissa value
     * \param  x2                   second abscissa value
     * \return                      effective mass ratio
     */
    inline double getMeffRatio(const double x1,
                               const double x2) const
    {
      const double M1 = interpolate(x1);
      const double M2 = interpolate(x2);
 
      const double ratio = M1 / (M2 > 0.0 ? M2 : MeffMin);
      
      return (range(x2) ? ratio : 1.0);
    }
    
 
    /**
     * Get effective mass ratio based on a scaling of the initial state energy.
     * 
     * \param  E0                   initial state energy [GeV]
     * \return                      effective mass ratio
     */
    inline double getMeffRatio1(const double E0) const
    {
      if (logE) {
        const double X = log10(E0);     
        return getMeffRatio(log10(ratio) + X, X);
      } else {
        return getMeffRatio(ratio * E0, E0);
      }
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the initial state energy.
     * 
     * \param  event                event
     * \return                      effective mass ratio
     */
    inline double getMeffRatio1(const Evt& event) const
    {
      const double E0 = getE0(event);
      
      return getMeffRatio1(E0);
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the primary neutrino energy.
     * 
     * \param  Enu                  primary neutrino energy [GeV]
     * \return                      effective mass ratio
     */
    inline double getMeffRatio2(const double Enu) const
    {
      if (logE) {
        const double X = log10(Enu);
        return getMeffRatio(log10(ratio) + X, X);
      } else {
        return getMeffRatio(ratio * Enu, Enu);  
      }
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the primary neutrino energy.
     * 
     * \param  event                event
     * \return                      effective mass ratio
     */
    inline double getMeffRatio2(const Evt& event) const
    {
      const Trk& primary = get_neutrino(event);
 
      const double Enu = primary.E;
 
      return getMeffRatio2(Enu);
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the total visible energy.
     * 
     * \param  Evis                 total visible energy [GeV]
     * \return                      effective mass ratio
     */
    inline double getMeffRatio3(const double Evis) const
    {
      if (logE) {
        const double X = log10(Evis);
        return getMeffRatio(log10(ratio) + X, X);
      } else {
        return getMeffRatio(ratio * Evis, Evis);        
      }
    }
    
 
    /**
     * Get effective mass ratio based on a scaling of the total visible energy.
     * 
     * \param  event                event
     * \return                      effective mass ratio
     */
    inline double getMeffRatio3(const Evt& event) const
    {
      const double Evis = JSIRENE::getVisibleEnergy(event, fiducialVolume);
 
      return getMeffRatio3(Evis);
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the leading leptonic contribution to the total visible energy.
     *
     * This function assumes that the hadronic component in the visible energy corresponds to the total visible energy\n
     * minus the leading leptonic contribution.
     * 
     * \param  Evis                 total visible energy          [GeV]
     * \param  EvisLL               leading lepton visible energy [GeV]
     * \return                      effective mass ratio
     */
    inline double getMeffRatio4(const double Evis,
                                const double EvisLL) const
    {
      if (logE) {
        return getMeffRatio(log10(Evis - (1-ratio) * EvisLL), log10(Evis));
      } else {
        return getMeffRatio(      Evis - (1-ratio) * EvisLL,        Evis);
      }
    }
    
 
    /**
     * Get effective mass ratio based on a scaling of the leading leptonic contribution to the total visible energy.
     *
     * This function assumes that the hadronic component in the visible energy corresponds to the total visible energy\n
     * minus the leading leptonic contribution.
     * 
     * \param  event                event
     * \return                      effective mass ratio
     */
    inline double getMeffRatio4(const Evt& event) const
    {
      const Trk& LL = get_leading_lepton(event);
 
      const double Evis   = JSIRENE::getVisibleEnergy(event, fiducialVolume);
      const double EvisLL = JSIRENE::getVisibleEnergy(LL,    fiducialVolume);
 
      return getMeffRatio4(Evis, EvisLL);
    }
 
 
    /**
     * Get effective mass ratio based on a scaling of the hadronic contribution to the total visible energy.
     *
     * This function assumes that the hadronic component in the visible energy corresponds to the total visible energy\n
     * minus the leading leptonic contribution.
     * 
     * \param  Evis                 total visible energy          [GeV]
     * \param  EvisLL               leading lepton visible energy [GeV]
     * \return                      effective mass ratio
     */
    inline double getMeffRatio5(const double Evis,
                                const double EvisLL) const
    {
      if (logE) {
        return getMeffRatio(log10(ratio * Evis + (1-ratio) * EvisLL), log10(Evis));
      } else {
        return getMeffRatio(      ratio * Evis + (1-ratio) * EvisLL,        Evis);
      }
    }
    
 
    /**
     * Get effective mass ratio based on a scaling of the hadronic contribution to the total visible energy.
     *
     * This function assumes that the hadronic component in the visible energy corresponds to the total visible energy\n
     * minus the leading leptonic contribution.
     * 
     * \param  event                event
     * \return                      effective mass ratio
     */
    inline double getMeffRatio5(const Evt& event) const
    {
      const Trk& LL = get_leading_lepton(event);
 
      const double Evis   = JSIRENE::getVisibleEnergy(event, fiducialVolume);
      const double EvisLL = JSIRENE::getVisibleEnergy(LL,    fiducialVolume);
 
      return getMeffRatio5(Evis, EvisLL);
    }
    
    
    /**
     * Get weighting factor for given event.
     *
     * \param  evt                  event
     * \return                      weighting factor
     */
    double getFactor(const Evt& evt) const override final
    {
      return (this->*mfp)(evt);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars               equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) override final
    {
      return JEvtWeightFactorMeffRatioHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars               equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JEvtWeightFactorMeffRatioHelper(*this, eqpars);      
    }
 
    
    /**
     * Read event-weight factor from input.
     *
     * \param  in                   input stream
     * \return                      input stream
     */
    std::istream& read(std::istream& in) override final
    {
      using namespace std;
      using namespace JPP;
 
      JStringStream is(in);
 
      if (getFileStatus(is.str().c_str())) {
        is.load();
      }
 
      JProperties properties = getProperties();
      is >> properties;
 
      check_validity();
 
      configure();
 
      return in;      
    }
    
    
  private:
    
 
    /** Type definition of pointer to member function for calculating effective mass ratios **/
    typedef double (JEvtWeightFactorMeffRatio::*pFunction)(const Evt&) const;
    
 
    /**
     * Auxiliary class for I/O of effective mass ratio factor.
     */
    struct JEvtWeightFactorMeffRatioHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  factor             constant event-weight factor
       * \param  eqpars             equation parameters 
       */
      template<class JEvtWeightFactorMeffRatio_t>
      JEvtWeightFactorMeffRatioHelper(JEvtWeightFactorMeffRatio_t& factor,
                                     const JEquationParameters&    eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JEvtWeightFactor::getTypeKey()] = "effective mass ratio";
 
        this->insert(gmake_property(factor.hMeff));
        this->insert(gmake_property(factor.option));
        this->insert(gmake_property(factor.ratio));
        this->insert(gmake_property(factor.range));
        this->insert(gmake_property(factor.fiducialVolume));
        this->insert(gmake_property(factor.logE));
      }
    };
 
    
    pFunction            mfp;            //!< Pointer to member method for calculating effective mass ratios
    
    std::unique_ptr<TH1> pMeff;          //!< Unique pointer to effective mass ratio histogram
    
    JRootObjectID        hMeff;          //!< Effective mass histogram OID
    int                  option;         //!< Effective mass ratio option
    double               ratio;          //!< Abscissa ratio
    JRange_t             range;          //!< Applicable range
    JCylinder3D          fiducialVolume; //!< Fiducial volume (for visible energy computation)
    bool                 logE;           //!< Toggle logarithmic abscissa
 
    const double         MeffMin;        //!< Minimum effective mass
  };
}
 
#endif