JEvtWeightToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTTOOLKIT__
#define __JAANET__JEVTWEIGHTTOOLKIT__

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

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"

#include "JLang/JException.hh"

#include "Jeep/JProperties.hh"

#include "JOscProb/JOscProbInterface.hh"

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"

#include "JAAnet/JEvtWeight.hh"
#include "JAAnet/JEvtWeightDAQ.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
 */

namespace JAANET {}
namespace JPP { using namespace JAANET; }

namespace JAANET {

  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 JEvtWeightDAQ());
      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;
    }
  };

    
  extern JEvtWeighter getEventWeighter;             //!< Function object for mapping header to event weighter.
}

#endif