RdVirtualStationNoiseImporter.cc

Go to the documentation of this file.

#include "RdVirtualStationNoiseImporter.h"

#include <evt/Event.h>
#include <evt/RadioSimulation.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/RandomEngine.h>
#include <utl/UTCDateTime.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 <TTree.h>
#include <TBranch.h>

#include <RadioFileIO.h>

#include <CLHEP/Random/RandFlat.h>

using namespace std;
using namespace evt;
using namespace det;
using namespace fwk;
using namespace revt;

namespace RdVirtualStationNoiseImporter {

fwk::VModule::ResultFlag
RdVirtualStationNoiseImporter::Init()
{

  utl::Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch("RdVirtualStationNoiseImporter");
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  INFODebug("");

  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("StationIDList").GetData(fStationIDList);
  topBranch.GetChild("RandomNoiseTimeIntervalStart").GetData(fRandomNoiseTimeIntervalStart);
  topBranch.GetChild("RandomNoiseTimeIntervalStop").GetData(fRandomNoiseTimeIntervalStop);
  topBranch.GetChild("CheckTimeWindow").GetData(fCheckTimeWindow);
  topBranch.GetChild("TimeWindow").GetData(fMaxTimeDifference);
  topBranch.GetChild("AllowTraceShifting").GetData(fTraceShifting);
  topBranch.GetChild("GainFactor").GetData(fGainFactor);
 
  topBranch.GetChild("ReplaceDataWithNoise").GetData(fReplaceDataWithNoise);

  topBranch.GetChild("TimeIntervalInFile").GetData(fTimeIntervalInFile);

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

  ostringstream msg;
  if (fStationIDList.empty()) {
    ERROR("Station list is empty, no noise can be added!");
  } else {
    msg << "\n\tNoise from the following stations is added: ";
    for (auto& id : fStationIDList)
      msg << id + 100 << " ";
    msg << std::endl;
    INFOFinal(msg);
  }


  // read in _additional_ bad period entries
  msg.str("");
  msg << "Additional Bad Period entries:";
  for (utl::Branch b = topBranch.GetFirstChild(); b; b = b.GetNextSibling()) {
    if (b.GetName() != "BadStation")
      continue;
    const size_t stationId = det::VManager::FindComponent<size_t>("id", b.GetAttributes());
    const string reason = det::VManager::FindComponent<string>("reason", b.GetAttributes());

    vector<utl::TimeStamp> range;
    b.GetData(range);
    fRejectMap.insert(make_pair(range[0], make_pair(range[1], stationId)));
    msg << "\n\tregistered: id = " << stationId << " [" << range[0] << ", " << range[1] << "]";

  }
  INFOFinal(msg);

  return eSuccess;
}


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

  utl::UTCDateTime eventTime(event.GetHeader().GetTime());
  const utl::TimeStamp origTimeStamp = event.GetHeader().GetTime();

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

  ostringstream infostr;
  ostringstream warnstr;
  ostringstream debugstr;

  string NoiseImportfilename;
  utl::TimeStamp noiseTimeStamp;
  utl::UTCDateTime noiseUTCDateTime;
  utl::UTCDateTime noiseFileStartingDay;
  long noiseTimeGPSSecond;

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

