earthmodel-service C++ API Reference

struct densityPathIntegrand

A helper structure for path integration of density.

Public Functions

inline densityPathIntegrand(const EarthParam &e, const I3Position &p, const I3Direction &d, IntegType integ_type)
inline double operator()(double dist) const

Evaluates the density at the appropriate point along the parameterized track.

Public Members

const EarthParam &medium

The medium in which the integration should be done.

I3Position pos

The base point for the path along which the integration is performed. Must be in earth-centered coordinates.

I3Direction dir

The direction along which the integration is performed. Must be in earth-centered coordinates.

IntegType intg_type
class EarthModelService : public I3ServiceBase
#include <EarthModelService.h>

Public Types

enum MediumType

Values:

enumerator INNERCORE
enumerator OUTERCORE
enumerator MANTLE
enumerator ROCK
enumerator ICE
enumerator AIR
enumerator VACUUM
enumerator LOWERMANTLE
enumerator UPPERMANTLE
enumerator LLSVP
enum IceCapType

Values:

enumerator NOICE
enumerator ICESHEET
enumerator SIMPLEICECAP
enum IntegType

Values:

enumerator PATH
enumerator RADIUS
enumerator CIRCLE
enumerator SPHERE
typedef std::map<MediumType, std::map<int, double>> MatRatioMap
typedef std::map<double, EarthParam> EarthParamMap

Public Functions

EarthModelService(const I3Context &c)

Constructor -.

EarthModelService(const std::string &name = "EarthModelService", const std::string &tablepath = "", const std::vector<std::string> &earthmodels = std::vector<std::string>(), const std::vector<std::string> &materialmodels = std::vector<std::string>(), const std::string &icecapname = "SimpleIceCap", double icecapangle = 20.0 * I3Units::degree, double detectordepth = 1948.0 * I3Units::m)

constructor for pybinding

Parameters:

  • name – [in] : (dummy) name to pass to I3ServiceBase

  • tablepath – [in] : path to table dir. blank: use default path

  • earthmodels – [in] : list of earthmodel files. You may add muliple files, and latter one overwrites former value. e.g. PREM_mmc + FLATCORE = PREM earth with flat density at core

  • materialmodels – [in] : list of material files. You may add muliple files, and latter one overwrites former value.

  • icecapname – [in] : type of icecap shape.

  • icecapangle – [in] valid only when SimpleIceCap is selected

  • detectordepth – [in] : depth of icecube center

virtual ~EarthModelService()

Destructor -

void Configure()

configure

const EarthParam &GetEarthParam(const I3Position &p_CE) const

Finds the material layer containing p_CE.

Parameters:

p_CE – the position for which to find the material layer, which must be specified in Earth-centered coordinates.

const double GetEarthDensityInCGS(const I3Position &p_CE) const

Computes the material density at the given point.

Parameters:

p_CE – the position at which the density is to be evaluated. This must be in Earth-centered coordinates.

Returns:

the density in g/cm^3

const double GetColumnDepthInCGS(const I3Position &from_posI3, const I3Position &to_posI3) const

This function calculates total column depth from from_posI3 to to_posI3.

GetColumnDepthInCGS

Parameters:

  • from_posI3 – [in] from position in I3 coordinate

  • to_posI3 – [in] to position in I3 coordinate

const double IntegrateDensityInCGS(const I3Position &from_posCE, const I3Position &to_posCE, IntegType intg_type = PATH) const

integrate density

IntegrateDensityInCGS

Parameters:

  • from_posCE – [in] from position in Earth centered coordinate

  • to_posCE – [in] to position in Earth centered coordinate

  • intg_type – [in] : integrate density with PATH : along a path from from_posCE to to_posCE. [g/cm2] RADIUS : similar to PATH but projected on radial direction [g/cm2] CIRCLE : 2*pi*r x RADIUS option [g/cm] SPHERE : 4*pi*r^2 x RADIUS option (volume mass)[g]

double DistanceForColumnDepthToPoint(const I3Position &to_posI3, const I3Direction &dirI3, double cDepth) const

Computes the distance along the given direction, ending at the given point, which must be traversed to accumulate the specified column depth.

