RdEventPreSelector.cc

Go to the documentation of this file.

#include "RdEventPreSelector.h"

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/TraceAlgorithm.h>
#include <utl/RadioTraceUtilities.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/EventTrigger.h>
#include <revt/StationRecData.h>
#include <revt/StationHeader.h>

#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h> // for SD zenith cut
#include <evt/ShowerFRecData.h> // for FD zenith cut

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

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>

#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string/split.hpp>


using namespace fwk;
using namespace utl;
using namespace evt;
using namespace revt;
using namespace fevt;
using std::vector;
using std::pair;
using std::set;
using std::string;
using std::ostringstream;
using std::fstream;
using std::min;
using std::max;


namespace RdEventPreSelector {

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

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("MinEventId").GetData(fMinEventId);
    topBranch.GetChild("MaxEventId").GetData(fMaxEventId);
    topBranch.GetChild("MinNumberOfStations").GetData(fMinNumberOfStations);
    topBranch.GetChild("MaxNumberOfStations").GetData(fMaxNumberOfStations);
    topBranch.GetChild("AllParticipatingStationIds").GetData(fAllParticipatingStationIds);
    topBranch.GetChild("AtLeastOneParticipatingStationIds").GetData(fAtLeastOneParticipatingStationIds);
    topBranch.GetChild("NotParticipatingStationIds").GetData(fNotParticipatingStationIds);
    topBranch.GetChild("UseTriggerInformation").GetData(fUseTriggerInformation);
    topBranch.GetChild("UseSelfTriggeredEvent").GetData(fUseSelfTriggeredEvent);
    topBranch.GetChild("UseExternallyTriggeredEvent").GetData(fUseExternallyTriggeredEvent);
    topBranch.GetChild("UseCalibrationTriggeredEvent").GetData(fUseCalibrationTriggeredEvent);
    topBranch.GetChild("UseScintillatorTriggeredEvent").GetData(fUseScintillatorTriggeredEvent);
    topBranch.GetChild("UseSDTriggeredEvent").GetData(fUseSDTriggeredEvent);
    topBranch.GetChild("UseGUITriggeredEvent").GetData(fUseGUITriggeredEvent);
    topBranch.GetChild("UseFDTriggeredEvent").GetData(fUseFDTriggeredEvent);
    topBranch.GetChild("UseHEATTriggeredEvent").GetData(fUseHEATTriggeredEvent);
    topBranch.GetChild("UseAERAletTriggeredEvent").GetData(fUseAERAletTriggeredEvent);
    topBranch.GetChild("UseAIRPLANETriggeredEvent").GetData(fUseAIRPLANETriggeredEvent);
    topBranch.GetChild("UsePeriodicTriggeredEvent").GetData(fUsePeriodicTriggeredEvent);
    topBranch.GetChild("UsePassThroughTriggeredEvent").GetData(fUsePassThroughTriggeredEvent);
    topBranch.GetChild("FirstAllowedEventDateTime").GetData(fAllowedTimeBegin);
    topBranch.GetChild("LastAllowedEventDateTime").GetData(fAllowedTimeStop);
    topBranch.GetChild("BreakLoopAtMaxDateTime").GetData(fBreakLoopAtMaxDateTime);
    topBranch.GetChild("MinSignalToNoise").GetData(fMinSignalToNoise);
    topBranch.GetChild("MinNumberOfStationsWithPulse").GetData(fMinNumberOfStationsWithPulse);
    topBranch.GetChild("MaxSDZenith").GetData(fMaxSDZenith);
    topBranch.GetChild("MinSDEnergy").GetData(fMinSDEnergy);
    topBranch.GetChild("AllowedTimeAroundEventTime").GetData(fAllowedTimeAroundEventTime);
    topBranch.GetChild("FDInput").GetData(fFDInput);
    topBranch.GetChild("MaxFDZenith").GetData(fMaxSDZenith);

    const auto vectorSelectedEventIds = topBranch.GetChild("SelectedEventIds").Get<vector<int>>();
    fSelectedEventIds = set<int>(vectorSelectedEventIds.begin(), vectorSelectedEventIds.end());

