AirflyFluorescenceModel.cc

Go to the documentation of this file.

#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/AirflyFluorescenceModel.h>

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


void
AirflyFluorescenceModel::Init()
{

  Branch topB = 
    CentralConfig::GetInstance()->GetTopBranch("AirflyFluorescenceModel");

  topB.GetChild("wavelength").GetData(fWavelength);
  topB.GetChild("relativeIntensity").GetData(fRelativeIntensity);
  topB.GetChild("relativeIntensityError").GetData(fRelativeIntensityError);
  topB.GetChild("pPrime").GetData(fPPrime);
  topB.GetChild("pPrimeErrorUncorrelated").GetData(fPPrimeErrorUncorrelated);
  topB.GetChild("pPrimeErrorCorrelated").GetData(fPPrimeErrorCorrelated);
  topB.GetChild("temperature0").GetData(fTemperature0);
  topB.GetChild("pressure0").GetData(fPressure0);
  topB.GetChild("yield337").GetData(fYield337);

  /* Extension of module from B. Keilhauer, Sep. 2009 according to 
     NIM A597 (2008) 50 and presentation from M. Bohacova at 6th Fl.
     Workshop in L'Aquila:
  */ 
  
  string tempParam, HumParam; 
  
  topB.GetChild("collisionalCrossSection").GetData(tempParam);
  if (tempParam == "eAIRFLYtemp")
    fTempParam = eAIRFLYtemp;
  else if (tempParam == "eNONE")
    fTempParam = eNoTempParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }
  
  topB.GetChild("humidity").GetData(HumParam);
  if (HumParam == "eAIRFLYhum")
    fHumParam = eAIRFLYhum;
  else if (HumParam == "eNONE")
    fHumParam = eNoHumParam;
  else {
    ERROR("Parametrization not implemented");
    throw utl::NonExistentComponentException();
  }
  
  //end of Extension

}


const
utl::TabulatedFunction&
AirflyFluorescenceModel::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);

  double fluorescenceYield = 0;

  /* Extension of module from B. Keilhauer, Sep. 2009 according to 
     NIM A597 (2008) 50 and presentation from M. Bohacova at 6th Fl.
     Workshop in L'Aquila; the given values are the measured ones 
     from Martina's presentation:
  */ 
  const int nWaveLength = 34;
  const double alphaAirfly[nWaveLength]={0.,0.,0.,0.,0.,-0.13,-0.19,0.,0.,0.,
                                         0.,0.,-0.35,0.,-0.38,-0.22,-0.35,0.,
                                         0.,-0.24,-0.17,-0.34,0.,0.,0.,-0.79,
                                         -0.2,-0.2,-0.37,0.,0.,0.,0.,-0.54};
  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.,0.,0.,0.,0.,0.,0.,0.,
                                       0.,0.,0.,0.};
  const double pPrimeHumAirfly[nWaveLength]={0.,0.,0.,0.,0.,1.2*100*pascal,
                                             1.1*100*pascal,0.,0.,0.,
                                             0.,0.,1.28*100*pascal,0.,
                                             1.5*100*pascal,1.27*100*pascal,
                                             1.3*100*pascal,0.,0.,
                                             1.3*100*pascal,
                                             1.1*100*pascal,1.4*100*pascal,
                                             0.,0.,0.,0.33*100*pascal,
                                             1.2*100*pascal,1.1*100*pascal,
                                             1.5*100*pascal,0.,
                                             0.,0.,0.,0.6*100*pascal};
  const double pPrimeHumNone[nWaveLength]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
                                           0.,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 == eAIRFLYtemp) ? alphaAirfly : alphaNone;
  const double * pPrimeHum =  (fHumParam == eAIRFLYhum) ? pPrimeHumAirfly : pPrimeHumNone;
  
  //end of Extension
  
  fFluorescenceSpectrum.Clear();
  for (unsigned int iwl = 0; iwl < fWavelength.size(); ++iwl) {
    
    const double wavelength = fWavelength[iwl];
    
    /* Extension of module from B. Keilhauer, Sep. 2009 according to 
       NIM A597 (2008) 50 and presentation from M. Bohacova at 6th Fl.
       Workshop in L'Aquila:
    */ 
    
    //    double Edep = GetdEdX0();
    
    if(pPrimeHum[iwl]==0.) {
      fluorescenceYield = density * fYield337 * fRelativeIntensity[iwl] * 0.01
        * (1 + fPressure0/fPPrime[iwl])
        / (1 + pressure/(fPPrime[iwl]*pow((fTemperature0/temperature), (alpha[iwl]-0.5))));
    }
    else {
      fluorescenceYield = density * fYield337 * fRelativeIntensity[iwl] * 0.01
        * ( (1 + fPressure0/fPPrime[iwl])
            / (1
               + (pressure/(fPPrime[iwl]*pow((fTemperature0/temperature), (alpha[iwl]-0.5))))
               - (vaporPressure/(fPPrime[iwl]*pow((fTemperature0/temperature), (alpha[iwl]-0.5))))
               + (vaporPressure/(pPrimeHum[iwl]*sqrt(temperature/fTemperature0))) 
               )
            );   
    } 
    
    //end of Extension
    
    fFluorescenceSpectrum.PushBack(wavelength, fluorescenceYield);
    
  }
  
  return fFluorescenceSpectrum;
}


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


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