FdLightCollectionEfficiency.cc

Go to the documentation of this file.

#include "FdLightCollectionEfficiency.h"
#include "SimMockEvent.h"
#include "SubModule.h"
#include "PixelToApertureLightConverter.h"
#include "Bootstrapper.h"
#include "TimeRangeCalculator.h"
#include "CleanupSentry.h"
#include "CorrectedApLightCalculator.h"

//#define NDEBUG
#undef NDEBUG
#include <assert.h>

// switch for writing debugging graphs to a ROOT file
//#define DEBUG
#undef DEBUG

#include <vector>
#include <map>
#include <string>

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

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

#include <utl/TabulatedFunction.h>

#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/TimeStamp.h>

#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/AugerCoordinateSystem.h>

#include <utl/Trace.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/VGaisserHillasParameter.h>
#include <evt/GaisserHillas2Parameter.h>
#include <evt/GaisserHillas4Parameter.h>

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

#include <det/Detector.h>
#include <atm/Atmosphere.h>

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

#include <utl/AugerException.h>
#include <utl/Vector.h>

#include <sstream>


// modules to be remote controlled
#include <FdReconstruction/FdProfileReconstructorKG/FdProfileReconstructor.h>
#include <FdSimulation/TelescopeSimulatorKG/TelescopeSimulator.h>
#include <FdSimulation/ShowerPhotonGeneratorOG/ShowerPhotonGenerator.h>
#include <FdSimulation/LightAtDiaphragmSimulatorKG/LightAtDiaphragmSimulator.h>
// just for the config MD5 hash generation for the calib. constants
#include <FdSimulation/FdElectronicsSimulatorOG/FdElectronicsSimulator.h>


// Still only for debugging, but then I'll get unused parameter errors. gcc--
//#ifdef DEBUG
#include <TFile.h>
#include <TH1D.h>
#include <TGraph.h>
#include <TGraphErrors.h>
//#endif

using namespace FdLightCollectionEfficiencyKG;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace det;
using namespace fdet;
using namespace std;


/*************************************************************************/
fwk::VModule::ResultFlag
FdLightCollectionEfficiency::Init()
{
  // Note: eCherDirect and eCherMieScattered will be automatically
  //       copied from eFluorDirect and eCherRayleighScattered
  //       efficiencies if missing!
  fLightComponents.push_back(FdConstants::eFluorDirect);
//  fLightComponents.push_back(FdConstants::eCherDirect);
  fLightComponents.push_back(FdConstants::eCherRayleighScattered);
//  fLightComponents.push_back(FdConstants::eCherMieScattered);

  // Read configuration
  CentralConfig* cc = CentralConfig::GetInstance();

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

  Branch profCalcB = topB.GetChild("profileCalculation");
  // steps in depth for calculation of the dEdX and Ne profiles
  profCalcB.GetChild("depthBinning").GetData(fProfileXBinning);

  Branch lightAtDiaB = topB.GetChild("lightAtDiaphragm");
  // factor the the number of bins in light@dia simulation
  lightAtDiaB.GetChild("binning").GetData(fLightAtDiaBinning);

  Branch lightCollB = topB.GetChild("lightCollection");
  lightCollB.GetChild("useZeta").GetData(fUseZeta);

  topB.GetChild("verbosity").GetData(fVerbosity);

  // Get the ShowerPhotonGeneratorOG overrides (ray tracing no. photons)
  Branch rayTracingB = topB.GetChild("rayTracing");
  rayTracingB.GetChild("minNRayTracePerBin").GetData(fMinNRayTracePerBin);
  rayTracingB.GetChild("maxNRayTracePerBin").GetData(fMaxNRayTracePerBin);
  rayTracingB.GetChild("extraRayTraceFactor").GetData(fExtraRayTraceFactor);
  rayTracingB.GetChild("nIterations").GetData(fNRayTracingIterations);
  rayTracingB.GetChild("targetUncertainty").GetData(fTargetUncertainty);
  rayTracingB.GetChild("minRelevantEfficiency").GetData(fMinRelevantEfficiency);
  rayTracingB.GetChild("maxIterations").GetData(fMaxIterations);
  rayTracingB.GetChild("targetApertureLightUncertaintyChange").GetData(fTargetWorstCaseLAtApUncertaintyChange);
  rayTracingB.GetChild("lowerLimitOnUncertainty").GetData(fLowerUncertaintyLimit);

  std::string stopCondition;
  rayTracingB.GetChild("stopCondition").GetData(stopCondition);
  if (stopCondition == "eFixed")
    fStopCondition = eFixed;
  else if (stopCondition == "eTargetUncertainty")
    fStopCondition = eTargetUncertainty;
  else if (stopCondition == "eSmallEnough")
    fStopCondition = eSmallEnough;
  else if (stopCondition == "eBootstrap")
    fStopCondition = eBootstrap;
  else
    throw utl::XMLParseException("Invalid stop conditon. No such enumeration value.");

  // TODO: flexible sub-module configuration and selection from XML

  // This is the list of all modules we "control"
  if (fSubModulesList.size() == 0) {
    fSubModulesList["FdProfileReconstructorKG"]    = SubModule("FdProfileReconstructorKG", fVerbosity);
    fSubModulesList["ShowerLightSimulatorKG"]      = SubModule("ShowerLightSimulatorKG", fVerbosity);
    fSubModulesList["LightAtDiaphragmSimulatorKG"] = SubModule("LightAtDiaphragmSimulatorKG", fVerbosity);
    fSubModulesList["ShowerPhotonGeneratorOG"]     = SubModule("ShowerPhotonGeneratorOG", fVerbosity);
    fSubModulesList["TelescopeSimulatorKG"]        = SubModule("TelescopeSimulatorKG", fVerbosity);
    fSubModulesList["FdElectronicsSimulatorOG"]    = SubModule("FdElectronicsSimulatorOG", fVerbosity);
  }

  // initialize the sub-modules
  for (map<string, SubModule>::iterator iModule = fSubModulesList.begin();
       iModule != fSubModulesList.end(); ++iModule) {
    SubModule& subModule = iModule->second;
    fwk::VModule::ResultFlag status = subModule.Init();
    if (status != eSuccess) {
      std::ostringstream msg;
      msg << "SubModule '" << subModule.GetName() << "' failed Init()" << endl;
      WARNING(msg);
      return status;
    }
  }

  if (fVerbosity > 0) {
    cout << "--[FdLightCollectionEfficiencyKG]--> Configuration:\n"
         << "  verbosity:                          " << fVerbosity << "\n"
         << "  simulated profile depth binning:    " << fProfileXBinning/(g/cm/cm) << "\n"
         << "  use only pixels within zeta:        " << (fUseZeta ? "true" : "false") << "\n"
         << "  number of ray-tracing iterations:   " << fNRayTracingIterations << "\n"
         << "  min no. ray traced photons per bin: " << fMinNRayTracePerBin << "\n"
         << "  max no. ray traced photons per bin: " << fMaxNRayTracePerBin << "\n"
         << "  extra-ray-tracing-factor:           " << fMaxNRayTracePerBin << "\n"
         << endl;
  }

  if (fVerbosity > -1) {
    if (!fUseZeta) {
      WARNING("Not using the light collection angle Zeta for calculating the light collection"
              " efficiency. If you are not debugging the FdLightCollectionEfficiencyKG module,"
              " then this is very likely an error in your configuration");
    }
  }

  return eSuccess;
}


/*************************************************************************/
fwk::VModule::ResultFlag
FdLightCollectionEfficiency::Run(evt::Event& event)
{
  if (!event.HasFEvent())
    return eSuccess;

  return ProcessEvent(event);
}

