RdChannelRiseTimeCalculator.cc

Go to the documentation of this file.

#include "RdChannelRiseTimeCalculator.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/Math.h>

#include <evt/Event.h>

#include <revt/REvent.h>
#include <revt/StationRecData.h>
#include <revt/Channel.h>
#include <revt/Station.h>

#include <vector>

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


namespace RdChannelRiseTimeCalculator {

  VModule::ResultFlag
  RdChannelRiseTimeCalculator::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelRiseTimeCalculator");
    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
  RdChannelRiseTimeCalculator::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;

      const double signalSearchWindowStart = station.GetRecData().GetParameter(eSignalSearchWindowStart);
      const double signalSearchWindowStop = station.GetRecData().GetParameter(eSignalSearchWindowStop);
      const double noiseWindowStart = station.GetRecData().GetParameter(eNoiseWindowStart);
      const double noiseWindowStop = station.GetRecData().GetParameter(eNoiseWindowStop);

      for (Station::ChannelIterator cIt = station.ChannelsBegin();
           cIt != station.ChannelsEnd(); ++cIt) {

        Channel& channel = *cIt;

        if (!channel.IsActive())
          continue;

        // Work on the time series of this channel
        ChannelTimeSeries& trace = channel.GetChannelTimeSeries();

        double noiseRMS = 0;
        Noisefinder(trace, noiseRMS, noiseWindowStart, noiseWindowStop);
        
        // Find the index of the pulse maximum position
        double peakAmplitude = 0;
        double peakTime = 0;
        double peakTimeError = 0;
        unsigned int maxPosition = 0;
        Pulsefinder(trace, peakAmplitude, peakTime, peakTimeError, signalSearchWindowStart, signalSearchWindowStop, maxPosition);
        
        if (!channel.HasRecData())
          channel.MakeRecData();
        channel.GetRecData().SetParameter(eChannelSignalToNoise, peakAmplitude / noiseRMS);

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

        // Construction of the cumulative function
        double cumulativeSum = 0;
        // unsigned int cumulativeVectorSize = riseTimeWindowFrequency;
        std::vector<double> cumulative;
        MakeCumulative(trace, startRiseTimeWindow, stopRiseTimeWindow, cumulative, cumulativeSum);

        // Normalize cumulative function
        std::vector<double> normalizedCumulative;
        for (const auto &value : cumulative)
          normalizedCumulative.push_back(value / cumulativeSum);

        // Rise time calculation
        double riseTime = 0;
        GetRiseTime(startRiseTimeWindow, stopRiseTimeWindow, normalizedCumulative, riseTime, samplingFrequency);
        channel.GetRecData().SetParameter(eRiseTime, riseTime);
      }
    }

    return eSuccess;
  }


  void
  RdChannelRiseTimeCalculator::Noisefinder(const revt::ChannelTimeSeries& channeltrace,
                                           double& rmsNoise, const double noiseWindowStart, const double noiseWindowStop)
     const
  {
    const unsigned int start = noiseWindowStart / channeltrace.GetBinning();
    const unsigned int stop = noiseWindowStop / channeltrace.GetBinning();
    rmsNoise = TraceAlgorithm::RootMeanSquare(channeltrace, start, stop);
  }


  void
  RdChannelRiseTimeCalculator::Pulsefinder(const revt::ChannelTimeSeries& channeltrace,
                                           double& peakAmplitude, double& peakTime, double& peakTimeError,
                                           const double signalSearchWindowStart, const double signalSearchWindowStop,
                                           unsigned int& sample)
    const
  {
    peakAmplitude = 0;
    unsigned int start = signalSearchWindowStart / channeltrace.GetBinning();
    unsigned int stop = signalSearchWindowStop / channeltrace.GetBinning();
    for (unsigned int i = start; i < stop; ++i) {
      if (channeltrace.At(i) > peakAmplitude) {
        peakAmplitude = channeltrace.At(i);
        peakTime = i * channeltrace.GetBinning();
        sample = i;
      }
    }
    //#warning Error of PeakTime is currently set to bin size
    peakTimeError = channeltrace.GetBinning(); // Binsize as error, has to be changed later
  }

  void
  RdChannelRiseTimeCalculator::MakeCumulative(const revt::ChannelTimeSeries& channeltrace,
                                              const unsigned int startRiseTimeWindow, const unsigned int stopRiseTimeWindow,
                                              std::vector<double>& cumulative, double &cumulativeSum)
      const
  {
    for (unsigned int k = startRiseTimeWindow; k < stopRiseTimeWindow + 1; ++k) {
      cumulativeSum += Sqr(channeltrace.At(k));
      cumulative.push_back(cumulativeSum);
    }
  }


  void
  RdChannelRiseTimeCalculator::GetRiseTime(const unsigned int startCumulative, const unsigned int stopCumulative,
                                           std::vector<double>& normalizedCumulative, double &riseTime,
                                           const double samplingFrequency)
      const
  {
    bool foundmin = false;
    bool foundmax = false;
    unsigned int minCumulative = 0;
    unsigned int maxCumulative = 0;
    double lowerCumulativeLimit = 0.2;
    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;
    }
    riseTime = (maxCumulative - minCumulative) / samplingFrequency;
  }

}