JEvtWeightFactorMupage.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTFACTORMUPAGE__
#define __JAANET__JEVTWEIGHTFACTORMUPAGE__

#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/Trk.hh"

#include "JLang/JException.hh"

#include "TFormula.h"

namespace {

  using JLANG::JValueOutOfRange;
  using JLANG::JNullPointerException;


  /**
   * Example object for reweighting mupage events according to a specifiable ROOT TFormula.
   *
   * Note: The ROOT TFormula may assume any number of parameters, but should be restricted to\n
   *       the physical variables listed among `JEvtWeightFactorMupageVariables`.\n
   *       These variables may be specified within the ROOT TFormula as 'x[<index>]',\n
   *       where <index> corresponds to the index of the desired physical variable within `JEvtWeightFactorMupageVariables`.
   */
  struct JEvtWeightFactorMupage :
    public TFormula
  {
    /**
     * Indices of reweighting variables for MUPAGE.
     */
    enum variables {
      MUON_MULTIPLICITY,  //!< Muon multiplicity      
      MEAN_ZENITH_ANGLE,  //!< Average cosine of zenith angle
      TOTAL_MUON_ENERGY,  //!< Total muon energy                [GeV]
      
      NUMBER_OF_VARIABLES //!< Number of reweighting variables; \n
                          //!< N.B.: This enum value needs to be specified last!
    }; 

    
    /**
     * Default constructor.
     */
    JEvtWeightFactorMupage() :
      TFormula()
    {}


    /**
     * Constructor.
     *
     * \param  name            name
     * \param  formula         formula
     */
    JEvtWeightFactorMupage(const char* name,
                           const char* formula) :
      TFormula(name, formula)
    {}

    
    /**
     * Get weighting factor for given event.
     *
     * \param  evt             event
     * \return                 weighting factor
     */
    double operator()(const Evt& evt) const
    {
      using namespace std;
      using namespace JPP;

      Double_t vars[NUMBER_OF_VARIABLES];

      for (vector<Trk>::const_iterator i = evt.mc_trks.cbegin(); i != evt.mc_trks.cend(); ++i) {
        
        if (is_muon(*i)) {

          vars[MUON_MULTIPLICITY] += 1;
          vars[MEAN_ZENITH_ANGLE] += i->dir.z / i->dir.len();
          vars[TOTAL_MUON_ENERGY] += i->E;
        }
      }

      if (vars[MUON_MULTIPLICITY] > 0) {

        vars[MEAN_ZENITH_ANGLE] /= vars[MUON_MULTIPLICITY];
          
        return this->DoEval(vars);

      } else {
        
        THROW(JNullPointerException, "JEvtWeightFactorMupage::operator(): No muon for event " << evt.id << '.' << endl);        
      }
    }
  };
}

#endif