RdStationPulseShapeRejector.cc

Go to the documentation of this file.

#include "RdStationPulseShapeRejector.h"

#include <cmath>

#include <fwk/CentralConfig.h>
#include <fwk/VModule.h>
#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <algorithm>

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

#include <evt/Event.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerRecData.h>
#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>

#include <TMath.h>

using namespace std;

namespace RdStationPulseShapeRejector {

  fwk::VModule::ResultFlag
  RdStationPulseShapeRejector::Init()
  {
    const auto topB = fwk::CentralConfig::GetInstance()->GetTopBranch("RdStationPulseShapeRejector");
    topB.GetChild("InfoLevel").GetData(fInfoLevel);
    topB.GetChild("checkT1Crossings").GetData(fCheckT1Crossings);
    topB.GetChild("checkT2Crossings").GetData(fCheckT2Crossings);
    topB.GetChild("T1").GetData(fT1);
    topB.GetChild("T2").GetData(fT2);
    topB.GetChild("T1MinSNR").GetData(fT1MinSNR);
    topB.GetChild("PulseRegion").GetData(fPulseRegion);
    topB.GetChild("TPer").GetData(fTPer);
    topB.GetChild("T2MinSNR").GetData(fT2MinSNR);
    topB.GetChild("NCMax").GetData(fNCMax);
    topB.GetChild("UseVerticalComponent").GetData(fUseVerticalComponent);
    return eSuccess;
  }

  fwk::VModule::ResultFlag
  RdStationPulseShapeRejector::Run(evt::Event& event)
  {
    stringstream sstr;
    revt::REvent& rEvent = event.GetREvent();
    for (revt::REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin(); sIt != rEvent.SignalStationsEnd();
         ++sIt) {
      fTotalStations++;
      revt::Station& station = *sIt;

      int traceMaxInt = FindTraceMaximum(station);
      bool t1CrossingDeselected = false;

      if (fCheckT1Crossings) {
        t1CrossingDeselected = !CheckT1Crossings(station, traceMaxInt);
        sstr.str("");
        sstr << "T1Crossing check returned " << t1CrossingDeselected;
        INFODebug(sstr.str());
        if (t1CrossingDeselected) {
          sstr.str("");
          sstr << "Rejecting station " << station.GetId();
          INFOFinal(sstr.str());

          station.SetRejectedReason(revt::ePulseShapeDeselected);
          fRejectedStations++;
        }
      }

      if ((!t1CrossingDeselected) && fCheckT2Crossings) {
        sstr.str("");
        sstr << "Checking station " << station.GetId();
        INFODebug(sstr.str());
        CheckT2Crossings(station, traceMaxInt);
      }
    }

    return eSuccess;
  }

  bool
  RdStationPulseShapeRejector::CheckT1Crossings(revt::Station& station, int peakInt)
  {
    stringstream sstr;
    sstr.str("");
    sstr << "RdStationPulseShapeRejector::CheckT1Crossings";
    INFODebug(sstr.str());

    revt::StationFFTDataContainer& efieldData = station.GetFFTDataContainer();
    const revt::StationTimeSeries& stationTrace = efieldData.GetConstTimeSeries();
    const double peakTime = station.GetRecData().GetParameter(revt::ePeakTimeMag) -
                            station.GetRecData().GetParameter(revt::eTraceStartTime);
    const double binning = stationTrace.GetBinning() / utl::nanosecond;
    const double peakMag = stationTrace[peakInt].GetMag();
    double prevMag = 0;
    double signalToNoise = 0;
    double prevSNR = 0;
    const revt::StationRecData& recStation = station.GetRecData();
    const double noise = recStation.GetParameter(revt::eNoise);
    int start = recStation.GetParameter(revt::eSignalSearchWindowStart) / binning;
    int stop = (peakTime - fPulseRegion) / utl::nanosecond / binning;

    for (int i = start; i < stop; i++) {
      signalToNoise = stationTrace[i].GetMag() / noise;
      if ((stationTrace[i].GetMag() > fT1 * peakMag && signalToNoise > sqrt(fT1MinSNR)) &&
          (prevMag < fT1 * peakMag || prevSNR < sqrt(fT1MinSNR))) {
        sstr.str("");
        sstr << "station rejected. peak located at " << peakInt
             << " with magnitude " << peakMag << " . Crossed at " << i
             << " with magnitude " << stationTrace[i].GetMag();
        INFODebug(sstr.str());
        return false;
      }

      prevMag = stationTrace[i].GetMag();
      prevSNR = signalToNoise;
    }

    prevMag = 0;
    prevSNR = 0;
    start = (peakTime + fPulseRegion) / utl::nanosecond / binning;
    stop = recStation.GetParameter(revt::eSignalSearchWindowStop) / binning;

    for (int i = start; i < stop; i++) {
      signalToNoise = stationTrace[i].GetMag() / noise;
      if ((stationTrace[i].GetMag() > fT1 * peakMag && signalToNoise > sqrt(fT1MinSNR)) &&
          (prevMag < fT1 * peakMag || prevSNR < sqrt(fT1MinSNR))) {
        sstr.str("");
        sstr << "station rejected. peak located at " << peakInt << " with magnitude "
             << peakMag << " . Crossed at " << i
             << " with magnitude " << stationTrace[i].GetMag();
        INFODebug(sstr.str());
        return false;
      }

      prevSNR = signalToNoise;
      prevMag = stationTrace[i].GetMag();
    }

    return true;
  }

