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 217 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 96 of file EarthModelBinned.cxx.

97{
98 SetDetPos(6368, 0, 0);
99 LoadModel(filename);
100 SetRemoveSmallPaths(false);
101 fLonError = 0;
102}
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 108 of file EarthModelBinned.cxx.

108{}

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 162 of file EarthModelBinned.cxx.

165{
166 fEarthBins.push_back(EarthBin(radius_out, radius_in, latitude, longitude,
167 density, zoa, index));
168}
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 408 of file EarthModelBinned.cxx.

409{
410 AddPathSegment(length, bin.density, bin.zoa, bin.index);
411}
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 122 of file EarthModelBinned.cxx.

123{
124 fnDepthBins = 0;
125 fnLonBins = 0;
126 fnLatBins = 0;
127 fEarthBins.clear();
128}
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 431 of file EarthModelBinned.cxx.

432{
433 // Find next lat bin boundary (bin center +/- 1/2-binwidth)
434 double NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
435 // Turn around if lat reaches +/- 90 degrees (or close enough)
436 if (abs(NewLat) - 0.5 * M_PI > -1.0e-14) {
437 if (L.maxreached) return -1; // Only turn around once
438 L.maxreached = true;
439 // cout << "Max latitude reached..." << endl;
440 L.sign = -L.sign;
441 NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
442 if (abs(NewLat) - 0.5 * M_PI > -1.0e-14)
443 return -1; // This would happen if there's only one latitude bin...
444 }
445
446 // Calculate parts of answer
447 double sinNewLat = sin(NewLat);
448 double NsinSqNewLat = -sinNewLat * sinNewLat;
449 double radicand = fC.maxSinSqLat +
450 NsinSqNewLat; // 1-[sin(NewLat)]^2-[sin(phi)cos(DetLat)]^2
451 // Turn around if the radicand is 0 (close enough) or less
452 if (radicand < 1.0e-14) {
453 if (L.maxreached) return -1; // Only turn around once
454 L.maxreached = true;
455 // cout << "Max latitude reached..." << endl;
456 L.sign = -L.sign;
457 NewLat = fEarthBins[cur_index].latitude + L.sign * L.dLat;
458 sinNewLat = sin(NewLat);
459 NsinSqNewLat = -sinNewLat * sinNewLat;
460 radicand = fC.maxSinSqLat + NsinSqNewLat;
461 if (radicand < 1.0e-14)
462 return -1; // The rest of the track is contained in this lat bin
463 }
464 double denominator = fC.gammaSq + NsinSqNewLat;
465
466 // Calculate distance from the detector at next lat bin boundary
467 double answer;
468 if (denominator == 0) {
469 answer = 0.5 * (fC.xlatextreme - fC.rDetSinDetLat / fC.gamma);
470 // There is a "possible" case where answer would be 0:
471 // if the detector is at the center of the Earth... just don't...
472 }
473 else {
474 answer = -(fC.rDetCosT * NsinSqNewLat + fC.rDetGammaSinDetLat +
475 L.sign * fC.rDetSinT * sinNewLat * sqrt(radicand)) /
476 denominator;
477 }
478
479 // Check to see if going to the edge of the bin jumped to an invalid
480 // part of the function (cosA*rDetCosDetLat/beta = x at max lat)
481 if (L.maxreached && answer >= fC.xlatextreme) {
482 return -1; // The rest of the track is contained in this lat bin
483 }
484
485 return answer;
486}
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 501 of file EarthModelBinned.cxx.

