ThresholdCalculator.cc

Go to the documentation of this file.

/*
  Module to calculate the threshold values for FD-Simulation
  \author Joachim Debatin
 */

#include "ThresholdCalculator.h"

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/TabulatedFunction.h>
#include <utl/MathConstants.h>
#include <utl/AugerUnits.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/Md5Sum.h>

#include <fwk/CentralConfig.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 <fevt/ChannelSimData.h>
#include <fevt/Channel.h>

#include <det/Detector.h>

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

#include <TFile.h>

#include <string>

using namespace ThresholdCalculatorKG;
using namespace std;
using namespace utl;
using namespace fwk;
using namespace fevt;
using namespace det;
using namespace fdet;

ThresholdCalculator::ThresholdCalculator() :
  fVerbosityLevel(0),
  fRootFileName("out.root")
{
}

ThresholdCalculator::~ThresholdCalculator()
{
}
  
VModule::ResultFlag
ThresholdCalculator::Init()
{

  CentralConfig* cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("ThresholdCalculator");

  fVerbosityLevel = 0;
  if (topB.GetChild("verbosityLevel"))
    topB.GetChild("verbosityLevel").GetData(fVerbosityLevel);

  //Get Npe from MeanPeFiller
  Branch topB2 = cc->GetTopBranch("MeanNpeFiller");
  if (topB2.GetChild("Npe")) {
    topB2.GetChild("Npe").GetData(fNpe);
  } 

  if (topB2.GetChild("Eye"))
    topB2.GetChild("Eye").GetData(fEye);

  if (topB2.GetChild("Tel"))
    topB2.GetChild("Tel").GetData(fTel);

  Detector::GetInstance().Update(UTCDateTime(2012,1,1,0,0,0).GetTimeStamp());
  const Camera& detCamera = Detector::GetInstance().GetFDetector().GetEye(fEye).GetTelescope(fTel).GetCamera();
  fBoxcarSumLength = detCamera.GetFLTBoxcarSumLength();

  //Get Baseline form ElectronicsSimulator
  Branch topB3 = cc->GetTopBranch("FdElectronicsSimulator");
  topB3.GetChild("baseline").GetData(fBaseline);

  if (topB.GetChild("RootFileName"))
    topB.GetChild("RootFileName").GetData(fRootFileName);
  fHistoFile = new TFile(fRootFileName.c_str(), "RECREATE");

  //Calculate histogram boundaries
  const double gain = detCamera.GetElectronicsGain();
  const double electronicVariance = detCamera.GetElectronicNoiseVariance();
  
  for (unsigned int i=0; i<fNpe.size(); ++i){
    ostringstream sampleName;
    sampleName << "BinSum" << i+1;
    
    int NSamplesMin = int(fBoxcarSumLength * fBaseline);
    
    //Factor of 3 is manually optimised. If yout get warnigs enlarge.
    int NSamplesMax = int(fBoxcarSumLength * (fBaseline + electronicVariance + 3 * sqrt(fNpe[i] * gain)) + .5);
    
    fSamples.push_back(new TH1D(sampleName.str().c_str(), sampleName.str().c_str(), NSamplesMax - NSamplesMin + 1, NSamplesMin-.5, NSamplesMax+.5));

    //Histogram maximum bin
    fSampleMaximum.push_back(NSamplesMax);

    fNumberOfTraces.push_back(0);
  }
  
  // info output
  ostringstream info;
  info << "Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n";

  if (fVerbosityLevel)
    info << "Verbosity level: " << fVerbosityLevel << '\n';

  if (fBoxcarSumLength)
    info << "Used NSamples value: " << fBoxcarSumLength << "\n";

  INFO(info);

  return eSuccess;
}


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

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

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

    fevt::Eye& eyeEvent = *iEye;

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

      //GetFDetector
      const fdet::FDetector& detFD = Detector::GetInstance().GetFDetector();
      fevt::Telescope& telEvent = *iTel;
      const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);

      //Needed to correct for scaling in ElectronicsSimulator
      const double normWavelength = detFD.GetReferenceLambda();

      //Loop pixels/channels
      for (unsigned int channelId=1; channelId<=detTel.GetLastPixelId(); ++channelId) {

        const fdet::Channel& detChannel = detTel.GetChannel(channelId);
        unsigned int pixelId = detChannel.GetPixelId();      

        if (!telEvent.HasChannel(channelId))
          continue;
        fevt::Channel& channelEvent = telEvent.GetChannel(channelId);

        if (!channelEvent.HasSimData())
          continue;
        fevt::ChannelSimData& channel_sim = channelEvent.GetSimData();

        if (!channel_sim.HasFADCTrace (fevt::FdConstants::eTotal))
          continue;

        if (!telEvent.HasPixel(pixelId))
          continue;
        fevt::Pixel& pixelEvent = telEvent.GetPixel(pixelId);

        if (pixelEvent.GetStatus() == ComponentSelector::eDeSelected) {
          continue;
        }

        if (!pixelEvent.HasSimData())
          continue;
        fevt::PixelSimData& simPixel = pixelEvent.GetSimData();

        const fdet::Pixel& detPixel = detFD.GetPixel(pixelEvent);
        const TabulatedFunction& qEff = detPixel.GetQEfficiency();
        const double qEffAtNormWavelength = qEff.Y(normWavelength);

        //Get meanPe cleaned of scaling in MeanPeFiller
        const double meanPe = simPixel.GetMeanBgPhotonFlux() * qEffAtNormWavelength;
        
        int meanPeBin = 0;
        for(unsigned int i = 0; i < fNpe.size(); ++i){
          if(abs(meanPe-fNpe[i])<.1){
            meanPeBin = i;
          }
        }

        //Calculate boxcarsums
        const TraceI& fadc_trace = channel_sim.GetFADCTrace(fevt::FdConstants::eTotal);
        vector<int> boxcar(fBoxcarSumLength, 0);
        fNumberOfTraces[meanPeBin]+=1;

        for (unsigned int bin=0; bin<fadc_trace.GetSize(); ++bin){
          // build boxcar
          for (unsigned int i=0; i<fBoxcarSumLength-1; ++i)
            boxcar[i] = boxcar[i+1];
          boxcar[fBoxcarSumLength-1] = fadc_trace [bin];
          int sum = 0;
          for (int i=fBoxcarSumLength-1; i >= 0; --i){
            if(boxcar[i]==0){
              sum = 0;
              break;
            }else{
              sum+=boxcar[i];
            }
          }
          if(sum){
            fSamples[meanPeBin]->Fill(sum);
            if(fVerbosityLevel>=1){
              if(sum > fSampleMaximum[meanPeBin]){
                ostringstream msg;
                msg << "BoxCarSum bigger than maximum of histogram. Calculation is wrong!! " << sum << " " << fSampleMaximum[meanPeBin];
                WARNING(msg);
              }
            }
          }
        }
      }
    }
  }
  return eSuccess;
}


VModule::ResultFlag
ThresholdCalculator::Finish()
{

  const fdet::Camera & detCamera = Detector::GetInstance().GetFDetector().GetEye(fEye).GetTelescope(fTel).GetCamera();

  //Cumulative histograms
  for (unsigned int i=0; i<fSamples.size(); ++i){
    fSummedSamples.push_back(fSamples[i]->GetCumulative(false));
  }

  //convert count to frequency
  for (unsigned int i=0; i<fSummedSamples.size(); ++i){
    fSummedSamples[i]->Scale(1. / (fNumberOfTraces[i] * detCamera.GetFADCTraceLength() * detCamera.GetFADCBinSize()));
  }

  //Create threshold signature
  const double FLTTriggerRate = detCamera.GetFLTTriggerRate();
  ostringstream ThresholdSignature;
  ThresholdSignature << "Threshold:" 
                     << " Gain = " << detCamera.GetElectronicsGain()
                     << " GainVariance = " << detCamera.GetGainVariance()
                     << " ElectronicNoise = " << detCamera.GetElectronicNoiseVariance()
                     << " CutoffFrequency = " << detCamera.GetCutoffFrequency()
                     << " FLTTriggerRate = " << detCamera.GetFLTTriggerRate()
                     << " FLTBoxcarSize = " << fBoxcarSumLength;

  if (fEye>0 && fEye<5){
    if (fTel>0 && fTel<7){
      ThresholdSignature << " Telescopetype = standard";
    }
  }else if (fEye==5){
    if (fTel>0 && fTel<4){
      ThresholdSignature << " Telescopetype = HEAT";
    }
  }
 
  //Start OutputFile
  ofstream thresholdXml("FSimulationThreshold.generated.xml.in");
  thresholdXml << "<?xml version=\"1.0\" encoding=\"iso-8859-1\"?>\n\n"
               << "<FSimulationThreshold xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n"
               << "                      xsi:noNamespaceSchemaLocation='@XMLSCHEMALOCATION@/FSimulationThreshold.xsd'>\n\n"
               << "  <Threshold signature=\"" << Md5Sum(ThresholdSignature.str()).GetHexDigest() << "\"> \n" 
               << "  <!--" << ThresholdSignature.str() <<"-->\n";

  //Calculate Threshold
  double xValues[fSummedSamples.size()];
  int yValues[fSummedSamples.size()];
  for (unsigned int i=0; i<fSummedSamples.size(); ++i){
    int x = 0;
    int& binnumber = x;
    const double dummy = fSummedSamples[i]->GetBinWithContent(FLTTriggerRate, binnumber, 0, fSummedSamples[i]->GetNbinsX(), 100);
    const int threshold = int(fSummedSamples[i]->GetBinCenter(binnumber) - fBoxcarSumLength * fBaseline - .5);

    if (fVerbosityLevel>=1){
      ostringstream info;
      info << "Result of Threshold Calculation for"
           << " meanPe = " << fNpe[i]
           << " Threshold = " << threshold
           << " Uncertainty = " << dummy
           << " Binnumber = " << x;
      INFO(info);
    }

    xValues[i] = fNpe[i];
    yValues[i] = threshold;
  }

  //Write threshold to xml
  thresholdXml << "    <x>";
  for(unsigned int i = 0; i < fSummedSamples.size(); ++i)
    thresholdXml << " " << xValues[i];
  thresholdXml << "</x>\n";

  thresholdXml << "    <y>";
  for(unsigned int i = 0; i < fSummedSamples.size(); ++i)
    thresholdXml << " " << yValues[i];
  thresholdXml << "</y>\n";

  thresholdXml << "  </Threshold>\n\n"
               << "</FSimulationThreshold>";

  fHistoFile->cd();
  fHistoFile->Write();
  fHistoFile->Close();

  return eSuccess;
}