    const auto vectorDeselectedEventIds = topBranch.GetChild("DeselectedEventIds").Get<vector<int>>();
    fDeselectedEventIds = set<int>(vectorDeselectedEventIds.begin(), vectorDeselectedEventIds.end());

    const auto vectorSelectedEventDates = topBranch.GetChild("SelectedEventDates").Get<vector<TimeStamp>>();
    fSelectedEventDates = set<TimeStamp>(vectorSelectedEventDates.begin(), vectorSelectedEventDates.end());

    ostringstream info;

    const auto vectorSelectedRunEventIds = topBranch.GetChild("SelectedRunEventIds").Get<vector<string>>();
    if (vectorSelectedRunEventIds.size() > 0) {
      info.str("");
      info << "only the following events will be selected:" << std::endl;
      for (string runEventId : vectorSelectedRunEventIds) {
        vector<string> splits;
        boost::split(splits, runEventId, boost::is_any_of("."));
        info << splits[0] << "." << splits[1] << ", ";
        fSelectedRunEventIds.emplace(boost::lexical_cast<int>(splits[0]),
                                     boost::lexical_cast<int>(splits[1]));
      }
      INFOFinal(info);
    }

    const auto vectorDeselectedRunEventIds = topBranch.GetChild("DeselectedRunEventIds").Get<vector<string>>();
    if (vectorDeselectedRunEventIds.size() > 0) {
      info.str("");
      info << "the following events will be rejected:" << std::endl;
      for (string runEventId : vectorDeselectedRunEventIds) {
        vector<string> splits;
        boost::split(splits, runEventId, boost::is_any_of("."));
        info << splits[0] << "." << splits[1] << ", ";
        fDeselectedRunEventIds.emplace(boost::lexical_cast<int>(splits[0]),
                                       boost::lexical_cast<int>(splits[1]));
      }
      INFOFinal(info);
    }

    string tmp = topBranch.GetChild("UsedEye").Get<string>();
    if (tmp == "CO")
      fWhichEye = 4;
    if (tmp == "HE")
      fWhichEye = 5;
    if (tmp == "HECO")
      fWhichEye = 6;

    info.str("");
    info <<" \nEvent id Minimum: "<< fMinEventId << " / Maximum: " << fMaxEventId << "\n";
    info <<"Number of station Minimum: " << fMinNumberOfStations << " / Maximum: " << fMaxNumberOfStations << "\n";
    if (!fUseTriggerInformation)
      info << "Will not use trigger information";
    else {
      info << "Will use the following trigger mask:";
      if (fUseSelfTriggeredEvent)
        info << "\n\tSelf Trigger" ;
      if (fUseExternallyTriggeredEvent)
        info << "\n\tExternal Trigger" ;
      if (fUseCalibrationTriggeredEvent)
        info << "\n\tCalibration Trigger" ;
      if (fUseScintillatorTriggeredEvent)
        info << "\n\tScintillator Trigger" ;
      if (fUseSDTriggeredEvent)
        info << "\n\tSD Trigger" ;
      if (fUseGUITriggeredEvent)
        info << "\n\tGUI Trigger" ;
      if (fUseFDTriggeredEvent)
        info << "\n\tFD Trigger" ;
      if (fUseHEATTriggeredEvent)
        info << "\n\tHEAT Trigger" ;
      if (fUseAERAletTriggeredEvent)
        info << "\n\tAERALET Trigger" ;
      if (fUseAIRPLANETriggeredEvent)
        info << "\n\tAIRPLANE Trigger" ;
      if (fUsePeriodicTriggeredEvent)
        info << "\n\tPeriodic Trigger" ;
      if (fUsePassThroughTriggeredEvent)
        info << "\n\tPassThrough Trigger" ;
    }
    INFOFinal(info);

