FD2ADST.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "FD2ADST.h"

#include "Offline2ADST.h"
#include "Config.h"
#include "FOVCalculator.h"
#include "FD2ADSTUtil.h"
#include "ConversionUtil.h"
#include "ErrorPropagation.h"
#include <FdPCGFData.h>

#include <adst/RecEvent.h>
#include <adst/FdApertureLight.h>
#include <adst/FdRecLevel.h>

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

#include <utl/AugerUnits.h>
#include <utl/CorrelationMatrix.h>
#include <utl/NumericalErrorPropagation.h>
#include <utl/ModifiedJulianDate.h>
#include <utl/PhysicalConstants.h>
#include <utl/Transformation.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/UTMPoint.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 <atm/Atmosphere.h>
#include <atm/AttenuationResult.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/LidarCloudDBModel.h>
#include <atm/ProfileResult.h>
#include <atm/GOESDB.h>

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

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/SortCriteria.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/GaisserHillas2Parameter.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/GaisserHillas6Parameter.h>
#include <evt/MultipleGaisserHillasParameters.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 <AugerEvent.h>
#include <EyeEvent.hh>
#include <EyeEventHeader.hh>

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

using namespace utl;
using namespace fwk;
using namespace std;
using namespace otoa;
using namespace otoa::fd;
using namespace fevt;
using sevt::SEvent;
using sevt::StationRecData;


FD2ADST::FD2ADST(const otoa::Offline2ADST& parent, otoa::Status& status) :
  fParent(parent),
  fStatus(status)
{ }


void
FD2ADST::Convert(const evt::Event& event, RecEvent& recEvent)
{
  for (fevt::FEvent::ConstEyeIterator iEye = event.GetFEvent().EyesBegin(ComponentSelector::eInDAQ),
       end = event.GetFEvent().EyesEnd(ComponentSelector::eInDAQ); iEye != end; ++iEye) {

    if (iEye->HasHeader() ||  // this requires eHasData and an eye trigger (event)
        !Config().DropUntriggeredMCProfiles()) {

      FDEvent fdEvent;
      fdEvent.SetRecLevel(eNoFdEvent);

      FillEyeHeader(*iEye, fdEvent);

      // ---- only for eyes without reconstruction ------
      bool haveSimTelData = false;
      if (iEye->HasRecData() ||
          !Config().DropUntriggeredMCProfiles()) {
        FillEyeSim(event, fdEvent);
        haveSimTelData = FillTelSimData(*iEye, fdEvent);
      }

      // ---- only for eyes without reconstruction ------
      if (iEye->HasRecData()) {

        if (otoa::fd::UsingMieAttenuationDatabase())
          recEvent.GetDetector().SetHasMieDatabase();
        if (otoa::fd::UsingGDASProfileDatabase())
          recEvent.GetDetector().SetHasGDASDatabase();

        FillEye(*iEye, fdEvent); // sets FdRecLevel

        FillTelRecData(*iEye, fdEvent);
        FillRecPixel(*iEye, fdEvent);
        FillFdRecStations(event, *iEye, fdEvent);

        if (Config().StoreCloudCameraData() > 0)
          FillCloudCameraDataBrief(*iEye, fdEvent);

        if (Config().StoreGOESData() > 0)
          FillCloudsBetweenEyeAndShower(*iEye, fdEvent);

        if (event.HasSEvent()) {
          FillSDEye(event,*iEye,fdEvent);
          FillHybridStations(event, *iEye, fdEvent.GetFdRecGeometry());
        }
      }

      if (iEye->HasHeader() || haveSimTelData)
        recEvent.AddEye(fdEvent);
    }

  }

  // Do the FD field of view calculations for all eyes
  if (Config().GetFieldOfViewOption()) {
    FOVCalculator fovCalc;
    if (Config().GetFieldOfViewOption() == 2)
      fovCalc.SetUseBGLoop(true);
    else
      fovCalc.SetUseBGLoop(false);
    fovCalc.FillFOVVariables(event.GetFEvent(), recEvent);
  }
}


void
FD2ADST::FillEyeHeader(const fevt::Eye& eye, FDEvent& theFd)
{
  // -------- eye header variables ---------
  const int eyeId = eye.GetId();
  theFd.SetEyeId((UShort_t) eyeId);

  if (eye.HasHeader()) {  // this requires a triggered eye and eHasData
    IncreaseFdRecLevel(theFd, eHasFdEvent);
    const fevt::EyeHeader& header = eye.GetHeader();
    theFd.SetRunId((UInt_t) header.GetRunNumber());
    theFd.SetEventId((UInt_t) header.GetEventNumber());
    const utl::TimeStamp& etime = header.GetTimeStamp();
    theFd.SetGPSSecond((UInt_t) etime.GetGPSSecond());
    theFd.SetGPSNanoSecond((UInt_t) etime.GetGPSNanoSecond());

    if (header.HasBadTimeCorrection()) {
      if (header.IsOfflineTimeCorrected())
        theFd.SetTimeCorrectionStatus(eOfflineCorrected);
      else
        theFd.SetTimeCorrectionStatus(eBadCorrection);
    } else
      theFd.SetTimeCorrectionStatus(eGoodCorrection);

    theFd.SetYYMMDD(TimeStamp2YYMMDD(etime));
    theFd.SetHHMMSS(TimeStamp2HHMMSS(etime));
    theFd.SetMoonCycle(TimeStamp2MoonCycle(etime));

  } else {
    // empty for now since all of the above are initialized to 0 by FDEvent
  }

  if (eye.HasTriggerData()) {
    const fevt::EyeTriggerData& eyeTrig = eye.GetTriggerData();
    theFd.SetT3Data(eyeTrig.GetT3SDPTheta(),
                    eyeTrig.GetT3SDPPhi(),
                    eyeTrig.GetT3NPixels(),
                    eyeTrig.GetT3AzimuthAtGround(),
                    eyeTrig.GetT3Time().GetGPSSecond(),
                    eyeTrig.GetT3Time().GetGPSNanoSecond());
  } else {
    // reducing printout
    //WARNING("FillEyeHeader --> This event has no T3-data. Not filling T3 data.");
  }

  theFd.SetTLT(CalcMirrorTLTMap(eye));
  theFd.SetTLTLabel(CalcMirrorTLTLabelMap(eye));

  theFd.SetMirrorTimeOffsets(CalcMirrorTimeOffsetMap(eye));
  const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(eyeId);
  const int nTel = detEye.GetLastTelescopeId() - detEye.GetFirstTelescopeId() + 1;

  unsigned long long evStatus = CalcMirrorsInEventBitField(eye);
  TBits evBits(7);
  for (int thisTelId = 1; thisTelId <= nTel; ++thisTelId) {
    if ((evStatus/(1ULL << (thisTelId-1))) % 2 > 0)
      evBits.SetBitNumber(thisTelId, true);
  }
  theFd.SetMirrorsInEvent(evBits);

  // ------------- Fill event class and type ----
  string evClass("unknown"), evType("unknown");
  if (eye.HasTriggerData())
    evClass = eye.GetTriggerData().GetT3Class();
  if (Status().IsMCEvent())
    evType = " Simulated Shower ";

  if (Status().GetEvent().HasRawEvent()) {
    const AugerEvent& rawEventRef = Status().GetEvent().GetRawEvent();
    AugerEvent* const rawEvent = (AugerEvent*)&rawEventRef;

    AugerEvent::EyeIterator iRawEye;
    for (iRawEye = rawEvent->EyesBegin(); iRawEye != rawEvent->EyesEnd(); ++iRawEye) {
      if (iRawEye->GetEventHeader()->GetEyeNo() == (unsigned int)eyeId)
        break;
    }

    if (iRawEye != rawEvent->EyesEnd()) {
      TEyeEventHeader& eyeheader = *iRawEye->GetEventHeader();
      evClass = string(eyeheader.GetVerboseEventClass(eyeheader.GetEventClass()));
      evType = string(eyeheader.GetVerboseEventType(eyeheader.GetEventType()));
    }
  }

  theFd.SetEventClassAndType(evClass, evType);
}


void
FD2ADST::FillEye(const fevt::Eye& eye, FDEvent& theFd)
{
  // Note: This isn't even called if the eye doesn't have RecData...

  const int eyeId = eye.GetId();
  const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(eyeId);
  const fevt::EyeRecData& eyeRec = eye.GetRecData();

  if (eyeRec.GetTimeFitNDof() > 0) {
    const evt::ShowerFRecData& frecData = eyeRec.GetFRecShower();
    // recalculate slant depth table that might have been
    // screwed up during error propagation or for multi-eye events
    INFO(" re-calculating slant tables...");
    const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
    try {
      atmo.InitSlantProfileModel(frecData.GetCorePosition(), frecData.GetAxis(), 10*g/cm2);
    } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException&) {
    }
  }

  if (!eye.HasHeader())  // this requires a triggered eye and eHasData
    return;

  IncreaseFdRecLevel(theFd, eHasFdEvent); // this is in principle already done...

  if (!eye.HasHeader())  // this requires a triggered eye and eHasData
    return;

  IncreaseFdRecLevel(theFd, eHasFdEvent); // this is in principle already done...

  // -----  Figure out triggered/pulsed pixel numbers ----------
  // Note: This is sort of duplicated from FillRecPixel. The reason
  //       is that we don't want to fill the BIG recpixel stuff just
  //       yet since we don't know whether we will use the FDEvent at all.
  //       Instead, we duplicate the stuff we need for the FdRecLevel here.
  unsigned int nTriggeredPixels = 0;
  unsigned int nPulsedPixels = 0;

  for (fevt::Eye::ConstTelescopeIterator teliter = eye.TelescopesBegin(ComponentSelector::eHasData),
       end = eye.TelescopesEnd(ComponentSelector::eHasData); teliter != end; ++teliter) {

    for (fevt::Telescope::ConstPixelIterator iPixel = teliter->PixelsBegin(ComponentSelector::eHasData),
         end = teliter->PixelsEnd(ComponentSelector::eHasData); iPixel != end; ++iPixel) {
      if (iPixel->HasTriggerData() && iPixel->GetTriggerData().IsTriggered())
        ++nTriggeredPixels;
    }

    for (fevt::EyeRecData::ConstPixelIterator pixiter = eye.GetRecData().PulsedPixelsBegin(),
         end = eye.GetRecData().PulsedPixelsEnd(); pixiter != end; ++pixiter)
      ++nPulsedPixels;
  }

  if (nTriggeredPixels <= 0)
    return;

  IncreaseFdRecLevel(theFd, eHasTriggeredPixels);

  if (nPulsedPixels <= 0)
    return;

  IncreaseFdRecLevel(theFd, eHasReconstructedPixels);
  //fnRecPixels = nPulsedPixels;

  // ----- End of duplicated stuff

  // -------- geometry variables ---------
  if (eyeRec.GetSDPFitNDof() <= 0)  // if no SDP fit
    return;

  IncreaseFdRecLevel(theFd, eHasSDP);

  FdRecGeometry& theGeometry = theFd.GetFdRecGeometry();
  FillSDP(detEye, eyeRec, theGeometry);

  if (eyeRec.GetTimeFitNDof() <= 0)  // if no time fit
    return;

  IncreaseFdRecLevel(theFd, eHasAxis);

  FillTimeFit(detEye, eyeRec, theGeometry);

  FillPCGFit(detEye, eyeRec, theGeometry);

//#warning checking hybrid stations was done here???
  // -------- SD related variables ---------
  //FillHybridStations(event, eye, theGeometry);

  // -------- high level geometry variables ---------
  FillFdCoreAxis(detEye, eye, theFd);

  // -------- shower profile variables ---------
  if (!FillEyeApertureLight(detEye, eye, theFd))
    return;

  IncreaseFdRecLevel(theFd, eHasAperturePhotons);

  if (!FillFdProfile(detEye, eye, theFd))
    return;

  // Note: RecLevel increased by FillFdProfile

  if (!FillGaisserHillas(detEye, eye, theFd))
    return;

  IncreaseFdRecLevel(theFd, eHasGHParameters);

  FdRecShower& fdRecShower = theFd.GetFdRecShower();
  fdRecShower.SetLinearProfileFitChi2(otoa::fd::LinearProfileFitChiSquare(fdRecShower));

  FillAtmosphericProfileVars(detEye, eye, theFd);

  if (!FillEnergy(detEye, eye, theFd))
    return;

  IncreaseFdRecLevel(theFd, eHasEnergy);
}


void
FD2ADST::FillSDP(const fdet::Eye& detEye, const fevt::EyeRecData& eyeRec, FdRecGeometry& theGeometry)
{
  utl::CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const utl::Vector& sdp = eyeRec.GetSDP();

  theGeometry.SetSDP(sdp.GetTheta(eyeCS), sdp.GetPhi(eyeCS));
  theGeometry.SetSDPThetaError(eyeRec.GetSDPThetaError());
  theGeometry.SetSDPPhiError(eyeRec.GetSDPPhiError());
  theGeometry.SetSDPChi2(eyeRec.GetSDPFitChiSquare());
  theGeometry.SetSDPNdF((UInt_t) eyeRec.GetSDPFitNDof());
  theGeometry.SetSDPThetaPhiCorrelation(eyeRec.GetSDPCorrThetaPhi());
}


void
FD2ADST::FillTimeFit(const fdet::Eye& /*detEye*/, const fevt::EyeRecData& eyeRec, FdRecGeometry& theGeometry)
{
  theGeometry.SetAxis(eyeRec.GetTZero(), eyeRec.GetRp(), eyeRec.GetChiZero());
  theGeometry.SetT0Error(eyeRec.GetTZeroError());
  theGeometry.SetRpError(eyeRec.GetRpError());
  theGeometry.SetChi0Error(eyeRec.GetChiZeroError());
  theGeometry.SetChi0T0Correlation(eyeRec.GetChi0TZeroCorrelation());
  theGeometry.SetChi0RpCorrelation(eyeRec.GetRpChi0Correlation());
  theGeometry.SetRpT0Correlation(eyeRec.GetRpTZeroCorrelation());
  theGeometry.SetTimeFitChi2(eyeRec.GetTimeFitChiSquare());
  theGeometry.SetTimeFitNdof(eyeRec.GetTimeFitNDof());
  theGeometry.SetNorthEastCorrelation(eyeRec.GetNorthEastCorrelation());
  theGeometry.SetNorthThetaCorrelation(eyeRec.GetNorthThetaCorrelation());
  theGeometry.SetNorthPhiCorrelation(eyeRec.GetNorthPhiCorrelation());
  theGeometry.SetNorthTCoreCorrelation(eyeRec.GetNorthTCoreCorrelation());
  theGeometry.SetEastThetaCorrelation(eyeRec.GetEastThetaCorrelation());
  theGeometry.SetEastPhiCorrelation(eyeRec.GetEastPhiCorrelation());
  theGeometry.SetEastTCoreCorrelation(eyeRec.GetEastTCoreCorrelation());
  theGeometry.SetThetaPhiCorrelation(eyeRec.GetThetaPhiCorrelation());
  theGeometry.SetThetaTCoreCorrelation(eyeRec.GetThetaTCoreCorrelation());
  theGeometry.SetPhiTCoreCorrelation(eyeRec.GetPhiTCoreCorrelation());
  theGeometry.SetAxisFitChi2(eyeRec.GetAxisFitChiSquare());
  theGeometry.SetAxisFitNdof(eyeRec.GetAxisFitNDof());
  theGeometry.SetNTimeFitPixels(eyeRec.GetNTimeFitPixels());

  // set the geometry rec level
  if (eyeRec.GetAxisFitChiSquare() <= 0) {
    if (eyeRec.GetFRecShower().GetStationIds().size() > 0)
      theGeometry.SetGeomRecLevel(eHybridGeometry);
    else
      theGeometry.SetGeomRecLevel(eMonoGeometry);
  } else {
    if (eyeRec.GetFRecShower().GetStationIds().size() > 0)
      theGeometry.SetGeomRecLevel(eAxisHybridGeometry);
    else
      theGeometry.SetGeomRecLevel(eAxisMonoGeometry);
  }

  // re-calculate the FD-only time fit Chi2 (for Bruce/Jose cuts)
  // clean up timefit pixels if in second iteration
  {
    double t_0 = eyeRec.GetTZero();
    double r_p = eyeRec.GetRp();
    double chi_0 = eyeRec.GetChiZero();
    double chisq = 0;
    int ndof = 0;
    for (auto iPixel = eyeRec.TimeFitPixelsBegin(); iPixel != eyeRec.TimeFitPixelsEnd(); ++iPixel) {
      const PixelRecData& pixrec = iPixel->GetRecData();
      const double t_i = pixrec.GetT_i().GetNanoSecond();
      const double t_i_Err = pixrec.GetT_iError().GetInterval();
      const double chi_i = pixrec.GetChi_i();

      if (t_i_Err > 0) {
        const double t_expected = t_0 + r_p/utl::kSpeedOfLight * std::tan(0.5*(chi_0-chi_i));
        const double tmp = t_i - t_expected;
        chisq += tmp*tmp / (t_i_Err*t_i_Err);
        ++ndof;
      } else
        cout << " zero error in pixel! " << t_i << endl;
    }

    theGeometry.SetTimeFitFDChi2(chisq);
  } // end scope of calculating FD only time fit chi2
}


