JWeightGSeaGen.hh

Go to the documentation of this file.

#ifndef __JAANET__JWEIGHTGSEAGEN__
#define __JAANET__JWEIGHTGSEAGEN__

#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/JWeightEvent.hh"

/**
 * \author bjung
 */

namespace JAANET {

  using JLANG::JClonable;
  using JLANG::JIndexOutOfRange;
  using JLANG::JValueOutOfRange;

  
  /**
   * Implementation of event weighing for GSeaGen data.
   */
  struct JWeightGSeaGen :
    public JClonable<JWeightEvent, JWeightGSeaGen>
  {
    /**
     * Default Constructor.
     */
    JWeightGSeaGen() : W(0.0)
    {}

    
    /** 
     * Constructor.
     *
     * \param header            header
     */
    JWeightGSeaGen(const JHead& header)
    {
      configure(header);
    }
    
    
    /**
     * Configuration.
     *
     * \param  header            header
     */
    virtual void configure(const JHead& header) override
    {
      if (check(header) && (header.tgen.numberOfSeconds  > 0.0 &&
                            header.genvol.numberOfEvents > 0.0)) {
        
        W = 1.0 / header.genvol.numberOfEvents / header.tgen.numberOfSeconds;
        
      } else {

        THROW(JValueOutOfRange, "JWeightGSeaGen::configure(): Provided header is inconsistent with GSeaGen.");
      }
    }


    /**
     * Check whether header is consistent with this event weighter.
     *
     * \param  header            header
     * \return                   true if consistent; else false
     */
    virtual bool check(const JHead& header) const override
    {
      return is_gseagen(header);
    }

    
    /**
     * Get weight of given event.
     *
     * \param  evt               event
     * \return                   weight        [Hz]
     */
    virtual double getWeight(const Evt& evt) const override
    {
      if (evt.w.size() > 2 && !(evt.w[2] < 0.0)) {
        
        return W * evt.w[2];
        
      } else {

        if (evt.w.size() < 3) {
          THROW(JIndexOutOfRange, "JWeightGSeaGen::getWeight(): w3-weight is empty.");
        } else {
          THROW(JValueOutOfRange, "JWeightGSeaGen::getWeight(): w3-weight is negative.");
        }
      }
    }
    
    
    /**
     * Get weight of given event.
     *
     * \param  evt               event
     * \param  flux              neutrino flux [m^-2 s^-1 sr^-1 GeV^-1]
     * \return                   weight        [Hz]
     */
    virtual double getWeight(const Evt&   evt,
                             const double flux) const override
    {
      if (evt.w.size() > 1 && !(evt.w[1] < 0.0 || flux < 0.0)) {
        
        return W * evt.w[1] * flux;
        
      } else {

        if (evt.w.size() < 2) {
          THROW(JIndexOutOfRange, "JWeightGSeaGen::getWeight(): w2-weight is empty.");
        } else {
          THROW(JValueOutOfRange, "JWeightGSeaGen::getWeight(): w2-weight or flux is negative");
        }
      }
    }

  private:
    double W;
  };
}

#endif