  void
  RdStationPulseShapeRejector::CheckT2Crossings(revt::Station& station, int peakInt)
  {
    stringstream sstr;
    sstr.str("");
    sstr << "CheckT2Crossings";
    INFODebug(sstr.str());

    revt::StationFFTDataContainer& efieldData = station.GetFFTDataContainer();
    const revt::StationTimeSeries& stationTrace = efieldData.GetConstTimeSeries();
    const double binning = stationTrace.GetBinning();
    const double peakMag = stationTrace[peakInt].GetMag();
    int start = peakInt - (fTPer / binning);
    int stop = peakInt + (fTPer / binning);
    int numberOfCrossings = 0;
    double prevMag = stationTrace[start].GetMag();
    const revt::StationRecData& recStation = station.GetRecData();
    double noise;
    if (fUseVerticalComponent) {
      noise = recStation.GetParameter(revt::eNoiseRmsMag);
    } else {
      noise = recStation.GetParameter(revt::eNoiseRmsNSEWPlane);
    }

    double signalToNoise = 0;
    double prevSNR = 0;
    double mag;
    for (int i = start; i < peakInt; i++) {
      if (fUseVerticalComponent) {
        mag = stationTrace[i].GetMag();
      } else {
        mag = sqrt(stationTrace[i][0] * stationTrace[i][0] + stationTrace[i][1] * stationTrace[i][1]);
      }

      signalToNoise = mag / noise;
      if ((mag > fT2 * peakMag && signalToNoise > sqrt(fT2MinSNR)) && (prevMag <= fT2 * peakMag || prevSNR < sqrt(fT2MinSNR))) {
        sstr.str("");
        sstr << "NC increased at index " << i << " and relative time " << (i - peakInt) * binning;
        INFODebug(sstr.str());
        numberOfCrossings++;
      }

      prevMag = mag;
      prevSNR = signalToNoise;
    }

    if (numberOfCrossings > fNCMax + 1) {    //Plus 1 because the pulse maximum itself will also be counted by NC
      sstr.str("");
      sstr << "Station rejected: " << station.GetId();
      INFOFinal(sstr.str());

      station.SetRejectedReason(revt::ePulseShapeDeselected);
      fRejectedStations++;
      return;
    } else {
      sstr.str("");
      sstr << "Station not rejected with NC=" << numberOfCrossings;
      INFODebug(sstr.str());
    }

    numberOfCrossings = 0;
    for (int i = peakInt; i < stop; i++) {
      if (fUseVerticalComponent) {
        mag = stationTrace[i].GetMag();
      } else {
        mag = sqrt(stationTrace[i][0] * stationTrace[i][0] + stationTrace[i][1] * stationTrace[i][1]);
      }

      signalToNoise = mag / noise;
      if ((mag > fT2 * peakMag && signalToNoise > sqrt(fT2MinSNR)) && (prevMag <= fT2 * peakMag || prevSNR < sqrt(fT2MinSNR))) {
        sstr.str("");
        sstr << "NC increased at index " << i << " and relative time " << (i - peakInt) * binning;
        INFODebug(sstr.str());
        numberOfCrossings++;
      }

      prevSNR = signalToNoise;
      prevMag = mag;
    }

    if (numberOfCrossings > fNCMax) {
      sstr.str("");
      sstr << "Station rejected: " << station.GetId();
      INFOFinal(sstr.str());

      station.SetRejectedReason(revt::ePulseShapeDeselected);
      fRejectedStations++;
      return;
    } else {
      sstr.str("");
      sstr << "Station not rejected with NC=" << numberOfCrossings;
      INFODebug(sstr.str());
    }
  }

  int
  RdStationPulseShapeRejector::FindTraceMaximum(revt::Station& station)
  {
    const revt::StationTimeSeries& stationTrace = station.GetFFTDataContainer().GetConstTimeSeries();
    const double peakTime = station.GetRecData().GetParameter(revt::ePeakTimeMag) -
                            station.GetRecData().GetParameter(revt::eTraceStartTime);
    const double binning = stationTrace.GetBinning();
    int peakInt = peakTime / binning;
    double peakMag = stationTrace[peakInt].GetMag();
    stringstream sstr;
    for (int i = peakInt - fPulseRegion; i < peakInt + fPulseRegion; i++) {
      if (stationTrace[i].GetMag() > peakMag) {
        sstr.str("");
        sstr << "Corrected peak from i=" << peakInt << " and P=" << peakMag
             << " to i=" << i << " t=" << i * binning << " and P=" << stationTrace[i].GetMag();
        INFODebug(sstr.str());

        peakMag = stationTrace[i].GetMag();
        peakInt = i;
      }
    }

    sstr.str("");
    sstr << "Peak is at index " << peakInt << " and time "
         << peakInt * binning + station.GetRecData().GetParameter(revt::eTraceStartTime);
    INFODebug(sstr.str());

    return peakInt;
  }

  fwk::VModule::ResultFlag
  RdStationPulseShapeRejector::Finish()
  {
    stringstream sstr;
    sstr.str("");
    sstr << "Rejected " << fRejectedStations << " stations out of a total of " << fTotalStations << " stations.";
    INFOFinal(sstr.str());
    return eSuccess;
  }

}