/*************************************************************************/
fwk::VModule::ResultFlag
FdLightCollectionEfficiency::ProcessEvent(evt::Event& realEvent)
{
  //const RunController& runCtrl = RunController::GetInstance();

  typedef map<int, ShowerFRecData*> FRecMap;
  FRecMap dataSets;
  const fwk::VModule::ResultFlag status = DoInitialReconstruction(realEvent, dataSets);
  if (status != eSuccess)
    return status;

  // Do simulation for each data set (each eye)
  //  FRecMap::iterator dataSetsIter = dataSets.begin(); // redundant definition
  for (FRecMap::iterator dataSetsIter = dataSets.begin(), end = dataSets.end();
       dataSetsIter != end; ++dataSetsIter)
  {
    const int eyeId = dataSetsIter->first;
    ShowerFRecData& frecData = *dataSetsIter->second;

    // Don't clone. Some modules may be irritated if the
    // event contains what they presume to be their own result before
    // they even start
    //evt::Event myEvent = realEvent;
    SimMockEvent myEvent(fProfileXBinning);
    myEvent.SetVerbosity(fVerbosity);

    bool atLeastOneTel = myEvent.PrepareEvent(eyeId, frecData, realEvent,
                                              realEvent.GetFEvent().GetEye(eyeId));
    if (!atLeastOneTel) {
      ostringstream msg;
      msg << "Skipping eye " << eyeId << ": Either not in DAQ or no good telescope.";
      INFO(msg);
      continue;
    }

    fwk::VModule::ResultFlag status;
    status = fSubModulesList.find("ShowerLightSimulatorKG")->second.Run(myEvent.GetEvent());
    if (status != eSuccess)
      return status;

    status = fSubModulesList.find("LightAtDiaphragmSimulatorKG")->second.Run(myEvent.GetEvent());
    if (status != eSuccess)
      return status;

    // Run LightAtDiaphragmSimulatorKG
    SubModule& lightAtDiaMod = fSubModulesList.find("LightAtDiaphragmSimulatorKG")->second;
    /*LightAtDiaphragmSimulatorKG::LightAtDiaphragmSimulator& lightAtDia =
      dynamic_cast<LightAtDiaphragmSimulatorKG::LightAtDiaphragmSimulator&>(lightAtDiaMod.GetModule());*/

#warning THIS DOES NOT WORK LIKE THIS ANY MORE!!!!!!!!!!!!
    //    const double lightAtDiaBinning = lightAtDia.GetBinning(); // remember
    //    lightAtDia.SetBinning(fLightAtDiaBinning); // overwrite
    status = lightAtDiaMod.Run(myEvent.GetEvent()); // RUN
    //    lightAtDia.SetBinning(lightAtDiaBinning); // reset

    if (status != eSuccess)
      return status;

    // Initialize cached converter for light-at-pixel to light-at-aperture conversion
    fLightConverter = new PixelToApertureLightConverter(det::Detector::GetInstance().GetFDetector(),
                                                        myEvent.GetEvent().GetFEvent().GetEye(eyeId, ComponentSelector::eExists));
    CleanupSentry<PixelToApertureLightConverter> converterGuard(fLightConverter); // this'll clean up before scope exit

    for (unsigned int iLightSource = 0; iLightSource < fLightComponents.size(); ++iLightSource)
    {
      const fevt::FdConstants::LightSource lightSource = fLightComponents[iLightSource];

      fTimeCorrectedApertureTraces = new std::map<unsigned int, utl::TraceD>();
      // this'll clean up before scope exit:
      CleanupSentry<std::map<unsigned int, utl::TraceD> > apTracesGuard(fTimeCorrectedApertureTraces);

      // Run photon-generation and ray tracing iteratively
      std::pair<VModule::ResultFlag, unsigned int> tmp
        = RunPhotonGenerationRayTracingLoop(myEvent.GetEvent(), realEvent,
                                            ( (fStopCondition==eTargetUncertainty
                                               ||fStopCondition==eSmallEnough)
                                              ? 1
                                              : (unsigned int)fNRayTracingIterations ),
                                            lightSource);
      status = tmp.first;
      const unsigned int nRayTraced = tmp.second;
      if (status != eSuccess)
        return status; // reporting done in RunPhotonGenerationRayTracingLoop

      // Note: We do not run FdElectronicsSimulatorOG at all,
      //       because it is only needed for the config hash calculation

      // The actual efficiency calculation...
      if (fStopCondition != eBootstrap)
        CalculateEfficiency(myEvent.GetEvent(), realEvent, eyeId,
                            nRayTraced, lightSource);
    } // end for each light source

#ifdef DEBUG
    TFile* theFile = NULL;
    if (fStopCondition == eBootstrap)
      theFile = new TFile("lightCollectionEfficiency.root", "UPDATE");
#endif /* DEBUG */

    // Dual purpose loop:
    // - Calculate lc-eff start/stop times
    // - Copy Fluorescence LCEff to direct Cherenkov LCEff
    //   and Rayleigh sc. Ch. light to Mie sc. Ch. light
    for (fevt::FEvent::EyeIterator iEye = realEvent.GetFEvent().EyesBegin(ComponentSelector::eHasData);
         iEye != realEvent.GetFEvent().EyesEnd(ComponentSelector::eHasData); ++iEye)
    {
      fevt::Eye& realEye = *iEye;
      for (fevt::Eye::TelescopeIterator iTel = realEye.TelescopesBegin(ComponentSelector::eHasData);
           iTel != realEye.TelescopesEnd(ComponentSelector::eHasData); ++iTel)
      {
        fevt::Telescope& realTel = *iTel;
        if (!realTel.HasRecData())
          continue;
        fevt::TelescopeRecData& recData = realTel.GetRecData();
        if (!recData.HasLightCollectionEfficiency())
          continue;

        utl::MultiTabulatedFunctionErrors& lcEff = recData.GetLightCollectionEfficiency();

    // lceff/aperture light time debugging couts...
    /*    const utl::TabulatedFunctionErrors& l = lcEff.GetTabulatedFunctionErrors(FdConstants::eFluorDirect);
        const utl::TabulatedFunctionErrors& ap = recData.GetLightFlux(FdConstants::eTotal);
        const list<pair<double, double> > ss = recData.GetSpotFarFromBorderTimeRanges();
        cout << "eye " << realEye.GetId() << " tel " << realTel.GetId() << endl;
        cout << "  AP start time:    " << ap.GetX(0) << endl;
        cout << "  lceff start time: " << l.GetX(0) << endl;
        cout << "  time ranges start times: ";
        for (list<pair<double, double> >::const_iterator it = ss.begin(), end=ss.end();
             it != end; ++it)
          cout << it->first << " ";
        cout << endl;
        cout << "  AP stop time:     " << ap.GetX(ap.GetNPoints()-1) << endl;
        cout << "  lceff stop time:  " << l.GetX(l.GetNPoints()-1) << endl;
        cout << "  time ranges stop times:  ";
        for (list<pair<double, double> >::const_iterator it = ss.begin(), end=ss.end();
             it != end; ++it)
          cout << it->second << " ";
        cout << endl;
    */


        // Update l@aperture validity start/stop times
        TimeRangeCalculator calc(recData, fVerbosity);
        calc.UpdateSpotFarFromBorderTimes(fMinRelevantEfficiency);

        // eFluorDirect => eCherDirect if necessary
        if (lcEff.HasLabel(FdConstants::eFluorDirect)
            && !lcEff.HasLabel(FdConstants::eCherDirect))
          lcEff.AddTabulatedFunctionErrors(lcEff.GetTabulatedFunctionErrors(FdConstants::eFluorDirect),
                                           FdConstants::eCherDirect);

        // eCherRayleighScattered => eCherMieScattered if necessary
        if (lcEff.HasLabel(FdConstants::eCherRayleighScattered)
            && !lcEff.HasLabel(FdConstants::eCherMieScattered))
          lcEff.AddTabulatedFunctionErrors(lcEff.GetTabulatedFunctionErrors(FdConstants::eCherRayleighScattered),
                                           FdConstants::eCherMieScattered);
#ifdef DEBUG
        if (fStopCondition == eBootstrap) {
          for (utl::MultiTabulatedFunctionErrors::Iterator componentIt = lcEff.Begin(), end = lcEff.End();
               componentIt != end; ++componentIt)
          {
            int component = componentIt->GetLabel();
            const utl::TabulatedFunctionErrors& thisLCEff = lcEff.GetTabulatedFunctionErrors(component);
            TGraphErrors* relHist = new TGraphErrors(thisLCEff.GetNPoints());
            for (unsigned int i = 0; i < thisLCEff.GetNPoints(); ++i) {
              relHist->SetPoint(i+1, thisLCEff.GetX(i), thisLCEff.GetY(i));
              relHist->SetPointError(i+1, thisLCEff.GetXErr(i), thisLCEff.GetYErr(i));
            }
            ostringstream relname;
            relname << "relHistFinal_e" << iEye->GetId() << "_t" << iTel->GetId() << "_c" << component;
            relHist->SetNameTitle(relname.str().c_str(), relname.str().c_str());
            relHist->Write();
          } // end foreach light component
        } // end if do bootstrap
#endif /* DEBUG */
      } // end foreach telescope
    } // end for eye

#ifdef DEBUG
    if (fStopCondition == eBootstrap)
      theFile->Close();
#endif /* DEBUG */

  } // end foreach data set

  return eSuccess;
}


