FdCalibrator.cc

Go to the documentation of this file.

#include "FdCalibrator.h"

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

#include <utl/config.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>

#include <evt/Event.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Header.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeTriggerData.h>
#include <fevt/Telescope.h>
#include <fevt/Channel.h>
#include <fevt/ChannelRecData.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <fevt/PixelTriggerData.h>
#include <fevt/TelescopeTriggerData.h>
#include <fevt/SLTData.h>

#include <det/Detector.h>

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

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>

#include <AugerEvent.h>
#include <EyeEvent.hh>
#include <EyeEventHeader.hh>
#include <EyeT3Data.hh>
#include <FadcData.hh>
#include <FdNumbering.hh>
#include <FDEventLibVersion.hh>

#include <sstream>
#include <cmath>
#include <iomanip>

using namespace FdCalibratorOG;
using namespace fevt;
using namespace fdet;
using namespace utl;
using namespace fwk;
using namespace std;


VModule::ResultFlag
FdCalibrator::Init()
{
  Branch topB = CentralConfig::GetInstance()->GetTopBranch("FdCalibrator");

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

  topB.GetChild("verbosity").GetData(fVerbosity);
  topB.GetChild("correctUncorrectedTimes").GetData(fCorrectTimes);
  topB.GetChild("firstMeanSlot").GetData(fFirstMeanSlot);
  topB.GetChild("lastMeanSlot").GetData(fLastMeanSlot);
  topB.GetChild("recoverSaturatedTraces").GetData(fRecoverSaturatedTraces);
  topB.GetChild("calibrationCorrection").GetData(fCalibrationCorrection);
  topB.GetChild("heatCalibrationEpochs").GetData(fHeatCaliEpochs);
  topB.GetChild("heatCalibrationCorrection").GetData(fHeatCalibration);

  if (fCalibrationCorrection.size() != 5) {
    ERROR( "Calibration Correction must have size of 5!");
    return eFailure;
  }

  // info output
  ostringstream info;
  info << "\n"
          " Version: " << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "         correct times: " << fCorrectTimes << "\n"
          "       first mean slot: " << fFirstMeanSlot << "\n"
          "        last mean slot: " << fLastMeanSlot << "\n"
          "    saturation recovery is " << (fRecoverSaturatedTraces ? "on" : "off") << "\n"
          "Calibration correction: ";
  for (unsigned int i = 0; i < fCalibrationCorrection.size(); ++i)
    info << fCalibrationCorrection[i] << ' ';
  INFO(info);

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

    if (fCalibrationCorrection[i] != 1 && (i+1) != 5) {
      ostringstream warn;
      warn << "Calibration correction of " << fCalibrationCorrection[i] << " "
              "is applied to eye: " << i+1 << ". "
              "ARE YOU SURE??";
      WARNING(warn);
    }

    if ((i+1) == 5 && fCalibrationCorrection[i] == 0)
      INFO("HEAT calbration correction will be esimated using a parametrization");
  }
  return eSuccess;
}


VModule::ResultFlag
FdCalibrator::Run(evt::Event& event)
{
  if (!event.HasRawEvent())
    return eContinueLoop;

  AugerEvent& rawEvent = event.GetRawEvent();

  if (rawEvent.EyesBegin() != rawEvent.EyesEnd())
    ++RunController::GetInstance().GetRunData().GetNamedCounters()["FdCalibrator/FdEvents"];

  for (AugerEvent::EyeIterator eyeIter = rawEvent.EyesBegin();
       eyeIter != rawEvent.EyesEnd(); ++eyeIter) {

    TEyeEvent& eyeEvent = *eyeIter;
    const unsigned int eyeNo = eyeEvent.GetEventHeader()->GetEyeNo();
#if FDEVENTLIB_VERSION_CODE < ModuleVersionCode(4, 0, 0)
    //const unsigned int eyeNo = eyeEvent.GetEventHeader()->GetEyeNo();
    if (eyeNo > 4) {
      ostringstream msg;
      msg << "Trying to read data with eye ids > 4 (HEAT) using an "
             "FDEventLib older than v4r0. Old FDEventLibs do not support "
             "HEAT data. I will skip this eye (id=" << eyeNo << ").";
      WARNING(msg);
      continue;
    }
#endif

    if (eyeNo == 5 && fCalibrationCorrection[eyeNo-1] == 0) {
      const unsigned int gpsSec =
        eyeEvent.GetEventHeader()->GetTimeStamp()->GetGPSSec();
      fCalibrationCorrection[eyeNo-1] = GetHeatCalibrationCorrection(gpsSec);
    }

    if (ApplyTimeCorrections(event, eyeEvent)) {
      FillTriggerData(event, eyeEvent);
      FillChannelRecData(event, eyeEvent);
      FillCalibratedPixels(event, eyeEvent);
      FillDAQInformation(event, eyeEvent);
    }

  }
  return eSuccess;
}


VModule::ResultFlag
FdCalibrator::Finish()
{
  return eSuccess;
}