void
FD2ADST::FillPCGFit(const fdet::Eye& /*dEye*/, const fevt::EyeRecData& eyeRec, FdRecGeometry& theGeometry)
{
  for (auto it = eyeRec.GetPCGF().begin(), end = eyeRec.GetPCGF().end(); it != end; ++it) {
    FdPCGFData pcgf;
    pcgf.SetChi0(it->GetChi0());
    pcgf.SetRp(it->GetRp());
    pcgf.SetT0(it->GetT0());
    switch (it->GetStatus()) {
    case PCGFData::eUnknown:
      pcgf.SetStatus(ePCGFUnknown);
      break;
    case PCGFData::eScan:
      pcgf.SetStatus(ePCGFScan);
      break;
    case PCGFData::eFit:
      pcgf.SetStatus(ePCGFFit);
      break;
    case PCGFData::ePreFit:
      pcgf.SetStatus(ePCGFInitial);
      break;
    case PCGFData::eFinal:
      pcgf.SetStatus(ePCGFFinal);
      break;
    }
    for (auto ichi2 = it->GetChi2().begin(), ichi2end = it->GetChi2().end(); ichi2 != ichi2end; ++ichi2)
      pcgf.AddChi2(ichi2->first, ichi2->second);

    for (auto ichi2 = it->GetNDof().begin(), ichi2end = it->GetNDof().end(); ichi2 != ichi2end; ++ichi2)
      pcgf.AddNDof(ichi2->first, ichi2->second);

    theGeometry.AddPCGFData(pcgf);
  }
}


void
FD2ADST::FillFdCoreAxis(const fdet::Eye& detEye, const fevt::Eye& eye, FDEvent& theFd)
{
  const evt::ShowerFRecData& showerFRec = eye.GetRecData().GetFRecShower();
  FdRecShower& theShower = theFd.GetFdRecShower();

  // -- core position --

  const utl::Point& core = showerFRec.GetCorePosition();
  const utl::Vector& axis = showerFRec.GetAxis();

  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
  const CoordinateSystemPtr coreLocalCS = fwk::LocalCoordinateSystem::Create(core);
  const CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();

  CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  Point eyePos = detEye.GetPosition();

  // core (first of all do the utm-exception check !!!!)
  double northing = 0;
  double easting = 0;
  double altitude = 0;
  try {
    utl::UTMPoint UTM(core, wgs84);
    northing = UTM.GetNorthing();
    easting = UTM.GetEasting();
  } catch (UTMPoint::UTMZoneException& zone_e) {
    ERROR ("UTMZoneException: fd rec shower invalid");
    return;
  } catch (UTMPoint::UTMException& e) {
    ERROR ("UTMException: fd rec shower invalid");
    return;
  }

  // check if altitude can be calculated -> NEED hottest station
  if (Status().GetEvent().HasSEvent()) {  // don't tell my mom I used the sneaky global Event
    const sevt::SEvent& sevent = Status().GetEvent().GetSEvent();

    int hottestId = 0;
    double hottestSignal = 0;
    bool first = true;
    const vector<unsigned short int>& hybridId = showerFRec.GetStationIds();
    for (auto iStation = hybridId.begin(), end = hybridId.end(); iStation != end; ++iStation) {
      if (!sevent.HasStation(*iStation))
        continue;

      const sevt::Station& station = sevent.GetStation(*iStation);
      if (!station.HasRecData())
        continue;

      const double signal = station.GetRecData().GetTotalSignal();

      if (first) {
        first = false;
        hottestSignal = signal;
        hottestId = *iStation;
      } else if (hottestSignal < signal) {
        hottestSignal = signal;
        hottestId = *iStation;
      }
    }

    // if there is a (hottest) hybrid station:
    if (!first) {
      utl::CoordinateSystemPtr stationCS =
        det::Detector::GetInstance().GetSDetector().GetStation(hottestId).GetLocalCoordinateSystem();

      // project FD core to ground
      const utl::Point coreStationCS = core - axis * core.GetZ(stationCS)/axis.GetZ(stationCS);
      northing = 0;
      easting = 0;
      altitude = 0;
      try {
        utl::UTMPoint UTM(coreStationCS, ReferenceEllipsoid::GetWGS84());
        northing = UTM.GetNorthing();
        easting = UTM.GetEasting();
        altitude = UTM.GetHeight();
      } catch (UTMPoint::UTMZoneException& zone_e) {
        ERROR("UTMZoneException: fd rec shower invalid");
      } catch (UTMPoint::UTMException& e) {
        ERROR("UTMException: fd rec shower invalid");
      }
    }
  }

  theShower.SetCoreUTMCS(TVector3(easting, northing, altitude));

  theShower.SetCoreSiteCS(ToTVector3(core, referenceCS, meter));

  theShower.SetCoreTime(showerFRec.GetCoreTime().GetGPSSecond(),
                        showerFRec.GetCoreTime().GetGPSNanoSecond());

  TVector3 axisCore(0,0,1);
  axisCore.SetTheta(axis.GetTheta(coreLocalCS));
  axisCore.SetPhi(axis.GetPhi(coreLocalCS));
  theShower.SetAxisCoreCS(axisCore);

  TVector3 axisSite(0,0,1);
  axisSite.SetTheta(axis.GetTheta(referenceCS));
  axisSite.SetPhi(axis.GetPhi(referenceCS));
  theShower.SetAxisSiteCS(axisSite);

  FillGeometricalUncertainties(detEye, eye, theFd);

  FillCelestialCoordinates(theFd.GetFdRecShower());
}


void
FD2ADST::FillGeometricalUncertainties(const fdet::Eye& detEye, const fevt::Eye& eye, FDEvent& theFd)
{
  // partly from FdGeomErrorsEstimator::FindErrors_B

  FdRecShower& theShower = theFd.GetFdRecShower();
  const fevt::EyeRecData& eyeRec = eye.GetRecData();
  CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();

  //const double t0 = eyeRec.GetTZero();
  const double chi0 = eyeRec.GetChiZero();
  const double rp = eyeRec.GetRp();

  //const double t0_err = eyeRec.GetTZeroError();
  const double chi0_err = eyeRec.GetChiZeroError();
  const double rp_err = eyeRec.GetRpError();

  const utl::AxialVector& sdp = eyeRec.GetSDP();
  const double theta_sdp = sdp.GetTheta(eyeCS);
  const double phi_sdp = sdp.GetPhi(eyeCS);

  const double theta_sdp_err = eyeRec.GetSDPThetaError();
  const double phi_sdp_err = eyeRec.GetSDPPhiError();

  // Retrieve correlation parameters
  const double corr_tp = eyeRec.GetSDPCorrThetaPhi();
  //const double corr_chi0t0 = eyeRec.GetChi0TZeroCorrelation();
  const double corr_chi0rp = eyeRec.GetRpChi0Correlation();
  //const double corr_rpt0 = eyeRec.GetRpTZeroCorrelation();

  // core error
  // INPUT-vector: Rp, Chi0, sdp-theta, sdp-phi
  std::vector<utl::Parameter> coreParameters;
  coreParameters.push_back(Parameter(rp, rp_err));
  coreParameters.push_back(Parameter(chi0, chi0_err));
  coreParameters.push_back(Parameter(theta_sdp, theta_sdp_err));
  coreParameters.push_back(Parameter(phi_sdp, phi_sdp_err));
  utl::CorrelationMatrix coreCorrelations(4);
  coreCorrelations.Set(0, 1, corr_chi0rp);
  coreCorrelations.Set(0, 2, 0);
  coreCorrelations.Set(0, 3, 0);
  coreCorrelations.Set(1, 2, 0);
  coreCorrelations.Set(1, 3, 0);
  coreCorrelations.Set(2, 3, corr_tp);
  otoa::err::CalculateFdCorePosition tmp(eye); // Temporary is necessary to work around 2001 compiler bug in Karlsruhe's ancient cluster
  utl::NumericalErrorPropagation coreErrors(tmp, coreParameters, coreCorrelations);

  theShower.SetCoreNorthingError(coreErrors.GetPropagatedErrors()[1]);
  theShower.SetCoreEastingError(coreErrors.GetPropagatedErrors()[0]);
  theShower.SetCoreNorthingEastingCorrelation(coreErrors.GetPropagatedCorrelations().Get(0,1));

  // axis error
  // INPUT-vector: Chi0, SDP-theta, SDP-phi
  // OUTPUT-vector: theta, phi in core-CS
  std::vector<utl::Parameter> axisParameter;
  axisParameter.push_back(Parameter(chi0, chi0_err));
  axisParameter.push_back(Parameter(theta_sdp, theta_sdp_err));
  axisParameter.push_back(Parameter(phi_sdp, phi_sdp_err));
  utl::CorrelationMatrix axisCorrelations(3);
  axisCorrelations.Set(0, 1, 0);
  axisCorrelations.Set(0, 2, 0);
  axisCorrelations.Set(1, 2, corr_tp);
  otoa::err::CalculateFdArrivalDirection tmp2(eye);  // Temporary is necessary to work around 2001 compiler bug in Karlsruhe's ancient cluster
  utl::NumericalErrorPropagation axisErrors(tmp2, axisParameter, axisCorrelations);
  //const vector<double>& out_e2 = axisErrors.GetPropagatedErrors();
  //const CorrelationMatrix& out_c2 = axisErrors.GetPropagatedCorrelations();

  theShower.SetZenithError(axisErrors.GetPropagatedErrors()[0]);
  theShower.SetAzimuthError(axisErrors.GetPropagatedErrors()[1]);
  theShower.SetZenithAzimuthCorrelation(axisErrors.GetPropagatedCorrelations().Get(0,1));
}


bool
FD2ADST::FillEyeApertureLight(const fdet::Eye& /*dEye*/, const fevt::Eye& eye, FDEvent& theFd)
{
  bool hasApertureLight = false;

  // -------- first check for existence of tel. reco light ------
  // Note: This is for the rec level!
  for (fevt::Eye::ConstTelescopeIterator teliter = eye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != eye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {
    if (!teliter->HasRecData())
      continue;
    const TelescopeRecData& oRecData = teliter->GetRecData();
    if (!oRecData.HasLightFlux(fevt::FdConstants::eTotal))
      continue;
    hasApertureLight = true;
    break;
  }

  const fevt::EyeRecData& eyeRec = eye.GetRecData();

  // ##### Getting Light at Aperture #########
  if (eyeRec.HasLightFlux(fevt::FdConstants::eTotal)) {
    hasApertureLight = true;

    const TabulatedFunctionErrors& diaPhotonTrace =
      eyeRec.GetLightFlux(fevt::FdConstants::eTotal);
    const double binSize = 100.*ns;
    const unsigned int time_size = diaPhotonTrace.GetNPoints();

    vector<double> dia_time(time_size);
    vector<double> dia_time_err(time_size);
    vector<double> dia_photons(time_size);
    vector<double> dia_photons_err(time_size);

    for (unsigned int i = 0; i < time_size; ++i) {
      dia_time[i] = diaPhotonTrace [i].X()/binSize;
      dia_time_err[i] = diaPhotonTrace.GetXErr(i)/binSize;

      if (dia_time_err[i] == 0)
        dia_time_err[i] = 1;

      dia_photons[i] = diaPhotonTrace[i].Y() / (2.*dia_time_err[i]);
      dia_photons_err[i] = diaPhotonTrace.GetYErr(i)/  (2.*dia_time_err[i]);
    }

    FdRecApertureLight& theProfile = theFd.GetFdRecApertureLight();
    theProfile.SetTimeBinning(binSize/ns);

//#warning Zeta is currently eye-global. That may become wrong in the face of telescope light at aperture.
    theProfile.SetZeta(eyeRec.GetZeta());
    theProfile.SetTime(dia_time, dia_time_err);
    theProfile.SetTotalLightAtAperture(dia_photons, dia_photons_err);
  }

  return hasApertureLight;
}


void
FD2ADST::FillCloudsBetweenEyeAndShower(const fevt::Eye& evtEye, FDEvent& theFd)
  const
{
  if (!evtEye.HasRecData() || !evtEye.GetRecData().HasFRecShower())
    return;

  try {
    const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
    const atm::GOESDB& goesDB = atmo.GetGOESDB();
    const evt::ShowerFRecData& shower = evtEye.GetRecData().GetFRecShower();

    const fdet::Eye& detEye =
      det::Detector::GetInstance().GetFDetector().GetEye(evtEye.GetId());

    const TimeStamp& eyeGPSTime = evtEye.GetHeader().GetTimeStamp();
    const TimeStamp eyeTraceStartTime = eyeGPSTime -
      TimeInterval(detEye.TelescopesBegin()->GetCamera().GetFADCTraceLength() *
                   detEye.TelescopesBegin()->GetCamera().GetFADCBinSize());

    // first loop over all TelRecData for multi-tel reconstruction results

    double maxFraction = -1;
    bool maxFractionFilled = false;

    for (auto iTel = evtEye.TelescopesBegin(fevt::ComponentSelector::eHasData);
         iTel != evtEye.TelescopesEnd(fevt::ComponentSelector::eHasData); ++iTel) {

      if (!iTel->HasRecData() ||
          !iTel->GetRecData().HasLightFlux(fevt::FdConstants::eTotal))
        continue;

      const Point& telescopePosition =
        detEye.GetTelescope(iTel->GetId()).GetPosition();

      const TabulatedFunctionErrors& photonTrace =
        iTel->GetRecData().GetLightFlux(fevt::FdConstants::eTotal);
      for (unsigned int timeBin = 0; timeBin < photonTrace.GetNPoints(); ++timeBin) {
        const TimeStamp timeAtTelescope = eyeTraceStartTime + photonTrace.GetX(timeBin);
        const Point pointAtShower = shower.CalculatePointOnShower(timeAtTelescope, telescopePosition);
        const double thisFrac =
          goesDB.GetMaximumCloudProbability(pointAtShower, telescopePosition);
        if (!maxFractionFilled || thisFrac > maxFraction)
          maxFraction = thisFrac;
        maxFractionFilled = true;
      }
    }

    if (maxFractionFilled) {
      theFd.GetFdRecShower().SetMaxCloudFractionBetweenEyeAndShower(maxFraction);
      return;
    }

    // if nothing filled yet, try the eye reconstruction
    const Point& eyePosition = detEye.GetPosition();
    if (evtEye.GetRecData().HasLightFlux(fevt::FdConstants::eTotal)) {
      const TabulatedFunctionErrors& photonTrace =
        evtEye.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);
      for (unsigned int timeBin = 0; timeBin < photonTrace.GetNPoints(); ++timeBin) {
        const TimeStamp timeAtEye = eyeTraceStartTime + photonTrace.GetX(timeBin);
        const Point pointAtShower = shower.CalculatePointOnShower(timeAtEye, eyePosition);
        const double thisFrac =
          goesDB.GetMaximumCloudProbability(pointAtShower, eyePosition);
        if (!maxFractionFilled || thisFrac > maxFraction)
          maxFraction = thisFrac;
        maxFractionFilled = true;
      }
    }

    theFd.GetFdRecShower().SetMaxCloudFractionBetweenEyeAndShower(maxFraction);

  } catch (utl::AugerException& augerException) {
    WARNING(augerException.GetExceptionName() + ", " + augerException.GetMessage());
  }
}