/*************************************************************************/
fwk::VModule::ResultFlag
FdLightCollectionEfficiency::Finish()
{
  utl::TabularStream timeTab = SubModule::TimingHeader();

  // finish the sub-modules
  for (map<string, SubModule>::iterator iModule = fSubModulesList.begin(), end = fSubModulesList.end();
       iModule != end; ++iModule) {
    SubModule& subModule = iModule->second;
    fwk::VModule::ResultFlag status = subModule.Finish();
    if (status != eSuccess) {
      cerr << "SubModule '" << subModule.GetName() << "' failed Finish()" << endl;
      return status;
    }
    if (fVerbosity >= 0)
      subModule.PrintTiming(timeTab);
  }

  // Print timing information
  if (fVerbosity >= 0)
    SubModule::FinishTimingPrintout(timeTab);

  return eSuccess;
}


/*************************************************************************/
bool
FdLightCollectionEfficiency::CalculateEfficiency(const evt::Event& simEvent, evt::Event& realEvent,
                                                 const int eyeId, const unsigned int nRayTracingIterations,
                                                 const fevt::FdConstants::LightSource lightSource)
{
  if (fVerbosity > 0)
    cout << "--[FdLightCollectionEfficiency]--> Calculating efficiency..." << endl;

  if (!simEvent.GetFEvent().HasEye(eyeId, ComponentSelector::eInDAQ)) {
    ostringstream msg;
    msg << "Efficiency calculation for eye " << eyeId
        << ": FEvent doesn't have that eye after simulation";
    ERROR(msg);
    return false;
  }

  const fevt::Eye& simEye = simEvent.GetFEvent().GetEye(eyeId, ComponentSelector::eInDAQ);
  fevt::Eye& realEye = realEvent.GetFEvent().GetEye(eyeId, ComponentSelector::eHasData);

  for (fevt::Eye::TelescopeIterator iTel = realEye.TelescopesBegin(ComponentSelector::eHasData);
       iTel != realEye.TelescopesEnd(ComponentSelector::eHasData); ++iTel) {
    const unsigned int telId = iTel->GetId();
    if (!realEye.HasTelescope(telId))
      continue;
    if (!realEye.GetTelescope(telId).HasRecData())
      continue;

    try {
      const fevt::Telescope& simTel = simEye.GetTelescope(iTel->GetId(), ComponentSelector::eInDAQ);
      if (!CalculateTelescopeEfficiency(simTel, *iTel, realEye, nRayTracingIterations, lightSource))
        continue;
    } catch (utl::NonExistentComponentException& e) {
      cout << "Caught exception: " << e.GetMessage() << " => Skipping." << endl;
      continue;
    }
  }

  return true;
}


/*************************************************************************/
bool
FdLightCollectionEfficiency::CalculateTelescopeEfficiency(const fevt::Telescope& simTel,
                                                          fevt::Telescope& realTel,
                                                          const fevt::Eye& realEye,
                                                          const unsigned int /*nRayTracingIterations*/,
                                                          const fevt::FdConstants::LightSource lightSource)
{
  if (!simTel.HasSimData()) {
    if (fVerbosity > 1)
      cout << "--[FdLightCollectionEfficiency]--> Skipping efficiency calculation for telescope "
           << simTel.GetId() << " of eye " << simTel.GetEyeId() << " (lack of sim data)\n";
    return false;
  }

  if (!realTel.HasRecData()
      || !realTel.GetRecData().HasLightFlux(fevt::FdConstants::eTotal))
  {
    if (fVerbosity > 1)
      cout << "--[FdLightCollectionEfficiency]--> Skipping efficiency calculation for telescope "
           << simTel.GetId() << " of eye " << simTel.GetEyeId() << " (lack of telescope total l@ap)\n";
    return false;
  }


  if (fVerbosity > 0) {
    cout << "--[FdLightCollectionEfficiency]--> Calculating efficiency for telescope "
         << simTel.GetId() << " of eye " << simTel.GetEyeId()
         << "...\n";
  }

  const utl::TimeStamp& tracesStartTime      = simTel.GetTracesStartTime();
  const fevt::TelescopeSimData& telSimData   = simTel.GetSimData();
  const utl::TimeStamp& simApertureStartTime = telSimData.GetPhotonsStartTime();

  if (!realTel.HasRecData()) {
    ostringstream msg;
    msg << "Telescope " << realTel.GetId() << " of eye " << realTel.GetEyeId()
        << " has no TelescopeRecData => Skipping.";
    WARNING(msg);
    return false;
  }

  fevt::TelescopeRecData& telRecData = realTel.GetRecData();

  // Get this telescope's time-corrected aperture light trace
  if (fTimeCorrectedApertureTraces->find(simTel.GetId())
      == fTimeCorrectedApertureTraces->end()) {
    ostringstream msg;
    msg << "Telescope " << realTel.GetId() << " of eye " << realTel.GetEyeId()
        << " has no corrected light at aperture trace!";
    ERROR(msg);
    return false;
  }
  const utl::TraceD& correctedApTrace = (*fTimeCorrectedApertureTraces)[simTel.GetId()];

  if (!telRecData.HasLightCollectionEfficiency())
    telRecData.MakeLightCollectionEfficiency();

#ifdef DEBUG
  TFile* const theFile = new TFile("lightCollectionEfficiency.root", "UPDATE");
  WARNING("ENABLED DEBUG OUTPUT TO ROOT FILE 'lightCollectionEfficiency.root'");
#endif

  const vector<vector<unsigned int> >& pixelsInZeta = telRecData.GetPixelsInZetaOverTime();
  const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(realEye.GetId());
  const fdet::Channel& detChannel = detEye.GetTelescope(realTel.GetId()).GetChannel(1);

  const utl::TimeStamp eyeTriggerTime =
    realEye.GetHeader().GetTimeStamp() -
    TimeInterval(detChannel.GetFADCTraceLength() * detChannel.GetFADCBinSize() * nanosecond);

  // debugging output only
  if (fVerbosity > 2) {
    const utl::TimeStamp simApLightStartTime =
      eyeTriggerTime + TimeInterval(simTel.GetTimeOffset()*nanosecond);
    const utl::TimeStamp realApLightStartTime =
      eyeTriggerTime + TimeInterval(realTel.GetTimeOffset()*nanosecond);

    cout << "================================================\n"
            "= FdLightCollectionEfficiency debugging output =\n"
            "================================================\n"
            "Event time stamp:                      " << realEye.GetHeader().GetTimeStamp() << "\n"
            "Eye Trigger time:                      " << eyeTriggerTime
         << "(GPS nano: " << eyeTriggerTime.GetGPSNanoSecond() << ")\n"
            "Real Aperture start time:              " << realApLightStartTime
         << "(GPS nano: " << realApLightStartTime.GetGPSNanoSecond() << ")\n"
            "Calc. sim. start time:                 " << simApLightStartTime
         << "(GPS nano: " << simApLightStartTime.GetGPSNanoSecond() << ")\n"
            "Sim. Aperture start time (from trace): " << simApertureStartTime
         << "(GPS nano: " << simApertureStartTime.GetGPSNanoSecond() << ")\n"
            "Sim. Traces   start time:              " << tracesStartTime
         << "(GPS nano: " << tracesStartTime.GetGPSNanoSecond() << ")\n"
            "Diff (sim-real)/ns: " << (simApertureStartTime-realApLightStartTime).GetNanoSecond() << "\n"
            "real TelTimeOff: " << realTel.GetTimeOffset() << "\n"
            "trace length: " << TimeInterval(detChannel.GetFADCTraceLength()
                                * detChannel.GetFADCBinSize() * nanosecond).GetNanoSecond() << endl;
  } // end debugging output only

  const double simTimeOffset = (simApertureStartTime - eyeTriggerTime).GetInterval() / ns;
  if (fVerbosity > 0)
    cout << "  Simulation time offset: " << simTimeOffset << "ns" << endl;

  const utl::TabulatedFunctionErrors& realApLightFlux = telRecData.GetLightFlux();

  atm::Atmosphere::EmissionMode emode;
  switch (lightSource) {
  case FdConstants::eFluorDirect:
    emode = atm::Atmosphere::eFluorescence;
    break;
  case FdConstants::eCherDirect:
  case FdConstants::eCherRayleighScattered:
  case FdConstants::eCherMieScattered:
    emode = atm::Atmosphere::eCerenkov;
    break;
  default:
    throw utl::NonExistentComponentException("Unsupported LightSource!");
  }

  const utl::TraceD totalWlTrace2 = CalculateTelescopeTraceSum(simTel, lightSource, emode);
  const utl::TraceD& totalWlTrace = correctedApTrace; // FIXME: Temporary override...
  cout << "Binning: " << totalWlTrace2.GetBinning() << " " << totalWlTrace.GetBinning() << endl;
  cout << "Size: "    << totalWlTrace2.GetSize() << " " << totalWlTrace.GetSize() << endl;
  cout << "Start: "    << totalWlTrace2.GetStart() << " " << totalWlTrace.GetStart() << endl;
  cout << "Stop: "    << totalWlTrace2.GetStop() << " " << totalWlTrace.GetStop() << endl;
  double sum1 = 0., sum2 = 0.;
  for (unsigned int i = 0; i< totalWlTrace.GetSize();++i) {
    sum1 +=totalWlTrace[i];
    sum2 += totalWlTrace2[i];
  }
  cout << "Sum:" << sum1 << " " << sum2 << endl;

  const double wlTraceBinSize = totalWlTrace.GetBinning();
  const double wlTraceBins = totalWlTrace.GetSize();
  if (wlTraceBins == 0) {
    if (fVerbosity > 0) {
      cout << "--[FdLightCollectionEfficiency]--> Skipping light component "
           << lightSource << " for this telescope. No light-at-aperture trace.\n";
    }
    return false;
  }

  const utl::TabulatedFunctionErrors totalPixelTrace =
    CalculatePixelTraceSum(simTel, lightSource, simTimeOffset, pixelsInZeta, realApLightFlux);

  const double pixTraceBins = totalPixelTrace.GetNPoints();
  if (pixTraceBins == 0) {
    if (fVerbosity > 0) {
      cout << "--[FdLightCollectionEfficiency]--> Skipping light component "
           << lightSource << " for this telescope. No pixel traces.\n";
    }
    return false;
  }

  if (fVerbosity > 1)  {
    cout << "Component " << lightSource << " wlTrace:    nBins="
         << wlTraceBins << ", binSize=" << wlTraceBinSize << "\n"
            "Component " << lightSource << " pixTrace:   nBins="
         << pixTraceBins << "\n"
            "Component " << lightSource << " nRayTraced: nBins="
         << simTel.GetSimData().GetRayTracedPhotonTrace().GetSize()
         << ", binSize=" << simTel.GetSimData().GetRayTracedPhotonTrace().GetBinning() << endl;
  }

  utl::TabulatedFunctionErrors telEff =
    CalcTraceDivision(totalPixelTrace, totalWlTrace,
                      simTel.GetSimData().GetRayTracedPhotonTrace(),
                      1);
                      //nRayTracingIterations);

#ifdef DEBUG
  theFile->cd();
  WriteDebugInfo(realEye.GetId(), simTel.GetId(), simTimeOffset, lightSource, totalPixelTrace, totalWlTrace, realApLightFlux, telEff);
#endif

  // Store the light collection efficiency for this light source
  // in the TelescopeSimData

  utl::MultiTabulatedFunctionErrors& lightCollEff = telRecData.GetLightCollectionEfficiency();
  if (lightCollEff.HasLabel(lightSource))
    lightCollEff.GetTabulatedFunctionErrors(lightSource) = telEff;
  else
    lightCollEff.AddTabulatedFunctionErrors(telEff, lightSource);

#ifdef DEBUG
  theFile->Close();
#endif

  return true;
}


