PremModel.cxx

Go to the documentation of this file.

 
//
// Implements an earth model with spherical shells
//
// jcoelho\@apc.in2p3.fr
 
#include <cmath>
#include <fstream>
#include <iostream>
 
#include "PremModel.h"
#include "prem_default.hpp"
 
using namespace std;
 
using namespace OscProb;
 
const double PremModel::DET_TOL = 0.2; // Tolerance in km
 
//.............................................................................
PremModel::PremModel(string filename) : fDetLayer(0)
{
  SetRemoveSmallPaths(false);
  SetDetPos(6368);
  LoadModel(filename);
}
 
//.............................................................................
PremModel::~PremModel() {}
 
//.............................................................................
void PremModel::CleanIdentical()
{
  int i = 0;
 
  while (i < int(fPremLayers.size()) - 1) {
    if (fPremLayers[i] == fPremLayers[i + 1]) {
      fPremLayers.erase(fPremLayers.begin() + i);
      continue;
    }
 
    i++;
  }
}
 
//.............................................................................
void PremModel::AddDetLayer()
{
  CleanIdentical();
 
  fDetLayer = fPremLayers.size() - 1;
 
  if (fPremLayers.size() && fDetRadius > fPremLayers.back().radius) {
    cerr << "WARNING: Detector must be inside Earth model." << endl
         << "WARNING: Adjusting detector radius from " << fDetRadius
         << " km to " << fPremLayers.back().radius << " km." << endl;
 
    fDetRadius = fPremLayers.back().radius;
 
    return;
  }
 
  for (int i = 0; i < fPremLayers.size(); i++) {
    double radius = fPremLayers[i].radius;
 
    // See if we passed the detector and decide whether
    // to create a special layer for it
    if (radius > fDetRadius - DET_TOL) {
      fDetLayer = i;
 
      // If detector is not near boundary, add a special layer
      if (radius > fDetRadius + DET_TOL) {
        PremLayer det_layer = fPremLayers[i];
        det_layer.radius    = fDetRadius;
        fPremLayers.insert(fPremLayers.begin() + fDetLayer, det_layer);
      }
      else if (radius != fDetRadius) {
        // update detector radius
        cerr << "WARNING: Adjusting detector radius from " << fDetRadius
             << " km to " << radius << " km." << endl;
        fDetRadius = radius;
      }
 
      break;
    }
  }
}
 
//.............................................................................
void PremModel::SetDetPos(double rad, double lat, double lon)
{
  SetDetectorCoordinates(rad, lat, lon);
 
  if (fPremLayers.size()) AddDetLayer();
}
 
//.............................................................................
vector<PremLayer> PremModel::GetPremLayers() { return fPremLayers; }
 
//.............................................................................
void PremModel::ClearModel()
{
  fDetLayer = 0;
  fPremLayers.clear();
}
 
//.............................................................................
void PremModel::AddLayer(double radius, double density, double zoa,
                         double layer)
{
  fPremLayers.push_back(PremLayer(radius, density, zoa, layer));
}
 
//.............................................................................
void PremModel::LoadModel(string filename)
{
  // Clear the current model
  ClearModel();
 
  // Use default if no file provided
  if (filename == "") { filename = PREM_DEFAULT; }
 
  // Open the file
  ifstream fin;
  fin.open(filename.c_str());
 
  if (!fin) {
    cerr << "ERROR: File " << filename << " not found!" << endl;
    return;
  }
 
  // Variables for storing table rows
  float radius, density, zoa, layer;
 
  // Keep track of previous radius
  double rprev = 0;
 
  // Loop over table rows
  while (fin >> radius >> density >> zoa >> layer) {
    // Radii must be ordered in model file
    if (radius <= rprev) {
      cerr
          << "ERROR: Radii are not sorted in increasing order in the model file"
          << endl;
      ClearModel();
      return;
    }
 
    // Add this layer to the model
    AddLayer(radius, density, zoa, layer);
  }
 
  AddDetLayer();
 
  // Set the maximum radius in the model
  fRadiusMax = fPremLayers.back().radius;
}
 
