RdChannelNoiseImporter_RD.cc

Go to the documentation of this file.

#include "RdChannelNoiseImporter_RD.h"

#include "TTree.h"
#include "TFile.h"
#include "TKey.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 <fwk/RandomEngineRegistry.h>
#include <fwk/CentralConfig.h>

#include <CLHEP/Random/RandFlat.h>


using namespace fwk;
using namespace revt;


namespace RdChannelNoiseImporter_RD {

  VModule::ResultFlag
  RdChannelNoiseImporter_RD::Init()
  {
    utl::Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelNoiseImporter_RD");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("NoiseFilePath").GetData(fNoiseFilePath);
    topBranch.GetChild("AddNoiseFromSameTimeToBothChannels").GetData(fAddNoiseFromSameTimeToBothChannels);
    topBranch.GetChild("AllowTraceShifting").GetData(fTraceShifting);
    topBranch.GetChild("ReplaceDataWithNoise").GetData(fReplaceDataWithNoise);

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

    return eSuccess;
  }


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

    ReadFromRootFile();
    if (fAddNoiseFromSameTimeToBothChannels) {
      if (fNoiseTraces.first.size() != fNoiseTraces.second.size() || fNoiseTraces.first.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.first.size() << " traces for channel 0 and "
            << fNoiseTraces.second.size() << " traces for channel 1"
            << "\n\tAddNoiseFromSameTimeToBothChannels : " << fAddNoiseFromSameTimeToBothChannels
            << "\n\tAllow trace shifting               : " << fTraceShifting;

    infostr << std::endl;
    INFO(infostr);

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

    // 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.first.size());
      unsigned int sampleIndexNoiseStart = 0;

      // Traces have a length 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.first[traceIndex].GetSize() / 25;
      }

      const rdet::Station& rDetStation = rDetector.GetStation(stationID);
      for (auto& channel : station.ChannelsRange()) {
        auto& 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;

        std::vector<ChannelADCTimeSeries> tempVec;
        if (channelID == 0) {
          tempVec = fNoiseTraces.first;
        } else {
          tempVec = fNoiseTraces.second;
        }

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

        // Get number of adc counts from RDetector
        const short bitDepth = rDetStation.GetChannel(channel.GetId()).GetBitDepth();
        const short numberOfADCCounts = pow(2., bitDepth) / 2;

        // Finally adding the noise event to the processed pedestal free event.
        const ChannelADCTimeSeries timeSeriesOrig = timeSeries;

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

          const int noiseADCValue = tempVec.at(traceIndex).At(sampleIndexNoise);

          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 > numberOfADCCounts - 1) {
            val = numberOfADCCounts - 1;
            infostr.str("");
            infostr << "ADC overflow (ADC = " << val << ", bin: " << i << ").\tStation "
                << channel.GetStationId() << ")\tnoise = " << noiseADCValue << "\tsignal = " << timeSeries[i];
            INFODebug(infostr);
          }

          if (val < -1 * numberOfADCCounts) {
            val = -1 * numberOfADCCounts;
            infostr.str("");
            infostr << "ADC overflow (ADC = " << val << ", bin: " << i << ").\tStation "
                << channel.GetStationId() << ")\tnoise = " << noiseADCValue << "\tsignal = " << timeSeries[i];
            INFODebug(infostr);
          }

          timeSeries[i] = val;
        }
      }
    }

    return eSuccess;
  }


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


  VModule::ResultFlag
  RdChannelNoiseImporter_RD::ReadFromRootFile()
  {
    std::ifstream noiseFile(fNoiseFilePath);
    if (!noiseFile) {
      ERROR("Given noise file not existent.");
      return eFailure;
    }

    TFile* fileIn = TFile::Open(fNoiseFilePath.c_str());
    TList* keys = fileIn->GetListOfKeys();
    int numTrees = 0;
    TKey* key;
    TIter next(keys);
    while ((key = (TKey*)next())) {
      TObject* obj = key->ReadObj();
      if (obj->InheritsFrom("TTree")) {
        numTrees++;
      }
    }

    for(int counter = 0; counter < numTrees -1; ++counter) {
      std::string treename = "Noisetree_" + std::to_string(counter);
      TTree* noiseTree = (TTree*)fileIn->Get(treename.c_str());

      RdChannelNoiseImporter_RD::EventN noiseEvent;

      TBranch* traceCh0Branch = noiseTree->GetBranch("traceCh0");
      traceCh0Branch->SetAddress(&noiseEvent.traceCh0);
      TBranch* traceCh1Branch = noiseTree->GetBranch("traceCh1");
      traceCh1Branch->SetAddress(&noiseEvent.traceCh1);
      /*
      // not available with daves noise files but needed in the future
      TBranch* stationIdBranch = noiseTree->GetBranch("stationId");
      stationIdBranch->SetAddress(&noiseEvent.stationId);
      TBranch* secBranch = noiseTree->GetBranch("sec");
      secBranch->SetAddress(&noiseEvent.sec);
      TBranch* nanosecBranch = noiseTree->GetBranch("nanosec");
      nanosecBranch->SetAddress(&noiseEvent.nanosec);
      TBranch* evtnoBranch = noiseTree->GetBranch("evtno");
      evtnoBranch->SetAddress(&noiseEvent.evtno);
      noiseTree->GetEntry(0);
      */

      ChannelADCTimeSeries noiseTraceCh0;
      ChannelADCTimeSeries noiseTraceCh1;

      for (const auto& bin : noiseEvent.traceCh0) {
        noiseTraceCh0.PushBack(bin);
      }
      for (const auto& bin : noiseEvent.traceCh1) {
        noiseTraceCh1.PushBack(bin);
      }

      fNoiseTraces.first.push_back(noiseTraceCh0);
      fNoiseTraces.second.push_back(noiseTraceCh1);
    }

    return eSuccess;
  }

}