AtmosphereUtilities.cc

Go to the documentation of this file.

#include <src/UniversalityConfig.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <det/Detector.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/AugerUnits.h>
#include <utl/UTCDateTime.h>
#include <utl/Branch.h>
#include <fwk/CentralConfig.h>
#include <src/PowerConfig.h>
#include <src/Functions.h>

using namespace std;
using namespace utl;


namespace un2 {

  /*void
  Print(const atm::ProfileResult& p,
        const string& label,
        const double xUnit, const double yUnit,
        const string& xLabel, const string& yLabel)
  {
    const auto& x = p.GetXValues();
    const auto& y = p.GetYValues();
    cout << label << '\n'
         << xLabel << '\t' << yLabel << '\n';
    for (unsigned int i = 0, n = x.size(); i < n; ++i)
      cout << x[i]/xUnit << '\t' << y[i]/yUnit << '\n';
  }

  // testing
  atm.InitSlantProfileModel(core, axis, 50*gpcm2);
  Print(atm.EvaluateSlantDepthVsDistance(), "\nEvaluateSlantDepthVsDistance:", meter, gpcm2, "Distance/m", "SlantDepth/(g/cm2)");
  Print(atm.EvaluateDistanceVsSlantDepth(), "\nEvaluateDistanceVsSlantDepth:", gpcm2, meter, "SlantDepth/(g/cm2)", "Distance/m");
  Print(atm.EvaluateHeightVsSlantDepth(), "\nEvaluateHeightVsSlantDepth:", gpcm2, meter, "SlantDepth/(g/cm2)", "Height/m");
  Print(atm.EvaluateHeightVsDistance(), "\nEvaluateHeightVsDistance:", meter, meter, "Distance/m", "Height/m");
  return 0;*/


  inline
  double
  Interpolate(const double dx, const double dy,
              const double x)
  {
    return dy * x / dx;
  }


  inline
  double
  Interpolate(const double x1, const double y1,
              const double x2, const double y2,
              const double x)
  {
    return y1 + Interpolate(x2 - x1, y2 - y1, x - x1);
  }


  inline
  double
  ExpInterpolate(const double x1, const double y1,
                 const double x2, const double y2,
                 const double x)
  {
    // assume log(y) is linear, ie y is exponential
    return exp(Interpolate(x1, log(y1), x2, log(y2), x));
  }


  /* Calculate the shower age
   *
   * differences to old method:
   * - non-abelean behaviour of calculations was taken into account
   * - using observation height gauged to what was actually used in simulations
   * - real atmospheric data can be used
   * - debug output
   * - units are now meters and radians */

  inline
  double
  CalcDXRealAtm(const bool useDL,
                const double spRadius, const double spPsi, const double groundHeight,
                const double xMax, const double theta, const int month)
  {
    // use the atmosphere configured in the bootstrap (preferably from FdRealDetConfig) to
    // do the full DX calculation. this is also used to create the atmosphere_parameters.xml
    // for faster analtic calculation.
    if (month<0)
      return 0;

    auto& detector = det::Detector::GetInstance();

    // make sure the atmosphere switches to the right month
    //detector.Update(UTCDateTime(2018, month, 15).GetTimeStamp());

    const auto& cs = detector.GetReferenceCoordinateSystem();  // ePampaAmarilla

    // referenceHeight is new effective "zero" height
    constexpr double paReferenceHeight = 1400;  // the UTM height of ePampaAmarilla origin
    const double paHeight = groundHeight - paReferenceHeight;  // universality height

    const Point core(0, 0, paHeight, cs);
    const Vector axis(1, theta, 0, cs, Vector::kSpherical);  // tilted towards x-axis

    const double gpcm2 = g/cm2;

    const auto& atm = detector.GetAtmosphere();

    atm.InitSlantProfileModel(core, axis, 0.5*gpcm2);

    // note: positive distance is downwards in direction of -axis
    const double projectedStationDistance = -spRadius * tan(theta) * cos(spPsi);

    if (useDL) {
      const auto& dVsX = atm.EvaluateDistanceVsSlantDepth();
      const double distanceToXMax = dVsX.Y(xMax * gpcm2) / meter;
      return projectedStationDistance - distanceToXMax;
    } else {
      const auto& xVsD = atm.EvaluateSlantDepthVsDistance();
      const double d = projectedStationDistance * meter;
      if (d < xVsD.MinX())
        // assume zero X
        return -xMax;
      else if (d < xVsD.MaxX()) {
        const double depthOfProjectedStation = xVsD.Y(d) / gpcm2;
        return depthOfProjectedStation - xMax;
      } else {
        // extrapolate below sea level
        const auto& dists = xVsD.GetXValues();
        const auto& depths = xVsD.GetYValues();
        const unsigned int n = dists.size();
        const double d1 = dists[n-2];
        const double d2 = dists[n-1];
        const double depth1 = depths[n-2];
        const double depth2 = depths[n-1];
        const double extrapolatedDepth = ExpInterpolate(d1, depth1, d2, depth2, d);
        return extrapolatedDepth / gpcm2 - xMax;
      }
    }
  }


  inline
  double
  CalcDXCorsikaAtm(const bool useDL,
                   const double spRadius, const double spPsi, const double groundHeight,
                   const double xMax, const double theta, const int month)
  {
    // use the maximum vertical depths from corsika profiles with the parametrised atmospheric
    // profiles from MonthlyAvgDB. results match the corsika atmosphere with ~1% precision and
    // better

    const UniversalityConfig& config = PowerConfig::GetInstance();
    fwk::CentralConfig& cc = *fwk::CentralConfig::GetInstance(/*config.GetBootstrap()*/);

    const auto& monthsDict = config.GetMonthsMapReverse();
    const auto m = monthsDict.at(month);
    auto mB = cc.GetTopBranch("AtmosphereParameters").GetChild("un2").GetChild(m);
    const auto hs = mB.GetChild("hs").Get<double>();
    const auto hs2 = mB.GetChild("hs2").Get<double>();
    // these are taken directly from the CORSIKA documentation
    const double xVert = config.GetGroundLevelDepths()[month-1];

    if (useDL) {

      const double hXmax = AnalyticHeightFromDepth(xMax, theta, xVert, hs, hs2);
      const double hProj = spRadius * cos(spPsi) * sin(theta);
      return (hXmax - hProj - 1400.) / cos(theta);

    } else {

      return AnalyticDX(spRadius, spPsi, theta, xMax, xVert, hs, hs2, groundHeight);

    }
  }


  double
  CalcDX(const bool realAtmosphere, const bool useDL,
         const double spRadius, const double spPsi, const double groundHeight,
         const double xMax, const double theta, const int month, const bool silent)
  {
    if (silent) {
      utl::ErrorLogger::GetInstance().SetStream(cerr); // redirect offline logging
      cerr.rdbuf(nullptr); // turn off offline output, which was directed to cerr
    }

    return (realAtmosphere) ?
      CalcDXRealAtm(useDL, spRadius, spPsi, groundHeight, xMax, theta, month) :
      CalcDXCorsikaAtm(useDL, spRadius, spPsi, groundHeight, xMax, theta, month);
  }

}