#ifndef __JAANET__JEVTWEIGHTFACTORMUPAGE__ #define __JAANET__JEVTWEIGHTFACTORMUPAGE__ #include "km3net-dataformat/offline/Evt.hh" #include "km3net-dataformat/offline/Trk.hh" #include "JLang/JException.hh" #include "JGeometry2D/JCircle2D.hh" #include "TFormula.h" /** * \author bjung */ namespace JAANET { using JLANG::JValueOutOfRange; using JLANG::JNullPointerException; /** * 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 : 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, //!< 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() : 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; vector muons; vector vars(NUMBER_OF_VARIABLES); for (vector::const_iterator i = evt.mc_trks.cbegin(); i != evt.mc_trks.cend(); ++i) { if (is_muon(*i)) { vars[MEAN_ZENITH_ANGLE] += i->dir.z / i->dir.len(); vars[TOTAL_MUON_ENERGY] += i->E; muons.push_back(getTrack(*i)); } } if (muons.size() > 0) { vars[MUON_MULTIPLICITY] = (Double_t) muons.size(); vars[MEAN_ZENITH_ANGLE] /= vars[MUON_MULTIPLICITY]; if (this->GetNdim() > (int) LATERAL_SPREAD) { JCircle2D circle = JCircle2D(muons.cbegin(), muons.cend()); // smallest enclosing circle vars[LATERAL_SPREAD] = (Double_t) circle.getRadius(); } return this->DoEval(&vars[0]); } else { THROW(JNullPointerException, "JEvtWeightFactorMupage::operator(): No muon for event " << evt.id << '.' << endl); } } }; } #endif