JEvtWeightToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTTOOLKIT__
#define __JAANET__JEVTWEIGHTTOOLKIT__

#include <map>
#include <vector>
#include <memory>

#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/JException.hh"
#include "JLang/JType.hh"

#include "Jeep/JProperties.hh"

#include "JOscProb/JOscProbInterface.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/JFlatFlux.hh"
#include "JAAnet/JPowerLawFlux.hh"
#include "JAAnet/JAtmosphericNeutrinoFlux.hh"
#include "JAAnet/JEvtWeightFactorMultiParticle.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 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.");
    }
  };


  /**
   * Auxiliary class for parsing a vector of neutrino PDG identifiers.
   */
  struct JNeutrinoTypeCollection
  {
    /**
     * Default constructor.
     */
    JNeutrinoTypeCollection() :
      identifiers()
    {}


    /**
     * Add identifier.
     *
     * \param  identifier      PDG identifier
     */
    void put(const int identifier)
    {
      if (abs(identifier) == TRACK_TYPE_NUE  ||
          abs(identifier) == TRACK_TYPE_NUMU ||
          abs(identifier) == TRACK_TYPE_NUTAU) {

        identifiers.push_back(identifier);
        
      } else {
        
        THROW(JValueOutOfRange, "JNeutrinoTypeCollection::put(): Invalid PDG identifier: " << identifier);
      }
    }


    /**
     * Get size of this collection.
     *
     * \return                 size of neutrino type collection.
     */
    size_t size() const
    {
      return identifiers.size();
    }


    /**
     * Get constant iterator to the first element of the collection.
     *
     * \return                 constant iterator to the first element of the collection
     */
    std::vector<int>::const_iterator cbegin() const noexcept
    {
      return identifiers.cbegin();
    }


    /**
     * Get constant iterator to the last element of the collection.
     *
     * \return                 constant iterator to the last element of the collection
     */
    std::vector<int>::const_iterator cend() const noexcept
    {
      return identifiers.cend();
    }
    

    /**
     * Stream input.
     *
     * \param  in              input stream
     * \param  collection      collection of neutrino PDG types
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JNeutrinoTypeCollection& collection)
    {
      for (int identifier; in >> identifier; ) {
        collection.put(identifier);
      }

      return in;
    }


    /**
     * Stream output.
     *
     * \param  out             output stream
     * \param  collection      collection of neutrino PDG types
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JNeutrinoTypeCollection& collection)
    {
      using namespace std;

      const vector<int>& identifiers = collection.identifiers;
      
      for (vector<int>::const_iterator i = identifiers.cbegin(); i != identifiers.cend(); ++i) {
        out << ' ' << *i;
      }
      
      return out;
    }

    
  private:
    
    std::vector<int> identifiers; //!< Container for identifiers
  };


  /**
   * Auxiliary class for parsing multiparticle fluxes.
   */
  struct JFluxMap :
    public JOscProbHelper
  {
    /**
     * Default constructor.
     */
    JFluxMap()
    {}


    /**
     * Get multiparticle flux function.
     *
     * \return                 multiparticle flux function
     */
    JFluxMultiParticle getMultiParticleFlux() const
    {
      using namespace std;
      using namespace JPP;

      JFluxMultiParticle multiFlux;

      for (map<int, JFlatFlux>::const_iterator i = flatFluxes.cbegin(); i != flatFluxes.cend(); ++i) {
        multiFlux.insert(i->first, make_fluxFunction(i->second));
      }

      for (map<int, JPowerLawFlux>::const_iterator i = powerLawFluxes.cbegin(); i != powerLawFluxes.cend(); ++i) {
        multiFlux.insert(i->first, make_fluxFunction(i->second));
      }

      const JOscProbInterface& oscProbInterface = getOscProbInterface();
        
      const JAtmosphericNeutrinoFlux atmFlux(oscProbInterface);
        
      const JEvtWeightFactorFunction<JAtmosphericNeutrinoFlux, JFlux> atmFluxFunction = make_fluxFunction(atmFlux);
        
      for (vector<int>::const_iterator i = atmosphericFluxes.cbegin(); i != atmosphericFluxes.cend(); ++i) {
        multiFlux.insert(*i, atmFluxFunction);
      }

      return multiFlux;
    }


    /**
     * Conversion operator.
     *
     * \return                 multiparticle flux function
     */
    operator JFluxMultiParticle() const
    {
      return getMultiParticleFlux();
    }


    /**
     * Stream input.
     *
     * \param  in              input stream
     * \param  fluxMap         flux map
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JFluxMap& fluxMap)
    {
      return getProperties(fluxMap).read(in);
    }

    
    /**
     * Stream output.
     *
     * \param  out             output stream
     * \param  fluxMap         flux map
     * \return                 output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JFluxMap& fluxMap)
    {
      return getProperties(fluxMap).write(out);
    }
    
    
    std::map<int, JFlatFlux>     flatFluxes;        //!< Uniform flux functions
    std::map<int, JPowerLawFlux> powerLawFluxes;    //!< Power-law flux functions  

    JNeutrinoTypeCollection      atmosphericFluxes; //!< Atmospheric neutrino flux functions

  private:
    /**
     * Get properties of this class.
     *
     * \param  object          flux map object
     * \return                 properties
     */
    template<class JFluxMap_t>
    static inline JProperties getProperties(JFluxMap_t& object)
    {
      JProperties properties;

      properties.insert(gmake_property(object.flatFluxes));
      properties.insert(gmake_property(object.powerLawFluxes));
      properties.insert(gmake_property(object.atmosphericFluxes));

      return properties;
    }
  };


  /**
   * 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