/*************************************************************************/
void
FdLightCollectionEfficiency::WriteDebugInfo(const int eyeId, const int simTelId, const double simTimeOffset,
                                            const FdConstants::LightSource component,
                                            const utl::TabulatedFunctionErrors& totalPixelTrace,
                                            const utl::TraceD& totalWlTrace,
                                            const utl::TabulatedFunctionErrors& realApLightFlux,
                                            const utl::TabulatedFunctionErrors& telEff)
{
  ostringstream pixname;
  pixname << "pixHist_e" << eyeId << "_t" << simTelId << "_c" << component;
  vector<double> pixHistX, pixHistY;
  for (unsigned int i = 0; i < totalPixelTrace.GetNPoints(); ++i) {
    pixHistX.push_back(totalPixelTrace.GetX(i) + simTimeOffset);
    pixHistY.push_back(totalPixelTrace.GetY(i));
  }
  TGraph* const pixHist =
    new TGraph(totalPixelTrace.GetNPoints(), &pixHistX.front(), &pixHistY.front());
  pixHist->SetNameTitle(pixname.str().c_str(), pixname.str().c_str());
  pixHist->Write();

  ostringstream apname;
  apname << "apHist_e" << eyeId << "_t" << simTelId << "_c" << component;
  TH1D* const apHist =
    new TH1D(apname.str().c_str(), apname.str().c_str(),
             totalWlTrace.GetSize(), simTimeOffset,
             simTimeOffset + totalWlTrace.GetBinning()*totalWlTrace.GetSize());
  for (unsigned int i = 1; i <= totalWlTrace.GetSize(); ++i)
    apHist->SetBinContent(i, totalWlTrace.At(i-1));
  apHist->Write();

  ostringstream realapname;
  realapname << "realApLight_e" << eyeId << "_t" << simTelId << "_c" << component;
  TGraph* const realAp =
    new TGraph(realApLightFlux.GetNPoints(),
               &realApLightFlux.XFront(), &realApLightFlux.YFront());
  realAp->SetNameTitle(realapname.str().c_str(), realapname.str().c_str());
  realAp->Write();

  ostringstream relname;
  relname << "relHist_e" << eyeId << "_t" << simTelId << "_c" << component;
  vector<double> relHistX, relHistXe, relHistY, relHistYe;
  for (unsigned int i = 0; i < telEff.GetNPoints(); ++i) {
    relHistX.push_back(telEff.GetX(i) + simTimeOffset);
    relHistXe.push_back(telEff.GetXErr(i));
    relHistY.push_back(telEff.GetY(i));
    relHistYe.push_back(telEff.GetYErr(i));
  }
  TGraphErrors* const relHist =
    new TGraphErrors(telEff.GetNPoints(), &relHistX.front(), &relHistY.front(),
                                          &relHistXe.front(), &relHistYe.front());
  relHist->SetNameTitle(relname.str().c_str(), relname.str().c_str());
  relHist->Write();
}

