InclinedAtmosphericProfile.cc

Go to the documentation of this file.

#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>
#include <atm/Atmosphere.h>

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <utl/CoordinateSystemPtr.h>

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/PhysicalConstants.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/GeometryUtilities.h>

#include <det/Detector.h>

#include <vector>
#include <sstream>

using namespace utl;
using namespace fwk;
using namespace atm;
using namespace std;


//#define DEBUG


InclinedAtmosphericProfile::InclinedAtmosphericProfile(const Point& core,
                                                       const Vector& axis,
                                                       const double deltaX)
{
  // axis points downward for upward going showers
  bool upward = (axis.GetZ(LocalCoordinateSystem::Create(core)) < 0);
  const Vector dir = -axis;
  const Vector dirInverse = axis;

#ifdef DEBUG
  INFO(string("Running ") + (upward ? "upward" : "downward") + " mode!");
#endif

  // If event is "skimming" the earth, we need to revert through atmosphere at vertex
  bool skimming = false;

  const auto& wgs84 = utl::ReferenceEllipsoid::Get(utl::ReferenceEllipsoid::eWGS84);

  auto& detector = det::Detector::GetInstance();
  auto& atmosphere = detector.GetAtmosphere();

  const auto& atmDepthVsHeight = atmosphere.EvaluateDepthVsHeight();
  const auto& atmHeightVsDepth = atmosphere.EvaluateHeightVsDepth();
  const auto& atmDensityVsHeight = atmosphere.EvaluateDensityVsHeight();

  // Safety margin to the cosmos .... and to earth
  const double atmSafetyMargin = 1*mm;
  const double earthSafetyMargin = 1*m;

  /* Get upper and lower limit for height of the atmosphere

     ---X-------------------- Top of Atmosphere => First
         \
          \
     ------X----------------- Actual altitude of core
     -------X---------------- Surface of WGS84 => Last */

  const double atmHeightMax = atmDepthVsHeight.MaxX() - atmSafetyMargin;
  const double atmHeightMin = atmDepthVsHeight.MinX() + earthSafetyMargin;
  const double atmHeightMin_density = atmDensityVsHeight.MinX() + earthSafetyMargin;
  const double minAtm = max(atmHeightMin, atmHeightMin_density);
  const double maxVerticalDepth = atmDepthVsHeight.Y(atmHeightMin);
  const double minVerticalDepth = atmDepthVsHeight.Y(atmHeightMax);

   // ---------------------- enter/leave atmosphere -------------------------

  const vector<Point> pAtm =
    Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), atmHeightMax, Line(core, dir));

  // Catch unphysical shower geometry - less than two intersections with
  // atmosphere at this point means no shower at all
  if (pAtm.size() < 2) {
    ostringstream err;
    err << "This shower geometry does not intersect the atmosphere!!!" << endl;
    ERROR(err);
    throw InclinedAtmosphereModelException(err.str());
  }

  // ---------------------- enter/leave earth ------------------------------

  const vector<Point> pEarth =
    Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), minAtm, Line(core, dir));
  // pEarth[0,1] points, where particle [leaves, enters] the earth's ellipsoid (order defind by direction of line)

  Point firstPoint;
  Point lastPoint;
  switch (pEarth.size()) {
  case 2:
    if (upward) {
      firstPoint = pEarth[0];  // Point, where particle enters the atmosphere (leaves earth)
      lastPoint = pAtm[0];     // Point, where particle leaves the atmosphere
    } else {
      firstPoint = pAtm[1];    // Point, where particle enters the atmosphere
      lastPoint = pEarth[1];   // Point, where particle enters earth
    }
    break;
  default:
    lastPoint = pAtm[0];       // Point, where particle leaves the atmosphere
    firstPoint = pAtm[1];      // Point, where particle enters the atmosphere

    // here we need to go through the atmosphere in the standard way because shower does not intersect ground
    upward = false;

#ifdef DEBUG
    cout << "Scraper!\n";
#endif
    break;
  }

#ifdef DEBUG
  const auto& csECEF = fwk::CoordinateSystemRegistry::Get("ECEF");
  cout << "Core\t" << core.GetCoordinates(csECEF) << "\n"
          "Destination\t" << lastPoint.GetCoordinates(csECEF) << "\n"
          "First\t" << firstPoint.GetCoordinates(csECEF) << endl;
