FdLaserEnergyReconstructor.cc

Go to the documentation of this file.

/*
  implementation FdLaserEnergyReconstructorKG

  \author Joachim Debatin

 */

#include "FdLaserEnergyReconstructor.h"

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/PhysicalConstants.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/MathConstants.h>
#include <utl/PairErr.h>

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


#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/FdConstants.h>
#include <fevt/Pixel.h>

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

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Camera.h>

#include <atm/Atmosphere.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>

#include <vector>
#include <cmath>

using namespace FdLaserEnergyReconstructorKG;
using namespace std;
using namespace utl;
using namespace fwk;
using namespace fevt;
using namespace evt;
using namespace det;


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

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

  if (!topB) {
    ERROR("Could not find branch FdLaserEnergyReconstructorKG");
    return eFailure;
  }

  topB.GetChild("laserwavelength").GetData(fLaserWaveLength);
  return eSuccess;

  topB.GetChild("pixelswithpeak").GetData(fPixelsWithPeak);
  return eSuccess;
}


fwk::VModule::ResultFlag
FdLaserEnergyReconstructor::Run(evt::Event& event)
{
  /*if (event.GetHeader().GetTime().GetGPSNanoSecond() < 249000000. || event.GetHeader().GetTime().GetGPSNanoSecond() > 252000000.) {
    const string msg = "Probably no CLF event";
    INFO(msg);
    return eContinueLoop;
  }*/

  if (!event.HasFEvent())
    return eSuccess;

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

  //loop eyes
  for (FEvent::EyeIterator eyeiter = event.GetFEvent().EyesBegin(ComponentSelector::eHasData);
       eyeiter != event.GetFEvent().EyesEnd(ComponentSelector::eHasData); ++eyeiter) {

    Eye& eye = *eyeiter;

    if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
      cout << " no recData for eye " << eye.GetId() << " --> skip! " << endl;
      continue;
    }

    GetGeometryData(eye);

    evt::ShowerFRecData& shower = eye.GetRecData().GetFRecShower();

    //init slant depth profile
    try {
      atmos.InitSlantProfileModel(fCore, fAxis, 10*g/cm2);
    } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
      ERROR(e.GetMessage());
      return eFailure;
    }

    const atm::ProfileResult& slantDepthProfile = atmos.EvaluateSlantDepthVsDistance();
    const double atmMinDistance = slantDepthProfile.MinX() + 1*utl::millimeter;
    const double atmMaxDistance = slantDepthProfile.MaxX() - 1*utl::millimeter;

    //reuse stuff to show Laser in the EventBrowser
    if (!shower.HasEnergyDeposit())
      shower.MakeEnergyDeposit();
    TabulatedFunctionErrors& dedxProfile = shower.GetEnergyDeposit();

    if (!shower.HasLongitudinalProfile())
      shower.MakeLongitudinalProfile();
    TabulatedFunctionErrors& longProfile = shower.GetLongitudinalProfile();

    //for average energy
    double sum = 0;
    double wsum = 0;

    //loop telescopes
    for (fevt::Eye::TelescopeIterator teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
         teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {

      const Telescope& tel = *teliter;

      const fdet::Telescope& detTel  = detector.GetFDetector().GetEye(eye).GetTelescope(tel.GetId());

      GetTelescopeData(tel, detTel);

       //This can be used to cut on clouds
      int pixels = 0;
      for (fevt::Telescope::ConstPixelIterator pixiter = tel.PixelsBegin(ComponentSelector::eHasData);
           pixiter != tel.PixelsEnd(ComponentSelector::eHasData); ++pixiter) {

        const Pixel& pix = *pixiter;
        if (pix.HasRecData() && pix.GetRecData().PulseIsFound()) {
          ++pixels;
        }
      }
      if (pixels < fPixelsWithPeak) {
        string msg = "Telescope with to few pixels";
        INFO(msg);
        continue;
      }

      const double eff = detTel.GetMeasuredRelativeEfficiency(fLaserWaveLength);

      if (!eye.GetRecData().HasLightFlux()) {
        cout << "No LightFlux for this Eye ---> skip" << endl;
        continue;
      }

      const utl::TabulatedFunctionErrors& FluxAtTelescope = tel.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);
      utl::ConstTabulatedFunctionErrIterator it = FluxAtTelescope.Begin();

      //loop timebins
      for (int i = 1; i <= int(fTraceLength); ++i) {

        const int midTimeBin = int(i * fBinSize - 0.5 * fBinSize);

        // check for FieldOfView by ApertureLightFinder
        const utl::PairErr& bin = *it;
        if (int(bin.X()) == midTimeBin) {

          if (!tel.GetRecData().IsSpotFarFromBorder(bin.X())) {
            ++it;
            continue;
          }

          const utl::Point pointAtLaserBeam = CalculateLaserBeamPosition(double(i));

          const double Coefficient = CalculateAtmosphereCoefficient(atmos, pointAtLaserBeam, eff, i);

          const double recEnergy = bin.Y() * utl::kPlanck * utl::kSpeedOfLight / fLaserWaveLength / Coefficient;
          const double recEnergyErr = bin.YErr() * utl::kPlanck * utl::kSpeedOfLight / fLaserWaveLength / Coefficient;

          //Get slant depth
          const double distanceToCore = (pointAtLaserBeam - fCore).GetMag();
          if (distanceToCore < atmMinDistance || distanceToCore > atmMaxDistance) {
            cerr << "Error - point outside atmosphere, d=" << distanceToCore << endl;
            continue;
          }

          const double X = slantDepthProfile.Y(distanceToCore);

          const double w = 1 / pow(recEnergyErr, 2);
          sum += w * recEnergy;
          wsum += w;

          //Fill profiles EventBrowser
          dedxProfile.PushBack(X, 0, recEnergy/g*cm2, recEnergyErr/g*cm2);
          longProfile.PushBack(X, 0, recEnergy, recEnergyErr);
          ++it;
        }

      }

    }
    //Calculate average Energy
    const double avEnergy = sum / wsum;
    const double avEnergyErr = sqrt(1 / wsum);

    if (std::isnan(avEnergy) || std::isnan(avEnergyErr) || avEnergy < 1 || avEnergyErr < 1) {
      cout << "Energy Reconstruction failed! --> skip" << endl;
      continue;
    }

    shower.SetTotalEnergy(avEnergy, avEnergyErr);
    //cout << " average energy " << avEnergy << endl;

    //calculate chi squared
    unsigned int ndf = 0;
    double chisquared = 0;
    bool setchisquared = false;
    if (longProfile.GetNPoints() > 0) {
      utl::TabulatedFunctionErrIterator it = longProfile.Begin();
      while (ndf < longProfile.GetNPoints()) {
        const utl::PairErr& bin = *it;
        chisquared += pow(((bin.Y() - avEnergy) / bin.YErr()), 2);
        ++ndf;
        ++it;
      }
      setchisquared = true;
    }

    //needed for EventBrowser to show energy versus slant depth
    if (!shower.HasGHParameters()) {
      evt::GaisserHillas4Parameter ghTmp;
      shower.MakeGHParameters(ghTmp);
    }
    GaisserHillas4Parameter& gh4 =
      dynamic_cast<GaisserHillas4Parameter&>(shower.GetGHParameters());

    gh4.SetShapeParameter(gh::eX0, -50000*g/cm2, 0);
    gh4.SetShapeParameter(gh::eLambda, 6000*g/cm2, 0);

    if (setchisquared) {
      gh4.SetChiSquare(chisquared, ndf-1);
    } else {
      gh4.SetChiSquare(1, 1);
    }

  }

  return eSuccess;
}


utl::Point
FdLaserEnergyReconstructor::CalculateLaserBeamPosition(const double bin)
{
  //Get absolute time for bin + half of a binSize for middle of timeBin
  const utl::TimeStamp binTime =
    fGPSTime - TimeInterval(((fTraceLength - bin) * fBinSize - (fBinSize / 2) + static_cast<double>(fTelTimeOffset)));

  //time between core and bin in nanoseconds
  const double deltaT = (binTime - fCoreTime).GetInterval();

  //GetVector pointing from Eye to Core
  const utl::Vector eyeToCore = fCore - fTelPos;

  const double factor = (eyeToCore * eyeToCore - utl::kSpeedOfLight2 * deltaT * deltaT)/(2 * (eyeToCore * fAxis - utl::kSpeedOfLight * deltaT));

  return fCore - factor * fAxis;
}


//Calculate Coefficient for Energy = c * Signal
double
FdLaserEnergyReconstructor::CalculateAtmosphereCoefficient(const atm::Atmosphere& atmos,
                                                           const utl::Point pointAtLaserBeam,
                                                           const double efficiency,
                                                           const double i)
{
  //Point A and B are the points at the laser beam for beginning and end of one time bin at the telescope
  const utl::Point pointAtLaserBeamA = CalculateLaserBeamPosition(i - 0.5);
  const utl::Point pointAtLaserBeamB = CalculateLaserBeamPosition(i + 0.5);
  const utl::Vector eyeToPointLaser = pointAtLaserBeam - fTelPos;
  const double angle = Angle(-eyeToPointLaser, fAxis);

  //Transmission Factors
  const double TRay1 = atmos.EvaluateRayleighAttenuation(fCore, pointAtLaserBeam, fLaserWaveLength);
  const double TRay2 = atmos.EvaluateRayleighAttenuation(pointAtLaserBeam, fTelPos, fLaserWaveLength);
  const double TMie1 = atmos.EvaluateMieAttenuation(fCore, pointAtLaserBeam, fLaserWaveLength);
  const double TMie2 = atmos.EvaluateMieAttenuation(pointAtLaserBeam, fTelPos, fLaserWaveLength);

  //Scattering Factors
  const double SRay = atmos.EvaluateRayleighScattering(pointAtLaserBeamA, pointAtLaserBeamB, angle, eyeToPointLaser.GetMag(), fLaserWaveLength);
  const double SMie = atmos.EvaluateMieScattering(pointAtLaserBeamA, pointAtLaserBeamB, angle, eyeToPointLaser.GetMag(), fLaserWaveLength);
  const double STot = SMie + SRay;

  const double Coefficient = TRay1 * TRay2 * TMie1 * TMie2 * STot * efficiency;

  return Coefficient;
}


void
FdLaserEnergyReconstructor::GetGeometryData(const fevt::Eye& eye)
{
  const evt::ShowerFRecData& shower = eye.GetRecData().GetFRecShower();

  fCore = shower.GetCorePosition();
  fAxis = shower.GetAxis();
  fCoreTime = shower.GetCoreTime();
  fGPSTime = eye.GetHeader().GetTimeStamp();
}


void
FdLaserEnergyReconstructor::GetTelescopeData(const fevt::Telescope& tel,
                                             const fdet::Telescope& detTel)
{
  fTelPos = detTel.GetPosition();
  fBinSize = detTel.GetCamera().GetFADCBinSize();
  fTraceLength = detTel.GetCamera().GetFADCTraceLength();
  fTelTimeOffset = tel.GetTimeOffset();
}


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