void
FdCalibrator::FillTriggerData(evt::Event& event, TEyeEvent& eyeEvent)
{
  // first the eye trigger data ...

  TEyeT3Data* const t3Data = eyeEvent.GetT3Data();
  TEyeEventHeader* const eyeHeader =  eyeEvent.GetEventHeader();
  const int eyeId = eyeEvent.GetEventHeader()->GetEyeNo();
  FEvent& fevent = event.GetFEvent();
  fevt::Eye& eye = fevent.GetEye(eyeId);

  if (t3Data) {

    if (!eye.HasTriggerData())
      eye.MakeTriggerData();
    fevt::EyeTriggerData& eyeTrig = eye.GetTriggerData();

    eyeTrig.SetT3Accepted(true);

    const TEyeEventHeader::EEventClass eventClass =
      eyeHeader->GetEventClass();
    const string eventClassString =
      string(TEyeEventHeader::GetVerboseEventClass(eventClass));
    eyeTrig.SetT3Class(eventClassString);

    try {

      const TimeStamp t3TimeStamp(int(fevent.GetHeader().GetTime().GetGPSSecond()),
                                  t3Data->GetTimeAtGround());
      eyeTrig.SetT3Time(t3TimeStamp);

    } catch (OutOfBoundException& e) {

      WARNING(e.GetMessage());
      WARNING("T3 time at ground seems to have an invalid value. "
              "Setting it to zero.");
      eyeTrig.SetT3Time(0);

    }
    eyeTrig.SetT3SDP(t3Data->GetSDPTheta(),
                     t3Data->GetSDPPhi(),
                     t3Data->GetAzimuthAtGround());
    eyeTrig.SetT3NPixels(t3Data->GetNPixels());
  }

  // ... then the mirror trigger data ...

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeId);
  TEyePixelData* eyePixelData = eyeEvent.GetPixelData();

  for (fdet::Eye::TelescopeIterator telIter = detEye.TelescopesBegin();
       telIter != detEye.TelescopesEnd(); ++telIter) {

    const unsigned int mirrorId = telIter->GetId();

    if (eyeHeader->IsMirrorDataPresent(mirrorId)) {

      if (!eye.HasTelescope(mirrorId))
        eye.MakeTelescope(mirrorId);

      fevt::Telescope& tel = eye.GetTelescope(mirrorId);

      if (fVerbosity) {
        ostringstream info;
        info << "Fill trigger-info for mirror: " << mirrorId;
        INFO(info);
      }

      if (tel.HasTriggerData()) {
        if (fVerbosity) {
          ostringstream warn;
          warn << " Mirror " << mirrorId << " has already trigger-data. Skipping.";
          WARNING(warn);
        }
        continue;
      }

      tel.MakeTriggerData();
      fevt::TelescopeTriggerData& telTrigger = tel.GetTriggerData();

      telTrigger.SetTLTAccepted(true);

      const fdet::Camera& theCam = telIter->GetCamera();
      const unsigned int sltNumberOfBins = (unsigned int)(theCam.GetSLTTraceLength());

      // fill telescope trigger data
      vector<fevt::SLTData> sltData;
      TMirrorPixelData* const miPixelData = eyePixelData->GetPixelData(mirrorId);

      for (unsigned int t2bin = 0; t2bin < sltNumberOfBins; ++t2bin) {
        fevt::SLTData sltDataCol(20);
        for (unsigned int col = 1; col <= 20; ++col) {
          TMirrorPixelData::PixelDataWord word = miPixelData->GetPixelData(col, t2bin);
          sltDataCol.SetSLTDataWord(col, PixelDataGetWord(word));
        }
        sltData.push_back(sltDataCol);
      }
      telTrigger.SetSLTData(sltData);

      // Filling of multiplicity
      // TODO: Multiplicity is always 1000 (also for HEAT!). Possibly move to FModelConfig?
      const unsigned int kNMultiplicityBins = 1000;
      const double kMultiplicityBinSize     = 100*ns;
      telTrigger.MakeMultiplicity(kNMultiplicityBins, kMultiplicityBinSize);
      TraceI& multiplicityTrace = telTrigger.GetMultiplicity();
      TMirrorPixelData::MultiplicityDataWord multiplicityData =
        miPixelData->GetMultiplicityData();

      ostringstream info;
      if (fVerbosity)
        info << "TLT multiplicity: ";
      for (unsigned int tbin = 0; tbin < kNMultiplicityBins; ++tbin) {
        info << multiplicityData[tbin] << ", ";
        multiplicityTrace[tbin] = multiplicityData[tbin];
      }
      if (fVerbosity)
        INFO(info);

    }
  }
}


