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/JClonable.hh"
#include "JLang/JException.hh"
#include "JLang/JPredicate.hh"
 
#include "JGeometry2D/JCircle2D.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JRotation3D.hh"
#include "JGeometry3D/JTrack3E.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtWeightFactorTFormula.hh"
 
 
/**
 * \author bjung
 */
 
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
 
  using JLANG::JClonable;
  
 
  /**
   * Implementation of reweighting factor for 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 `JEvtWeightFactorMupage::variables`.\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 `JEvtWeightFactorMupage::variables`.
   */
  struct JEvtWeightFactorMupage final:
    public JClonable<JEvtWeightFactorTFormula, JEvtWeightFactorMupage>
  {
    /**
     * Indices of reweighting variables for MUPAGE.
     */
    enum variables {
      MUON_MULTIPLICITY,  //!< Muon multiplicity      
      MEAN_ZENITH_ANGLE,  //!< Average cosine of zenith angle
      TOTAL_MUON_ENERGY,  //!< Muon bundle total energy         [GeV]
      LATERAL_SPREAD,     //!< Muon bundle lateral spread       [m]
      
      NUMBER_OF_VARIABLES //!< Number of reweighting variables; \n
                          //!< N.B.\ This enum value needs to be specified last!
    }; 
 
    
    /**
     * Default constructor.
     */
    JEvtWeightFactorMupage()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  name            name
     * \param  formula         formula
     */
    JEvtWeightFactorMupage(const char* const name,
                           const char* const formula)
    {
      getFormula().SetName(name);
      
      compile(formula);
    }
 
 
    /**
     * Check whether this formula is valid.
     *
     * \return                 true if valid; else false
     */
    bool is_valid() const override final
    {
      const int N = getFormula().GetNpar();
 
      return N >= 0 && N <= (int) NUMBER_OF_VARIABLES;
    }
 
    
    /**
     * Get weighting factor for given event.
     *
     * \param  evt             event
     * \return                 weighting factor
     */
    double getFactor(const Evt& evt) const override final
    {
      using namespace std;
      using namespace JPP;
 
      vector<JTrack3E> muons;
      
      for (vector<Trk>::const_iterator i = evt.mc_trks.cbegin(); i != evt.mc_trks.cend(); ++i) {
        if (is_muon(*i) && is_finalstate(*i)) {
          muons.push_back(getTrack(*i));
        }
      }
 
      if (muons.empty()) {
        THROW(JValueOutOfRange, "JEvtWeightFactorMupage::operator(): No muon for event " << evt.id << '.' << endl);
      }
 
      Double_t  vars[NUMBER_OF_VARIABLES] = { 0.0 };
 
      double    Etot = 0.0;
      JVector3D dir;
 
      for (vector<JTrack3E>::const_iterator i = muons.begin(); i != muons.end(); ++i) {
        Etot += i->getE();
        dir  += i->getDirection();
      }
        
      const JDirection3D D(dir);
      const JRotation3D  R(D);
 
      for (vector<JTrack3E>::iterator i = muons.begin(); i != muons.end(); ++i) {
        i->rotate(R);
      }
 
      vector<Trk>::const_iterator iBundle = find_if(evt.mc_trks.cbegin(), evt.mc_trks.cend(),
                                                    make_predicate(&Trk::status, TRK_ST_MUONBUNDLE));
 
      vars[MUON_MULTIPLICITY] = (double) muons.size();
      vars[TOTAL_MUON_ENERGY] = (iBundle != evt.mc_trks.cend() ? iBundle->E : Etot);
      vars[MEAN_ZENITH_ANGLE] = D.getDZ();
 
      if (this->GetNdim() > (int) LATERAL_SPREAD) {
 
        const JCircle2D circle = JCircle2D(muons.cbegin(), muons.cend()); // smallest enclosing circle
 
        vars[LATERAL_SPREAD] = circle.getRadius();
      }
          
      return getFormula().EvalPar(&vars[0]);
    }
  };
}
 
#endif