JEvtWeightToolkit.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTTOOLKIT__
#define __JAANET__JEVTWEIGHTTOOLKIT__

#include <map>
#include <vector>

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

#include "JLang/JException.hh"
#include "JLang/JSharedPointer.hh"

#include "Jeep/JProperties.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 {
  
  /**
   * Look-up table for event weighters.
   */
  struct JEvtWeighter :
    public std::vector< JLANG::JSharedPointer<JEvtWeight> >
  {
    /**
     * Constructor
     */
    JEvtWeighter()
    {
      this->push_back(new JEvtWeightGSeaGen());
      this->push_back(new JEvtWeightKM3BUU());
      this->push_back(new JEvtWeightCorsika());
      this->push_back(new JEvtWeightMupage());
      this->push_back(new JEvtWeightDAQ());
      this->push_back(new JEvtWeightMiscellaneous());
    }

      
    /**
     * Get event weighter corresponding to given header.
     *
     * \param  header          header
     * \return                 event weighter
     */
    const JEvtWeight& operator()(const JHead& header) const
    {
      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 :
    std::vector<int>
  {
    /**
     * Default constructor.
     */
    JNeutrinoTypeCollection()
    {}


    /**
     * 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 type; in >> type; ) {
        
        if (abs(type) == TRACK_TYPE_NUE  ||
            abs(type) == TRACK_TYPE_NUMU ||
            abs(type) == TRACK_TYPE_NUTAU) {
          
          collection.push_back(type);
          
        } else {

          THROW(JValueOutOfRange, "JNeutrinoTypeCollection::operator>>(): Invalid particle type: " << type);
        }
      }

      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)
    {
      for (std::vector<int>::const_iterator i = collection.cbegin(); i != collection.cend(); ++i) {
        out << ' ' << *i;
      }
      
      return out;
    }
  };


  /**
   * Auxiliary class for parsing multiparticle fluxes.
   */
  struct JFluxMapParser :
    public JEEP::JProperties,
    public JOscProbHelper
  {
    /**
     * Constructor.
     */
    JFluxMapParser()
    {
      this->put("flat",        flatFluxes);
      this->put("powerlaw",    powerLawFluxes);
      this->put("atmospheric", atmosphericFluxes);      
    }


    /**
     * 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));
      }

      if (is_valid()) {

        const JOscProb& oscProbInterface = *(this->get());
        
        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  fluxMapParser   flux map parser
     * \return                 input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JFluxMapParser& fluxMapParser)
    {
      using namespace std;
      using namespace JPP;
      
      JStringStream is(in);

      if (getFileStatus(is.str().c_str())) {
        is.load();
      }

      is >> static_cast<JProperties&>(fluxMapParser);

      return in;
    }

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

    JNeutrinoTypeCollection      atmosphericFluxes;   //!< Atmospheric neutrino flux functions
  };

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

#endif