RdStationRiseTimeCalculator.cc

Go to the documentation of this file.

#include "RdStationRiseTimeCalculator.h"

#include <fwk/CentralConfig.h>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Trace.h>
#include <utl/Math.h>

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

#include <rdet/RDetector.h>
#include <rdet/Station.h>

using namespace revt;
using namespace fwk;
using namespace utl;
using namespace std;


namespace RdStationRiseTimeCalculator {

  VModule::ResultFlag RdStationRiseTimeCalculator::Init()
  {
    // Initialize your module here. This method
    // is called once at the beginning of the run.
    // The eSuccess flag indicates the method ended
    // successfully.  For other possible return types,
    // see the VModule documentation.

    INFO("RdStationRiseTimeCalculator::Init()");

    // Read in the configurations of the xml file
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationRiseTimeCalculator");
    topBranch.GetChild("RiseTimeWindow").GetData(fRiseTimeWindow);

    if (fRiseTimeWindow > 2000*nanosecond) {
      ERROR("this window should be narrower: short window around the main pulse");
      return eFailure;
    }

    return eSuccess;
  }


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

    REvent& rEvent = event.GetREvent();

    // loop through stations and for every station through every channel
    for (REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
         sIt != rEvent.CandidateStationsEnd(); ++sIt) {

       Station& station = *sIt;

      if (!station.HasRecData()) {
        WARNING("Station has no RecData! Please call RdEventInitializer first!");
        return eFailure;
      }

      //double SignalSearchWindowStart = station.GetRecData().GetParameter(eSignalSearchWindowStart);
      //double SignalSearchWindowStop = station.GetRecData().GetParameter(eSignalSearchWindowStop);

      // get station component with max SNR
      const double SNR[3] = {
        station.GetRecData().GetParameter(ePeakAmplitudeNS)/station.GetRecData().GetParameter(eNoiseRmsNS),
        station.GetRecData().GetParameter(ePeakAmplitudeEW)/station.GetRecData().GetParameter(eNoiseRmsEW),
        station.GetRecData().GetParameter(ePeakAmplitudeV)/station.GetRecData().GetParameter(eNoiseRmsV)
      };

      int chwithmaxSNR = 0;
      double tempSNRmax = 0;
      for (int k = 0; k < 3; ++k) {
        if (SNR[k] > tempSNRmax) {
          tempSNRmax = SNR[k];
          chwithmaxSNR = k;
        }
      }

      const double relTime = station.GetRecData().GetParameter(eTraceStartTime);

      // get trace with max SNR
      const StationTimeSeries& stationtrace = station.GetStationTimeSeries();

      unsigned int maxPosition = 0;
      if (chwithmaxSNR == 0)
        maxPosition = station.GetRecData().GetParameter(ePeakTimeNS) - relTime;
      else if (chwithmaxSNR == 1)
        maxPosition = station.GetRecData().GetParameter(ePeakTimeEW) - relTime;
      else if (chwithmaxSNR == 2)
        maxPosition = station.GetRecData().GetParameter(ePeakTimeV) - relTime;

      const double samplingFrequency = 1. / stationtrace.GetBinning();

      // define the window used for rise time calculation
      unsigned int startRiseTimeWindow = 0;
      unsigned int stopRiseTimeWindow = 0;
      if (maxPosition) {
        startRiseTimeWindow = (maxPosition - (fRiseTimeWindow/2.)) * samplingFrequency;
        stopRiseTimeWindow = (maxPosition + (fRiseTimeWindow/2.)) * samplingFrequency;
        // work-around when startRiseTimeWindow moves out of the trace
        //#warning JN: this check should probably be done somewhere else
        if (maxPosition < fRiseTimeWindow/2.) {
          WARNING("RiseTimeWindowStart outside of trace, RiseTimeWindow adjusted to start of trace");
          startRiseTimeWindow = 0;
          stopRiseTimeWindow = fRiseTimeWindow * samplingFrequency;
        }
      }

      // construction of the cumulative function
      double cumulativeSum = 0;
      const unsigned int cumulativeVectorSize = static_cast<unsigned int>(fRiseTimeWindow * samplingFrequency);
      vector<double> cumulative(cumulativeVectorSize);
      MakeCumulative(stationtrace, startRiseTimeWindow, stopRiseTimeWindow, cumulative, cumulativeSum, chwithmaxSNR);

      // normalize cumulative function
      vector<double> normalizedCumulative(cumulativeVectorSize);
      for (unsigned int k = startRiseTimeWindow; k < stopRiseTimeWindow + 1; ++k)
        normalizedCumulative.at(k - startRiseTimeWindow) = cumulative.at(k - startRiseTimeWindow) / cumulativeSum;

      // rise time calculation
      unsigned int minCumulative = 0;
      unsigned int maxCumulative = 0;
      GetRiseTime(startRiseTimeWindow, stopRiseTimeWindow, normalizedCumulative, minCumulative, maxCumulative);
      const double riseTime = (maxCumulative - minCumulative) / samplingFrequency;

      station.GetRecData().SetParameter(eStationRiseTime, riseTime);

    }

    return eSuccess;
  }


  void
  RdStationRiseTimeCalculator::MakeCumulative(const revt::StationTimeSeries& stationtrace,
                                              const unsigned int startRiseTimeWindow, const unsigned int stopRiseTimeWindow,
                                              vector<double>& cumulative, double& cumulativeSum, const unsigned int chwithmaxSNR)
    const
  {
    cumulativeSum = Sqr(stationtrace.At(startRiseTimeWindow)[chwithmaxSNR]);
    cumulative.at(0) = cumulativeSum;
    for (unsigned int k = startRiseTimeWindow + 1; k < stopRiseTimeWindow + 1; ++k) {
      const double tmpSum = Sqr(stationtrace.At(k)[chwithmaxSNR]);
      cumulativeSum += tmpSum;
      cumulative.at(k - startRiseTimeWindow) = cumulative.at(k - startRiseTimeWindow - 1) + tmpSum;
    }
  }


  void
  RdStationRiseTimeCalculator::GetRiseTime(const unsigned int startCumulative, const unsigned int stopCumulative,
                                           vector<double>& normalizedCumulative, unsigned int& minCumulative, unsigned int& maxCumulative)
    const
  {
    bool foundmin = false;
    bool foundmax = false;
    const double lowerCumulativeLimit = 0.2;
    const double upperCumulativeLimit = 0.6;

    for (unsigned int k = startCumulative; k < stopCumulative+1; ++k) {
      if (!foundmin) {
        if (normalizedCumulative.at(k - startCumulative) < lowerCumulativeLimit)
          continue;
        else {
          minCumulative = k;
          foundmin = true;
        }
      }
      if (!foundmax) {
        if (normalizedCumulative.at(k - startCumulative) < upperCumulativeLimit)
          continue;
        else {
          maxCumulative = k;
          foundmax = true;
        }
      }
      if (foundmin && foundmax)
        break;
    }
  }

}