CorrectedApLightCalculator.cc

Go to the documentation of this file.

#include "CorrectedApLightCalculator.h"

#include <vector>
#include <list>
#include <limits>
#include <sstream>
#include <iostream>

#include <evt/Event.h>

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

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

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>

using namespace FdLightCollectionEfficiencyKG;
using namespace utl;
using namespace evt;
using namespace fevt;
using namespace std;


/*************************************************************************/
CorrectedApLightCalculator::CorrectedApLightCalculator(const fevt::TelescopeSimData& telSimData,
                                                       utl::TraceD& outputApTrace,
                                                       //const unsigned int eyeId, // unused. LN.
                                                       const unsigned int telId,
                                                       const int verbosity) :
                                                       fTelSimData(telSimData),
                                                       fOutputApTrace(outputApTrace),
                                                       //fEyeId(eyeId), // unused. LN.
                                                       fTelId(telId),
                                                       fVerbosity(verbosity)
{
}

/*************************************************************************/
void
CorrectedApLightCalculator::Run()
{
  if (fVerbosity > 1) {
    ostringstream msg;
    msg << "Calculating time-corrected aperture light trace for telescope "
        << fTelId << " ...";
    INFO(msg);
  }

  const utl::TraceI& rtPhotTrace = fTelSimData.GetRayTracedPhotonTrace();

  // initialize
  if (fOutputApTrace.GetSize() != rtPhotTrace.GetSize()) {
    fOutputApTrace.SetBinning(rtPhotTrace.GetBinning());
    fOutputApTrace.SetStart(rtPhotTrace.GetStart());
    fOutputApTrace.SetStop(rtPhotTrace.GetStop());

    for (unsigned int i = 0; i < rtPhotTrace.GetSize(); ++i)
      fOutputApTrace.PushBack(0.);
  }
  //const fdet::FDetector& theFDet = det::Detector::GetInstance().GetFDetector();
  //const fdet::Telescope& detTel  = theFDet.GetEye(fEyeId).GetTelescope(fTelId);
  //const double diaArea = detTel.GetDiaphragmArea();

  const int telTraceNBins = fOutputApTrace.GetSize();
  const double traceBinWidth = fOutputApTrace.GetBinning();

  double sumOutOfBounds = 0.;
  double sumInBounds = 0.;

  // Loop over all photons and add them to the corrected aperture
  // light trace
  for (fevt::TelescopeSimData::ConstPhotonIterator iPhoton = fTelSimData.PhotonsBegin(),
       end = fTelSimData.PhotonsEnd(); iPhoton != end; ++iPhoton)
  {
    const utl::Photon& p = *iPhoton;
    const double time = p.GetTime().GetInterval();
    const int bin = int(time/traceBinWidth);

    if (bin < 0 || bin >= int(telTraceNBins)) { // FIXME think about this
      if (fVerbosity > 2) {
        ostringstream err;
        err << "Bin out of range. bin=" << bin
            << " weight=" << p.GetWeight();
        WARNING(err);
      }
      sumOutOfBounds += p.GetWeight();
      continue;
    }
    sumInBounds += p.GetWeight();

    const double photonWeight = p.GetWeight();
    //fOutputApTrace[bin] += photonWeight * diaArea;
    fOutputApTrace[bin] += photonWeight;
  } // end for photons

  if (fVerbosity > 1) {
    ostringstream msg;
    msg << "  For this ray-tracing iteration: Total sum of photon weights in range: " << sumInBounds
        << " out of range: " << sumOutOfBounds << endl;
    INFO(msg);
  }
}


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