    string eventTimesFilename = topBranch.GetChild("EventTimesFilename").Get<string>();
    if (eventTimesFilename != "") {
      fstream fin;
      char line[256];
      fin.open(eventTimesFilename.c_str(), std::ios::in);  // opening data
      fin.getline(line, 256); // read/skip first line

      info.str("");
      info << "loop through file " << eventTimesFilename;
      INFODebug(info);

      if (!fin.good()) {
        WARNING("File could not be opened");
        return eFailure;
      }

      while(fin.good()) {
        unsigned long runId, eventId, GPSSecond, GPSNanoSecond;
        fin >> runId >> eventId >> GPSSecond >> GPSNanoSecond;
        fEventGPSSeconds.push_back(GPSSecond);
        fin.getline(line, 256); // read rest of the lines
        info.str("");
        info << runId << "\t" << eventId << "\t" << GPSSecond << "\t"
             << GPSNanoSecond;
        INFODebug(info);
      }

      fin.close();
      std::sort(fEventGPSSeconds.begin(), fEventGPSSeconds.end());
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdEventPreSelector::Run(evt::Event& event)
  {
    ++fNumberOfEvents;

    // Check if there is an radio event at all
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      ++fNRejected;
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();
    ostringstream info;

    // select only event around specific times
    if (fEventGPSSeconds.size() > 0) {
      const unsigned long currentGPSSecond = rEvent.GetHeader().GetTime().GetGPSSecond();
      // find entry that is closest to current event time
      unsigned long earlierEventTime = 0;
      unsigned long laterEventTime = 10000000000;

      const auto low = std::lower_bound(fEventGPSSeconds.begin(), fEventGPSSeconds.end(), currentGPSSecond);
      const auto up = std::upper_bound(fEventGPSSeconds.begin(), fEventGPSSeconds.end(), currentGPSSecond);

      if (low != fEventGPSSeconds.end()) {
        earlierEventTime = fEventGPSSeconds[(low - fEventGPSSeconds.begin()) - 1];
      }

      if (up != fEventGPSSeconds.end()) {
        laterEventTime = *up;;
      }

      if ( !(((currentGPSSecond - earlierEventTime) * utl::second < fAllowedTimeAroundEventTime) ||
             ((laterEventTime - currentGPSSecond) * utl::second < fAllowedTimeAroundEventTime)) ) {
        info.str("");
        info << "event is larger then " << fAllowedTimeAroundEventTime / utl::minute
             << " min away from any given event time" << " (" << currentGPSSecond << " , "
             << earlierEventTime << " - " << laterEventTime << ")";
        INFODebug(info);
        return eContinueLoop;
      }
    } // end select only event around specific times

    // check if there is an SD event and apply max zenith cut
    if (event.HasRecShower() && event.GetRecShower().HasSRecShower()) {
      const det::Detector& detector = det::Detector::GetInstance();
      const utl::CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();
      const utl::Vector axis = event.GetRecShower().GetSRecShower().GetAxis();
      const double zenith = axis.GetTheta(referenceCS);
      const double energy = event.GetRecShower().GetSRecShower().GetEnergy();

      if (zenith > fMaxSDZenith) {
        info.str("");
        info << "SD reconstructed zenith angle is greater than " << fMaxSDZenith / utl::deg
             << "deg, event will be skipped";
        INFOFinal(info.str());

        fNRejected++;
        return eContinueLoop;
      }

      if (zenith < fMinSDZenith) {
        info.str("");
        info << "SD reconstructed zenith angle is smaller than " << fMinSDZenith / utl::deg
             << "deg, event will be skipped";
        INFOFinal(info.str());

        fNRejected++;
        return eContinueLoop;
      }

      if (energy < fMinSDEnergy){
          info.str("");
          info << "SD reconstructed energy is smaller than " <<fMinSDEnergy /utl::eV
               << "eV, event will be skipped";
          INFOFinal(info.str());

        fNRejected++;
        return eContinueLoop;
      }
    }

    // if FD is selected as input: check if there is an FD event, apply max zenith cut and test selected eye
    if (fFDInput) {
      if (!event.HasFEvent()) {
        ERROR("No FEvent available but FD input specified in the settings!");
        return eContinueLoop;
      }

      FEvent& fdEvent = event.GetFEvent();
      for (auto eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
           eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIter) {
        if (eyeIter->GetId() != fWhichEye)
          continue; // continue until selected eye
        const fevt::Eye& eye = *eyeIter;

        if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
          ERROR("Selected FD eye has no RecData or FRecShower but FD input specified in the settings!");
          return eContinueLoop;
        }

        const ShowerFRecData& shower = eye.GetRecData().GetFRecShower();
        const utl::Vector axis = shower.GetAxis();
        const det::Detector& detector = det::Detector::GetInstance();
        const utl::CoordinateSystemPtr referenceCS = detector.GetReferenceCoordinateSystem();
        const double zenith = axis.GetTheta(referenceCS);
        if (zenith > fMaxFDZenith) {
          info.str("");
          info << "FD reconstructed zenith angle is greater than " << fMaxFDZenith / utl::deg
               << "deg, event will be skipped";
          INFOFinal(info);
          fNRejected++;
          return eContinueLoop;
        }
      }
    }

    // check preselected time range
    const TimeStamp& eventts = rEvent.GetHeader().GetTime();
    if (eventts < fAllowedTimeBegin) {
      info.str("");
      info << "Event ID " << rEvent.GetHeader().GetId()
           << " comes before minimum allowed Time";
      INFOFinal(info);

      fNRejected++;
      return eContinueLoop;
    }

    if (eventts > fAllowedTimeStop) {
      info.str("");
      info << "Event ID " << rEvent.GetHeader().GetId()
           << " comes after maximum allowed Time";
      INFOFinal(info.str());

      fNRejected++;
      if (fBreakLoopAtMaxDateTime)
        return eBreakLoop;

      return eContinueLoop;
    }

    // check for specific (de)selected event ids
    if (fSelectedEventIds.size() > 0 &&
        fSelectedEventIds.find(rEvent.GetHeader().GetId()) == fSelectedEventIds.end()) {
      info.str("");
      info << "The current event ID " << rEvent.GetHeader().GetId()
           << " is not selected";
      INFOFinal(info);
      ++fNRejected;
      return eContinueLoop;
    }

    if (fDeselectedEventIds.size() > 0 &&
        fDeselectedEventIds.find(rEvent.GetHeader().GetId()) != fDeselectedEventIds.end()) {
      info.str("");
      info << "The current event ID " << rEvent.GetHeader().GetId()
           << " is deselected";
      INFOFinal(info);
      ++fNRejectedEventRejectionList;
      ++fNRejected;
      return eContinueLoop;
    }

    if (fSelectedEventDates.size() > 0 &&
        fSelectedEventDates.find(rEvent.GetHeader().GetTime()) == fSelectedEventDates.end()) {
      info.str("");
      info << "The current event date " << rEvent.GetHeader().GetTime()
           << " is not selected";
      INFOFinal(info.str());
      ++fNRejected;
      return eContinueLoop;
    }

    const pair<int, int> eventIdentifier (rEvent.GetHeader().GetRunNumber(), rEvent.GetHeader().GetId());
    if (fSelectedRunEventIds.size() > 0 &&
        fSelectedRunEventIds.find(eventIdentifier) == fSelectedRunEventIds.end()) {
      info.str("");
      info << "The current event " << eventIdentifier.first << "." << eventIdentifier.second
           << " is not selected";
      INFOFinal(info);
      ++fNRejected;
      return eContinueLoop;
    }

    if (fDeselectedRunEventIds.size() > 0 &&
        fDeselectedRunEventIds.find(eventIdentifier) != fDeselectedRunEventIds.end()) {
      info.str("");
      info << "The current event ID " << eventIdentifier.first << "." << eventIdentifier.second
           << " is deselected";
      INFOFinal(info);
      ++fNRejectedEventRejectionList;
      ++fNRejected;
      return eContinueLoop;
    }

    // check trigger flags
    if (fUseTriggerInformation && !CheckTrigger(rEvent.GetTrigger())) {
      INFOFinal("Event does not fulfill selected trigger conditions.");
      ++fNRejectedTriggerConditions;
      ++fNRejected;
      return eContinueLoop;
    }

    // check if number of stations contained in the event passes the cuts
    if (rEvent.GetNumberOfStations() < fMinNumberOfStations) {
      info.str("");
      info << "Event must contain at least " << fMinNumberOfStations
           << " stations but does contain only " << rEvent.GetNumberOfStations()
           << " stations.";
      INFOFinal(info);
      ++fNRejectedMinNumberOfStations;
      ++fNRejected;

      if (rEvent.GetNumberOfStations() > fMaxNumberOfAvailableStationsInLargestEvents)
        fMaxNumberOfAvailableStationsInLargestEvents = rEvent.GetNumberOfStations();
      return eContinueLoop;
    }

    if (fMaxNumberOfStations > 0 && rEvent.GetNumberOfStations() > fMaxNumberOfStations) {
      info.str("");
      info << "Event must contain not more than " << fMaxNumberOfStations
           << " stations but does contain " << rEvent.GetNumberOfStations()
           << " stations.";
      INFOFinal(info);
      ++fNRejected;
      return eContinueLoop;
    }

    // check if cuts for event IDs are passed
    if (rEvent.GetHeader().GetId() < fMinEventId) {
      info.str("");
      info << "Event ID must be at least " << fMinEventId << " but it is "
           << rEvent.GetHeader().GetId() << ".";
      INFOFinal(info);
      ++fNRejected;
      return eContinueLoop;
    }

    if (fMaxEventId > 0 && rEvent.GetHeader().GetId() > fMaxEventId) {
      info.str("");
      info << "Event ID must be at most " << fMaxEventId << " but it is "
           << rEvent.GetHeader().GetId() << ".";
      INFOFinal(info);
      ++fNRejected;
      return eContinueLoop;
    }

    // look if all required station ids take part
    for (int stationId : fAllParticipatingStationIds){
      if (!rEvent.HasStation(stationId)) {
        info.str("");
        info << "Required station with ID " << stationId
             << " does not take part in this event.";
        INFOFinal(info);
        ++fNRejected;
        return eContinueLoop;
      }
    }

    // look if one of the required stations take part
    if (fAtLeastOneParticipatingStationIds.size() > 0 ) {
      bool OneTakesPlace = false;
      for (int stationId : fAtLeastOneParticipatingStationIds) {
        if (rEvent.HasStation(stationId)) {
          OneTakesPlace = true;
          break;
        }
      }

      if (!OneTakesPlace) {
        INFOFinal("None of the required station take part in this event.");
        ++fNRejected;
        return eContinueLoop;
      }
    }

    // look if any not-allowed station ids take part
    for (int stationId : fNotParticipatingStationIds)
      if (rEvent.HasStation(stationId)) {
        info.str("");
        info << "Not allowed station with ID " << stationId
             << " does take part in this event.";
        INFOFinal(info);
        ++fNRejectedStationRejectionList;
        ++fNRejected;
        return eContinueLoop;
      }

    // check if a minimum number of stations fulfills the SNR condition (has pulses in the ADC trace)
    if (fMinSignalToNoise > 0){
      const unsigned int numberOfStationsWithPulseFound = CountStationWithPulse(event);
      if (numberOfStationsWithPulseFound < fMinNumberOfStationsWithPulse) {
        info.str("");
        info << "Event must contain at least " << fMinNumberOfStationsWithPulse
             << " stations with ADC pulse fulfilling the condition on an SNR cut of "
             << fMinSignalToNoise << ", but it only contains: "
             << numberOfStationsWithPulseFound << " stations.";
        INFOFinal(info);
        ++fNRejectedMinNumberOfStationsWithPulse;
        ++fNRejected;
        return eContinueLoop;
      }
    }

    // if all cuts are passed, return eSuccess at the end
    ++fNAccepted;
    return eSuccess;
  }


