AirFluorescenceModel.cc

Go to the documentation of this file.

#include <sstream>
#include <fstream>
#include <utl/TabulatedFunction.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/PhysicalConstants.h>
#include <fwk/CentralConfig.h>
#include <det/Detector.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <atm/AirFluorescenceModel.h>

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


void
AirFluorescenceModel::Init()
{
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("AirFluorescenceModel");

  string dataSetSelection;
  topB.GetChild("DataSetSelection").GetData(dataSetSelection);

  Branch dataSetBranch;

//#warning Kakimoto and Keilhauer data sets are not fully implemented

  if (dataSetSelection == "eAuger")
    dataSetBranch = topB.GetChild("AugerFluorescenceYieldDataSet");
  else if (dataSetSelection == "eNagano")
    dataSetBranch = topB.GetChild("NaganoFluorescenceYieldDataSet");
  else if (dataSetSelection == "eKakimoto")
    dataSetBranch = topB.GetChild("KakimotoFluorescenceYieldDataSet");
  else if (dataSetSelection == "eKeilhauer")
    dataSetBranch = topB.GetChild("KeilhauerFluorescenceYieldDataSet");
  else if (dataSetSelection == "eUserDefined")
    dataSetBranch = topB.GetChild("UserDefinedFluorescenceYieldDataSet");
  else {
    ERROR("DataSet not implemented");
    throw utl::NonExistentComponentException();
  }

  dataSetBranch.GetChild("wavelength").GetData(fWavelength);
  dataSetBranch.GetChild("relativeIntensity").GetData(fRelativeIntensity);
  dataSetBranch.GetChild("relativeIntensityError").GetData(fRelativeIntensityError);
  dataSetBranch.GetChild("pPrime").GetData(fPPrime);
  dataSetBranch.GetChild("pPrimeErrorUncorrelated").GetData(fPPrimeErrorUncorrelated);
  dataSetBranch.GetChild("pPrimeErrorCorrelated").GetData(fPPrimeErrorCorrelated);

  const string temp = dataSetBranch.GetChild("collisionalCrossSection").Get<string>();
  if (temp == "eAIRtemp")
    fTempParam = eAIRtemp;
  else if (temp == "eNONE")
    fTempParam = eNoTempParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }
  dataSetBranch.GetChild("alpha").GetData(fAlpha);

  const string water = dataSetBranch.GetChild("collisionalCrossSectionWater").Get<string>();
  if (water == "eAIRWatertemp")
    fTempWaterParam = eAIRWatertemp;
  else if (water == "eNONE")
    fTempWaterParam = eNoWaterTempParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }

  dataSetBranch.GetChild("alphaWater").GetData(fAlphaWater);

  const string hum = dataSetBranch.GetChild("humidity").Get<string>();
  if (hum == "eAIRhum")
    fHumParam = eAIRhum;
  else if (hum == "eNONE")
    fHumParam = eNoHumParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }

  dataSetBranch.GetChild("pPrimeWater").GetData(fPPrimeHum);

  dataSetBranch.GetChild("temperature0").GetData(fTemperature0);
  dataSetBranch.GetChild("pressure0").GetData(fPressure0);
  dataSetBranch.GetChild("yield337").GetData(fYield337);

  ostringstream info;
  info << "Version:\n"
          "Parameters:\n"
          "  DataSet:         " << dataSetSelection << "\n"
          "  yield337:        " << fYield337/(1/MeV) << " [1/MeV]\n"
          "  ref temperature: " << fTemperature0/kelvin << " [kelvin]\n"
          "  ref pressure:    " << fPressure0/100/pascal << " [100*pascal]\n"
          "  temp dependence: " << temp << "\n"
          "  hum dependence:  " << hum << "\n"
          "  water temp dep:  " << water << '\n';
  INFO(info);
}


const
utl::TabulatedFunction&
AirFluorescenceModel::EvaluateFluorescenceYield(const double heightAboveSeaLevel)
  const
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  const ProfileResult& tempProfile = atmo.EvaluateTemperatureVsHeight();
  const ProfileResult& densityProfile = atmo.EvaluateDensityVsHeight();
  const ProfileResult& pressureProfile = atmo.EvaluatePressureVsHeight();
  const ProfileResult& pvaporProfile = atmo.EvaluateVaporPressureVsHeight();

  const double temperature = tempProfile.Y(heightAboveSeaLevel);
  const double density = densityProfile.Y(heightAboveSeaLevel);
  const double pressure = pressureProfile.Y(heightAboveSeaLevel);
  const double vaporPressure = pvaporProfile.Y(heightAboveSeaLevel);

  fFluorescenceSpectrum.Clear();
  for (unsigned int iwl = 0; iwl < fWavelength.size(); ++iwl) {

    const double current_alpha = (fTempParam == eAIRtemp) ? fAlpha[iwl] : 0;
    const double current_alpha_w = (fTempWaterParam == eAIRWatertemp) ? fAlphaWater[iwl] : 0;
    const double current_pprime_w = (fHumParam == eAIRhum)  ? fPPrimeHum[iwl] : 0;

    double pprime = fPPrime[iwl] * pow(fTemperature0/temperature, current_alpha - 0.5);
    // P'(lambda,T) = P'(lambda,T_0)(T_0/T)^(alpha-1/2)

    const double pprime_w = current_pprime_w*pow(fTemperature0/temperature, current_alpha_w - 0.5);
    /* In a near future, alpha values for water might be available. Presently alpha_w = 0 for all wavelengths.
       P'_w(lambda,T) = P'_w(lambda,T_0)(T_0/T)^(alpha_w-1/2) */

    if (fHumParam == eAIRhum) {
      if (pprime_w) {
        pprime = (1/pprime) * (1 - vaporPressure/pressure) + vaporPressure / (pressure*pprime_w);
        pprime = 1 / pprime;
      }
    }

    // Y(lambda,P,T) = Y(lambda,P_0,T_0) * (1 + P_0/(P'(lambda,T_0))) / (1+P/(P'(lambda,T)))
    const double fluorescenceYield =
      (fYield337 * fRelativeIntensity[iwl] * 0.01) *
      (1 + fPressure0/fPPrime[iwl]) / (1 + pressure/pprime) *
      density;

    fFluorescenceSpectrum.PushBack(fWavelength[iwl], fluorescenceYield);

  }

  return fFluorescenceSpectrum;
}


double
AirFluorescenceModel::GetdEdX0()
  const
{
  // Energy loss evaluated for an electron of 0.85 MeV
  return 1;  //.658*MeV/(g/cm/cm);
}