/*************************************************************************/
utl::TraceD
FdLightCollectionEfficiency::CalculateTelescopeTraceSum(const fevt::Telescope& tel,
                                                        const fevt::FdConstants::LightSource lightSource,
                                                        const atm::Atmosphere::EmissionMode /*lightType*/)
{
  const fevt::TelescopeSimData& telSimData = tel.GetSimData();
  utl::TraceD totalWlTrace;
  if (!telSimData.HasPhotonTrace(lightSource))
    return totalWlTrace;

  const fdet::FDetector& fdetector = det::Detector::GetInstance().GetFDetector();
  const fdet::Telescope& detTel = fdetector.GetEye(tel.GetEyeId()).GetTelescope(tel.GetId());

  //const atm::Atmosphere& atmo = det::Detector::GetInstance().GetAtmosphere();
  //const vector<double>& wavelengths = atmo.GetWavelengths(lightType);

  const double refWl = fdetector.GetReferenceLambda();
  //const TabulatedFunction& modelEfficiency = detTel.GetModelRelativeEfficiency(refWl);
  const TabulatedFunction& measuredEfficiency = detTel.GetMeasuredRelativeEfficiency();
  //const double refModelEfficiency = modelEfficiency.Y(refWl);
  //const double refMeasuredEfficiency = measuredEfficiency.Y(refWl);

  // Note: Taking the Q-Eff from a single pixel may be a hack, but should *currently*
  //       be reliable: Here, the pixQEff is only used for determining the min/max
  //       supported wavelength. This whole min/max logi mimics what is done elsewhere
  //       in the same context, so we have to do the same in order to remain consistent.
  const TabulatedFunction& pixelQEff = detTel.GetPixel(1).GetQEfficiency();
  //const double refPixelQEff = pixelQEff.Y(refWl);

  // The efficiency TabulatedFunctions don't cover arbitrary wavelengths
  const double minWl =
    std::max(std::max(detTel.GetModelMinWavelength(), measuredEfficiency.GetX(0)), pixelQEff.GetX(0));
  const double maxWl =
    std::min(std::min(detTel.GetModelMaxWavelength(), measuredEfficiency.GetX(measuredEfficiency.GetNPoints()-1)),
                      pixelQEff.GetX( pixelQEff.GetNPoints() - 1));

  bool firstWl = true;

  const double diaArea = detTel.GetDiaphragmArea();
  //cout << "diaArea: " << diaArea << endl;

  for (TelescopeSimData::ConstPhotonTraceIterator traceIter = telSimData.PhotonTracesBegin(lightSource),
       end = telSimData.PhotonTracesEnd(lightSource);
       traceIter != end; ++traceIter) {
    const utl::TraceD& wlTrace = traceIter->GetTrace();
    //cout << "Trace " << traceIter->GetLabel() << " length: " << wlTrace.GetSize() << endl;
    // Note: overriding the wavelengths with the reference wavelength since
    // we manually override the wavelength sampling in the ShowerPhotonGenerator
    // TODO understand why this is necessary?
    const double wl = refWl;//wavelengths[traceIter->GetLabel()];

    // Replicate the behaviour of the TelescopeSimulatorKG
    if (wl < minWl || wl > maxWl) {
      //cout << "Skipping wl. min="<<minWl << " wl=" << wl << " max=" << maxWl << endl;
      continue;
    }

    if (firstWl) {
      // Set the binning/start/stop
      totalWlTrace.SetBinning(wlTrace.GetBinning());
      totalWlTrace.SetStart(wlTrace.GetStart());
      totalWlTrace.SetStop(wlTrace.GetStop());

      for (unsigned int i = 0; i < wlTrace.GetSize(); ++i)
        totalWlTrace.PushBack(0);

      firstWl = false;
    } else {
      // Note: This is likely not necessary as the start/end times should be the same for all
      //       wl's, but let's be defensive for now, it doesn't hurt.
      //       Also, the trace length should be quite the same.
      if (wlTrace.GetStart() < totalWlTrace.GetStart())
        totalWlTrace.SetStart(wlTrace.GetStart());
      if (wlTrace.GetStop() > totalWlTrace.GetStop())
        totalWlTrace.SetStop(wlTrace.GetStop());
      if (wlTrace.GetSize() != totalWlTrace.GetSize())
        throw OutOfBoundException("The wavelength-specific traces differ in their length!");
    }

    for (unsigned int i = 0; i < wlTrace.GetSize(); ++i)
      totalWlTrace[i] += wlTrace[i] * diaArea;

  } // end for each wavelength

  //cout << "total wl trace size=" << totalWlTrace.GetSize() << endl;
  return totalWlTrace;
}


/*************************************************************************/
utl::TabulatedFunctionErrors
FdLightCollectionEfficiency::CalculatePixelTraceSum(const fevt::Telescope& tel,
                                                    const fevt::FdConstants::LightSource component,
                                                    const double simTimeOffset,
                                                    const vector<vector<unsigned int> >& pixelsInZeta,
                                                    const utl::TabulatedFunctionErrors& realApLight)
{
  utl::TabulatedFunctionErrors totalPixelTrace;

  const fdet::FDetector& theFDet = det::Detector::GetInstance().GetFDetector();
  const fdet::Telescope& detTel = theFDet.GetTelescope( tel );
  const unsigned int telId = tel.GetId();

  const unsigned int nPoints = realApLight.GetNPoints();
  for (unsigned int iTimeBin = 0; iTimeBin < nPoints; ++iTimeBin) {
    const vector<unsigned int>& thisTimeZetaPixels = pixelsInZeta[iTimeBin];
    const double timeBinCenter = realApLight.GetX(iTimeBin) - simTimeOffset;
    const double dtHalf = realApLight.GetXErr(iTimeBin); // half time bin width
    double timeBinSignal = 0.;
    double timeBinSignalVar = 0.;

    const unsigned int nPixels =
      fUseZeta ? thisTimeZetaPixels.size() : detTel.GetLastPixelId()-detTel.GetFirstPixelId()+1;

    for (unsigned int iPixId = 0; iPixId < nPixels; ++iPixId) {
      const unsigned int pixId = fUseZeta ? thisTimeZetaPixels[iPixId] : iPixId + 1;
      if (!tel.HasPixel(pixId, ComponentSelector::eInDAQ))
        continue;

      // Fetch or calculate pixel efficiency modification
      double pixEffMod = fLightConverter->GetConversionConstant(telId, pixId);

      const fevt::Pixel& pix = tel.GetPixel(pixId);
      const fevt::PixelSimData& pixSim = pix.GetSimData();
      if (!pixSim.HasPhotonTrace(component)) {
//#warning The light-at-pixel information should be available at this point. Skipping the light component seems wrong.
        //cout << "PixelSimData for pixel " << pix.GetId() << " of telescope " << tel.GetId() << " does not have light component " << component << ". Skipping." << endl;
        continue;
      }
      const utl::TraceD& pixTrace   = pixSim.GetPhotonTrace(component);
      const double pixSignalInBin =
        CalcTraceBinContent<double>(pixTrace, timeBinCenter - dtHalf, timeBinCenter + dtHalf);
      timeBinSignal += pixSignalInBin * pixEffMod;

      const utl::TraceD& pixWSTrace = pixSim.GetPhotonWeightSquareTrace(component);
      const double pixSignalVarInBin =
        CalcTraceBinContent<double>(pixWSTrace, timeBinCenter - dtHalf, timeBinCenter + dtHalf);
      timeBinSignalVar += pixSignalVarInBin * pixEffMod*pixEffMod;
    } // end for each zeta pixel

    totalPixelTrace.PushBack(timeBinCenter, dtHalf, timeBinSignal, sqrt(timeBinSignalVar));
  } // end for each aperture time bin

  return totalPixelTrace;
}


