RdStationAssociator.cc

Go to the documentation of this file.

#include "RdStationAssociator.h"

#include <fwk/CentralConfig.h>

#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/Minou.h>
#include <utl/PhysicalConstants.h>
#include <utl/Point.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/Trace.h>
#include <utl/UTMPoint.h>
#include <utl/TraceAlgorithm.h>

#include <rdet/RDetector.h>
#include <rdet/Station.h>
#include <rdet/Channel.h>
#include <rdet/RStationListManager.h>
#include <rdet/RSimulationStationListManager.h>

#include <evt/RadioSimulation.h>
#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <revt/REvent.h>
#include <revt/EventTrigger.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationTriggerData.h>
#include <revt/StationSimData.h>

#include <sevt/SEvent.h>
#include <sevt/StationTriggerData.h>


using namespace std;
using namespace evt;
using namespace det;
using namespace utl;
using namespace fwk;


namespace RdStationAssociator {

  VModule::ResultFlag
  RdStationAssociator::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationAssociator");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("AddVirtualStations").GetData(fAddVirtualStations);
    topBranch.GetChild("FirstStationId").GetData(fFirstStationId);
    topBranch.GetChild("ExtrapolateTheTailOfThePulses").GetData(fExtrapolateTheTailOfThePulses);
    topBranch.GetChild("MaximumAllowedDistance").GetData(fMaximumDistance);
    topBranch.GetChild("TracePaddingAtFront").GetData(fTracePaddingAtFront);
    topBranch.GetChild("MinimumTraceLength").GetData(fMinimumTraceLength);
    topBranch.GetChild("IncludedStationIds").GetData(fIncludedStationIds);
    topBranch.GetChild("ExcludedStationIds").GetData(fExcludedStationIds);
    topBranch.GetChild("AddNoSignalStations").GetData(fAddNoSignalStations);
    topBranch.GetChild("AddEmptyTriggeredRDStations").GetData(fAddEmptyTriggeredRDStations);
    topBranch.GetChild("MaximumDistanceOfNoSignalStations").GetData(fMaximumDistanceOfNoSignalStations);
    topBranch.GetChild("CreateAllSDStationsWithRDTrace").GetData(fCreateAllSDStationsWithRDTrace);

    ostringstream info;
    info << "\n\tMaximum allowed distance "<< fMaximumDistance/meter << " m\n"
         << "\tAdding time series data at front by " << fTracePaddingAtFront / nanosecond
         << " ns and at end to a total length of at least " << fMinimumTraceLength / nanosecond << " ns.\n"
         << "\tAdd virtual stations: " << fAddVirtualStations << "(" << fFirstStationId << ")\n";

    if (fAddNoSignalStations) {
      if (fAddEmptyTriggeredRDStations)
        info << "\tAdd empty traces for triggered but not simulated RD stations\n";
      else {
        info << "\tAdd empty traces for all stations within " << fMaximumDistanceOfNoSignalStations / utl::km
             << " km around the shower axis.\n";
      }
    }