#endif

  /* Now we'll loop along the shower in steps of deltaX_incl =: deltaX
     whilst deltaX_vert = deltaX_incl*cos(theta)
     and cos(theta) is the angle between z-axis in the P_i-frame.
     From these we'll get a vertical height which gives P_i+1 as
     intersection from shower axis with the atmospheric ellipsoid
     Giving us the distance between P_core and P_i+1.*/

  Point iPoint = firstPoint;
  Point nextPoint;
  double cosTheta = 0;
  try {
    cosTheta = dirInverse.GetCosTheta(LocalCoordinateSystem::Create(UTMPoint(firstPoint, wgs84)));
  } catch (UTMPoint::UTMException& e) {
    ostringstream err;
    err << "InclinedAtmosphericProfile: un-physical local point ";
    ERROR(err);
    throw InclinedAtmosphereModelException(err.str());
  }

  // Temporary variables for calculation
  double verticalDeltaX = 0;
  double tmpHeight = UTMPoint(firstPoint, wgs84).GetHeight();
  double tmpDepth = atmDepthVsHeight.Y(tmpHeight);
  vector<Point> tmpPoints;

  // Set vertical start points
  double tmpSlantDepth = upward ? cosTheta*(tmpDepth - atmDepthVsHeight.Y(0)) : tmpDepth;

  // Result vectors
  vector<double> verticalHeight;
  vector<double> slantDepth;
  vector<double> distanceToImpactPoint;

  // Store initial values as first entry
  verticalHeight.push_back(atmHeightVsDepth.Y(tmpDepth));
  slantDepth.push_back(log(tmpSlantDepth)); // log(slantDepth)!
  distanceToImpactPoint.push_back(0 - (firstPoint - core).GetMag());
  // Pi0 at first

  // Start loop
  while (true) {

    // Double-check whether everything made sense up to this point
    try {
      cosTheta = dirInverse.GetCosTheta(LocalCoordinateSystem::Create(UTMPoint(iPoint, wgs84)));
    } catch (UTMPoint::UTMException& e) {
      ostringstream err;
      err << "InclinedAtmosphericProfile: un-physical local point ";
      ERROR(err.str());
      throw InclinedAtmosphereModelException(err.str());
    }

    verticalDeltaX = deltaX * cosTheta;
    tmpDepth += verticalDeltaX;

    // Break loop as soon as we reach the last point
    if (maxVerticalDepth <= tmpDepth || tmpDepth <= minVerticalDepth) {
      tmpHeight = UTMPoint(lastPoint, wgs84).GetHeight();
      verticalDeltaX = atmDepthVsHeight.Y(tmpHeight) - (tmpDepth - verticalDeltaX);
      tmpSlantDepth += verticalDeltaX / cosTheta;

      // Save vertical height of new point
      verticalHeight.push_back(tmpHeight);
      // Save slant depth at new point
      slantDepth.push_back(log(tmpSlantDepth)); // log(slantDepth)!
      // Save distance to impact point
      distanceToImpactPoint.push_back((lastPoint - core).GetMag());

      break;
    } else {
      tmpSlantDepth += deltaX;
      tmpHeight = atmHeightVsDepth.Y(tmpDepth);

      // Calculate new point
      tmpPoints = Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), tmpHeight, Line(core, dir));

      if (tmpPoints.size() < 2) {
        if (!skimming) {
          // For showers with no earth intersection
          tmpDepth -= verticalDeltaX;
          tmpSlantDepth -= deltaX;
          tmpHeight = atmHeightVsDepth.Y(tmpDepth);
          tmpPoints = Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), tmpHeight, Line(core, dir));
          tmpSlantDepth += atmDensityVsHeight.Y(tmpHeight) * (tmpPoints[0] - tmpPoints[1]).GetMag();
          skimming = true;
        } else {
          // Beak, if we have reached space in this step! Should not happen...
          INFO("Reached space!");
          break;
        }
      }

      // Get correct next point
      nextPoint = skimming ? (upward ? tmpPoints[1] : tmpPoints[0]) : (upward ? tmpPoints[0] : tmpPoints[1]);

      // Save vertical height of new point
      verticalHeight.push_back(tmpHeight);
      // Save slant depth at new point
      slantDepth.push_back(log(tmpSlantDepth)); // log(slantDepth)!
      // Save distance to impact point
      distanceToImpactPoint.push_back((nextPoint - firstPoint).GetMag() - (core - firstPoint).GetMag());

      // Next step
      iPoint = nextPoint;
    }
  } // End of loop

  // Check for empty vectors
  if (slantDepth.size() <= 1) {
    ostringstream err;
    err << "No points calculated. Table is empty. " << endl;
    ERROR(err);
    throw InclinedAtmosphereModelException(err.str());
  }

#ifdef DEBUG  
  const double testCosTheta = dirirection.GetCosTheta(LocalCoordinateSystem::Create(UTMPoint(core, wgs84)));
  cout << "CosTheta\tSlantDepth\tverticalHeight\tDistance\n";
  for (unsigned int i = 0; i < verticalHeight.size(); ++i) {
    cout << -testCosTheta << '\t'
         << exp(slantDepth[i])/g*cm2 << '\t'
         << verticalHeight[i]/m << '\t'
         << distanceToImpactPoint[i]/m << endl;
  }
