OscProb
OscProb::EarthModelBinned Class Reference

Implements an earth model with depth/latitude/longitude bins. More...

#include <EarthModelBinned.h>

Inheritance diagram for OscProb::EarthModelBinned:
OscProb::EarthModelBase

Classes

struct  LatBinInfo
 A struct holding information about upcoming latitude bin crossings along the neutrino's trajectory. More...
 
struct  LonBinInfo
 A struct holding information about upcoming longitude bin crossings along the neutrino's trajectory. More...
 

Public Member Functions

 EarthModelBinned (std::string filename="")
 Constructor. More...
 
virtual ~EarthModelBinned ()
 Destructor. More...
 
void SetDetPos (double rad, double lat=0, double lon=0)
 
int FillPath (double cosT, double phi=0)
 
virtual void LoadModel (std::string filename)
 Load an earth model from a file. More...
 
virtual std::vector< EarthBinGetEarthBins ()
 
virtual void SetRegionZoA (int index, double zoa)
 Set Z/A of all bins with specified region index. More...
 
virtual double GetRegionZoA (int index)
 Get Z/A of all bins with specified region index. More...
 
virtual void ScaleRegionDensity (int index, double scalingfactor)
 
virtual std::vector< NuPathGetNuPath ()
 
virtual std::vector< NuPathGetMergedPaths (double prec=0.25)
 Get merged path sequence in a vector. More...
 
virtual double GetTotalL (double cosT)
 Get the total baseline for a given cosTheta. More...
 
virtual double GetCosT (double L)
 Get the cosTheta for a given total baseline. More...
 
virtual void SetRemoveSmallPaths (bool rp=true)
 Set tag to remove small paths. More...
 

Protected Member Functions

virtual void ClearModel ()
 Clear the earth model information. More...
 
virtual void AddBin (double radius_out, double radius_in, double latitude, double longitude, double density, double zoa, double layer)
 Add a bin to the model (angles in degrees) More...
 
virtual void AddPath (double length, EarthBin bin)
 Add a path segment to the sequence. More...
 
virtual double DetDistForNextLatBin (int cur_index, LatBinInfo &L)
 
virtual double DetDistForNextLonBin (double prev_lon, LonBinInfo &L)
 
virtual void RecordLatLonBinCrossings (double detDist_nextDbin, double &DetDist, int &index, LatBinInfo &latI, LonBinInfo &lonI)
 
virtual int LonBinIndex (double longitude)
 Find lon bin index containing longitude. More...
 
 ClassDef (EarthModelBinned, 1)
 
virtual void SetDetectorCoordinates (double rad, double lat, double lon)
 
virtual void AddPathSegment (double length, double density, double zoa, int index)
 Add a path segment to the sequence. More...
 
 ClassDef (EarthModelBase, 1)
 

Protected Attributes

TrajConstants fC
 Useful constants for the trajectory. More...
 
std::vector< EarthBinfEarthBins
 The bins in the earth model. More...
 
int fnDepthBins
 Total number of depth bins. More...
 
int fnLonBins
 Total number of longitude bins. More...
 
int fnLatBins
 Total number of latitude bins. More...
 
double fInvLonBinWidth
 1/binwidth for each longitude bin More...
 
double fInvLatBinWidth
 1/binwidth for each latitude bin More...
 
double fHalfLonBinWidth
 Half-width of each longitude bin. More...
 
double fHalfLatBinWidth
 Half-width of each latitude bin. More...
 
int fmaxRegIndex
 Largest region index in model. More...
 
std::string fErrorMessage_LonInfo
 
int fLonError
 
std::vector< NuPathfNuPath
 The current neutrino path sequence. More...
 
double fRadiusMax
 Maximum radius in Earth model (in km) More...
 
double fDetRadius
 The radius where the detector sits (in km) More...
 
double fDetLat
 The latitude (in rad) where the detector sits. More...
 
double fDetLon
 The longitude (in rad) where the detector sits. More...
 
bool fRemoveSmallPaths
 Tag whether to merge small paths. More...
 

Detailed Description

This class implements a 3D (asymmetric) model of the earth using EarthBin's to store bins with different properties.

The class is then able to produce path sequences through the Earth as a function of the azimuthal angle and cosine of the zenith angle with respect to the detector.

The detector can be positioned at any point within the Earth, and the path sequences will take into account the fact that some layers are above the detector.

Using indices to specify various regions (which must have a constant Z/A value), the Z/A value for a region can be changed, and the density throughout the region can be scaled by a constant factor.

By default this implements the model stored in EarthTables/earth_binned_default.txt with the detector at the bottom of the ocean layer (radius = 6368 km) where the prime meridian intersects the equator.

This class inherits from EarthModelBase and can be saved in ROOT files.

Author
rpestes@apc.in2p3.fr

Definition at line 201 of file EarthModelBinned.h.

Constructor & Destructor Documentation

◆ EarthModelBinned()

EarthModelBinned::EarthModelBinned ( std::string  filename = "")

Constructor.

By default this implements the model stored in EarthTables/earth_binned_default.txt with the detector 6368 km from the center of the Earth, having a latitude of 0 deg N and a longitude of 0 deg E.

Parameters
filename- The txt file containing a table of earth layers

Definition at line 121 of file EarthModelBinned.cxx.

122{
123 SetDetPos(6368, 0, 0);
124 LoadModel(filename);
125 SetRemoveSmallPaths(false);
126 fLonError = 0;
127}
virtual void SetRemoveSmallPaths(bool rp=true)
Set tag to remove small paths.
virtual void LoadModel(std::string filename)
Load an earth model from a file.
void SetDetPos(double rad, double lat=0, double lon=0)

References fLonError, LoadModel(), SetDetPos(), and OscProb::EarthModelBase::SetRemoveSmallPaths().

◆ ~EarthModelBinned()

EarthModelBinned::~EarthModelBinned ( )
virtual

Nothing to clean.

Definition at line 133 of file EarthModelBinned.cxx.

133{}

Member Function Documentation

◆ AddBin()

void EarthModelBinned::AddBin ( double  radius_out,
double  radius_in,
double  latitude,
double  longitude,
double  density,
double  zoa,
double  index 
)
protectedvirtual

Add a bin to the earth model.

Parameters
radius_out- The outer depth of the bin in km
radius_in- The inner depth of the bin in km
latitude- The latitude of the bin center in degrees
longitude- The longitude of the bin center in degrees
density- The density of the matter in the bin in g/cm^3
zoa- The effective Z/A value of the matter in the bin
index- Region index

Definition at line 187 of file EarthModelBinned.cxx.

190{
191 fEarthBins.push_back(EarthBin(radius_out, radius_in, latitude, longitude,
192 density, zoa, index));
193}
std::vector< EarthBin > fEarthBins
The bins in the earth model.

References fEarthBins.

Referenced by LoadModel().

◆ AddPath()

void EarthModelBinned::AddPath ( double  length,
EarthBin  bin 
)
protectedvirtual

Add a path segment to the sequence.

For a given EarthBin, adds a path of a given length in that material

Parameters
length- The length of the path segment in km
bin- The bin we are crossing

Definition at line 433 of file EarthModelBinned.cxx.

434{
435 AddPathSegment(length, bin.density, bin.zoa, bin.index);
436}
virtual void AddPathSegment(double length, double density, double zoa, int index)
Add a path segment to the sequence.
double zoa
The effective Z/A value of the matter in the bin.
int index
Region index.
double density
The density of the matter in the bin in g/cm^3.

References OscProb::EarthModelBase::AddPathSegment(), OscProb::EarthBin::density, OscProb::EarthBin::index, and OscProb::EarthBin::zoa.

Referenced by FillPath(), and RecordLatLonBinCrossings().

◆ AddPathSegment()

void EarthModelBase::AddPathSegment ( double  length,
double  density,
double  zoa,
int  index 
)
protectedvirtualinherited

Add a path segment to the sequence.

For a given EarthBin, adds a path of a given length in that material

Parameters
length- The length of the path segment in km
density- The density along the path segment
zoa- Z/A along the path segment
index- Index for the matter along the path segment

Definition at line 65 of file EarthModelBase.cxx.

67{
68 fNuPath.push_back(NuPath(length, density, zoa, index));
69}
std::vector< NuPath > fNuPath
The current neutrino path sequence.
A struct representing a neutrino path segment.
Definition: NuPath.h:34

References OscProb::EarthModelBase::fNuPath.

Referenced by AddPath(), and OscProb::PremModel::AddPath().

◆ ClassDef() [1/2]

OscProb::EarthModelBase::ClassDef ( EarthModelBase  ,
 
)
protectedinherited

◆ ClassDef() [2/2]

OscProb::EarthModelBinned::ClassDef ( EarthModelBinned  ,
 
)
protected

◆ ClearModel()

void EarthModelBinned::ClearModel ( )
protectedvirtual

Clear the earth model information.

Definition at line 147 of file EarthModelBinned.cxx.

148{
149 fnDepthBins = 0;
150 fnLonBins = 0;
151 fnLatBins = 0;
152 fEarthBins.clear();
153}
int fnDepthBins
Total number of depth bins.
int fnLatBins
Total number of latitude bins.
int fnLonBins
Total number of longitude bins.

References fEarthBins, fnDepthBins, fnLatBins, and fnLonBins.

Referenced by LoadModel().

◆ DetDistForNextLatBin()

double EarthModelBinned::DetDistForNextLatBin ( int  cur_index,
LatBinInfo L 
)
protectedvirtual

Calculate the detector distance at the edge of the current lat bin along the neutrino's trajectory

Calculate the distance to the detector along a neutrino trajectory, specified by cosT and the azimuthal angle, at the edge of the current latitude bin in the direction specified by dir. Returns a negative value if the next latitude bin will not be reached before the neutrino reaches the detector. Flips L.sign and changes L.maxreached to true if the extreme latitude is within the cur_index Earth bin.

In an attempt to make this calculation more robust, this function turns the latitude bin iteration around when sin^2(lat) gets within 10^(-14) of its extreme limit.

Parameters
cur_index- Index of current Earth bin
L- Information for the lat bin crossings
Returns
Calculated detector distance

Definition at line 456 of file EarthModelBinned.cxx.

457{
458 // Find next lat bin boundary (bin center +/- 1/2-binwidth)
459 double NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
460 // Turn around if lat reaches +/- 90 degrees (or close enough)
461 if (abs(NewLat) - 0.5 * M_PI > -1.0e-14) {
462 if (L.maxreached) return -1; // Only turn around once
463 L.maxreached = true;
464 // cout << "Max latitude reached..." << endl;
465 L.sign = -L.sign;
466 NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
467 if (abs(NewLat) - 0.5 * M_PI > -1.0e-14)
468 return -1; // This would happen if there's only one latitude bin...
469 }
470
471 // Calculate parts of answer
472 double sinNewLat = sin(NewLat);
473 double NsinSqNewLat = -sinNewLat * sinNewLat;
474 double radicand = fC.maxSinSqLat +
475 NsinSqNewLat; // 1-[sin(NewLat)]^2-[sin(phi)cos(DetLat)]^2
476 // Turn around if the radicand is 0 (close enough) or less
477 if (radicand < 1.0e-14) {
478 if (L.maxreached) return -1; // Only turn around once
479 L.maxreached = true;
480 // cout << "Max latitude reached..." << endl;
481 L.sign = -L.sign;
482 NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
483 sinNewLat = sin(NewLat);
484 NsinSqNewLat = -sinNewLat * sinNewLat;
485 radicand = fC.maxSinSqLat + NsinSqNewLat;
486 if (radicand < 1.0e-14)
487 return -1; // The rest of the track is contained in this lat bin
488 }
489 double denominator = fC.gammaSq + NsinSqNewLat;
490
491 // Calculate distance from the detector at next lat bin boundary
492 double answer;
493 if (denominator == 0) {
494 answer = 0.5 * (fC.xlatextreme - fC.rDetSinDetLat / fC.gamma);
495 // There is a "possible" case where answer would be 0:
496 // if the detector is at the center of the Earth... just don't...
497 }
498 else {
499 answer = -(fC.rDetCosT * NsinSqNewLat + fC.rDetGammaSinDetLat +
500 L.sign * fC.rDetSinT * sinNewLat * sqrt(radicand)) /
501 denominator;
502 }
503
504 // Check to see if going to the edge of the bin jumped to an invalid
505 // part of the function (cosA*rDetCosDetLat/beta = x at max lat)
506 if (L.maxreached && answer >= fC.xlatextreme) {
507 return -1; // The rest of the track is contained in this lat bin
508 }
509
510 return answer;
511}
TrajConstants fC
Useful constants for the trajectory.
double gammaSq
[sin(T)*cos(phi)*cos(DetLat)+cosT*sin(DetLat)]^2
double rDetSinT
rDet*sin(T)
double rDetCosT
rDet*cosT
double rDetGammaSinDetLat
rDet*sin(DetLat)*[sin(T)*cos(phi)*cos(DetLat)+cosT*sin(DetLat)]
double gamma
sin(T)*cos(phi)*cos(DetLat)+cosT*sin(DetLat)
double maxSinSqLat
1 - [sin(phi)*cos(DetLat)]^2
double rDetSinDetLat
rDet*sin(DetLat)

References OscProb::EarthModelBinned::LatBinInfo::dLat, fC, fEarthBins, OscProb::TrajConstants::gamma, OscProb::TrajConstants::gammaSq, OscProb::EarthModelBinned::LatBinInfo::maxreached, OscProb::TrajConstants::maxSinSqLat, OscProb::TrajConstants::rDetCosT, OscProb::TrajConstants::rDetGammaSinDetLat, OscProb::TrajConstants::rDetSinDetLat, OscProb::TrajConstants::rDetSinT, OscProb::EarthModelBinned::LatBinInfo::sign, and OscProb::TrajConstants::xlatextreme.

Referenced by FillPath(), and RecordLatLonBinCrossings().

◆ DetDistForNextLonBin()

double EarthModelBinned::DetDistForNextLonBin ( double  prev_lon,
LonBinInfo L 
)
protectedvirtual

Calculate the detector distance at the edge of the current lon bin along the neutrino's trajectory and increment L.bin

Calculate the distance to the detector at the edge of the current longitude bin along a neutrino trajectory specified by cosT and the azimuthal angle, setting the index for the next longitude bin in the process. The maximum/minimum longitude are expressed around the 0-2pi boundary. So, if the track spans that boundary, the actual allowed range for lon is [min_lon,2pi)U[0,max_lon].

Parameters
prev_lon- Value at center of initial longitude bin
L- Information about longitude bin crossings
Returns
Calculated detector distance

Definition at line 526 of file EarthModelBinned.cxx.

527{
528 double twoPi = 2 * M_PI;
529 // Increment to next bin in direction dLon
530 double lon = prev_lon + L.dLon;
531 if (lon < 0 || lon >= twoPi) { lon -= floor(lon / twoPi) * twoPi; }
532 if (L.dLon > 0) {
533 L.nextBin = L.bin + 1;
534 if (L.nextBin >= fnLonBins) {
535 // cross from 2pi to 0
536 L.nextBin = 0;
537 }
538 }
539 else {
540 L.nextBin = L.bin - 1;
541 if (L.nextBin < 0) {
542 // cross from 0 to 2pi
543 L.nextBin = fnLonBins - 1;
544 }
545 }
546
547 // Check if lon is outside of range covered by neutrino's trajectory
548 if (L.max < L.min) {
549 if (0 > lon - L.min && L.max - lon < 0) return -1;
550 }
551 else {
552 if (0 > lon - L.min || L.max - lon < 0) return -1;
553 }
554
555 // Parts of the answer
556 double cosNewLon = cos(lon);
557 double sinNewLon = sin(lon);
558 double expression = fC.sinDetLon * cosNewLon - sinNewLon * fC.cosDetLon;
559 double denominator = fC.alpha * expression +
560 fC.sinTsinAcosDetLon * cosNewLon +
561 fC.sinTsinAsinDetLon * sinNewLon;
562
563 // Calculate & return detector distance
564 if (denominator == 0) { // Is there ever a case when this is true???
565 // Error message will be printed at the end of the FillPath function
566 L.err_message.assign("Longitude along neutrino trajectory = " +
567 std::to_string(lon));
568 L.error = -1;
569 return -1;
570 }
571 return fC.rDetCosDetLat * expression / denominator;
572}
double cosDetLon
cos(DetLon)
double alpha
sin(T)*cos(phi)*sin(DetLat)-cosT*cos(DetLat)
double sinDetLon
sin(DetLon)
double sinTsinAsinDetLon
sin(T)*sin(phi)*sin(DetLon)
double rDetCosDetLat
rDet*cos(DetLat)
double sinTsinAcosDetLon
sin(T)*sin(phi)*cos(DetLon)

References OscProb::TrajConstants::alpha, OscProb::EarthModelBinned::LonBinInfo::bin, OscProb::TrajConstants::cosDetLon, OscProb::EarthModelBinned::LonBinInfo::dLon, OscProb::EarthModelBinned::LonBinInfo::err_message, OscProb::EarthModelBinned::LonBinInfo::error, fC, fnLonBins, OscProb::EarthModelBinned::LonBinInfo::max, OscProb::EarthModelBinned::LonBinInfo::min, OscProb::EarthModelBinned::LonBinInfo::nextBin, OscProb::TrajConstants::rDetCosDetLat, OscProb::TrajConstants::sinDetLon, OscProb::TrajConstants::sinTsinAcosDetLon, and OscProb::TrajConstants::sinTsinAsinDetLon.

Referenced by FillPath(), and RecordLatLonBinCrossings().

◆ FillPath()

int EarthModelBinned::FillPath ( double  cosT,
double  phi = 0 
)
virtual

Fill the path sequence in a vector (phi in degrees)

Fill the path sequence in a vector.