  VModule::ResultFlag
  RdEventPreSelector::Finish()
  {
    ostringstream sstr;
    sstr << fNRejected << " events rejected (total number of events = " << fNumberOfEvents << ")"
         << "\n\tEvent does not fulfill trigger conditions:\t" << fNRejectedTriggerConditions
         << "\n\tEvent in EventRejectionList:\t" << fNRejectedEventRejectionList
         << "\n\tEvent has not enough stations:\t" << fNRejectedMinNumberOfStations
         << " (of those rejected stations, no event had more than "
         << fMaxNumberOfAvailableStationsInLargestEvents << " stations.)"
         << "\n\tEvent has Station in StationRejectionList:\t" << fNRejectedStationRejectionList;
    INFO(sstr);
    return eSuccess;
  }


  unsigned int
  RdEventPreSelector::CountStationWithPulse(evt::Event& event)
    const
  {
    const det::Detector& zedet = det::Detector::GetInstance();
    const REvent& rEvent = event.GetREvent();

    unsigned int numberOfStationsWithPulseFound = 0;
    bool PulseFound = false;

    // loop over stations and for every station over each channel
    for (const auto& station : rEvent.StationsRange()) {

      // RRecStation event should have been generated by RdEventInitializer
      if (!station.HasRecData()) {
        ERROR("Station has no RecData! Please call RdEventInitializer first!");
      }
      const auto& recStation = station.GetRecData();

      // loop over the channels
      for (const auto& channel : station.ChannelsRange()) {
        if (!channel.IsActive())
          continue;

        const auto& ChannelDet =
          zedet.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId());

        // read the trace of the channel
        const auto& trace = channel.GetChannelTimeSeries();

        // check if the content of each trace is correct
        const double bitDepth = ChannelDet.GetBitDepth();
        const double samplingFrequency = ChannelDet.GetSamplingFrequency();
        const double traceBinning = trace.GetBinning();
        const double traceMax = TraceAlgorithm::Max(trace, 0, trace.GetSize()-1);
        const double traceMin = TraceAlgorithm::Min(trace, 0, trace.GetSize()-1);

        if (!CloseTo(1.0 / traceBinning, samplingFrequency, 1e-6)) {
          ostringstream warning;
          warning << "ADC Sampling Frequency of " << samplingFrequency / MHz
                  << " MHz does not match sampling frequency of data, which is "
                  << 1.0/traceBinning / MHz << " MHz. Aborting."
                  << " Was processing Station : " << station.GetId() << " Channel "
                  << channel.GetId();
          ERROR(warning);
        }

        if ((traceMin < 0) || (traceMax > pow(2, bitDepth)-1)) {
          ostringstream warning;
          warning << "ADC trace of Channel " << channel.GetId()
                  << ", Station " << station.GetId()
                  << " has illegal entries. Aborting.";
          ERROR(warning);
          continue;
        }

        if (!traceMin || traceMax == pow(2, bitDepth)-1) {
          ostringstream warning;
          warning << "ADC trace of Channel " << channel.GetId()
                  << ", Station " << station.GetId()
                  << " has saturated samples!";
          INFOIntermediate(warning);
        }

        double NoiseWindowStart = recStation.GetParameter(eNoiseWindowStart);
        double NoiseWindowStop = recStation.GetParameter(eNoiseWindowStop);
        double SignalSearchWindowStart = recStation.GetParameter(eSignalSearchWindowStart);
        double SignalSearchWindowStop = recStation.GetParameter(eSignalSearchWindowStop);

        double MeanNoise = 0;
        double RMSNoise = 0;

        unsigned int sample = 0;

        double PeakAmplitude = 0;
        double PeakTime = 0;
        double PeakTimeError = 0;

        double TraceStop = (trace.GetSize() - 1) * trace.GetBinning();

        // clip windows at trace boundaries
        NoiseWindowStart = min(max(0.0, NoiseWindowStart), TraceStop);
        NoiseWindowStop = max(min(NoiseWindowStop, TraceStop), NoiseWindowStart);

        SignalSearchWindowStart = min(max(0.0, SignalSearchWindowStart), TraceStop);
        SignalSearchWindowStop = max(min(SignalSearchWindowStop, TraceStop),
              SignalSearchWindowStart);

        // calculate noise and find pulse
        utl::RadioTraceUtilities::Noisefinder(trace, RMSNoise, MeanNoise, NoiseWindowStart, NoiseWindowStop);
        utl::RadioTraceUtilities::Pulsefinder(trace, PeakAmplitude, PeakTime, PeakTimeError,
                                              SignalSearchWindowStart, SignalSearchWindowStop, sample);

        PeakAmplitude = PeakAmplitude - MeanNoise;
        RMSNoise = sqrt(Sqr(RMSNoise) - Sqr(MeanNoise));

        const double SignalToNoise = PeakAmplitude / RMSNoise;

        if (SignalToNoise > fMinSignalToNoise)
          PulseFound = true;

      } //end of loop over the channels

    if (PulseFound)
      ++numberOfStationsWithPulseFound;

    } //end of loop over stations

