RdStationSimPulseFinder.cc

Go to the documentation of this file.

#include "RdStationSimPulseFinder.h"

#include <fwk/CentralConfig.h>
#include <fwk/MagneticFieldModel.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/FFTDataContainer.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/RadioTraceUtilities.h>

#include <evt/Event.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSimData.h>

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

#include <TVector3.h>


using namespace utl;


namespace RdStationSimPulseFinder {

  fwk::VModule::ResultFlag
  RdStationSimPulseFinder::Init()
  {
    Branch topBranch = fwk::CentralConfig::GetInstance()->GetTopBranch("RdStationSimPulseFinder");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("LowerCutoffFrequency").GetData(fLowerCutoffFrequency);
    topBranch.GetChild("UpperCutoffFrequency").GetData(fUpperCutoffFrequency);
    topBranch.GetChild("ApplyBandpassFilter").GetData(fApplyBandpassFilter);
    topBranch.GetChild("UseHilbertEnvelope").GetData(fUseHilbertEnvelope);
    topBranch.GetChild("SignalIntegrationWindowSize").GetData(fSignalIntegrationWindowSize);
    topBranch.GetChild("ReferenceShower").GetData(fReferenceShower);

    std::ostringstream info;
    info << "\n\tApply bandpass filter: " << fApplyBandpassFilter
         << " (" << fLowerCutoffFrequency / MHz << " - " << fUpperCutoffFrequency / MHz << ") MHz"
         << "\n\tUse hilbert envelope: " << fUseHilbertEnvelope
         << "\n\tSize of the signal integration window: " << fSignalIntegrationWindowSize / ns << " ns"
         << "\n\tUse reference direction for coorindate transformation from: " << fReferenceShower << "\n";
    INFOFinal(info);

    return eSuccess;
  }


  fwk::VModule::ResultFlag
  RdStationSimPulseFinder::Run(evt::Event& event)
  {
    auto& rEvent = event.GetREvent();
    const TimeStamp& rEventTime = rEvent.GetHeader().GetTime();

    std::ostringstream info;

    for (auto& station : rEvent.StationsRange()) {
      info.str("");
      info << "Getting StationRecData for station " << station.GetId();
      INFODebug(info);

      // create has sim data when not done yet
      if (!station.HasSimData()) {
        station.MakeSimData();
      }

      auto& stationSimData = station.GetSimData();
      auto efieldData = station.GetFFTDataContainer();
      auto& spectrum = efieldData.GetFrequencySpectrum();

      // apply the bandpass filter
      if (fApplyBandpassFilter) {
        for (unsigned int i = 0; i < spectrum.GetSize(); ++i) {
          if ((efieldData.GetFrequencyOfBin(i) < fLowerCutoffFrequency) ||
              (efieldData.GetFrequencyOfBin(i) > fUpperCutoffFrequency)) {
            spectrum[i][0] = std::complex<double>(0., 0.);
            spectrum[i][1] = std::complex<double>(0., 0.);
            spectrum[i][2] = std::complex<double>(0., 0.);
          }
        }
      }

      // get trace without envelope
      const revt::StationTimeSeries& stationTrace = efieldData.GetConstTimeSeries();

      // apply envelope
      revt::StationTimeSeries stationTraceForPeakFinding;
      if (fUseHilbertEnvelope) {
        INFODebug("Applying Hilbert envelope");

        auto efieldDataEnvelope = efieldData;
        FFTDataContainerAlgorithm::HilbertEnvelope(efieldDataEnvelope);
        stationTraceForPeakFinding = efieldDataEnvelope.GetConstTimeSeries();
      }
      else {
        stationTraceForPeakFinding = stationTrace;
      }

      // searching for pulse
      unsigned int peakInt = 0;
      double peakAmplitude = 0;
      double peakTime = 0;
      double peakTimeError = 0;
      const double windowstart = 0;
      const double windowstop = stationTraceForPeakFinding.GetStop() * stationTraceForPeakFinding.GetBinning();
      INFODebug("Finding pulse ...");
      RadioTraceUtilities::Pulsefinder(stationTraceForPeakFinding, peakAmplitude, peakTime, peakTimeError,
                                       windowstart, windowstop, peakInt);

      info.str("");
      info << "Setting SimPulseTime of station " << station.GetId() << " to " << peakTime
           << " (sample " << peakInt << ")" << " and SimPulseAmplitude to " << peakAmplitude;
      INFOIntermediate(info);

      if (!peakAmplitude) {
        info.str("");
        info << "No sim pulse found for station " << station.GetId();
        INFOIntermediate(info);
        continue;
      }

      const TimeInterval traceStartTime = station.GetRawTraceStartTime() - rEventTime;
      stationSimData.SetParameter(revt::eTraceStartTime, traceStartTime);

      stationSimData.SetParameter(revt::eSignal, peakAmplitude);
      stationSimData.SetParameter(revt::eSignalTime, peakTime + traceStartTime);

      // calculating E-field vector
      double maxPeak = 0;
      unsigned int maxPol = 0;
      for (unsigned int j = 0; j <= 2; ++j) {
        for (unsigned int i = 0; i < stationTrace.GetSize(); ++i) {
          if (fabs(stationTrace[i][j]) > maxPeak) {
            maxPeak = fabs(stationTrace[i][j]);
            maxPol = j;
          }
        }
      }

      double tMaxAmp = 0;
      unsigned int tMaxPos = 0;
      double peakSearchWindow = 20. * ns;
      for (unsigned int i = peakInt + peakSearchWindow / stationTrace.GetBinning();
           i > (peakInt - peakSearchWindow / stationTrace.GetBinning()); --i) {
        if (fabs(stationTrace[i][maxPol]) > tMaxAmp) {
          tMaxAmp = fabs(stationTrace[i][maxPol]);
          tMaxPos = i;
        }
      }

      // first we calculated the FWHM
      unsigned int FWHMlowerBound = 0;
      unsigned int FWHMupperBound = 0;
      for (unsigned int i = tMaxPos; i > 0; --i) {
        if (stationTraceForPeakFinding[i].GetMag() < peakAmplitude * .5) {
          FWHMlowerBound = i;
          break;
        }
      }

      if (FWHMlowerBound == 0) {
        INFODebug("FWHMLowerBound is 0. This may mean that something went wrong");
      }

      for (unsigned int i = tMaxPos; i < stationTrace.GetSize(); ++i) {
        if (stationTraceForPeakFinding[i].GetMag() < peakAmplitude * .5) {
          FWHMupperBound = i;
          break;
        }
      }

      if (FWHMupperBound == stationTrace.GetSize()) {
        INFODebug("FWHMupperBound is the last entry in the station trace. This may mean that something went wrong");
      }

      // now we know over which interval to integrate. Continue with calculating the average E-field vector
      TVector3 initialGuess(stationTrace[peakInt][0], stationTrace[peakInt][1], stationTrace[peakInt][2]);
      TVector3 meanEField(0, 0, 0);
      for (unsigned int i = FWHMlowerBound; i < FWHMupperBound; ++i) {
        TVector3 currentEField(stationTrace[i][0], stationTrace[i][1], stationTrace[i][2]);
        if (currentEField.Angle(initialGuess) > 90. * degree) {
          currentEField *= -1;
        }
        meanEField += currentEField;
      }
      meanEField.SetMag(peakAmplitude);

      info.str("");
      info << "Sim E-field vector for this station is (" << meanEField[0] << ","
           << meanEField[1] << "," << meanEField[2] << ").";
      INFODebug(info);

      stationSimData.SetParameter(revt::eEFieldVectorEW, meanEField[0]);
      stationSimData.SetParameter(revt::eEFieldVectorNS, meanEField[1]);
      stationSimData.SetParameter(revt::eEFieldVectorV, meanEField[2]);

      // get reference geometry
      Point coreSite;
      Vector showerAxis;
      if (fReferenceShower == "Reference") {
        const evt::ShowerRRecData& recShower = event.GetRecShower().GetRRecShower();
        coreSite = recShower.GetReferenceCorePosition(event);
        showerAxis = recShower.GetReferenceAxis(event);
      }
      if (fReferenceShower == "SD") {
        const evt::ShowerSRecData& recShower = event.GetRecShower().GetSRecShower();
        coreSite = recShower.GetCorePosition();
        showerAxis = recShower.GetAxis();
      }
      if (fReferenceShower == "FD") {
        const evt::ShowerFRecData& recShower = event.GetRecShower().GetFRecShower();
        coreSite = recShower.GetCorePosition();
        showerAxis = recShower.GetAxis();
      }
      if (fReferenceShower == "MC") {
        const evt::ShowerSimData& recShower = event.GetSimShower();
        coreSite = recShower.GetPosition();
        showerAxis = recShower.GetDirection();
      }

      // init trace container
      revt::ChannelTimeSeries stationChannel;
      stationChannel.SetBinning(stationTrace.GetBinning());

      // transformation in vxB vxvxB system
      const auto coreCS = fwk::LocalCoordinateSystem::Create(coreSite);
      const auto localMagneticField = fwk::MagneticFieldModel::instance().GetMagneticFieldVector(
          coreSite, rEventTime.GetGPSSecond());
      const auto csTrans = RadioGeometryUtilities(showerAxis, coreCS, localMagneticField);
      TraceV3D traceShowerPlane = csTrans.GetTraceInShowerPlaneVxB(stationTrace);

      // loop over polarisations
      double integratedSignal = 0;
      for (unsigned int i = 0; i <= 2; ++i) {
        integratedSignal = 0;
        stationChannel.Clear();

        if (i == 0) {
          for (const auto& element : stationTrace) {
            stationChannel.PushBack(element.GetMag());
          }
        } else {
          for (revt::StationTimeSeries::SizeType j = 0; j < traceShowerPlane.GetSize(); ++j) {
            stationChannel.PushBack(traceShowerPlane[j][i - 1]);
          }
        }

        // calculation of energy fluence
        double integratedNoisePowerFixedWidthStartTime = 0;
        double integratedNoisePowerFixedWidthStopTime = 0;
        RadioTraceUtilities::PulseFixedWindowIntegrator(stationChannel, peakInt, fSignalIntegrationWindowSize,
            integratedSignal, integratedNoisePowerFixedWidthStartTime, integratedNoisePowerFixedWidthStopTime, true);
        const double energyFluence = integratedSignal * kConversionRadioSignalToEnergyFluence;

        if (i == 0) {
          stationSimData.SetParameter(revt::eSignalEnergyFluenceMag, energyFluence);
        } else if (i == 1) {
          stationSimData.SetParameter(revt::eSignalEnergyFluenceVxB, energyFluence);
        } else {
          stationSimData.SetParameter(revt::eSignalEnergyFluenceVxVxB, energyFluence);
        }
      }
    }

    return eSuccess;
  }


  fwk::VModule::ResultFlag
  RdStationSimPulseFinder::Finish()
  {
    return eSuccess;
  }

}