RdChannelNoiseImporter_AERA.cc

Go to the documentation of this file.

#include "RdChannelNoiseImporter_AERA.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 <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/RandomEngine.h>
#include <utl/TimeStamp.h>

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

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

#include <TKey.h>
#include <TFile.h>

#include <CLHEP/Random/RandFlat.h>

namespace RdChannelNoiseImporter_AERA {


  fwk::VModule::ResultFlag
  RdChannelNoiseImporter_AERA::Init()
  {
    INFODebug("");

    const utl::Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch("RdChannelNoiseImporter_AERA");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("NoiseFileSelection").GetData(fNoiseFileSelection);
    topBranch.GetChild("NoiseFilePath").GetData(fNoiseFilePath);
    topBranch.GetChild("FilenamePrefix").GetData(fFilenamePrefix);
    topBranch.GetChild("FilenameSuffix").GetData(fFilenameSuffix);
    topBranch.GetChild("ExpectSubDirectories").GetData(fExpectSubDirectories);
    topBranch.GetChild("MaxTimeDifference").GetData(fMaxTimeDifference);
    topBranch.GetChild("NoiseEvtNr").GetData(fNoiseEvtNr);
    topBranch.GetChild("RejectStationsWithoutNoiseInformation").GetData(fRejectStationsWithoutNoiseInformation);
    topBranch.GetChild("ReplaceDataWithNoise").GetData(fReplaceDataWithNoise);
    topBranch.GetChild("RandomNoiseTimeIntervalStart").GetData(fRandomNoiseTimeIntervalStart);
    topBranch.GetChild("RandomNoiseTimeIntervalStop").GetData(fRandomNoiseTimeIntervalStop);

    //Add a / at the end of the path if missing and if the file is searched automatically in a directory
    if (fNoiseFileSelection == "Automatically") {
      if (fNoiseFilePath.compare(fNoiseFilePath.size() - 1, fNoiseFilePath.size(), "/") != 0)
        fNoiseFilePath.append("/");
    }

    topBranch.GetChild("EvtSelInNoiseFile").GetData(fEvtSelInNoiseFile);
    topBranch.GetChild("TimeIntervalInFile").GetData(fTimeIntervalInFile);
    fRandomEngine = &fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector).GetEngine();

