JEvtWeightToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTTOOLKIT__
#define __JAANET__JEVTWEIGHTTOOLKIT__
 
#include <vector>
 
#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<std::shared_ptr<JEvtWeight> >
  {
    /**
     * 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<JEvtWeightGSeaGen>& 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<JEvtWeightKM3BUU>& 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<JEvtWeightGenhen>& 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<const JEvtWeightGSeaGen*>(&weighter) != NULL) { return getVolume(JType<JEvtWeightGSeaGen>(), evt); }
    if (dynamic_cast<const JEvtWeightGenhen*> (&weighter) != NULL) { return getVolume(JType<JEvtWeightGenhen> (), evt); }
    if (dynamic_cast<const JEvtWeightKM3BUU*> (&weighter) != NULL) { return getVolume(JType<JEvtWeightKM3BUU> (), evt); }
 
    THROW(JCastException, "No valid weighter.");
  }
    
  extern JEvtWeighter getEventWeighter;             //!< Function object for mapping header to event weighter.
}
 
#endif