void
FdCalibrator::FillChannelRecData(evt::Event& event, TEyeEvent& eyeEvent)
{
  TEyePixelList* const eyePixelList = eyeEvent.GetPixelList();
  TEyeFADCData* const eyeFADCData = eyeEvent.GetFADCData();
  const unsigned int npixels = eyePixelList->GetNumPixels();

  const int eyeId = eyeEvent.GetEventHeader()->GetEyeNo();
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeId);

  FEvent& fdEvent = event.GetFEvent();
  fevt::Eye& eye = fdEvent.GetEye(eyeId);

  for (unsigned int iPixel = 0; iPixel < npixels; ++iPixel) {

    if (fVerbosity) {
      ostringstream info;
      info << "copying trace for pixel: " << iPixel;
      INFO(info);
    }

    FdUtil::Fd::FdPixelNumber pixelNumber = eyePixelList->GetPixel(iPixel);

    const unsigned int mirrorId = FdUtil::Fd::GetEyeMirrorNo(pixelNumber);
    const unsigned int mirrorChannelId = FdUtil::Fd::GetMirrorPixelNo(pixelNumber);

    if (!eye.HasTelescope(mirrorId))
      eye.MakeTelescope(mirrorId);
    fevt::Telescope& tel = eye.GetTelescope(mirrorId);

    // copy fadc trace
    const TFADCData* const fadcData = eyeFADCData->GetFADCData(pixelNumber);

    // Protect yourself against the stupid non-existent traces...
    if (!fadcData)
      continue;

    TFADCData::FADCDataWord fadcword = fadcData->GetFADCTrace();

    const fdet::Telescope& detTel = detEye.GetTelescope(mirrorId);
    const fdet::Camera& theCam = detTel.GetCamera();

    const double fadcBinSize = theCam.GetFADCBinSize();
    const unsigned int fadcBinLength = (unsigned int)(theCam.GetFADCTraceLength());
    const unsigned int fadcDynamicRange = (unsigned int)(theCam.GetADCDynamicRange());
    const int adcOverFlow = (1 << fadcDynamicRange) - 1;

    // copy the Raw trace into an utl::Trace

    utl::TraceI rawFADCTrace(fadcBinLength, fadcBinSize);

    const unsigned int startBin = fadcData->GetTraceStartBin();
    const unsigned int endBin = fadcData->GetTraceEndBin();

    // sanity check, MAM 09/23/04
    if (startBin > 200)
      continue;

    rawFADCTrace.SetStart(startBin);
    rawFADCTrace.SetStop(endBin);

    // loop on fadc bin
    bool isSaturated = false;
    bool firstSatBinSet = false;
    int firstSatBin = 0;
    int lastSatBin = 0;

    ostringstream info;
    info << "trace data: ";
    for (unsigned int pos = startBin; pos <= endBin; ++pos) {

      rawFADCTrace[pos] = int(FADCDataWordGetData(&fadcword[pos]));

      if (rawFADCTrace[pos] >= adcOverFlow) {
        isSaturated = true;
        if (!firstSatBinSet) {
          firstSatBinSet = true;
          firstSatBin = pos;
        }
        lastSatBin = pos;
      }
      if (fVerbosity) {
        info << rawFADCTrace[pos] << ", ";
      }
    }
    if (fVerbosity)
      INFO(info);

    // find baseline in first portion of trace
    double baseline = 0;
    try {
      baseline = TraceAlgorithm::Mean(rawFADCTrace, fFirstMeanSlot, fLastMeanSlot);
    } catch (utl::InvalidTraceBoundException& ex) {
      ERROR(ex.GetMessage());
      continue;
    }

    if (!tel.HasChannel(mirrorChannelId))
      tel.MakeChannel(mirrorChannelId);
    fevt::Channel& thisChannel = tel.GetChannel(mirrorChannelId);

    if (!thisChannel.HasRecData())
      thisChannel.MakeRecData();

    ChannelRecData& channelRecData = thisChannel.GetRecData();

    channelRecData.SetBaseLine(baseline);

    if (!channelRecData.HasFADCTrace(FdConstants::eTotal))
      channelRecData.MakeFADCTrace(rawFADCTrace, FdConstants::eTotal);
    else
      WARNING("channel recdata exists! Not overwriting....");

    if (isSaturated) {
      channelRecData.SetChannelSaturated();
      channelRecData.SetFirstSatBin(firstSatBin);
      channelRecData.SetLastSatBin(lastSatBin);
    }
  }
}


