JEvtWeightCorsika.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTCORSIKA__
#define __JAANET__JEVTWEIGHTCORSIKA__
 
#include "km3net-dataformat/definitions/weightlist.hh"
#include "km3net-dataformat/definitions/w2list_gseagen.hh"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
#include "JLang/JException.hh"
#include "JLang/JClonable.hh"
 
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JEvtWeight.hh"
#include "JAAnet/JEvtWeightFactorHelper.hh"
#include "JAAnet/JFlux.hh"
#include "JAAnet/JFluxHelper.hh"
 
 
/**
 * \author bjung
 */
 
namespace JAANET {}
namespace JPP { using namespace JAANET; }
 
namespace JAANET {
  
  using JLANG::JClonable;
 
  
  /**
   * Implementation of event weighting for Corsika data.
   */
  struct JEvtWeightCorsika final :
    public JFluxHelper,
    public JClonable<JEvtWeight, JEvtWeightCorsika>
  {
    /**
     * Default Constructor.
     */
    JEvtWeightCorsika() :
      JFluxHelper()
    {}
 
 
    /** 
     * Constructor.
     *
     * \param  header            header
     */    
    JEvtWeightCorsika(const JHead& header) :
      JFluxHelper()
    {
      configure(header);
    }
 
 
    /** 
     * Constructor.
     *
     * \param  header            header
     * \param  flux              flux
     */
    JEvtWeightCorsika(const JHead& header,
                      const JFlux& flux)
    {
      configure(header);
      JFluxHelper::configure(flux);
    }
 
 
    /**
     * Get name.
     *
     * \return                   name
     */
    const char* const getName() const override final
    {
      return "Corsika";
    }
    
    
    /**
     * Weight configuration.
     *
     * \param  header            header
     */
    void configure(const JHead& header) override final
    {
      using namespace JPP;
      
      if (check(header)) {
 
        setNormalisation(1.0 / header.genvol.numberOfEvents);
        
      } else {
 
        THROW(JValueOutOfRange, "JEvtWeightCorsika::configure(): Provided header is inconsistent with Corsika.");
      }
    }
 
 
    /**
     * Check whether header is consistent with this event weighter.
     *
     * \param  header            header
     * \return                   true if consistent; else false
     */
    bool check(const JHead& header) const override final
    {
      return is_corsika(header);
    }
 
    
    /**
     * Get weight of given event.
     *
     * \param  evt               event
     * \return                   weight        [Hz]
     */
    double getWeight(const Evt& evt) const override final
    {
      using namespace JPP;
      
      const double Nretries = evt.w2list[W2LIST_GSEAGEN_N_RETRIES];
      
      if (static_cast<const JFluxHelper&>(*this) &&
          evt.w.size() > WEIGHTLIST_DIFFERENTIAL_EVENT_RATE) {
          
        return evt.w[WEIGHTLIST_DIFFERENTIAL_EVENT_RATE] * getFactor(evt) * getNormalisation(evt) / (1 + Nretries);
        
      } else if (evt.w.size() > WEIGHTLIST_RESCALED_EVENT_RATE) {
        
        return evt.w[WEIGHTLIST_RESCALED_EVENT_RATE];
 
      } else if (evt.w.size() > WEIGHTLIST_EVENT_RATE) {
 
        return evt.w[WEIGHTLIST_EVENT_RATE] * getNormalisation(evt) / (1 + Nretries);
        
      } else {
 
        if (!static_cast<const JFluxHelper&>(*this)) {
          THROW(JNullPointerException, "JEvtWeightCorsika::getWeight(): Unspecified flux function.");
        } else {
          THROW(JIndexOutOfRange,      "JEvtWeightCorsika::getWeight(): Empty " << (evt.w.size() < 3 ? "w2-" : "w3-") << "weight.");
        }
      }
    }
  };
}
 
#endif