This will start at the upper-most layer (where the neutrino's path started) and find path lengths to the next bin crossed. When reaching the inner-most depth in the given direction, the paths move back to outer layers until hitting the detector. The neutrino direction is specified as the zenith angle and azimuthal angle of the vector pointing from the detector to where the neutrino entered the Earth.

The path sequence is stored as an attribute and can be retrieved with the function GetNuPath.

Note: This function assumes uniform bin sizes for the latitude bins and the longitude bins.

Parameters
cosT- The cosine of the zenith angle of the neutrino direction
phi- The azimuthal angle (in degrees from North) of the neutrino direction
Returns
The number of path segments in the sequence

Implements OscProb::EarthModelBase.

Definition at line 691 of file EarthModelBinned.cxx.

692{
693 // Clear current path sequence
694 fNuPath.clear();
695
696 // Do nothing if cosine is unphysical
697 if (fabs(cosT) > 1) return 0;
698
699 // Calculate useful combinations of variables
700 double Az = phi / 180.0 * M_PI; // convert from degrees to radians
701 fC.UpdateNuAngles(cosT, Az);
702 fC.Recalculate();
703 double rDetSqSinSqT = pow(fC.rDetSinT, 2);
704 double distEdgeToMin =
705 sqrt(fRadiusMax * fRadiusMax -
706 rDetSqSinSqT); // distance from point where neutrino enters to
707 // minimum depth (neglecting stopping at the detector)
708 double baseline = distEdgeToMin - fC.rDetCosT; // total baseline
709
710 // Define the maximum depth (minimum radius) along the path (detector
711 // depth/radius, if cosT is non-negative)
712 double minRadius = fDetRadius;
713 if (cosT < 0) { minRadius = fC.rDetSinT; }
714
715 // Starting latitude
716 /* cout << "Pieces of init_lat calculation:" << endl
717 << "\trDetSinT = " << fC.rDetSinT << endl
718 << "\tbeta = " << fC.beta << endl
719 << "\tgamma = " << fC.gamma << endl
720 << "\tdistEdgeToMin = " << distEdgeToMin << endl
721 << "\tRadiusMax = " << fRadiusMax << endl; */
722 double init_lat = asin((fC.rDetSinT * fC.beta + fC.gamma * distEdgeToMin) /
723 fRadiusMax); // in radians
724
725 // Starting longitude
726 double init_lon = atan2(
728 baseline * (fC.alpha * fC.sinDetLon + fC.sinTsinAcosDetLon),
730 baseline *
732 fC.alpha * fC.cosDetLon)); // in radians (between -M_PI and M_PI)
733 if (init_lon < 0) { init_lon += 2 * M_PI; }
734
735 // Get initial Earth bin and # of bins per depth layer
736 LatBinInfo latinfo;
737 LonBinInfo loninfo;
738 int binsPerDepth = fEarthBins.size() / fnDepthBins;
739 int init_lonBin = LonBinIndex(init_lon);
740 int init_latBin = floor((init_lat + M_PI / 2.0) * fInvLatBinWidth);
741 int index =
742 fEarthBins.size() - binsPerDepth + init_lonBin * fnLatBins + init_latBin;
743
744 latinfo.bin = init_latBin;
745 loninfo.bin = init_lonBin;
746 // cout << "Neutrino starting at (lat,lon)=(" << init_lat << "rad," <<
747 // init_lon << "rad) in Earth bin " << index << " = (r" <<
748 // floor(index/binsPerDepth) << ", lon" << loninfo.bin << ", lat" <<
749 // latinfo.bin << ")" << endl;
750
751 /* Find detector distances for 1st latitude/longitude bin changes */
752 latinfo.dLat = fHalfLatBinWidth; // change in lat from bin center to next bin
753 // (0 if no more than 1 lat bin change)
754 // Condition calculation for starting sign in x(lat) equation
755 latinfo.sign = 1; //+ => lat incr., - => lat decr. (with decr. x)
756 latinfo.maxreached =
757 false; // this changes to true if past the lat func transition
758 if (fC.beta < 0) latinfo.sign = -1;
759 if (fC.beta == 0) {
760 if (fC.cosA > 0) latinfo.sign = -1;
761 }
762 else if (baseline < fC.xlatextreme) {
763 latinfo.sign = -(latinfo.sign);
764 latinfo.maxreached = true;
765 }
766 // Condition calculation for incr./decr. Longitude (with decr. x)
767 loninfo.dLon =
768 fHalfLonBinWidth; // change in lon from bin center to next bin (includes
769 // direction; 0 if no more than 1 lon bin change)
770 loninfo.min = init_lon; // initial lon if lon is inc.
771 loninfo.max = fDetLon; // detector lon if lon is inc.
772 if (fC.sinA < 0) {
773 loninfo.dLon = -(loninfo.dLon);
774 loninfo.min = loninfo.max;
775 loninfo.max = init_lon;
776 // Switch max/min if dec., instead of inc.
777 loninfo.max = init_lon;
778 loninfo.min = fDetLon;
779 }
780 // Remainder of detDist (and nextBin) calculations
781 if (fC.sinT == 0) {
782 // both latitude and logitude flip at center of the Earth
783 loninfo.detDist_nextBin =
784 -fC.rDetCosT; // formula is divided by cosT, but cosT = +-1
785 loninfo.nextBin = LonBinIndex(fDetLon);
786 loninfo.dLon = 0;
787 latinfo.detDist_nextBin = loninfo.detDist_nextBin;
788 latinfo.nextBin = floor((fDetLat + 0.5 * M_PI) * fInvLatBinWidth);
789 latinfo.dLat = 0;
790 }
791 else {
792 if (fC.sinA == 0) {
793 if (fC.alpha == 0) { // crosses "at" infinity
794 loninfo.detDist_nextBin = -1;
795 }
796 else { // longitude flips at center of the Earth
797 loninfo.detDist_nextBin = fC.rDetCosDetLat / fC.alpha;
798 loninfo.nextBin = LonBinIndex(fDetLon);
799 }
800 loninfo.dLon = 0;
801 }
802 else if (fC.cosDetLat == 0) {
803 loninfo.detDist_nextBin = -1;
804 }
805 else {
806 // x(lon) calculation
807 loninfo.detDist_nextBin = DetDistForNextLonBin(
808 fEarthBins[index].longitude, loninfo); // also sets loninfo.nextBin
809 }
810 // x(lat) calculation
811 latinfo.detDist_nextBin = DetDistForNextLatBin(index, latinfo);
812 latinfo.nextBin = init_latBin + latinfo.sign;
813 }
814
815 // Get initial Detector Distance
816 double prev_DetDistance = baseline;
817
818 /* Start at the top layer of Earth model and go down to minimum (not including
819 * minimum), looping over depth bins */
820 int dBin = 0; // depth bin index
821 for (dBin = 0; dBin < fnDepthBins - 1; dBin++) {
822 // Check if minimum radius is contained in depth bin
823 double r_in = fEarthBins[index].radius_in; // inner-most radius
824 if (r_in < minRadius) break;
825
826 // Distance from detector at inner-most radius
827 double DetDistance = -fC.rDetCosT + sqrt(r_in * r_in - rDetSqSinSqT);
828
829 // Check for latitude/longitude bin crossings
830 RecordLatLonBinCrossings(DetDistance, prev_DetDistance, index, latinfo,
831 loninfo);
832
833 // Add Segment to next depth bin to path
834 AddPath(prev_DetDistance - DetDistance, fEarthBins[index]);
835
836 // Reset variables for next depth bin
837 prev_DetDistance = DetDistance;
838 index -= binsPerDepth;
839 // cout << "Depth bin crossing at x=" << DetDistance << "km...now in
840 // Earth bin " << index << " = (r" << floor(index/binsPerDepth) << ",
841 // lon" << loninfo.bin << ", lat" << latinfo.bin << ")" << endl;
842 }
843
844 /* Do minimum depth bin and then go up to the detector */
845 int index2 = index;
846 for (index2 = index; fEarthBins[index2].radius_out < fDetRadius;
847 index2 += binsPerDepth) {
848 double r_out = fEarthBins[index2].radius_out; // outer-most radius
849
850 // Distance from detector at outer-most radius
851 double DetDistance = -fC.rDetCosT - sqrt(r_out * r_out - rDetSqSinSqT);
852
853 // Check for latitude/longitude bin crossings
854 RecordLatLonBinCrossings(DetDistance, prev_DetDistance, index2, latinfo,
855 loninfo);
856
857 // Add Segment in depth bin to path
858 AddPath(prev_DetDistance - DetDistance, fEarthBins[index2]);
859
860 // Reset variables for "previous" depth bin
861 prev_DetDistance = DetDistance;
862 // cout << "Depth bin crossing at x=" << DetDistance << "km...now in
863 // Earth bin " << index2 << " = (r" <<
864 // floor((index2+binsPerDepth)/binsPerDepth) << ", lon" << loninfo.bin <<
865 // ", lat" << latinfo.bin << ")" << endl;
866 }
867
868 /* Do path to detector */
869 // Check for latitude/longitude bin crossings
870 RecordLatLonBinCrossings(0, prev_DetDistance, index2, latinfo, loninfo);
871 // Add the path segment within the detector bin
872 AddPath(prev_DetDistance, fEarthBins[index2]);
873
874 // Longitude calculation error message, if needed
875 if (loninfo.error < 0) {
876 cerr
877 << "ERROR: Oops... It turns out that I was wrong about the denominator"
878 << " of the x(long) equation never being 0 apart from the cases where "
879 << "sin(Zenith) = 0, sin(Azimuthal) = 0, or cos(Detector Longitude) = "
880 << "0."
881 << " Send the following message to rpestes@apc.in2p3.fr:" << endl
882 << endl
883 << "You were wrong about the denominator of x(long)... "
884 << "Here are the values of the variables that broke it:" << endl
885 << "\tDetector Coordinates (r,lat,lon) = (" << fDetRadius << ", "
886 << fDetLat << ", " << fDetLon << ")" << endl
887 << "\tcos(Zenith Angle) = " << cosT << endl
888 << "\tAzimuthal Angle = " << phi << " deg" << endl
889 << "\t" << loninfo.err_message << endl
890 << "Please fix this ASAP!" << endl;
891 return -1;
892 }
893
894 // Return the number of path segments
895 return fNuPath.size();
896}
double fRadiusMax
Maximum radius in Earth model (in km)
double fDetLat
The latitude (in rad) where the detector sits.
double fDetLon
The longitude (in rad) where the detector sits.
double fDetRadius
The radius where the detector sits (in km)
virtual double DetDistForNextLonBin(double prev_lon, LonBinInfo &L)
virtual void AddPath(double length, EarthBin bin)
Add a path segment to the sequence.
virtual void RecordLatLonBinCrossings(double detDist_nextDbin, double &DetDist, int &index, LatBinInfo &latI, LonBinInfo &lonI)
double fHalfLonBinWidth
Half-width of each longitude bin.
double fHalfLatBinWidth
Half-width of each latitude bin.
double fInvLatBinWidth
1/binwidth for each latitude bin
virtual int LonBinIndex(double longitude)
Find lon bin index containing longitude.
virtual double DetDistForNextLatBin(int cur_index, LatBinInfo &L)
double beta
sin(T)*sin(DetLat)-cos(phi)*cosT*cos(DetLat)
void UpdateNuAngles(double cosTheta, double phi)
double cosDetLat
cos(DetLat)
double sinT
sin(T) = sqrt(sinSqT)