void
FdCalibrator::FillCalibratedPixels(evt::Event& event, TEyeEvent& eyeEvent)
{
  TEyePixelList* const eyePixelList = eyeEvent.GetPixelList();
  TEyeFADCData* const eyeFADCData = eyeEvent.GetFADCData();
  TEyePixelData* const eyePixelData = eyeEvent.GetPixelData();
  const unsigned int npixels = eyePixelList->GetNumPixels();

  const int eyeId = eyeEvent.GetEventHeader()->GetEyeNo();
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeId);

  FEvent& fdEvent = event.GetFEvent();
  fevt::Eye& eye = fdEvent.GetEye(eyeId);

  ostringstream info;
  info << "Eye " << eyeId << ": nPixels=" << npixels;
  int countSaturated = 0;
  int countErrors = 0;
  int countVirtuals = 0;
  int countNoCalib = 0;
  int countCalibrated = 0;

  for (unsigned int iPixel = 0; iPixel < npixels; ++iPixel) {

    FdUtil::Fd::FdPixelNumber pixelNumber = eyePixelList->GetPixel(iPixel);

    const unsigned int mirrorId = FdUtil::Fd::GetEyeMirrorNo(pixelNumber);
    const unsigned int mirrorChannelId = FdUtil::Fd::GetMirrorPixelNo(pixelNumber);

    if (!eye.HasTelescope(mirrorId))
      eye.MakeTelescope(mirrorId);
    fevt::Telescope& tel = eye.GetTelescope(mirrorId);

    const fdet::Telescope& detTel = detEye.GetTelescope(mirrorId);
    const fdet::Camera& detCamera = detTel.GetCamera();
    const fdet::Channel& thisChannel = detTel.GetChannel(mirrorChannelId);

    // virtual channels must be treated differently
    if (thisChannel.IsVirtual()) {
      ++countVirtuals;
      continue;
    }

    // copy fadc trace

    const TFADCData* const fadcData = eyeFADCData->GetFADCData(pixelNumber);

    // Protect yourself against the stupid non-existent traces...
    if (!fadcData) {
      ++countErrors;
      continue;
    }

    TFADCData::FADCDataWord fadcword = fadcData->GetFADCTrace();

    const double fadcBinSize = detCamera.GetFADCBinSize();
    const unsigned int fadcBinLength = (unsigned int)(detCamera.GetFADCTraceLength());
    const unsigned int fadcDynamicRange = (unsigned int)(detCamera.GetADCDynamicRange());
    const int adcOverFlow = (1 << fadcDynamicRange) - 1;

    // copy the Raw trace into an utl::Trace

    utl::TraceI rawFADCTrace(fadcBinLength, fadcBinSize);
    utl::TraceB rawFLTtrace(fadcBinLength, fadcBinSize);

    const unsigned int startBin = fadcData->GetTraceStartBin();
    const unsigned int endBin = fadcData->GetTraceEndBin();

    // sanity check, MAM 09/23/04
    if (startBin > 200)
      continue;

    rawFADCTrace.SetStart(startBin);
    rawFADCTrace.SetStop(endBin);
    rawFLTtrace.SetStart(startBin);
    rawFLTtrace.SetStop(endBin);

    bool hasSaturated = false;

    // loop on fadc bin
    for (unsigned int pos = startBin; pos <= endBin; ++pos) {

#if FDEVENTLIB_VERSION_CODE >= ModuleVersionCode(4, 0, 5)
      rawFLTtrace[pos] = bool(FADCDataWordGetPixelTrigger(&fadcword[pos]));
#else
      rawFLTtrace[pos] = bool(FADCDataWordPixelTrigger(&fadcword[pos]));
#endif
      rawFADCTrace[pos] = int(FADCDataWordGetData(&fadcword[pos]));

      if (rawFADCTrace[pos] >= adcOverFlow)
        hasSaturated = true;
    }

    // Correct for channel/pixel mismatch
    const fdet::Channel& detChannel = detTel.GetChannel(mirrorChannelId);
    unsigned int pixelId = detChannel.GetPixelId();
    if (mirrorChannelId != pixelId) {
      ostringstream info;
      info << "remapped pixel during calibration: channel=" << mirrorChannelId
           << " -> pixel=" << pixelId;
      INFO(info);
    }

    if (!tel.HasPixel(pixelId))
      tel.MakePixel(pixelId);
    fevt::Pixel& pixel = tel.GetPixel(pixelId);

    if (hasSaturated) {
      ++countSaturated;
      pixel.SetHighGainSaturation();
      if (fRecoverSaturatedTraces)
        RecoverSaturatedTrace(rawFADCTrace, pixel, tel);
    } else
      pixel.SetNoSaturation();

    if (!pixel.HasTriggerData())
      pixel.MakeTriggerData();

    PixelTriggerData& triggerdata = pixel.GetTriggerData();

    triggerdata.SetThreshold(fadcData->GetActualThreshold());
    triggerdata.SetPreviousThreshold(fadcData->GetPreviousThreshold());
    triggerdata.SetThresholdChanged();
    triggerdata.SetMean(fadcData->GetMean());
    triggerdata.SetVariance(fadcData->GetVariance());
    triggerdata.SetRate(fadcData->GetActualRate());

    triggerdata.MakeFLTTrace(rawFLTtrace);

    if (!ApplyCalibration(rawFADCTrace, pixel)) {
      pixel.SetStatus(ComponentSelector::eBadCalibration);
      ++countNoCalib;
      continue;
    }
    ++countCalibrated;

    // RU Wed May 18 10:08:27 CEST 2005
    // the FirstTriggeredTimeBin info was missing in event structure
    if (pixel.HasRecData()) {

      const TMirrorPixelData* const miPixelData =
        eyePixelData->GetPixelData(mirrorId);

      if (!miPixelData) {
        ERROR("No TMirrorPixelData available in TEyePixelData!");
        continue;
      }

      // note TMirrorPixelData::GetFirstTriggeredTimeBin(FdUtil::Fd::FdPixelNumber, unsigned int, unsigned int)
      // does currently (v4.1.2 not work correctly)
      const int triggerpos =
        miPixelData->GetFirstTriggeredTimeBin(FdUtil::Fd::GetColumnNo(FdUtil::Fd::GetPixelNumber(pixelNumber)),
                                              FdUtil::Fd::GetRowNo(FdUtil::Fd::GetPixelNumber(pixelNumber)));
      pixel.GetRecData().SetFirstTriggeredTimeBin(triggerpos);

    }
  }
  info << ", nCal=" << countCalibrated << ", nSat=" << countSaturated << ", nNoCal=" << countNoCalib
       << ", nErr=" << countErrors << ", nVirt=" << countVirtuals;
  INFO(info);
}


