#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 (is_valid()) {
      
	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