JVisibleEnergyToolkit.hh

Go to the documentation of this file.

#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 "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::R_EARTH_KM;
  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 (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  = (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 JCylinder3D can(JCircle2D(JVector2D(), R_EARTH_KM * 1e3), -R_EARTH_KM * 1e3, R_EARTH_KM * 1e3);

    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