void
FdCalibrator::FillDAQInformation(evt::Event& event, TEyeEvent& eyeEvent)
{
  FEvent& fevent = event.GetFEvent();

  TEyeEventHeader* const eyeHeader = eyeEvent.GetEventHeader();
  const int eyeId = eyeHeader->GetEyeNo();

  if (!fevent.HasEye(eyeId, fevt::ComponentSelector::eInDAQ))
    fevent.MakeEye(eyeId, fevt::ComponentSelector::eInDAQ);
  fevt::Eye& eyeDAQ = fevent.GetEye(eyeId, fevt::ComponentSelector::eInDAQ);

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeDAQ);

  for (fdet::Eye::TelescopeIterator iTel = detEye.TelescopesBegin();
       iTel != detEye.TelescopesEnd(); ++iTel) {

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

    if (eyeHeader->IsMirrorPresent(telId) &&
        !eyeDAQ.HasTelescope(telId, fevt::ComponentSelector::eInDAQ))
      eyeDAQ.MakeTelescope(telId, fevt::ComponentSelector::eInDAQ);

  }
}


bool
FdCalibrator::ApplyTimeCorrections(evt::Event& event, TEyeEvent& eyeEvent)
{
  if (!event.HasFEvent())
    event.MakeFEvent();
  FEvent& fevent = event.GetFEvent();

  TEyeEventHeader* const eyeHeader =  eyeEvent.GetEventHeader();
  const int eyeId = eyeHeader->GetEyeNo();

  if (!fevent.HasEye(eyeId))
    fevent.MakeEye(eyeId);

  // correct trigger time if real data and not already done online
  if (eyeHeader->GetEventType() != TEyeEventHeader::kSimulated) {

#if FDEVENTLIB_VERSION_CODE >= ModuleVersionCode(2, 6, 0)
    TEyeEventHeader::ETimeCorrectionQuality corrQuality =
      eyeHeader->GetTimeCorrectionQuality();
    const bool badCorrection = (corrQuality < TEyeEventHeader::kFullQuality);
#else
    Float_t t10, offset;
    eyeHeader->GetTimeCorrectionParameters(&t10, &offset);
    const bool badCorrection = (t10 < 50 || t10 > 150);
#endif

    if (badCorrection) {

      EyeHeader& headerEye = fevent.GetEye(eyeId).GetHeader();
      headerEye.SetBadTimeCorrection(true);

      if (fCorrectTimes) {
        CorrectFDTime(eyeHeader);

        const unsigned int sec = eyeHeader->GetTimeStamp()->GetGPSSec();
        const unsigned int nsec = eyeHeader->GetTimeStamp()->GetNanoSec();

        const utl::TimeStamp tstamp(sec,nsec);
        Header& headerFEvent = fevent.GetHeader();

        headerEye.SetOfflineTimeCorrected(true);  // set a flag
        headerEye.SetTimeStamp(tstamp); // set the corrected time
        headerFEvent.SetTime(tstamp);
      }
    }
  }

  // check for time-zero events (bug in split-event merging...)
  FdRoot::TTimeStamp uninitializedTimeStamp(FdRoot::TTimeStamp::FromGPS(0, 0));
  if (eyeHeader->GetTimeStamp()->GetUnixSec() <= uninitializedTimeStamp.GetUnixSec()) {

    ostringstream info;
    info << "correcting GPS second=" <<  eyeHeader->GetTimeStamp()->GetSec()
         << " for eye ID="
         << eyeId;

    AdjustMirrorTimes(eyeHeader, fevent, fevent.GetEye(eyeId));

    info << " corrected GPS time: "
         <<  eyeHeader->GetTimeStamp()->GetSec() << '\n';

    // check if still bad ...
    if (eyeHeader->GetTimeStamp()->GetSec() <= uninitializedTimeStamp.GetUnixSec()) {
      info << "  ==> GPS correction failed - skipping eye!";
      WARNING(info);
      return false;
    }
  }

  // correction for the SD- FD time offset for this eye
  // Definition: SD_FD_Offset = SD_time - FD_time
  //           ===> SD_time = FD_time + SD_FD_Offset
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();

  const double SD_FD_Offset = detFD.GetEye(eyeId).GetSDTimeOffset().GetInterval();

  unsigned int sec = eyeHeader->GetTimeStamp()->GetGPSSec();
  const unsigned int nsec = eyeHeader->GetTimeStamp()->GetNanoSec() + int(SD_FD_Offset);

  // correction for wrong leap second in start of 2006
  // (note correction is from 4 Jan 2006 - 23 Jan 2006)
  if (sec >= 820374477 && sec < 822020847)
    ++sec;

  const utl::TimeStamp tstamp(sec, nsec);
  EyeHeader& headerEye = fevent.GetEye(eyeId).GetHeader();
  Header& headerFEvent = fevent.GetHeader();

  headerEye.SetTimeStamp(tstamp); // set the corrected time
  headerFEvent.SetTime(tstamp);

  return true;
}


