RdStationEFieldVectorCalculator.cc

Go to the documentation of this file.

#include "RdStationEFieldVectorCalculator.h"

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

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

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

#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>

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


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


namespace RdStationEFieldVectorCalculator {

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

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("meanEFieldVector").GetData(fMeanEFieldVector);
    topBranch.GetChild("trace").GetData(fTrace);
    topBranch.GetChild("peakSearch").GetData(fPeakSearch);
    topBranch.GetChild("peakSearchWindow").GetData(fPeakSearchWindowWidth);
    topBranch.GetChild("ChargeExcessStrength").GetData(fChargeExcessStrength);

    return eSuccess;
  }


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

    // check if needed reconstructed quantities exist
    if (!event.GetRecShower().HasRRecShower()) {
      WARNING("Event does not have RRecShower, nothing will be done!");
      return eSuccess;
    }

    ostringstream info;

    REvent& rEvent = event.GetREvent();
    auto& rRecShower = event.GetRecShower().GetRRecShower();
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    const Vector showerAxis =rRecShower.GetReferenceAxis(event);
    const Point showerCore = rRecShower.GetReferenceCorePosition(event);
    const CoordinateSystemPtr coordinateSystem = LocalCoordinateSystem::Create(showerCore);
    Vector magneticFieldAxis = rRecShower.GetMagneticFieldVector();

    // use magnetic field from simulations
    if (event.HasSimShower()) {
      INFOFinal("Using magnetic field from simulated air shower.");
      const double magFieldAzimuth = event.GetSimShower().GetMagneticFieldAzimuth();
      const double magFieldZenith = event.GetSimShower().GetMagneticFieldZenith();
      magneticFieldAxis = Vector(1, magFieldZenith, magFieldAzimuth, coordinateSystem,
                                 Vector::kSpherical);
    }

    // loop through stations
    for (auto& station : rEvent.StationsRange()) {
      if (!station.HasRecData()) {
        WARNING("No RecData found");
        return eContinueLoop;
      }

      const Point stationPosition = rDetector.GetStation(station).GetPosition();

      // Read the time trace of this station
      // getting a local copy of the electric field data
      const auto& efieldData = station.GetFFTDataContainer();
      const auto& stationtrace = efieldData.GetConstTimeSeries();
      fPeakSearchWindow = static_cast<int>(fPeakSearchWindowWidth / stationtrace.GetBinning());

      info.str("");
      info << "PeakSearchWindow = " << fPeakSearchWindowWidth << " ns = "
           << fPeakSearchWindow << " bins";
      INFOIntermediate(info);

      StationRecData& recStation = station.GetRecData();

      if (!(station.GetRecData().HasParameter(ePeakAmplitudeEW)
          && station.GetRecData().HasParameter(ePeakAmplitudeNS)
          && station.GetRecData().HasParameter(ePeakAmplitudeV))) {
        WARNING("ePeakAmplitude(EW,NS,V) was not set");
        continue;
      }

      const array<double, 3> peakAmp = {recStation.GetParameter(ePeakAmplitudeEW),
                                        recStation.GetParameter(ePeakAmplitudeNS),
                                        recStation.GetParameter(ePeakAmplitudeV)
                                       };

      // the polarization with the maximum amplitude or the one specified in the xml
      unsigned int maxPol = 0;
      if (fTrace == 3) {
        auto iteratorOfMax = max_element(peakAmp.begin(), peakAmp.end());
        maxPol = distance(peakAmp.begin(), iteratorOfMax);
      }
      else {
        maxPol = fTrace;  // use polarization defined in xml
      }

      // determine peak position of combo envelope
      const double tPeakTimeEnvelope =
        (recStation.GetParameter(eSignalTime) - recStation.GetParameter(eTraceStartTime));
      const long unsigned int tPeakPos = tPeakTimeEnvelope / stationtrace.GetBinning();

      long unsigned int peakPos = 0;

      if (fPeakSearch != 0) {
        info.str("");
        info << "Station " << station.GetId() << ": fPeakSearch";
        INFOIntermediate(info);

        double tMaxAmp = 0;
        unsigned int tMaxPos = 0;

        // search for maximum before the maximum of the combo envelope
        if (fPeakSearch == -1) {
          for (auto i = tPeakPos; i > (tPeakPos - fPeakSearchWindow); --i) {
            if (fabs(stationtrace[i][maxPol]) > tMaxAmp) {
              tMaxAmp = fabs(stationtrace[i][maxPol]);
              tMaxPos = i;
            }
          }
          peakPos = tMaxPos;

          info.str("");
          info << "Station " << station.GetId() << ": searching " << fPeakSearchWindow
               << "ns before the maximum of the combo hilbert envelope. MaxCombo at "
               << tPeakPos << "ns and max of polarization " << maxPol << " at " << peakPos;
          INFOIntermediate(info);
        }

        // search for maximum after the maximum of the combo envelope
        if (fPeakSearch == 1) {
          for (auto i = tPeakPos; i < (tPeakPos + fPeakSearchWindow); ++i) {
            if (fabs(stationtrace[i][maxPol]) > tMaxAmp) {
              tMaxAmp = fabs(stationtrace[i][maxPol]);
              tMaxPos = i;
            }
          }
          peakPos = tMaxPos;

          info.str("");
          info << "Station " << station.GetId() << ": searching " << fPeakSearchWindow
               << "ns after the maximum of the combo hilbert envelope. MaxCombo at "
               << tPeakPos << "ns and max of polarization " << maxPol << " at " << peakPos;
          INFODebug(info);
        }
      }

      info.str("");
      info << "Station " << station.GetId() << ": Found peak position at " << peakPos << " = "
           << peakPos * stationtrace.GetBinning() << "ns in polarization " << maxPol;
      INFOIntermediate(info);

      // set electric field vector to maximum amplitudes of each efield component.
      // This is used as first guess in the efield vector reconstruction (averaging)
      array<double, 3> EVector = {0, 0, 0};
      for (unsigned int iPol = 0; iPol < 3; ++iPol) {
        EVector[iPol] = stationtrace[peakPos][iPol];
      }

      // values are meaningless if not fMeanEFieldVector?
      double angleToLorentzForce = 0;
      double angleToExpectedEField = 0;

      if (fMeanEFieldVector) {
        // Get HilbertEnvelope of Efield Traces
        auto efieldEnvelope = station.GetFFTDataContainer();
        FFTDataContainerAlgorithm::HilbertEnvelope(efieldEnvelope);
        const auto& stationEnvelope = efieldEnvelope.GetConstTimeSeries();

        // determine FWHM of HilbertEnvelope 
        const double maxSignal = recStation.GetParameter(eSignal);
        const pair<double, double> fwhm = GetFWHM(stationEnvelope, maxSignal, tPeakPos);

        const TVector3 firstGuess(EVector[0], EVector[1], EVector[2]);
        TVector3 meanEField = GetMeanEfieldInRange(stationtrace, fwhm, firstGuess);

        const Vector vCurrentEField =
          Normalized(Vector(meanEField.x(), meanEField.y(), meanEField.z(), coordinateSystem));

        const double meanAngleToLorentzForce =
          RadioGeometryUtilities::GetAngleToLorentzVector(vCurrentEField, showerAxis, magneticFieldAxis);

        const double meanAngleToExpectedEField = RadioGeometryUtilities::GetAngleToEFieldExpectation(
            vCurrentEField, showerCore, showerAxis, stationPosition, magneticFieldAxis, fChargeExcessStrength);

        // set magnitude of the mean efield vector to the maximum of the hilbert combo trace
        meanEField.SetMag(maxSignal);

        // save mean EFieldvector
        EVector[0] = meanEField.x();
        EVector[1] = meanEField.y();
        EVector[2] = meanEField.z();
        angleToLorentzForce = NormalizeAngle(meanAngleToLorentzForce);
        angleToExpectedEField = NormalizeAngle(meanAngleToExpectedEField);
      }

      info.str("");
      info << "Station " << station.GetId() << ": Calculated EFieldvector to ("
           << EVector[0] << ", " << EVector[1] << ", " << EVector[2] << ")";
      INFODebug(info);

      recStation.SetParameter(eEFieldVectorEW, EVector[0]);
      recStation.SetParameter(eEFieldVectorNS, EVector[1]);
      recStation.SetParameter(eEFieldVectorV, EVector[2]);

      const double angleToLorentzForceError =
        GetAngleToLorentzForceError(recStation.GetParameter(eSignalToNoise));

      recStation.SetParameter(eAngleToLorentzVector, angleToLorentzForce);
      recStation.SetParameterError(eAngleToLorentzVector, angleToLorentzForceError);
      recStation.SetParameter(eAngleToExpectedEFieldVector, angleToExpectedEField);
      recStation.SetParameterError(eAngleToExpectedEFieldVector, angleToLorentzForceError);

      info.str("");
      info << "Station " << station.GetId() << ": "
           << "Angle to Lorentz force vector " << round(angleToLorentzForce / degree * 10.) / 10.
           << " +/- " << round(angleToLorentzForceError / degree  * 10.) / 10. << " deg, "
           << "Angle to expected eField vector " <<  round(angleToExpectedEField / degree  * 10.) / 10.
           << " +/- " << round(angleToLorentzForceError / degree * 10.) / 10. << " deg "
           << "with a charge excess strength of " << fChargeExcessStrength;
      INFOIntermediate(info);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationEFieldVectorCalculator::Finish()
  {
    return eSuccess;
  }


  pair<double, double>
  RdStationEFieldVectorCalculator::GetFWHM(const StationTimeSeries& stationEnvelope, const double maxSignal,
                                           const double tPeakPos)
  {
    ostringstream info;

    unsigned int FWHMPosLeft = 0;
    for (auto i = tPeakPos; i > 0; --i) {
      if (stationEnvelope[i].GetMag() < maxSignal / 2) {
        FWHMPosLeft = i;
        break;
      }
    }

    unsigned int FWHMPosRight = 0;
    for (auto i = tPeakPos; i < stationEnvelope.GetSize(); ++i) {
      if (stationEnvelope[i].GetMag() < maxSignal / 2) {
        FWHMPosRight = i;
        break;
      }
    }

    info.str("");
    info << "FWHM: " << FWHMPosLeft << " = "
         << FWHMPosLeft * stationEnvelope.GetBinning() << "ns to " << FWHMPosRight
         << " = " << FWHMPosRight * stationEnvelope.GetBinning() << "ns";
    INFODebug(info);

    return std::make_pair(FWHMPosLeft, FWHMPosRight);
  }


  TVector3
  RdStationEFieldVectorCalculator::GetMeanEfieldInRange(const StationTimeSeries& stationTrace,
                                                        const pair<double, double>& range,
                                                        const TVector3& firstGuess)
  {
    ostringstream info;

    // waiting for cpp17
    double FWHMPosLeft, FWHMPosRight;
    tie(FWHMPosLeft, FWHMPosRight) = range;

    TVector3 meanEField;
    TVector3 currentEField;
    double sumOfWeights = 0;

    // loop through trace from FWHMPosLeft to FWHMPosRight and compute mean EFieldvector
    for (auto i = FWHMPosLeft; i < FWHMPosRight; ++i) {
      currentEField = TVector3(stationTrace[i][0], stationTrace[i][1], stationTrace[i][2]);

      if (currentEField.Angle(firstGuess) > 90 * degree) {
        info.str("");
        info << "i = " << i << " = " << i * stationTrace.GetBinning()
             << "ns: vector flipped";
        INFODebug(info);
        currentEField *= -1;
      }

      meanEField += currentEField;
      sumOfWeights += currentEField.Mag();
    }
    meanEField *= 1. / sumOfWeights;

    return meanEField;
  }


  double
  RdStationEFieldVectorCalculator::GetAngleToLorentzForceError(double snr)
  {
    const double sigma1 = (23.6 / sqrt(snr) - 4.3 / snr + 41.8 / pow(snr, 1.5)) * degree;
    const double sigma2 = 1 * deg;  // 1 degree alignment uncertainty
    return sqrt(Sqr(sigma1) + Sqr(sigma2));
  }


  double
  RdStationEFieldVectorCalculator::NormalizeAngle(const double x)
  {
    return (x > kPi / 2) ? kPi - x : x;
  }

}