#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::const_iterator __begin, std::vector::const_iterator __end, const JCylinder3D& can = getMaximumContainmentVolume()) { using namespace std; double Evis = 0.0; for (vector::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::const_iterator __begin, std::vector::const_iterator __end, const JCylinder3D& can = getMaximumContainmentVolume()) { using namespace std; Vec Evis(0.0, 0.0, 0.0); for (vector::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