bool
FdCalibrator::ApplyCalibration(const utl::TraceI& rawtrace, fevt::Pixel& pixel)
  const
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Pixel& detPixel = detFD.GetPixel(pixel);

  const int telId = detPixel.GetTelescopeId();
  const int eyeId = detPixel.GetEyeId();
  const int pixelId = detPixel.GetId();

  // Get the calibration constant for this channel
  const int channelId = detPixel.GetChannelId();
  const fdet::Pixel& channel =
    detFD.GetEye(eyeId).GetTelescope(telId).GetPixel(channelId);

  const TabulatedFunction* calibration = nullptr;
  try {
    calibration = &channel.GetEndToEndCalibrationConstant();
  } catch (utl::NonExistentComponentException& e)
  { }

  if (!calibration || !calibration->GetNPoints()) {
    ostringstream err;
    err << "Could not find calibration for eye " << eyeId
        << " tel " << telId << " pix " << pixelId;
    ERROR(err);
    return false;
  }

  // check calibration status
  if (channel.GetStatus() == fdet::Pixel::eBadCalibration) {
    ostringstream warn;
    warn << " bad calibration data for eye " << eyeId
         << " tel " << telId << " pix " << pixelId << " channel " << channelId
         << " status = " << detPixel.GetStatus();
    WARNING(warn);
    return false;
  } else if (channel.GetStatus() != fdet::Pixel::eGood) {
    ostringstream err;
    err << " unknown calibration status for eye " << eyeId
        << " tel " << telId << " pix " << pixelId << " channel " << channelId
        << " status = " << detPixel.GetStatus();
    ERROR(err);
    throw DataNotFoundInDBException(err.str());
  }

  const double bestCalibConst = channel.GetEndToEndCalibrationAtReferenceWavelength();

  // find baseline in first portion of trace
  double baseline = 0;
  try {
    baseline = TraceAlgorithm::Mean(rawtrace, fFirstMeanSlot, fLastMeanSlot);
  } catch (utl::InvalidTraceBoundException& ex) {
    ERROR(ex.GetMessage());
    return false;
  }

  // subtract baseline
  const TraceD baselineSubtractedTrace = (rawtrace - baseline);

  if (!pixel.HasRecData())
    pixel.MakeRecData();
  PixelRecData& recdata = pixel.GetRecData();

  // apply calibration constant
  if (recdata.HasPhotonTrace()) {
    WARNING("Pixel already calibrated ! Skipping!");
    return true;
  }

  recdata.MakePhotonTrace(baselineSubtractedTrace * bestCalibConst * fCalibrationCorrection[eyeId-1]);

  const TraceD& photontrace = recdata.GetPhotonTrace();

  // get rms
  const double rms =
    TraceAlgorithm::StandardDeviation(photontrace, fFirstMeanSlot, fLastMeanSlot);

  recdata.SetBaseline(0); // by definition
  recdata.SetRMS(rms);

  return true;
}


void
FdCalibrator::CorrectFDTime(TEyeEventHeader* const eyeHeader)
{
  const int eyeId = eyeHeader->GetEyeNo();
  det::Detector& theDet=det::Detector::GetInstance();
  const fdet::FDetector& detFD = theDet.GetFDetector();

  // code below is needed, if this is the first event to be processed (--> detector
  // not yet updated)
  const utl::TimeStamp detTime = theDet.GetTime();

  if (!detTime) {
    const unsigned int sec = eyeHeader->GetTimeStamp()->GetGPSSec();
    const unsigned int nsec = eyeHeader->GetTimeStamp()->GetNanoSec();
    const utl::TimeStamp tstamp(sec, nsec);
    theDet.Update(tstamp);
  }

  // get average correction value for this month
  const fdet::Eye& detEye = detFD.GetEye(eyeId);
  const utl::TabulatedFunction& timeCorrection =
    detEye.GetTimeCorrectionFactors();

  // there may be no values in FTelescopeList (warning printed in Eye.cc)
  if (timeCorrection.GetNPoints() < 1)
    return;

  int year, month, day;
  eyeHeader->GetTimeStamp()->GetDate(kTRUE, 0, &year, &month, &day);
  double iMonth = month + (year - 2004) * 12;

  // some out of bound checks ...
  if (iMonth < timeCorrection.GetX(0)) {
    INFO("Warning - request for time correction before 2004");
    iMonth = timeCorrection.GetX(0);
  } else if (iMonth > timeCorrection.GetX(timeCorrection.GetNPoints() - 1)) {
    INFO("Warning - request for time correction ouside range");
    iMonth = timeCorrection.GetX(timeCorrection.GetNPoints() - 1);
  }

  // recalculate mirror offsets and event timestamp
  const float t_10 = (1. - timeCorrection.Y(iMonth)) * 100.;
  const float offset = 4.*utl::ns;  // average

  cerr << " CorrectFDTime(): " << year << '/' << month << '/' << day << "\n"
          "              --> field: " << iMonth << " value: " << t_10 << "\n"
          "        before  : " << eyeHeader->GetTimeStamp()->GetGPSSec()
       << " s" << eyeHeader->GetTimeStamp()->GetNanoSec() << " ns\n";

  eyeHeader->SetTimeCorrectionParameters(t_10, offset);
  eyeHeader->AdjustMirrorTimes();

  cerr << "        after   : " << eyeHeader->GetTimeStamp()->GetGPSSec()
       << " s" << eyeHeader->GetTimeStamp()->GetNanoSec() << " ns" << endl;
}