    INFOFinal(info);
    return eSuccess;
  }


  VModule::ResultFlag
  RdStationAssociator::Run(Event& theevent)
  {
    ostringstream out;

    if (fAddEmptyTriggeredRDStations && !theevent.HasSEvent()) {
      ERROR("SEvent not available but required!");
      return eFailure;
    }

    auto& simshow = theevent.GetSimShower();
    auto& rsim = simshow.GetRadioSimulation();
    auto corsys = simshow.GetLocalCoordinateSystem();
    auto pulsecorsys = rsim.GetLocalCoordinateSystem();

    const auto& eventtime = simshow.GetTimeStamp();
    const auto& det = Detector::GetInstance();
    const auto& radioDet = det.GetRDetector();

    // Add generic stations to detector describtion at the exact positions of the simulated pulses
    if (fAddVirtualStations) {
      RdStationAssociator::AddVirtualStations(theevent, pulsecorsys);
    }

    // check if the event time has an uninitialized (illegal) value
    if (eventtime == TimeStamp(0, 0)) {
      WARNING("RdStationAssociator: Event time is not available. Event is skipped!");
      return eContinueLoop;
    }

    if (!theevent.HasREvent())
      theevent.MakeREvent();

    revt::REvent& therev = theevent.GetREvent();

    if (!rsim.GoToFirstSimRadioPulse()) {
      WARNING("RdStationAssociator: No SimRadioPulse in RadioSimulation. Event is skipped!");
      return eContinueLoop;
    }

    therev.GetHeader().SetId(rsim.GetEventNumber());
    therev.GetHeader().SetRunNumber(rsim.GetRunNumber());
    therev.GetHeader().SetTime(eventtime);

    /* Add Trigger Information to the Event. By default Event are assumed as Externally triggered
       Once somebody implements a module for RdTrigger simulation in OffLine this part
       should become obsolete */
    therev.MakeTrigger();
    therev.GetTrigger().SetExternalTrigger(true);
    // End of Trigger Information part

    if (fInfoLevel >= eInfoDebug) {
      out.str("");
      for (const auto& detStation : radioDet.StationsRange()) {
        auto position = detStation.GetPosition();
        out << "Station " << detStation.GetId() << ": " << position.GetX(corsys) << ", "
            << position.GetY(corsys) << ", "  << position.GetZ(corsys) << endl;
      }
      INFODebug(out);
    }

    bool ok;
    // association between RDetector station and sim. pulse
    do {
      const SimRadioPulse& srpulse = rsim.GetNextSimRadioPulse(ok);

      // leave the loop if there were no further SimRadioPulses
      if (!ok)
        break;

      srpulse.SetLocalCoordinateSystem(pulsecorsys);

      // check if there is a station in the close vicinity of the SimRadioPulse
      const utl::Point loc = srpulse.GetLocation();
      const int statID = FindClosestStationFromPoint(loc, radioDet, 1);

      if (statID != -1) {
        out.str("");
        out << "associated station " << statID << " ("
            << (radioDet.GetStation(statID).GetPosition() - loc).GetMag()  / meter << " m)";
        INFODebug(out);
      } else {
        out.str("");
        out << "Did not find any station at position ("
            << loc.GetX(corsys) << ", " << loc.GetY(corsys) << ", " << loc.GetZ(corsys) << ")";
        INFOFinal(out);

        // if there was no station corresponding to the current SimRadioPulse, go to next SimRadioPulse
        continue;
      }

      // if station id is listed in fExcludedStationIds, skip this simulated radio pulse
      // (it will not be re-used for the second-closest antenna station!)
      if (!fExcludedStationIds.empty())
        if (find(fExcludedStationIds.begin(), fExcludedStationIds.end(), statID) != fExcludedStationIds.end()) {
          out.str("");
          out << "Station " << statID << " not associated as it is in the exclude list.";
          INFOFinal(out);
          continue;
        }

      // if there is a list of fIncludedStationIds skip this simulated radio pulse if station id is not in the list
      if (!fIncludedStationIds.empty()) {
        if (find(fIncludedStationIds.begin(), fIncludedStationIds.end(), statID) == fIncludedStationIds.end()) {
          out.str("");
          out << "Station " << statID << " not associated as it is not in the include list.";
          INFOFinal(out);
          continue;
        }
      }

      // If a station does not have a valid detector describtion == not at least two AntennaHighGain channel, skip it
      // (will not be added to the event). It can happen that such stations are simulated for real events while
      // their real conterpart is not taken data.
      const auto& station = radioDet.GetStation(statID);
      const unsigned int numberOfValidConfiguredChannel = count_if(station.ChannelsBegin(), station.ChannelsEnd(),
        [&](auto const& channel) { return channel.GetChannelType() == "AntennaHighGain"; }
      );

      if (numberOfValidConfiguredChannel < 2) {
        out.str("");
        out << "Station " << statID << " is not associated with a pulse as"
            << " its description is incomplete! (Does not has two channel with type: \"AntennaHighGain\"";
        WARNING(out);
        continue;
      }

      // Create all SD stations already here if you activated the RD triggering.
      // Otherwise you can get errors during SD Reconstruction
      if (fCreateAllSDStationsWithRDTrace) {
        const int sStID = statID - radioDet.GetRdSdStationIdLink();

        sevt::SEvent& sEvent = theevent.GetSEvent();
        if (!sEvent.HasStation(sStID)) {
          sEvent.MakeStation(sStID);
          sEvent.GetStation(sStID).MakeSimData();
        }
      }

      // otherwise, fill SimRadioPulse into the associated station
      therev.MakeStation(statID); // Create a new Station, already checks if the station exists
      auto& stat = therev.GetStation(statID);

      // create data structures
      stat.MakeGPSData();     // should check before if exists?
      stat.MakeTriggerData(); // should check before if exists?
      stat.MakeSimData();     // should check before if exists?

      const auto& simtimeseries = srpulse.GetEfieldTimeSeries();

      // maybe have to check if there is data already and throw an exception?
      auto& statTimeSeries = stat.GetStationTimeSeries();

      // fill the Station SimData TimeSeries with the (unpadded) simulated time series data
      stat.GetSimData().SetSimulatedTrace(simtimeseries);

      // copy over the simtimeseries to the Station TimeSeries, padding appropriately at the front
      Vector3D v3DZero;
      v3DZero = 0.0, 0.0, 0.0;

      const double binning = simtimeseries.GetBinning();
      statTimeSeries.SetBinning(binning);
      long numpaddedsamples = static_cast<long>(fTracePaddingAtFront / binning + 0.5);

      // set the start time of the station trace, taking into account the front padding appropriately
      stat.SetRawTraceStartTime(therev.GetHeader().GetTime()
        + TimeInterval(srpulse.GetStartTime() - numpaddedsamples * binning));

      // pad the station time series trace at the front
      for (long i = 0; i < numpaddedsamples; ++i)
        statTimeSeries.PushBack(v3DZero);

      out.str("");
      out << "sim time series has size " << simtimeseries.GetSize();
      INFODebug(out);

      // copy over the simulated trace
      for (const auto& el : simtimeseries)
        statTimeSeries.PushBack(el);

      /* Extrapolating the tail of the sim. pulse (for each polarisation) migh be nessecary because
         the pulse has not converged to 0 in the time window it was simulated. This introduces a cut-off when
         padded with zeros at the end. Bandpass filtered, this cutoff can look like a regular pulse. With the
         extrapolation the pulse amplitude can be brought close to 0 before padding */
      if (fExtrapolateTheTailOfThePulses) {
        unsigned int fitRangeStart = statTimeSeries.GetSize() - int((50 * utl::ns) / binning);

        vector<double> vPowerOrig{0, 0, 0}; // to compare with extrapolated trace
        vector<vector<double>> extraPulse;

        for (unsigned int pol = 0; pol < 3; pol++) {
          vector<pair<double, double>> fitData;

          for (unsigned int sample = 0; sample < statTimeSeries.GetSize(); ++sample) {
            vPowerOrig.at(pol) += Sqr(statTimeSeries[sample][pol]);
            if (sample >= fitRangeStart)
              fitData.emplace_back(sample * binning, move(statTimeSeries[sample][pol]));
          }

          // Initalize parameters (WARNING: Take care of correct ordering)
          vector<utl::Minou::ParameterDef> parameters;
          // name, value, step, min, max, fixed
          parameters.emplace_back("a", fitData[0].second, 1e-6, -1, 1, false);
          parameters.emplace_back("b", -5e-3, 1e-3, -1, -1e-3, false);

          ExpoTail eT(parameters);
          eT.Setx0(fitRangeStart * binning);
          eT.SetData(fitData);
          utl::Minou::Minimizer<ExpoTail> min(eT);
          min.Minimize();

          const auto result = min.GetParameters();
          if (result.at(1) >= -1.00001e-3)
            ERROR("Extrapolating the radio pulse, slope at bondary");

          vector<double> extraPulsePol;
          for (int sample = 0; sample < int(2000 * utl::ns / binning); ++sample) {
            double time = (statTimeSeries.GetSize() + sample) * binning;
            extraPulsePol.push_back(eT.GetY(time, result));
          }
          extraPulse.push_back(extraPulsePol);
        }

        Vector3D v3D;
        for (int sample = 0; sample < int(2000 * utl::ns / binning); ++sample) {
          v3D = extraPulse.at(0).at(sample), extraPulse.at(1).at(sample), extraPulse.at(2).at(sample);
          statTimeSeries.PushBack(v3D);
        }

        for (unsigned int pol = 0; pol < 3; pol++) {
          double powerExt = 0;
          for (const auto& el : statTimeSeries)
            powerExt += Sqr(el[pol]);

          if ((fabs(powerExt / vPowerOrig.at(pol) - 1) > 0.1) &&
              (vPowerOrig.at(pol) * binning * kConversionRadioSignalToEnergyFluence > 1 * utl::eV / utl::m2)) {
            out.str("");
            out << "Due to the extrapolation of the pulse tail, the power (i.e., sum over squared amplitudes) "
                << "changed by more than 10%. Please be cautious. Station: "
                << statID << ", Polarisation: " << pol << ": (extrapolated/original) "
                <<  powerExt << " / " << vPowerOrig.at(pol) << " = " << powerExt / vPowerOrig.at(pol);
            WARNING(out);
          }
        }
      }

      /* Check amplitudes at the end of the pulse, the amplitudes should be close to zero.
         See also the comment above for fExtrapolateTheTailOfThePulses == True */
      Trace<double> statTimeSeriesMag;
      statTimeSeriesMag.SetBinning(statTimeSeries.GetBinning());
      for (const auto& el : statTimeSeries)
        statTimeSeriesMag.PushBack(el.GetMag());

      const double noiseRMS = TraceAlgorithm::RootMeanSquare(statTimeSeriesMag,
        statTimeSeriesMag.GetSize() - int(20 * utl::ns / binning), statTimeSeriesMag.GetSize() - 1);
      if (noiseRMS > 50 * utl::microvolt) {
        out.str("");
        out << "Found a RMS amplitude of " << noiseRMS / utl::microvolt
            << " muV/m in the last 20ns of the simulated pulse. "
            << "Such a cutoff might introduces a second pulse. "
            << "Condider using fExtrapolateTheTailOfThePulses";
        WARNING(out);
      }

      out.str("");
      out << "station time series has size " << statTimeSeries.GetSize();
      INFODebug(out);

      // now pad the station trace at the end up to the given minimum trace length
      while ((fMinimumTraceLength - statTimeSeries.GetSize() * binning) > 1e-6)
        statTimeSeries.PushBack(v3DZero);

      // if number of samples is uneven after padding, pad one more sample
      // (otherwise it will again be clipped at the next FFT)
      if (statTimeSeries.GetSize() % 2 == 1)
        statTimeSeries.PushBack(v3DZero);

    } while (true);

    // add empty stations with empty trace for stations that do not have been simulated but where
    // no signal is expected -> trace of zeros can be added
    if (fAddNoSignalStations) {
      unsigned numberOfEmptyButTriggeredRDStations = 0;

      for (const auto& station : radioDet.StationsRange()) {
        const unsigned int rdId = station.GetId();

        if (!fExcludedStationIds.empty()) {
          if (find(fExcludedStationIds.begin(), fExcludedStationIds.end(), station.GetId()) != fExcludedStationIds.end())
            continue;
        }

        if (!therev.HasStation(rdId)) {  // check is station has no data yet
          utl::Point positionStat = station.GetPosition();
          utl::Vector mcAxis = -1 * simshow.GetDirection();
          const Point& mcCore = simshow.GetPosition();

          const double dist = RadioGeometryUtilities::GetDistanceToAxis(
            mcAxis, mcCore, positionStat);

          if (fAddEmptyTriggeredRDStations) {
            const int sStID = rdId - radioDet.GetRdSdStationIdLink();

            const sevt::SEvent& sEvent = theevent.GetSEvent();
            if (!sEvent.HasStation(sStID) || !sEvent.GetStation(sStID).HasTriggerData())
              continue; // skip this station

            const auto& trig = sEvent.GetStation(sStID).GetTriggerData();
            if (!trig.IsT2() && !trig.IsT1())
              // This is not equal to station.IsSilent() but this is desired here
              continue; // skip this station

            // the station is being added
            numberOfEmptyButTriggeredRDStations++;
          } else {
            // decide upon distance
            if (dist > fMaximumDistanceOfNoSignalStations)
              continue; // skip this station
          }

          therev.MakeStation(rdId);  // Create a new Station, already checks if the station exists
          revt::Station& stat = therev.GetStation(rdId);

          // create data structures
          stat.MakeGPSData();  // should check before if exists?
          stat.MakeTriggerData();  // should check before if exists?
          stat.MakeSimData();  // should check before if exists?

          // maybe have to check if there is data already and throw an exception?
          auto& statTimeSeries = stat.GetStationTimeSeries();

          Vector3D v3DZero;
          v3DZero = 0.0, 0.0, 0.0;

          // hardcoded binning of 1ns, will anyway be resampled to detector sampling
          const double binning = 1 * utl::ns;
          statTimeSeries.SetBinning(binning);

          out.str("");
          out << "adding zero trace to station " << rdId;
          INFODebug(out);

          for(unsigned int i=0; i < (fMinimumTraceLength / binning); ++i) {
            statTimeSeries.PushBack(v3DZero);
          }

          // if number of samples is uneven after padding, pad one more sample
          // (otherwise it will again be clipped at the next FFT)
          if (statTimeSeries.GetSize() % 2  == 1) {
            statTimeSeries.PushBack(v3DZero);
          }

          // Calculating expected pulse time if event had a pulse. the trace start time is then set,
          // so that this time is within the trace
          const CoordinateSystemPtr coreCS = LocalCoordinateSystem::Create(mcCore);

          // calculate Rd station distance to Sd core
          const Vector VecDistanceToSdCore = positionStat - mcCore;

          // project the vector on shower direction
          const double projVec = -mcAxis.GetX(coreCS) * VecDistanceToSdCore.GetX(coreCS)
              - mcAxis.GetY(coreCS) * VecDistanceToSdCore.GetY(coreCS)
              - mcAxis.GetZ(coreCS) * VecDistanceToSdCore.GetZ(coreCS);

          // calculate expected pulse arrival time
          const TimeInterval relativePulseTime = projVec / (kSpeedOfLight) - fTracePaddingAtFront;

          out.str("");
          out << "Set relative pulse time of noSignalStation to " << relativePulseTime
              << " relative to rd event time";
          INFODebug(out);

          const TimeStamp traceStartTime = therev.GetHeader().GetTime() + relativePulseTime;
          stat.SetRawTraceStartTime(traceStartTime);

          out.str("");
          out << "created empty trace with " << statTimeSeries.GetSize() << " samples";
          INFODebug(out);
        }
      }

      INFOIntermediate(out);

      if (numberOfEmptyButTriggeredRDStations) {
        out.str("");
        out << numberOfEmptyButTriggeredRDStations
            << " empty traces have been add for triggered but not simulated RD stations.";
        INFO(out);
      }

    }

    return eSuccess;
  }


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


  int
  RdStationAssociator::FindClosestStationFromPoint(const utl::Point& pt,
                                                   const rdet::RDetector& rDet,
                                                   double maxDistanceFactor)
  {
    double dist = 0; // distance between point and station
    double foundmindist = fMaximumDistance * maxDistanceFactor;
    int stid = -1;

    for (const auto& station : rDet.StationsRange()) {
      const utl::Point Position = station.GetPosition();

      dist = (Position - pt).GetMag();
      if (dist < foundmindist) {
        stid = station.GetId(); // make current station the closest one
        foundmindist = dist;
      }
    }

    return stid;
  }


  void
  RdStationAssociator::AddVirtualStations(evt::Event& event, utl::CoordinateSystemPtr corsys)
   const
  {
    // initialize starshaped station with id (without proper position)
    INFODebug("Use RdStationAssociator to create generic (starshaped) stations.");
    ostringstream out;

    Detector& det = Detector::GetInstance();

    rdet::RSimulationStationListManager& listManager =
        dynamic_cast<rdet::RSimulationStationListManager&>
        (det.GetRManagerRegister().GetManager("RSimulationStationListManager"));

    if (!event.HasSimShower()) {
      const string err = "Attempt to generate radio array with starshaped pattern failed because "
                         "no simulated shower is present in the event.";
      throw MissingEventDataException(err);
    }

    const string fEllipsoid_string = "WGS84";
    const string fBand = "H";
    const int fZone = 19;

    const ReferenceEllipsoid ellipsoid = ReferenceEllipsoid::Get(
        ReferenceEllipsoid::GetEllipsoidIDFromString(fEllipsoid_string));

    RadioSimulation& radiosim = event.GetSimShower().GetRadioSimulation();
    const int numberOfAntennas = radiosim.GetNumPulses();
    unsigned int nstat = 0;
    // starting station id, increment ids from that on
    for (int stationId = fFirstStationId; stationId < fFirstStationId + numberOfAntennas; stationId++) {

      rdet::RStationListManager::StationData station;
      station.fSource = "XML";

      station.fId = stationId;
      station.fCommissionTime = "2000-01-01T00:00:00Z";
      station.fDecommissionTime = "2030-01-01T00:00:00Z";
      station.fNChannels = 2;
      station.fFirstChannelId = 0;
      station.fLastChannelId = 1;
      station.fEllipsoid = fEllipsoid_string;
      station.fBand = fBand;
      station.fZone = fZone;

      station.fIsVirtual = 1;

      // get index of pulse -> accosciation between arbitary station and pulse
      // the RdStationAssociator will make sure that everything is fine
      // however once a station id is accoiated with a pulse index this should not change!
      const int pulseIndex = stationId - fFirstStationId;

      if (radiosim.HasSimPulseByIndex(pulseIndex)) {
        // if the station was initialized as starshaped station
        const SimRadioPulse& srpulse = radiosim.GetSimPulseByIndex(pulseIndex);  // get specific SimRadioPulse
        srpulse.SetLocalCoordinateSystem(corsys);
        const Point pos = srpulse.GetLocation();
        const UTMPoint stationUTM(pos, ellipsoid);

        station.fNorthing = stationUTM.GetNorthing();
        station.fEasting = stationUTM.GetEasting();
        station.fAltitude = stationUTM.GetHeight();
        station.fName = srpulse.GetAntennaName();

        listManager.AddVirtualStation(station);
        nstat++;
      } else {
        ERROR("Something went wrong during reading of virtual stations");
      }

    }

    // Force update of the RDetector after virtual stations were added to RStationListManager
    det.Update(det.GetTime(), false, false, true);

    out.str("");
    out << "Initialized " << nstat << " virtual stations.";
    INFO(out);
  }

}