    return numberOfStationsWithPulseFound;
  }


  bool
  RdEventPreSelector::CheckTrigger(const EventTrigger& evtrig)
  {
    INFOIntermediate("RdEventPreSelector::CheckTrigger");
    if (evtrig.IsSelfTrigger())
      INFODebug("IS SELF TRIGGER");
    if (evtrig.IsExternalTrigger())
      INFODebug("IS EXTERNAL TRIGGER");
    if (evtrig.IsCalibrationTrigger())
      INFODebug("IS CALIBRATION TRIGGER");
    if (evtrig.IsScintillatorTrigger())
      INFODebug("IS SCINTILLATOR TRIGGER");
    if (evtrig.IsSDTrigger())
      INFODebug("IS SD TRIGGER");
    if (evtrig.IsGUITrigger())
      INFODebug("IS GUI TRIGGER");
    if (evtrig.IsFDTrigger())
      INFODebug("IS FD TRIGGER");
    if (evtrig.IsHEATTrigger())
      INFODebug("IS HEAT TRIGGER");
    if (evtrig.IsAERAletTrigger())
      INFODebug("IS AERALET TRIGGER");
    if (evtrig.IsAIRPLANETrigger())
      INFODebug("IS AIRPLANE TRIGGER");
    if (evtrig.IsPeriodicTrigger())
      INFODebug("IS PERIODIC TRIGGER");
    if (evtrig.IsPassThroughTrigger())
      INFODebug("IS PASSTHROUGH TRIGGER");

    if (fUseSelfTriggeredEvent && evtrig.IsSelfTrigger())
      return true;
    if (fUseExternallyTriggeredEvent && evtrig.IsExternalTrigger())
      return true;
    if (fUseCalibrationTriggeredEvent && evtrig.IsCalibrationTrigger())
      return true;
    if (fUseScintillatorTriggeredEvent && evtrig.IsScintillatorTrigger())
      return true;
    if (fUseSDTriggeredEvent && evtrig.IsSDTrigger())
      return true;
    if (fUseGUITriggeredEvent && evtrig.IsGUITrigger())
      return true;
    if (fUseFDTriggeredEvent && evtrig.IsFDTrigger())
      return true;
    if (fUseHEATTriggeredEvent && evtrig.IsHEATTrigger())
      return true;
    if (fUseAERAletTriggeredEvent && evtrig.IsAERAletTrigger())
      return true;
    if (fUseAIRPLANETriggeredEvent && evtrig.IsAIRPLANETrigger())
      return true;
    if (fUsePeriodicTriggeredEvent && evtrig.IsPeriodicTrigger())
      return true;
    if (fUsePassThroughTriggeredEvent && evtrig.IsPassThroughTrigger())
      return true;

    return false;
  }


  bool
  RdEventPreSelector::CloseTo(const double a, const double b, const double tolerance)
    const
  {
    return fabs(a - b) < tolerance;
  }

}