Keilhauer2008FluorescenceModel.cc

Go to the documentation of this file.

#include <utl/TabulatedFunction.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <fwk/CentralConfig.h>
#include <det/Detector.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <atm/Keilhauer2008FluorescenceModel.h>

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


void
Keilhauer2008FluorescenceModel::Init()
{
  Branch topB = 
    CentralConfig::GetInstance()->GetTopBranch("Keilhauer2008FluorescenceModel");
  
  string tempParam, HumParam; 

  topB.GetChild("collisionalCrossSection").GetData(tempParam);
  if (tempParam == "eAIRFLY")
    fTempParam = eAIRFLY;
  else if (tempParam == "eNONE")
    fTempParam = eNoTempParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }
  
  topB.GetChild("humidity").GetData(HumParam);
  if (HumParam == "eMorozov")
    fHumParam = eMorozov;
  else if (HumParam == "eTilo")
    fHumParam = eTilo;
  else if (HumParam == "eNONE")
    fHumParam = eNoHumParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }
  
  const int nWaveLength = 23;
  double wavelength[nWaveLength] = 
  {
    3117*angstrom, 3136*angstrom,3159*angstrom,3285*angstrom,
    3309*angstrom, 3339*angstrom,3371*angstrom,3469*angstrom,
    3500*angstrom, 3537*angstrom,3577*angstrom,3672*angstrom,
    3711*angstrom, 3756*angstrom,3805*angstrom,3894*angstrom,
    3914*angstrom, 3943*angstrom,3998*angstrom,4059*angstrom,
    4141*angstrom, 4201*angstrom,4270*angstrom
  }; 

  for (int i = 0; i < nWaveLength; ++i)
    fWavelength.push_back(wavelength[i]);
}


