SimShowerProfileModel.cc

Go to the documentation of this file.

#include "SimShowerProfileModel.h"

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/TabulatedFunction.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/UTMPoint.h>

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/LocalCoordinateSystem.h>

#include <atm/ProfileResult.h>

#include <evt/Event.h>
#include <evt/Header.h>
#include <evt/ShowerSimData.h>
#include <evt/RadioSimulation.h>
#include <evt/AtmosphereParameters.h>

#include <iostream>
#include <iomanip>
#include <limits>
#include <string>
#include <sstream>
#include <vector>

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


void
SimShowerProfileModel::Init()
{
  INFO("Using SimShowerProfileModel");
}


void
SimShowerProfileModel::LazyInit()
  const
{
  const Event& event = RunController::GetInstance().GetCurrentEvent();

  if (fInitialized && fCurrentEvent == event.GetHeader().GetId())
    return;

  if (!(event.HasSimShower() && event.GetSimShower().HasAtmosphereParameters()))
    throw NoDataForModelException("SimShowerProfileModel requires that the event has a simulated shower "
                                  "with AtmosphereParameters. One of the two was not found.");

  const ShowerSimData& shower = event.GetSimShower();

  if (!shower.HasGeometry())
    throw MissingEventDataException("The shower geometry is not defined yet.");

  fCurrentEvent = event.GetHeader().GetId();
  fInitialized = true;

  // this is for paranoid people (like me): The shower coordinate
  // system might now point "upward" like an AugerCorrdinateSystem
  // in all possible cases. Thus, use the original reference frame of
  // the air shower simulation program (ie CORSIKA)
  const CoordinateSystemPtr coordSys = shower.GetGroundParticleCoordinateSystem();
  const Point pCORSIKA(0, 0, 0, coordSys);
  const Vector zCORSIKA(0, 0, 1, coordSys);
  const CoordinateSystemPtr localSys = LocalCoordinateSystem::Create(pCORSIKA);
  const Vector zLocal(0, 0, 1, localSys);
  const double cosZenith = zCORSIKA.GetCosTheta(localSys);

  const double hObservationLevel = UTMPoint(pCORSIKA, ReferenceEllipsoid::eWGS84).GetHeight();

  const AtmosphereParameters& params = shower.GetAtmosphereParameters();
  fHlay = params.GetLayerAltitudes();
  fAatm = params.GetA();
  fBatm = params.GetB();
  fCatm = params.GetC();

  double refractiveIndexAtSeaLevel = 0;
  if (shower.HasRadioSimulation()) {
    const RadioSimulation& radioSim = shower.GetRadioSimulation();
    refractiveIndexAtSeaLevel = radioSim.GetRefractiveIndexAtSeaLevel();
    ostringstream info;
    info << "Event has RadioSimulation, use refractive index at sea level from .reas file: n = "
         << std::setprecision(10) << refractiveIndexAtSeaLevel;
    INFO(info);
  } else {
    refractiveIndexAtSeaLevel = kRefractiveIndexSeaLevel;
  }

  // make tabulated functions
  TabulatedFunctionErrors tabRIVsHeight;
  TabulatedFunctionErrors tabLogPressureVsHeight;
  TabulatedFunctionErrors tabLogVaporPressureVsHeight;
  TabulatedFunctionErrors tabTemperatureVsHeight;
  TabulatedFunctionErrors tabLogDepthVsHeight;
  TabulatedFunctionErrors tabLogDensityVsHeight;
  TabulatedFunctionErrors tabHeightVsLogDepth;

  const double minLog = kLn10 * numeric_limits<double>::min_exponent10;
  const double densityAtSeaLevel = GetDensityAtHeight(0);

  vector<double> heights;
  vector<double> depths;

  /* Loop just over layer altitudes turned out to be insufficent.
     With a more finer step-size results are much more accurate. */
  const double heightStep = 10 * meter;
  double height = -heightStep;
  while (height < fHlay.back()) {
    height += heightStep;

    if (height > fHlay.back())
      height = fHlay.back();

    const double heightCORSIKA = height;
    const double heightOffline = cosZenith * (heightCORSIKA - hObservationLevel) + hObservationLevel;

    const double depth = GetDepthAtHeight(heightCORSIKA);
    const double pressure = depth * atmosphere / kOverburdenSeaLevel;
    const double density = GetDensityAtHeight(heightCORSIKA);

    // by eye average, from GAP-2009-037
    const double avgrelhumidity = 27 * percent;

    // assuming dry, clean atmosphere
    // exploiting the molar ideal gas law

    const double temperature = pressure * kDryAirMolarMass / kMolarGasConstant / density;

    const double pSat = SaturationVaporPressure(temperature);
    // from GDASProfileModel.cc
    const double pVapor = avgrelhumidity * pSat;

    double refracIndex = 1;
    if (shower.HasRadioSimulation()) {
      // Gladstone-Dale law analogous to CoREAS
      refracIndex = RefractionIndex::GladstoneDale(density, densityAtSeaLevel, refractiveIndexAtSeaLevel);
    } else {
      // analogous to GDASProfileModel.cc
      refracIndex = RefractionIndex::LorentzLorentz(depth);
    }

    tabLogDepthVsHeight.PushBack(heightOffline, 0, log(depth), minLog);

    tabLogPressureVsHeight.PushBack(heightOffline, 0, log(pressure), minLog);

    tabLogDensityVsHeight.PushBack(heightOffline, 0, log(density), minLog);

    tabRIVsHeight.PushBack(heightOffline, 0., log(refracIndex), minLog);

    tabLogVaporPressureVsHeight.PushBack(heightOffline, 0, log(pVapor), minLog);

    tabTemperatureVsHeight.PushBack(heightOffline, 0, temperature, 0);

    heights.push_back(heightOffline);
    depths.push_back(depth);
  }

  // have to fill tabulated functions in ascending order or ordinate (tabulated function should be fixed..)
  for (int i = heights.size() - 1; i >= 0; --i)
    tabHeightVsLogDepth.PushBack(log(depths[i]), minLog, heights[i], 0);

  // create profile results
  *fTabLogPressureVsHeight =
    ProfileResult(tabLogPressureVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabLogVaporPressureVsHeight =
    ProfileResult(tabLogVaporPressureVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabTemperatureVsHeight =
    ProfileResult(tabTemperatureVsHeight, ProfileResult::eLinear, ProfileResult::eLinear);

  *fTabLogDensityVsHeight =
    ProfileResult(tabLogDensityVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabLogDepthVsHeight =
    ProfileResult(tabLogDepthVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabHeightVsLogDepth =
    ProfileResult(tabHeightVsLogDepth, ProfileResult::eLog, ProfileResult::eLinear);

  *fTabRIVsHeight =
     ProfileResult(tabRIVsHeight, ProfileResult::eLinear, ProfileResult::eLog);
}


const ProfileResult&
SimShowerProfileModel::EvaluatePressureVsHeight()
  const
{
  LazyInit();
  return *fTabLogPressureVsHeight;
}


const ProfileResult&
SimShowerProfileModel::EvaluateVaporPressureVsHeight()
  const
{
  LazyInit();
  return *fTabLogVaporPressureVsHeight;
}


const ProfileResult&
SimShowerProfileModel::EvaluateTemperatureVsHeight()
  const
{
  LazyInit();
  return *fTabTemperatureVsHeight;
}


const ProfileResult&
SimShowerProfileModel::EvaluateDensityVsHeight()
  const
{
  LazyInit();
  return *fTabLogDensityVsHeight;
}


const ProfileResult&
SimShowerProfileModel::EvaluateDepthVsHeight()
  const
{
  LazyInit();
  return *fTabLogDepthVsHeight;
}


const ProfileResult&
SimShowerProfileModel::EvaluateHeightVsDepth()
  const
{
  LazyInit();
  return *fTabHeightVsLogDepth;
}


const ProfileResult&
SimShowerProfileModel::EvaluateRefractionIndexVsHeight()
  const
{
  LazyInit();
  return *fTabRIVsHeight;
}


double
SimShowerProfileModel::GetDensityAtHeight(const double height)
  const
{
  const int n = fHlay.size();
  for (int l = 0; l < n-1; ++l)
    if (fHlay[l] <= height && height < fHlay[l+1])
      return fBatm[l] * exp(-height/fCatm[l]) / fCatm[l];
  if (height >= fHlay[n-1])
    return fBatm[n-1] / fCatm[n-1];
  return 0;
}


double
SimShowerProfileModel::GetDepthAtHeight(const double height)
  const
{
  int section = 0;
  const int n = fHlay.size();
  for (int l = 0; l < n-1; ++l) {
    if (fHlay[l] <= height && height < fHlay[l+1]) {
      section = l;
      break;
    }
  }
  if (height >= fHlay[n-1])
    section = n-1;

  const double hc = height / fCatm[section];
  if (section < 4)
    return fAatm[section] + fBatm[section] * exp(-hc);
  return fAatm[section] - fBatm[section] * hc;
}