502{
503 double twoPi = 2 * M_PI;
504 // Increment to next bin in direction dLon
505 double lon = prev_lon + L.dLon;
506 if (lon < 0 || lon >= twoPi) { lon -= floor(lon / twoPi) * twoPi; }
507 if (L.dLon > 0) {
508 L.nextBin = L.bin + 1;
509 if (L.nextBin >= fnLonBins) {
510 // cross from 2pi to 0
511 L.nextBin = 0;
512 }
513 }
514 else {
515 L.nextBin = L.bin - 1;
516 if (L.nextBin < 0) {
517 // cross from 0 to 2pi
518 L.nextBin = fnLonBins - 1;
519 }
520 }
521
522 // Check if lon is outside of range covered by neutrino's trajectory
523 if (L.max < L.min) {
524 if (0 > lon - L.min && L.max - lon < 0) return -1;
525 }
526 else {
527 if (0 > lon - L.min || L.max - lon < 0) return -1;
528 }
529
530 // Parts of the answer
531 double cosNewLon = cos(lon);
532 double sinNewLon = sin(lon);
533 double expression = fC.sinDetLon * cosNewLon - sinNewLon * fC.cosDetLon;
534 double denominator = fC.alpha * expression +
535 fC.sinTsinAcosDetLon * cosNewLon +
536 fC.sinTsinAsinDetLon * sinNewLon;
537
538 // Calculate & return detector distance
539 if (denominator == 0) { // Is there ever a case when this is true???
540 // Error message will be printed at the end of the FillPath function
541 L.err_message.assign("Longitude along neutrino trajectory = " +
542 std::to_string(lon));
543 L.error = -1;
544 return -1;
545 }
546 return fC.rDetCosDetLat * expression / denominator;
547}
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 666 of file EarthModelBinned.cxx.

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

