#ifndef __JSIRENE__JVISIBLEENERGYTOOLKIT__
#define __JSIRENE__JVISIBLEENERGYTOOLKIT__

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

#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JParticleTypes.hh"
#include "JAAnet/JPDB.hh"

#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JCircle2D.hh"

#include "JGeometry3D/JCylinder3D.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JGeane.hh"

#include "JSirene/pythia.hh"


/**
 * \file
 *
 * Auxiliary methods for evaluating visible energies.
 * \author bjung
 */
namespace JSIRENE {}
namespace JPP { using namespace JSIRENE; }

namespace JSIRENE {

  using JAANET::JPDB;
  using JAANET::is_muon;
  using JAANET::getAxis;
  using JAANET::getPosition;  
  using JAANET::getKineticEnergy;

  using JGEOMETRY2D::JVector2D;
  using JGEOMETRY2D::JCircle2D;
  
  using JGEOMETRY3D::JCylinder3D;

  using JPHYSICS::geanc;
  using JPHYSICS::gWater;
  using JPHYSICS::MASS_MUON;
  using JPHYSICS::JGeaneWater;
  using JPHYSICS::getSinThetaC;
  
  
  /**
   * Get the visible energy of a track.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  track         track
   * \param  can           detector can  
   * \return               visible energy [GeV]
   */
  inline Vec getVisibleEnergy(const Trk&         track,
			      const JCylinder3D& can) {
    using namespace std;
    using namespace JPP;

    double Evis = 0.0;
    
    if (track.is_finalstate()) {
      if (is_muon(track)) {

	// Determine muon pathlength inside detector [m]

	const JCylinder3D::intersection_type& intersection = can.getIntersection(getAxis(track));

	const double Lmuon = gWater.getX(track.E, MASS_MUON / getSinThetaC()); 
	const double Leff  = (intersection.first < 0.0        ? 
			      min(Lmuon, intersection.second) : 
			      min(Lmuon, intersection.second) - intersection.first);

	// Determine visible energy deposition [GeV]
      
	const double dEb = gWater.getEb(track.E, Leff);
	const double dEc = Leff / geanc();
      
	Evis = dEb + dEc;
	  
      } else if (!is_neutrino(track) && JPDB::getInstance().hasPDG(track.type)) {

	Evis = pythia(track.type, getKineticEnergy(track));
      }
    }

    return Evis * track.dir / track.dir.len();
  }


  /**
   * Get the visible energy of a track, assuming an infinite detector volume.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  track         track
   */  
  inline Vec getVisibleEnergy(const Trk& track)
  {
    static const double Rearth = 6386 * 1e3; // Radius of the Earth [m]
    
    static const JCylinder3D can(JCircle2D(JVector2D(), Rearth), Rearth, Rearth);

    return getVisibleEnergy(track, can);
  }


  /**
   * Get the visible energy of an event.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  evt           event
   * \param  can           detector can  
   * \return               visible energy [GeV]
   */
  inline Vec getVisibleEnergy(const Evt&         evt,
			      const JCylinder3D& can)
  {
    Vec Evis(0.0, 0.0, 0.0);    

    for (std::vector<Trk>::const_iterator track = evt.mc_trks.begin(); track != evt.mc_trks.end(); ++track) {
      Evis += getVisibleEnergy(*track, can);
    }

    return Evis;
  }
    

  /**
   * Get the visible energy of an event, assuming an infinite detector volume.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  evt           event
   * \return               visible energy [GeV]
   */
  inline Vec getVisibleEnergy(const Evt& evt)
  {
    Vec Evis(0.0, 0.0, 0.0);

    for (std::vector<Trk>::const_iterator track = evt.mc_trks.cbegin(); track != evt.mc_trks.cend(); ++track) {
      Evis += getVisibleEnergy(*track);
    }
    
    return Evis;
  }
}

#endif