RdChannelGalacticBackgroundCalibrator.cc

Go to the documentation of this file.

#include "RdChannelGalacticBackgroundCalibrator.h"
#include <utl/ErrorLogger.h>

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/StationRecData.h>
#include <revt/StationTriggerData.h>
#include <revt/StationGPSData.h>

#include <rdet/RDetector.h>
#include <rdet/Channel.h>

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>



using namespace std;
using namespace utl;
using namespace fwk;
using namespace revt;
using namespace rdet;


namespace RdChannelGalacticBackgroundCalibrator {

  VModule::ResultFlag
  RdChannelGalacticBackgroundCalibrator::Init()
  {
    INFO("UserModule::Init()");

    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelGalacticBackgroundCalibrator");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("InputFileName").GetData(fInputFileName);
    topBranch.GetChild("InterpolateFilterGaps").GetData(fInterpolateFilterGaps);
    topBranch.GetChild("RejectUncalibratedStations").GetData(fRejectUncalibratedStations);

    // open input file and fill calibration constants in vectors
    fstream finput;
    finput.open(fInputFileName.c_str(), ios::in);
    if(!finput.good()) {
      ERROR("Input file could not be opened");
      return eFailure;
    }
    while(finput.good()) {
      std::vector<double> vTmp;
      int inputStationId, inputChannel;
      double inputFrequency, inputConstant;
      finput >> inputStationId >> inputChannel >> inputFrequency >> inputConstant;
      v_channelfrequencies[inputStationId][inputChannel-1].push_back(inputFrequency);
      v_channelconstants[inputStationId][inputChannel-1].push_back(inputConstant);
    }
    finput.close();

    // linearly interpolate filter gaps if option is selected
    if (fInterpolateFilterGaps) {
      for (int jj=0; jj<200; jj++) {
        for (int ichannel=0; ichannel<2; ichannel++) {
          for (unsigned int ibin=1; ibin<=v_channelconstants[jj][ichannel].size(); ibin++) {
            if (v_channelconstants[jj][ichannel][ibin] == 0.0 &&
                v_channelconstants[jj][ichannel][ibin-1] != 0.0) {
              int j=0;
              while (v_channelconstants[jj][ichannel][j+ibin] == 0 &&
                     (j+ibin)<v_channelconstants[jj][ichannel].size()-1)
                j++;
              if ((j+ibin)!=v_channelconstants[jj][ichannel].size()-1) {
                double loweredge = v_channelconstants[jj][ichannel][ibin-1];
                double upperedge = v_channelconstants[jj][ichannel][j+ibin];
                for (int ifill=1; ifill<=j;ifill++) {
                  v_channelconstants[jj][ichannel][ibin+ifill-1] =
                    loweredge + (upperedge-loweredge)/(j+1) * ifill;
                }
              }
            }
          }
        }
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelGalacticBackgroundCalibrator::Run(evt::Event& event)
  {
    INFO("RdChannelGalacticBackgroundCalibrator::Run()");

    REvent& rEvent = event.GetREvent();
    std::ostringstream info;

    for (auto& station : rEvent.CandidateStationsRange()) {
      const int fStationId = station.GetId() - 100;

      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        const int fChannelId = channel.GetId();

        if(v_channelfrequencies[fStationId][fChannelId-1].size() == 0) {
          info << "Station ID " << fStationId+100 << " with channel "
               << fChannelId << " not in calibration list!";

          if (fRejectUncalibratedStations) {
            info << " Station rejected!";
            // probably a new rejection flag should be created for this
            station.SetRejectedReason(revt::eManuallyRejected);
          }
          WARNING(info);
          continue;
        }

        const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
        const rdet::Channel& detChannel = rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());
        const double lowerDesignFreq = detChannel.GetDesignLowerFreq();
        const double upperDesignFreq = detChannel.GetDesignUpperFreq();

        revt::ChannelFrequencySpectrum& spectrum = channel.GetChannelFrequencySpectrum();

        const int lowerBin = (int)(lowerDesignFreq / spectrum.GetBinning());
        const int upperBin = (int)(upperDesignFreq / spectrum.GetBinning());
        const unsigned int nSpecBins = upperBin - lowerBin;

        const double lowerFreq = channel.GetFrequencyOfBin(lowerBin+1) / MHz;
        const double upperFreq = channel.GetFrequencyOfBin(upperBin+1) / MHz;

        const double spectrumBinning = spectrum.GetBinning();
        const unsigned int spectrumSize = spectrum.GetSize();

        double currentFreq = 0;
        double nearestFreqDist = 0;
        int jjNearestFreq = 0;

        if (nSpecBins != v_channelfrequencies[fStationId][fChannelId-1].size()) {
          info << "The calibration constants (" << v_channelfrequencies[fStationId][fChannelId-1].size()
               << ") of station ID " << fStationId+100 << " with channel " << fChannelId
               << " do not match the number of spectral bins (" << nSpecBins << ").";
          WARNING(info);
        }

        // loop to DIVIDE spectral elements with calibration factors
        // the nearest frequency in the calibration constant list is applied
        for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrumSize; i++) {
          currentFreq=channel.GetFrequencyOfBin(i) / MHz;
          if (currentFreq >= lowerFreq && currentFreq < upperFreq) {
            nearestFreqDist = upperFreq;
            jjNearestFreq = 0;
            for (unsigned int jj=0; jj<v_channelfrequencies[fStationId][fChannelId-1].size(); jj++) {
              if (abs(v_channelfrequencies[fStationId][fChannelId-1][jj] -
                  (currentFreq + 0.5*spectrumBinning/MHz)) < nearestFreqDist) {
                nearestFreqDist = abs(v_channelfrequencies[fStationId][fChannelId-1][jj] -
                    (currentFreq + 0.5*spectrumBinning/MHz));
                jjNearestFreq = jj;
              }
            }
            info << fixed << "Station ID " << fStationId + 100 << " Channel "
                 << fChannelId << " has a calibration constant of "
                 << v_channelconstants[fStationId][fChannelId-1][jjNearestFreq]
                 << " applied to frequency " << currentFreq + 0.5*spectrumBinning/MHz
                 << " MHz." << endl;
            INFODebug(info);

            if (v_channelconstants[fStationId][fChannelId-1][jjNearestFreq] == 0.0)
              // set frequency bin to zero if calibration constant is zero (gap)
              spectrum[i] = 0.0;
            else
              spectrum[i] /= v_channelconstants[fStationId][fChannelId-1][jjNearestFreq];
          }
        }

      } // channel loop
    } // station loop
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelGalacticBackgroundCalibrator::Finish()
  {
    INFO("RdChannelGalacticBackgroundCalibrator::Finish()");

    return eSuccess;
  }
}