// --> basically a copy of the TEyeEventHeader function
void
FdCalibrator::AdjustMirrorTimes(TEyeEventHeader* const eyeHeader,
                                fevt::FEvent& theFEvent,
                                fevt::Eye& theEye)
{
  FdRoot::TTimeStamp first(time_t(INT_MAX), 0);
  UInt_t first_mirror = 0;

  // --- find the earliest mirror

  for (UInt_t i = 0; i < FdUtil::Fd::kEYE_NMIRRORS; ++i) {
    const UInt_t j = i + FdUtil::Fd::kEYE_FIRST_MIRROR;
    if (eyeHeader->IsMirrorDataPresent(j)) {

      const FdRoot::TTimeStamp mirrorTimeStamp =
        *eyeHeader->GetRawMirrorTimeStamp(j);

      cerr << "\t FdCalibrator::AdjustMirrorTimes() -- mirror "
           << j
           << ", raw mirror time stamp "
           << mirrorTimeStamp << endl;

      if (mirrorTimeStamp < first) {
        first = mirrorTimeStamp;
        first_mirror = i;
      }
    }
  }

  // --- correct the time of the first mirror

  UInt_t sec_out, nsec_out;

  eyeHeader->DoTimeCorrection(first.GetSec(), first.GetNanoSec(), &sec_out, &nsec_out);

  // TTimeStamp has internally Unix time format

  FdRoot::TTimeStamp corrected(FdRoot::TTimeStamp::FromGPS(sec_out, nsec_out));
  eyeHeader->SetTimeStamp(corrected);

  const unsigned int sec = eyeHeader->GetTimeStamp()->GetGPSSec();
  const unsigned int nsec = eyeHeader->GetTimeStamp()->GetNanoSec();

  utl::TimeStamp tstamp(sec, nsec);
  theEye.GetHeader().SetTimeStamp(tstamp);
  theEye.GetHeader().SetOfflineTimeCorrected(true);
  fevt::Telescope& telescope = theEye.GetTelescope(first_mirror + FdUtil::Fd::kEYE_FIRST_MIRROR);
  telescope.SetTimeOffset(0);

  theFEvent.GetHeader().SetTime(tstamp);

  // --- adjust now the offsets for all mirrors which are present

  for (UInt_t i = 0; i < FdUtil::Fd::kEYE_NMIRRORS; ++i) {

    const UInt_t j = i + FdUtil::Fd::kEYE_FIRST_MIRROR;

    if (eyeHeader->IsMirrorDataPresent(j) && i != first_mirror) {

      const FdRoot::TTimeStamp mirrorTimeStamp =
        *eyeHeader->GetRawMirrorTimeStamp(j);
      eyeHeader->DoTimeCorrection(mirrorTimeStamp.GetSec(),
                                  mirrorTimeStamp.GetNanoSec(), &sec_out, &nsec_out);

      // now calculate offset, might be in next second
      UInt_t mirrorTimeOffset;

      if (sec_out > UInt_t(eyeHeader->GetTimeStamp()->GetSec()))
        mirrorTimeOffset = 1000000000 + nsec_out - eyeHeader->GetTimeStamp()->GetNanoSec();
      else
        mirrorTimeOffset = nsec_out - eyeHeader->GetTimeStamp()->GetNanoSec();
      fevt::Telescope& telescope = theEye.GetTelescope(j);
      telescope.SetTimeOffset(mirrorTimeOffset);

    }
  }
}


