FdProfileFinder.cc

Go to the documentation of this file.

#include "FdProfileFinder.h"

#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Vector.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Reader.h>

#include <utl/TabulatedFunctionErrors.h>
#include <utl/UTMPoint.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeRecData.h>

#include <det/Detector.h>
#include <atm/Atmosphere.h>
#include <atm/AttenuationResult.h>
#include <atm/ProfileResult.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>
#include <fdet/Channel.h>

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


using namespace FdProfileFinderOG;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace det;
using namespace fdet;
using namespace std;
using namespace atm;


FdProfileFinder::FdProfileFinder()
{
}

fwk::VModule::ResultFlag
FdProfileFinder::Init()
{
  CentralConfig* cc = CentralConfig::GetInstance();

  Branch topB = cc->GetTopBranch("FdProfileFinder");
  topB.GetChild("energyCutoff").GetData(fEnergyCutoff);

  return eSuccess;
}

fwk::VModule::ResultFlag
FdProfileFinder::Finish()
{
  return eSuccess;
}

fwk::VModule::ResultFlag
FdProfileFinder::Run(evt::Event& event)
{
  if (!event.HasFEvent())
    return eSuccess;

  FEvent::EyeIterator eyeiter;
  for (eyeiter= event.GetFEvent().EyesBegin();
       eyeiter != event.GetFEvent().EyesEnd();
       ++eyeiter)
  {                    // loop on showers
    fevt::Eye& eye = *eyeiter;

    if (! eye.GetRecData().HasLightFlux())
      continue;

    CoordinateSystemPtr eyeCS =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetEyeCoordinateSystem();

    Point eyepos =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();

    if (!eye.HasRecData()) {
      INFO("Eye has no geometry info");
      continue;
    }

    // retrieve geometrical parameters
    fAxis = eye.GetRecData().GetFRecShower().GetAxis();

    fT0 = eye.GetRecData().GetTZero();
    fChi0 =   eye.GetRecData().GetChiZero();
    fRp =     eye.GetRecData().GetRp();

    // use geometry as seen by this eye.
    // The vector  fCore_Eye_vec needs to be in the
    // horizontal plane of the eye
    //

    Vector vertical(0,0,1,eyeCS);
    Vector sdp = eye.GetRecData().GetSDP();
    Vector horizontalInSDP = cross(sdp,vertical);
    horizontalInSDP.Normalize();
    double core_distance = fRp / sin( kPi - fChi0);

    fCore_Eye_vec =  core_distance * horizontalInSDP;
    fCore = eyepos + fCore_Eye_vec;

    LightAtApertureToSize(eye);

  } // for eye
  return eSuccess;
}




void
FdProfileFinder::LightAtApertureToSize(fevt::Eye& eye)
{
  const ReferenceEllipsoid& wgs84 =
      ReferenceEllipsoid::Get (ReferenceEllipsoid::eWGS84);

  if (!eye.GetRecData().HasFRecShower())
    eye.GetRecData().MakeFRecShower();
  evt::ShowerFRecData& frecdata = eye.GetRecData().GetFRecShower();

  // Book energy cutoff
  frecdata.SetEnergyCutoff(fEnergyCutoff);

  double showerMax = 800*g/cm/cm;
  if (frecdata.HasGHParameters())
    showerMax = frecdata.GetGHParameters().GetXMax();

  CoordinateSystemPtr coreCS=
      fwk::LocalCoordinateSystem::Create(fCore);

  const Atmosphere& atm = det::Detector::GetInstance().GetAtmosphere();
  double dEdX0 = atm.GetdEdX0();

  ProfileResult dvh = atm.EvaluateDepthVsHeight();
  ProfileResult tvh = atm.EvaluateTemperatureVsHeight();
  ProfileResult denvh = atm.EvaluateDensityVsHeight();
  double atmHeightMin = max( tvh.MinX(), max(dvh.MinX(), denvh.MinX()));
  double atmHeightMax = min( tvh.MaxX(), max(dvh.MaxX(), denvh.MaxX()));

  if (!eye.GetRecData().HasLightFlux(fevt::FdConstants::eFluorTotal) ) return ;

  const TabulatedFunctionErrors& lightflux =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eFluorTotal);

  TabulatedFunctionErrors::ConstIteratorErr iter;

  // BRD Using the model efficiency of the first telescope as the eye
  // efficiency
  fevt::Eye::TelescopeIterator itTel = eye.TelescopesBegin();
  const fdet::Telescope& dettel =
    det::Detector::GetInstance().GetFDetector().GetTelescope(*itTel);
  const double refLambda =
    det::Detector::GetInstance().GetFDetector().GetReferenceLambda();
  // detector efficiency at reference wavelength
  const double detector_eff_ref = dettel.GetModelMeanEfficiency(dettel.GetConfigSignature("TelescopeSimulatorKG"), refLambda);

  // clear long profile if it exists
  if (frecdata.HasLongitudinalProfile())
    frecdata.GetLongitudinalProfile().Clear();
  if (frecdata.HasFluorescencePhotons())
    frecdata.GetFluorescencePhotons().Clear();

  for (iter = lightflux.Begin(); // loop on light profile at diaphragm
       iter != lightflux.End();
       ++iter)
  {
    double time      =     iter->X();
    double dtime     =     iter->XErr();
    double ph_at_dia =     iter->Y();
    double ph_at_dia_err = iter->YErr();

    double chi_i = fChi0 - 2* atan ( kSpeedOfLight/fRp *(time - fT0) );
    double l = fCore_Eye_vec.GetMag() * sin (chi_i) / sin (fChi0 - chi_i);
    if (l<0) l*=-1;

    // vector detector-shower at t=time
    Vector chi_i_vec = fCore_Eye_vec + l * fAxis;

    // point where the shower is at t=time
    Point chi_i_point = fCore + l * fAxis;

    // height of shower at t=time
    // note: we have to trasform to a UTM point to know
    // the z of the point above the ellipsoid
    double height = wgs84.PointToLatitudeLongitudeHeight(chi_i_point).get<2>();
    if (height>atmHeightMax || height<atmHeightMin) {
        //longitudinalprofile.PushBack (0,0,0,0);
        //photons_at_track.PushBack (0,0,0,0);
        continue;
    }

    // create a profile if none exists yet
    if (!frecdata.HasLongitudinalProfile() )
      frecdata.MakeLongitudinalProfile();
    TabulatedFunctionErrors& longitudinalprofile =
      frecdata.GetLongitudinalProfile() ;
    if (!frecdata.HasFluorescencePhotons())
      frecdata.MakeFluorescencePhotons();
    TabulatedFunctionErrors& photons_at_track=
      frecdata.GetFluorescencePhotons();

    double time1 = time/ns - dtime/ns;
    double time2 = time/ns + dtime/ns;

    double chi_1 = fChi0 - 2* atan ( kSpeedOfLight/fRp *(time1 - fT0) );
    double l1 = fCore_Eye_vec.GetMag() * sin (chi_1) / sin (fChi0 - chi_1);
    double chi_2 = fChi0 - 2* atan ( kSpeedOfLight/fRp *(time2 - fT0) );
    double l2 = fCore_Eye_vec.GetMag() * sin (chi_2) / sin (fChi0 - chi_2);

    // traveled track length
    double deltaL = abs (l2 - l1);

    // vertdepth
    double verticalDepth = dvh.Y(height);

    // slant depth
    double Xdepth =  verticalDepth / abs (cos(fAxis.GetTheta(coreCS)));

    double solidangle = 1 / (4.0 * kPi* chi_i_vec.GetMag2());

    // Evaluate fluorescence yield for that height
    const TabulatedFunction fluoyield_tab =
      atm.EvaluateFluorescenceYield(height);

    std::vector<double> relfluorescenceYield_vec( fluoyield_tab.GetNPoints());
    std::vector<double> waveLength_vec( fluoyield_tab.GetNPoints());

    // copy the fluo yield into the vector
    copy(fluoyield_tab.YBegin(),fluoyield_tab.YEnd(),
         relfluorescenceYield_vec.begin());
    // normalize the spectrum to 1
    double fluoyield_normalization = fluoyield_tab.SumY();

    double dEdX = EnergyDeposit(Xdepth, showerMax, fEnergyCutoff);

    // dEdX = 2;

    double totalyield = fluoyield_normalization * dEdX / dEdX0;

    for ( vector<double>::iterator it= relfluorescenceYield_vec.begin();
          it != relfluorescenceYield_vec.end();
          ++it)
      (*it) /= fluoyield_normalization;

    // copy the wavelength
    copy(fluoyield_tab.XBegin(), fluoyield_tab.XEnd(),
         waveLength_vec.begin());

    // Rayleigh attenuation
    Point eyepos =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();
    Point trackpos = eyepos + chi_i_vec;


    AttenuationResult rAttenuation =
      atm.EvaluateRayleighAttenuation(eyepos,trackpos,waveLength_vec);

    const utl::TabulatedFunctionErrors& rayleighAtt =
      rAttenuation.GetTransmissionFactor();

    // Mie Attenuation
    AttenuationResult mAttenuation =
      atm.EvaluateMieAttenuation(eyepos,trackpos,waveLength_vec);

    const utl::TabulatedFunctionErrors& mieAtt =
      mAttenuation.GetTransmissionFactor();

    // calculation of shower size
    double n_e = 0; // size
    double L =   0; // light at track
    double convert = 0; // conversion factor

    // wavelength dependent part
    for (unsigned int i = 0; i<waveLength_vec.size(); i++){
      convert += relfluorescenceYield_vec[i]* dettel.GetModelMeanEfficiency(dettel.GetConfigSignature("TelescopeSimulatorKG"), waveLength_vec[i])
        * rayleighAtt.Y(waveLength_vec[i]) *  mieAtt.Y(waveLength_vec[i]);
    } // for wavelength

    // non wavelength dependent part
    n_e =   ph_at_dia * detector_eff_ref/solidangle/deltaL/totalyield/convert;
    L =   ph_at_dia * detector_eff_ref/solidangle/convert;
    // n_e error: only error on number of photons at diaphragm considered
    double n_e_err= ph_at_dia_err* n_e / ph_at_dia;

    longitudinalprofile.PushBack(Xdepth,0,n_e,n_e_err);

    // L error: only error on number of photons at diaphragm considered
    double L_err= ph_at_dia_err* L / ph_at_dia;

    photons_at_track.PushBack(Xdepth,0,L,L_err);

  }// loop on light profile
}

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