/*************************************************************************/
fwk::VModule::ResultFlag
FdLightCollectionEfficiency::DoInitialReconstruction(evt::Event& event, std::map<int, evt::ShowerFRecData*>& dataSets)
{
  SubModule& profRecMod = fSubModulesList.find("FdProfileReconstructorKG")->second;

  FdProfileReconstructorKG::FdProfileReconstructor& reconstructor =
    dynamic_cast<FdProfileReconstructorKG::FdProfileReconstructor&>(profRecMod.GetModule());

  const bool propGeo = reconstructor.GetPropagateGeometryErrors();
  const bool propAtm = reconstructor.GetPropagateAtmUncertainties();
  reconstructor.SetPropagateGeometryErrors(false);
  reconstructor.SetPropagateAtmUncertainties(false);

  fwk::VModule::ResultFlag status = profRecMod.Run(event);

  reconstructor.SetPropagateGeometryErrors(propGeo);
  reconstructor.SetPropagateAtmUncertainties(propAtm);

  // Now, do some minor validation to make sure the event has been reconstructed
  // in a *remotely* physical way.

  // At the same time, create a map of eye-id => ShowerFRecData as the set
  // of data to operate on down the road
  for (fevt::FEvent::EyeIterator iEye = event.GetFEvent().EyesBegin(fevt::ComponentSelector::eHasData),
       end = event.GetFEvent().EyesEnd(fevt::ComponentSelector::eHasData);
       iEye != end; ++iEye)
  {
    // Hierarchically check all the necessary structures
    if (!iEye->HasRecData()) {
      if (fVerbosity > 0)
        cout << "    Skipping eye " << iEye->GetId() << ": no EyeRecData.\n";
      continue;
    }

    EyeRecData& eyeRecData = iEye->GetRecData();
    if (!eyeRecData.HasFRecShower()) {
      if (fVerbosity > 0)
        cout << "    Skipping eye " << iEye->GetId() << ": no ShowerFRecData.\n";
      continue;
    }

    ShowerFRecData& frecData = eyeRecData.GetFRecShower();
    if (!frecData.HasGHParameters() ||
        !frecData.HasLongitudinalProfile() ||
        !frecData.HasEnergyDeposit()) {
      if (fVerbosity > 0)
        cout << "    Skipping eye " << iEye->GetId() << ": no profile reconstruction.\n";
      continue;
    }

    const evt::VGaisserHillasParameter& ghParam = frecData.GetGHParameters();
    const evt::GaisserHillas4Parameter* gh4 = dynamic_cast<const evt::GaisserHillas4Parameter*>(&ghParam);
    if (!gh4) {
      if (fVerbosity > 0)
        cout << "    Skipping eye " << iEye->GetId() << ": no profile GH parameters after profile reconstruction.\n";
      continue;
    }

    if (gh4->GetXMax() < 0. || gh4->GetNMax() < 0.) {
      if (fVerbosity > 0)
        cout << "    Skipping eye " << iEye->GetId() << ": unphysical profile reconstruction.\n";
      continue;
    }

    dataSets[iEye->GetId()] = &frecData;
  } // end of eye loop

  return status;
}


/*************************************************************************/
std::pair<fwk::VModule::ResultFlag, unsigned int>
FdLightCollectionEfficiency::RunPhotonGenerationRayTracingLoop(evt::Event& event,
                                                               evt::Event& realEvent,
                                                               unsigned int nRayTracingIterations,
                                                               const fevt::FdConstants::LightSource lightSource)
{
  fevt::FEvent& theFEvent = event.GetFEvent();

  // Prepare ShowerPhotonGeneratorOG configuration
  SubModule& photGenMod = fSubModulesList.find("ShowerPhotonGeneratorOG")->second;
  ShowerPhotonGeneratorOG::ShowerPhotonGenerator& photGen =
    dynamic_cast<ShowerPhotonGeneratorOG::ShowerPhotonGenerator&>(photGenMod.GetModule());

  const unsigned int maxNRayTraceOld = photGen.GetMaxNRayTrace(); // remember
  const unsigned int minNRayTraceOld = photGen.GetMinNRayTrace(); // remember
  const double extraFactor = photGen.GetExtraRayTraceFactor(); // remember
  const int sourceSelection = photGen.GetLightSourceSelection(); // remember
  const bool useOnlyReferenceWl = photGen.GetUseOnlyReferenceWavelength(); // remember
  photGen.SetMaxNRayTrace(fMaxNRayTracePerBin); // overwrite
  photGen.SetMinNRayTrace(fMinNRayTracePerBin); // overwrite
  photGen.SetUseOnlyReferenceWavelength(true); // overwrite
  photGen.SetExtraRayTraceFactor(fExtraRayTraceFactor); // overwrite
  photGen.SetLightSourceSelection(lightSource); // overwrite

  // Prepare TelescopeSimulatorKG configuration
  SubModule& telSimMod = fSubModulesList.find("TelescopeSimulatorKG")->second;
  TelescopeSimulatorKG::TelescopeSimulator& telSim =
    dynamic_cast<TelescopeSimulatorKG::TelescopeSimulator&>(telSimMod.GetModule());

  const bool storeLightAtPixels = telSim.StoreLightComponentsAtPixels(); // remember
  const bool telSimVerbosity = telSim.GetVerbosity(); // remember

  telSim.SetStoreLightComponentsAtPixels(true); // overwrite
  telSim.SetVerbosity(1); // overwrite

  if (fVerbosity > 0)
    cout << "\n  Performing at least " << nRayTracingIterations
         << " ray-tracing iteration(s):" << endl;

  typedef std::list<utl::TabulatedFunctionErrors> TelLCEfficiencies;
  typedef std::map<unsigned int, TelLCEfficiencies> EyeLCEfficiences;
  typedef std::map<unsigned int, EyeLCEfficiences> EventLCEfficiencies;
  EventLCEfficiencies lcEfficiencies;

  fwk::VModule::ResultFlag status = eFailure;
  for (unsigned int iRayTrace = 1; iRayTrace <= (unsigned int)nRayTracingIterations; ++iRayTrace) {

    bool needRepeat = false;

    if (fVerbosity > 0)
      cout << "\n  Ray-tracing iteration " << iRayTrace << " ..." << endl;

//#warning Temporary measure until I can reproduce the ShowerPhotonGenerator exceptions
    try { // FIXME this is a temporary measure in order not to lose a whole job for absurd events
      status = photGenMod.Run(event); // RUN
      if (status != eSuccess)
        break;
    }
    catch (...) {
      ostringstream err;
      err << "ShowerPhotonGenerator threw an exception."
             " Skipping the event for lack of better options...";
      ERROR(err);
      status = eContinueLoop;
      break;
    }

    // Note: we won't catch exceptions from the TelescopeSimulator.
    //       There is no reason to throw any even for a fairly malformed
    //       event.
    status = telSimMod.Run(event); // RUN
    if (status != eSuccess)
      break;

    // Now, calculate the sum of photon weights at aperture (but using their
    // 3D-shifted time), then empty the photon traces
    // (TelescopeSimulatorKG doesn't do the latter, currently)

    // loop eyes
    for (fevt::FEvent::EyeIterator iEye = theFEvent.EyesBegin(fevt::ComponentSelector::eInDAQ),
         end = theFEvent.EyesEnd(fevt::ComponentSelector::eInDAQ);
         iEye != end; ++iEye)
    {
      if (iEye->GetStatus() == fevt::ComponentSelector::eDeSelected)
        continue;

      fevt::Eye& eyeEvent = *iEye;

      const unsigned int eyeId = eyeEvent.GetId();

      // loop telescopes for calculating time-corrected aperture light trace
      for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ),
           end = eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
           iTel != end; ++iTel)
      {
        const unsigned int telId = iTel->GetId();

        fevt::TelescopeSimData& simData = iTel->GetSimData();
        if (simData.HasRayTracedPhotonTrace()) {
          // Get the corrected aperture light trace
          utl::TraceD& apTrace = (*fTimeCorrectedApertureTraces)[telId];
          // Add photons to the corrected ap. light trace
          CorrectedApLightCalculator calc(simData,
                                          apTrace,
                                          //eyeId, // unused. LN.
                                          telId,
                                          fVerbosity);
          calc.Run();
        }
      } // end for telescope

      if (fStopCondition != eFixed)
        CalculateEfficiency(event, realEvent, eyeId, iRayTrace, lightSource);

      EyeLCEfficiences& eyeLCEff = lcEfficiencies[eyeId];

      // loop telescopes
      for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ),
           end = eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
           iTel != end; ++iTel)
      {
        const unsigned int telId = iTel->GetId();

        fevt::TelescopeSimData& simData = iTel->GetSimData();
        simData.ClearPhotons();

        if (fStopCondition == eBootstrap) {
          // If we're bootstrapping, we calculate the efficiency for each iteration,
          // then use the list of efficiencies to arrive at a final results.
          // That means we need to clean up the intermediate results or we'll get
          // a wrong end result.
          if (simData.HasRayTracedPhotonTrace())
            simData.ClearRayTracedPhotonTrace();

          // Clear the simulated photon traces @ pixel
          for (fevt::Telescope::PixelIterator iPix = iTel->PixelsBegin(fevt::ComponentSelector::eInDAQ),
               end = iTel->PixelsEnd(fevt::ComponentSelector::eInDAQ);
               iPix != end; ++iPix)
          {
            iPix->GetSimData().ClearPhotonTraces();
            iPix->GetSimData().ClearPhotonWeightSquareTraces();
          }

          TelLCEfficiencies& telLCEff = eyeLCEff[telId];
          if (!realEvent.GetFEvent().HasEye(eyeId, ComponentSelector::eHasData))
            continue;

          fevt::Eye& reye = realEvent.GetFEvent().GetEye(eyeId, ComponentSelector::eHasData);
          if (!reye.HasTelescope(telId, ComponentSelector::eHasData))
            continue;

          fevt::Telescope& rtel = reye.GetTelescope(telId, ComponentSelector::eHasData);
          if (!rtel.HasRecData())
            continue;

          TelescopeRecData& recData = reye.GetTelescope(telId, ComponentSelector::eHasData).GetRecData();
          telLCEff.push_back( recData.GetLightCollectionEfficiency().GetTabulatedFunctionErrors(lightSource) );

          // Check whether we arrived at the end of the mandatory number of
          // iterations. If so, recalculate the net efficiency
          if (iRayTrace == nRayTracingIterations) {
            cout << "Bootstrapping uncertainties..." << endl;
            Bootstrapper bs(telLCEff, fVerbosity);
            // Repeat if we haven't met the uncertainty target yet
            if (!needRepeat &&
                !bs.MaxRelUncertaintyBelowThreshold(fTargetUncertainty, fMinRelevantEfficiency)) {
              needRepeat = true;
              continue;
            }
            utl::MultiTabulatedFunctionErrors& mult = recData.GetLightCollectionEfficiency();
            mult.GetTabulatedFunctionErrors(lightSource) = bs.GetMean();
          }
        } // end if using bootstrap method
      } // end for telescopes
    } // end for eyes


    // If we're aiming for a certain upper limit on the uncertainty of the efficiency,
    // we predict the number of photons remaining to be simulated
    if ((fStopCondition == eTargetUncertainty
         || fStopCondition == eSmallEnough)
        && iRayTrace == 1)
    {
      double targetRelUncertainty = 0.;
      if (fStopCondition == eSmallEnough) {
        // Get uncertainty from l@aperture errors
        targetRelUncertainty =
          CalculateTargetRelUncertainty(realEvent.GetFEvent(),
                                        event.GetFEvent(), lightSource,
                                        fMinRelevantEfficiency,
                                        fTargetWorstCaseLAtApUncertaintyChange);
        // Make sure we don't go crazy on the uncertainty...
        if (targetRelUncertainty < fLowerUncertaintyLimit)
          targetRelUncertainty = fLowerUncertaintyLimit;
      }
      else
        targetRelUncertainty = fTargetUncertainty;

      const double additionalIter =
          CalculateAdditionalIterationsToMeetTarget(realEvent.GetFEvent(),
                                                    event.GetFEvent(), lightSource,
                                                    targetRelUncertainty, fMinRelevantEfficiency);
      if (additionalIter > 0.) {
        needRepeat = true;
        unsigned int addIter = (unsigned int)additionalIter;
        addIter = (additionalIter - (double)addIter > 0.) ? addIter+1 : addIter; // ceil
        ostringstream msg;
        msg << "Determined that at least " << additionalIter
            << " ray-tracing iterations are necessary to reach uncertainty target."
            << endl;
        INFO(msg);
        nRayTracingIterations += addIter-1; // one iteration added by virtue of needRepeat
      }
      else
        INFO("Met uncertainty target for this light component for all bins.");
    } // end if in eTargetUncertainty mode and at first iteration


    // Now check whether we need to keep going according to the uncertainty target and the
    // calculated uncertainties from bootstrapping.
    if (needRepeat) {
      if (iRayTrace == fMaxIterations) {
        cout << "Would need an additional iteration according to stop "
             << "criterion, but we hit the hard upper limit. Stopping iteration." << endl;
        nRayTracingIterations = iRayTrace;
        break;
      }
      else {
        cout << "Performing additional iteration to meet uncertainty target..." << endl;
        nRayTracingIterations++;
      }
    }
    else if (fStopCondition != eFixed && iRayTrace == fMaxIterations) {
      cout << "Would need an additional iteration according to stop "
           << "criterion, but we hit the hard upper limit. Stopping iteration." << endl;
      nRayTracingIterations = iRayTrace;
      break;
    }
  } // end for nRayTrace iterations

  // Reset ShowerPhotonGeneratorOG configuration
  photGen.SetMaxNRayTrace(maxNRayTraceOld); // reset
  photGen.SetMinNRayTrace(minNRayTraceOld); // reset
  photGen.SetExtraRayTraceFactor(extraFactor); // reset
  photGen.SetUseOnlyReferenceWavelength(useOnlyReferenceWl); // reset
  if (sourceSelection < 0)
    photGen.ResetLightSourceSelection(); // reset
  else
    photGen.SetLightSourceSelection((fevt::FdConstants::LightSource)sourceSelection); // reset

  // Reset TelescopeSimulatorKG configuration
  telSim.SetStoreLightComponentsAtPixels(storeLightAtPixels); // reset
  telSim.SetVerbosity(telSimVerbosity); // reset

  return std::make_pair(status, nRayTracingIterations);
}


