FdBackgroundSimulator.cc

Go to the documentation of this file.

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Trace.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/Vector.h>

#include <fwk/CentralConfig.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>
#include <fdet/Camera.h>
#include <fdet/Mirror.h>
#include <fdet/Corrector.h>
#include <fdet/Filter.h>
#include <fdet/Eye.h>


#include <evt/Event.h>
#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/PixelSimData.h>
#include <fevt/Pixel.h>

#include "FdBackgroundSimulator.h"

using namespace FdBackgroundSimulatorOG;
using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fevt;


VModule::ResultFlag FdBackgroundSimulator::Init(){

  CentralConfig* cc = CentralConfig::GetInstance();

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

  fVerbosity = 0;
  if (topB.GetChild("verbosity"))
    topB.GetChild("verbosity").GetData(fVerbosity);

  string tmp;
  topB.GetChild("bgMode").GetData(tmp);
  if (tmp=="table") {
    fMode = eTable;
  } else if (tmp=="parametric") {
    fMode = eParametric;
  } else if (tmp=="bgLoop") {
    fMode = eBGLoop;
  } else {
    ostringstream err;
    err << "unknown mode for bg model: " << tmp;
    ERROR(err.str());
    return eFailure;
  }

  topB.GetChild("backgroundFlux").GetData(fBackgroundFlux);

  if (fBackgroundFlux.size()!=440) {
    ERROR("Dimensions of the the vector with the photon background is different from 440. ");
    return eFailure;
  }

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "      bg photon model: ";
  switch(fMode) {
      case eTable: info << "XML datacard tables \n"; break;
      case eBGLoop: info << "bgLoop database \n"; break;
      case eParametric: info << "sky-variance parametrisation \n"; break;
      default: info << "UNKNOWN \n"; break;
  }
  INFO(info);

  return eSuccess;

} // end of Init


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

  if (!event.HasFEvent()) {
    ERROR ("Event has no FEvent.");
    return eFailure;
  }

  FEvent&  fEvent = event.GetFEvent();

  for (FEvent::EyeIterator iEye = fEvent.EyesBegin(ComponentSelector::eInDAQ);
       iEye != fEvent.EyesEnd(ComponentSelector::eInDAQ);
       ++iEye){

    if (iEye->GetStatus()==fevt::ComponentSelector::eDeSelected)
      continue;

    fevt::Eye& eyeEvent = *iEye;

    for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(ComponentSelector::eInDAQ);
         iTel != eyeEvent.TelescopesEnd(ComponentSelector::eInDAQ);
         ++iTel) {


      fevt::Telescope& telEvent = *iTel;

      for (fevt::Telescope::PixelIterator iPixel = telEvent.PixelsBegin(ComponentSelector::eInDAQ);
           iPixel != telEvent.PixelsEnd(ComponentSelector::eInDAQ);
           ++iPixel) {

        if (iPixel->HasSimData()) {
          if (fMode == eTable)
            AddBackground(telEvent, event);
          else
            AddBackgroundFromVariance(telEvent, event);

          break; // move to next tel
        }

      }// end loop over Pixels

    }// end loop over Telescopes

  }// end loop over Eyes

  return eSuccess;

}// end of Run



VModule::ResultFlag FdBackgroundSimulator::Finish(){
  return eSuccess;
}// end of Finish


bool FdBackgroundSimulator::AddBackground(fevt::Telescope& telEvent,
                                          evt::Event& /*event*/)
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const double normWavelength = detFD.GetReferenceLambda();

  unsigned int telId = telEvent.GetId();
  unsigned int eyeId = telEvent.GetEyeId();

  if (fVerbosity > 0) {
    ostringstream info;
    info << "Loading photon background to eye=" << eyeId
         << " mirror=" << telId;
    INFO (info);
  }

  const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);
  const fdet::Camera& detCamera = detTel.GetCamera();

  for (unsigned int pixelId = detTel.GetFirstPixelId();
       pixelId <= detTel.GetLastPixelId();
       ++pixelId) {

    if (!telEvent.HasPixel(pixelId))
      telEvent.MakePixel(pixelId);

    fevt::Pixel& pix = telEvent.GetPixel(pixelId);
    const fdet::Pixel& detPixel = detTel.GetPixel(pixelId);

    const double pixelAperture =
      detTel.GetDiaphragmArea() * detPixel.GetSolidAngle();

    if (!pix.HasSimData())
      pix.MakeSimData();

    fevt::PixelSimData& simdata = pix.GetSimData();

    if (simdata.HasPhotonTrace(fevt::FdConstants::eBackground)) {
      ostringstream warn;
      warn << "The pixel: " << pixelId << " already has a background trace.";
      INFO (warn);
      continue;
    }

    const double timeBinSize = detCamera.GetFADCBinSize();

    //use diaPhotonToPE with the correct config of the TelescopeSimulator used
    const string& configSignature = telEvent.GetSimData().GetConfigSignature();
    const double diaPhotonToPE =
      detPixel.GetDiaPhoton2PEFactor(normWavelength, configSignature);
    const double opticalEfficiency =
      diaPhotonToPE / detPixel.GetQEfficiency().Y(normWavelength);

    const double flux = fBackgroundFlux[pixelId-1];
    const double bgPhotons = flux * pixelAperture * timeBinSize * opticalEfficiency;

    simdata.SetMeanBgPhotonFlux(bgPhotons);

    if (fVerbosity > 0) {
      ostringstream info;
      info << "bg-TABLE: Eye n." << eyeId << " mirror n." << telId
           << " pixel n." << pixelId << " bgPhotons " << bgPhotons
           << " OpticalFactor -> " << opticalEfficiency;
      INFO(info);
    }

  } // end loop over Pixels

  return true;

} // end of AddBackground