#endif

  // Currently no uncertainties ...
  vector<double> verticalHeightErrors(slantDepth.size());
  vector<double> slantDepthErrors(slantDepth.size());
  vector<double> distanceErrors(slantDepth.size());

  // reverse -----------------------------
  vector<double> slantDepthReverse(slantDepth.size());
  vector<double> verticalHeightReverse(slantDepth.size());

  reverse_copy(slantDepth.begin(), slantDepth.end(), slantDepthReverse.begin());
  reverse_copy(verticalHeight.begin(), verticalHeight.end(), verticalHeightReverse.begin());

  const TabulatedFunctionErrors tabSlantDist(distanceToImpactPoint, distanceErrors,
                                             slantDepth, slantDepthErrors);

  fSlantDepthVsDistance =
    new ProfileResult(tabSlantDist, ProfileResult::eLinear, ProfileResult::eLog);

  const TabulatedFunctionErrors tabDistSlant(slantDepth, slantDepthErrors,
                                             distanceToImpactPoint, distanceErrors);

  fDistanceVsSlantDepth =
    new ProfileResult(tabDistSlant, ProfileResult::eLog, ProfileResult::eLinear);

  const TabulatedFunctionErrors tabHeightSlant(slantDepthReverse, slantDepthErrors,
                                               verticalHeightReverse, verticalHeightErrors);

  fHeightVsSlantDepth =
    new ProfileResult(tabHeightSlant, ProfileResult::eLog, ProfileResult::eLinear);

  const TabulatedFunctionErrors tabHeightDistance(distanceToImpactPoint, distanceErrors,
                                                  verticalHeight, verticalHeightErrors);

  fHeightVsDistance =
    new ProfileResult(tabHeightDistance, ProfileResult::eLinear, ProfileResult::eLinear);
}


InclinedAtmosphericProfile::~InclinedAtmosphericProfile()
{
  delete fSlantDepthVsDistance;
  delete fDistanceVsSlantDepth;
  delete fHeightVsSlantDepth;
  delete fHeightVsDistance;
}


const atm::ProfileResult&
InclinedAtmosphericProfile::EvaluateSlantDepthVsDistance()
  const
{
  if (fSlantDepthVsDistance)
    return *fSlantDepthVsDistance;
  else
    throw InclinedAtmosphereModelException("First Initialize the InclinedAtmosphericProfile Model !!");
}


const atm::ProfileResult&
InclinedAtmosphericProfile::EvaluateDistanceVsSlantDepth()
  const
{
  if (fDistanceVsSlantDepth)
    return *fDistanceVsSlantDepth;
  else
    throw InclinedAtmosphereModelException("First Initialize the InclinedAtmosphericProfile Model !!");
}


const atm::ProfileResult&
InclinedAtmosphericProfile::EvaluateHeightVsSlantDepth()
  const
{
  if (fHeightVsSlantDepth)
    return *fHeightVsSlantDepth;
  else
    throw InclinedAtmosphereModelException("First Initialize the InclinedAtmosphericProfile Model !!");
}


const atm::ProfileResult&
InclinedAtmosphericProfile::EvaluateHeightVsDistance()
  const
{
  if (fHeightVsDistance)
    return *fHeightVsDistance;
  else
    throw InclinedAtmosphereModelException("First Initialize the InclinedAtmosphericProfile Model !!");
}


//############# Member functions ######################

double
InclinedAtmosphericProfile::IntegratedGrammage(const utl::Point& pStart,
                                               const utl::Point& pStop,
                                               const double delta)
{
  const auto& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  auto& detector = det::Detector::GetInstance();
  auto& atmosphere = detector.GetAtmosphere();
  const ProfileResult& atmDepthVsHeight = atmosphere.EvaluateDepthVsHeight();

  const double atmHeightMin = atmDepthVsHeight.MinX();
  const double atmHeightMax = atmDepthVsHeight.MaxX();

  Vector dir(pStop - pStart);
  const double length = dir.GetMag();
  dir.Normalize();

  double depth = 0;
  double distance = 0;

  while (distance < length) {

    const double delta_now = (length - distance < delta ? length - distance : delta);
    const Point p1 = pStart + distance * dir;
    distance += delta_now;
    const Point p2 = pStart + distance * dir;

    double cosLocalTheta = 0;
    try {
      const UTMPoint utm(p1, wgs84);
      cosLocalTheta = dir.GetCosTheta(LocalCoordinateSystem::Create(utm));
    } catch (UTMPoint::UTMException& e) {
      continue;
    }

    const double heightP1 = wgs84.PointToLatitudeLongitudeHeight(p1).get<2>();
    const double heightP2 = wgs84.PointToLatitudeLongitudeHeight(p2).get<2>();

    if (heightP1 <= atmHeightMin || heightP1 >= atmHeightMax ||
        heightP2 <= atmHeightMin || heightP2 >= atmHeightMax) {
      // add nothing
      continue;
    }

    const double depthP1 = atmDepthVsHeight.Y(heightP1);
    const double depthP2 = atmDepthVsHeight.Y(heightP2);

    depth += (depthP1 - depthP2) / cosLocalTheta;

  }

  return depth;
}