void
FD2ADST::FillCloudCameraDataBrief(const fevt::Eye& iEye, FDEvent& theFd)
{
  float maxCloudFraction = -1;
  float maxCloudFractionLidar = -1;
  float maxCloudFractionHighElevation = -1;
  float maxCloudFractionLidarHighElevation = -1;
  const fdet::FDetector& theFDet = det::Detector::GetInstance().GetFDetector();
  const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();

  for (auto iTel = iEye.TelescopesBegin(fevt::ComponentSelector::eHasData);
       iTel != iEye.TelescopesEnd(fevt::ComponentSelector::eHasData); ++iTel) {

    const unsigned int physicalEyeId =
      theFDet.GetTelescope(*iTel).GetParentPhysicalEyeId();

    utl::CoordinateSystemPtr eyeCS =
      theFDet.GetEye(physicalEyeId).GetEyeCoordinateSystem();

    bool eyeHasCloudData = true;

    if (!iTel->HasRecData())
      continue;

    if (!iEye.GetRecData().HasFRecShower())
      continue;

    const evt::ShowerFRecData& shower = iEye.GetRecData().GetFRecShower();
    const Vector& axis = shower.GetAxis();
    const Point& core = shower.GetCorePosition();

    if (!iEye.HasHeader())
      continue;

    const TimeStamp& coreTime = shower.GetCoreTime();
    const TimeStamp& eyeGPSTime = iEye.GetHeader().GetTimeStamp();
//#warning hardcoded trace time length!
    const TimeStamp eyeTriggerTime = eyeGPSTime-TimeInterval(1000*100*utl::ns);

    // for distance from shower core to emission point (aDist)
    // law of cosines: bDist^2 = aDist^2 + cDist^2 - 2*aDist*cDist*cos(beta)
    // cDist = distance(telescope, core)
    // bDist-aDist = (telTime-coreTime)/c = delta (see below)
    // telTime= TOF emission point to telescope
    const Vector coreTelescopeVec = theFDet.GetTelescope(*iTel).GetPosition() - core;
    const double cDist = coreTelescopeVec.GetMag();
    const double cosBeta = CosAngle(axis, coreTelescopeVec);

    if (!iTel->GetRecData().HasLightFlux(fevt::FdConstants::eTotal))
      continue;

    const TabulatedFunctionErrors& photonTrace =
      iTel->GetRecData().GetLightFlux(fevt::FdConstants::eTotal);

    //pixels used for aperture light per time bin (returns at least empty vec)
    const vector<vector<unsigned int> >& aperturePixels =
      iTel->GetRecData().GetPixelsInZetaOverTime();

    if (aperturePixels.size() != photonTrace.GetNPoints()) {
      // this should never happen
      ostringstream errMsg;
      errMsg << "inconsistent trace/zetaPixels: N(trace) = "
             << photonTrace.GetNPoints() << " N(zetaPixels)="
             << aperturePixels.size() << endl;
      ERROR(errMsg.str());
      throw NonExistentComponentException(errMsg.str());
    }

    //loop over all time bins
    for (unsigned int timeBin = 0; timeBin < aperturePixels.size(); timeBin++) {

      if (!eyeHasCloudData)
        break;

      const double t = photonTrace.GetX(timeBin);

      TimeStamp TimeAtTelescope = eyeTriggerTime + t;
      const double delta =
        (TimeAtTelescope-coreTime).GetInterval()*kSpeedOfLight;
      const double aDist =
        -(cDist*cDist-delta*delta)/(2*(delta+cDist*cosBeta));
      const Point pointOnShower = core-aDist*axis;

      //loop over all pixel in time bin
      for (auto zetaIter = aperturePixels[timeBin].begin();
           zetaIter != aperturePixels[timeBin].end(); ++zetaIter) {

        if (!eyeHasCloudData)
          break;

        try {
          const fevt::Pixel& evtPixel = iTel->GetPixel(*zetaIter);
          const fdet::Pixel& detectorPix = theFDet.GetPixel(evtPixel);
          const double pixelElevation =
            kPiOnTwo - detectorPix.GetDirection().GetTheta(eyeCS);
          const double highElevation = pixelElevation > 5.5*degree;

          //goto next eye and fill with "unknown"=-100
          if (!detectorPix.HasCloudFraction()) {
            eyeHasCloudData = false;
            break;
          }

          //now let us compare with max cloud coverage
          float cloudFraction = detectorPix.GetCloudFraction();
          if (cloudFraction > maxCloudFraction)
            maxCloudFraction = cloudFraction;

          if (highElevation && cloudFraction > maxCloudFractionHighElevation)
            maxCloudFractionHighElevation = cloudFraction;

          try {
            cloudFraction =
              atmo.EvaluateCloudCoverage(detectorPix, pointOnShower).GetCoverage();
            if (cloudFraction > maxCloudFractionLidar)
              maxCloudFractionLidar = cloudFraction;
            if (highElevation && cloudFraction > maxCloudFractionLidarHighElevation)
              maxCloudFractionLidarHighElevation = cloudFraction;
          } catch (NonExistentComponentException&) {
            //ignore
          } catch (NoDataForModelException&) {
            //ignore
          }
        } catch (NonExistentComponentException&) {
          const fevt::Pixel& evtPixel = iTel->GetPixel(*zetaIter);
          ostringstream errMsg;
          errMsg << "could not resolve fevt::Pixel to fdet::Pixel, "
                 << "eye " << evtPixel.GetEyeId() << ", tel " << evtPixel.GetTelescopeId()
                 << ", pix " << evtPixel.GetId();
          ERROR(errMsg);
          FdRecApertureLight& theApLight = theFd.GetFdRecApertureLight();
          theApLight.SetMaxCloudCameraFractionInZeta(-1.);
          theApLight.SetMaxCloudCameraFractionWithLidarInZeta(-1.);
          theApLight.SetMaxCloudCameraFractionInZetaHighElev(-1.);
          theApLight.SetMaxCloudCameraFractionWithLidarInZetaHighElev(-1.);
          return;
        }
      }
    }
  }
  // write maxCloudFraction to adst class
  FdRecApertureLight& theApLight = theFd.GetFdRecApertureLight();
  theApLight.SetMaxCloudCameraFractionInZeta(maxCloudFraction);
  theApLight.SetMaxCloudCameraFractionWithLidarInZeta(maxCloudFractionLidar);
  theApLight.SetMaxCloudCameraFractionInZetaHighElev(maxCloudFractionHighElevation);
  theApLight.SetMaxCloudCameraFractionWithLidarInZetaHighElev(maxCloudFractionLidarHighElevation);
}


bool
FD2ADST::FillFdProfile(const fdet::Eye& /*dEye*/, const fevt::Eye& eye, FDEvent& theFd)
{
  const fevt::EyeRecData& eyeRec = eye.GetRecData();
  const evt::ShowerFRecData& showerFRec = eyeRec.GetFRecShower();

  if (!showerFRec.HasLongitudinalProfile())
    return false;

  // ####### Getting reconstructed light contributions at diaphragm #####
  FdRecApertureLight& theProfile = theFd.GetFdRecApertureLight();
  const vector<double>& dia_time = theProfile.GetTime();
  const vector<double>& dia_time_err = theProfile.GetTimeError();

  const utl::TabulatedFunctionErrors& diaPhotonTrace =
    eyeRec.GetLightFlux(fevt::FdConstants::eTotal);

  unsigned int time_size = diaPhotonTrace.GetNPoints();
  vector<double> dia_fluor(time_size);
  vector<double> dia_fluor_err(time_size);
  vector<double> dia_direct_cher(time_size);
  vector<double> dia_direct_cher_err(time_size);
  vector<double> dia_mie_cher(time_size);
  vector<double> dia_mie_cher_err(time_size);
  vector<double> dia_ray_cher(time_size);
  vector<double> dia_ray_cher_err(time_size);
  vector<double> dia_mult_cher(time_size);
  vector<double> dia_mult_fluor(time_size);

  if (eyeRec.HasLightFlux(fevt::FdConstants::eFluorDirect)) {
    const TabulatedFunctionErrors& directFluorescence =
      eyeRec.GetLightFlux(fevt::FdConstants::eFluorDirect);
    std::copy(directFluorescence.YBegin(), directFluorescence.YEnd(), &dia_fluor.front());
  }

  if (eyeRec.HasLightFlux(fevt::FdConstants::eCherDirect)) {
    const TabulatedFunctionErrors& directCherenkov =
      eyeRec.GetLightFlux(fevt::FdConstants::eCherDirect);
    std::copy(directCherenkov.YBegin(), directCherenkov.YEnd(), &dia_direct_cher.front());
  }

  if (eyeRec.HasLightFlux(fevt::FdConstants::eCherMieScattered)) {
    const TabulatedFunctionErrors& mieCherenkov =
      eyeRec.GetLightFlux(fevt::FdConstants::eCherMieScattered);
    std::copy(mieCherenkov.YBegin(), mieCherenkov.YEnd(), &dia_mie_cher.front());
  }

  if (eyeRec.HasLightFlux(fevt::FdConstants::eCherRayleighScattered)) {
    const TabulatedFunctionErrors& rayCherenkov =
      eyeRec.GetLightFlux(fevt::FdConstants::eCherRayleighScattered);
    std::copy(rayCherenkov.YBegin(), rayCherenkov.YEnd(), &dia_ray_cher.front());
  }

  if (eyeRec.HasLightFlux(fevt::FdConstants::eCherMultScattered)) {
    const TabulatedFunctionErrors& rayCherenkov =
      eyeRec.GetLightFlux(fevt::FdConstants::eCherMultScattered);
    std::copy(rayCherenkov.YBegin(), rayCherenkov.YEnd(), &dia_mult_cher.front());
  }

  if (eyeRec.HasLightFlux(fevt::FdConstants::eFluorMultScattered)) {
    const TabulatedFunctionErrors& tmp =
      eyeRec.GetLightFlux(fevt::FdConstants::eFluorMultScattered);
    std::copy(tmp.YBegin(), tmp.YEnd(), &dia_mult_fluor.front());
  }

  for (unsigned int i = 0, n = dia_time.size(); i < n; ++i) {
    const double timeBin = 2 * dia_time_err[i];
    dia_fluor[i] /= timeBin;
    dia_direct_cher[i] /= timeBin;
    dia_mie_cher[i] /= timeBin;
    dia_ray_cher[i] /= timeBin;
    dia_mult_fluor[i] /= timeBin;
    dia_mult_cher[i] /= timeBin;
  }

  theProfile.SetFluorLightAtAperture(dia_fluor);
  theProfile.SetCherLightAtAperture(dia_direct_cher);
  theProfile.SetMieCherLightAtAperture(dia_mie_cher);
  theProfile.SetRayCherLightAtAperture(dia_ray_cher);
  theProfile.SetMultScatCherLightAtAperture(dia_mult_cher);
  theProfile.SetMultScatFluorLightAtAperture(dia_mult_fluor);

  // Calculate Cherenkov fraction
  double totalLight = 0;
  double totalCherenkov = 0;
  for (unsigned int i = 0, n = dia_time.size(); i < n; ++i) {
    double timeBin = 2 * dia_time_err[i];
    totalLight += (dia_fluor[i] + dia_mult_fluor[i] +
                   dia_direct_cher[i] + dia_mie_cher[i] +
                   dia_ray_cher[i] + dia_mult_cher[i]) * timeBin;
    totalCherenkov += (dia_direct_cher[i] + dia_mie_cher[i] + dia_ray_cher[i] +
                       dia_mult_cher[i]) * timeBin;
  }

  FdRecShower& theShower = theFd.GetFdRecShower();
  if (totalCherenkov > 0 && totalLight > 0)
    theProfile.SetCherenkovFraction(100 * totalCherenkov / totalLight);
  else
    theProfile.SetCherenkovFraction(-1);

  //cout << " Fill rec viewing angle" << endl;
  FillViewingAngle(&theProfile, &(theFd.GetFdRecGeometry()), theShower.GetXmaxChi());

  // ######  Getting the Profile (Electrons vs Depths)  #######

  const TabulatedFunctionErrors& electronTrace = showerFRec.GetLongitudinalProfile();
  const unsigned int esize = electronTrace.GetNPoints();

  vector<double> depth(electronTrace.XBegin(), electronTrace.XEnd());
  vector<double> depth_err(electronTrace.XErrBegin(), electronTrace.XErrEnd());
  vector<double> electrons(electronTrace.YBegin(), electronTrace.YEnd());
  vector<double> electrons_err(electronTrace.YErrBegin(), electronTrace.YErrEnd());

  bool hasdEdX = false;
  if (esize > 0) {
    for (unsigned int i = 0; i < esize; ++i) {
      depth[i]     /= g/cm2;
      depth_err[i] /= g/cm2;
    }

    theShower.SetDepth(depth, depth_err);
    theShower.SetElectrons(electrons, electrons_err);

    IncreaseFdRecLevel(theFd, eHasLongitudinalProfile);

    if (showerFRec.HasEnergyDeposit()) {  // reusing the electrons vectors!

      const TabulatedFunctionErrors& dedxTrace = showerFRec.GetEnergyDeposit();
      const unsigned int dedx_size = dedxTrace.GetNPoints();

      if (dedx_size != esize) {
        ERROR("FillFdProfile - Error dedx_size!=size (?)");
        return false;
      }

      std::copy(dedxTrace.YBegin(),dedxTrace.YEnd(), &electrons.front());
      std::copy(dedxTrace.YErrBegin(),dedxTrace.YErrEnd(), &electrons_err.front());

      for (unsigned int i = 0; i < esize; ++i) {
        electrons[i]     /= PeV/(g/cm2);
        electrons_err[i] /= PeV/(g/cm2);
      }

      theShower.SetEnergyDeposit(electrons, electrons_err);
      IncreaseFdRecLevel(theFd, eHasdEdXProfile);
      hasdEdX = true;
    }
  }

  return hasdEdX;
}