  // Reading in a noise-event and checking if requested stations are present.
  bool fileOK = false;
  RadioFileIO* RadioFileIORead = nullptr;
  TFile* Readfile = nullptr;
  TTree* IOLibTree = nullptr;
  TBranch* IObranch = nullptr;
  // TBranch* evtSecBranch;
  vector<int> closeEventIds;
  vector<pair<int,vector<int>>> EventListWithStations;
  int tries = 0;
  while (!fileOK) {
    tries += 1;

    if (tries > 1 && fNoiseFileSelection == "Manually") {
      ERROR("Manually file sleection do not try a second time... aborting!");
      return eFailure;
    }
    
    // If the Noise file is selected based on a random event time
    if (fNoiseFileSelection == "Randomly") {
      const long noiseTimeIntervalStartGPSSecond = fRandomNoiseTimeIntervalStart.GetGPSSecond();
      const long noiseTimeIntervalStopGPSSecond = fRandomNoiseTimeIntervalStop.GetGPSSecond();
      noiseTimeGPSSecond =
        CLHEP::RandFlat::shootInt(fRandomEngine, noiseTimeIntervalStartGPSSecond, noiseTimeIntervalStopGPSSecond);
      
      const unsigned long secondsPerDay = 24 * 60 * 60;
      noiseTimeStamp = utl::TimeStamp(noiseTimeGPSSecond);
      
      // already here check with BP DB so that on does not have to open the file.
      det.Update(noiseTimeStamp);  // update for BadPeriod DP
      TimeRangeMap::IdCollection badIds = fRejectMap.GetBadStations(noiseTimeStamp);  // manual BP
      debugstr.str("");
      bool dont_skip = false;
      for (auto it = fStationIDList.begin() ; it != fStationIDList.end(); ++it) {
        if (!rDetector.GetBadStationReason(*it) && !(std::find(badIds.begin(), badIds.end(), *it) != badIds.end())) 
          dont_skip = true; // if at least one station is not in BP Database.
      }
      if (!dont_skip) {
        INFOIntermediate("Random time stamp for noise import is rejected, because of Bad Period. "
                         "Look for next file... ");
        continue;
      }

      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;
      
      if (!GetNoiseFileNameAtTime(noiseUTCDateTime, noiseFileStartingDay, NoiseImportfilename)) {
        if (tries > 9) {
          ERROR("Failed to fetch noise file too many times. Continue loop...");
          return eContinueLoop;
        }
        infostr.str("");
        infostr << "Attempt " << tries << ": Failed to import noise from file " << NoiseImportfilename;
        infostr << ". please check if the path and suffix/prefix are set properly.\nTrying next random time";
        INFOFinal(infostr);
        continue;
      }
      infostr << "Importing noise from file " << NoiseImportfilename;
      INFOIntermediate(infostr);
    }

    RadioFileIORead = new RadioFileIO();
    Readfile = new TFile(NoiseImportfilename.c_str(), "READ");
    IOLibTree = (TTree*)Readfile->Get("AERAIOTree");
    IObranch = IOLibTree->GetBranch("AERAIO");
    // evtSecBranch = IOLibTree->GetBranch("vAERAevents.EVseconds");
    IObranch->SetAddress(&RadioFileIORead);

    unsigned int NumberEvents = IOLibTree->GetEntries();

    // Search in a given interval around event time, noise-event is randomly selected from candidates
    // get list of event numbers that are within the allowed time interval
    vector<int> notBadStations;
    unsigned int stationsInTotal = 0;
    vector<int> infoList;
    EventListWithStations.clear();
    debugstr << "Number of events: " << NumberEvents << endl
             << "TimeStamp:" << noiseTimeStamp << endl;
    INFODebug(debugstr);

    // determine number of traces which are needed
    evt::ShowerSimData& simshow = event.GetSimShower();
    evt::RadioSimulation& rsim = simshow.GetRadioSimulation();
    
    // Number of requiered noise traces, reduces when allow to shift traces 
    unsigned int requiered_number = rsim.GetNumPulses();
    if (fTraceShifting)
      requiered_number = (unsigned int)requiered_number / 6 + 1;
    
    // loop over all events in selected file
    int idx = -1;
    // if fCheckTimeWindow is false only ask 5 events for data within one file
    // if no stations where found go the next file.
    unsigned int n_tried_events_in_file = 0;  
    do {
      ++idx;

      debugstr.str("");
      warnstr.str("");

      // find the index to get the event
      int used_index = 0;
      if (fCheckTimeWindow) {
        used_index = idx;  // here could be a smarter code, such that it is not iterated to much
      } else {
        // this is done so that not always the same time (beginning of file) is used for fCheckTimeWindow = 0
        used_index = CLHEP::RandFlat::shootInt(fRandomEngine, 0, NumberEvents - 1);
        n_tried_events_in_file++;
        if (n_tried_events_in_file > 5 && stationsInTotal == 0)
          break;
      }

      // get event
      IOLibTree->GetEntry(used_index);
      AERAevent* const aeraevent = RadioFileIORead->getAeraRadioFileEvent();
      if (!aeraevent->getAeraEventTriggerIsPeriodic())
        continue;

      if (fCheckTimeWindow) {
        utl::TimeStamp probeTimeStamp = utl::TimeStamp(aeraevent->getAeraEventSeconds());
        if (noiseTimeStamp < probeTimeStamp
        || std::abs(probeTimeStamp - noiseTimeStamp) > fMaxTimeDifference){
          debugstr << "Event " << used_index << ": (" << probeTimeStamp << ") is NOT in the specified range. continue. ";
          INFODebug(debugstr);
          continue;
        }

        debugstr << "event " << used_index << " is in the specified range. "
                 << "Now check how many stations are found (also check BP database). ";
        INFODebug(debugstr); 
      }

      if (EventHasStations(fStationIDList, aeraevent, rDetector).size() == 0) {
        INFODebug("No matched stations in event.. try next");
        continue;   
      }

      //check for bad period (again, for specific time)
      notBadStations.clear();
      TimeRangeMap::IdCollection badIds = fRejectMap.GetBadStations(noiseTimeStamp);  // manual BP
      debugstr.str("");
      for (auto it = fStationIDList.begin() ; it != fStationIDList.end(); ++it) {
        auto reason = rDetector.GetBadStationReason(*it);  // from DB
        if (reason) {
          debugstr << "\n\tStation rejceted: " << *it << ", reason: " << reason;
          continue; // skip if station is in BP Database.
        }
        if(!(std::find(badIds.begin(), badIds.end(), *it) != badIds.end())) {
          notBadStations.push_back(*it);
          debugstr << "\n\tStation accepted: " << *it;
        }
      }
      INFODebug(debugstr);

      // check
      if (notBadStations.empty()) {
        INFOFinal("all stations in event are bad.. try next");
        continue;
      }

      vector<int> matchedIDs = EventHasStations(notBadStations, aeraevent, rDetector);
      if (matchedIDs.size() == 0) {
        INFODebug("No matched good stations in event.. try next");
        continue;   
      }

      closeEventIds.push_back(used_index);
      EventListWithStations.emplace_back(used_index, matchedIDs);
      stationsInTotal += matchedIDs.size();
    
    } while (stationsInTotal < requiered_number and (unsigned int)idx < NumberEvents);  // still stop after one full iteration

    // check for existing noise trace
    if (closeEventIds.size() == 0 || stationsInTotal < requiered_number) {
      warnstr << "To few stations found within " << fMaxTimeDifference / utl::minute
              << " minutes has been found. (Checked time: " << fCheckTimeWindow << ")" << endl 
              << "Or no stations were found (in total: " 
              << stationsInTotal << "). File will be skipped";
      delete Readfile;
      delete RadioFileIORead;
      if (fNoiseFileSelection == "Manually") {
        ERROR("The manually selected noise file doesn't contain the defined stations...aborting!");
        return eFailure;
      }
      WARNING(warnstr);
      continue;
    }

    infostr.str("");
    infostr << "Found " << stationsInTotal << " stations in " << closeEventIds.size() 
            << " events. (If Shifting is used, each trace is used 6 times)";
    INFOFinal(infostr);

    fileOK = true;
  }