References AddPath(), OscProb::TrajConstants::alpha, OscProb::TrajConstants::beta, OscProb::EarthModelBinned::LatBinInfo::bin, OscProb::EarthModelBinned::LonBinInfo::bin, OscProb::TrajConstants::cosA, OscProb::TrajConstants::cosDetLat, OscProb::TrajConstants::cosDetLon, OscProb::EarthModelBinned::LatBinInfo::detDist_nextBin, OscProb::EarthModelBinned::LonBinInfo::detDist_nextBin, DetDistForNextLatBin(), DetDistForNextLonBin(), OscProb::EarthModelBinned::LatBinInfo::dLat, OscProb::EarthModelBinned::LonBinInfo::dLon, OscProb::EarthModelBinned::LonBinInfo::err_message, OscProb::EarthModelBinned::LonBinInfo::error, fC, OscProb::EarthModelBase::fDetLat, OscProb::EarthModelBase::fDetLon, OscProb::EarthModelBase::fDetRadius, fEarthBins, fHalfLatBinWidth, fHalfLonBinWidth, fInvLatBinWidth, fnDepthBins, fnLatBins, OscProb::EarthModelBase::fNuPath, OscProb::EarthModelBase::fRadiusMax, OscProb::TrajConstants::gamma, LonBinIndex(), OscProb::EarthModelBinned::LonBinInfo::max, OscProb::EarthModelBinned::LatBinInfo::maxreached, OscProb::EarthModelBinned::LonBinInfo::min, OscProb::EarthModelBinned::LatBinInfo::nextBin, OscProb::EarthModelBinned::LonBinInfo::nextBin, OscProb::TrajConstants::rDetCosDetLat, OscProb::TrajConstants::rDetCosT, OscProb::TrajConstants::rDetSinT, OscProb::TrajConstants::Recalculate(), RecordLatLonBinCrossings(), OscProb::EarthModelBinned::LatBinInfo::sign, OscProb::TrajConstants::sinA, OscProb::TrajConstants::sinDetLon, OscProb::TrajConstants::sinT, OscProb::TrajConstants::sinTsinAcosDetLon, OscProb::TrajConstants::sinTsinAsinDetLon, OscProb::TrajConstants::UpdateNuAngles(), and OscProb::TrajConstants::xlatextreme.

◆ GetCosT()

double EarthModelBase::GetCosT ( double  L)
virtualinherited

Get the cosTheta for a given total baseline.

Given a baseline, find the direction of the neutrino. This could be useful for experiments with fixed baselines for example.

The baseline must be within the range of possible values in this earth model. Will return vertical neutrinos otherwise.

Parameters
L- The total baseline of the neutrino

Definition at line 99 of file EarthModelBase.cxx.

100{
101 if (L < fRadiusMax - fDetRadius) return 1;
102 if (L > fRadiusMax + fDetRadius) return -1;
103
104 return (fRadiusMax * fRadiusMax - fDetRadius * fDetRadius - L * L) /
105 (2 * fDetRadius * L);
106}

References OscProb::EarthModelBase::fDetRadius, and OscProb::EarthModelBase::fRadiusMax.

◆ GetEarthBins()

vector< EarthBin > EarthModelBinned::GetEarthBins ( )
virtual

Get the set of earth layers

Get the set of Earth bins

This returns the set of bins for this Earth model.

Definition at line 141 of file EarthModelBinned.cxx.

141{ return fEarthBins; }

References fEarthBins.

◆ GetMergedPaths()

vector< NuPath > EarthModelBase::GetMergedPaths ( double  prec = 0.25)
virtualinherited

Merge similar paths to reduce number of steps

This method will merge consecutive paths and take their averages until it finds a large enough gap to start a new merged path.

The merged paths will be returned, and the original detailed path will not be changed and will stay stored as an attribute.

Parameters
prec- The precision to merge paths in g/cm^3
Returns
The vector of merged path segments

Definition at line 121 of file EarthModelBase.cxx.

122{
123 // The output vector
124 vector<NuPath> mergedPath;
125
126 // Start with the first path
127 NuPath path = fNuPath[0];
128
129 // Track the total length
130 double totL = 0;
131
132 // Loop over all paths starting from second
133 for (int i = 1; i < fNuPath.size(); i++) {
134 // If this path electron density is beyond the tolerance
135 if (fabs(path.density * path.zoa - fNuPath[i].density * fNuPath[i].zoa) >
136 prec * path.zoa) {
137 // Add merged path to vector
138 mergedPath.push_back(path);
139
140 // Set this path as new merged path
141 path = fNuPath[i];
142 }
143 else { // If path is within tolerance
144
145 // Merge the path with current merged path
146 path = AvgPath(path, fNuPath[i]);
147 }
148
149 // Increment total length
150 totL += fNuPath[i].length;
151
152 } // End of loop over paths
153
154 // Add the final merged path to vector
155 mergedPath.push_back(path);
156
157 // If tag is true, remove small paths
158 if (fRemoveSmallPaths) {
159 // Start at first path
160 int k = 0;
161
162 // While not at the end of vector
163 while (k + 1 < mergedPath.size()) {
164 // If length is less than 1% of total
165 if (mergedPath[k].length < 0.01 * totL) {
166 // Merge path with following path
167 mergedPath = MergePaths(mergedPath, k, k + 1);
168 }
169 // If path is long enough skip it
170 else
171 k++;
172
173 } // End of while loop
174
175 } // End of if statement
176
177 // return the merged vector
178 return mergedPath;
179}
bool fRemoveSmallPaths
Tag whether to merge small paths.
std::vector< NuPath > MergePaths(std::vector< NuPath > &inputPath, int j, int k)
Merge paths j and k in vector.
Definition: NuPath.cxx:89
NuPath AvgPath(NuPath &p1, NuPath &p2)
Get the average of two paths.
Definition: NuPath.cxx:27
double density
The density of the path segment in g/cm^3.
Definition: NuPath.h:79
double zoa
The effective Z/A value of the path segment.
Definition: NuPath.h:80

