JEvtWeightFactorGSeaGen.hh

Go to the documentation of this file.

#ifndef __JAANET__JEVTWEIGHTFACTORGSEAGEN__
#define __JAANET__JEVTWEIGHTFACTORGSEAGEN__

#include "km3net-dataformat/definitions/w2list_gseagen.hh"

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

#include "JLang/JManip.hh"
#include "JLang/JToken.hh"
#include "JLang/JException.hh"

#include "JAAnet/JFlux.hh"
#include "JAAnet/JEvtWeightFactorHelper.hh"

#include "JGizmo/JGizmoToolkit.hh"

#include "TFormula.h"


/**
 * \author bjung
 */

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

namespace JAANET {

  /**
   * Implementation of reweighting factor for simulated neutrino interactions according to a specifiable ROOT TFormula.
   *
   * Note: The ROOT TFormula may assume any number of parameters, but should be restricted to\n
   *       the physical variables listed among `JEvtWeightFactorGSeaGen::variables`.\n
   *       These variables may be specified within the ROOT TFormula as 'x[<index>]',\n
   *       where <index> corresponds to the index of the desired physical variable within `JEvtWeightFactorGSeaGen::variables`.
   */
  struct JEvtWeightFactorGSeaGen :
    public JFluxHelper,
    public TFormula
  {
    /**
     * Indices of reweighting variables for GSeaGen.
     */
    enum variables {
      COSTH,              //!< Cosine zenith angle
      ENERGY,             //!< Initial state energy
      BJORKEN_X,          //!< Björken-x (= fractional momentum carried by the struck nucleon)
      INELASTICITY,       //!< Inelasticity (= Björken-y)
      NUMBER_OF_VARIABLES //!< Number of reweighting variables; \n
                          //!< N.B.\ This enum value needs to be specified last!
    }; 

    
    /**
     * Default constructor.
     */
    JEvtWeightFactorGSeaGen() :
      JFluxHelper(),
      TFormula()
    {}


    /**
     * Constructor.
     *
     * \param  flux            flux function
     * \param  name            name
     * \param  formula         formula
     */
    JEvtWeightFactorGSeaGen(const JFlux& flux,
                            const char*  name,
                            const char*  formula) :
      JFluxHelper(flux),
      TFormula(name, formula)
    {}

    
    /**
     * Get weighting factor for given event.
     *
     * \param  evt             event
     * \return                 weighting factor
     */
    double operator()(const Evt& evt) const
    {
      using namespace std;
      using namespace JPP;

      if (static_cast<const JFluxHelper&>(*this)) {
      
        Double_t vars[NUMBER_OF_VARIABLES] = { 0.0 };
        
        const Trk&   neutrino = get_neutrino(evt);
        const double flux     = getFactor(evt);
      
        vars[COSTH]        = neutrino.dir.z;
        vars[ENERGY]       = getE0(evt);
        vars[BJORKEN_X]    = evt.w2list[W2LIST_GSEAGEN_BX];
        vars[INELASTICITY] = evt.w2list[W2LIST_GSEAGEN_BY];
      
        return this->DoEval(&vars[0]) * flux;
        
      } else {

        THROW(JNullPointerException, "JEvtWeightFactorGSeaGen::operator(): Unspecified flux function.");
      }
    }


    /**
     * Stream input.
     *
     * \param  in                     input stream
     * \param  object                 gSeaGen event-weight factor
     * \return                        input stream
     */
    inline friend std::istream& operator>>(std::istream&            in,
                                           JEvtWeightFactorGSeaGen& object)
    {
      using namespace std;
      using namespace JPP;

      typedef JToken<';'>   JToken_t;

      JToken_t formula;
      JToken_t equation;
      
      in >> formula;

      object.Clear();
      object.Compile(formula.c_str());
      
      for (int count = 0; count < object.GetNpar() && in >> equation; ++count) {

        const int    index = getParameter(equation);
        const double value = getValue(equation, 0);

        object.SetParameter(index, value);
      }
      
      return in;
    }

      
      /**
       * Stream output.
       *
       * \param  out                    output stream
       * \param  object                 gSeaGen event-weight factor
       * \return                        output stream
       */
      inline friend std::ostream& operator<<(std::ostream&                  out,
                                             const JEvtWeightFactorGSeaGen& object)
      {
        using namespace std;
        using namespace JPP;

        out << object.GetExpFormula() << ';';
        
        for (int i = 0; i < object.GetNpar(); ++i) {
          out << " p" << i << " = " << SCIENTIFIC(6,3) << object.GetParameter(i) << ';';
        }
        
        return out;
      }
  };
}

#endif