116{ 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.

◆ 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 347 of file EarthModelBinned.cxx.

348{
349 if (index > fmaxRegIndex) {
350 cerr << "ERROR: Index " << index << " does not exist in model "
351 << "(max index = " << fmaxRegIndex << "). Returning 0" << endl;
352 return 0;
353 }
354
355 int nbins = fEarthBins.size();
356
357 for (int i = 0; i < nbins; i++) {
358 if (fEarthBins[i].index != index)
359 continue;
360
361 else
362 return fEarthBins[i].zoa;
363 }
364
365 // End of vector reached without finding index
366 cerr << "ERROR: Index " << index << " not found! Returning 0" << endl;
367 return 0;
368}
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 189 of file EarthModelBinned.cxx.

190{
191 // Clear the current model
192 ClearModel();
193 double EarthRadius = 6371.0;
194
195 // Use default if no file provided
196 if (filename == "") { filename = PREM3D_DEFAULT; }
197 cout << "Loading Earth model table from " << filename << "..." << endl;
198
199 // Open the file
200 ifstream fin;
201 fin.open(filename.c_str());
202 if (!fin) {
203 cerr << "ERROR: File " << filename << " not found!" << endl;
204 return;
205 }
206
207 // Variables for storing table rows
208 float depth, latitude, longitude, density, zoa, index;
209
210 // Keep track of previous depth/latitude/longitude
211 double radius_prev = 0;
212 double lat_prev = -90.0;
213 double lon_prev = 0;
214
215 // Hold onto previous different depth, so we have depth maximum for bin
216 double radius_min = radius_prev;
217
218 // Initialize max region index
219 fmaxRegIndex = -1;
220
221 // Loop over table rows
222 while (fin >> longitude >> latitude >> depth >> density >> zoa >> index) {
223 double outer_radius = EarthRadius - depth;
224
225 // Move minimum depth (and reset latitude/longitude) if depth has changed
226 // from previous one
227 if (outer_radius > radius_prev) {
228 radius_min = radius_prev;
229 lat_prev = -90.0;
230 lon_prev = 0.0;
231 fnDepthBins++;
232 }
233 else if (outer_radius <
234 radius_prev) { // Depths must be ordered in model file
235 cerr << "ERROR: Depths are not sorted in decreasing order in the model "
236 "file "
237 << "(or depth is greater than " << EarthRadius << " km)." << endl;
238 ClearModel();
239 return;
240 }
241 else if (longitude > lon_prev) { // Reset latitude if longitude has changed
242 // from previous one
243 lat_prev = -90.0;
244 fnLonBins++;
245 }
246 else if (longitude < lon_prev) { // Longitudes must be ordered in model file
247 // (after depths)
248 cerr << "ERROR: Longitudes are not sorted in increasing order (after "
249 "depths) in the model file."
250 << endl;
251 ClearModel();
252 return;
253 }
254 else if (latitude < lat_prev) { // Latitudes must be ordered in model file
255 // (after longitudes)
256 cerr << "ERROR: Latitudes are not sorted in increasing order (after "
257 "longitudes) in the model file."
258 << endl;
259 ClearModel();
260 return;
261 }
262
263 // Add this bin to the model
264 AddBin(outer_radius, radius_min, latitude, longitude, density, zoa, index);
265
266 // Check region index
267 if (index > fmaxRegIndex) { fmaxRegIndex = index; }
268
269 // Set previous coordinates for next bin
270 radius_prev = outer_radius;
271 lat_prev = latitude;
272 lon_prev = longitude;
273 }
274
275 // Set Outermost Radius
276 fRadiusMax = radius_prev;
277
278 // Modify number of longitude bins to how many there actually are (instead of
279 // how many times we changed the longitude bin without changing the depth bin
280 // = nLonBins-1 times per depth bin)
282 fInvLonBinWidth = 0.5 * fnLonBins / M_PI; // in radians^(-1)
283 fHalfLonBinWidth = M_PI / fnLonBins; // in radians
284
285 // Calculate number of latitude bins
287 fInvLatBinWidth = fnLatBins / M_PI; // in radians^(-1)
288 fHalfLatBinWidth = 0.5 * M_PI / fnLatBins; // in radians
289
290 // Error if first latitude bin center deviates from about half of the
291 // calculated latitude bin width
292 if (abs(fHalfLatBinWidth - (fEarthBins[0].latitude + 0.5 * M_PI)) >
293 fHalfLatBinWidth * 0.01) { // Warn if first latitude bin center isn't
294 // about half of the latitude bin width
295 cerr << "ERROR: Latitude of 1st bin center (" << fEarthBins[0].latitude
296 << ") is not in the middle of the first bin, whose size should be "
297 << 2 * fHalfLatBinWidth << " rad, as calculated from having "
298 << fnLatBins << " latitude bins." << endl;
299 ClearModel();
300 return;
301 }
302
303 cout << "\t...done (" << fEarthBins.size() << " bins: " << fnDepthBins
304 << " depth, " << fnLatBins << " latitude, and " << fnLonBins
305 << " longitude)." << endl;
306}
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 630 of file EarthModelBinned.cxx.

631{
632 // Find difference between longitude and beginning of first lon bin
633 double londiff = longitude - (fEarthBins[0].longitude - fHalfLonBinWidth);
634
635 // Shift londiff to be between 0 and 2pi
636 if (londiff < 0 || londiff >= 2 * M_PI)
637 londiff -= floor(londiff / (2 * M_PI));
638
639 // Number of binwidths after from first lon bin
640 return floor(londiff * fInvLonBinWidth);
641}

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 561 of file EarthModelBinned.cxx.

565{
566 // Keep going as long as there is a lat or lon bin crossing before
567 // the next depth bin crossing
568 while (detDist_nextDbin < latI.detDist_nextBin ||
569 detDist_nextDbin < lonI.detDist_nextBin) {
570 if (latI.detDist_nextBin >
571 lonI.detDist_nextBin) { // if next lat bin crossing is before next lon
572 // bin crossing
573 // Add segment to next lat bin to path
574 AddPath(DetDist - latI.detDist_nextBin, fEarthBins[index]);
575 double latitude_crossing =
576 fEarthBins[index].latitude + latI.sign * latI.dLat;
577
578 // Move lat bins
579 DetDist = latI.detDist_nextBin;
580 index += latI.nextBin - latI.bin;
581 latI.bin = latI.nextBin;
582 // cout << "Latitude bin crossing at x=" << DetDist << "km,lat=" <<
583 // latitude_crossing << "rad...now in Earth bin " << index << " = (r"
584 // << floor(index/(fnLonBins*fnLatBins)) << ", lon" << lonI.bin << ",
585 // lat" << latI.bin << ")" << endl;
586
587 // Find next lat bin
588 if (latI.dLat ==
589 0) { // Only 1 lat bin crossing along entire neutrino trajectory
590 latI.detDist_nextBin = -1;
591 }
592 else {
593 latI.detDist_nextBin = DetDistForNextLatBin(index, latI);
594 latI.nextBin = latI.bin + latI.sign;
595 }
596 }
597 else { // if next lon bin crossing is before (or simultaneous with) next lat
598 // bin corssing
599 // Add segment to next lon bin to path
600 AddPath(DetDist - lonI.detDist_nextBin, fEarthBins[index]);
601
602 // Move lon bins
603 DetDist = lonI.detDist_nextBin;
604 index += (lonI.nextBin - lonI.bin) * fnLatBins;
605 lonI.bin = lonI.nextBin;
606 // cout << "Longitude bin crossing at x=" << DetDist << "km...now in
607 // Earth bin " << index << " = (r" <<
608 // floor(index/(fnLonBins*fnLatBins)) << ", lon" << lonI.bin << ",
609 // lat" << latI.bin << ")" << endl;
610
611 // Find next lon bin
612 if (lonI.dLon ==
613 0) { // Only 1 lon bin crossing along entire neutrino trajectory
614 lonI.detDist_nextBin = -1;
615 }
616 else {
617 lonI.detDist_nextBin =
618 DetDistForNextLonBin(fEarthBins[index].longitude, lonI);
619 }
620 }
621 }
622}

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 378 of file EarthModelBinned.cxx.

379{
380 if (index > fmaxRegIndex) {
381 cerr << "WARNING: Index " << index << " does not exist in model "
382 << "(max index = " << fmaxRegIndex << ")." << endl
383 << "Doing nothing." << endl;
384 return;
385 }
386
387 int nbins = fEarthBins.size();
388
389 // Loop over all bins and change the ones
390 // with the given index
391 for (int i = 0; i < nbins; i++) {
392 if (fEarthBins[i].index != index) continue;
393
394 double density_init = fEarthBins[i].density;
395 fEarthBins[i].density = factor * density_init;
396 }
397}

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 144 of file EarthModelBinned.cxx.

145{
146 SetDetectorCoordinates(rad, lat, lon);
148}
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 318 of file EarthModelBinned.cxx.

319{
320 if (index > fmaxRegIndex) {
321 cerr << "WARNING: Index " << index << " does not exist in model "
322 << "(max index = " << fmaxRegIndex << ")." << endl
323 << "Doing nothing." << endl;
324 return;
325 }
326
327 int nbins = fEarthBins.size();
328
329 // Loop over all bins and change the ones
330 // with the given index
331 for (int i = 0; i < nbins; i++) {
332 if (fEarthBins[i].index != index) continue;
333
334 fEarthBins[i].zoa = zoa;
335 }
336}

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 334 of file EarthModelBinned.h.

◆ fHalfLatBinWidth

double OscProb::EarthModelBinned::fHalfLatBinWidth
protected

Definition at line 329 of file EarthModelBinned.h.

Referenced by FillPath(), and LoadModel().

◆ fHalfLonBinWidth

double OscProb::EarthModelBinned::fHalfLonBinWidth
protected

Definition at line 328 of file EarthModelBinned.h.

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

◆ fInvLatBinWidth

double OscProb::EarthModelBinned::fInvLatBinWidth
protected

Definition at line 327 of file EarthModelBinned.h.

Referenced by FillPath(), and LoadModel().

◆ fInvLonBinWidth

double OscProb::EarthModelBinned::fInvLonBinWidth
protected

Definition at line 326 of file EarthModelBinned.h.

Referenced by LoadModel(), and LonBinIndex().

◆ fLonError

int OscProb::EarthModelBinned::fLonError
protected

Definition at line 336 of file EarthModelBinned.h.

Referenced by EarthModelBinned().

◆ fmaxRegIndex

int OscProb::EarthModelBinned::fmaxRegIndex
protected

Definition at line 331 of file EarthModelBinned.h.

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

◆ fnDepthBins

int OscProb::EarthModelBinned::fnDepthBins
protected

Definition at line 323 of file EarthModelBinned.h.

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

◆ fnLatBins

int OscProb::EarthModelBinned::fnLatBins
protected

Definition at line 325 of file EarthModelBinned.h.

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

◆ fnLonBins

int OscProb::EarthModelBinned::fnLonBins
protected

Definition at line 324 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: