RdChannelResponseIncorporator.cc

Go to the documentation of this file.

#include "RdChannelResponseIncorporator.h"

#include <fwk/CentralConfig.h>
#include <fwk/RandomEngineRegistry.h>

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

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

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

#include <CLHEP/Random/RandGauss.h>


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

using CLHEP::RandGauss;


namespace RdChannelResponseIncorporator {


  VModule::ResultFlag
  RdChannelResponseIncorporator::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelResponseIncorporator");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

    topBranch.GetChild("NumResponsesToCache").GetData(fNumResponsesToCache);
    topBranch.GetChild("ForwardResponseOnFirstCall").GetData(fForwardResponseOnFirstCall);
    topBranch.GetChild("OverrideInverseBandLimits").GetData(fInverseBandLimitsOverride);
    topBranch.GetChild("InverseLowerFrequencyLimit").GetData(fInverseLowerFrequency);
    topBranch.GetChild("InverseUpperFrequencyLimit").GetData(fInverseUpperFrequency);
    topBranch.GetChild("OverrideForwardBandLimits").GetData(fForwardBandLimitsOverride);
    topBranch.GetChild("ForwardLowerFrequencyLimit").GetData(fForwardLowerFrequency);
    topBranch.GetChild("ForwardUpperFrequencyLimit").GetData(fForwardUpperFrequency);

    topBranch.GetChild("SmearChannelTraceToMimicInaccurateHardwareDescription").GetData(
      fSmearChannelTraceToMimicInaccurateHardwareDescription);

    topBranch.GetChild("RelativeGaussianUcertainty").GetData(fRelativeGaussianUcertainty);
    if (fSmearChannelTraceToMimicInaccurateHardwareDescription)
      fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

    ostringstream info;
    info << "Read config:"
         << "\n\t ForwardResponseOnFirstCall                            : "
         << fForwardResponseOnFirstCall
         << "\n\t SmearChannelTraceToMimicInaccurateHardwareDescription : "
         << fSmearChannelTraceToMimicInaccurateHardwareDescription
         << "\n\t RelativeGaussianUcertainty                            : "
         << fRelativeGaussianUcertainty << "\n";
    INFO(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelResponseIncorporator::Run(evt::Event& event)
  {
    // Check if there are events at all
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    det::Detector& thedet = det::Detector::GetInstance();
    REvent& rEvent = event.GetREvent();

    ostringstream message;
    message << "Run first time: " << (rEvent.GetChannelResponseApplied() ? "no" : "yes");
    INFO(message);

    // loop through stations and for every station through every channel
    for (auto& station : rEvent.StationsRange()) {
      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive()) {
          continue;
        }

        const rdet::Channel& channeldet =
            thedet.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId());
        ResponseMap channelresponsemap = channeldet.GetResponseMap();

        const TabulatedFunctionComplexLgAmpPhase& channelresponse = GetCachedOverallResponse(channelresponsemap);

        // for later check if in frequency range
        const double designLowerfrequency = channeldet.GetDesignLowerFreq();
        const double designUpperfrequency = channeldet.GetDesignUpperFreq();


        if (!fForwardResponseOnFirstCall) {
          // module is just called once in module sequence, i.e. this is a data pipeline
          // -> apply inverse response
          message.str("");
          message << "Applying inverse detector response to channel "
                  << channel.GetId() << " of station " << station.GetId();
          INFOIntermediate(message);

          if (fInverseBandLimitsOverride)
            ApplyResponse(channel, channelresponse, true, fInverseLowerFrequency, fInverseUpperFrequency);
          else
            ApplyResponse(channel, channelresponse, true, designLowerfrequency, designUpperfrequency);
        } else {
          // module is called more than once in module sequence; now check whether this is first or later call
          if (rEvent.GetChannelResponseApplied()) {
            // this is (at least) the second call of the module for this event
            message.str("");
            message << "Later call: applying inverse response to channel "
                    << channel.GetId() << " of station " << station.GetId();
            INFODebug(message);

            if (fInverseBandLimitsOverride)
              ApplyResponse(channel, channelresponse, true, fInverseLowerFrequency, fInverseUpperFrequency);
            else
              ApplyResponse(channel, channelresponse, true, designLowerfrequency, designUpperfrequency);
          } else {
            // this is the first call of the module for this event
            message.str("");
            message << "First call: applying forward response to channel "
                    << channel.GetId() << " of station " << station.GetId();
            INFODebug(message);

            if (fForwardBandLimitsOverride)
              ApplyResponse(channel, channelresponse, false, fForwardLowerFrequency, fForwardUpperFrequency);
            else
              ApplyResponse(channel, channelresponse, false, designLowerfrequency, designUpperfrequency);

            if (fSmearChannelTraceToMimicInaccurateHardwareDescription) {
              revt::ChannelTimeSeries& trace = channel.GetChannelTimeSeries();
              const double smearFactor = RandGauss::shoot(&fRandomEngine->GetEngine(), 1,
                                                          fRelativeGaussianUcertainty);

              message.str("");
              message << "Station " << station.GetId() << " channel " << channel.GetId()
                      << " smeared with " << smearFactor;
              INFOIntermediate(message);

              for (auto& val : trace)
                val *= smearFactor;
            }
          }
        }
      }
    }

    // note down application of channel response in event
    rEvent.SetChannelResponseApplied(true);

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelResponseIncorporator::Finish()
  {
    ostringstream msg;
    msg << "RdChannelResponseIncorporator cache hit fraction was "
        << static_cast<double>(fCacheHits) / static_cast<double>(fCacheAsks)
        << " for a cache with " << fNumResponsesToCache << " entries.";
    INFO(msg);

    return eSuccess;
  }


  const
  TabulatedFunctionComplexLgAmpPhase&
  RdChannelResponseIncorporator::GetCachedOverallResponse(const ResponseMap& responsemap) const
  {
    ++fCacheAsks;
    vector<CacheRecord>::iterator leastusedentry = fCacheRegistry.begin();
    for (vector<CacheRecord>::iterator i = fCacheRegistry.begin(); i != fCacheRegistry.end(); ++i) {
      if (i->get<0>() == responsemap) { // is the response map in the CacheRegistry equal to the specified one?
        ++fCacheHits;      // then register cache hit
        ++i->get<2>();     // count hit to cache entry
        return i->get<1>(); // return the found cache entry
      }
      // keep track of the least-used cache entry
      if (i->get<2>() < leastusedentry->get<2>())
        leastusedentry = i;
    }

    // if we arrive here, then the cache did not contain the appropriate entry
    if (fCacheRegistry.size() == fNumResponsesToCache)  // if cache registry is full
      fCacheRegistry.erase(leastusedentry);             // make room for the new entry
    fCacheRegistry.emplace_back(responsemap, CalculateOverallResponse(responsemap), 0);
    return fCacheRegistry.back().get<1>();              // return the freshly generated entry
  }


  TabulatedFunctionComplexLgAmpPhase
  RdChannelResponseIncorporator::CalculateOverallResponse(const ResponseMap& responsemap) const
  {
    TabulatedFunctionComplexLgAmpPhase totalresponse;
    const rdet::RDetector& rdetector = det::Detector::GetInstance().GetRDetector();

    // treat first entry of response map manually
    ResponseMap::const_iterator beginentry = responsemap.Begin();
    // retrieve the response profile of the first entry
    totalresponse = rdetector.GetHardwareResponseProfile(beginentry->first);
    totalresponse.Pow(beginentry->second);  // scale it with the correct weight
    ++beginentry;

    // now accumulate the rest of the response map
    for (auto it = beginentry; it != responsemap.End(); ++it) {
      TabulatedFunctionComplexLgAmpPhase tempresponse =
        rdetector.GetHardwareResponseProfile(it->first);
      tempresponse.Pow(it->second);
      totalresponse *= tempresponse;
    }

    return totalresponse;
  }


  void
  RdChannelResponseIncorporator::ApplyResponse(Channel& channel, const TabulatedFunctionComplexLgAmpPhase& response,
                                               bool inverse, double lowerfreq, double upperfreq) const
  {
    ChannelFrequencySpectrum& channelspectrum = channel.GetChannelFrequencySpectrum();
    for (ChannelFrequencySpectrum::SizeType i = 0; i < channelspectrum.GetSize(); ++i) {
      const double freq = channel.GetFrequencyOfBin(i);
      if ((freq < lowerfreq) || (freq > upperfreq)) {
        channelspectrum[i] = complex<double>(0.0, 0.0); // bin is outside design specifications, set it to zero
      } else {
        if (inverse)
          channelspectrum[i] /= response.Y(freq).GetComplex();
        else
          channelspectrum[i] *= response.Y(freq).GetComplex();
      }
    }
  }

}