/*************************************************************************/
double
FdLightCollectionEfficiency::CalculateAdditionalIterationsToMeetTarget(const fevt::FEvent& theRealFEvent,
                                                                       const fevt::FEvent& theSimFEvent,
                                                                       const fevt::FdConstants::LightSource& lightSource,
                                                                       const double targetRelUncertainty,
                                                                       const double minRelevantEfficiency)
{
  double minNIterationsRequired = 0.;
  double photonsRayTracedAtMinimum = 0;
  double efficiencyAtMinimum = 0;
  double efficiencyUncertaintyAtMinimum = 0;
  // loop eyes
  for (fevt::FEvent::ConstEyeIterator iEye = theSimFEvent.EyesBegin(fevt::ComponentSelector::eInDAQ),
       end = theSimFEvent.EyesEnd(fevt::ComponentSelector::eInDAQ);
       iEye != end; ++iEye)
  {
    if (iEye->GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;
    if (!theRealFEvent.HasEye(iEye->GetId(), fevt::ComponentSelector::eHasData))
      continue;

    const fevt::Eye& eyeEvent = *iEye;
    const fevt::Eye& realEye = theRealFEvent.GetEye(iEye->GetId(), fevt::ComponentSelector::eHasData);

    // loop telescopes
    for (fevt::Eye::ConstTelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ),
         end = eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
         iTel != end; ++iTel)
    {
      if (!iTel->HasSimData())
        continue;
      const fevt::TelescopeSimData& simData = iTel->GetSimData();

      if (!realEye.HasTelescope(iTel->GetId(), fevt::ComponentSelector::eHasData))
        continue;
      const fevt::Telescope& realTel = realEye.GetTelescope(iTel->GetId(), fevt::ComponentSelector::eHasData);

      if (!realTel.HasRecData())
        continue;
      const TelescopeRecData& recData = realTel.GetRecData();

      if (!recData.HasLightCollectionEfficiency()
          || !recData.GetLightCollectionEfficiency().HasLabel(lightSource))
        continue;

      const utl::TabulatedFunctionErrors& lcEff = recData.GetLightCollectionEfficiency().GetTabulatedFunctionErrors(lightSource);

      if (!simData.HasRayTracedPhotonTrace())
        continue;
      const utl::TraceI& nRayTraced = simData.GetRayTracedPhotonTrace();

      const unsigned int nPoints = lcEff.GetNPoints();
      for (unsigned int i = 0; i < nPoints; ++i) {
        if (lcEff.GetY(i) >= minRelevantEfficiency) {
          if (nRayTraced[i] == 0)
            continue;
          const double yerr = lcEff.GetYErr(i);
          const double nPhotNow = nRayTraced[i];
          if (yerr > targetRelUncertainty) {
            // sigma > sigma_t => sigma/sqrt(n) = sigma_t => n = (sigma/sigma_t)^2
            // s_mean_now = s / sqrt(n_now); s_mean_target = s / sqrt(n_target)
            // => s_mean_target * sqrt(n_target) = s_mean_now * sqrt(n_now)
            // => n_target = (s_mean_now / s_mean_target)^2 * n_now
            // Note: targetRelUncertainty is a relative uncertainty!
            const double nTarget = pow(yerr/(targetRelUncertainty*lcEff.GetY(i)), 2.) * nPhotNow;
            const double thisMinNPhot = nTarget - nPhotNow;
            const double thisMinNIter = thisMinNPhot / (double)nPhotNow; // one iteration yields nPhotNow photons in this bin
            if (thisMinNIter > minNIterationsRequired) {
              minNIterationsRequired         = thisMinNIter;
              photonsRayTracedAtMinimum      = nPhotNow;
              efficiencyAtMinimum            = lcEff.GetY(i);
              efficiencyUncertaintyAtMinimum = yerr;
            }
          }
        }
      } // end foreach lceff bin

    } // end for telescopes
  } // end for eyes

  if (fVerbosity > 0) {
    ostringstream msg;
    msg << "Calculated that we'll need " << minNIterationsRequired
        << " more iterations to meet uncertainty target.";
    if (fVerbosity > 1) {
      msg << " In the bin of the minimum, we ray-traced " << photonsRayTracedAtMinimum
          << " so far, therefore, we need "
          << (minNIterationsRequired+1)*photonsRayTracedAtMinimum
          << " photons in total to meet target. In this bin, the efficiency was: "
          << efficiencyAtMinimum << " +/- " << efficiencyUncertaintyAtMinimum;
    }
    INFO(msg);
  }

  return minNIterationsRequired;
}