//.............................................................................
void PremModel::SetLayerZoA(int layer, double zoa)
{
  int nlayers = fPremLayers.size();
 
  // Loop over all layers and change the ones
  // with the given index
  for (int i = 0; i < nlayers; i++) {
    if (fPremLayers[i].layer != layer) continue;
 
    fPremLayers[i].zoa = zoa;
  }
}
 
//.............................................................................
double PremModel::GetLayerZoA(int layer)
{
  int nlayers = fPremLayers.size();
 
  // This check assumes the layer types are in increasing order
  if (layer > fPremLayers.back().layer) {
    cerr << "ERROR: Not that many layer types. Returning 0" << endl;
    return 0;
  }
 
  for (int i = 0; i < nlayers; i++) {
    if (fPremLayers[i].layer != layer)
      continue;
 
    else
      return fPremLayers[i].zoa;
  }
 
  // End of vector reached without finding input type
  cerr << "ERROR: layer type not found. Returning 0" << endl;
  return 0;
}
 
//.............................................................................
void PremModel::SetTopLayerSize(double thick)
{
  if (thick <= 0) {
    cerr << "WARNING: Layer thickness should be positive. Do nothing." << endl;
    return;
  }
 
  int nlayers = fPremLayers.size();
 
  if (nlayers < 1) {
    cerr << "WARNING: PremModel has no layers. Do nothing." << endl;
    return;
  }
 
  double bottomRadius;
  if (nlayers > 1)
    bottomRadius = fPremLayers[nlayers - 2].radius;
  else
    bottomRadius = 0;
 
  double topRadius = bottomRadius + thick;
 
  fPremLayers[nlayers - 1].radius = topRadius;
 
  // Update max radius
  fRadiusMax = fPremLayers.back().radius;
}
 
//.............................................................................
void PremModel::AddPath(double length, PremLayer pl)
{
  AddPathSegment(length, pl.density, pl.zoa, pl.layer);
}
 
//.............................................................................
int PremModel::FillPath(double cosT, double phi)
{
  // Clear current path sequence
  fNuPath.clear();
 
  // Do nothing if cosine is unphysical
  if (fabs(cosT) > 1) return 0;
 
  // Define the minimum path radius
  double minR   = fDetRadius * sqrt(1 - cosT * cosT);
  double minRsq = minR * minR;
 
  // Set the top layer index
  int toplayer = fPremLayers.size() - 1;
 
  // Find the inner-most crossed layer
  int minlayer = 0;
  while (fPremLayers[minlayer].radius < minR) minlayer++;
 
  // Compute the number of path segments needed
  int nsteps = toplayer - fDetLayer;
  if (cosT < 0) nsteps += 2 * (fDetLayer - minlayer) + 1;
 
  // Start at the top layer and go down
  int layer = toplayer;
  int dl    = -1;
 
  // Loop over all path segments
  for (int i = 0; i < nsteps; i++) {
    // Get square of the path length between this layer's
    // outer radius and  inner-most radius
    double L1 = pow(fPremLayers[layer].radius, 2) - minRsq;
 
    // If L1 is negative, outer radius is not crossed.
    // This only happens if detector is at the top layer and the
    // neutrino is coming from above.
    if (L1 < 0) return true;
 
    // Get square of the path length between this layer's
    // inner radius and  inner-most radius
    double L2 = -minRsq;
    if (layer > 0) L2 += pow(fPremLayers[layer - 1].radius, 2);
 
    // If L2 is negative, inner radius is not crossed,
    // so set this as the minimum layer.
    bool ismin = (L2 <= 0 && cosT < 0);
 
    // Store the path segment length
    double dL;
 
    // If it's the minimum layer, connect two outer radius
    // crossing points. If not, compute difference between
    // inner and outer radius crossings.
    if (ismin)
      dL = 2 * sqrt(L1);
    else if (L2 >= 0)
      dL = sqrt(L1) - sqrt(L2);
    else
      dL = sqrt(L1); // This should never actually happen,
                     // but protect in case L2 is slightly
                     // negative when it should be zero.
                     // e.g. arriving at the detector from above.
 
    // Add this path segment to the sequence
    AddPath(dL, fPremLayers[layer]);
 
    // If we reached the inner-most layer,
    // start moving up again.
    if (ismin) dl = 1;
 
    // Move to next layer
    layer += dl;
  }
 
  // Return the number of path segments
  return fNuPath.size();
}