Parameters:

  • to_posI3 – the endpoint of the path

  • dirI3 – the direction of the path

  • cDepth – the column depth required in g/cm^2

Returns:

the distance in meters

double DistanceForColumnDepthFromPoint(const I3Position &from_posI3, const I3Direction &dirI3, double cDepth) const

Computes the distance along the given direction, starting at the given point, which must be traversed to accumulate the specified column depth.

Parameters:

  • from_posI3 – the starting point of the path

  • dirI3 – the direction of the path

  • cDepth – the column depth required in g/cm^2

Returns:

the distance in meters

const double GetLeptonRangeInMeterFrom(const I3Position &posI3, const I3Direction &dirI3, double energy, bool isTau = false, EarthModelCalculator::LeptonRangeOption opt = EarthModelCalculator::DEFAULT, double scale = 1.0) const

This function calculates MuonRange [m.w.e] and convert it to distance [m] with given particle info, earth model and the start position (o)

GetLeptonRangeInMeterFrom 

      d[m]
|————-&#8212;| <————-&#8212;o dir start pos

Parameters:

  • posI3 – [in] start position to calculate range

  • dirI3 – [in] direction where the particle is moving to (particle direction)

  • energy – [in] particle energy

  • isTau – [in] if set the lepton will be treated as a tau, otherwsie as a muon.

  • opt, scale – [in] Used by EarthModelCalculator::GetLeptonRange() function. Leave it as defaut unless you understand well about the parameters.

Returns:

distance d[m]

const double GetLeptonRangeInMeterTo(const I3Position &posI3, const I3Direction &dirI3, double energy, bool isTau = false, EarthModelCalculator::LeptonRangeOption opt = EarthModelCalculator::DEFAULT, double scale = 1.0) const

This function calculates MuonRange [m.w.e] and convert it to distance [m] with given particle info, earth model and the end position (o)

GetLeptonRangeInMeterTo 

      d[m]
|————–&#8212;| o<————-&#8212;| end pos dir

Parameters:

  • posI3 – [in] end position of the conversion

  • dirI3 – [in] direction where the particle is moving to (particle direction)

  • energy – [in] particle energy

  • isTau – [in] if set the lepton will be treated as a tau, otherwsie as a muon.

  • opt, scale – [in] Used by EarthModelCalculator::GetLeptonRange() function. Leave it as defaut unless you understand well about the parameters.

Returns:

distance d[m]

double DistanceToNextBoundaryCrossing(const I3Position &posCE, const I3Direction &dirCE, bool &exitsAtmosphere = ignored_bool) const

Computes the disance to the next boundary between material layers along the given track.

Parameters:

  • posCE – the starting point of the track in EarthCenter coordinate

  • dirCE – the direction of the track in EarthCenter coordinate

  • exitsAtmosphere – [out] whether the next boundary crossing will be the outer surface of the atmosphere. This parameter is optional and may be omitted if this information is not of interest

const MediumType GetMedium(const I3Position &p_CE) const

Get medium type at a given point. YOU MUST USE Earth-centered coordinate for argument.

inline const MatRatioMap GetMatRatioMap() const

Get material ratio map (for genie)

const std::map<int, double> &GetMatRatioMap(MediumType med) const
double GetMatRatio(MediumType med, int id) const
inline const MatRatioMap GetPNERatioMap() const

Get ratio map of proton, neutron and electron (for nugen)

const std::map<int, double> &GetPNERatioMap(MediumType med) const
double GetPNERatio(MediumType med, int id) const
const double GetDistanceFromEarthEntranceToDetector(double zenith_rad) const

Get track path length from Earth entrance point of track to detector center

const std::vector<double> GetDistanceFromSphereSurfaceToDetector(double zenrad, double radius, double det_z) const
std::string PrintEarthParams() const

Printing Earth Params

inline boost::python::list GetEarthParamsList()

Get Earth Params.

const double GetPREM(double r) const

Earth model function. CAUTION : Return unit is [g/cm3] !

const I3Position GetEarthCoordPosFromDetCoordPos(const I3Position &p) const

convert I3Position in detector coordinate to Earth Center Coordinate