const
utl::TabulatedFunction&
Keilhauer2008FluorescenceModel::EvaluateFluorescenceYield(const double heightAboveSeaLevel)
  const
{

  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const ProfileResult& pressureProfile = atmo.EvaluatePressureVsHeight();
  const ProfileResult& tempProfile     = atmo.EvaluateTemperatureVsHeight();
  const ProfileResult& densityProfile  = atmo.EvaluateDensityVsHeight();
  const ProfileResult& pvaporProfile   = atmo.EvaluateVaporPressureVsHeight(); 
  const double pressure      = pressureProfile.Y(heightAboveSeaLevel);
  const double temperature   = tempProfile.Y(heightAboveSeaLevel);
  const double density       = densityProfile.Y(heightAboveSeaLevel);
  const double vaporPressure = pvaporProfile.Y(heightAboveSeaLevel);

  fFluorescenceSpectrum.Clear();

  const unsigned int nWaveLength=23;
  const double sigma_nn[nWaveLength]={1.2e-19*m2 ,8.8e-20*m2 ,3.77e-20*m2,1.2e-19*m2,
                                      8.8e-20*m2 ,3.77e-20*m2,1.82e-20*m2,1.2e-19*m2,
                                      8.8e-20*m2 ,3.77e-20*m2,1.82e-20*m2,1.2e-19*m2,
                                      8.8e-20*m2 ,3.77e-20*m2,1.82e-20*m2,1.2e-19*m2,
                                      4.37e-19*m2,8.8e-20*m2 ,3.77e-20*m2,1.82e-20*m2,
                                      1.2e-19*m2 ,8.8e-20*m2 ,3.77e-20*m2};
  
  const double sigma_no[nWaveLength]={8.e-19*m2 ,7.e-19*m2,5.e-19*m2,8.e-19*m2,7.0e-19*m2,
                                      5.e-19*m2,2.1e-19*m2,8.e-19*m2,7.e-19*m2,5.0e-19*m2,
                                      2.1e-19*m2,8.e-19*m2,7.e-19*m2,5.e-19*m2,2.1e-19*m2,
                                      8.e-19*m2,13.e-19*m2,7.e-19*m2,5.e-19*m2,2.1e-19*m2,
                                      8.e-19*m2 ,7.e-19*m2,5.e-19*m2};
  

  const double eff[nWaveLength]={0.005*percent*percent,0.029*percent,0.05*percent,0.0154*percent,0.002*percent,
                                 0.0041*percent,0.082*percent,0.0063*percent,0.004*percent,0.029*percent,
                                 0.0615*percent,0.0046*percent,0.01*percent,0.0271*percent,0.0213*percent,
                                 0.003*percent,0.33*percent,0.0064*percent,0.016*percent,0.0067*percent,
                                 0.0017*percent,0.0046*percent,0.0035*percent};
  
  const double tau[nWaveLength]={6.65e-8*s,4.45e-8*s,4.17e-8*s,6.65e-8*s,4.45e-8*s,4.17e-8*s,
                                 4.17e-8*s,6.65e-8*s,4.45e-8*s,4.17e-8*s,4.17e-8*s,6.65e-8*s,
                                 4.45e-8*s,4.17e-8*s,4.17e-8*s,6.65e-8*s,6.58e-8*s,4.45e-8*s,
                                 4.17e-8*s,4.17e-8*s,6.65e-8*s,4.45e-8*s,4.17e-8*s};
  
  
  const double alphaAirfly[nWaveLength]={0.,-0.09,-0.21,0.,-0.09,-0.21,-0.36,0.,-0.09,-0.21,-0.36,0.,
                                         -0.09,-0.21,-0.36,0.,-0.8,-0.09,-0.21,-0.36,0.,-0.09,-0.21};
  const double alphaNone[nWaveLength]={0.,0.,0.,0.,0.,0.,0.,0.,0.,
                                       0.,0.,0.,0.,0.,0.,0.,0.,0.,
                                       0.,0.,0.,0.,0.};
  
  const double * alpha =  (fTempParam == eAIRFLY) ? alphaAirfly : alphaNone;
  
  const double sigma_nvMorozov[nWaveLength]={0.*m2,0.*m2,8.91e-19*m2,0.*m2,0.*m2,8.91e-19*m2,
                                             9.45e-19*m2,0.*m2,0.*m2,8.91e-19*m2,9.45e-19*m2,
                                             0.*m2,0.*m2,8.91e-19*m2,9.45e-19*m2,0.*m2,0.*m2,
                                             0.*m2,8.91e-19*m2,9.45e-19*m2,0.*m2,0.*m2,8.91e-19*m2};
  const double sigma_nvTilo[nWaveLength]={0.*m2,0.*m2,7.71e-19*m2,0.*m2,0.*m2,7.71e-19*m2,
                                          7.25e-19*m2,0.*m2,0.*m2,7.71e-19*m2,7.25e-19*m2,
                                          0.*m2,0.*m2,7.71e-19*m2,7.25e-19*m2,0.*m2,0.*m2,
                                          0.*m2,7.71e-19*m2,7.25e-19*m2,0.*m2,0.*m2,7.71e-19*m2};
  const double sigma_nvNone[nWaveLength]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
                                          0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};

  const double * sigma_nv = (fHumParam == eMorozov) ?
    sigma_nvMorozov :
    (fHumParam == eTilo) ?
    sigma_nvTilo : sigma_nvNone;

  // Get water vapor pressure, then calculate the H20, N2, O2 volume fractions.
  const double volV = vaporPressure / pressure;
  const double volN = (1. - volV) * kN2AirFraction;
  const double volO = (1. - volV) * kO2AirFraction;
  
  const double preFactor = pressure*kAvogadro / (kMolarGasConstant*temperature)
                           * sqrt(kBoltzmann*temperature * kAvogadro / kPi);

  for (unsigned int iwl = 0; iwl < nWaveLength; ++iwl)
  {
    const double snn = sigma_nn[iwl] * pow(293 * kelvin, -alpha[iwl]);
    const double sno = sigma_no[iwl] * pow(293 * kelvin, -alpha[iwl]);
    
    const double nnQuenching = 4 * volN * snn * pow(temperature, alpha[iwl])
      * sqrt(1/kN2MolarMass);
    const double noQuenching = 2 * volO * sno * pow(temperature, alpha[iwl])
      * sqrt(2 * (1 / kN2MolarMass + 1 / kO2MolarMass));
    const double nvQuenching = 2 * volV * sigma_nv[iwl]
      * sqrt(2 * (1 / kN2MolarMass + 1 / kH2OMolarMass));
    const double pOverpPrime = tau[iwl] * preFactor 
      * (nnQuenching + noQuenching + nvQuenching);

    const double deff = eff[iwl] / (1 + pOverpPrime);
    const double dscn = deff * fWavelength[iwl] / (kPlanck * kSpeedOfLight);

    const double fluorescenceYield = dscn * GetdEdX0() * density;
    
    fFluorescenceSpectrum.PushBack(fWavelength[iwl], fluorescenceYield);
  }

  return fFluorescenceSpectrum;
}


double
Keilhauer2008FluorescenceModel::GetdEdX0()
  const
{
  // Energy loss evaluated for an electron of .85 MeV
  return 1.675*MeV/(g/cm/cm);  
}

const std::vector<double>&
Keilhauer2008FluorescenceModel::GetWavelengths() const{
  
  return fWavelength;

}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. -k"
// End: