FD2ADSTUtil.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "FD2ADSTUtil.h"
#include "ConversionUtil.h"

#include <adst/RecEvent.h>

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

#include <det/VManager.h>
#include <det/Detector.h>
#include <det/ManagerRegister.h>

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

#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/RiseTime1000.h>

#include <evt/GaisserHillas4Parameter.h>
#include <evt/GaisserHillas6Parameter.h>

#include <fevt/FEvent.h>
#include <fevt/Header.h>
#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/EyeTriggerData.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeTriggerData.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <fevt/FdComponentSelector.h>

#include <atm/MeasuredDBMieModel.h>
#include <atm/GDASProfileModel.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>

#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/ReferenceEllipsoid.h>

//#warning FDEventLib/IoAuger should not be used in Offline2ADST converter
#include <AugerEvent.h>
#include <EyeEvent.hh>
#include <EyeEventHeader.hh>

#include <iostream>
#include <string>
#include <vector>

using namespace std;
using namespace otoa;
using namespace otoa::fd;
using namespace fevt;
using namespace utl;
using namespace atm;


std::map<int, bool>
otoa::fd::CalcMirrorTLTMap(const fevt::Eye& eye)
{
  std::map<int, bool> hasTLT;

  for (auto teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {
    const int telID = teliter->GetId();
    if (teliter->HasTriggerData()) {
      const fevt::TelescopeTriggerData& telTrig = teliter->GetTriggerData();
      hasTLT[telID] = telTrig.IsTLTAccepted();
    } else
      hasTLT[telID] = false;
  }
  return hasTLT;
}


std::map<int, std::string>
otoa::fd::CalcMirrorTLTLabelMap(const fevt::Eye& eye)
{
  std::map<int, std::string> label;

  for (auto teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {
    const int telID = teliter->GetId();
    if (teliter->HasTriggerData()) {
      const fevt::TelescopeTriggerData& telTrig = teliter->GetTriggerData();
      label[telID] = telTrig.IsTLTAccepted() ? telTrig.GetTLTLabel() : " --- ";
    } else
      label[telID] = " --- ";
  }
  return label;
}


unsigned long long
otoa::fd::CalcMirrorsInEventBitField(const fevt::Eye& eye)
{
  unsigned long long telBits = 0;
  bool telHasData = false;
  for (auto teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {
    telHasData = false;
    for (auto iPixel = teliter->PixelsBegin(ComponentSelector::eHasData);
         iPixel != teliter->PixelsEnd(ComponentSelector::eHasData); ++iPixel) {
      if (!iPixel->HasRecData())
        continue;
      telHasData = true;
      break;
    }
    if (telHasData) {
      const int bit = teliter->GetId() - 1;
      telBits += 1ULL << bit;
    }
  }
  return telBits;
}


unsigned long long
otoa::fd::CalcMirrorsInDAQBitField(const fevt::Eye& eye)
{
  unsigned long long telBits = 0;
  for (auto iTel =  eye.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
       iTel != eye.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {
    const int bit = iTel->GetId() - 1;
    telBits += 1ULL << bit;
  }
  return telBits;
}


std::map<UShort_t, Int_t>
otoa::fd::CalcMirrorTimeOffsetMap(const fevt::Eye& eye)
{
  std::map<UShort_t, Int_t> timeOffset;
  for (auto teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter)
    timeOffset[teliter->GetId()] = teliter->GetTimeOffset();

  return timeOffset;
}


bool
otoa::fd::UsingMieAttenuationDatabase()
{
  // first check if Mie DB data is actually used ...
  utl::Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("Atmosphere");
  utl::Branch mieModel = topB.GetChild("MieModel");
  if (!mieModel)
    return false;
  string mieModelName;
  mieModel.GetData(mieModelName);

  // ... then check if there is data inside
  if (mieModelName == "Super" || mieModelName == "MeasuredDB") {
    atm::MeasuredDBMieModel measuredMieModel;
    return measuredMieModel.HasData();
  }
  return false;
}


bool
otoa::fd::UsingGDASProfileDatabase()
{
  // first check if GDAS data is actually used ...
  utl::Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("Atmosphere");
  utl::Branch profileModel = topB.GetChild("ProfileModel");
  if (!profileModel)
    return false;
  string profileModelName;
  profileModel.GetData(profileModelName);

  // ... then check if there is data inside
  if (profileModelName == "Super" || profileModelName == "GDAS") {
    atm::GDASProfileModel gdasProfileModel;
    return gdasProfileModel.HasData();
  }
  return false;
}


double
otoa::fd::ConvertXToChi(const double x, const fdet::Eye& detEye,
                        const Point& core, const Vector& axis, const Vector& sdp)
{
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

  try {
    const double deltaX = 10*g/cm2;
    InclinedAtmosphericProfile incModel(core, axis, deltaX);

    // use the inclined atmosphere model
    const ProfileResult& distVsSlantDepth = incModel.EvaluateDistanceVsSlantDepth();

    double distance = 0;
    if (distVsSlantDepth.MinX() < x && x < distVsSlantDepth.MaxX())
      distance = -distVsSlantDepth.Y(x);  // DEF: p = core + distance * axis

    const utl::Point& eyePos = detEye.GetPosition();

    const utl::Point& xPoint = core + distance*axis;

    const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
    const Vector vertical(0, 0, 1, eyeCS);
    const Vector horizontal = cross(sdp, vertical);

    double angle = Angle((xPoint-eyePos), horizontal);

    return angle;
  } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    // do nothing
  } catch (OutOfBoundException& e) {
    // do nothing
  }

  return 0;
}


double
otoa::fd::ConvertXToChi(const double x, const fdet::Telescope& detTel,
                        const Point& core, const Vector& axis, const Vector& sdp)
{
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

  try {
    const double deltaX = 10*g/cm2;
    InclinedAtmosphericProfile incModel(core, axis, deltaX);

    // use the inclined atmosphere model
    const ProfileResult& distVsSlantDepth = incModel.EvaluateDistanceVsSlantDepth();
    //const ProfileResult &heightVsSlantDepth = incModel.EvaluateHeightVsSlantDepth();

    double distance = 0;
    if (distVsSlantDepth.MinX() < x && x < distVsSlantDepth.MaxX())
      distance = -distVsSlantDepth.Y(x); // DEF: p = core + distance * axis

    const utl::Point& telPos = detTel.GetPosition();

    const utl::Point& xPoint = core + distance*axis;

    CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(telPos);
    // Note that the following is wrong for telescopes:
    //const CoordinateSystemPtr localCS = detTel.GetTelescopeeCoordinateSystem();
    const Vector vertical(0, 0, 1, localCS);
    const Vector horizontal = cross(sdp, vertical);

    double angle = Angle((xPoint-telPos), horizontal);

    return angle;
  } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    // do nothing
  } catch (OutOfBoundException& e) {
    // do nothing
  }

  return 0;
}


double
otoa::fd::LinearProfileFitChiSquare(const FdRecShower& recShower)
{
  const vector<double>& edep = recShower.GetEnergyDeposit();
  const vector<double>& edep_err = recShower.GetEnergyDepositError();
  const vector<double>& ne = recShower.GetElectrons();
  const vector<double>& ne_err = recShower.GetElectronsError();

  const vector<double>& y = edep.empty() ? ne : edep;
  const vector<double>& ye = edep_err.empty() ? ne_err : edep_err;

  double a1 = 0;
  double a2 = 0;
  double chi2 = 0;

  otoa::LinearFit(recShower.GetDepth(), y, ye, a1, a2, chi2);

  return std::isnan(chi2) ? 0 : chi2;
}


double
otoa::fd::GetCalibrationConstant(const fevt::Pixel& evtpixel)
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Pixel& detPixel = detFD.GetPixel(evtpixel);

  const int telid = detPixel.GetTelescopeId();
  const int eyeid = detPixel.GetEyeId();
  const int pixelid = detPixel.GetId();

  // Get the calibration constant for this channel
  const TabulatedFunction& calibration = detPixel.GetEndToEndCalibrationConstant();

  if (!calibration.GetNPoints()) {
    ostringstream err;
    err << "Could not find calibration for eye " << eyeid
        << " tel " << telid << " pix " << pixelid;
    ERROR(err);
    return false;
  }

  // check calibration status
  if (detPixel.GetStatus() != fdet::Pixel::eGood) {
    ostringstream warn;
    warn << " bad calibration data for eye " << eyeid
         << " tel " << telid << " pix " << pixelid
         << " status = " << detPixel.GetStatus();
    WARNING(warn);
    return false;
  }

  const double bestCalibConst =
    detPixel.GetEndToEndCalibrationAtReferenceWavelength();

  return bestCalibConst;
}