#ifndef __JAANET__JEVTWEIGHTTOOLKIT__ #define __JAANET__JEVTWEIGHTTOOLKIT__ #include #include "km3net-dataformat/definitions/weightlist.hh" #include "km3net-dataformat/definitions/w2list_gseagen.hh" #include "km3net-dataformat/definitions/w2list_km3buu.hh" #include "km3net-dataformat/definitions/w2list_genhen.hh" #include "km3net-dataformat/offline/Head.hh" #include "km3net-dataformat/offline/Evt.hh" #include "JLang/JType.hh" #include "JLang/JException.hh" #include "JPhysics/JConstants.hh" #include "JAAnet/JHead.hh" #include "JAAnet/JHeadToolkit.hh" #include "JAAnet/JEvtWeight.hh" #include "JAAnet/JEvtWeightDAQ.hh" #include "JAAnet/JEvtWeightGenhen.hh" #include "JAAnet/JEvtWeightMupage.hh" #include "JAAnet/JEvtWeightCorsika.hh" #include "JAAnet/JEvtWeightGSeaGen.hh" #include "JAAnet/JEvtWeightKM3BUU.hh" #include "JAAnet/JEvtWeightMiscellaneous.hh" #include "JAAnet/JParticleTypes.hh" #include "JAAnet/JAAnetToolkit.hh" /** * \author bjung, mdejong */ namespace JAANET {} namespace JPP { using namespace JAANET; } namespace JAANET { using JLANG::JType; using JLANG::JCastException; using JEEP::JProperties; /** * Look-up table for event weighters. */ struct JEvtWeighter : public std::vector > { /** * Constructor */ JEvtWeighter() { this->emplace_back(new JEvtWeightGSeaGen()); this->emplace_back(new JEvtWeightKM3BUU()); this->emplace_back(new JEvtWeightCorsika()); this->emplace_back(new JEvtWeightMupage()); this->emplace_back(new JEvtWeightGenhen()); this->emplace_back(new JEvtWeightDAQ()); // must be after Monte Carlo weighters this->emplace_back(new JEvtWeightMiscellaneous()); } /** * Get event weighter corresponding to given header. * * \param header header * \return event weighter */ const JEvtWeight& operator()(const JHead& header) const { using namespace std; using namespace JPP; for (const_iterator i = this->begin(); i != this->end(); ++i) { if ((*i)->check(header)) { return *(*i); } } THROW(JValueOutOfRange, "JEvtWeighter::operator(): No event weighter found for given header:" << endl << header); } }; /** * Get volume of given event according given weighter. * * \param type type * \param evt event * \return volume [m^3 Gev sr] */ inline double getVolume(const JType& type, const Evt& evt) { return (evt.w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] * evt.w2list[W2LIST_GSEAGEN_WATER_INT_LEN] / evt.w2list[W2LIST_GSEAGEN_P_EARTH]); } /** * Get volume of given event according given weighter. * * \param type type * \param evt event * \return volume [m^3 Gev sr] */ inline double getVolume(const JType& type, const Evt& evt) { return (evt.w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] * evt.w2list[W2LIST_KM3BUU_WATER_INT_LEN] / evt.w2list[W2LIST_KM3BUU_P_EARTH]); } /** * Get volume of given event according given weighter. * * \param type type * \param evt event * \return volume [m^3 Gev sr] */ inline double getVolume(const JType& type, const Evt& evt) { using namespace JPHYSICS; const double l_int = 1.0e-6 * NUCLEON_MOLAR_MASS / (evt.w2list[W2LIST_GENHEN_SIG] * DENSITY_SEA_WATER * AVOGADRO); const double year = 60*60*24*365; // [s] return (evt.w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] * l_int / year / evt.w2list[W2LIST_GENHEN_P_EARTH]); } /** * Get volume of given event according given weighter. * * The return value should be normalised to * - number of generated events; * - solid angle covered by the generation; and * - bin width of the histogram. * * \param weighter weighter * \param evt event * \return volume [m^3 GeV sr] */ inline double getVolume(const JEvtWeight& weighter, const Evt& evt) { if (dynamic_cast(&weighter) != NULL) { return getVolume(JType(), evt); } if (dynamic_cast (&weighter) != NULL) { return getVolume(JType (), evt); } if (dynamic_cast (&weighter) != NULL) { return getVolume(JType (), evt); } THROW(JCastException, "No valid weighter."); } extern JEvtWeighter getEventWeighter; //!< Function object for mapping header to event weighter. } #endif