    std::ostringstream info;
    info << "\n\t Reject stations without noise information (or try to find matching station with noise information): "
         << fRejectStationsWithoutNoiseInformation << std::endl;
    INFO(info);
    return (RdChannelNoiseImporter_AERA::checkSetup());
  }


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

    utl::UTCDateTime eventTime(event.GetHeader().GetTime());
    const utl::TimeStamp eventTimeStamp = event.GetHeader().GetTime();
    det::Detector& det = det::Detector::GetInstance();
    revt::REvent& rEvent = event.GetREvent();

    std::ostringstream sstr;
    sstr.str("");

    utl::TimeStamp noiseTimeStamp;

    // find noise file based on module config
    std::string NoiseImportfilename = GetNoiseFileName(eventTime, noiseTimeStamp);

    // Check if at selection method was successful
    if (NoiseImportfilename == "None") {
      return eContinueLoop;
    }

    // Reading in a noise-event from selected file
    RadioFileIO* RadioFileIORead = new RadioFileIO();
    TFile* const Readfile = new TFile(NoiseImportfilename.c_str(), "READ");
    TTree* const IOLibTree = (TTree*)Readfile->Get("AERAIOTree");
    TBranch* const IObranch = IOLibTree->GetBranch("AERAIO");
    TBranch* const evtSecBranch = IOLibTree->GetBranch("vAERAevents.EVseconds");
    IObranch->SetAddress(&RadioFileIORead);
    unsigned int NumberEvents = IOLibTree->GetEntries();
    std::map<int, int> IDtoNumber; //maps rootFileStationNumber on LS_ID

    // from here on different selection methods are applied, which are ordered hierarchically
    // in case several methods are selected, the last one used will determine the event number
    const int EvtNum = GetNoiseEventNumber(event, RadioFileIORead, evtSecBranch, IOLibTree, NumberEvents, eventTime);

    // Check if at selection method was successfully applied
    if (EvtNum == -1) {
      ERROR("No noise event selected. Can not continue with import procedure. Skipping event!");
      return eContinueLoop;
    }

    // Fetching the whole selected noise-event from the root-file
    IOLibTree->GetEntry(EvtNum);
    AERAevent* const aeraevent = RadioFileIORead->getAeraRadioFileEvent();
    // update to detector time to noise event time
    if (fNoiseFileSelection == "RandomlyInTimeInterval") {
      noiseTimeStamp.SetGPSTime(aeraevent->getAeraEventSeconds(), 0);
      det.Update(noiseTimeStamp);
    }

    // Printing some general information
    sstr << "General info concerning the noise-event selection:\n\tGPSSecond of the chosen noise-event: "
         << aeraevent->getAeraEventSeconds() << "\n\tGPSSecond of the processed event: "
         << eventTime.GetTimeStamp().GetGPSSecond() << '\n';
    INFOIntermediate(sstr);

    // Checking the time difference between processed data-event and the selected noise-event,
    // if noise-event not set fixed or random by xml-config.
    if (fEvtSelInNoiseFile == "ByEstimationFromEvtNr" || fEvtSelInNoiseFile == "ByTimestamp" ||
        fEvtSelInNoiseFile == "RandomlyInTimeWindow") {
      const unsigned int noiseEvtNanoSec = aeraevent->getAeraEventNanoSeconds();
      const unsigned int noiseEvtSec = aeraevent->getAeraEventSeconds();
      const utl::TimeStamp noiseEvtTimestamp(noiseEvtSec, noiseEvtNanoSec);
      if (std::abs(noiseEvtTimestamp - eventTimeStamp) > fMaxTimeDifference) {
        std::ostringstream warning;
        warning <<
          "The difference between the time of the estimated noise-event "
          "time and the currently processed event exceeds " << fMaxTimeDifference / utl::minute <<
          " minutes. Skipping event.";
        WARNING(warning);
        return eContinueLoop;
      }
    }

    // Mapping the indices of the noise event stations within the root container
    // to their station number for later use.
    for (int i = 0; i < aeraevent->getAeraEventLsCount(); ++i) {
      aerarootiostation::Station* const station = RadioFileIORead->getAeraRadioFileEvent()->getAeraEventStation(i);
      IDtoNumber.insert(std::pair<int, int>(station->getAeraStationLsID(), i));
    }

    // Getting rDetector after time might have been updated
    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();
      std::map<int, int>::iterator idIt;
      int rootID = 0;
      //set to true if the current station has no noise trace and the noise trace from a different station is used
      bool usingReplacementStation;
      // this array stores the channel association of the replacement station. This is necessary as numbering of the
      // channels might not be the same. However, in most cases channel 1 of the current station corresponds also to
      // channel 1 of the replacement station
      short replacementChannelId[4] = { -1, -1, -1, -1 };

      // Check if the event-station exists in the noise-file and reject it if fRejectStationsWithoutNoiseInformation == true.
      if ((idIt = IDtoNumber.find(stationID - 100)) == IDtoNumber.end()) { // if not in data:
        if (fRejectStationsWithoutNoiseInformation) {
          sstr.str("");
          sstr << "Did not find Station " << stationID << " in event number " << EvtNum << " of noise-file " <<
            NoiseImportfilename << "---> Excluding station!";
          INFOFinal(sstr);
          station.SetExcludedReason(revt::eNoNoiseDataAvailable);
          continue;
        } else {
          sstr.str("");
          sstr << "No noise information available for station " << stationID
               << ". Searching for other station to use its noise";
          INFOIntermediate(sstr);

          // if that bool is turned true noise is added otherwise the station is excluded.
          // in the following the loop is continued if it turns out that the tested candidate is not suitable
          // (and the bool is then not switched)
          bool suitableStationFound = false;
          // create randomized list of all LS ids of the root file
          std::vector<int>::size_type n = aeraevent->getAeraEventLsCount();
          std::vector<int> rootIds;
          rootIds.reserve(n);
          while (rootIds.size() < n) {
             int randInt = CLHEP::RandFlat::shootInt(fRandomEngine, 0, n);
             if (std::find(rootIds.begin(), rootIds.end(), randInt) == rootIds.end()) {
                rootIds.push_back(randInt);
             }
          }

          for (auto tmpRootID : rootIds) {
            aerarootiostation::Station* const rootIoReplacementStation =
              RadioFileIORead->getAeraRadioFileEvent()->getAeraEventStation(tmpRootID);
            //FIXME: Station Ids in RadioFile and RDetector differ by 100. Do not compare station Ids but rather station names (AERA_XXX)
            const int tmpReplacementStationID = rootIoReplacementStation->getAeraStationLsID() + 100;
            sstr.str("");
            sstr << "checking station " << tmpReplacementStationID
                 << " if it is a suitable replacement (this corresponds to station number "
                 << tmpRootID << " in rootio file)";
            INFODebug(sstr);
            const std::vector<int> stationIdVector = rDetector.GetFullStationList();
            // first check if station present in the data is also available in the detector description
            if (std::find(stationIdVector.begin(), stationIdVector.end(),
                          tmpReplacementStationID) == stationIdVector.end()) {
              sstr.str("");
              sstr << "Station " << tmpReplacementStationID
                   << " is present in data but not in the detector description";
              INFOIntermediate(sstr);
              continue;
            }

            // checking BP database
            unsigned long long int reason = rDetector.GetBadStationReason(tmpReplacementStationID);
            if (reason) {
              sstr.str("");
              sstr << "Station " << tmpReplacementStationID << " is excluded for bad period reason: " << reason;
              INFOFinal(sstr);
              continue;
            }

            const rdet::Station& replacementStation = rDetector.GetStation(tmpReplacementStationID);
            bool stationIsSuitable = true;

            // loop through the channels of the station and check if the station in the noise file has the same properties
            for (auto chIt = rDetector.GetStation(station.GetId()).ChannelsBegin();
                 chIt != rDetector.GetStation(station.GetId()).ChannelsEnd(); ++chIt) {
              sstr.str("");
              sstr << "checking channel " << chIt->GetId();
              INFODebug(sstr);
              if (!station.HasChannel(chIt->GetId())) { // if channel is not present in data, skip to next channel
                continue;
              }

              const revt::Channel& channel = station.GetChannel(chIt->GetId());
              if (!channel.IsActive()) {
                sstr.str("");
                sstr << "channel " << chIt->GetId() << " is not active, continue with next channel";
                INFODebug(sstr);
                continue;
              }

              const rdet::Channel& detChannel = *chIt;
              bool similarChannelFound = false; //used to break the loop when we found a channel

              for (auto replacementChIt = replacementStation.ChannelsBegin();
                   (replacementChIt != replacementStation.ChannelsEnd()) && !similarChannelFound; ++replacementChIt) {
                const rdet::Channel& replacementChannel = *replacementChIt;
                similarChannelFound = CheckMatchingChannel(detChannel, replacementChannel);

                if (similarChannelFound) {
                  replacementChannelId[chIt->GetId() - 1] = replacementChannel.GetId();
                  sstr.str("");
                  sstr << "found replacement for channel " << chIt->GetId() << ": it is channel "
                       << replacementChannel.GetId() << " of station " << replacementStation.GetId() << " ("
                       << tmpReplacementStationID << ")" << "\n " << detChannel.GetAntennaTypeName() << " "
                       << replacementChannel.GetAntennaTypeName();
                  INFOIntermediate(sstr);
                }
              }

              if (!similarChannelFound) {
                  stationIsSuitable = false;
              }
              
            }

            if (stationIsSuitable) {
              suitableStationFound = true;
              sstr.str("");
              sstr << "Suitable noise station was found for station " << stationID << ". Will use station "
                   << tmpReplacementStationID << " in the noise file";
              INFOIntermediate(sstr);
              rootID = tmpRootID;
              break;
            } else {
              sstr.str("");
              sstr << "Station " << tmpReplacementStationID << " is not a suitable replacement for station "
                   << stationID;
              INFODebug(sstr);
            }
          }

          if (!suitableStationFound) {
            sstr.str("");
            sstr << "Exclude station " << stationID << "!" << "\n\t It does not exist in event number " << EvtNum
                 << " of noise-file " << NoiseImportfilename
                 << ". Tried to use noise information from another station but no suitable station was found.";
            INFOFinal(sstr);
            station.SetExcludedReason(revt::eNoNoiseDataAvailable);
            fNumberOfRejectedStationsWithoutNoise++;
            continue;
          } else {
            usingReplacementStation = true;
          }

        } // end looking for replacement station

      } else {
        sstr.str("");
        sstr << "Found station " << station.GetId() << " in noise data.";
        // checking BP database
        unsigned long long int reason = rDetector.GetBadStationReason(station.GetId());
        if (reason) {
          station.SetExcludedReason(revt::eNoNoiseDataAvailable);
          sstr << " But was excluded because for bad period reason: " << reason;
          INFOFinal(sstr);
          fNumberOfRejectedStationsWithoutNoise++;
          continue;
        }

        INFOIntermediate(sstr);
        rootID = idIt->second;
        usingReplacementStation = false;
      }

      // Storing which channels of the station contribute to the noise-event.
      const int noiseChMask =
        RadioFileIORead->getAeraRadioFileEvent()->getAeraEventStation(rootID)->getAeraStationChannelMask();
      const bool noiseChannelExists[4] = {bool(noiseChMask & 1), bool(noiseChMask & 2),
                                          bool(noiseChMask & 4), bool(noiseChMask & 8)};
      sstr << "ChannelMask from noise-station" << noiseChMask << ":" << noiseChannelExists[0]
           << noiseChannelExists[1] << noiseChannelExists[2] << noiseChannelExists[3];
      INFODebug(sstr);

      for (auto& channel : station.ChannelsRange()) {
        // Checking if the channel of the processed station contribute to the event.
        // If not, continue with the next channel.
        if (!station.HasChannel(channel.GetId()) || !channel.IsActive())
          continue;

        revt::ChannelADCTimeSeries& timeSeries = channel.GetChannelADCTimeSeries();

        // Get number of adc counts from RDetector
        const short bitDepth = det.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId()).GetBitDepth();
        const short numberOfADCCounts = pow(2., bitDepth) - 1;
        // Checking if the noise of this channel can be added to the event. If not, reject channel.
        unsigned int noiseSeriesLength = 0;
        int channelIdInNoiseFile = channel.GetId();

        if (usingReplacementStation) {
          channelIdInNoiseFile = replacementChannelId[channel.GetId() - 1];
          sstr.str("");
          sstr << "original channel id " << channel.GetId() << " cooresponds to channel id " << channelIdInNoiseFile
               << " in noise file";
          INFODebug(sstr);
        }

        if (!noiseChannelExists[channelIdInNoiseFile - 1]) {
          channel.SetNotActive();
          sstr.str("");
          sstr << "The channel " << channel.GetId() << " of station " << stationID
               << " exists for data but not for noise. ---> Setting channel not active!";
          INFOIntermediate(sstr);
          continue;
        }

        noiseSeriesLength = RadioFileIORead->getAeraRadioFileEvent()->getAeraEventStation(rootID)->getAeraStationADC(
          channelIdInNoiseFile - 1)->GetADCValuesSize();
        if (noiseSeriesLength <= timeSeries.GetSize()) {
          channel.SetNotActive();
          sstr.str("");
          sstr << "The length " << noiseSeriesLength << " of the noise timeseries in channel " << channel.GetId()
               << " is shorter then the length " << timeSeries.GetSize() << " of data timeseries in station "
               << stationID << "---> Setting channel not active!";
          WARNING(sstr);
          continue;
        }

        // Finally adding the noise event to the processed pedestal free event.
        const revt::ChannelADCTimeSeries timeSeriesOrig = timeSeries;
        revt::ChannelADCTimeSeries pedestalFree = timeSeries;
        //Removing pedestal of event time-series, in order to add the pedestal-free event to a noise-trace.
        RemovePedestal(pedestalFree);
        for (unsigned int i = 0; i < timeSeries.GetSize(); ++i) {
          const int noiseADCValue =
            RadioFileIORead->getAeraRadioFileEvent()->getAeraEventStation(rootID)->getAeraStationADC(
              channelIdInNoiseFile - 1)->GetADCValues(i);
          // Decide if noise should be added to simulated data or replace it to create pure noise events
          int val = fReplaceDataWithNoise ? noiseADCValue : pedestalFree[i] + noiseADCValue;
          // Need to be done correctly, and info should be taken from detector geometry
          // WHY, it seems to be done properly??
          if (val > numberOfADCCounts) {
            sstr.str("");
            sstr << "ADC overflow (ADC = " << val << ", \tbin " << i << " , \tstation " << channel.GetStationId()
                 << ")\tnoise = " << noiseADCValue << "\tsignal = " << timeSeries[i] << "\tsignal without pedestal = "
                 << pedestalFree[i];
            WARNING(sstr);
            val = numberOfADCCounts;
          }

          if (val < 0) {
            sstr.str("");
            sstr << "ADC underflow (ADC = " << val << ", \tbin " << i << " , \tstation " << channel.GetStationId()
                 << ")\tnoise = " << noiseADCValue << "\tsignal = " << timeSeries[i] << "\tsignal without pedestal = "
                 << pedestalFree[i];
            WARNING(sstr);
            val = 0;
          }

          timeSeries[i] = val;
          sstr.str("");
          sstr << "Timeseries orig: " << timeSeriesOrig[i] << "    Timeseries after NoiseImporter " << timeSeries[i];
          INFODebug(sstr);
        }
      } // for Channel

      int notActiveChannels = 0;
      for (auto& channel : station.ChannelsRange()) {
        if (!station.HasChannel(channel.GetId()) || !channel.IsActive()) {
          ++notActiveChannels;
        }
      }

      if (notActiveChannels > 2) {
        sstr.str("");
        sstr << "Station " << stationID << " has more than two not active channels " << "---> Exclude station!";
        WARNING(sstr);
        station.SetExcludedReason(revt::eHardwareMismatch);
      }
    } // for Station

    // Some clean up:
    delete Readfile;
    delete RadioFileIORead;

    // in case we changed the detector time when adding the noise
    // (happens when the option <RandomlyInTimeInterval> is chosen)
    // we now change it back to its original value
    if (fNoiseFileSelection == "RandomlyInTimeInterval") {
      det.Update(eventTimeStamp);
    }

    return eSuccess;
  }


  fwk::VModule::ResultFlag
  RdChannelNoiseImporter_AERA::Finish()
  {
    std::ostringstream out;
    out << "\n\tNumber of stations rejected due to missing noise: " << fNumberOfRejectedStationsWithoutNoise;
    INFO(out);
    return eSuccess;
  }


  void
  RdChannelNoiseImporter_AERA::RemovePedestal(revt::ChannelADCTimeSeries& timeSeries)
  {
    const double mean = utl::TraceAlgorithm::Mean(timeSeries, 0, timeSeries.GetSize() - 1);
    const int intMean = int(std::abs(mean));
    for (auto it = timeSeries.Begin(); it != timeSeries.End(); ++it)
      *it -= intMean;
  }


  void
  RdChannelNoiseImporter_AERA::RemovePedestal(std::vector<int>& timeSeries)
  {
    const double sum = std::accumulate(timeSeries.begin(), timeSeries.end(), 0);
    const int length = timeSeries.size();
    const int intMean = int(std::abs(sum)/length);

    for (auto& e : timeSeries) {
      e -= intMean;
    }
  }


  fwk::VModule::ResultFlag
  RdChannelNoiseImporter_AERA::checkSetup()
  {
    // Check if a proper path to the noise file or directory is given.
    if (fNoiseFilePath == "/path/to/noise/file/") {
      std::ostringstream error;
      error << "The <NoiseFilePath> has not been set in the initfile";
      ERROR(error);
      return eFailure;
    }

    return eSuccess;
  }


  bool
  RdChannelNoiseImporter_AERA::GetNoiseFileNameAtTime(const utl::UTCDateTime& dateTime,
                                                      const utl::UTCDateTime& startingDay,
                                                      std::string& NoiseImportfilename)
  {
    std::ostringstream filename;
    std::ostringstream subdirectory;
    if (fExpectSubDirectories) {
      subdirectory << startingDay.GetYear() << std::setfill('0') << "/" << std::setw(2) << startingDay.GetMonth() << "/";
      if ((fTimeIntervalInFile == "Hour")) {
        subdirectory << std::setfill('0') << std::setw(2) << startingDay.GetDay() << "/";
      }
    }

    filename << fFilenamePrefix;
    filename << startingDay.GetYear() << "_" << std::setfill('0') << std::setw(2) << startingDay.GetMonth() << "_"
             << std::setfill('0') << std::setw(2) << startingDay.GetDay();
    if ((fTimeIntervalInFile == "Hour")) {
      filename << "_" << std::setfill('0') << std::setw(2) << dateTime.GetHour();
    }

    filename << fFilenameSuffix << ".root";
    NoiseImportfilename = fNoiseFilePath + subdirectory.str() + filename.str();
    // Checking if the noise-file exists.
    std::ifstream checkfile(NoiseImportfilename.c_str());
    if (!checkfile) {
      std::ostringstream warning;
      warning << "Cannot find the corresponding noise-file to the event from " << dateTime.GetYear()
              << std::setfill('0') << std::setw(2) << dateTime.GetMonth() << std::setfill('0') << std::setw(2)
              << dateTime.GetDay() << " at " << NoiseImportfilename << "." << "/n"
              << "Check if the desired file exists at the path specified by the variable " "\"NoiseFilePath\" in the initfile!";
      WARNING(warning);
      return false;
    } else {
      std::ostringstream info;
      info << "Reading from noise-file " << NoiseImportfilename << "." << std::endl;
      INFOFinal(info.str());
    }

    return true;
  }


  int
  RdChannelNoiseImporter_AERA::GetNoiseEventNumber(evt::Event& event, RadioFileIO* RadioFileIORead,
                                                   TBranch* const evtSecBranch, TTree* const IOLibTree,
                                                   unsigned int NumberEvents, utl::UTCDateTime eventTime)
  {
    int EvtNum = -1;
    std::ostringstream sstr;

    // Select noise-event by event number specified by xml-config
    if (fEvtSelInNoiseFile == "ByNrInFile") {
      EvtNum = fNoiseEvtNr;
      sstr.str("");
      sstr << "Noise-event number chosen in init-file is " << EvtNum << ".";
      INFOIntermediate(sstr);
    }

    // Select a random noise-event
    if (fEvtSelInNoiseFile == "Randomly") {
      EvtNum = CLHEP::RandFlat::shootInt(fRandomEngine, NumberEvents - 1);
      // Fetching the whole selected noise-event from the root-file
      IOLibTree->GetEntry(EvtNum);
      sstr.str("");
      sstr << "Randomly chosen noise-event number is " << EvtNum << " out of " << NumberEvents << " events.";
      INFOIntermediate(sstr);
    }

    // Estimate the noise-event according to the data-event
    if (fEvtSelInNoiseFile == "ByEstimationFromEvtNr") {
      if (fTimeIntervalInFile == "Day") {
        int hourOffset = 0;
        if (eventTime.GetHour() < 12) {
          hourOffset = 12;
        } else {
          hourOffset = -12;
        }
        EvtNum = int(((eventTime.GetHour() + hourOffset) * 60. + eventTime.GetMinute()) * (NumberEvents / 24. / 60.));
      }

      if (fTimeIntervalInFile == "Hour") {
        EvtNum = int(eventTime.GetMinute() * (NumberEvents / 60.));
      }

      sstr.str("");
      sstr << "Estimated noise-event number by the processed data-event time is " << EvtNum << ".";
      INFOIntermediate(sstr);
    }

    // Select the noise-event succeeding the processed data-event in time
    if (fEvtSelInNoiseFile == "ByTimestamp") {
      for (unsigned int EvtCounter = 0; EvtCounter < NumberEvents; ++EvtCounter) {
        evtSecBranch->GetEntry(EvtCounter);
        // Be CAREFUL: Only the event-second of aeraevent is filled and updated in this scope.
        // The other members of aeraevent (like ID,...) contain unknown/initial values. Do not call them!
        AERAevent* const aeraevent = RadioFileIORead->getAeraRadioFileEvent();
        const unsigned int noiseEvtSec = aeraevent->getAeraEventSeconds();
        const utl::TimeStamp eventTimestamp = event.GetHeader().GetTime();
        const double diff = utl::TimeStamp(noiseEvtSec) - eventTimestamp;
        if (diff > 0) {
          EvtNum = EvtCounter;
          break;
        }

        if (EvtCounter == NumberEvents - 1) {
          EvtNum = EvtCounter;
        }
      }

      sstr.str("");
      sstr << "Noise-event number " << EvtNum << " succeeds the processed data-event in time.";
      INFOIntermediate(sstr);
    }

    // Search in a given interval around event time, noise-event is randomly selected from candidates
    if (fEvtSelInNoiseFile == "RandomlyInTimeWindow") {
      // get list of event numbers that are within the allowed time interval
      std::vector<int> closeEventIds;
      for (unsigned int i = 0; i < NumberEvents; ++i) {
        evtSecBranch->GetEntry(i);
        // Be CAREFUL: Only the event-second of aeraevent is filled and updated in this scope.
        // The other members of aeraevent (like ID,...) contain unknown/initial values. Do not call them!
        AERAevent* const aeraevent = RadioFileIORead->getAeraRadioFileEvent();
        const unsigned int noiseGPSSeconds = aeraevent->getAeraEventSeconds();
        const utl::TimeStamp eventTimestamp = event.GetHeader().GetTime();
        const double diff = std::abs(eventTimestamp - utl::TimeStamp(noiseGPSSeconds));
        if (diff < fMaxTimeDifference) {
          closeEventIds.push_back(i);
        }
      }

      // check for existing noise trace
      if (closeEventIds.size() == 0) {
        sstr << "No noise trace within " << fMaxTimeDifference / utl::minute
             << " minutes has been found. Event will be skipped";
        INFOIntermediate(sstr);
        return -1;
      }

      // Fetch a random noise-event from the the closeEventIds-vector
      int randVecPos = CLHEP::RandFlat::shootInt(fRandomEngine, closeEventIds.size() - 1);
      EvtNum = closeEventIds[randVecPos];
    }

    return EvtNum;
  }


  std::string
  RdChannelNoiseImporter_AERA::GetNoiseFileName(const utl::UTCDateTime& eventTime, utl::TimeStamp& noiseTimeStamp)
  {
    std::ostringstream sstr;

    utl::UTCDateTime noiseFileStartingDay;
    utl::UTCDateTime noiseUTCDateTime;
    std::string NoiseImportfilename;

    if (fNoiseFileSelection == "RandomlyInTimeInterval") {
      const long noiseTimeIntervalStartGPSSecond = fRandomNoiseTimeIntervalStart.GetGPSSecond();
      const long noiseTimeIntervalStopGPSSecond = fRandomNoiseTimeIntervalStop.GetGPSSecond();
      int tries = 0;
      do {
        if (tries > 10) {
          sstr << "Tried more than 10 times to import noise from a randomly chosen file. Skip Event...";
          INFO(sstr);
          return "None";
        }

        long noiseTimeGPSSecond = CLHEP::RandFlat::shootInt(fRandomEngine, noiseTimeIntervalStartGPSSecond,
                                                            noiseTimeIntervalStopGPSSecond);
        const unsigned long secondsPerDay = 24 * 60 * 60;
        noiseTimeStamp = utl::TimeStamp(noiseTimeGPSSecond);
        noiseUTCDateTime = utl::UTCDateTime(noiseTimeStamp);
        const utl::TimeStamp noiseTimeStampMinus24h = utl::TimeStamp(noiseTimeGPSSecond - secondsPerDay);
        const utl::UTCDateTime noiseUTCDateTimeMinus24h = utl::UTCDateTime(noiseTimeStampMinus24h);
        noiseFileStartingDay = (noiseUTCDateTime.GetHour() < 12) ? noiseUTCDateTimeMinus24h : noiseUTCDateTime;
        tries++;
      } while (!GetNoiseFileNameAtTime(noiseUTCDateTime, noiseFileStartingDay, NoiseImportfilename));

      sstr << "Importing noise from file " << NoiseImportfilename;
      INFOIntermediate(sstr);
    }

    //If only a directory is given in init-file: Search for and check a suiting noise-file.
    if (fNoiseFileSelection == "Automatically") {
      const unsigned long secondsPerDay = 24 * 60 * 60;
      const utl::TimeStamp eventTimeStampMinus24h(eventTime.GetTimeStamp().GetGPSSecond() - secondsPerDay,
                                                  eventTime.GetTimeStamp().GetGPSNanoSecond());
      const utl::UTCDateTime eventTimeMinus24h(eventTimeStampMinus24h);

      // Taking into account, that daily recorded background-files
      // contain data from 12 noon of the starting day until 12 noon
      // of the following day.
      bool eventBefore12amUTC = false;
      if (eventTime.GetHour() < 12) {
        eventBefore12amUTC = true;
      }

      //  Using the file one day before if the processed event has taken place before 12 noon UTC.
      // But only if the files contain noise-data of a whole day.
      if (eventBefore12amUTC && (fTimeIntervalInFile == "Day")) {
        noiseFileStartingDay = eventTimeMinus24h;
      } else {
        noiseFileStartingDay = eventTime;
      }

      // Composing the path to the noise-file.
      if (!GetNoiseFileNameAtTime(eventTime, noiseFileStartingDay, NoiseImportfilename)) {
        return "None";
      }
    }

    // If the noise file is given explicitly in the init-file:
    if (fNoiseFileSelection == "Manually") {
      NoiseImportfilename = fNoiseFilePath;
    }

    return NoiseImportfilename;
  }


  bool
  RdChannelNoiseImporter_AERA::CheckMatchingChannel(const rdet::Channel& detChannel, const rdet::Channel& replacementChannel)
  {
    // check if orientation of the channel is the same
    if (std::abs(detChannel.GetOrientationZenith() - replacementChannel.GetOrientationZenith()) > 3 * utl::degree ||
        std::abs(detChannel.GetOrientationAzimuth() - replacementChannel.GetOrientationAzimuth()) > 3 * utl::degree)
      return false;

    // check is channel and antenna type are the same
    if ((detChannel.GetChannelType() != replacementChannel.GetChannelType()) ||
        (detChannel.GetAntennaTypeName() != replacementChannel.GetAntennaTypeName()))
      return false;

    // check if sampling frequency is the same
    if (std::abs(detChannel.GetSamplingFrequency() - replacementChannel.GetSamplingFrequency()) > 0.01)
      return false;

    // check if ADC is the same
    if (detChannel.GetBitDepth() != replacementChannel.GetBitDepth())
      return false;

    // replacement found
    return true;
  }

}