  // Mapping the indices of the noise event stations within the root container
  // to their station number for later use.
  unsigned int iterationCounter = 0;
  unsigned int internalStationCounter = 0;
  // Adding noise to the data, depending on if the station exists in both: noise-file and data.
  for (auto sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
    revt::Station& station = *sIt;
    const int stationID = station.GetId();
    unsigned int rootID = 0;
    unsigned int eventsRead = fTraceShifting ? (int)iterationCounter / 6 : iterationCounter;
    unsigned int nTimesUsed = fTraceShifting ? iterationCounter % 6: iterationCounter;
    
    // Fetching the whole selected noise-event from the root-file
    IOLibTree->GetEntry(EventListWithStations[eventsRead].first);
    AERAevent* const aeraevent = RadioFileIORead->getAeraRadioFileEvent();

    debugstr.str("");
    debugstr << "\n\tIteration " << iterationCounter << ", Event: " << eventsRead
             << " (is event no. " << closeEventIds[eventsRead] << ")"
             << " Traces were used " << nTimesUsed << " times."
             << "\n\tImporting Noise for station " << stationID
             << ", fetching aera event " << EventListWithStations[eventsRead].first;
    INFODebug(debugstr);
    
    // Printing some general information
    // std::ostringstream generalInfo;
    // generalInfo << "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(generalInfo);

    // get stations
    vector<int> stationsInEvent = EventListWithStations[eventsRead].second;

    // 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 };  
    debugstr.str("");
    const int tmpRootID = stationsInEvent[internalStationCounter];
    debugstr << "\n\tevent station counter " << internalStationCounter << " (index) of " << stationsInEvent.size() << " (size)"
             << "\n\tEvent: " << eventsRead << " of " << EventListWithStations.size();
    INFODebug(debugstr);


    aerarootiostation::Station* const rootIoReplacementStation = aeraevent->getAeraEventStation(tmpRootID);
    const int tmpReplacementStationID = rootIoReplacementStation->getAeraStationLsID() + 100;  //FIXME: Station Ids in RadioFile and RDetector differ by 100. Do not compare station Ids but rather station names (AERA_XXX).

    debugstr << "\n\tChecking station " << tmpReplacementStationID
             << " if it is a suitable replacement (this corresponds to station number "
             << tmpRootID << " in rootio file)";
    INFODebug(debugstr);

    // for next iteration
    if (stationsInEvent.size() > internalStationCounter + 1) {  // next stations
      internalStationCounter += 1;
    }
    else { // next event (resp. same event shifted)
      internalStationCounter = 0;
      iterationCounter += 1;
    }
    
    if (eventsRead > EventListWithStations.size()-2) {
      WARNING("Repeat loop ...");
      iterationCounter = 0; //repeat if not enough stations are present....
    }

    // first check if station present in the data is also available in the detector description
    const std::vector<int> stationIdVector = rDetector.GetFullStationList();
    if (std::find(stationIdVector.begin(), stationIdVector.end(), tmpReplacementStationID) == stationIdVector.end()) {
      infostr.str("");
      infostr << "Station " << tmpReplacementStationID
           << " is present in data but not in the detector description";
      INFOIntermediate(infostr);
      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 (rdet::Station::ChannelIterator chIt = rDetector.GetStation(station.GetId()).ChannelsBegin();
         chIt != rDetector.GetStation(station.GetId()).ChannelsEnd(); ++chIt) {
      debugstr.str("");
      debugstr << "checking channel " << chIt->GetId();
      INFODebug(debugstr);
      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()) {
        debugstr.str("");
        debugstr << "channel " << chIt->GetId() << " is not active, continue with next channel";
        INFODebug(debugstr);
        continue;
      }

      const rdet::Channel& detChannel = *chIt;
      bool similarChannelFound = false;  //used to break the loop when we found a channel
      for (rdet::Station::ChannelIterator replacementChIt = replacementStation.ChannelsBegin();
           (replacementChIt != replacementStation.ChannelsEnd()) && !similarChannelFound;
           ++replacementChIt) {
        const rdet::Channel& replacementChannel = *replacementChIt;
        // 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) {
          // check is channel and antenna type are the same
          if ((detChannel.GetChannelType() == replacementChannel.GetChannelType()) &&
          (detChannel.GetAntennaTypeName() == replacementChannel.GetAntennaTypeName())) {
            // check if sampling frequency is the same
            if (std::abs(detChannel.GetSamplingFrequency() - replacementChannel.GetSamplingFrequency()) < 0.01) {
              // check if ADC is the same
              if (detChannel.GetBitDepth() == replacementChannel.GetBitDepth()) {
                similarChannelFound = true;
                replacementChannelId[chIt->GetId() - 1] = replacementChannel.GetId();
                infostr.str("");
                infostr << "found replacement for channel " << chIt->GetId()
                 << ": it is channel " << replacementChannel.GetId() << " of station "
                 << replacementStation.GetId() << " (" << tmpReplacementStationID << ")";
                INFOIntermediate(infostr);
              }
            }
          }
        }
      }

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

    if (stationIsSuitable) {
      infostr.str("");
      infostr << "Suitable noise station was found for station " << stationID
          << ". Will use station " << tmpReplacementStationID << " in the noise file";
      INFOIntermediate(infostr);
      rootID = tmpRootID;
    } else {
      debugstr.str("");
      debugstr << "Station " << tmpReplacementStationID
           << " is not a suitable replacement for station " << stationID
           << " I think this shouldn't have happened... ";
      INFODebug(debugstr);
      continue;
    }
    // end looking for replacement station

    // Storing which channels of the station contribute to the noise-event.
    const int noiseChMask = aeraevent->getAeraEventStation(rootID)->getAeraStationChannelMask();
    const bool noiseChannelExists[4] = {
      bool(noiseChMask & 1),
      bool(noiseChMask & 2),
      bool(noiseChMask & 4),
      bool(noiseChMask & 8)
    };
    

    /* the channel numbering of the replacement station might be different
     * -> resort the  */
//    bool tmpNoiseChannelExists[4] = noiseChannelExists;
//    noiseChannelExists[replacementChannelId[0]-1] = tmpNoiseChannelExists[0];
//    noiseChannelExists[replacementChannelId[1]-1] = tmpNoiseChannelExists[1];
//    noiseChannelExists[replacementChannelId[2]-1] = tmpNoiseChannelExists[2];
//    noiseChannelExists[replacementChannelId[3]-1] = tmpNoiseChannelExists[3];

    std::ostringstream debug;
    debug << "ChannelMask from noise-station" << noiseChMask << ":" << noiseChannelExists[0]
          << noiseChannelExists[1] << noiseChannelExists[2] << noiseChannelExists[3];
    INFODebug(debug);
    
    //if (!usingReplacementStation) {
    for (revt::Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {

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

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

      // 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 = cIt->GetId();
      channelIdInNoiseFile = replacementChannelId[cIt->GetId() - 1];
      debug.str("");
      debug << "original channel id " << cIt->GetId() << " cooresponds to channel id "
            << channelIdInNoiseFile << " in noise file";
      INFODebug(debug);
      
      if (!noiseChannelExists[channelIdInNoiseFile - 1]) {
        cIt->SetNotActive();
        infostr.str("");
        infostr << "The channel " << cIt->GetId() << " of station " << stationID
                << " exists for data but not for noise. ---> Setting channel not active!";
        INFOIntermediate(infostr);
        continue;
      }

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

      // Finally adding the noise event to the processed pedestal free event.
      const revt::ChannelADCTimeSeries timeSeriesOrig = timeSeries;
      revt::ChannelADCTimeSeries pedestalFree = timeSeries;
      for (unsigned int i = 0; i < timeSeries.GetSize(); ++i) {
        //Removing pedestal of event time-series, in order to add the pedestal-free event to a noise-trace.
        RemovePedestal(pedestalFree);
        
        // get (shifted) adc count
        unsigned int j = fTraceShifting ? (i + nTimesUsed * noiseSeriesLength / 6) % noiseSeriesLength : i;
        const int tmpADCValue = 
            aeraevent->getAeraEventStation(rootID)->getAeraStationADC(channelIdInNoiseFile - 1)->GetADCValues(j);
        
        // 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)));

        // 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) {
          warnstr.str("");
          warnstr << "ADC overflow (ADC = " << val << ", \tbin " << i << " , (shifted-> bin is actually: " << j << ").\tstation "
              << cIt->GetStationId() << ")\tnoise = " << noiseADCValue << "\tsignal = "
              << timeSeries[i] << "\tsignal without pedestal = " << pedestalFree[i];
          WARNING(warnstr);
          val = numberOfADCCounts;
        }
        
        if (val < 0) {
          warnstr.str("");
          warnstr << "ADC underflow (ADC = " << val << ", \tbin " << i << " , (shifted-> bin is actually: " << j << ".\tstation "
              << cIt->GetStationId() << ")\tnoise = " << noiseADCValue << "\tsignal = "
              << timeSeries[i] << "\tsignal without pedestal = " << pedestalFree[i];
          WARNING(warnstr);
          val = 0;
        }
        timeSeries[i] = val;

      }
    }  // for Channel
    
    if (fTraceShifting)
      nTimesUsed += 1;
    
    int NotActiveChannels = 0;
    for (revt::Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {
      revt::Channel& channel = *cIt;
      if (!station.HasChannel(cIt->GetId()) || !channel.IsActive()) {
        ++NotActiveChannels;
      }
    }

    if (NotActiveChannels > 2) {
      warnstr.str("");
      warnstr << "Station " << stationID << " has more than two not active channels "
              << "---> Rejecting whole station!";
      WARNING(warnstr);
      station.SetRejectedReason(revt::eManuallyRejected);
    }
  } // for Station

  // Some clean up:
  delete Readfile;
  delete RadioFileIORead;
  
  // in case we changed the detector time when adding the noise
  // we now change it back to its original value
  det.Update(origTimeStamp);


  return eSuccess;
}


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

  INFO(out);
  
  return eSuccess;
}


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


void
RdVirtualStationNoiseImporter::RemovePedestal(std::vector<int>& timeSeries)
{
  double mean = 0;
  unsigned int N = 0;
  for (std::vector<int>::iterator it = timeSeries.begin(); it != timeSeries.end(); ++it) {
    mean += *it;
    ++N;
  }
  mean /= N;
  const int int_mean = int(std::abs(mean));
  for (std::vector<int>::iterator it = timeSeries.begin(); it != timeSeries.end(); ++it)
    *it -= int_mean;
}


bool
RdVirtualStationNoiseImporter::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();
  }
  if (fFilenameSuffix != "0")
      filename << fFilenameSuffix;
  filename << ".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 << "."
            << std::endl
            << "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);
  }
  return true;
}

vector<int>
RdVirtualStationNoiseImporter::EventHasStations(vector<int> stationID, AERAevent* const aeraevent, const rdet::RDetector& rDetector)
{
  ostringstream debugstr;
  std::vector<int> matchedID;

  // loop over all LS (stations) 
  for (int i = 0; i < aeraevent->getAeraEventLsCount(); ++i) {
    const int tmpRootID = i;
    aerarootiostation::Station* const rootIoReplacementStation = aeraevent->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;  

    for (std::vector<int>::iterator matchedIdIt = stationID.begin(); matchedIdIt != stationID.end(); ++matchedIdIt) {
      // match
      if (*matchedIdIt == tmpReplacementStationID - 100) {
        
        //check
        if (rootIoReplacementStation->getNumberOfAvailableAERAStationADCs() == 0
            || rootIoReplacementStation->getAeraStationStatus() == 0) {
          debugstr.str("");
          debugstr << "number of ADS is " << rootIoReplacementStation->getNumberOfAvailableAERAStationADCs()
              << "\nStation Status is " << rootIoReplacementStation->getAeraStationStatus() << " --> skip " << endl;
          INFODebug(debugstr);
          continue;
        }
        
        /*
        cout << "matched ID: " << *matchedIdIt << endl;
        cout << "station in Ev.: " << tmpRootID << endl;
        cout << "station trigger is periodic = " << rootIoReplacementStation->getAeraStationTriggerIsPeriodic() << endl;
        cout << "station status = " << rootIoReplacementStation->getAeraStationStatus() << endl;
        cout << "station message ADCcount = " << rootIoReplacementStation->getAeraStationMessageADCcount() << endl;
        cout << "# ADCs: " << rootIoReplacementStation->getNumberOfAvailableAERAStationADCs() << endl;
        cout << "ADC values ch 0: " << rootIoReplacementStation->getAeraStationADC(0)->GetADCValuesSize() << endl;
        cout << "ADC values ch 1: " << rootIoReplacementStation->getAeraStationADC(1)->GetADCValuesSize() << endl;
        cout << "ADC values ch 2: " << rootIoReplacementStation->getAeraStationADC(2)->GetADCValuesSize() << endl;
        cout << "ADC values ch 3: " << rootIoReplacementStation->getAeraStationADC(3)->GetADCValuesSize() << endl;
        cout << "now checking compatibitlity of stations..." << endl;
        */

        // first check if station present in the data is also available in the detector description
        const std::vector<int> stationIdVector = rDetector.GetFullStationList();
        if (std::find(stationIdVector.begin(), stationIdVector.end(), tmpReplacementStationID) == stationIdVector.end()) {
          debugstr.str("");
          debugstr << "Station " << tmpReplacementStationID << " is present in data but not in the detector description";
          INFODebug(debugstr);
          continue;
        }

        debugstr << "Matched station "<< tmpReplacementStationID << " (root id: " << tmpRootID << ")";

        // const rdet::Station& replacementStation = rDetector.GetStation(tmpReplacementStationID);
        // if (stationIsCompatible(testStation, replacementStation)) {
        //   if (!rootIoReplacementStation->getAeraStationStatus()) {
        //     cout << "station has no active channel... skip! " << endl;
        //     continue;
        //   }
        //   matchedID.push_back(tmpRootID);
        // }

        INFODebug(debugstr);
        matchedID.push_back(tmpRootID);
      }
    }
  }
  return matchedID;
}

