FdNoiseAdder.cc

Go to the documentation of this file.

#include "FdNoiseAdder.h"
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/Reader.h>
#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>
#include <utl/RandomEngine.h>
#include <CLHEP/Random/RandGauss.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 <evt/Event.h>
#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <fevt/EyeHeader.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <utl/TraceAlgorithm.h>

#include <sstream>

using CLHEP::RandGauss;
using namespace std;
using namespace utl;
using namespace fwk;
using namespace evt;
using namespace fevt;


namespace FdNoiseAdderKG {


  double
  GetCalibConst(const fevt::Pixel& pixel) {
    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();
    // Get the calibration constant for this channel
    const int channelId = detPixel.GetChannelId();
    const fdet::Pixel& channel =
      detFD.GetEye(eyeId).GetTelescope(telId).GetPixel(channelId);
    return channel.GetEndToEndCalibrationAtReferenceWavelength();
  }

  double
  CalcMeanVariance(const fevt::Eye& eye)
  {
      double meanVariance = 0;
      int nVariance = 0;
      for (fevt::Eye::ConstTelescopeIterator telIter =
             eye.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
           telIter != eye.TelescopesEnd(fevt::ComponentSelector::eInDAQ);
           ++telIter) {
        for (fevt::Telescope::ConstPixelIterator pixIter =
               telIter->PixelsBegin(ComponentSelector::eHasData);
             pixIter != telIter->PixelsEnd(ComponentSelector::eHasData);
             ++pixIter) {
          const fevt::PixelRecData& pixRec = pixIter->GetRecData();
          meanVariance += pow(pixRec.GetRMS() / GetCalibConst(*pixIter),2);
          ++nVariance;
        }
      }
      meanVariance /= nVariance;
      return meanVariance;
  }


  FdNoiseAdder::FdNoiseAdder()
  {
  }


  FdNoiseAdder::~FdNoiseAdder()
  {
  }


  VModule::ResultFlag
  FdNoiseAdder::Init()
  {

    CentralConfig* cc = CentralConfig::GetInstance();

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

    const double lgTargetVariance =
      topB.GetChild("lgTargetVariance").Get<double>();
    fTargetVariance = pow(10, lgTargetVariance);

    return eSuccess;
  }


  VModule::ResultFlag
  FdNoiseAdder::Run(evt::Event& event)
  {

    utl::RandomEngine& theRandom =
      RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::ePhysics);

    FEvent& fdEvent = event.GetFEvent();

    for (FEvent::EyeIterator eyeIter =
           fdEvent.EyesBegin(ComponentSelector::eHasData);
         eyeIter!=fdEvent.EyesEnd(ComponentSelector::eHasData);
         ++eyeIter) {

      const double meanVariance = CalcMeanVariance(*eyeIter);
      ostringstream msg;
      msg << " Eye " << eyeIter->GetId() << " <Var> = " << meanVariance;
      INFO(msg.str());
      if (meanVariance > fTargetVariance) {
        INFO(" target < actual variance ");
        continue;
      }
      const double toAdd = fTargetVariance - meanVariance;
      for (fevt::Eye::TelescopeIterator telIter =
             eyeIter->TelescopesBegin(fevt::ComponentSelector::eInDAQ);
           telIter != eyeIter->TelescopesEnd(fevt::ComponentSelector::eInDAQ);
           ++telIter) {
        for (fevt::Telescope::PixelIterator pixIter =
               telIter->PixelsBegin(ComponentSelector::eHasData);
             pixIter != telIter->PixelsEnd(ComponentSelector::eHasData);
             ++pixIter) {
          const double calConst = GetCalibConst(*pixIter);
          fevt::PixelRecData& pixRec = pixIter->GetRecData();
          TraceD& trace = pixRec.GetPhotonTrace();
          for (unsigned int i =  trace.GetStart(); i <= trace.GetStop(); ++i) {
            const double adcOrig = trace[i] / calConst;
            const double adcNoise = RandGauss::shoot(&theRandom.GetEngine(),
                                                     0, sqrt(toAdd));
            trace[i] = (adcOrig + adcNoise) * calConst;
          }

          const double rms =
            TraceAlgorithm::StandardDeviation(trace, 5, 200);
          pixRec.SetRMS(rms);
        }
      }

      const double finalVariance = CalcMeanVariance(*eyeIter);
      msg.str("");
      msg << " Eye " << eyeIter->GetId() << " final <Var> = "
           << finalVariance;
      INFO(msg.str());

    }
    return eSuccess;

  }

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

}