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)
      INTERACTION_TYPE,    //!< Interaction channel type
      CURRENT_TYPE,        //!< Weak current type (CC or NC)
      XSEC,                //!< Exclusive total cross-section of the interaction
      XSEC_MEAN,           //!< Average interaction cross-section per nucleon along neutrino path [m2]
      XSEC_DIFFERENTIAL,   //!< Differential cross-section of the interaction
      XSEC_WATER,          //!< Inclusive cross-section in water
      INT_LENGTH_WATER,    //!< Interaction length in pure water
      COLUMN_DEPTH,        //!< Column density [m.w.e]
      P_EARTH,             //!< Earth transmission probability
      P_SCALE,             //!< GENIE ineraction probability scale
      TARGET_A,            //!< Number of nucleons in the target
      TARGET_Z,            //!< Number of protons  in the target
      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];
        vars[INTERACTION_TYPE]  = evt.w2list[W2LIST_GSEAGEN_ICHAN];
        vars[CURRENT_TYPE]      = evt.w2list[W2LIST_GSEAGEN_CC];
        vars[XSEC]              = evt.w2list[W2LIST_GSEAGEN_XSEC];
        vars[XSEC_MEAN]         = evt.w2list[W2LIST_GSEAGEN_XSEC_MEAN];
        vars[XSEC_DIFFERENTIAL] = evt.w2list[W2LIST_GSEAGEN_DXSEC];
        vars[XSEC_WATER]        = evt.w2list[W2LIST_GSEAGEN_WATERXSEC];
        vars[INT_LENGTH_WATER]  = evt.w2list[W2LIST_GSEAGEN_WATER_INT_LEN];
        vars[COLUMN_DEPTH]      = evt.w2list[W2LIST_GSEAGEN_COLUMN_DEPTH];
        vars[P_EARTH]           = evt.w2list[W2LIST_GSEAGEN_P_EARTH];
        vars[P_SCALE]           = evt.w2list[W2LIST_GSEAGEN_P_SCALE];
        vars[TARGET_A]          = evt.w2list[W2LIST_GSEAGEN_TARGETA];
        vars[TARGET_Z]          = evt.w2list[W2LIST_GSEAGEN_TARGETZ];

        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