EarthModelBase.cxx

Go to the documentation of this file.

/*************************************************************************
 * EarthModelBase.cxx
 * Base class for all Earth model clases in OscProb
 * (created by Rebekah Pestes based on PremModel.cxx)
 *************************************************************************/
 
#include <cmath>
#include <iostream>
 
#include "EarthModelBase.h"
 
using namespace std;
 
using namespace OscProb;
 
//.............................................................................
void EarthModelBase::SetDetectorCoordinates(double rad, double lat, double lon)
{
  // Force radius to be non-negative
  if (rad < 0) {
    cerr << "WARNING: Negative radius detected. Setting to absolute value."
         << endl;
    rad = -rad;
  }
  fDetRadius = rad;
 
  // Force latitude to be between -90 and 90 deg
  lat -= floor((lat + 90.0) / 360.0) * 360.0;
  if (lat > 90) { lat = 180.0 - lat; }
  fDetLat = lat / 180.0 * M_PI; // convert to radians
 
  // Force longitude to be between 0 and 360 deg
  lon -= floor(lon / 360.0) * 360.0;
  fDetLon = lon / 180.0 * M_PI; // convert to radians
}
 
//.............................................................................
vector<NuPath> EarthModelBase::GetNuPath() { return fNuPath; }
 
//.............................................................................
void EarthModelBase::AddPathSegment(double length, double density, double zoa,
                                    int index)
{
  fNuPath.push_back(NuPath(length, density, zoa, index));
}
 
//.............................................................................
double EarthModelBase::GetTotalL(double cosT)
{
  if (fabs(cosT) > 1) return 0;
 
  double sinsqrT = 1 - cosT * cosT;
 
  return -fDetRadius * cosT +
         sqrt(fRadiusMax * fRadiusMax - fDetRadius * fDetRadius * sinsqrT);
}
 
//.............................................................................
double EarthModelBase::GetCosT(double L)
{
  if (L < fRadiusMax - fDetRadius) return 1;
  if (L > fRadiusMax + fDetRadius) return -1;
 
  return (fRadiusMax * fRadiusMax - fDetRadius * fDetRadius - L * L) /
         (2 * fDetRadius * L);
}
 
//.............................................................................
vector<NuPath> EarthModelBase::GetMergedPaths(double prec)
{
  // The output vector
  vector<NuPath> mergedPath;
 
  // Start with the first path
  NuPath path = fNuPath[0];
 
  // Track the total length
  double totL = 0;
 
  // Loop over all paths starting from second
  for (int i = 1; i < fNuPath.size(); i++) {
    // If this path electron density is beyond the tolerance
    if (fabs(path.density * path.zoa - fNuPath[i].density * fNuPath[i].zoa) >
        prec * path.zoa) {
      // Add merged path to vector
      mergedPath.push_back(path);
 
      // Set this path as new merged path
      path = fNuPath[i];
    }
    else { // If path is within tolerance
 
      // Merge the path with current merged path
      path = AvgPath(path, fNuPath[i]);
    }
 
    // Increment total length
    totL += fNuPath[i].length;
 
  } // End of loop over paths
 
  // Add the final merged path to vector
  mergedPath.push_back(path);
 
  // If tag is true, remove small paths
  if (fRemoveSmallPaths) {
    // Start at first path
    int k = 0;
 
    // While not at the end of vector
    while (k + 1 < mergedPath.size()) {
      // If length is less than 1% of total
      if (mergedPath[k].length < 0.01 * totL) {
        // Merge path with following path
        mergedPath = MergePaths(mergedPath, k, k + 1);
      }
      // If path is long enough skip it
      else
        k++;
 
    } // End of while loop
 
  } // End of if statement
 
  // return the merged vector
  return mergedPath;
}
 
//.............................................................................
void EarthModelBase::SetRemoveSmallPaths(bool rp) { fRemoveSmallPaths = rp; }