References OscProb::AvgPath(), OscProb::NuPath::density, OscProb::EarthModelBase::fNuPath, OscProb::EarthModelBase::fRemoveSmallPaths, OscProb::MergePaths(), and OscProb::NuPath::zoa.

◆ GetNuPath()

vector< NuPath > EarthModelBase::GetNuPath ( )
virtualinherited

Get the current neutrino path sequence

Get the current neutrino path sequence

The path needs to be filled for a given cosTheta before calling this function.

Definition at line 52 of file EarthModelBase.cxx.

52{ return fNuPath; }

References OscProb::EarthModelBase::fNuPath.

Referenced by OscProb::PMNS_Avg::AvgAlgo(), OscProb::PMNS_Avg::AvgAlgoCosT(), and OscProb::PMNS_Avg::ExtrapolationProbCosT().

◆ GetRegionZoA()

double EarthModelBinned::GetRegionZoA ( int  index)
virtual

Get the effective Z/A value for all bins in region specified by index, e.g. all outer-core layers. (Assumes that all bins of given index have same Z/A value)

Parameters
index- The region index
Returns
Z/A corresponding to index

Definition at line 372 of file EarthModelBinned.cxx.

373{
374 if (index > fmaxRegIndex) {
375 cerr << "ERROR: Index " << index << " does not exist in model "
376 << "(max index = " << fmaxRegIndex << "). Returning 0" << endl;
377 return 0;
378 }
379
380 int nbins = fEarthBins.size();
381
382 for (int i = 0; i < nbins; i++) {
383 if (fEarthBins[i].index != index)
384 continue;
385
386 else
387 return fEarthBins[i].zoa;
388 }
389
390 // End of vector reached without finding index
391 cerr << "ERROR: Index " << index << " not found! Returning 0" << endl;
392 return 0;
393}
int fmaxRegIndex
Largest region index in model.

References fEarthBins, and fmaxRegIndex.

◆ GetTotalL()

double EarthModelBase::GetTotalL ( double  cosT)
virtualinherited

Get the total baseline for a given cosTheta.

Parameters
cosT- The cosine of the neutrino direction

Definition at line 77 of file EarthModelBase.cxx.

78{
79 if (fabs(cosT) > 1) return 0;
80
81 double sinsqrT = 1 - cosT * cosT;
82
83 return -fDetRadius * cosT +
84 sqrt(fRadiusMax * fRadiusMax - fDetRadius * fDetRadius * sinsqrT);
85}

References OscProb::EarthModelBase::fDetRadius, and OscProb::EarthModelBase::fRadiusMax.

◆ LoadModel()

void EarthModelBinned::LoadModel ( std::string  filename)
virtual

Load an earth model from a file.

By default it loads the model stored in EarthTables/earth_binned_default.txt

The row format for the model table needs to be: longitude (deg) latitude (deg) outer depth (km) density (g/cm^3?) Z/A region_index The data in the table must be in order of decreasing depth, followed by increasing longitude, and then increasing latitude. Latitude bin widths and longitude bin widths must each be constant. In each bin of one coordinate, the bins centers for the other coordinates must be the same. This assumes that the center of the Earth is at a depth of 6371km.

Parameters
filename- The txt file containing a table of earth layers

Definition at line 214 of file EarthModelBinned.cxx.

215{
216 // Clear the current model
217 ClearModel();
218 double EarthRadius = 6371.0;
219
220 // Use default if no file provided
221 if (filename == "") { filename = PREM3D_DEFAULT; }
222 cout << "Loading Earth model table from " << filename << "..." << endl;
223
224 // Open the file
225 ifstream fin;
226 fin.open(filename.c_str());
227 if (!fin) {
228 cerr << "ERROR: File " << filename << " not found!" << endl;
229 return;
230 }
231
232 // Variables for storing table rows
233 float depth, latitude, longitude, density, zoa, index;
234
235 // Keep track of previous depth/latitude/longitude
236 double radius_prev = 0;
237 double lat_prev = -90.0;
238 double lon_prev = 0;
239
240 // Hold onto previous different depth, so we have depth maximum for bin
241 double radius_min = radius_prev;
242
243 // Initialize max region index
244 fmaxRegIndex = -1;
245
246 // Loop over table rows
247 while (fin >> longitude >> latitude >> depth >> density >> zoa >> index) {
248 double outer_radius = EarthRadius - depth;
249
250 // Move minimum depth (and reset latitude/longitude) if depth has changed
251 // from previous one
252 if (outer_radius > radius_prev) {
253 radius_min = radius_prev;
254 lat_prev = -90.0;
255 lon_prev = 0.0;
256 fnDepthBins++;
257 }
258 else if (outer_radius <
259 radius_prev) { // Depths must be ordered in model file
260 cerr << "ERROR: Depths are not sorted in decreasing order in the model "
261 "file "
262 << "(or depth is greater than " << EarthRadius << " km)." << endl;
263 ClearModel();
264 return;
265 }
266 else if (longitude > lon_prev) { // Reset latitude if longitude has changed
267 // from previous one
268 lat_prev = -90.0;
269 fnLonBins++;
270 }
271 else if (longitude < lon_prev) { // Longitudes must be ordered in model file
272 // (after depths)
273 cerr << "ERROR: Longitudes are not sorted in increasing order (after "
274 "depths) in the model file."
275 << endl;
276 ClearModel();
277 return;
278 }
279 else if (latitude < lat_prev) { // Latitudes must be ordered in model file
280 // (after longitudes)
281 cerr << "ERROR: Latitudes are not sorted in increasing order (after "
282 "longitudes) in the model file."
283 << endl;
284 ClearModel();
285 return;
286 }
287
288 // Add this bin to the model
289 AddBin(outer_radius, radius_min, latitude, longitude, density, zoa, index);
290
291 // Check region index
292 if (index > fmaxRegIndex) { fmaxRegIndex = index; }
293
294 // Set previous coordinates for next bin
295 radius_prev = outer_radius;
296 lat_prev = latitude;
297 lon_prev = longitude;
298 }
299
300 // Set Outermost Radius
301 fRadiusMax = radius_prev;
302
303 // Modify number of longitude bins to how many there actually are (instead of
304 // how many times we changed the longitude bin without changing the depth bin
305 // = nLonBins-1 times per depth bin)
307 fInvLonBinWidth = 0.5 * fnLonBins / M_PI; // in radians^(-1)
308 fHalfLonBinWidth = M_PI / fnLonBins; // in radians
309
310 // Calculate number of latitude bins
312 fInvLatBinWidth = fnLatBins / M_PI; // in radians^(-1)
313 fHalfLatBinWidth = 0.5 * M_PI / fnLatBins; // in radians
314
315 // Error if first latitude bin center deviates from about half of the
316 // calculated latitude bin width
317 if (abs(fHalfLatBinWidth - (fEarthBins[0].latitude + 0.5 * M_PI)) >
318 fHalfLatBinWidth * 0.01) { // Warn if first latitude bin center isn't
319 // about half of the latitude bin width
320 cerr << "ERROR: Latitude of 1st bin center (" << fEarthBins[0].latitude
321 << ") is not in the middle of the first bin, whose size should be "
322 << 2 * fHalfLatBinWidth << " rad, as calculated from having "
323 << fnLatBins << " latitude bins." << endl;
324 ClearModel();
325 return;
326 }
327
328 cout << "\t...done (" << fEarthBins.size() << " bins: " << fnDepthBins
329 << " depth, " << fnLatBins << " latitude, and " << fnLonBins
330 << " longitude)." << endl;
331}
virtual void ClearModel()
Clear the earth model information.
double fInvLonBinWidth
1/binwidth for each longitude bin
virtual void AddBin(double radius_out, double radius_in, double latitude, double longitude, double density, double zoa, double layer)
Add a bin to the model (angles in degrees)

References AddBin(), ClearModel(), fEarthBins, fHalfLatBinWidth, fHalfLonBinWidth, fInvLatBinWidth, fInvLonBinWidth, fmaxRegIndex, fnDepthBins, fnLatBins, fnLonBins, and OscProb::EarthModelBase::fRadiusMax.

Referenced by EarthModelBinned().

◆ LonBinIndex()

int EarthModelBinned::LonBinIndex ( double  longitude)
protectedvirtual

Find and return longitude bin index that contains longitude.

Parameters
longitude- longitude for which bin is found

Definition at line 655 of file EarthModelBinned.cxx.

656{
657 // Find difference between longitude and beginning of first lon bin
658 double londiff = longitude - (fEarthBins[0].longitude - fHalfLonBinWidth);
659
660 // Shift londiff to be between 0 and 2pi
661 if (londiff < 0 || londiff >= 2 * M_PI)
662 londiff -= floor(londiff / (2 * M_PI));
663
664 // Number of binwidths after from first lon bin
665 return floor(londiff * fInvLonBinWidth);
666}

References fEarthBins, fHalfLonBinWidth, and fInvLonBinWidth.

Referenced by FillPath().

◆ RecordLatLonBinCrossings()

void EarthModelBinned::RecordLatLonBinCrossings ( double  detDist_nextDbin,
double &  DetDist,
int &  index,
LatBinInfo latI,
LonBinInfo lonI 
)
protectedvirtual

Record path segments for each latitude/longitude bin crossed before reaching detDist_nextDbin

Record all the path segments as the neutrino crosses into new latitude and longitude bins until it reaches the next depth bin.

Parameters
detDist_nextDbin- Distance from the detector at which neutrino passes into next depth bin
DetDist- Current distance from the detector
index- Index of current Earth bin
latI- Information about lat bin crossings
lonI- Information about lon bin crossings

Definition at line 586 of file EarthModelBinned.cxx.

590{
591 // Keep going as long as there is a lat or lon bin crossing before
592 // the next depth bin crossing
593 while (detDist_nextDbin < latI.detDist_nextBin ||
594 detDist_nextDbin < lonI.detDist_nextBin) {
595 if (latI.detDist_nextBin >
596 lonI.detDist_nextBin) { // if next lat bin crossing is before next lon
597 // bin crossing
598 // Add segment to next lat bin to path
599 AddPath(DetDist - latI.detDist_nextBin, fEarthBins[index]);
600 double latitude_crossing =
601 fEarthBins[index].latitude + latI.sign * latI.dLat;
602
603 // Move lat bins
604 DetDist = latI.detDist_nextBin;
605 index += latI.nextBin - latI.bin;
606 latI.bin = latI.nextBin;
607 // cout << "Latitude bin crossing at x=" << DetDist << "km,lat=" <<
608 // latitude_crossing << "rad...now in Earth bin " << index << " = (r"
609 // << floor(index/(fnLonBins*fnLatBins)) << ", lon" << lonI.bin << ",
610 // lat" << latI.bin << ")" << endl;
611
612 // Find next lat bin
613 if (latI.dLat ==
614 0) { // Only 1 lat bin crossing along entire neutrino trajectory
615 latI.detDist_nextBin = -1;
616 }
617 else {
618 latI.detDist_nextBin = DetDistForNextLatBin(index, latI);
619 latI.nextBin = latI.bin + latI.sign;
620 }
621 }
622 else { // if next lon bin crossing is before (or simultaneous with) next lat
623 // bin corssing
624 // Add segment to next lon bin to path
625 AddPath(DetDist - lonI.detDist_nextBin, fEarthBins[index]);
626
627 // Move lon bins
628 DetDist = lonI.detDist_nextBin;
629 index += (lonI.nextBin - lonI.bin) * fnLatBins;
630 lonI.bin = lonI.nextBin;
631 // cout << "Longitude bin crossing at x=" << DetDist << "km...now in
632 // Earth bin " << index << " = (r" <<
633 // floor(index/(fnLonBins*fnLatBins)) << ", lon" << lonI.bin << ",
634 // lat" << latI.bin << ")" << endl;
635
636 // Find next lon bin
637 if (lonI.dLon ==
638 0) { // Only 1 lon bin crossing along entire neutrino trajectory
639 lonI.detDist_nextBin = -1;
640 }
641 else {
642 lonI.detDist_nextBin =
643 DetDistForNextLonBin(fEarthBins[index].longitude, lonI);
644 }
645 }
646 }
647}

References AddPath(), OscProb::EarthModelBinned::LatBinInfo::bin, OscProb::EarthModelBinned::LonBinInfo::bin, OscProb::EarthModelBinned::LatBinInfo::detDist_nextBin, OscProb::EarthModelBinned::LonBinInfo::detDist_nextBin, DetDistForNextLatBin(), DetDistForNextLonBin(), OscProb::EarthModelBinned::LatBinInfo::dLat, OscProb::EarthModelBinned::LonBinInfo::dLon, fEarthBins, fnLatBins, OscProb::EarthModelBinned::LatBinInfo::nextBin, OscProb::EarthModelBinned::LonBinInfo::nextBin, and OscProb::EarthModelBinned::LatBinInfo::sign.

Referenced by FillPath().

◆ ScaleRegionDensity()

void EarthModelBinned::ScaleRegionDensity ( int  index,
double  factor 
)
virtual

Set Z/A of all bins with specified region index

Scale density by scaling factor for all bins in region specified by given index.

Parameters
index- The region index
factor- The scaling factor for changing the density

Definition at line 403 of file EarthModelBinned.cxx.

404{
405 if (index > fmaxRegIndex) {
406 cerr << "WARNING: Index " << index << " does not exist in model "
407 << "(max index = " << fmaxRegIndex << ")." << endl
408 << "Doing nothing." << endl;
409 return;
410 }
411
412 int nbins = fEarthBins.size();
413
414 // Loop over all bins and change the ones
415 // with the given index
416 for (int i = 0; i < nbins; i++) {
417 if (fEarthBins[i].index != index) continue;
418
419 double density_init = fEarthBins[i].density;
420 fEarthBins[i].density = factor * density_init;
421 }
422}

References fEarthBins, and fmaxRegIndex.

◆ SetDetectorCoordinates()

void EarthModelBase::SetDetectorCoordinates ( double  rad,
double  lat,
double  lon 
)
protectedvirtualinherited