const I3Position GetDetCoordPosFromEarthCoordPos(const I3Position &p) const

convert I3Position in Earth Center Coordinate to Detector Coordinate

const I3Direction GetEarthCoordDirFromDetCoordDir(const I3Direction &p) const

convert I3Direction in Detector Coordinate to Earth Center Coordinate

const I3Direction GetDetCoordDirFromEarthCoordDir(const I3Direction &p) const

convert I3Direction in Earth Center Coordinate to Detector Coordinate

inline void SetPath(const std::string &s)

Set path of parameter files

void SetEarthModel(const std::vector<std::string> &s)

Set name of Earth Model

void SetMaterialModel(const std::vector<std::string> &s)

Set Material Model

void SetIceCapTypeString(std::string s)

Set IceCap type

  • ”NoIce” … no ice at all

  • ”IceSheet” … use ice sheet wrapps entirely the Earth

  • ”SimpleIceCap” … use simple spherical dorm ice

void SetIceCapSimpleAngle(double cap_angle)

Set open angle of ice cap

void SetDetectorDepth(double d)

Set detector depth

void SetDetectorXY(double x, double y)

Set detector XY

inline const double GetDetectorDepth() const

Get Detector depth.

inline const I3Position &GetDetectorPosInEarthCoord() const

Get Detector coordinate origin with respect to Earth Center Coordinate

inline const std::string &GetPath() const

get path of parameter files

inline const std::string &GetIceCapTypeString() const

Get open angle of ice cap

inline const double GetIceCapSimpleAngle() const

Get open angle of ice cap

const double GetBoundary(const std::string &s) const

Get Boundary [m]

inline const double GetMohoBoundary() const

Get MohoBoundary [m]

inline const double GetRockIceBoundary() const

Get Rock-Ice boundary [m]

inline const double GetIceAirBoundary() const

Get Ice-Air boundary [m] If you have simple cap ice, this value is most far Ice-Air boundary from the center of the Earth.

inline const double GetAtmoRadius() const

Get Radius of atmosphare [m]

const double RadiusToCosZen(double r) const

Convert radius to coszen of tangent line of the radius. This function always return 1 if the radius exceeds distance from center of the Earth to IceCube center

void Init()

init function

Public Static Functions

static std::string GetMediumTypeString(MediumType m)
static MediumType ConvertMediumTypeString(const std::string &s)
static double GetEarthDensityInCGS(const EarthParam &ep, const I3Position &p_CE)

Computes the material density in the given layer.

This function is useful when multiple density queries are needed which are known to all be within the same layer, since the layer can be cached and reused for each query. It is the user’s responsibility to ensure that ep is actually the correct material layer for p_CE (most likey by ensuring that it is the result of a call to GetEarthParam(p_CE)).

Parameters:

  • ep – the material layer in which the density should be computed

  • p_CE – the position at which the density is to be evaluated. This must be in Earth-centered coordinates.

Returns:

the density in g/cm^3

Private Functions

const double GetLeptonRangeInMeter(double energy, const I3Position &posI3, const I3Direction &dirI3, bool isTau = false, bool isReverse = false, EarthModelCalculator::LeptonRangeOption opt = EarthModelCalculator::DEFAULT, double scale = 1.0) const

This function converts MuonRange in m.w.e to m with given earth model and geometry of the track.

GetLeptonRangeInMeter

Parameters:

  • energy – [in] particle energy

  • posI3 – [in] anchor point of the range which is normally the beginning, but will be the end if isReverse is set to true.

  • dirI3 – [in] the direction the particle is traveling.

  • isTau – [in] if set the lepton will be treated as a tau, otherwsie as a muon.

  • isReverse – [in] if true, it set pos as end point and calculate m.w.e or length in opposite direction of the track. If you want to know how far from a muon can reach to the detector, set pos as the entrance of icecube and set isReverse = true.

  • opt, scale – [in] Used by EarthModelCalculator::GetLeptonRange() function. Leave it as defaut unless you understand well about the parameters.

Returns:

distance[m]

const double FlatDepthCalculator(const I3Position &frompos_CE, const I3Position &topos_CE, double density, IntegType intg_type) const

A utility function to get depth for integration.

void GetAZ(int pdgcode, int &np, int &nn)

convert material pdg code to number of protons and neutrons

bool CheckIntersectionWithLayer(const I3Position &posCE, const I3Direction &dirCE, EarthParamMap::const_iterator ep, EarthParamMap::const_iterator &closestBoundary, double &closestDist, bool &willExit) const

Computes the next point at which a track intersects the given material layer. This function exists as a helper for DistanceToNextBoundaryCrossing, and should probably not be used from other places.

Parameters:

  • posCE – the starting point of the track in EarthCenter coordinate

  • dirCE – the direction of the track in EarthCenter coordinate

  • ep – the layer to check for intersection

  • closestBoundary – [out] if the distance to ep is smaller than closestDist ep will be copied to this variable

  • closestDist – [out] the closest boundary crossing yet found, which will be updated if the crossing point for this boundary is closer

  • willExit – [out] whether the next crossing will involve the track leaving this layer (ep)

SET_LOGGER ("EarthModelService")

Private Members

std::string fPath_

path to data file

std::vector<std::string> fEarthModelStrings_

name of EarthModel data file

std::vector<std::string> fMatRatioStrings_

name of Material Ratio Map data file

double fMohoBoundary_

name of EarthModel data file

double fRockIceBoundary_
double fIceAirBoundary_
double fAtmoRadius_
double fDetDepth_
I3Position fDetPos_
EarthModelService::IceCapType fIceCapType_
std::string fIceCapTypeString_
double fIceCapSimpleAngle_
double fIceCapSimpleRadius_
double fIceCapSimpleZshift_
EarthParamMap fEarthParams_
EarthParamMap fIceParams_
MatRatioMap fMatRatioMap_

map of materials and their WEIGHT ratio in unit volume. used by Genie format : map<MediumType, map<int(pdg nuclear code), double> > example : [Ice] [1000080160][0.8881016] // O [1000010010][0.1118983] // H2

MatRatioMap fPNERatioMap_

map of light particles and their NUMBER ratio in unit volume used by NuGen or may be by oscillation calculation format : map<MediumType, map<int(pdg code), double> > example : [Ice] [2212][0.5555555555] // proton [2112][0.4444444444] // neutron [11][0.5555555555] // electron, the weight is called Ye value This map will be filled automatically when you set material file.

Private Static Attributes

static bool ignored_bool

A dummy variable used as a sink for ignored results.

THis variable’s value should be considered undefined and it should never be read.

class EarthParam
#include <EarthModelService.h>

A representation of one spherical shell of material.

Public Functions

inline EarthParam()
inline EarthParam(const EarthParam &p)
inline const double GetDensity(double x) const

Evaluates the density within this shell

Parameters:

x – the radial coordinate in meters

Returns:

the material density in g/cm^3

inline std::string PrintDensity() const

print density in density file format

Public Members

double fUpperRadius_

The outer radius of this shell.

Note that the inner radius is not stored here; it is either zero or the outer radius of the next shell inwards.

double fZOffset_

The offset (towards to south pole) of the center of this shell from the center of the earth.

std::string fBoundaryName_

The name of the upper boundary of this shell

MediumType fMediumType_

The type of material from which this shell is composed

std::vector<double> fParams_

The coefficients of the polynomial which gives the density as a function of radius in this shell

template<typename FuncType>
struct TrapezoidalIntegrator
#include <EarthModelCalculator.h>

An object which performs 1-D intgration of a function using the trapezoid rule, and allows the approximation to be refined incrementally.

Public Functions

inline TrapezoidalIntegrator(const FuncType &f, double a, double b)
Parameters:

  • f – the function to be integrated

  • a – the lower bound of integration

  • b – the upper bound of integration

inline double integrate(unsigned int detail)

Get the integral approximation, updating it with higher detail if necessary

Parameters:

detail – how finely to approximate the integral. A detail level of n requires 1+2^n function evaluations, but reuses any evaluations already performed when lower detail levels were calculated.

inline unsigned int getDetail() const

Get the current detail level of the ingral approximation.

Private Functions

inline void update()

Add one level of detail to the integral approximation

Private Members

const FuncType &f

The function being integrated

double a

The lower bound of integration

double b

The upper bound of integration

unsigned int currentDetail

Counter expressing the number of times the integral approximation has been refined

double value

The current approximation of the integral

namespace earthmodel

This is a namespace which provides a collection of stand-alone functions that calculate various geometrical information for particle propagation of the Earth.

This is a namespace which provides constants of particle like muon mass etc.

Functions

I3_POINTER_TYPEDEFS(EarthModelService)

Variables

static double TOLERANCE = 1e-6
namespace EarthModelCalculator

Enums

enum LeptonRangeOption

Values:

enumerator DEFAULT
enumerator LEGACY
enumerator NUSIM

Functions

double GetImpactParameter(const I3Position &p0, const I3Direction &d, double &t, I3Position &p)

calculate impact parameter with rewpect to origin and scalar value t that fulfills p = p0 + t*d, where p0 is start position of a track, d is direction of the track, and p is most closest position on a track from origin. this function should work in any Cartecian coordinate

Parameters:

  • p0 – [in] particle start position

  • d – [in] particle direction (unit vector)

  • t – [out] distance from p0 to p

  • p – [out] most closest position on a track from origin

Returns:

impact parameter, distance from origin to p

int GetIntersectionsWithSphere(const I3Position &pos, const I3Direction &dir, double r, I3Position &startPos, I3Position &endPos)

This function returns intersection-positions between a track and a sphere with radius r. Note that the origin of track position and direction must be at the center of the sphere. Return positions are in same coordinate as input. If there is only one intersection, it will be stored in both output parameters.

Parameters:

  • pos – [in] track position

  • dir – [in] track direction (unit vector)

  • r – [in] radius

  • startPos – [out] track-entering-position to the sphere

  • endPos – [out] track-exitting-position from the sphere

Returns:

number of intersections

int GetDistsToIntersectionsWithSphere(const I3Position &pos, const I3Direction &dir, double r, double &enterdist, double &exitdist)

wrapper function of GetIntersectionsWithSphere

Parameters:

  • pos – [in] track position

  • dir – [in] track direction (unit vector)

  • r – [in] radius

  • enterdist – [out] distance from pos to track-entering-position negative value indicates behind the pos

  • exitdist – [out] distance from pos to track-exiting-position negative value indicates behind the pos

Returns:

number of intersections

double GetLeptonRange(double particle_energy, bool isTau = false, LeptonRangeOption option = DEFAULT, double scale = 1)

Returns muon range in m.w.e. If you need surviving length [m] of muon/tau, use EarthModelService::GetLeptonRangeInMeter().

Now MuonRange calculation offers three options:

DEFAULT -> use Dima’s fitting function and optimized parameter confirmed on Mar. 29, 2011

LEGACY -> use legacy nugen equation and parameter obtained study with mmc, using ice medium

NUSIM -> use Gary’s fitting function (used in NUSIM)

scale gives a scaling factor for MuonRange in Meter. See GetLeptonRangeInMeter().

[Dima’s equation]

muon and tau: R = (1/b)*log[(bE/a)+1] for the equation see arXiv:hep-ph/0407075, section 6.2 P16.

muon: (legacy nugen parameter) b~3.4*10^-4 [1/m] a~2*10^-1 [GeV/m] R_mu ~ 3000 * log(1.0 + 1.7*10^-3*E[GeV]) [mwe]

tau: (see: https://web.archive.org/web/20100630110917/http://www.ice.phys.psu.edu/~toale/icecube/nutau/nutau-mc.html) b~2.6*10^-5 [1/m] a~1.5*10^3 [GeV/m] R_tau ~ 38000 * log(1.0 + 1.8*10^-8*E[GeV]) [mwe]

[Gary’s equation] (muon only)

if (energy < 1e3) { R_mu = energy/0.002; } else if (energy < 1e4) { R_mu = 1.0e5 * (3.5 + 9.0 * (TMath::Log10(energy) - 3.0)); } else { R_mu = 1.0e5 * (12.0 + 6.0 * (TMath::Log10(energy) - 4.0)); }

Returns:

range [m.w.e]

double ColumnDepthCGStoMWE(double cdep_CGS)

unit conversion from g/cm2 to m.w.e

double MWEtoColumnDepthCGS(double range_MWE)

unit conversion from m.w.e to g/cm2

namespace Integration

Functions

template<typename FuncType>
double rombergIntegrate(const FuncType &func, double a, double b, double tol = 1e-6)

Performs a fifth order Romberg integration of a function to a chosen tolerance.

This routine is rather simplistic and not suitable for complicated functions, particularly not ones with discontinuities, but it is very fast for smooth functions.

Parameters:

  • func – the function to integrate

  • a – the lower bound of integration

  • b – the upper bound of integration

  • tol – the (absolute) tolerance on the error of the integral

namespace ParticleConstants

Variables

static const double m_mu = 1.05658E-1 * I3Units::GeV

muon rest mass

static const double m_tau = 1.777 * I3Units::GeV

tau rest mass

namespace std

STL namespace.

file EarthModelCalculator.cxx
#include “earthmodel-service/EarthModelCalculator.h”
#include “icetray/I3Units.h”
#include “dataclasses/I3Constants.h”
file EarthModelCalculator.h
#include <cmath>
#include “dataclasses/I3Position.h”
#include “dataclasses/I3Direction.h”

copyright (C) 2004 the icecube collaboration

Rcs

EarthModelCalculator.h

Version

Rcs

1.16

Date

Rcs

2012-03-13 10:40:52 -0500 (Tue, 13 Mar 2012)

Author

Kotoyo Hoshina hoshina@icecube.wisc.edu

file EarthModelService.cxx
#include <fstream>
#include <sstream>
#include <iostream>
#include “icetray/I3Units.h”
#include “dataclasses/I3Constants.h”
#include “earthmodel-service/EarthModelService.h”
#include “earthmodel-service/ParticleConstants.h”
#include “icetray/I3SingleServiceFactory.h”

EarthModelService manages density profile of the Earth.

Copyright (C) 2005 The IceCube Collaboration

Rcs

Version

Rcs
Date

Rcs
Author

Kotoyo Hoshina kotoyo.hoshina@icecube.wisc.edu

Typedefs

typedef I3SingleServiceFactory<EarthModelService> I3EarthModelServiceFactory

Functions

I3_SERVICE_FACTORY(I3EarthModelServiceFactory)

Variables

static const double PREM_EARTH_RADIUS = 6371000
static const int CHAR_BUF_SIZE = 8196
static const double M_TO_CM = (I3Units::m / I3Units::cm)
static const double CM_TO_M = (I3Units::cm / I3Units::m)
file EarthModelService.h
#include “icetray/I3ServiceBase.h”
#include “icetray/I3SingleServiceFactory.h”
#include “dataclasses/I3Position.h”
#include “dataclasses/I3Direction.h”
#include “earthmodel-service/EarthModelCalculator.h”
#include <sstream>
#include <string>
#include <vector>
#include <boost/foreach.hpp>
#include <boost/python.hpp>

A class for managing Earth’s density profile and ice geometry.

It also offers utility functions that may be useful for propagators.

  • GetDensityInCGS(pos) function

  • GetMedium(pos) function (currently returns ROCK or ICE)

  • Convert position & direction in Earth-centered coordinate to Detector Coordinate, and vice vasa

This class requires a crust data file (ascii format). see resources/earthparam/PREM_mmc.dat for details

CAUTION : The internal units of the module are [m] for length and [g/cm^3] for density, as in the data format mentioned above. The unit conversion is done inside the module so that you need to set density etc. with IceCube units.

(c) the IceCube Collaboration

Author

Kotoyo Hoshina (hoshina@icecube.wisc.edu)

Version

Rcs
Date

Rcs

file ParticleConstants.h
#include “icetray/I3Units.h”

copyright (C) 2004 the icecube collaboration

Rcs

Version

Rcs

1.16

Date

Rcs

2012-03-13 10:40:52 -0500 (Tue, 13 Mar 2012)

Author

K.Hoshina

dir earthmodel-service
dir earthmodel-service
dir earthmodel-service
dir icetray
dir private
dir public