void
FdCalibrator::RecoverSaturatedTrace(utl::TraceI& rawFADCTrace,
                                    fevt::Pixel& pixel,
                                    fevt::Telescope& telEvent)
{
  const unsigned int eyeId = telEvent.GetEyeId();
  const unsigned int telId = telEvent.GetId();

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& detEye = detFD.GetEye(eyeId);
  const fdet::Telescope& detTel = detEye.GetTelescope(telId);

  const fdet::Camera& detCamera = detTel.GetCamera();
  const unsigned int fadcDynamicRange = (unsigned int) detCamera.GetADCDynamicRange();
  const int adcOverFlow = (1 << fadcDynamicRange) - 1;

  const fdet::Channel& detChannel = detTel.GetChannel(pixel);
  const unsigned int channelId = detChannel.GetId();
  const double channelGain = detChannel.GetElectronicsGain();
  const unsigned int virtualChannelId = detChannel.GetVirtualChannelId();
  const double virtualChannelGain = detTel.GetChannel(virtualChannelId).GetElectronicsGain();

  ChannelRecData& channelRecData = telEvent.GetChannel(channelId).GetRecData();

  const unsigned int startBin = channelRecData.GetFirstSatBin();
  const unsigned int endBin = channelRecData.GetLastSatBin();

  // prepare list of channel data sharing the same virtual channel
  const ChannelRecData* vChannelRecData = nullptr;
  vector<const ChannelRecData*> otherChannels;
  vector<double> electronicGains;

  for (fevt::Telescope::ConstChannelIterator chanIter = telEvent.ChannelsBegin();
       chanIter != telEvent.ChannelsEnd(); ++chanIter) {

    const unsigned int thisChannelId = chanIter->GetId();
    if (!chanIter->HasRecData())
      continue;

    const ChannelRecData& recData = chanIter->GetRecData();
    const fdet::Channel& thisDetChannel = detTel.GetChannel(thisChannelId);

    if (thisDetChannel.IsVirtual()) {
      if (thisChannelId == virtualChannelId) { // virtual channel for de-saturation
        if (recData.ChannelIsSaturated()) {
          const unsigned int t1 = recData.GetFirstSatBin();
          const unsigned int t2 = recData.GetLastSatBin();
          if ((startBin <= t1 && t1 <= endBin) ||
              (startBin <= t2 && t2 <= endBin)) {
            WARNING("trace recovery failed because virtual channel"
                    " is saturated!");
            pixel.SetLowGainSaturation();
            return;
          }
        }
        vChannelRecData = &recData;
      }

    } else if (thisChannelId != channelId) {

      const unsigned int thisVirtualChannelId =
        detTel.GetChannel(thisChannelId).GetVirtualChannelId();

      if (thisVirtualChannelId == virtualChannelId) { // other pixels in same virtual channel
        if (recData.ChannelIsSaturated()) {
          const unsigned int t1 = recData.GetFirstSatBin();
          const unsigned int t2 = recData.GetLastSatBin();
          if ((startBin <= t1 && t1 <= endBin) ||
              (startBin <= t2 && t2 <= endBin)) {
            WARNING("more than one saturated channels in same time range!!");
            return;
          }
        }
        otherChannels.push_back(&recData);
        electronicGains.push_back(thisDetChannel.GetElectronicsGain());
      }
    }
  }

  if (!vChannelRecData) {
    ERROR("Missing virtual channel data!!");
    pixel.SetLowGainSaturation();
    return;
  }

  const TraceI& virtualTrace = vChannelRecData->GetFADCTrace(FdConstants::eTotal);

  // correct trace according to Eq. (1) in GAP-2006-104

  for (unsigned int pos = startBin; pos <= endBin; ++pos) {
    const int adcValue = rawFADCTrace[pos];
    if (adcValue >= adcOverFlow) {
      double otherChannelADCSum = 0;
      for (unsigned int i = 0; i < otherChannels.size(); ++i) {
        const TraceI& thisTrace = otherChannels[i]->GetFADCTrace(FdConstants::eTotal);
        const double adcValue = thisTrace[pos] - otherChannels[i]->GetBaseLine();
        otherChannelADCSum += adcValue/electronicGains[i];
      }
      const double vADCBaseLineSub = virtualTrace[pos] - vChannelRecData->GetBaseLine();
      const double recoveredSignal =
        (vADCBaseLineSub/virtualChannelGain - otherChannelADCSum) * channelGain;
      rawFADCTrace[pos] = int(recoveredSignal + 0.5);
      channelRecData.AddDeSaturatedBin(pos);
    }
  }
  pixel.SetSaturationRecovered();
}


double
FdCalibrator::GetHeatCalibrationCorrection(const unsigned int gpsSecond)
  const
{
  int thisBin = -1;
  for (unsigned int i = 0; i < fHeatCaliEpochs.size(); ++i) {
    if (fHeatCaliEpochs[i] > gpsSecond) {
      thisBin = i - 1;
      break;
    }
  }

  int lastBin = fHeatCaliEpochs.size() - 2;

  std::vector<double> p;
  if (thisBin == -1) {
    for (unsigned int i = 0; i < 3; ++i)
      p.push_back(fHeatCalibration[3*lastBin + i]);
  } else {
    for (unsigned int i = 0; i < 3; ++i)
      p.push_back(fHeatCalibration[3*thisBin + i]);
  }

  /*
  cout <<"thisBin: " << thisBin <<"\n";
  cout << "gpsSecond " << gpsSecond << "\n";
  cout << p[0] << " " << p[1] << " " << p[2] <<"\n";
  */
  unsigned int lastEpoch = fHeatCaliEpochs[lastBin+1];
  int x;

  if (thisBin > -1)    // means that  gpsSecond is within Epoch range
    x = gpsSecond - 1016048411;
  else                 // means that  gpsSecond is greater than Epoch range
    x = lastEpoch - 1016048411;

  const double HEATcali =
    p[0] + p[1] * x + p[2] * pow(x, 2);

  return HEATcali;
}