/*************************************************************************/
double
FdLightCollectionEfficiency::CalculateTargetRelUncertainty(const fevt::FEvent& theRealFEvent,
                                                           const fevt::FEvent& theSimFEvent,
                                                           const fevt::FdConstants::LightSource& lightSource,
                                                           const double minRelevantEfficiency,
                                                           const double maxUncertaintyChangeFraction)
{
  unsigned int worstTel = 0;
  unsigned int worstBin = 0;
  double worstTime      = 0.;
  double worstCaseRelUncertainty = 1.e9;
  // loop eyes
  for (fevt::FEvent::ConstEyeIterator iEye = theSimFEvent.EyesBegin(fevt::ComponentSelector::eInDAQ),
       end = theSimFEvent.EyesEnd(fevt::ComponentSelector::eInDAQ);
       iEye != end; ++iEye)
  {
    if (iEye->GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;
    if (!theRealFEvent.HasEye(iEye->GetId(), fevt::ComponentSelector::eHasData))
      continue;

    const fevt::Eye& eyeEvent = *iEye;
    const fevt::Eye& realEye = theRealFEvent.GetEye(eyeEvent.GetId(), fevt::ComponentSelector::eHasData);

    // loop telescopes
    for (fevt::Eye::ConstTelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ),
         end = eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
         iTel != end; ++iTel)
    {
      if (!iTel->HasSimData())
        continue;

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

      if (!realEye.HasTelescope(telId, fevt::ComponentSelector::eHasData))
        continue;
      const fevt::Telescope& realTel = realEye.GetTelescope(telId, fevt::ComponentSelector::eHasData);

      if (!realTel.HasRecData())
        continue;
      const TelescopeRecData& recData = realTel.GetRecData();

      if (!recData.HasLightFlux(fevt::FdConstants::eTotal)
          || !recData.HasLightCollectionEfficiency()
          || !recData.GetLightCollectionEfficiency().HasLabel(lightSource))
        continue;

      const utl::TabulatedFunctionErrors& lcEff = recData.GetLightCollectionEfficiency().GetTabulatedFunctionErrors(lightSource);
      const utl::TabulatedFunctionErrors& lAtAp = recData.GetLightFlux(fevt::FdConstants::eTotal); // FIXME eTotal light is a worst-case?

      const unsigned int nPoints = lcEff.GetNPoints();
      assert(nPoints == lAtAp.GetNPoints());

      for (unsigned int i = 0; i < nPoints; ++i) {
        const double eff = lcEff.GetY(i);
        if (eff >= minRelevantEfficiency) {
          const double light    = fabs(lAtAp.GetY(i));
          const double lightErr = lAtAp.GetYErr(i);
          const double requiredRelEff = sqrt( pow(maxUncertaintyChangeFraction+1, 2)-1 ) * lightErr / light;
          if (requiredRelEff < worstCaseRelUncertainty) {
            worstTel                = telId;
            worstBin                = i;
            worstTime               = lAtAp.GetX(i);
            worstCaseRelUncertainty = requiredRelEff;
              cout << telId << " " << i << " " << maxUncertaintyChangeFraction << " "
              << lightErr << " " << light << " " << lightErr/light << " " << requiredRelEff
              << endl;
          }
        }
      } // end foreach lceff bin

    } // end for telescopes
  } // end for eyes

  if (fVerbosity > 0) {
    ostringstream msg;
    msg << "Calculated worst-case required relative uncertainty on efficiency is "
        << worstCaseRelUncertainty << ". ";
    if (worstCaseRelUncertainty < fLowerUncertaintyLimit) {
      msg << " It is smaller than the lower bound for the rel. unc. ("
          << fLowerUncertaintyLimit << ") and thus, we'll use the latter instead. ";
    }
    msg << "This will be the basis for calculating the "
           "required number of iterations.";
    if (fVerbosity > 1) {
      msg << " The most stringent requirement was found in tel. " << worstTel
          << " for bin " << worstBin << " at time " << worstTime << ".";
    }
    INFO(msg);
  }

  return worstCaseRelUncertainty;
}


/*************************************************************************/
utl::TabulatedFunctionErrors
FdLightCollectionEfficiency::CalcTraceDivision(const utl::TabulatedFunctionErrors& pixelTrace,
                                               const utl::TraceD& apertureTrace,
                                               const utl::TraceI& /*nRayTracedTrace*/,
                                               const double nRayTracingIterations)
{
  const unsigned int nBins = pixelTrace.GetNPoints();
  vector<double> vx(nBins), vxerr(nBins), vy(nBins), vyerr(nBins)/*, n(nBins)*/;

  for (unsigned int i = 0; i < nBins; ++i) {
    const double binCenter = pixelTrace.GetX(i);
    const double halfBinWidth = pixelTrace.GetXErr(i);
    const double binStart = binCenter - halfBinWidth;
    const double binEnd = binCenter + halfBinWidth;

    const double pixelSignal    = pixelTrace.GetY(i);
    const double pixelSignalErr = pixelTrace.GetYErr(i);
    assert(isfinite(pixelSignal));
    assert(isfinite(pixelSignalErr));

    const double apertureSignal = CalcTraceBinContent(apertureTrace, binStart, binEnd) * nRayTracingIterations;
    assert(isfinite(apertureSignal));
    //const double nRayTraced = CalcTraceBinContent(nRayTracedTrace, binStart, binEnd);

    const double x = binCenter;
    if (apertureSignal >= 1.e-9) {
      const double y = pixelSignal / apertureSignal;
      assert(isfinite(y));
      vy[i]    = y;
      vyerr[i] = pixelSignalErr / apertureSignal;
    }
    else {
      vy[i]    = 0.; // default to zero efficiency if there was no light@aperture
      vyerr[i] = 1.e3;
    }

    vx[i]    = x;
    vxerr[i] = halfBinWidth;
    //n[i]     = nRayTraced;
    //efficiency.PushBack(x, xerr, 0., 0.);
  } // end foreach bin

  utl::TabulatedFunctionErrors division;
  for (unsigned int i = 0; i < nBins; ++i) {
    const double yerr = vyerr[i];
    assert(isfinite(vx[i]));
    assert(isfinite(vxerr[i]));
    assert(isfinite(vy[i]));
    assert(isfinite(yerr));
    division.PushBack(vx[i], vxerr[i], vy[i], yerr);
  }

  return division;
}

#undef DEBUG

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: