#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[\]',\n * where \ corresponds to the index of the desired physical variable within `JEvtWeightFactorMupage::variables`. */ struct JEvtWeightFactorMupage final: public JClonable { /** * 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 muons; for (vector::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::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::iterator i = muons.begin(); i != muons.end(); ++i) { i->rotate(R); } vector::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