JEvtWeightFactorTriggerEfficiencyRatio.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTFACTORTRIGGEREFFICIENCYRATIO__
#define __JAANET__JEVTWEIGHTFACTORTRIGGEREFFICIENCYRATIO__
 
#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&);
  double getVisibleEnergyLeadingLepton(const Trk&, const JCylinder3D&);  
  double getVisibleEnergyLeadingLepton(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 trigger efficiency ratios.
   */
  struct JEvtWeightFactorTriggerEfficiencyRatio final :
    public JClonable<JEvtWeightFactor, JEvtWeightFactorTriggerEfficiencyRatio>
  {
    typedef JRange<double>      JRange_t;
    
    /**
     * Indices of options for calculating trigger efficiency ratios.\n
     * The option names correspond to the variable which is used to compute the trigger efficiency ratio.
     */
    enum JTriggerEfficiencyRatioOption {
      ENERGY_INITIAL_STATE,      //!< Compute trigger efficiency ratio based on scaling of the initial state energy
      ENERGY_NEUTRINO,           //!< Compute trigger efficiency ratio based on scaling of theprimary neutrino energy
      EVIS,                      //!< Compute trigger efficiency ratio based on scaling of the total visible energy
      EVIS_LEADING_LEPTON_CONTR, //!< Compute trigger efficiency ratio based on scaling of the visible energy contribution of the leading lepton
      EVIS_HADRONIC_CONTR,       //!< Compute trigger efficiency 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.
     */
    JEvtWeightFactorTriggerEfficiencyRatio() :
      pTE           (),
      hTE           ("", "default"),
      option        (ENERGY_NEUTRINO),
      ratio         (0.0),
      range         (0.0, 4.0),
      fiducialVolume(getMaximumContainmentVolume()),
      logE          (true)
    {
      pTE = std::make_unique<TH1D>(hTE.getObjectName(), "", 10, range.getLowerLimit(), range.getUpperLimit());
      
      for (Int_t i = 1; i <= pTE->GetNbinsX(); ++i) {
        pTE->SetBinContent(i, 1.0);
      }
 
      pTE->SetDirectory(0);
 
      check_validity();
      configure();
    }
 
    
    /**
     * Constructor.
     *
     * \param  objectID             trigger efficiency 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
     */
    JEvtWeightFactorTriggerEfficiencyRatio(const JRootObjectID&                objectID,
                                           const JTriggerEfficiencyRatioOption option,
                                           const double                        ratio,
                                           const JRange_t&                     range,
                                           const bool                          logE,
                                           const JCylinder3D&                  fiducialVolume) :
      pTE           (),
      hTE           (objectID),
      option        (option),
      ratio         (ratio),
      range         (range),
      fiducialVolume(fiducialVolume),
      logE          (logE)
    {
      check_validity();
      configure();
    }
 
 
    /**
     * Copy constructor.
     *
     * \param  factor               trigger efficiency ratio factor
     */
    JEvtWeightFactorTriggerEfficiencyRatio(const JEvtWeightFactorTriggerEfficiencyRatio& factor) :
      pTE           (),
      hTE           (factor.hTE),
      option        (factor.option),
      ratio         (factor.ratio),
      range         (factor.range),
      fiducialVolume(factor.fiducialVolume),
      logE          (factor.logE)
    {
      pTE.reset((TH1*) factor.pTE->Clone());
      
      pTE->SetDirectory(0);
 
      check_validity();
      configure();
    }
 
 
    /**
     * Configure trigger efficiency ratio factor.
     */
    void configure()
    {
      using namespace std;
      using namespace JPP;
 
      if (hTE.is_valid() && (!pTE || hTE.getObjectName() != string(pTE->GetName()))) {
 
        if (!getFileStatus(hTE.getFilename().c_str())) {
          THROW(JFileOpenException,    "JEvtWeightFactorTriggerEfficiencyRatio::configure(): Could not open file " << hTE.getFilename());
        }
        
        TFile f(hTE.getFilename().c_str(), "read");
        
        TH1* h = (TH1*) f.Get(hTE.getObjectName());
 
        if (h != NULL) {
          pTE.reset(h);
          pTE->SetDirectory(0);
        } else {
          THROW(JNullPointerException, "JEvtWeightFactorTriggerEfficiencyRatio::configure(): Could not extract histogram " << hTE);
        }
        
        f.Close();
      }
 
      switch (option) { 
      case (int) ENERGY_INITIAL_STATE :
        mfp = &JEvtWeightFactorTriggerEfficiencyRatio::getTriggerEfficiencyRatio1;
        break;
      case (int) ENERGY_NEUTRINO:
        mfp = &JEvtWeightFactorTriggerEfficiencyRatio::getTriggerEfficiencyRatio2;
        break;
      case (int) EVIS:
        mfp = &JEvtWeightFactorTriggerEfficiencyRatio::getTriggerEfficiencyRatio3;
        break;
      case (int) EVIS_LEADING_LEPTON_CONTR:
        mfp = &JEvtWeightFactorTriggerEfficiencyRatio::getTriggerEfficiencyRatio4;
        break;
      case (int) EVIS_HADRONIC_CONTR:
        mfp = &JEvtWeightFactorTriggerEfficiencyRatio::getTriggerEfficiencyRatio5;
        break;
      default:
        THROW(JValueOutOfRange, "JEvtWeightFactorRatio::configure(): Invalid option " << option);
      }
    }
 
 
    /**
     * Retrieve trigger efficiency histogram.
     *
     * \return                      trigger efficiency histogram
     */
    const TH1& getHistogram() const
    {
      if (pTE) {
        return *pTE;
      } else {
        THROW(JLANG::JNullPointerException, "JEvtWeightFactorTriggerEfficiencyRatio::getHistogram(): Trigger efficiency 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 trigger efficiency ratio weight factor is valid.
     *
     * \return                      true if valid; else false
     */
    bool is_valid() const override final
    {
      return (pTE                                 &&
              option >= ENERGY_INITIAL_STATE      &&
              option <  NUMBER_OF_OPTIONS         &&
              ratio  >= 0.0                       &&
              range.is_valid()                    &&
              fiducialVolume.getVolume() >= 0.0);
    }
 
 
    /**
     * Perform linear inter- or extrapolation of trigger efficiency for given abscissa value.
     *
     * \param  x                    abscissa value
     * \return                      trigger efficiency
     */
    inline double interpolate(const double x) const
    {
      using namespace std;
 
      const Int_t    Nbins = pTE->GetNbinsX();
 
      //!< Ensure that only given user-range is used for interpolation      
      const Double_t X  = range.constrain(x); 
      
      const Int_t    i0 = max(1, min(Nbins, pTE->FindBin(X)));
      
      const Double_t x0 = pTE->GetBinCenter (i0);
      const Double_t y0 = pTE->GetBinContent(i0);
      
      const Int_t    i1 = max(2, min(Nbins-1, (x > x0 ? i0+1 : i0-1)));
      
      const Double_t x1 = pTE->GetBinCenter (i1);
      const Double_t y1 = pTE->GetBinContent(i1);
      
      return max(0.0, y0 + (x-x0) * (y1-y0)/(x1-x0));
    }
 
 
    /**
     * Get trigger efficiency ratio based on a scaling of the initial state energy.
     * 
     * \param  x1                   first  abscissa value
     * \param  x2                   second abscissa value
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio(const double x1,
                                            const double x2) const
    {
      static const double Xmin = -3.0; //!< Minimum logarithmic energy (= 1 MeV)
      
      double X1 = x1;
      double X2 = x2;
      
      if (logE) {
        X1 = x1 > 0.0 ? log10(x1) : Xmin;
        X2 = x2 > 0.0 ? log10(x2) : Xmin;
      }
      
      const double M1 = interpolate(X1);
      const double M2 = interpolate(X2);
      
      return M2 > 0.0 ? M1/M2 : 1.0;
    }
    
    
    /**
     * Get trigger efficiency ratio based on a scaling of the initial state energy.
     * 
     * \param  E0                   initial state energy [GeV]
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio1(const double E0) const
    {
      return getTriggerEfficiencyRatio(ratio * E0, E0);
    }
    
    
    /**
     * Get trigger efficiency ratio based on a scaling of the initial state energy.
     * 
     * \param  event                event
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio1(const Evt& event) const
    {
      const double E0 = getE0(event);
      
      return getTriggerEfficiencyRatio1(E0);
    }
 
 
    /**
     * Get trigger efficiency ratio based on a scaling of the primary neutrino energy.
     * 
     * \param  Enu                  primary neutrino energy [GeV]
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio2(const double Enu) const
    {
      return getTriggerEfficiencyRatio(ratio * Enu, Enu);
    }
 
 
    /**
     * Get trigger efficiency ratio based on a scaling of the primary neutrino energy.
     * 
     * \param  event                event
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio2(const Evt& event) const
    {
      const Trk& primary = get_neutrino(event);
 
      const double Enu = primary.E;
 
      return getTriggerEfficiencyRatio2(Enu);
    }
 
 
    /**
     * Get trigger efficiency ratio based on a scaling of the total visible energy.
     * 
     * \param  Evis                 total visible energy [GeV]
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio3(const double Evis) const
    {
      return getTriggerEfficiencyRatio(ratio * Evis, Evis);
    }
    
 
    /**
     * Get trigger efficiency ratio based on a scaling of the total visible energy.
     * 
     * \param  event                event
     * \return                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio3(const Evt& event) const
    {
      const double Evis = JSIRENE::getVisibleEnergy(event, fiducialVolume);
 
      return getTriggerEfficiencyRatio3(Evis);
    }
 
 
    /**
     * Get trigger efficiency 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                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio4(const double Evis,
                                             const double EvisLL) const
    {
      return getTriggerEfficiencyRatio(Evis - (1-ratio)*EvisLL, Evis);
    }
    
 
    /**
     * Get trigger efficiency 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                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio4(const Evt& event) const
    {
      const double Evis   = JSIRENE::getVisibleEnergy             (event, fiducialVolume);
      const double EvisLL = JSIRENE::getVisibleEnergyLeadingLepton(event, fiducialVolume);
 
      return getTriggerEfficiencyRatio4(Evis, EvisLL);
    }
 
 
    /**
     * Get trigger efficiency 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                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio5(const double Evis,
                                             const double EvisLL) const
    {
      return getTriggerEfficiencyRatio(ratio*Evis + (1-ratio)*EvisLL, Evis);
    }
    
 
    /**
     * Get trigger efficiency 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                      trigger efficiency ratio
     */
    inline double getTriggerEfficiencyRatio5(const Evt& event) const
    {
      const double Evis   = JSIRENE::getVisibleEnergy             (event, fiducialVolume);
      const double EvisLL = JSIRENE::getVisibleEnergyLeadingLepton(event, fiducialVolume);
 
      return getTriggerEfficiencyRatio5(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 JEvtWeightFactorTriggerEfficiencyRatioHelper(*this, eqpars);
    }
 
 
    /**
     * Get properties of this class.
     *
     * \param  eqpars               equation parameters
     */
    JProperties getProperties(const JEquationParameters& eqpars = JEvtWeightFactor::getEquationParameters()) const override final
    {
      return JEvtWeightFactorTriggerEfficiencyRatioHelper(*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 trigger efficiency ratios **/
    typedef double (JEvtWeightFactorTriggerEfficiencyRatio::*pFunction)(const Evt&) const;
    
 
    /**
     * Auxiliary class for I/O of trigger efficiency ratio factor.
     */
    struct JEvtWeightFactorTriggerEfficiencyRatioHelper :
      public JProperties
    {
      /**
       * Constructor.
       *
       * \param  factor             constant event-weight factor
       * \param  eqpars             equation parameters 
       */
      template<class JEvtWeightFactorTriggerEfficiencyRatio_t>
      JEvtWeightFactorTriggerEfficiencyRatioHelper(JEvtWeightFactorTriggerEfficiencyRatio_t& factor,
                                                   const JEquationParameters&   eqpars) :
        JProperties(eqpars, 1)
      {
        (*this)[JEvtWeightFactor::getTypeKey()] = "trigger efficiency ratio";
 
        this->insert(gmake_property(factor.hTE));
        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 trigger efficiency ratios
    
    std::unique_ptr<TH1> pTE;            //!< Unique pointer to trigger efficiency ratio histogram
    
    JRootObjectID        hTE;            //!< Trigger efficiency histogram OID
    int                  option;         //!< Trigger efficiency ratio option
    double               ratio;          //!< Abscissa ratio
    JRange_t             range;          //!< Applicable range
    JCylinder3D          fiducialVolume; //!< Fiducial volume (for visible energy computation)
    bool                 logE;           //!< Toggle logarithmic abscissa
  };
}
 
#endif