bool
FD2ADST::FillGaisserHillas(const fdet::Eye& detEye, const fevt::Eye& eye, FDEvent& theFd)
{
  // reconstructed Gaisser-Hillas profile
  const fevt::EyeRecData& eyeRec = eye.GetRecData();
  const evt::ShowerFRecData& showerFRec = eyeRec.GetFRecShower();
  if (!showerFRec.HasGHParameters())
    return false;

  const double nmax     = showerFRec.GetGHParameters().GetNMax();
  const double nmaxErr  = showerFRec.GetGHParameters().GetNMaxError();
  const double xmax     = showerFRec.GetGHParameters().GetXMax();
  const double chisq    = showerFRec.GetGHParameters().GetChiSquare();
  const double ndof     = showerFRec.GetGHParameters().GetNdof();
  const double rhoNmaxXmax =
    showerFRec.GetGHParameters().GetNMaxXMaxCorrelation();

  double xzero = 0;
  double xzeroErr = 0;
  double lam = 0;
  double lamErr = 0;
  double fwhm = 0;
  double fwhmErr = 0;
  double asym = 0;
  double asymErr = 0;
  double uspL = 0;
  double uspLErr = 0;
  double uspR = 0;
  double uspRErr = 0;
  double rhoNMaxLambda = 0;
  double rhoNMaxX0 = 0;
  double rhoXMaxLambda = 0;
  double rhoXMaxX0 = 0;
  double rhoLambdaX0 = 0;

  EGaisserHillasType ghType = eGHUndefinedType;

  if (showerFRec.HasGHParameters<evt::GaisserHillas2Parameter>()) {

    ghType = eGHTwoPar;

  } else if (showerFRec.HasGHParameters<evt::GaisserHillas4Parameter>()) {

    const evt::GaisserHillas4Parameter& gh4 = showerFRec.GetGHParameters<evt::GaisserHillas4Parameter>();
    using namespace evt::gh;
    switch (gh4.GetFunctionType()) {
    case evt::gh::eClassic:
      ghType = eGHFourParClassic;
      break;
    case evt::gh::eWidth:
      ghType = eGHFourParWidth;
      break;
    case evt::gh::eUSP:
      ghType = eGHFourParUSP;
      break;
    default:
      {
        ostringstream errMsg;
        errMsg << " no ADST enum for GH type "
               << evt::gh::GetFunctionTypeName(gh4.GetFunctionType());
        throw DoesNotComputeException(errMsg.str());
      }
    }

    xzero = gh4.GetShapeParameter(eX0) / (g/cm2);
    xzeroErr = gh4.GetShapeParameterError(eX0) / (g/cm2);
    fwhm = gh4.GetShapeParameter(eFWHM) / (g/cm2);
    fwhmErr = gh4.GetShapeParameterError(eFWHM) / (g/cm2);
    asym = gh4.GetShapeParameter(eAsym);
    asymErr = gh4.GetShapeParameterError(eAsym);
    uspL = gh4.GetShapeParameter(eUSP_L) / (g/cm2);
    uspLErr = gh4.GetShapeParameterError(eUSP_L) / (g/cm2);
    uspR = gh4.GetShapeParameter(eUSP_R);
    uspRErr = gh4.GetShapeParameterError(eUSP_R);

    lam = gh4.GetShapeParameter(eLambda) / (g/cm2);
    // some OG module does not set lambda, but use kLambdaGH
    if (lam == 0)
      lam = kLambdaGH / (g/cm2);
    lamErr = gh4.GetShapeParameterError(eLambda) / (g/cm2);

    if (gh4.GetFunctionType() == evt::gh::eClassic) {
      rhoNMaxLambda = gh4.GetCorrelationNMaxShapeParameter(eLambda);
      rhoNMaxX0 = gh4.GetCorrelationNMaxShapeParameter(eX0);
      rhoXMaxLambda = gh4.GetCorrelationXMaxShapeParameter(eLambda);
      rhoXMaxX0 = gh4.GetCorrelationXMaxShapeParameter(eX0);
      rhoLambdaX0 = gh4.GetCorrelationShapeParameters();
    }

  } else if (showerFRec.HasGHParameters<evt::GaisserHillas6Parameter>()) {  // this should never happen, not used in reco.
    const evt::GaisserHillas6Parameter& gh6 = showerFRec.GetGHParameters<evt::GaisserHillas6Parameter>();
    WARNING("6 Parameter GaisserHillas-fit? Are you serious?");
    xzero = gh6.GetXZero() / (g/cm2);
    xzeroErr = gh6.GetXZeroError() / (g/cm2);
  }

  // Calculate Xmax Chi, etc
  const FdRecGeometry& recGeo = theFd.GetFdRecGeometry();
  const CoordinateSystemPtr referenceCS =
    det::Detector::GetInstance().GetReferenceCoordinateSystem();
  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const Vector recSDP(1, recGeo.GetSDPTheta(), recGeo.GetSDPPhi(), eyeCS, Vector::kSpherical);

  // calculate chi(xmax)
  FdRecShower& theShower = theFd.GetFdRecShower();
  theShower.SetGHType(ghType);

  const TVector3& coreSiteCS = theShower.GetCoreSiteCS();
  Point recCore(coreSiteCS.X()*m, coreSiteCS.Y()*m, coreSiteCS.Z()*m, referenceCS);
  const TVector3& axisSiteCS = theShower.GetAxisSiteCS();
  const Vector recAxis(1, axisSiteCS.Theta(), axisSiteCS.Phi(), referenceCS, Vector::kSpherical);

  const double xMaxChi = otoa::fd::ConvertXToChi(xmax, detEye, recCore, recAxis, recSDP);

  if (ndof > 0) {
    const double gcm2 = g/cm2;
    FdRecShower& theShower = theFd.GetFdRecShower();

    const double conv = (PeV/gcm2);
    theShower.SetdEdXmax(nmax / conv);
    theShower.SetdEdXError(nmaxErr / conv);

    theShower.SetXmax(xmax / gcm2, xMaxChi);

    theShower.SetXmaxError(showerFRec.GetXmaxError() / gcm2, 0.);
    if (showerFRec.HasXmaxError(evt::ShowerFRecData::eAtmospheric))
      theShower.SetXmaxAtmError(showerFRec.GetXmaxError(evt::ShowerFRecData::eAtmospheric) / gcm2);
    theShower.SetX0(xzero);
    theShower.SetX0Error(xzeroErr);
    theShower.SetLambda(lam);
    theShower.SetLambdaError(lamErr);
    theShower.SetFWHM(fwhm);
    theShower.SetFWHMError(fwhmErr);
    theShower.SetAsymmetry(asym);
    theShower.SetAsymmetryError(asymErr);
    theShower.SetUspL(uspL);
    theShower.SetUspLError(uspLErr);
    theShower.SetUspR(uspR);
    theShower.SetUspRError(uspRErr);
    theShower.SetGHChi2(chisq);
    theShower.SetGHNdf(int(ndof));

    theShower.SetdEdXmaxXMaxCorrelation(rhoNmaxXmax);
    theShower.SetdEdXmaxLambdaCorrelation(rhoNMaxLambda);
    theShower.SetdEdXmaxX0Correlation(rhoNMaxX0);
    theShower.SetXMaxLambdaCorrelation(rhoXMaxLambda);
    theShower.SetXMaxX0Correlation(rhoXMaxX0);
    theShower.SetLambdaX0Correlation(rhoLambdaX0);

    // Set multiple GH parameters if they exist
    if (showerFRec.HasMultipleGHParameters()) {
      const evt::MultipleGaisserHillasParameters& mghPars =
        showerFRec.GetMultipleGHParameters();
      if (mghPars.GetVirtualParameters().size() == 2) {
        const double xmax1 = (*(mghPars.GetVirtualParameters()[0])).GetXMax();
        const double xmax2 = (*(mghPars.GetVirtualParameters()[1])).GetXMax();
        const double mgChi2 = mghPars.GetChi2Improvement();

        theShower.SetDGHXmax1(xmax1 / gcm2);
        theShower.SetDGHXmax2(xmax2 / gcm2);
        theShower.SetDGHChi2Improvement(mgChi2);
      }
    }

    return true;
  }

  return false;
}


void
FD2ADST::FillAtmosphericProfileVars(const fdet::Eye& detEye, const fevt::Eye& eye, FDEvent& theFd)
{
  // reconstructed Gaisser-Hillas profile
  const fevt::EyeRecData& eyeRec = eye.GetRecData();
  const evt::ShowerFRecData& showerFRec = eyeRec.GetFRecShower();
  FdRecShower& recShower = theFd.GetFdRecShower();
  const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

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

  const utl::Point& core = showerFRec.GetCorePosition();
  const utl::Vector& axis = showerFRec.GetAxis();
  const double xmax = showerFRec.GetGHParameters().GetXMax();

  double distCoreXmax = 0;
  double xMaxHeight = 0;
  double xMaxHeightAboveSeaLevel = 0;
  try {
    try {
      // use the inclined atmosphere model
      const atm::ProfileResult& distVsSlantDepth = atmo.EvaluateDistanceVsSlantDepth();
      const atm::ProfileResult& heightVsSlantDepth = atmo.EvaluateHeightVsSlantDepth();

      xMaxHeight = (heightVsSlantDepth.MaxX() > xmax)
        ? heightVsSlantDepth.Y(xmax) : 0.;

      // DEF: p = core + distance * axis
      distCoreXmax = (distVsSlantDepth.MaxX() > xmax)
        ? -distVsSlantDepth.Y (xmax) : 0.;

      xMaxHeightAboveSeaLevel = xMaxHeight; // includes core height!
    } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
      // use the flat atmosphere geometry
      const atm::ProfileResult& heightVsDepth = atmo.EvaluateHeightVsDepth();
      const double fdTheta = recShower.GetZenith();
      const double Xvert = xmax * cos(fdTheta);
      xMaxHeightAboveSeaLevel = heightVsDepth.Y(Xvert);
      xMaxHeight = xMaxHeightAboveSeaLevel - wgs84.PointToLatitudeLongitudeHeight(core).get<2>();
      distCoreXmax = xMaxHeight / cos(fdTheta);
    }
  } catch (OutOfBoundException& e) {
    // ignore
  }

  const Point pointXmax = core + axis * distCoreXmax;
  recShower.SetCoreXmaxDist(distCoreXmax);

  double distXmax = (pointXmax - eyePos).GetMag();
  // The following can happen if theta > 180 deg ->
  // heightVsDepth(-xx)= IsNan but no OutOfBoundException
  if (std::isnan(distXmax))
    distXmax = 0;
  recShower.SetDistXmax(distXmax);

  recShower.SetHeightXmax(xMaxHeightAboveSeaLevel);

  // ------calculate the heights corresponding to the depth vector
  const vector<double>& depthVec = recShower.GetDepth();

  vector<double> xHeightAboveSeaLevel;
  try {
    try {
      // use the inclined atmosphere model
      const atm::ProfileResult& heightVsSlantDepth = atmo.EvaluateHeightVsSlantDepth();
      for (unsigned int iDepth = 0; iDepth < depthVec.size(); ++iDepth) {
        const double depth = depthVec[iDepth] * g/cm2;
        double heightAboveSeaLevel;
        if (heightVsSlantDepth.MaxX() > depth && heightVsSlantDepth.MinX() < depth)
          heightAboveSeaLevel = heightVsSlantDepth.Y(depth); // includes core height!
        else
          heightAboveSeaLevel = 0;
        xHeightAboveSeaLevel.push_back(heightAboveSeaLevel);
      }
    } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
      xHeightAboveSeaLevel.clear();
      // use the flat atmosphere geometry
      const atm::ProfileResult& heightVsDepth = atmo.EvaluateHeightVsDepth();
      const double fdTheta = recShower.GetZenith();
      for (unsigned int iDepth = 0; iDepth < depthVec.size(); ++iDepth) {
        const double depth = depthVec[iDepth] * g/cm2;
        const double Xvert = depth*cos(fdTheta);
        const double heightAboveSeaLevel = heightVsDepth.Y(Xvert);
        xHeightAboveSeaLevel.push_back(heightAboveSeaLevel);
      }
    }
  } catch (OutOfBoundException& e) {/*ignore*/}

  recShower.SetHeight(xHeightAboveSeaLevel);

  std::vector<double> lambda;
  const double refLambda = det::Detector::GetInstance().GetFDetector().GetReferenceLambda();
  lambda.push_back(refLambda);
  const atm::AttenuationResult mCAtt = atmo.EvaluateMieAttenuation(eyePos, pointXmax, lambda);
  const atm::AttenuationResult rCAtt = atmo.EvaluateRayleighAttenuation(eyePos, pointXmax, lambda);
  const utl::TabulatedFunctionErrors& rayCAtt = rCAtt.GetTransmissionFactor();
  const utl::TabulatedFunctionErrors& mieCAtt = mCAtt.GetTransmissionFactor();

  double attXmax = mieCAtt.GetY(0)*rayCAtt.GetY(0);
  if (std::isnan(attXmax))
    attXmax = 0;
  double mieCAttXmax = mieCAtt.GetY(0);
  if (std::isnan(mieCAttXmax))
    mieCAttXmax = 0;

  recShower.SetAttXmax(attXmax);
  recShower.SetMieAttXmax(mieCAttXmax);

  // ------------- Calculate Mie/Rayleigh correction
  const double dX = 25*g/cm2;
  const std::vector<double> referenceWaveLength(1, refLambda);

  const TabulatedFunctionErrors& dedxTrace = showerFRec.GetEnergyDeposit();
  const unsigned int dedx_size = dedxTrace.GetNPoints();
  if (dedx_size < 1)
    return;

  const double depth1 = dedxTrace.GetX(0);
  const double depth2 = dedxTrace.GetX(dedx_size-1);
  const double minDepth = fmin(depth1, depth2);
  const double maxDepth = fmax(depth1, depth2);
  const evt::VGaisserHillasParameter& ghParameters = showerFRec.GetGHParameters();

  double rayleighCorrection = 0;
  double mieCorrection      = 0;
  try {
    const atm::ProfileResult& distVsSlantDepth = atmo.EvaluateDistanceVsSlantDepth();

    double profileSum       = 0;
    double profileMieAttSum = 0;
    double profileRayAttSum = 0;

    for (double currentDepth = maxDepth; currentDepth > minDepth; currentDepth -= dX) {

      if (currentDepth > distVsSlantDepth.MinX() &&
          currentDepth < distVsSlantDepth.MaxX()) {

        const double dEdX = ghParameters.Eval(currentDepth);
        const double distance = -distVsSlantDepth.Y(currentDepth);
        const Point pointOnShower = core + distance*axis;

        const atm::AttenuationResult mCAtt = atmo.EvaluateMieAttenuation(eyePos, pointOnShower, referenceWaveLength);
        const utl::TabulatedFunctionErrors& mieCAtt = mCAtt.GetTransmissionFactor();
        const atm::AttenuationResult rCAtt = atmo.EvaluateRayleighAttenuation(eyePos, pointOnShower, referenceWaveLength);
        const utl::TabulatedFunctionErrors& rayCAtt = rCAtt.GetTransmissionFactor();

        profileSum       += dEdX;
        profileMieAttSum += dEdX * mieCAtt.GetY(0);
        profileRayAttSum += dEdX * rayCAtt.GetY(0);
      }

    }

    if (profileMieAttSum > 0)
      mieCorrection = profileSum / profileMieAttSum;
    if (profileRayAttSum > 0)
      rayleighCorrection = profileSum / profileRayAttSum;
  } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) {
    ERROR("calculation of transmission factors not implemented for flat atmosphere...");
    mieCorrection = 0;
    rayleighCorrection = 0;
  }
  recShower.SetMieCorrection(mieCorrection);
  recShower.SetRayleighCorrection(rayleighCorrection);

  // -------- geometric variables useful for cuts ---------
  if (!depthVec.empty()) {
    recShower.SetXTrackMin(depthVec.front(), 0.); // FIXME second argument==chi. WTF?
    recShower.SetXTrackMax(depthVec.back());
//#warning Maybe this needs to be made aware that the time-at-diaphragm data might only be available per telescope in a future reconstruction
    const vector<double>& aptime = theFd.GetFdRecApertureLight().GetTime();
    if (aptime.size() > 1)
      recShower.SetTimTrackObs(aptime.back() - aptime.front());
    else
      recShower.SetTimTrackObs(0.);
    // angular track length is set in FillRecPixel...
  }
}


