RdChannelGalacticCalibrator.cc

Go to the documentation of this file.

#include "RdChannelGalacticCalibrator.h"

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/StationGPSData.h>

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

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

#include <fstream>


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


namespace RdChannelGalacticCalibrator {

  VModule::ResultFlag
  RdChannelGalacticCalibrator::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelGalacticCalibrator");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("NameConstantsFile").GetData(fNameConstantsFile);
    topBranch.GetChild("RejectIfCalibrationFailed").GetData(fRejectIfCalibrationFailed);

    GetCalibrationConstants();
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelGalacticCalibrator::Run(evt::Event& event)
  {
    REvent& rEvent = event.GetREvent();
    ostringstream info;

    for (auto& station : rEvent.CandidateStationsRange()) {

      const int stationId = station.GetId() - 100;
      const unsigned long gpsSec = station.GetGPSData().GetSecond();
      int month = 0;
      int year = 0;
      GetDataEvent(gpsSec, year, month);

      info.str("");
      info << "Run galactic calibration for station " << stationId
           << " in " << month << "/" << year << ".";
      INFODebug(info);

      bool stationRejected = false;

      for (auto& channel : station.ChannelsRange()) {
        if (!channel.IsActive() || stationRejected)
          continue;

        const int channelId = channel.GetId();
        const auto& rDetector = det::Detector::GetInstance().GetRDetector();
        const auto& detChannel = rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());

        auto& spectrum = channel.GetChannelFrequencySpectrum();
        const double lowerDesignFreq = detChannel.GetDesignLowerFreq();
        const double upperDesignFreq = detChannel.GetDesignUpperFreq();
        const int lowerBin = int(lowerDesignFreq / spectrum.GetBinning());
        const int upperBin = int(upperDesignFreq / spectrum.GetBinning());
        const double lowerFreq = channel.GetFrequencyOfBin(lowerBin + 1) / MHz;
        const double upperFreq = channel.GetFrequencyOfBin(upperBin + 1) / MHz;

        bool frequencyFilled[100] = { false };

        const int indexMonth = GetMonthIndex(month, year);

        // save corrections factors in channel parameter storage
        if (!channel.HasRecData()) {
          channel.MakeRecData();
        }

        for (ChannelFrequencySpectrum::SizeType f = 0; f < spectrum.GetSize(); ++f) {
          const double currentFreq = channel.GetFrequencyOfBin(f) / MHz;
          if (currentFreq > upperFreq || currentFreq < lowerFreq) {
            continue;
          }

          calibrationMap channelMap;
          if (channelId == 1) {
            channelMap = fCalibrationConstantChannel1[stationId];
          }
          if (channelId == 2) {
            channelMap = fCalibrationConstantChannel2[stationId];
          }

          const int intFreq = int(currentFreq);
          const auto it = channelMap.find(make_pair(indexMonth, intFreq));
          if (it == channelMap.end()) {
            info.str("");
            info << "no calibration constant in map for station " << stationId
                 << ", channel " << channel.GetId() << " at a frequency of "
                 << intFreq << "MHz.";

            if (fRejectIfCalibrationFailed) {
              info << " Rejecting station " << stationId << ".";
              INFOFinal(info);

              // does not fit perfectly but best we have
              station.SetRejectedReason(eLNANotCalibrated);

              // break the current loop over frequencies and do not iterate other channels of same station
              stationRejected = true;
              break;
            }

            INFOIntermediate(info);

            // do no correct current frequency but continue correcting other frequencies
            continue;
          }

          const double C0 = it->second;
          spectrum[f] /= C0;

          if (!frequencyFilled[intFreq]) {
            const auto par = GetParameterName(intFreq);
            channel.GetRecData().SetParameter(par, C0);
            frequencyFilled[intFreq] = true;
          }
        }
      }
    }

    return eSuccess;
  }


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


  // function to calculate LST
  void
  RdChannelGalacticCalibrator::GetDataEvent(const int gps,int& year, int& month)
  {
    const time_t utc = gps + 315964819 + 16; // GPS UTC offset, Leap Seconds have to be updated
    const tm* const t = gmtime(&utc);
    year = t->tm_year + 1900;
    month = t->tm_mon + 1;
  }


  void
  RdChannelGalacticCalibrator::GetCalibrationConstants()
  {
    ifstream obj;
    obj.open(fNameConstantsFile.c_str());

    if (!obj.is_open()) {
      const string error = "File with calibration constants not found. The trace cannot be corrected.";
      ERROR(error);
      throw utl::AugerException(error);
    }

    ostringstream info;
    info << "file found : " << fNameConstantsFile;
    INFO(info);

    while (!obj.eof()) {
      int channel, year, month, antenna, frequency;
      double calConstant;
      obj >> channel >> year >> month >> antenna >> frequency >> calConstant;

      const int indexMonth = GetMonthIndex(month, year);
      if (channel == 1) {
        fCalibrationConstantChannel1[antenna].insert({{indexMonth, frequency}, calConstant});
      }
      if (channel == 2){
        fCalibrationConstantChannel2[antenna].insert({{indexMonth, frequency}, calConstant});
      }
    }
  }


  unsigned int
  RdChannelGalacticCalibrator::GetMonthIndex(int month, int year)
  {
    return 12 * (year - fInitYear) + month;
  }


  revt::ChannelRRecDataQuantities
  RdChannelGalacticCalibrator::GetParameterName(unsigned int intFreq)
  {
    if (intFreq == 30)
      return eChannelGalacticCalibrationFactor30MHz;
    if (intFreq == 31)
      return eChannelGalacticCalibrationFactor31MHz;
    if (intFreq == 32)
      return eChannelGalacticCalibrationFactor32MHz;
    if (intFreq == 33)
      return eChannelGalacticCalibrationFactor33MHz;
    if (intFreq == 34)
      return eChannelGalacticCalibrationFactor34MHz;
    if (intFreq == 35)
      return eChannelGalacticCalibrationFactor35MHz;
    if (intFreq == 36)
      return eChannelGalacticCalibrationFactor36MHz;
    if (intFreq == 37)
      return eChannelGalacticCalibrationFactor37MHz;
    if (intFreq == 38)
      return eChannelGalacticCalibrationFactor38MHz;
    if (intFreq == 39)
      return eChannelGalacticCalibrationFactor39MHz;
    if (intFreq == 40)
      return eChannelGalacticCalibrationFactor40MHz;
    if (intFreq == 41)
      return eChannelGalacticCalibrationFactor41MHz;
    if (intFreq == 42)
      return eChannelGalacticCalibrationFactor42MHz;
    if (intFreq == 43)
      return eChannelGalacticCalibrationFactor43MHz;
    if (intFreq == 44)
      return eChannelGalacticCalibrationFactor44MHz;
    if (intFreq == 45)
      return eChannelGalacticCalibrationFactor45MHz;
    if (intFreq == 46)
      return eChannelGalacticCalibrationFactor46MHz;
    if (intFreq == 47)
      return eChannelGalacticCalibrationFactor47MHz;
    if (intFreq == 48)
      return eChannelGalacticCalibrationFactor48MHz;
    if (intFreq == 49)
      return eChannelGalacticCalibrationFactor49MHz;
    if (intFreq == 50)
      return eChannelGalacticCalibrationFactor50MHz;
    if (intFreq == 51)
      return eChannelGalacticCalibrationFactor51MHz;
    if (intFreq == 52)
      return eChannelGalacticCalibrationFactor52MHz;
    if (intFreq == 53)
      return eChannelGalacticCalibrationFactor53MHz;
    if (intFreq == 54)
      return eChannelGalacticCalibrationFactor54MHz;
    if (intFreq == 55)
      return eChannelGalacticCalibrationFactor55MHz;
    if (intFreq == 56)
      return eChannelGalacticCalibrationFactor56MHz;
    if (intFreq == 57)
      return eChannelGalacticCalibrationFactor57MHz;
    if (intFreq == 58)
      return eChannelGalacticCalibrationFactor58MHz;
    if (intFreq == 59)
      return eChannelGalacticCalibrationFactor59MHz;
    if (intFreq == 60)
      return eChannelGalacticCalibrationFactor60MHz;
    if (intFreq == 61)
      return eChannelGalacticCalibrationFactor61MHz;
    if (intFreq == 62)
      return eChannelGalacticCalibrationFactor62MHz;
    if (intFreq == 63)
      return eChannelGalacticCalibrationFactor63MHz;
    if (intFreq == 64)
      return eChannelGalacticCalibrationFactor64MHz;
    if (intFreq == 65)
      return eChannelGalacticCalibrationFactor65MHz;
    if (intFreq == 66)
      return eChannelGalacticCalibrationFactor66MHz;
    if (intFreq == 67)
      return eChannelGalacticCalibrationFactor67MHz;
    if (intFreq == 68)
      return eChannelGalacticCalibrationFactor68MHz;
    if (intFreq == 69)
      return eChannelGalacticCalibrationFactor69MHz;
    if (intFreq == 70)
      return eChannelGalacticCalibrationFactor70MHz;
    if (intFreq == 71)
      return eChannelGalacticCalibrationFactor71MHz;
    if (intFreq == 72)
      return eChannelGalacticCalibrationFactor72MHz;
    if (intFreq == 73)
      return eChannelGalacticCalibrationFactor73MHz;
    if (intFreq == 74)
      return eChannelGalacticCalibrationFactor74MHz;
    if (intFreq == 75)
      return eChannelGalacticCalibrationFactor75MHz;
    if (intFreq == 76)
      return eChannelGalacticCalibrationFactor76MHz;
    if (intFreq == 77)
      return eChannelGalacticCalibrationFactor77MHz;
    if (intFreq == 78)
      return eChannelGalacticCalibrationFactor78MHz;
    if (intFreq == 79)
      return eChannelGalacticCalibrationFactor79MHz;
    if (intFreq == 80)
      return eChannelGalacticCalibrationFactor80MHz;

    throw logic_error("no corresponding parameter found");
  }

}