FdCherenkovFinder.cc

Go to the documentation of this file.

#include <utl/ErrorLogger.h>
#include <utl/Vector.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Reader.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/ScatteringResult.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 <utl/Vector.h>
#include <utl/Point.h>

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

#include "FdCherenkovFinder.h"

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


FdCherenkovFinder::FdCherenkovFinder(){}

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

fwk::VModule::ResultFlag FdCherenkovFinder::Init(){

  fwk::CentralConfig* cc = fwk::CentralConfig::GetInstance();

  Branch topB =
    cc->GetTopBranch("FdCherenkovFinder");

  topB.GetChild("x0FOV").GetData(fX0);
  topB.GetChild("xMaxFOV").GetData(fXMax);
  topB.GetChild("stepX").GetData(fStep);
  topB.GetChild("size").GetData(fSize);

  fIter = 0;

  return eSuccess;

}// end of Init


fwk::VModule::ResultFlag FdCherenkovFinder::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.HasRecData()) continue;
    fevt::EyeRecData& eyerecdata = eye.GetRecData();

    if (!eyerecdata.HasFRecShower()) continue; // check added on May 6, 2005 (MAM)
    if (!eyerecdata.GetFRecShower().HasLongitudinalProfile() )continue;
    if (eyerecdata.GetFRecShower().GetLongitudinalProfile().
        GetNPoints()<1 ) continue;
    if (!eyerecdata.GetFRecShower().HasGHParameters() )continue;
    // GH fit failed !
    if (eyerecdata.GetFRecShower().GetGHParameters().GetXMax() < 0) continue;
    if (!eyerecdata.HasLightFlux(fevt::FdConstants::eTotal) )continue;
    if (!eyerecdata.HasLightFlux(fevt::FdConstants::eFluorTotal) )continue;

    fRecShower = &eyerecdata.GetFRecShower();

    // Loop control
    if(fIter >= 10)
      fIter = 0;

    fAxis = fRecShower->GetAxis();
    fCore = fRecShower->GetCorePosition();


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

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


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

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

    fCoreEyeVec =  fCoreDistance*horizontalInSDP;
    fCore = fEyePos + fCoreEyeVec;

    // Note: Since the measured efficiency is really a telescope property (albeit
    //       currently the same for all telescopes), this shouldn't simply grab the
    //       efficiency of the first telescope. --Steffen 18.5.09
    fEfficiencyCalibration =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetTelescope(1).GetMeasuredRelativeEfficiency();

    fWavelength.resize(fEfficiencyCalibration.GetNPoints());

    std::copy( fEfficiencyCalibration.XBegin(),
               fEfficiencyCalibration.XEnd(),
               fWavelength.begin()) ;

    CherenkovFinderAtAperture(eye);

    } // end of loop over eyes

  ++fIter;

  return eSuccess;

} // end of Run


fwk::VModule::ResultFlag
FdCherenkovFinder::CherenkovFinderAtAperture(fevt::Eye& eye) {

  const ReferenceEllipsoid& wgs84 =
      ReferenceEllipsoid::Get (ReferenceEllipsoid::eWGS84);


  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherDirect) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherDirect);

  TabulatedFunctionErrors& directCherenkovFlux =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherDirect);

  directCherenkovFlux.Clear();

  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherRayleighScattered) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherRayleighScattered);

  TabulatedFunctionErrors& rayScattCherenkovFlux =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherRayleighScattered);

  rayScattCherenkovFlux.Clear();

  if(!eye.GetRecData().HasLightFlux(fevt::FdConstants::eCherMieScattered) )
    eye.GetRecData().MakeLightFlux(fevt::FdConstants::eCherMieScattered);

  TabulatedFunctionErrors& mieScattCherenkovFlux =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eCherMieScattered);

  mieScattCherenkovFlux.Clear();


  TabulatedFunction& longitudinalProfile =
    fRecShower->GetLongitudinalProfile() ;

  double ghXMax = fRecShower->GetGHParameters().GetXMax();

  if (ghXMax < 0) return eSuccess;

  const Atmosphere& atm = det::Detector::GetInstance().GetAtmosphere();
  ProfileResult dvh = atm.EvaluateDepthVsHeight();
  double atmHeightMin = dvh.MinX();
  double atmHeightMax = dvh.MaxX();

  vector<double> cherenkovAtTrack = InitialCherenkov(eye);
  const TabulatedFunctionErrors& lightFlux =
    eye.GetRecData().GetLightFlux(fevt::FdConstants::eFluorTotal);


  CoordinateSystemPtr siteCS =
    det::Detector::GetInstance().GetSiteCoordinateSystem();

  int n = 0;

  double dirCkv      = 0.0;
  double rayScattCkv = 0.0;
  double mieScattCkv = 0.0;

  vector<double> vecDirCkv(fWavelength.size());
  vector<double> vecRayScattCkv(fWavelength.size());
  vector<double> vecMieScattCkv(fWavelength.size());

  const double lambdaRef  =
    det::Detector::GetInstance().GetFDetector().GetReferenceLambda();
  const double epsilonRef = fEfficiencyCalibration.Y(lambdaRef);

  TabulatedFunctionErrors::ConstIteratorErr iter;

  // loop on light profile at diaphragm
  for(iter = lightFlux.Begin(); iter != lightFlux.End(); ++iter) {

    double time  = iter->X();
    double dtime = iter->XErr();

    double chiAv = fChi0 - 2.* atan(kSpeedOfLight/fRp*(time - fT0));
    double l = fCoreEyeVec.GetMag()*sin(chiAv)/sin(fChi0 - chiAv);

    double time1 = time/ns - dtime/ns;
    double chi1 = fChi0 - 2.0*atan(kSpeedOfLight/fRp*(time1 - fT0));
    double l1 = fCoreEyeVec.GetMag()*sin(chi1)/sin(fChi0 - chi1);

    double time2 = time/ns + dtime/ns;
    double chi2 = fChi0 - 2.0*atan(kSpeedOfLight/fRp *(time2 - fT0));
    double l2 = fCoreEyeVec.GetMag()*sin(chi2)/sin(fChi0 - chi2);

    // Vector detector-shower at t=time
    Vector chiAvVec = fCoreEyeVec + l*fAxis;

    Point pos1   = fCore + l1*fAxis;
    Point pos2   = fCore + l2*fAxis;

    double pos1_z = (wgs84.PointToLatitudeLongitudeHeight(pos1)).get<2>();
    double pos2_z = (wgs84.PointToLatitudeLongitudeHeight(pos2)).get<2>();
    if ((pos1_z>atmHeightMax || pos1_z<atmHeightMin) ||
        (pos2_z>atmHeightMax || pos2_z<atmHeightMin)) {
        directCherenkovFlux.PushBack(time,0,0,0);
        rayScattCherenkovFlux.PushBack(time,0,0,0);
        mieScattCherenkovFlux.PushBack(time,0,0,0);
        continue;
    }

    // Average position along the track.
    double xpos = (pos1.GetX(siteCS) + pos2.GetX(siteCS))/2.0;
    double ypos = (pos1.GetY(siteCS) + pos2.GetY(siteCS))/2.0;
    double zpos = (pos1.GetZ(siteCS) + pos2.GetZ(siteCS))/2.0;

    Point avepos(xpos,ypos,zpos,siteCS);

    double xMax = fXMax;

    if(fIter != 0)
      xMax = ghXMax;

    const double slantDepth = longitudinalProfile.GetX(n);
    const double showerAge = ShowerAge(slantDepth, xMax);

    // Rayleigh attenuation in a atmospheric depth bin along the track
    AttenuationResult rAttenuation =
      atm.EvaluateRayleighAttenuation(pos1,pos2,fWavelength);

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

    // Mie attenuation in a atmospheric depth bin along the track
    AttenuationResult mAttenuation =
      atm.EvaluateMieAttenuation(pos1,pos2,fWavelength);

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

    // Integrated Cherenkov at track.
    const TabulatedFunction CherenkovProd =
      atm.EvaluateCherenkovPhotons(pos1,pos2,showerAge);

    for(unsigned int w = 0; w < fWavelength.size(); w++) {

     cherenkovAtTrack[w] *= mieAtt.GetY(w)*rayleighAtt.GetY(w);
     if (longitudinalProfile.GetY(n)>0.0)
       cherenkovAtTrack[w] += CherenkovProd.GetY(w)*longitudinalProfile.GetY(n);

    }

    // Distance between average point to the eye.
    double disteye = (avepos - fEyePos).GetMag();

    // Emisson angle from the average point to the eye.
    double angle    = fChi0 - chiAv;

    const utl::TabulatedFunction directCherenkov =
      atm.EvaluateCherenkovDirect(pos1,pos2,fEyePos,showerAge);


    // Rayleigh and Mie Attenuation in a atmospheric depth bin from shower
    // track to eye
    AttenuationResult rAttTrackEye =
      atm.EvaluateRayleighAttenuation(avepos,fEyePos,fWavelength);

    const utl::TabulatedFunctionErrors& rayleighAttTrackEye =
      rAttTrackEye.GetTransmissionFactor();

    AttenuationResult mAttTrackEye =
      atm.EvaluateMieAttenuation(avepos,fEyePos,fWavelength);

    const utl::TabulatedFunctionErrors& mieAttTrackEye =
      mAttTrackEye.GetTransmissionFactor();

    dirCkv = 0.0;

    for(unsigned int i=0; i<fWavelength.size(); i++) {

      vecDirCkv[i] = 0.0;
      if(longitudinalProfile.GetY(n)>0.0) {

        vecDirCkv[i]  = directCherenkov.GetY(i)*longitudinalProfile.GetY(n);
        vecDirCkv[i] *= mieAttTrackEye.GetY(i)*rayleighAttTrackEye.GetY(i);

        dirCkv += vecDirCkv[i]*fEfficiencyCalibration[i].Y()/epsilonRef;
      }

    }


    // Evaluate Scattering Mie and Rayleigh from shower track to eye
    ScatteringResult rScattering =
      atm.EvaluateRayleighScattering(pos1,pos2,angle,disteye,fWavelength);

    const utl::TabulatedFunctionErrors& rayleighScattTrackEye =
      rScattering.GetScatteringFactor();

    ScatteringResult mScattering =
      atm.EvaluateMieScattering(pos1,pos2,angle,disteye,fWavelength);

    const utl::TabulatedFunctionErrors& mieScattTrackEye =
      mScattering.GetScatteringFactor();

    // Scatttered Cherenkov light
    rayScattCkv  = 0.0;
    mieScattCkv  = 0.0;

    for(unsigned int i=0; i<fWavelength.size(); i++) {

      vecRayScattCkv[i] = 0.0;
      vecMieScattCkv[i] = 0.0;

      vecRayScattCkv[i] = cherenkovAtTrack[i]
                                 *rayleighScattTrackEye.GetY(i);
      vecMieScattCkv[i] = cherenkovAtTrack[i]
                                 *mieScattTrackEye.GetY(i);

      vecRayScattCkv[i] *= mieAttTrackEye.GetY(i)
                                     *rayleighAttTrackEye.GetY(i);
      vecMieScattCkv[i] *= mieAttTrackEye.GetY(i)
                                     *rayleighAttTrackEye.GetY(i);

      rayScattCkv += vecRayScattCkv[i]
                                *fEfficiencyCalibration[i].Y()/epsilonRef;
      mieScattCkv += vecMieScattCkv[i]
                                *fEfficiencyCalibration[i].Y()/epsilonRef;

    }

    directCherenkovFlux.PushBack(time,0,dirCkv,0);
    rayScattCherenkovFlux.PushBack(time,0,rayScattCkv,0);
    mieScattCherenkovFlux.PushBack(time,0,mieScattCkv,0);

    n++;

  } // end loop on light profile

  return eSuccess;
}


//BRD 28/4/05 new algorithm for calculating l1 and l2
vector<double>
FdCherenkovFinder::InitialCherenkov(fevt::Eye& /*eye*/)
  const
{
  const TabulatedFunction& longitudinalProfile =
    fRecShower->GetLongitudinalProfile() ;

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

  cout << "- " << longitudinalProfile.GetNPoints() << endl;

  double xInit    = longitudinalProfile.GetX(0);
  double sizeInit = longitudinalProfile.GetY(0);

  cout << "-- " << xInit << " " << sizeInit << endl;

  double ghXMax = fRecShower->GetGHParameters().GetXMax();

  int xBin = int((xInit - fX0)/fStep);
  if (xBin<0) xBin=0;

  double x1        = 0.0;
  double x2        = 0.0;
  double xMed      = 0.0;
  double xMax      = 0.0;
  double size      = 0.0;
  double normSize  = 0.0;

  vector<double> cherenkov;
  for(unsigned int w=0; w<fWavelength.size(); w++)
    cherenkov.push_back(0.);

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

  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));


  double cosZenithAngle = fAxis.GetCosTheta(coreCS);

  ProfileResult hvd = atm.EvaluateHeightVsDepth();
  double atmDepthMin = hvd.MinX();
  double atmDepthMax = hvd.MaxX();

  for(int i=0; i<xBin; i++) {

    //slant depths
    x1 = fX0 + (double(i) - 1.0)*fStep;
    x2 = fX0 + double(i)*fStep;

    xMed = (x1 + x2)*0.5;

    size = GaisserHillas(xMed,fX0,fXMax,fSize,kLambdaGH);
    normSize = GaisserHillas(xInit,fX0,fXMax,fSize,kLambdaGH);
    size *= sizeInit/normSize;

    xMax = fXMax;

    if(fIter)
      xMax = ghXMax;

    double depth1=x1*cosZenithAngle;
    double depth2=x2*cosZenithAngle;
    if( depth1<atmDepthMin ||  depth1>atmDepthMax ||
        depth2<atmDepthMin ||  depth2>atmDepthMax ) {
      continue;
    }


    double h1 = hvd.Y(depth1);
    double h2 = hvd.Y(depth2);

    //for defined heights within the molecular profile model
    if (h1 !=0. && h2 !=0.) {

      //vertical height above core

      h1 -= (wgs84.PointToLatitudeLongitudeHeight(fCore)).get<2>();
      h2 -= (wgs84.PointToLatitudeLongitudeHeight(fCore)).get<2>();

      double l1 = h1/cosZenithAngle;
      double l2 = h2/cosZenithAngle;


      Point pos1 = fCore + l1*fAxis;
      Point pos2 = fCore + l2*fAxis;

      const double showerAge = ShowerAge(xMed, xMax);

      TabulatedFunction cherenkovProd =
        atm.EvaluateCherenkovPhotons(pos1,pos2,showerAge);


      // Rayleigh Attenuation
      AttenuationResult rAttenuation =
        atm.EvaluateRayleighAttenuation(pos1,pos2,fWavelength);

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

      // Mie Attenuation
      AttenuationResult mAttenuation =
        atm.EvaluateMieAttenuation(pos1,pos2,fWavelength);

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

      for(unsigned int w=0; w<fWavelength.size(); w++) {

        cherenkov.at(w) *= mieAtt.GetY(w)*rayleighAtt.GetY(w);
        cherenkov.at(w) += cherenkovProd.GetY(w)*size;
      }
    } // for defined heights

  }

  return cherenkov;

} // end of InitialCherenkov


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