Set the coordinates of the detector (rad = radius in km, lat/lon in deg)

Set the coordinates of the detector: radius in km, latitude in degrees, longitude in degrees

Parameters
rad- The distance from the detector to the Earth's center in km
lat- The latitude of the detector in deg N (between -90 and 90)
lon- The longitude of the detector in deg E (between 0 and 360)

Definition at line 25 of file EarthModelBase.cxx.

26{
27 // Force radius to be non-negative
28 if (rad < 0) {
29 cerr << "WARNING: Negative radius detected. Setting to absolute value."
30 << endl;
31 rad = -rad;
32 }
33 fDetRadius = rad;
34
35 // Force latitude to be between -90 and 90 deg
36 lat -= floor((lat + 90.0) / 360.0) * 360.0;
37 if (lat > 90) { lat = 180.0 - lat; }
38 fDetLat = lat / 180.0 * M_PI; // convert to radians
39
40 // Force longitude to be between 0 and 360 deg
41 lon -= floor(lon / 360.0) * 360.0;
42 fDetLon = lon / 180.0 * M_PI; // convert to radians
43}

References OscProb::EarthModelBase::fDetLat, OscProb::EarthModelBase::fDetLon, and OscProb::EarthModelBase::fDetRadius.

Referenced by SetDetPos(), and OscProb::PremModel::SetDetPos().

◆ SetDetPos()

void EarthModelBinned::SetDetPos ( double  rad,
double  lat = 0,
double  lon = 0 
)
virtual

Set the detector position (rad = radius in km, lat/lon in deg)

Set the coordinates of the detector: radius in km, latitude in degrees, longitude in degrees Also, updates the detector-coordinate-dependent parts of the trajectory constants.

Note: If lat = -90 (or 90), lon (or -lon) defines where the azimuthal angle for the neutrino's trajectory is 0.

Parameters
rad- The distance from the detector to the Earth's center in km
lat- The latitude of the detector in deg N (between -90 and 90)
lon- The longitude of the detector in deg E (between 0 and 360)

Implements OscProb::EarthModelBase.

Definition at line 169 of file EarthModelBinned.cxx.

170{
171 SetDetectorCoordinates(rad, lat, lon);
173}
virtual void SetDetectorCoordinates(double rad, double lat, double lon)
void UpdateDetPos(double rDet, double DetLat, double DetLon)

References fC, OscProb::EarthModelBase::fDetLat, OscProb::EarthModelBase::fDetLon, OscProb::EarthModelBase::fDetRadius, OscProb::EarthModelBase::SetDetectorCoordinates(), and OscProb::TrajConstants::UpdateDetPos().

Referenced by EarthModelBinned().

◆ SetRegionZoA()

void EarthModelBinned::SetRegionZoA ( int  index,
double  zoa 
)
virtual

Set the effective Z/A value for all bins with given region index.

Use this to change the Z/A of indexed bin, e.g. all outer-core layers

Parameters
index- The region index
zoa- The effective Z/A value to use

Definition at line 343 of file EarthModelBinned.cxx.

344{
345 if (index > fmaxRegIndex) {
346 cerr << "WARNING: Index " << index << " does not exist in model "
347 << "(max index = " << fmaxRegIndex << ")." << endl
348 << "Doing nothing." << endl;
349 return;
350 }
351
352 int nbins = fEarthBins.size();
353
354 // Loop over all bins and change the ones
355 // with the given index
356 for (int i = 0; i < nbins; i++) {
357 if (fEarthBins[i].index != index) continue;
358
359 fEarthBins[i].zoa = zoa;
360 }
361}

References fEarthBins, and fmaxRegIndex.

◆ SetRemoveSmallPaths()

void EarthModelBase::SetRemoveSmallPaths ( bool  rp = true)
virtualinherited

Set the boolean to tag whether to remove small paths when merging Small is defined as <1% of the total baseline

Parameters
rp- Boolean value to set

Definition at line 188 of file EarthModelBase.cxx.

188{ fRemoveSmallPaths = rp; }

References OscProb::EarthModelBase::fRemoveSmallPaths.

Referenced by EarthModelBinned(), and OscProb::PremModel::PremModel().

Member Data Documentation

◆ fC

TrajConstants OscProb::EarthModelBinned::fC
protected

◆ fDetLat

double OscProb::EarthModelBase::fDetLat
protectedinherited

◆ fDetLon

double OscProb::EarthModelBase::fDetLon
protectedinherited

◆ fDetRadius

◆ fEarthBins

std::vector<EarthBin> OscProb::EarthModelBinned::fEarthBins
protected

◆ fErrorMessage_LonInfo

std::string OscProb::EarthModelBinned::fErrorMessage_LonInfo
protected

Part of error message containing the longitude information

Definition at line 318 of file EarthModelBinned.h.

◆ fHalfLatBinWidth

double OscProb::EarthModelBinned::fHalfLatBinWidth
protected

Definition at line 313 of file EarthModelBinned.h.

Referenced by FillPath(), and LoadModel().

◆ fHalfLonBinWidth

double OscProb::EarthModelBinned::fHalfLonBinWidth
protected

Definition at line 312 of file EarthModelBinned.h.

Referenced by FillPath(), LoadModel(), and LonBinIndex().

◆ fInvLatBinWidth

double OscProb::EarthModelBinned::fInvLatBinWidth
protected

Definition at line 311 of file EarthModelBinned.h.

Referenced by FillPath(), and LoadModel().

◆ fInvLonBinWidth

double OscProb::EarthModelBinned::fInvLonBinWidth
protected

Definition at line 310 of file EarthModelBinned.h.

Referenced by LoadModel(), and LonBinIndex().

◆ fLonError

int OscProb::EarthModelBinned::fLonError
protected

Definition at line 320 of file EarthModelBinned.h.

Referenced by EarthModelBinned().

◆ fmaxRegIndex

int OscProb::EarthModelBinned::fmaxRegIndex
protected

Definition at line 315 of file EarthModelBinned.h.

Referenced by GetRegionZoA(), LoadModel(), ScaleRegionDensity(), and SetRegionZoA().

◆ fnDepthBins

int OscProb::EarthModelBinned::fnDepthBins
protected

Definition at line 307 of file EarthModelBinned.h.

Referenced by ClearModel(), FillPath(), and LoadModel().

◆ fnLatBins

int OscProb::EarthModelBinned::fnLatBins
protected

Definition at line 309 of file EarthModelBinned.h.

Referenced by ClearModel(), FillPath(), LoadModel(), and RecordLatLonBinCrossings().

◆ fnLonBins

int OscProb::EarthModelBinned::fnLonBins
protected

Definition at line 308 of file EarthModelBinned.h.

Referenced by ClearModel(), DetDistForNextLonBin(), and LoadModel().

◆ fNuPath

std::vector<NuPath> OscProb::EarthModelBase::fNuPath
protectedinherited

◆ fRadiusMax

double OscProb::EarthModelBase::fRadiusMax
protectedinherited

◆ fRemoveSmallPaths

bool OscProb::EarthModelBase::fRemoveSmallPaths
protectedinherited

The documentation for this class was generated from the following files: