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 "km3net-dataformat/offline/io_ascii.hh"
 
#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JCircle2D.hh"
 
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
 
#include "JPhysics/JConstants.hh"
#include "JPhysics/JDeltaRays.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.\n
   *
   * Note: The optional parameter `can` is used only when the given track\n
   *       corresponds to a muon
   *
   * \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)) {
 
      const bool isMuon = is_muon(track);
      
      if (isMuon) {
 
        // 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 * JDeltaRays::getEnergyLossFromMuon(Emidpoint);
      
        Evis = dEb + dEc + dEd;
          
      } else if (JPDB::getInstance().hasPDG(track.type) &&
                 can.is_inside(getPosition(track))) {
        
        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);
  }
  
  
  /**
   * Get visible energy of the leading lepton of a neutrino interaction.
   *
   * \param  event         event
   * \param  can           containment volume
   * \return               visible energy of leading lepton [GeV]
   */
  inline double getVisibleEnergyLeadingLepton(const Evt&         event,
                                              const JCylinder3D& can = getMaximumContainmentVolume())
  {
    using namespace std;
    using namespace JPP;
    
    double Evis = 0.0;
 
    const Trk& leading_lepton = get_leading_lepton(event);
 
    if (is_finalstate(leading_lepton)) {
      Evis = getVisibleEnergy(leading_lepton, can);
    } else { // For tau-leptons
      for (vector<Trk>::const_iterator track = event.mc_trks.cbegin(); track != event.mc_trks.cend(); ++track) {
        if (is_finalstate(*track) && track->mother_id == leading_lepton.id) { 
          Evis += getVisibleEnergy(*track, can);
        }
      }
    }
    
    return Evis;
  }
 
 
  /**
   * Get visible energy vector of the leading lepton of a neutrino interaction.
   *
   * \param  event         event
   * \param  can           containment volume
   * \return               visible energy vector of leading lepton [GeV]
   */
  inline Vec getVisibleEnergyVectorLeadingLepton(const Evt&         event,
                                                 const JCylinder3D& can = getMaximumContainmentVolume())
  {
    using namespace std;
    using namespace JPP;
 
    Vec Evis(0.0, 0.0, 0.0);        
    
    const Trk& leading_lepton = get_leading_lepton(event);
 
    if (is_finalstate(leading_lepton)) {
      Evis = getVisibleEnergyVector(leading_lepton, can);
    } else { // For tau-leptons
      for (vector<Trk>::const_iterator track = event.mc_trks.cbegin(); track != event.mc_trks.cend(); ++track) {
        if (is_finalstate(*track) && track->mother_id == leading_lepton.id) { 
          Evis += getVisibleEnergyVector(*track, can);
        }
      }
    }
    
    return Evis;
  }
}
 
#endif