Kakimoto1996FluorescenceModel.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/Vector.h>
#include <utl/PhysicalConstants.h>
#include <utl/Transformation.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/AugerException.h>
#include <det/Detector.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>
#include <atm/Kakimoto1996FluorescenceModel.h>
#include <fwk/CentralConfig.h>
#include <sstream>

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


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

  topB.GetChild("wavelength").GetData(fWavelength);
  topB.GetChild("fluorescenceYield").GetData(fFluorescenceYield);

  // Finding the bin that constains the 391.4nm line
  bool found = false;
  const double line = 391.4*nanometer;
  for (int i = 0; i < int(fWavelength.size()-1); ++i)
    if (fWavelength[i] < line && line < fWavelength[i+1]) {
      found = true;
      // Find closest bin
      if ((line - fWavelength[i]) < (fWavelength[i+1] - line)) {
        fBin391 = i;
        break;
      } else {
        fBin391 = i + 1;
        break;
      }
    }

  if (!found) {
    string msg("Supplied Fluorescence spectrum misses wavelength for 391.4 nm");
    FATAL(msg);
    throw NoDataForModelException(msg);
  }

  // Computing the total light
  fTotalLight = 0;
  for (int i = 0; i < int(fFluorescenceYield.size()); ++i)
    fTotalLight += fFluorescenceYield[i];
}


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


const utl::TabulatedFunction&
Kakimoto1996FluorescenceModel::EvaluateFluorescenceYield(const double heightAboveSeaLevel)
  const
{
  const Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  const ProfileResult& tempProfile    = atmo.EvaluateTemperatureVsHeight();
  const ProfileResult& densityProfile = atmo.EvaluateDensityVsHeight();
  const double temperature = tempProfile.Y(heightAboveSeaLevel);
  const double density     = densityProfile.Y(heightAboveSeaLevel);

  const double rho = density / (g/cm3);

  TabulatedFunction seaLevelFluorescenceSpectrum(fWavelength,
                                                 fFluorescenceYield);

  // Experimental parameters
  const double a1 = 0.929071;
  const double b1 = 1850.0;
  const double a2 = 0.574145;
  const double b2 = 6500.0;
  const double f  = 1.0439e-5;

  // Total light except from bin with the 391nm line
  const double k2 = fFluorescenceYield[fBin391];
  const double k = fTotalLight - k2;

  const double waveLengthBin391 = fWavelength[fBin391];

  fFluorescenceSpectrum.Clear();
  for (TabulatedFunction::Iterator it = seaLevelFluorescenceSpectrum.Begin();
       it != seaLevelFluorescenceSpectrum.End(); ++it) {

    const double wavelength = it->X();
    const double seaLevelFluorescenceYield = it->Y();

    double fluorescenceYield;
    if (wavelength != waveLengthBin391)
      fluorescenceYield = seaLevelFluorescenceYield*rho*a1 /
        ((1 + rho*b1*sqrt(temperature))*k*f);
    else
      fluorescenceYield = seaLevelFluorescenceYield*rho*a2 /
        ((1 + rho*b2*sqrt(temperature))*k2*f);

    fFluorescenceSpectrum.PushBack(wavelength, fluorescenceYield);
  }

  return fFluorescenceSpectrum;
}


double Kakimoto1996FluorescenceModel::GetdEdX0() const {

    // Energy loss evaluated for an electron of 1.4 MeV at standard
    // temperature and pressure.
    // Value taken from the ESTAR program of NIST (2004)
    // http://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html
    return 1.658*MeV/(g/cm/cm);
}

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