/*
bool
RdVirtualStationNoiseImporter::IsCompatible(const revt::Station& testStation, const rdet::Station& replacementStation)
{
  std::ostringstream infostr;
  short replacementChannelId[4] = { -1, -1, -1, -1 };

  for (rdet::Station::ChannelIterator chIt = replacementStation.ChannelsBegin();
       chIt != replacementStation.ChannelsEnd(); ++chIt) {
    infostr.str("");
    infostr << "checking channel " << chIt->GetId();
    INFODebug(infostr);
    if (!testStation.HasChannel(chIt->GetId())) {  // if channel is not present in data, skip to next channel
      continue;
    }
    const revt::Channel& channel = testStation.GetChannel(chIt->GetId());
    if (!channel.IsActive()) {
      infostr.str("");
      infostr << "channel " << chIt->GetId() << " is not active, continue with next channel";
      INFODebug(infostr);
      continue;
    }
    const rdet::Channel& detChannel = *chIt;
    bool similarChannelFound = false;  //used to break the loop when we found a channel
    for (rdet::Station::ChannelIterator replacementChIt = replacementStation.ChannelsBegin();
         (replacementChIt != replacementStation.ChannelsEnd()) && !similarChannelFound;
         ++replacementChIt) {
      const rdet::Channel& replacementChannel = *replacementChIt;
      // 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) {
        // check is channel and antenna type are the same
        if ((detChannel.GetChannelType() == replacementChannel.GetChannelType()) &&
            (detChannel.GetAntennaTypeName() == replacementChannel.GetAntennaTypeName())) {
          // check if sampling frequency is the same
          if (std::abs(detChannel.GetSamplingFrequency() - replacementChannel.GetSamplingFrequency()) < 0.01) {
            // check if ADC is the same
            if (detChannel.GetBitDepth() == replacementChannel.GetBitDepth()) {
              similarChannelFound = true;
              replacementChannelId[chIt->GetId() - 1] = replacementChannel.GetId();
              infostr.str("");
              infostr << "found replacement for channel " << chIt->GetId()
               << ": it is channel " << replacementChannel.GetId() << " of station "
               << replacementStation.GetId() << endl;
              INFOIntermediate(infostr);
            }
          }
        }
      }
    }

  }
}
*/

}