#ifndef __JSIRENE__JVISIBLEENERGYTOOLKIT__
#define __JSIRENE__JVISIBLEENERGYTOOLKIT__

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

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

#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"

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

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

#include "JSirene/pythia.hh"


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

namespace JSIRENE {
  
  using JGEOMETRY3D::JCylinder3D;

  
  /**
   * Auxiliary function to retrieve the maximum cylindrical containment volume.
   *
   * \return               maximum cylindrical containment volume
   */
  inline const JCylinder3D getMaximumContainmentVolume()
  {
    using namespace JPP;

    const double R_Earth = R_EARTH_KM * 1e3; // m
    
    const JVector2D center(0.0, 0.0);
    const JCircle2D circle(center, R_Earth);
    
    return JCylinder3D(circle, -R_Earth, R_Earth);
  }


  /**
   * 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]
   */
  double getVisibleEnergy(const Trk&         track,
			  const JCylinder3D& can = getMaximumContainmentVolume())
  {
    using namespace std;
    using namespace JPP;

    double Evis = 0.0;
    
    if (is_finalstate(track)) {
      
      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  = (min(Lmuon, max(intersection.second, 0.0)) -
			      min(Lmuon, max(intersection.first,  0.0)));
				  

	// Determine visible energy deposition [GeV]

	const double Emidpoint = gWater.getE(track.E, Lmuon/2.0);
      
	const double dEb = gWater.getEb(track.E, Leff);
	const double dEc = Leff / geanc();
	const double dEd = Leff * getDeltaRaysFromMuon(Emidpoint);
      
	Evis = dEb + dEc + dEd;
	  
      } else if (!is_neutrino(track) && JPDB::getInstance().hasPDG(track.type)) {

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

    return Evis;
  }
  
  
  /**
   * Get the visible energy vector of a track.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  track         track
   * \param  can           detector can  
   * \return               visible energy vector [GeV]
   */
  inline Vec getVisibleEnergyVector(const Trk&         track,
				    const JCylinder3D& can = getMaximumContainmentVolume()) {
    return getVisibleEnergy(track, can) * track.dir;
  }


  /**
   * Get the visible energy of a given range of tracks.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  __begin       start of track data
   * \param  __end         end of track data
   * \param  can           detector can  
   * \return               visible energy [GeV]
   */
  inline double getVisibleEnergy(std::vector<Trk>::const_iterator __begin,
				 std::vector<Trk>::const_iterator __end,
				 const JCylinder3D&               can = getMaximumContainmentVolume())
  {
    using namespace std;
    
    double Evis = 0.0;
    
    for (vector<Trk>::const_iterator track = __begin; track != __end; ++track) {
      Evis += getVisibleEnergy(*track, can);
    }
    
    return Evis;
  }


  /**
   * Get the visible energy vector of a given range of tracks.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  __begin       start of track data
   * \param  __end         end of track data
   * \param  can           detector can  
   * \return               visible energy vector [GeV]
   */
  inline Vec getVisibleEnergyVector(std::vector<Trk>::const_iterator __begin,
				    std::vector<Trk>::const_iterator __end,
				    const JCylinder3D&               can = getMaximumContainmentVolume())
  {
    using namespace std;
    
    Vec Evis(0.0, 0.0, 0.0);    

    for (vector<Trk>::const_iterator track = __begin; track != __end; ++track) {
      Evis += getVisibleEnergyVector(*track, can);
    }

    return Evis;
  }


  /**
   * Get the visible energy vector of an event.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  evt           event
   * \param  can           detector can  
   * \return               visible energy [GeV]
   */
  inline double getVisibleEnergy(const Evt&         evt,
				 const JCylinder3D& can  = getMaximumContainmentVolume())
  {
    return getVisibleEnergy(evt.mc_trks.begin(), evt.mc_trks.end(), can);
  }
  

  /**
   * Get the visible energy vector of an event.\n
   * This method accounts for muon radiative energy losses.
   *
   * \param  evt           event
   * \param  can           detector can  
   * \return               visible energy vector [GeV]
   */
  inline Vec getVisibleEnergyVector(const Evt&         evt,
				    const JCylinder3D& can  = getMaximumContainmentVolume())
  {
    return getVisibleEnergyVector(evt.mc_trks.begin(), evt.mc_trks.end(), can);
  }
}

#endif