bool FdBackgroundSimulator::AddBackgroundFromVariance(fevt::Telescope& telEvent, evt::Event& /*event*/)
{
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);
  const fdet::Eye& detEye = detFD.GetEye(detTel.GetParentPhysicalEyeId());
  const fdet::Camera& detCamera = detTel.GetCamera();
  const double timeBinSize = detCamera.GetFADCBinSize();
  const CoordinateSystemPtr& eyeCS = detEye.GetEyeCoordinateSystem();


  // see IEEE Trans.Nucl.Sci.50:1208-1213,2003
  //
  // comment: here we increase the variance by analogElectronicsFactor
  //          to account for later reduction by Gaussian convolution
  //          in FdElectronicsSimulator...
  //
  const double cutoffFrequency   = detCamera.GetCutoffFrequency();
  const double noiseEqBdw        = 0.5*cutoffFrequency*sqrt(kPi/log(2.0));
  const double samplingFrequency = 1./timeBinSize;
  const double analogElectronicsFactor = 2. * noiseEqBdw / samplingFrequency;


  const double normWavelength = detFD.GetReferenceLambda();

  const unsigned int telId = telEvent.GetId();
  const unsigned int eyeId = telEvent.GetEyeId();

  if (fVerbosity > 0) {
    ostringstream info;
    info << "Loading photon background to eye=" << eyeId
         << " mirror=" << telId;
    INFO (info);
  }

  for (unsigned int pixelId = detTel.GetFirstPixelId();
       pixelId <= detTel.GetLastPixelId();
       ++pixelId) {

    if (!telEvent.HasPixel(pixelId))
      telEvent.MakePixel(pixelId);

    fevt::Pixel& pixel = telEvent.GetPixel(pixelId);
    const fdet::Channel& detChannel = detFD.GetChannel(pixel);
    const fdet::Pixel& detPixel = detFD.GetPixel(pixel);

    if (!pixel.HasSimData())
      pixel.MakeSimData();

    fevt::PixelSimData& simdata = pixel.GetSimData();

    if (simdata.HasPhotonTrace(fevt::FdConstants::eBackground)) {
      ostringstream warn;
      warn << "The pixel: " << pixelId << " already has a background trace.";
      INFO (warn);
      continue;
    }

    const double VARelec = detChannel.GetElectronicNoiseVariance(); // [ADC counts**2]

    double VARsky = 0;
    if(fMode==eBGLoop) {
      VARsky = (detChannel.GetADCVariance() - VARelec);
    }
    else {
      const double theta =
        detTel.GetPixel(pixelId).GetDirection().GetTheta(eyeCS);
      const double VARskyStandardASTel = GetParametricSkyVariance(theta);
      VARsky = VARskyStandardASTel * timeBinSize / (100*ns);
    }
    VARsky /= analogElectronicsFactor;

    if (VARsky<=0) {
      simdata.SetMeanBgPhotonFlux(0.0);
      pixel.SetStatus(ComponentSelector::eDeSelected);
      cout << " pixel de-selected: " << pixelId << endl;
      continue;
    }

    const TabulatedFunction& qEff = detPixel.GetQEfficiency();
    const double absGain = detChannel.GetElectronicsGain(); // [ ADC per PE ]

    // see GAP 2004-072 for conversion

    const double Vg = detChannel.GetGainVariance();

    //double gain = 1; // todo
    const double VARskyPE = VARsky / (absGain*absGain);
    const double nPE = VARskyPE / (1+Vg);

    const double bgPhotons = nPE / qEff.Y(normWavelength);

    if (fVerbosity > 0) {
      const double theta = detTel.GetPixel(pixelId).GetDirection().GetTheta(eyeCS);
      ostringstream info;
      if (fMode==eBGLoop) info << " bg-BGLoop: " ;
      else                info << " bg-Prmtrc: ";
      info << " Eye n." << eyeId << " mirror n." << telId
           << " pixel n." << pixelId << " pixel zenith "
           << theta/deg << " var_sky " << VARsky << " bgPhotons " << bgPhotons;
      INFO(info);
    }

    simdata.SetMeanBgPhotonFlux(bgPhotons);

  } // end loop over Pixels

  return eSuccess;

} // end of AddParametricBkg





// variance [ADC^2] of standard Auger South telescopes, see GAP 2006-090
double FdBackgroundSimulator::GetParametricSkyVariance(double theta) const {

  const double f = 0.501;
  const double k = 0.3862;
  const double bl = 14.27;

  double V =  1 / sqrt (1 - 0.96*sin(theta)*sin(theta) );
  double expon = -.4*k*(V-1);
  double var = bl * ( 1 - f + f * V  ) * pow(10., expon);

  return var;
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "cd $AUGER_BASE && make"
// End: