RdChannelASCIINoiseImporterRD.cc

Go to the documentation of this file.

#include "RdChannelASCIINoiseImporterRD.h"

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Station.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>
#include <utl/RandomEngine.h>
#include <utl/UTCDateTime.h>
#include <utl/TimeStamp.h>

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

#include <boost/algorithm/string.hpp>

#include <vector>
#include <string>
#include <iostream>
#include <fstream>

#include <CLHEP/Random/RandFlat.h>


namespace RdChannelASCIINoiseImporterRD {

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

    topBranch.GetChild("NoiseFilePathCh0").GetData(fNoiseFilePathCh0);
    topBranch.GetChild("NoiseFilePathCh1").GetData(fNoiseFilePathCh1);
    topBranch.GetChild("AddNoiseFromSameTimeToBothChannels").GetData(fAddNoiseFromSameTimeToBothChannels);
    topBranch.GetChild("CheckTimeWindow").GetData(fCheckTimeWindow);
    topBranch.GetChild("TimeWindow").GetData(fMaxTimeDifference);
    topBranch.GetChild("AllowTraceShifting").GetData(fTraceShifting);
    topBranch.GetChild("SkipTrailingColumnsInNoiseFile").GetData(fSkipTrailingColumnsInNoiseFile);
    topBranch.GetChild("GainFactor").GetData(fGainFactor);

    topBranch.GetChild("ReplaceDataWithNoise").GetData(fReplaceDataWithNoise);

    fRandomEngine =
      &fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector).GetEngine();

    ReadNoiseFile(fNoiseFilePathCh0, 0);
    ReadNoiseFile(fNoiseFilePathCh1, 1);

    if (fAddNoiseFromSameTimeToBothChannels) {
      if (fNoiseTraces[0].size() != fNoiseTraces[1].size() || fNoiseTraces[0].size() == 0) {
        ERROR("Imported noise for the two channels are not same size or size is 0!");
        return eFailure;
      }
    }

    std::ostringstream infostr;
    infostr << "\n\tRead " << fNoiseTraces[0].size() << " traces for channel 0 and "
            << fNoiseTraces[1].size() << " traces for channel 1"
            << "\n\tAddNoiseFromSameTimeToBothChannels : " << fAddNoiseFromSameTimeToBothChannels
            << "\n\tAllow trace shifting               : " << fTraceShifting
            << "\n\tCheck time window (not implmented) : " << fCheckTimeWindow;

    infostr << std::endl;
    INFO(infostr);
    return eSuccess;
  }


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

    det::Detector& det = det::Detector::GetInstance();
    revt::REvent& rEvent = event.GetREvent();
    const rdet::RDetector& rDetector = det.GetRDetector();

    std::ostringstream infostr;
    std::ostringstream warnstr;
    std::ostringstream debugstr;

    // Adding noise to the data, depending on if the station exists in both: noise-file and data.
    for (auto& station : rEvent.StationsRange()) {
      const int stationID = station.GetId();

      // choose same trace (index) for both channel
      unsigned int traceIndex = CLHEP::RandFlat::shootInt(fRandomEngine, 0, fNoiseTraces[0].size());
      unsigned int sampleIndexNoiseStart = 0;

      // Traces have a lenght of 2048 (samples) / 250 MHz = 8192 ns
      // Shift traces by multiple of 100ns, i.e., 25 sample
      // This does not prevent the trace end to be within the signal window
      if (fTraceShifting && fAddNoiseFromSameTimeToBothChannels) {
        sampleIndexNoiseStart = CLHEP::RandFlat::shootInt(fRandomEngine, 0, 25) *
                                fNoiseTraces[0][traceIndex].size() / 25;
      }

      const rdet::Station& rDetStation = rDetector.GetStation(stationID);
      for (auto& channel : station.ChannelsRange()) {
        revt::ChannelADCTimeSeries& timeSeries = channel.GetChannelADCTimeSeries();
        const unsigned int channelID = channel.GetId();

        // Checking if the channel of the processed station contribute to the event.
        // If not, continue with the next channel.
        if (!station.HasChannel(channelID) || !channel.IsActive())
          continue;

        if (!fAddNoiseFromSameTimeToBothChannels) {
          // choose trace and shift for each channel independently
          traceIndex = CLHEP::RandFlat::shootInt(fRandomEngine, 0, fNoiseTraces[channelID].size());
          if (fTraceShifting) {
            sampleIndexNoiseStart = CLHEP::RandFlat::shootInt(fRandomEngine, 0, 25) *
                                    fNoiseTraces[channelID][traceIndex].size() / 25;
          }
        }

        // Get number of adc counts from RDetector
        const short bitDepth = rDetStation.GetChannel(channel.GetId()).GetBitDepth();
 //       const short numberOfADCCounts = pow(2., bitDepth) - 1;
        double minVoltage = rDetStation.GetChannel(channel.GetId()).GetMinVoltage();
        double maxVoltage = rDetStation.GetChannel(channel.GetId()).GetMaxVoltage();
        double voltagePerCount = (maxVoltage - minVoltage) / pow(2, bitDepth);
        
        // Finally adding the noise event to the processed pedestal free event.
        const revt::ChannelADCTimeSeries timeSeriesOrig = timeSeries;

        const unsigned int traceLength = fNoiseTraces[channelID][traceIndex].size();
        for (unsigned int i = 0; i < timeSeries.GetSize(); ++i) {
          // get (shifted) adc count
          unsigned int sampleIndexNoise = (sampleIndexNoiseStart + i) % traceLength;

          const int noiseADCValue = fNoiseTraces.at(channelID).at(traceIndex).at(sampleIndexNoise);

          // // account for gain correction (if wanted)
          // const double mean = utl::TraceAlgorithm::Mean(timeSeriesOrig, 0, timeSeriesOrig.GetSize() - 1);
          // const int int_mean = int(std::abs(mean));
          // const int noiseADCValue = int(tmpADCValue * pow(10.0, fGainFactor / 20) + int_mean * (1 - pow(10.0, fGainFactor/ 20)));

          int val = timeSeries[i] + noiseADCValue;

          if (fReplaceDataWithNoise)
            val = noiseADCValue;

          // ADC count can not be out of range. This channel has to be flagged later.
          if (val > maxVoltage/voltagePerCount -1) {
            val = maxVoltage/voltagePerCount -1;
            infostr.str("");
            infostr << "ADC overflow (ADC = " << val << ", bin: " << i << ").\tStation "
                << channel.GetStationId() << ")\tnoise = " << noiseADCValue << "\tsignal = " << timeSeries[i];
            INFODebug(infostr);
          }
          if (val < minVoltage/voltagePerCount) {
            val = minVoltage/voltagePerCount;
            infostr.str("");
            infostr << "ADC underflow (ADC = " << val << ", bin: " << i << ").\tStation "
                    << channel.GetStationId() << ")\tnoise = " << noiseADCValue << "\tsignal = " << timeSeries[i];
            INFODebug(infostr);
          }

          timeSeries[i] = val;
        }
      }
    }

    return eSuccess;
  }

  fwk::VModule::ResultFlag
  RdChannelASCIINoiseImporterRD::Finish()
  {
    return eSuccess;
  }

  void
  RdChannelASCIINoiseImporterRD::ReadNoiseFile(const std::string fileName, const unsigned int channelId) {
    std::string line;
    std::ifstream noiseFile(fileName);
    if (!noiseFile)
      ERROR("Could not imported noise from asci file.");

    std::vector<std::vector<int>> noiseTracesChannel;
    if (noiseFile.is_open())
    {
      while (getline(noiseFile, line))
      {
        if (line[0] == '#') // skip any comments
          continue;

        std::vector<std::string> lineSplited;
        boost::split(lineSplited, line, [](char c) { return c == ','; });
        std::vector<int> noiseTrace;
        // assumes that first 2 entries of each row is something else than trace
        for (unsigned int column = fSkipTrailingColumnsInNoiseFile; column < lineSplited.size(); column++)
          noiseTrace.push_back(std::stoi(lineSplited[column]));

        noiseTracesChannel.push_back(noiseTrace);
      }
      noiseFile.close();
    }
    fNoiseTraces.at(channelId) = noiseTracesChannel;
  }

}