bool
FD2ADST::FillTelSimData(const fevt::Eye& eye, FDEvent& fd)
{
  const int eyeId = fd.GetEyeId();

  det::Detector& detector = det::Detector::GetInstance();
  const fdet::FDetector& detFd = detector.GetFDetector();
  const fdet::Eye& detEye = detFd.GetEye(eyeId);
  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();

  const double minWavelength = detFd.GetModelMinWavelength();
  const double maxWavelength = detFd.GetModelMaxWavelength();

  const atm::Atmosphere& atmo = detector.GetAtmosphere();
  const vector<double>& WLengthFluo = atmo.GetWavelengths(atm::Atmosphere::eFluorescence);
  const vector<double>& WLengthCkov = atmo.GetWavelengths(atm::Atmosphere::eCerenkov);
  unsigned int NWLengthFluo = WLengthFluo.size();
  unsigned int NWLengthCkov = WLengthCkov.size();

  const bool hasEyeTrigger = eye.HasHeader();
  TimeStamp eyeTriggerTime;
  if (hasEyeTrigger) {
    const TimeStamp eyeGPSTime = eye.GetHeader().GetTimeStamp();
    eyeTriggerTime = eyeGPSTime - TimeInterval(1000*100*utl::ns);
  }

  // ----------------  MC light profile -------------------

  bool oneTel = false;

  // select if only save just "telescopes WITH reco. data", or "ALL simulated telescopes IN FOV"
  ComponentSelector::Status select = ComponentSelector::eHasData;
  if (!Config().DropUntriggeredMCTelescopes())
    select = ComponentSelector::eInDAQ;

  for (auto teliter = eye.TelescopesBegin(select); teliter != eye.TelescopesEnd(select); ++teliter) {

    const int telId = teliter->GetId();

    double totalGenLight = 0;
    double totalGenCherenkov = 0;

    if (!teliter->HasSimData())
      continue;
    const fevt::TelescopeSimData& simTel = teliter->GetSimData();

    if (!simTel.HasDistanceTrace())
      continue;
    const TraceD& distanceTrace = simTel.GetDistanceTrace();

    const bool hasNoPixels = (teliter->PixelsBegin(select) == teliter->PixelsEnd(select));
    if (hasNoPixels)
      continue;

    const uint nBins = distanceTrace.GetSize();
    double binWidth = distanceTrace.GetBinning(); // [ns] this can be changed via xml by user
    const double telStartTime = hasEyeTrigger ? (simTel.GetPhotonsStartTime()-eyeTriggerTime).GetInterval() : 0;

    vector<double> timeVec(nBins);
    vector<double> lightDiaTotal(nBins);
    vector<double> lightDiaFluo(nBins);
    vector<double> lightDiaCkovD(nBins);
    vector<double> lightDiaCkovM(nBins);
    vector<double> lightDiaCkovR(nBins);
#ifdef PRIM_CHER
    vector<double> lightDiaCkovPrimD(nBins);
    vector<double> lightDiaCkovPrimM(nBins);
    vector<double> lightDiaCkovPrimR(nBins);
    bool HasCkovPrimD = false;
    bool HasCkovPrimM = false;
    bool HasCkovPrimR = false;
#endif
    bool HasTotal = false;
    bool HasFluo = false;
    bool HasCkovD = false;
    bool HasCkovM = false;
    bool HasCkovR = false;

    const fdet::Telescope& detTel = detEye.GetTelescope (telId);
    //const Point& telPos = detTel.GetPosition();
    const CoordinateSystemPtr telCS = detTel.GetTelescopeCoordinateSystem();
    // const double diaphragmArea = detTel.GetDiaphragmArea(); // unused -- silence warning
    // const double refLambda = detFd.GetReferenceLambda(); // unused -- silence warning

    for (uint iBin = 0; iBin < nBins; ++iBin) {

      timeVec[iBin] = (telStartTime + double(0.5 + iBin) * binWidth) / binWidth; // bin center in units of [bins]

      // Loop fluorescence wavelengths
      for (unsigned int iwl = 0; iwl < NWLengthFluo; ++iwl) {

        const double wavelength = WLengthFluo[iwl];
        if (wavelength<minWavelength || wavelength>maxWavelength)
          continue;

        const double norm = detTel.GetModelRelativeEfficiency(wavelength);

        if (simTel.HasPhotonTrace(fevt::FdConstants::eFluorDirect, iwl)) {
          const TraceD& photonTrace =
            simTel.GetPhotonTrace (fevt::FdConstants::eFluorDirect, iwl);
          const double vBin = photonTrace[iBin] * norm;
          lightDiaFluo[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          totalGenLight += vBin;
          HasFluo = true;
          HasTotal = true;
        }
      }

      // Loop Cherenkov wavelengths
      for (unsigned int iwl = 0; iwl < NWLengthCkov; ++iwl) {

        const double wavelength = WLengthCkov[iwl];
        if (wavelength<minWavelength || wavelength>maxWavelength)
          continue;

        const double norm = detTel.GetModelRelativeEfficiency(wavelength);

        if (simTel.HasPhotonTrace(fevt::FdConstants::eCherDirect, iwl)) {
          const TraceD& photonTrace =
            simTel.GetPhotonTrace (fevt::FdConstants::eCherDirect, iwl);
          const double vBin = iBin<photonTrace.GetSize() ? photonTrace[iBin] * norm : 0;
          lightDiaCkovD[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          totalGenLight += vBin;
          totalGenCherenkov += vBin;
          HasCkovD = true;
          HasTotal = true;
        }

        if (simTel.HasPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl)) {
          const TraceD& photonTrace =
            simTel.GetPhotonTrace(fevt::FdConstants::eCherMieScattered, iwl);
          const double vBin = photonTrace[iBin] * norm;
          lightDiaCkovM[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          totalGenLight += vBin;
          totalGenCherenkov += vBin;
          HasCkovM = true;
          HasTotal = true;
        }

        if (simTel.HasPhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl)) {
          const TraceD& photonTrace =
            simTel.GetPhotonTrace(fevt::FdConstants::eCherRayleighScattered, iwl);
          const double vBin = photonTrace[iBin] * norm;
          lightDiaCkovR[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          totalGenLight += vBin;
          totalGenCherenkov += vBin;
          HasCkovR = true;
          HasTotal = true;
        }

#ifdef PRIM_CHER
        if (simTel.HasPhotonTrace(fevt::FdConstants::ePrimaryCherDirect, iwl)) {
          const TraceD& photonTrace = simTel.GetPhotonTrace(fevt::FdConstants::ePrimaryCherDirect, iwl);
          const double vBin = photonTrace[iBin] * norm;
          lightDiaCkovPrimD[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          HasCkovPrimD = true;
          HasTotal = true;
        }

        if (simTel.HasPhotonTrace(fevt::FdConstants::ePrimaryCherMieScattered, iwl)) {
          const TraceD& photonTrace = simTel.GetPhotonTrace(fevt::FdConstants::ePrimaryCherMieScattered, iwl);
          const double vBin = photonTrace[iBin] * norm;
          lightDiaCkovPrimM[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          HasCkovPrimM = true;
          HasTotal = true;
        }

        if (simTel.HasPhotonTrace(fevt::FdConstants::ePrimaryCherRayleighScattered, iwl)) {
          const TraceD& photonTrace = simTel.GetPhotonTrace(fevt::FdConstants::ePrimaryCherRayleighScattered, iwl);
          const double vBin = photonTrace[iBin] * norm;
          lightDiaCkovPrimR[iBin] += vBin;
          lightDiaTotal[iBin] += vBin;
          HasCkovPrimR = true;
          HasTotal = true;
        }
#endif
      }
    }

    if (totalGenLight <= 0)  // no light here....
      continue;

    // rebin traces if wished
    if (Config().GetRebinSimTelescopeTraces() > 0) {

      /* The simulated light binning can be very fine.
         For display purposes enlarge bin width */
      int nCombine = max(1, int(Config().GetRebinSimTelescopeTraces() / binWidth));
      {
        const int nOrg = timeVec.size();
        const int nReduced = max(2, nOrg/nCombine); // require minimum 2 bins leftover
        nCombine = max(1, int(double(nOrg)/double(nReduced)));
      }

      binWidth *= nCombine;

      int iCombine = 0;
      double c1 = 0;
      double c2 = 0;
      double c3 = 0;
      double c4 = 0;
      double c5 = 0;
      double c6 = 0;
      vector<double> time_comb;
      vector<double> fluo_comb;
      vector<double> dirC_comb;
      vector<double> mieC_comb;
      vector<double> rayC_comb;
      vector<double> total_comb;
      for (unsigned int i = 0; i < timeVec.size(); ++i) {
        c1 = (telStartTime + double(0.5 + time_comb.size()) * binWidth) / binWidth; // bin center in units of [bins]
        c2 += lightDiaFluo[i];
        c3 += lightDiaCkovD[i];
        c4 += lightDiaCkovM[i];
        c5 += lightDiaCkovR[i];
        c6 += lightDiaTotal[i];
        if (++iCombine >= nCombine) {
          time_comb.push_back(c1);
          fluo_comb.push_back(c2);
          dirC_comb.push_back(c3);
          mieC_comb.push_back(c4);
          rayC_comb.push_back(c5);
          total_comb.push_back(c6);
          c1 = c2 = c3 = c4 = c5 = c6 = 0;
          iCombine = 0;
        }
      }
      // the LAST (incomplete) bin !
      if (iCombine > 0) {
        time_comb.push_back(c1);
        fluo_comb.push_back(c2);
        dirC_comb.push_back(c3);
        mieC_comb.push_back(c4);
        rayC_comb.push_back(c5);
        total_comb.push_back(c6);
      }
      timeVec = time_comb;
      lightDiaFluo = fluo_comb;
      lightDiaCkovD = dirC_comb;
      lightDiaCkovM = mieC_comb;
      lightDiaCkovR = rayC_comb;
      lightDiaTotal = total_comb;

    }

    oneTel = true;
    if (!fd.HasTelescopeData(telId))
      fd.MakeTelescopeData(telId);

    FdTelescopeData& telSimD = fd.GetTelescopeData(telId);

    FdGenApertureLight& fdProfile = telSimD.GetGenApertureLight();
    fdProfile.SetTimeBinning(binWidth/ns);

    if (HasTotal)
      fdProfile.SetTotalLightAtAperture(lightDiaTotal);
    if (HasFluo)
      fdProfile.SetFluorLightAtAperture(lightDiaFluo);
    if (HasCkovD)
      fdProfile.SetCherLightAtAperture(lightDiaCkovD);
    if (HasCkovR)
      fdProfile.SetRayCherLightAtAperture(lightDiaCkovR);
    if (HasCkovM)
      fdProfile.SetMieCherLightAtAperture(lightDiaCkovM);
#ifdef PRIM_CHER
    if (HasCkovPrimD)
      fdProfile.SetCherPrimaryLightAtAperture(lightDiaCkovPrimD);
    if (HasCkovPrimR)
      fdProfile.SetRayCherPrimaryLightAtAperture(lightDiaCkovPrimR);
    if (HasCkovPrimM)
      fdProfile.SetMieCherPrimaryLightAtAperture(lightDiaCkovPrimM);
#endif

    if (HasTotal || HasFluo || HasCkovD || HasCkovR || HasCkovM)
      fdProfile.SetTime(timeVec);

    //cout << " Fill gen viewing angle" << endl;
    FillViewingAngle((FdApertureLight*)&(telSimD.GetGenApertureLight()), (FdGeometry*)&(fd.GetGenGeometry()), fd.GetGenShower().GetXmaxChi());

    telSimD.GetGenApertureLight().SetCherenkovFraction(0);
    if (totalGenCherenkov > 0 && totalGenLight > 0)
      telSimD.GetGenApertureLight().SetCherenkovFraction(100.*totalGenCherenkov/totalGenLight);

  }

  return oneTel;
}


void
FD2ADST::FillTelRecData(const fevt::Eye& theEye, FDEvent& theFd)
{
  const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(theEye);

  for (auto teliter = theEye.TelescopesBegin(ComponentSelector::eHasData);
       teliter != theEye.TelescopesEnd(ComponentSelector::eHasData); ++teliter) {

    const unsigned int telId = teliter->GetId();

    if (!teliter->HasRecData())
      continue;

    const fevt::TelescopeRecData& oRecData = teliter->GetRecData();
    if (!theFd.HasTelescopeData(telId))
      theFd.MakeTelescopeData(telId);
    FdTelescopeData& telRecD = theFd.GetTelescopeData(telId);

    // Fill FdRecGeometry
    // FIXME fill it completely, reduce code duplication?
    FdRecGeometry& recGeo = telRecD.GetRecGeometry();

    const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
    const double SDPphi   = oRecData.GetSDP().GetPhi(eyeCS);
    const double SDPtheta = oRecData.GetSDP().GetTheta(eyeCS);
    recGeo.SetSDP(SDPtheta, SDPphi);

    recGeo.SetAxis(oRecData.GetTZero(), oRecData.GetRp(), oRecData.GetChiZero());
    recGeo.SetT0Error(oRecData.GetTZeroError());
    recGeo.SetRpError(oRecData.GetRpError());
    recGeo.SetChi0Error(oRecData.GetChiZeroError());
    recGeo.SetChi0T0Correlation(oRecData.GetChi0TZeroCorrelation());
    recGeo.SetChi0RpCorrelation(oRecData.GetRpChi0Correlation());
    recGeo.SetRpT0Correlation(oRecData.GetRpTZeroCorrelation());
    recGeo.SetTimeFitChi2(oRecData.GetTimeFitChiSquare());
    recGeo.SetTimeFitNdof(oRecData.GetTimeFitNDof());

    // Fill FdRecApertureLight
    FdRecApertureLight& recLight = telRecD.GetRecApertureLight();

    recLight.SetZeta(oRecData.GetZeta());

    if (!oRecData.HasLightFlux(fevt::FdConstants::eTotal))
      continue;

    const TabulatedFunctionErrors& diaPhotonTrace =
      oRecData.GetLightFlux(fevt::FdConstants::eTotal);

    const double binSize = 100*ns;
    const unsigned int time_size = diaPhotonTrace.GetNPoints();
    vector<double> diaTime(time_size), diaTimeErr(time_size),
      diaPhotons(time_size), diaPhotonsErr(time_size);

    recLight.SetTimeBinning(binSize/ns);

    // Note: We need to find the first bin that was included in the
    //       profile fit because that's where our information about
    //       the light components starts.
    // Note2: This assumes the light components all start at the same bin!
    TabulatedFunctionErrors directFluorescence;
    if (oRecData.HasLightFlux(fevt::FdConstants::eFluorDirect))
      directFluorescence = oRecData.GetLightFlux(fevt::FdConstants::eFluorDirect);
    unsigned int firstComponentBin = 0;
    const double firstComponentBinTime = (directFluorescence.GetNPoints() == 0) ? 0 : directFluorescence.GetX(0);

    // Fill l@aperture and the time vector, find first component bin
    for (unsigned int i = 0; i < time_size; ++i) {
      if (diaPhotonTrace[i].X() <= firstComponentBinTime)
        firstComponentBin = i;

      diaTime[i]    = diaPhotonTrace[i].X() / binSize;
      diaTimeErr[i] = diaPhotonTrace.GetXErr(i) / binSize;

      if (diaTimeErr[i] == 0)
        diaTimeErr[i] = 1;

      diaPhotons[i]    = diaPhotonTrace[i].Y() / (2. * diaTimeErr[i]);
      diaPhotonsErr[i] = diaPhotonTrace.GetYErr(i) / (2. * diaTimeErr[i]);
    }

    recLight.SetTime(diaTime, diaTimeErr);
    recLight.SetTotalLightAtAperture(diaPhotons, diaPhotonsErr);

    // With the below test, we assume that if there's no fluorescence
    // light trace, there aren't any others either. The reverse
    // is not assumed.
    if (!oRecData.HasLightFlux(fevt::FdConstants::eFluorDirect))
      continue;

    vector<double> dia_fluor(time_size);
    vector<double> dia_fluor_err(time_size);
    vector<double> dia_direct_cher(time_size);
    vector<double> dia_direct_cher_err(time_size);
    vector<double> dia_mie_cher(time_size);
    vector<double> dia_mie_cher_err(time_size);
    vector<double> dia_ray_cher(time_size);
    vector<double> dia_ray_cher_err(time_size);
    vector<double> dia_mult_cher(time_size);
    vector<double> dia_mult_fluor(time_size);

    // Fetch/write direct fluorescence light
    if (oRecData.HasLightFlux(fevt::FdConstants::eFluorDirect)) {
      const TabulatedFunctionErrors& directFluorescence =
        oRecData.GetLightFlux(fevt::FdConstants::eFluorDirect);
      std::copy(directFluorescence.YBegin(),
                directFluorescence.YEnd(),
                dia_fluor.begin() + firstComponentBin);
      std::copy(directFluorescence.YErrBegin(),
                directFluorescence.YErrEnd(),
                dia_fluor_err.begin() + firstComponentBin);
    }

    // Fetch/write direct Cherenkov light
    if (oRecData.HasLightFlux(fevt::FdConstants::eCherDirect)) {
      const TabulatedFunctionErrors& directCherenkov =
        oRecData.GetLightFlux(fevt::FdConstants::eCherDirect);
      std::copy(directCherenkov.YBegin(),
                directCherenkov.YEnd(),
                dia_direct_cher.begin() + firstComponentBin);
      std::copy(directCherenkov.YErrBegin(),
                directCherenkov.YErrEnd(),
                dia_direct_cher_err.begin() + firstComponentBin);
    }

    // Fetch/write Mie scattered Cherenkov light
    if (oRecData.HasLightFlux(fevt::FdConstants::eCherMieScattered)) {
      const TabulatedFunctionErrors& mieCherenkov =
        oRecData.GetLightFlux(fevt::FdConstants::eCherMieScattered);
      std::copy(mieCherenkov.YBegin(),
                mieCherenkov.YEnd(),
                dia_mie_cher.begin() + firstComponentBin);
      std::copy(mieCherenkov.YErrBegin(),
                mieCherenkov.YErrEnd(),
                dia_mie_cher_err.begin() + firstComponentBin);
    }

    // Fetch/write Rayleigh scattered Cherenkov light
    if (oRecData.HasLightFlux(fevt::FdConstants::eCherRayleighScattered)) {
      const TabulatedFunctionErrors& rayCherenkov =
        oRecData.GetLightFlux(fevt::FdConstants::eCherRayleighScattered);
      std::copy(rayCherenkov.YBegin(),
                rayCherenkov.YEnd(),
                dia_ray_cher.begin() + firstComponentBin);
      std::copy(rayCherenkov.YErrBegin(),
                rayCherenkov.YErrEnd(),
                dia_ray_cher_err.begin() + firstComponentBin);
    }

    // Fetch/write multiply scattered Cherenkov light
    if (oRecData.HasLightFlux(fevt::FdConstants::eCherMultScattered)) {
      const TabulatedFunctionErrors& multCherenkov =
        oRecData.GetLightFlux(fevt::FdConstants::eCherMultScattered);
      std::copy(multCherenkov.YBegin(),
                multCherenkov.YEnd(),
                dia_mult_cher.begin() + firstComponentBin);
    }

    // Fetch/write multiply scattered fluorescence light
    if (oRecData.HasLightFlux(fevt::FdConstants::eFluorMultScattered)) {
      const TabulatedFunctionErrors& multFl =
        oRecData.GetLightFlux(fevt::FdConstants::eFluorMultScattered);
      std::copy(multFl.YBegin(),
                multFl.YEnd(),
                dia_mult_fluor.begin() + firstComponentBin);
    }

    for (unsigned int i = 0; i< diaPhotons.size(); ++i) {
      // diaTimeErr has the half bin width
      const double timeBin = 2 * diaTimeErr[i];
      dia_fluor[i]           /= timeBin;
      dia_fluor_err[i]       /= timeBin;
      dia_direct_cher[i]     /= timeBin;
      dia_direct_cher_err[i] /= timeBin;
      dia_mie_cher[i]        /= timeBin;
      dia_mie_cher_err[i]    /= timeBin;
      dia_ray_cher[i]        /= timeBin;
      dia_ray_cher_err[i]    /= timeBin;
      dia_mult_fluor[i]      /= timeBin;
      dia_mult_cher[i]       /= timeBin;
    }

    recLight.SetFluorLightAtAperture(dia_fluor);
    recLight.SetCherLightAtAperture(dia_direct_cher);
    recLight.SetMieCherLightAtAperture(dia_mie_cher);
    recLight.SetRayCherLightAtAperture(dia_ray_cher);
    recLight.SetMultScatCherLightAtAperture(dia_mult_cher);
    recLight.SetMultScatFluorLightAtAperture(dia_mult_fluor);

    // Set the light collection efficiencies for the most prominent light components
    if (Config().StoreLCEfficiency() && oRecData.HasLightCollectionEfficiency()) {
      const utl::MultiTabulatedFunctionErrors& lcEff = oRecData.GetLightCollectionEfficiency();

      if (lcEff.HasLabel(fevt::FdConstants::eFluorDirect)) {
        const TabulatedFunctionErrors& effTab =
          lcEff.GetTabulatedFunctionErrors(fevt::FdConstants::eFluorDirect);
        const vector<Double_t> eff(effTab.YBegin(), effTab.YEnd());
        recLight.SetFluorLightCollEfficiency(eff);
        const vector<Double_t> effErr(effTab.YErrBegin(), effTab.YErrEnd());
        recLight.SetFluorLightCollEfficiencyErrors(effErr);
      }

      if (lcEff.HasLabel(fevt::FdConstants::eCherDirect)) {
        const TabulatedFunctionErrors& effTab =
          lcEff.GetTabulatedFunctionErrors(fevt::FdConstants::eCherDirect);
        const vector<Double_t> eff(effTab.YBegin(), effTab.YEnd());
        recLight.SetCherLightCollEfficiency(eff);
        const vector<Double_t> effErr(effTab.YErrBegin(), effTab.YErrEnd());
        recLight.SetCherLightCollEfficiencyErrors(effErr);
      }

      if (lcEff.HasLabel(fevt::FdConstants::eCherMieScattered)) {
        const TabulatedFunctionErrors& effTab =
          lcEff.GetTabulatedFunctionErrors(fevt::FdConstants::eCherMieScattered);
        const vector<Double_t> eff(effTab.YBegin(), effTab.YEnd());
        recLight.SetMieCherLightCollEfficiency(eff);
        const vector<Double_t> effErr(effTab.YErrBegin(), effTab.YErrEnd());
        recLight.SetMieCherLightCollEfficiencyErrors(effErr);
      }

      if (lcEff.HasLabel(fevt::FdConstants::eCherRayleighScattered)) {
        const TabulatedFunctionErrors& effTab =
          lcEff.GetTabulatedFunctionErrors(fevt::FdConstants::eCherRayleighScattered);
        const vector<Double_t> eff(effTab.YBegin(), effTab.YEnd());
        recLight.SetRayCherLightCollEfficiency(eff);
        const vector<Double_t> effErr(effTab.YErrBegin(), effTab.YErrEnd());
        recLight.SetRayCherLightCollEfficiencyErrors(effErr);
      }
    }

  }
}


bool
FD2ADST::FillEnergy(const fdet::Eye& /*dEye*/, const fevt::Eye& eye, FDEvent& theFd)
{
  FdRecShower& theShower = theFd.GetFdRecShower();
  const evt::ShowerFRecData& showerFRec = eye.GetRecData().GetFRecShower();

  bool retval = false;
  const double enTot = showerFRec.GetTotalEnergy();
  if (isfinite(enTot) &&
      enTot > 0 &&
      isfinite(showerFRec.GetTotalEnergyError()) &&
      isfinite(showerFRec.GetEmEnergy()) &&
      isfinite(showerFRec.GetEmEnergyError())) {
    theShower.SetCalorimetricEnergy(showerFRec.GetEmEnergy());
    theShower.SetCalorimetricEnergyError(showerFRec.GetEmEnergyError());
    theShower.SetEnergy(enTot);
    theShower.SetEnergyError(showerFRec.GetTotalEnergyError());
    if (showerFRec.HasTotalEnergyError(evt::ShowerFRecData::eAtmospheric))
      theShower.SetEnergyAtmError(showerFRec.GetTotalEnergyError(evt::ShowerFRecData::eAtmospheric));
    retval = true;
  }

  return retval;
}


void
FD2ADST::IncreaseFdRecLevel(FDEvent& fdevent, EFdRecLevel reclevel)
{
  if (fdevent.GetRecLevel() < reclevel)
    fdevent.SetRecLevel(reclevel);
}


void
FD2ADST::FillFdRecStations(const evt::Event& event, const fevt::Eye& eye, FDEvent& adstFDEvent)
  const
{
  if (!event.HasSEvent())
    return;
  const sevt::SEvent& sevent = event.GetSEvent(); // this gets used

  if (!eye.HasRecData())
    return;
  const fevt::EyeRecData& eyeRec = eye.GetRecData();

  // the following is RecDataProvider::HasFdReconstruction
  if (!(eyeRec.HasFRecShower() && eye.HasHeader()))
    return;
  const evt::ShowerFRecData& showerfrec = eyeRec.GetFRecShower(); // this gets used
  const fevt::EyeHeader& eyeHeader = eye.GetHeader(); // this gets used once, for the time

  // ----------- preparation ---------------
  // needed for age calculations
  const atm::ProfileResult& depthProfile = det::Detector::GetInstance().GetAtmosphere().EvaluateDepthVsHeight();
  const double atmHeightMin = depthProfile.MinX();
  const double atmHeightMax = depthProfile.MaxX();

  det::Detector& Det = det::Detector::GetInstance();
  const sdet::SDetector& theSDetector = Det.GetSDetector();
  CoordinateSystemPtr EyeCS = Det.GetFDetector().GetEye(eye).GetEyeCoordinateSystem();
  CoordinateSystemPtr CS = Det.GetFDetector().GetEye(eye).GetLocalCoordinateSystem();
  Point eyePos = Det.GetFDetector().GetEye(eye).GetPosition();

  utl::Vector vertical(0,0,1, CS);
  utl::Vector SDP = eyeRec.GetSDP();
  const AxialVector horizontal = Normalized(Cross(SDP, vertical));

  const double Chi0 = eyeRec.GetChiZero();
  const double Rp = eyeRec.GetRp();
  const double T0 = eyeRec.GetTZero();

  const Vector eyeCoreVec = Rp / sin(kPi - Chi0) * horizontal;
  eyeCoreVec.TransformTo(CS) ;

  Transformation rot(Transformation::Rotation(-Chi0, SDP, CS));

  Vector axis(rot*horizontal);  // apply rotation
  if (axis.GetTheta(CS) > (kPi/2))
    axis *= -1;  // reflect if necessary
  axis.Normalize();
  double cZen = CosAngle(axis, vertical);

  // eye timing
  const TimeStamp eyeTriggerTime =
    eyeHeader.GetTimeStamp() - TimeInterval (1e5);

  const_cast<sevt::SEvent&>(sevent).SortStations(sevt::ByDecreasingSignal());

  for (auto sditer = sevent.StationsBegin(); sditer != sevent.StationsEnd(); ++sditer) {

    if (!sditer->HasRecData())
      continue;

    try {
      theSDetector.GetStation(*sditer);
    } catch (AugerException& e) {
      ERROR(e.GetMessage());
      continue;
    }

    const sevt::StationRecData& recdata = sditer->GetRecData();

    if (sditer->IsSilent())
      continue;

    const unsigned int thisStationId = sditer->GetId();

    FdRecStation recStation;

    // ------ set all basic station properties ----
    recStation.SetId(sditer->GetId());

    const double time = recdata.GetSignalStartTime().GetGPSNanoSecond();
    //#warning ---- This does not work for events with GPS second crossing (nsec is not enough!!) ---
    recStation.SetTime(time/100, 0);
    recStation.SetTotalSignal(recdata.GetTotalSignal(), recdata.GetTotalSignalError());

    if (sditer->IsCandidate())
      recStation.SetCandidate();
    if (sditer->IsRejected())
      recStation.SetAccidental();

    // station trigger info
    if (sditer->HasTriggerData()) {

      const sevt::StationTriggerData& trigStation = sditer->GetTriggerData();

      if (trigStation.IsTimeOverThreshold())
        recStation.SetToT();
      else if (trigStation.IsTimeOverThresholdDeconvoluted())
        recStation.SetToTd();
      else if (trigStation.IsMultiplicityOfPositiveSteps())
        recStation.SetMoPS();
      else if (trigStation.IsT2Threshold())
        recStation.SetT2Threshold();
      else if (trigStation.IsT1Threshold())
        recStation.SetT1Threshold();
      else if (trigStation.IsRandom())
        recStation.SetRandom();
      else if (trigStation.IsRDThreshold())
        recStation.SetRDThreshold();

      recStation.SetPLDVersion(trigStation.GetPLDVersion());
      recStation.SetPLDTimeOffset(trigStation.GetPLDTimeOffset());
    }

    if (sditer->HasCalibData())
      recStation.SetCalibrationVersion(sditer->GetCalibData().GetVersion());

    // check if station was used in hybrid fit
    const vector<unsigned short int>& hybridStations = showerfrec.GetStationIds();
    for (unsigned int k = 0; k < hybridStations.size(); ++k) {
      if (hybridStations[k] == thisStationId) {
        recStation.SetHybrid();
        break;
      }
    }

    // --- time at eye ------

    const sdet::Station& currentSt = theSDetector.GetStation(*sditer);
    const Point& SdPos = currentSt.GetPosition();
    const Vector eyeTankVec(SdPos.GetX(CS), SdPos.GetY(CS), SdPos.GetZ(CS), CS) ;
    const Vector coreTankVec = eyeTankVec - eyeCoreVec;

    const double azimuth = eyeTankVec.GetPhi(EyeCS);
    const double elevation = kPi/2 - eyeTankVec.GetTheta(EyeCS);

    recStation.SetAzimuthAndElevation(azimuth, elevation);

    const double tankProj = coreTankVec * axis;
    const Vector tankAxisVec = coreTankVec - (tankProj * axis);
    const Vector tankInAxisVec = eyeCoreVec + tankProj * axis;
    const double rho = tankAxisVec.GetMag();  // perp. distance to axis

    // ---- trigger time at station ---

    TimeStamp tank_time = recdata.GetSignalStartTime();

    const double t_i = ((tank_time + TimeInterval (tankInAxisVec.GetMag()/kSpeedOfLight))  // add travel time to Eye
                        - eyeTriggerTime).GetInterval() / 100;

    // angle along axis in SDP
    double chi_i = Angle(tankInAxisVec,horizontal);
    if (tankInAxisVec.GetZ(CS) < 0)
      chi_i *= -1;
    recStation.SetChi_i(chi_i);
    recStation.SetTimeEye(t_i);

    // time residual to axis fit
    const double beta = Chi0 - kPi/2 - chi_i;
    const double t_fit = (T0 + Rp/kSpeedOfLight * (tan (beta)+1./cos (beta)))/100.;
    const double t_residual = std::abs (t_fit - t_i);

    recStation.SetTimeResidual(t_residual);
    recStation.SetDistanceToAxis (rho);

    // ----- shower age of tank ----------

    recStation.SetAge(-1);
    recStation.SetSlantDepth(-1/g*cm2);

    if (adstFDEvent.GetRecLevel() >= eHasLongitudinalProfile) {

      const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();

      const double TankOverSeaLevel =
        wgs84.PointToLatitudeLongitudeHeight(SdPos).get<2>();

      if (TankOverSeaLevel>atmHeightMin && TankOverSeaLevel<atmHeightMax) {
        const double DepthTank = depthProfile.Y (TankOverSeaLevel);
        const double P_USSTD = 1013.272684 * pow(1 - 2.255771988e-5*TankOverSeaLevel/m, 5.255876);  // hPa

        Point tankInAxis = eyePos + tankInAxisVec;
        const double TankHeightAxis =
          wgs84.PointToLatitudeLongitudeHeight(tankInAxis).get<2>();

        double Xtank = 1013.272684 * pow(1 - 2.255771988e-5*TankHeightAxis/m, 5.255876) /*hPa*/
          * (DepthTank/P_USSTD);
        Xtank /= cZen;

        const double Xmax = adstFDEvent.GetFdRecShower().GetXmax();
        if (Xmax)
          recStation.SetAge(FdRecShower::ShowerAge(Xtank, Xmax));

        recStation.SetSlantDepth(Xtank/g*cm2);
      }
    }

    adstFDEvent.AddStation(recStation);

  }
}


void
FD2ADST::FillRecPixel(const fevt::Eye& eye, FDEvent& theFd)
{
  const fevt::EyeRecData& eyerec = eye.GetRecData();
  FdRecGeometry& theGeometry = theFd.GetFdRecGeometry();
  FdRecShower& theShower = theFd.GetFdRecShower();

  // ------ pixel stuff ------
  int nPixel = 0;
  int nTrigPixel = 0;
  int nRecPixel = 0;
  int nTracePixel = 0;
  int nAxisPixel = 0;

  double meanPixelRMS = 0;
  double meanADCRMS = 0;
  vector<unsigned short int> pixelID;
  vector<unsigned short int> pixelStatus;
  vector<double> pixelTime;
  vector<double> pixelTimeErr;
  vector<double> pixelChi;
  vector<double> pixelCharge;
  vector<int> pixelThreshold;

  vector<FDSaturationStatus> pixelSaturationStatus;

  vector<vector<double>> pixelTrace;
  vector<vector<double>> spotRecPixelTrace;
  vector<map<int,vector<double>>> pixelSimTraces;
  vector<int> pixelSimTraceOffsets;

  vector<int> pulseStart;
  vector<int> pulseStop;
  vector<double> pixelRMS;
  map<int, double> calibConst;
  map<int, int> LookUp;

  UShort_t istat = eBackGroundPix;

  bool hasFdPixel = false;
  bool hasFdTracePixel = false;

  vector<int> evPixelID;
  vector<int> mirrorTimeOffset;
  vector<bool> triggeredPixel;
  int nRawPix = 0;
  int nTrigPix = 0;

  double photonRMSSum = 0;
  double ADCRMSSum = 0;

  double angularTrackLength = 0;

  const double kTraceBinning = 100*ns;

  // FIXME
  // This implies that the number of pixels in each telescope is the same
  // I think this assumption is all over RecDataWriter and ADST, so I won't
  // fix it now since the assumption will hold for a long time still --Steffen

  const int nDetectorPixelsPerTel =
    TelescopeGeometry::GetMaxNumberOfPixelsPerCamera();

  for (auto iTel = eye.TelescopesBegin(ComponentSelector::eHasData); iTel != eye.TelescopesEnd(ComponentSelector::eHasData); ++iTel) {
    for (auto iPixel = iTel->PixelsBegin(ComponentSelector::eHasData); iPixel != iTel->PixelsEnd(ComponentSelector::eHasData); ++iPixel) {
      if (!iPixel->HasRecData())
        continue;
      const int pixID = (iPixel->GetId()-1) +
        (iPixel->GetTelescopeId()-1)*nDetectorPixelsPerTel;

      evPixelID.push_back(pixID);
      LookUp[pixID] = nRawPix;
      calibConst[nRawPix] = otoa::fd::GetCalibrationConstant(*iPixel);

      const double photonRMS = iPixel->GetRecData().GetRMS();
      photonRMSSum += photonRMS;

      // undo the calibration to get the ADC-RMS
      double ADCRMS = 0;
      if (calibConst[nRawPix] != 0)
        ADCRMS = photonRMS / calibConst[nRawPix];
      ADCRMSSum += ADCRMS;

      mirrorTimeOffset.push_back(iTel->GetTimeOffset());

      ++nRawPix;

      int firstBin = iPixel->GetRecData().GetFirstTriggeredTimeBin();
      if (firstBin > -1) {
        triggeredPixel.push_back(true);
        ++nTrigPix;
      } else {
        triggeredPixel.push_back(false);
      }
    }
  }

  // store raw pixels
  if (nRawPix > 0)
    hasFdPixel = true;

  if (nRawPix > 0) {
    nPixel = nRawPix;
    nTrigPixel = nTrigPix;
    pixelID.resize(nRawPix);
    pixelStatus.resize(nRawPix);
    pixelTime.resize(nRawPix);
    pixelTimeErr.resize(nRawPix);
    pixelChi.resize(nRawPix);
    pixelCharge.resize(nRawPix);
    pixelSaturationStatus.resize(nRawPix);
    pixelThreshold.resize(nRawPix);

    meanPixelRMS = photonRMSSum / nRawPix;
    meanADCRMS = ADCRMSSum / nRawPix;

    for (int pix = 0; pix < nRawPix; ++pix) {
      if (!triggeredPixel[pix])
        istat = eBackGroundPix;
      else
        istat = eTriggeredPix;

      pixelID[pix]               = evPixelID[pix];
      pixelStatus[pix]           = istat;
      pixelTime[pix]             = 0;
      pixelTimeErr[pix]          = 0;
      pixelChi[pix]              = 0;
      pixelCharge[pix]           = 0;
      pixelThreshold[pix]        = 0;
      pixelSaturationStatus[pix] = eNoSaturation;

      // calculate geometric chi_i for all pixel, IF time fit exists.
      // Later overwrite with reconstructed chi_i for the relevant pixels
      if (eyerec.GetNSDPPixels() > 0)  {
        const int thisTelId = pixelID[pix] / nDetectorPixelsPerTel + 1;
        const int thisPixelId = pixelID[pix] - (thisTelId-1)*nDetectorPixelsPerTel + 1;
        const fdet::Eye& eyeDet = det::Detector::GetInstance().GetFDetector().GetEye(eye.GetId());
        const Vector& pixeldir = eyeDet.GetTelescope(thisTelId).GetPixel(thisPixelId).GetDirection();
        const CoordinateSystemPtr eyeCS = eyeDet.GetEyeCoordinateSystem();
        const Vector vertical(0, 0, 1, eyeCS);
        const utl::Vector SDP = eyerec.GetSDP();
        const Vector horizontalWithinSDP = Normalized(Cross(SDP, vertical));
        // subtract component parallel to sdp
        const Vector chi_i_vector = Normalized(pixeldir - (SDP * pixeldir) * SDP);
        //  Estimating chi_i
        const double chi_i = Angle(chi_i_vector, horizontalWithinSDP);
        pixelChi[pix] = chi_i;
      }
    }
  } else {
    cerr << " RecDataWriter::ProvideFdRecPixel - no triggered pixels???" << endl;
    return;
  }

  // retrieve reconstructed (pulse found) pixels
  if (eyerec.GetNPulsedPixels() > 0) {
    istat = ePulseRecPix;

    int NRecPix = 0;
    for (auto pixiter = eyerec.PulsedPixelsBegin(); pixiter != eyerec.PulsedPixelsEnd(); ++pixiter) {
      const fevt::Pixel& evtpixel = *pixiter;
      int pixNum = evtpixel.GetId();
      int miID = evtpixel.GetTelescopeId();

      int pixID = (miID-1)*nDetectorPixelsPerTel+pixNum-1;

      if (LookUp.count(pixID) == 0) {
        ostringstream err;
        err << " RecDataWriter::ProvideFdRecPixel - empty lookup table slot"
            << pixID << endl;
        throw OutOfBoundException(err.str());
      }

      int pixSerial = LookUp[pixID];

      const PixelRecData& pixrec = evtpixel.GetRecData();
      double time = 0;
      double time_err = 0;
      if (eyerec.GetNTimeFitPixels() > 0) {
        time = pixrec.GetT_i().GetNanoSecond() / kTraceBinning;
        time_err = pixrec.GetT_iError().GetNanoSecond() / kTraceBinning;
      } else {
        Double_t TelescopeTimeOffset = mirrorTimeOffset[pixSerial] / kTraceBinning;
        time = pixrec.GetCentroid() + TelescopeTimeOffset;
        time_err = pixrec.GetCentroidError();
      }

      const double charge = pixrec.GetTotalCharge();
      const double chi = pixrec.GetChi_i();

      int threshold = 0;
      if (evtpixel.HasTriggerData()) {
        const PixelTriggerData& pixtrig = evtpixel.GetTriggerData();
        threshold = pixtrig.GetThreshold();
      }

      pixelStatus[pixSerial] = istat;
      pixelTime[pixSerial] = time;
      pixelTimeErr[pixSerial] = time_err;
      pixelCharge[pixSerial] = charge;
      pixelChi[pixSerial] = chi;
      pixelThreshold[pixSerial] = threshold;

      if (evtpixel.IsHighGainSaturation())
        pixelSaturationStatus[pixSerial] = eHighGainSat;
      else if (evtpixel.IsLowGainSaturation())
        pixelSaturationStatus[pixSerial] = eLowGainSat;
      else if (evtpixel.IsSaturationRecovered())
        pixelSaturationStatus[pixSerial] = eRecovered;

      ++NRecPix;

    }

    if (!NRecPix)
      return;

    nRecPixel = NRecPix;
  }

  // mark pixels used in SDP fit and calculate RMS to SDP
  if (eyerec.GetNSDPPixels() > 0)  {

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

    const Vector sdp = eyerec.GetSDP();

    istat = eSDPRecPix;
    int NSDPPix = 0;
    double RMS = 0;

    for (auto pixiter = eyerec.SDPPixelsBegin(); pixiter != eyerec.SDPPixelsEnd(); ++pixiter) {

      const fevt::Pixel& evtpixel = *pixiter;
      const fdet::Pixel& detpixel =
        det::Detector::GetInstance().GetFDetector().GetPixel(evtpixel);

      const Vector dir = detpixel.GetDirection();
      dir.TransformTo(eyeCS);

      const Vector SDP = eyerec.GetSDP();;

      const double d = kPi/2 - Angle(dir, SDP);
      RMS += (d*d);

      const int pixNum = evtpixel.GetId();
      const int miID = evtpixel.GetTelescopeId();

      const int pixID = (miID-1) * nDetectorPixelsPerTel+pixNum-1;

      if (LookUp.count(pixID) == 0) {
        cerr << " RecDataWriter::FillFdRecPixel - empty lookup table slot" << endl;
        exit(1);
      }
      const int pixSerial = LookUp[pixID];

      pixelStatus[pixSerial] = istat;
      ++NSDPPix;
    }

    RMS = (NSDPPix > 0 ? sqrt(RMS/(double)NSDPPix) : 0.);
  }

  // mark pixels used in time fit
  if (eyerec.GetNTimeFitPixels() > 0)  {

    bool start = true;
    double MinChi = 0;
    double MaxChi = 0;

    int NAxisPix = 0;
    istat = eTimeFitPix;
    for (auto pixiter = eyerec.TimeFitPixelsBegin(); pixiter != eyerec.TimeFitPixelsEnd(); ++pixiter) {
      const fevt::Pixel& evtpixel = *pixiter;
      const int pixNum = evtpixel.GetId();
      const int miID = evtpixel.GetTelescopeId();

      const int pixID = (miID-1) * nDetectorPixelsPerTel+pixNum-1;

      if (LookUp.count(pixID) == 0) {
        cerr << " RecDataWriter::FillFdRecPixel - empty lookup table slot" << endl;
        exit(1);
      }
      const int pixSerial = LookUp[pixID];

      if (start || pixelChi[pixSerial] > MaxChi || start)
        MaxChi = pixelChi[pixSerial];
      if (start || pixelChi[pixSerial] < MinChi || start)
        MinChi = pixelChi[pixSerial];

      start = false;

      pixelStatus[pixSerial] = istat;

      ++NAxisPix;
    }

    nAxisPixel = NAxisPix;
    if (nAxisPixel > 0) {
      angularTrackLength = fabs(MaxChi - MinChi);
      theShower.SetAngTrackObs(angularTrackLength);
    }
  }

  // Pixel traces
  if (Config().StoreFDTraces() > 0) {

    nTracePixel = nPixel;
    const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();

    int tracePixel = 0;

    pulseStart.resize(nTracePixel);
    pulseStop.resize(nTracePixel);
    pixelRMS.resize(nTracePixel);
    pixelTrace.resize(nTracePixel);
    spotRecPixelTrace.resize(nTracePixel);
    pixelSimTraces.resize(nTracePixel);
    pixelSimTraceOffsets.resize(nTracePixel);

    for (auto iTel = eye.TelescopesBegin(ComponentSelector::eHasData); iTel != eye.TelescopesEnd(ComponentSelector::eHasData); ++iTel) {
      const fdet::Camera& detCamera = detFD.GetTelescope(*iTel).GetCamera();

      for (auto iPixel = iTel->PixelsBegin(ComponentSelector::eHasData); iPixel != iTel->PixelsEnd(ComponentSelector::eHasData); ++iPixel) {
        const fevt::Pixel& traceEvtPixel = *iPixel;

        if (!traceEvtPixel.HasRecData())
          continue;

        const int tracePixNo = traceEvtPixel.GetId();
        const int traceMiID = traceEvtPixel.GetTelescopeId();

        const int tracePixId = (traceMiID-1)*nDetectorPixelsPerTel+tracePixNo-1;

        if (LookUp.count(tracePixId) == 0) {
          ostringstream err;
          err << " RecDataProvider::FillFdRecPixel - empty lookup table slot" << tracePixId << endl;
          throw OutOfBoundException(err.str());
        }

        const int tracePixSerial = LookUp.find(tracePixId)->second;

        const PixelRecData& tracePixRec = traceEvtPixel.GetRecData();

        pixelTrace[tracePixSerial].assign(tracePixRec.GetPhotonTrace().Begin(), tracePixRec.GetPhotonTrace().End());

        //spotRecPixelTrace[tracePixSerial].resize(size);
        // to be uncommented as soon as SpotReconstruction becomes available
        //if (tracePixRec.HasPhotonTrace(fevt::FdConstants::eReconstructed))
        //  spotRecPixelTrace[tracePixSerial].assign(tracePixRec.GetPhotonTrace(fevt::FdConstants::eReconstructed).Begin(), tracePixRec.GetPhotonTrace(fevt::FdConstants::eReconstructed).End());
        //  ++spotRecTracePixel;
        //}

        if (Config().StoreMCTraces() && traceEvtPixel.HasSimData()) {
          for (auto iSource = traceEvtPixel.GetSimData().PhotonTracesBegin();
               iSource != traceEvtPixel.GetSimData().PhotonTracesEnd(); ++iSource) {

            FdApertureLight::ELightComponent sourceLabel = FdApertureLight::eUnknown;
            switch ((FdConstants::LightSource)iSource->GetLabel()) {
            //case eTotal:
              //eFluorTotal,
            case FdConstants::eFluorDirect:
              sourceLabel = FdApertureLight::eFluorescence;
              break;
              //eFluorRayleighScattered, // Rayleigh scattered
              //eFluorMieScattered,      // Mie scattered
              //eFluorMultScattered,     // Combination of Rayleigh and
              //                         // Mie multiple scattered
            case FdConstants::eCherDirect:
              {
                // ostringstream inf;
                // inf << "fd2adst cherdirect " << tracePixNo << " " << traceMiID << " " << tracePixSerial;
                // INFO(inf);
              }
              sourceLabel = FdApertureLight::eDirectCherenkov;
              break;
            case FdConstants::eCherRayleighScattered:
              sourceLabel = FdApertureLight::eRayleighCherenkov;
              break;
            case FdConstants::eCherMieScattered:
              sourceLabel = FdApertureLight::eMieCherenkov;
              break;
              //ePrimaryCherDirect,
              //ePrimaryCherRayleighScattered,
              //ePrimaryCherMieScattered,
              //case eLaserRayleighScattered:
              //case eLaserMieScattered:
              //            eBackground:             // from the starry night
              //eReconstructed,          // result of spot reconstruction
              //eCherMultScattered,
              //eUserDefined1 = 32,
              //eUserDefined2,
              //eLastSource = eUserDefined2
            default:
              continue;
              break;
            }

            pixelSimTraces[tracePixSerial][sourceLabel].assign(iSource->GetTrace().Begin(), iSource->GetTrace().End());
            pixelSimTraceOffsets[tracePixSerial] = detCamera.GetSLTTriggerBin() - iTel->GetSimData().GetSltTimeShift();
          }
        }

        pulseStart[tracePixSerial] = tracePixRec.GetPulseStart();
        pulseStop[tracePixSerial] = tracePixRec.GetPulseStop();
        pixelRMS[tracePixSerial] = tracePixRec.GetRMS();
        ++tracePixel;

      }
    }
    nTracePixel = tracePixel;

    if (nTracePixel > 0)
      hasFdTracePixel = true;

  }

  // ########################  SAVE PIXELS ############################

  if (hasFdPixel) {

    FdRecPixel& thePixel = theFd.GetFdRecPixel();

    if (Config().StoreAllPixels())
      thePixel.SetNumberOfPixels(nPixel);
    else
      thePixel.SetNumberOfPixels(nTrigPixel);

    thePixel.SetNumberOfTriggeredPixels(nTrigPixel);
    thePixel.SetNumberOfPulsedPixels(nRecPixel);
    if (Config().StoreFDTraces() > 0 && hasFdTracePixel)
      thePixel.SetNumberOfTracePixels(nTracePixel);
    else
      thePixel.SetNumberOfTracePixels(0);

    thePixel.SetMeanPixelRMS(meanPixelRMS);
    thePixel.SetMeanADCRMS(meanADCRMS);

    int pixSerial = 0;
    for (int pix = 0; pix < nPixel; ++pix) {

      if (Config().StoreAllPixels() || (EPixelStatus) pixelStatus[pix] >= eTriggeredPix) {
        thePixel.SetID(pixSerial, pixelID[pix]);
        thePixel.SetStatus(pixSerial, (EPixelStatus)pixelStatus[pix]);
        thePixel.SetTime(pixSerial, pixelTime[pix]);
        thePixel.SetTimeErr(pixSerial, pixelTimeErr[pix]);
        thePixel.SetCharge(pixSerial, pixelCharge[pix]);
        thePixel.SetChi(pixSerial, pixelChi[pix]);
        thePixel.SetThreshold (pixSerial, pixelThreshold[pix]);
        if (pixelSaturationStatus[pix] == eHighGainSat)
          thePixel.SetHighGainSaturated(pixSerial);
        else if (pixelSaturationStatus[pix] == eLowGainSat)
          thePixel.SetLowGainSaturated(pixSerial);
        else if (pixelSaturationStatus[pix] == eRecovered)
          thePixel.SetSaturationRecovered(pixSerial);

        if (Config().StoreFDTraces() > 0 && hasFdTracePixel) {
          thePixel.SetDataTrace(pixSerial, pixelTrace[pix]);
          thePixel.SetPulseStart(pixSerial, pulseStart[pix]);
          thePixel.SetPulseStop(pixSerial, pulseStop[pix]);
          thePixel.SetRMS(pixSerial, pixelRMS[pix]);
          thePixel.SetCalibrationConstant(pixSerial, calibConst[pix]);
          if (Config().StoreMCTraces()) {
            for (auto iSource = pixelSimTraces[pix].begin(); iSource != pixelSimTraces[pix].end(); ++iSource) {
              thePixel.SetSimTrace(iSource->second, pixSerial, (FdApertureLight::ELightComponent)iSource->first, pixelSimTraceOffsets[pix]);
            }
          }
        }

        ++pixSerial;
      }
    }
  }

  if (!pixelChi.empty()) {
    double a1;
    double a2;
    double chi2;
    double linearTimeFitChi2 = 0;
    otoa::LinearFit(pixelChi, pixelTime, pixelTimeErr, a1, a2, chi2);

    if (!TMath::IsNaN(chi2))
      linearTimeFitChi2 = chi2;

    theGeometry.SetLinearTimeFitChi2(linearTimeFitChi2);
  }
}


void
FD2ADST::FillEyeSim(const evt::Event& event, FDEvent& fd)
{
  if (!event.HasSimShower() || !event.GetSimShower().HasGeometry())
    return;

  const int eyeId = fd.GetEyeId();

  const evt::ShowerSimData& sim = event.GetSimShower();

  det::Detector& detector = det::Detector::GetInstance();
  const fdet::FDetector& detFd = detector.GetFDetector();
  const fdet::Eye& detEye = detFd.GetEye(eyeId);
  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const Point& eyePos = detEye.GetPosition();

  // ----------------  MC geometry -------------------
  FdGenGeometry& fdGeometry = fd.GetGenGeometry();
  TimeStamp eyeTriggerTime;

  bool hasEyeTrigger = false;
  if (event.GetFEvent().HasEye(eyeId, ComponentSelector::eInDAQ)) {
    hasEyeTrigger = true;
    const fevt::Eye& eye = event.GetFEvent().GetEye(eyeId, ComponentSelector::eInDAQ);

    bool hasEyeTrigger = eye.HasHeader();
    if (hasEyeTrigger) {
      eyeTriggerTime = eye.GetHeader().GetTimeStamp() - utl::TimeInterval (1000.*100.*ns);
    } else {
      hasEyeTrigger = false;
    }
  }

  const Point& showerCore =  sim.GetPosition();
  const TimeStamp& timeAtCore = sim.GetTimeStamp();
  const Vector& showerAxisReal = sim.GetDirection(); // REAL SHOWER MOVEMENT
  Vector showerAxis = showerAxisReal;
  bool upward = true;
  if (showerAxis.GetTheta(eyeCS) > 90*deg) {
    showerAxis *= -1;  // reco-def (ALWAYS POINTING UPWARDS)
    upward = false;
  }

  const Vector eyeToCore = showerCore - eyePos;
  Vector sdp = Cross(showerAxis, eyeToCore);
  const double rp = sdp.GetMag() * (upward ? -1 : 1);
  sdp.Normalize();

  const Vector vertical(0, 0, 1, eyeCS);
  const utl::Vector horizontalInSDP = Normalized(Cross(sdp, vertical));

  const double chi0 = Angle(horizontalInSDP, showerAxis) - (upward ? 180*deg : 0);

  const double distCoreEye = eyeToCore.GetMag();
  const double beta = Angle(-showerAxisReal, eyeToCore);
  const double T0 = (hasEyeTrigger ? ((timeAtCore + TimeInterval(distCoreEye/kSpeedOfLight * cos(beta)))
                                      - eyeTriggerTime).GetInterval()
                     : distCoreEye/kSpeedOfLight * cos(beta));
  fdGeometry.SetSDP(sdp.GetTheta(eyeCS), sdp.GetPhi(eyeCS));
  fdGeometry.SetAxis(T0, rp, chi0);

  // --- shower profile parameters seen by eye -----
  FdGenShower& fdShower = fd.GetGenShower();

  if (sim.HasGHParameters()) {
    const double XmaxSim = sim.GetGHParameters().GetXMax();
    fdShower.SetXmaxChi(otoa::fd::ConvertXToChi(XmaxSim, detEye,
                                                showerCore, showerAxis, sdp));
  }
  fdShower.SetX1Chi(otoa::fd::ConvertXToChi(sim.GetXFirst(), detEye,
                                            showerCore, showerAxis, sdp));
}


void
FD2ADST::FillHybridStations(const evt::Event& event, const fevt::Eye& eye, FdRecGeometry& theGeometry)
{
  //TODO: one could pass only eyeRec, but then this check has to be made before.
  const fevt::EyeRecData& eyeRec = eye.GetRecData();

  if (!(eyeRec.HasFRecShower() && eye.HasHeader()))
    return;

  if (!event.HasSEvent())
    return;

  const SEvent& sevent = event.GetSEvent();

  det::Detector& Det = det::Detector::GetInstance();
  const sdet::SDetector& theSDetector = Det.GetSDetector();
  CoordinateSystemPtr EyeCS = Det.GetFDetector().GetEye(eye).GetEyeCoordinateSystem();
  CoordinateSystemPtr CS = Det.GetFDetector().GetEye(eye).GetLocalCoordinateSystem();
  Point eyePos = Det.GetFDetector().GetEye(eye).GetPosition();

  utl::Vector vertical(0,0,1, CS);
  utl::Vector SDP = eyeRec.GetSDP();
  const AxialVector horizontal = Normalized(Cross(SDP, vertical));

  double Chi0 = eyeRec.GetChiZero();
  double Rp = eyeRec.GetRp();

  Vector eyeCoreVec = Rp / sin(kPi - Chi0) * horizontal;
  eyeCoreVec.TransformTo(CS) ;

  Transformation rot (Transformation::Rotation(-Chi0, SDP, CS));

  Vector axis(rot*horizontal);                // apply rotation
  if (axis.GetTheta(CS) > (kPi/2))
    axis *= -1; // reflect if necessary
  axis.Normalize();

  bool first = true;
  int hottestId = 0;
  double hottestSignal = 0;
  double hottestRho = 0;
  double hottestSDPdist = 0;
  for (auto sditer = sevent.StationsBegin(); sditer != sevent.StationsEnd(); ++sditer) {

    if (!sditer->HasRecData())
      continue;

    const StationRecData& recdata = sditer->GetRecData();

    if (sditer->IsSilent())
      continue;

    unsigned int thisStationId = sditer->GetId();
    double signal = recdata.GetTotalSignal();

    const sdet::Station& currentSt = theSDetector.GetStation(*sditer);
    const Point& SdPos = currentSt.GetPosition();
    Vector eyeTankVec(SdPos.GetX(CS), SdPos.GetY(CS), SdPos.GetZ(CS), CS);
    Vector coreTankVec = eyeTankVec - eyeCoreVec;

    double tankProj = (coreTankVec * axis);
    Vector tankAxisVec = coreTankVec - (tankProj * axis);
    Vector tankInAxisVec = eyeCoreVec + tankProj * axis;
    double rho = tankAxisVec.GetMag();  // perp. distance to axis

    // x_0 = (x0,y0,z0) point
    // D=n^ . x_0 + p AND  p=0 (origin at eye n^ is SDP at eye -> plane in Hessian form)
    double sdp_dist = SDP * eyeTankVec;
    sdp_dist = std::abs(sdp_dist);

    // check if station was used in hybrid fit
    const vector<unsigned short int>& hybridStations = eyeRec.GetFRecShower().GetStationIds();
    for (unsigned int k = 0; k < hybridStations.size(); ++k) {
      if (hybridStations[k]==thisStationId) {
        if (first) {
          first = false;
          hottestId = thisStationId;
          hottestSignal = signal;
          hottestSDPdist = sdp_dist;
          hottestRho = rho;
        } else {
          if (signal>hottestSignal) {
            hottestSignal = signal;
            hottestRho = rho;
            hottestSDPdist = sdp_dist;
            hottestId = thisStationId;
          }
        }
        break;
      }
    }
  }

  const unsigned int nStations = eyeRec.GetFRecShower().GetStationIds().size();
  if (nStations > 0) {
    theGeometry.SetNHybridStations(nStations);
    theGeometry.SetHottestStation(hottestId);
    theGeometry.SetStationAxisDistance(hottestRho);
    theGeometry.SetSDFDTimeOffset(eyeRec.GetFRecShower().GetSDTimeResidual());
    theGeometry.SetStationSDPDistance(hottestSDPdist);
  }
}


void
FD2ADST::FillViewingAngle(FdApertureLight* apLight, FdGeometry* geo, const double XmaxChi)
{
  // this will calculate the min and mean viewing angles for profile
  // #warning FillViewingAngle must be made aware of that the time-at-diaphragm data might only be available per telescope.

  const double rp = geo->GetRp();
  const double chi0 = geo->GetChi0();
  const double t0 = geo->GetT0() / apLight->GetTimeBinning();

  apLight->SetXmaxAngle(chi0 - XmaxChi);

  const vector<double>& time = apLight->GetTime();  // in units of time bins
  if (time.empty())  // make extra sure...
    return;

  // convert times to angle towards shower axis
  double minang = 0;
  double maxang = 0;
  double meanang = 0;

  const unsigned int nPoints = time.size();

  if (nPoints > 1) {
    const double c = 0.299792458 * apLight->GetTimeBinning(); // speed of light (m pro time-bin)

    bool first = true;
    for (unsigned int i = 0; i < nPoints; ++i) {
      const double B = (time[i]-t0) * c/rp;
      const double chii = chi0 - atan(B)*2;
      const double angl = chi0 - chii;

      if (first) {
        first = false;
        minang = angl;
        maxang = angl;
      } else {
        if (angl < minang)
          minang = angl;
        if (angl > maxang)
          maxang = angl;
      }
      meanang += angl;
    }

    meanang /= (double)nPoints;
  }

  // FdRecShower& recShower = theFd.GetFdRecShower();
  apLight->SetMinAngle(minang);
  apLight->SetMaxAngle(maxang);
  apLight->SetMeanAngle(meanang);
}


void
FD2ADST::FillSDEye(const evt::Event& theEvent, const fevt::Eye& theEye, FDEvent& theFd)
{
  SdRecGeometry& sdGeom = theFd.GetSdRecGeometry();

  if (theEvent.HasRecShower() &&
      theEvent.GetRecShower().HasSRecShower()) {
    const evt::ShowerSRecData& sdRecShower =
      theEvent.GetRecShower().GetSRecShower();

    const Point& showerCore = sdRecShower.GetCorePosition();
    const Vector& showerAxis = sdRecShower.GetAxis();

    // "SDP" as seen by SD
    utl::Point eyePos = det::Detector::GetInstance().GetFDetector().GetEye(theEye).GetPosition();

    utl::Vector eyeToShower = showerCore - eyePos;

    utl::Vector sdp = Normalized(Cross(showerAxis, eyeToShower));

    CoordinateSystemPtr eyeCS = det::Detector::GetInstance().GetFDetector().GetEye(theEye).GetEyeCoordinateSystem();
    sdGeom.SetSDP(sdp.GetTheta(eyeCS),sdp.GetPhi(eyeCS));

    // T0, Chi0 and Rp as seen by SD

    TimeStamp timeAtCore = sdRecShower.GetCoreTime();

    if (timeAtCore == TimeStamp())
      return;

    utl::TimeStamp eyeTriggerTime = theEye.GetHeader().GetTimeStamp() - utl::TimeInterval(1000*100*ns);
    const Vector vertical(0,0,1, eyeCS);
    const Vector horizontalInSDP = Normalized(Cross(sdp, vertical));  // horizontal

    double chi0 = Angle(horizontalInSDP, showerAxis);
    double Rp = Cross(showerAxis, (showerCore - eyePos)).GetMag() / showerAxis.GetMag();
    utl::Vector vecCoreEye = showerCore - eyePos;
    double distCoreEye = vecCoreEye.GetMag();
    double beta = Angle(-showerAxis, vecCoreEye);
    double T0 = ((timeAtCore - TimeInterval(distCoreEye/kSpeedOfLight * cos (beta))) - eyeTriggerTime).GetInterval();

    sdGeom.SetAxis(T0,Rp,chi0);
  }
}