JMultiParticleFlux.hh

Go to the documentation of this file.

#ifndef __JAANET__JMULTIPARTICLEFLUX__
#define __JAANET__JMULTIPARTICLEFLUX__

#include <map>

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

#include "JLang/JClonable.hh"
#include "JLang/JPredicate.hh"
#include "JLang/JException.hh"

#include "JAAnet/JFlux.hh"
#include "JAAnet/JFluxHelper.hh"
#include "JAAnet/JAAnetToolkit.hh"

/**
 * \author bjung
 */

namespace JAANET {
  
  using JLANG::JClonable;
  using JLANG::JValueOutOfRange;

  /**
   * Implementation of flux function for multiple particle types.
   */
  struct JMultiParticleFlux :
    public JClonable<JFlux, JMultiParticleFlux>,
    public std::map<int, JFluxHelper>
  {
    typedef std::map<int, JFluxHelper>      map_type;
    
    /**
     * Default constructor.
     */
    JMultiParticleFlux() :
      map_type()
    {}

    
    /**
     * Constructor.
     *
     * \param  type                   particle PDG code
     * \param  flux                   flux function
     * \param  args                   remaining pairs of PDG codes and flux functions
     */
    template<class ...Args>
    JMultiParticleFlux(const int      type,
                       const JFlux&   flux,
                       const Args& ...args)
    {
      insert(type, flux, args...);
    }


    /**
     * Insert pair of particle code and flux function.
     *
     * \param  type                    particle PDG code
     * \param  flux                    flux function
     * \param  args                    remaining pairs of PDG codes and flux functions
     */ 
    void insert(const int    type,
                const JFlux& flux)
    {
      map_type::insert(std::make_pair(type, JFluxHelper(flux)));
    }
    

    /**
     * Insert pairs of particle codes and flux functions.
     *
     * \param  type                    particle PDG code
     * \param  flux                    flux function
     * \param  args                    remaining pairs of PDG codes and flux functions
     */ 
    template<class ...Args>
    void insert(const int      type,
                const JFlux&   flux,
                const Args& ...args)
    {
      insert(type, flux);
      insert(args...);
    }


    /**
     * Get flux of given event.
     *
     * \param  evt               event
     * \return                   flux         [GeV * m^-2 * sr^-1 * s^-1]
     */
    double getFlux(const Evt& evt) const override
    {
      using namespace std;
      using namespace JPP;

      vector<Trk>::const_iterator primary = find_if(evt.mc_trks.cbegin(),
                                                    evt.mc_trks.cend(),
                                                    make_predicate(&Trk::is_primary, true));

      if (primary != evt.mc_trks.cend()) {
        
        return this->at(primary->type).getFlux(evt);
        
      } else if (has_neutrino(evt)) {
        
        const Trk& neutrino = get_neutrino(evt);

        return this->at(neutrino.type).getFlux(evt);
        
      } else {
        
        THROW(JValueOutOfRange, "JMultiParticleFlux::getFlux(): No primary track found for event " << evt.id << ".");
      }
    }
  };
}

#endif