RdAntennaStationToChannelConverter.cc

Go to the documentation of this file.

#include "RdAntennaStationToChannelConverter.h"

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

#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Vector.h>
#include <utl/Point.h>

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

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

#include <det/Detector.h>

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

#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

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


namespace RdAntennaStationToChannelConverter {

VModule::ResultFlag
RdAntennaStationToChannelConverter::Init()
{
  Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdAntennaStationToChannelConverter");
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  INFODebug("");

  topBranch.GetChild("InterpolationMode").GetData(fInterpolationMode);
  std::string tmpstring = topBranch.GetChild("Direction").Get<std::string>();
  fDirectionPerStation = tmpstring == "LineOfSightStationShowerMaximum";

  topBranch.GetChild("CalculateDistanceToXmaxOnTheFly").GetData(fCalculateDistanceToXmaxOnTheFly);


  std::ostringstream info;
  info << "\n\tInterplation mode : " << fInterpolationMode
       << "\n\tDirection used    : " << tmpstring << "\n ";
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
RdAntennaStationToChannelConverter::Run(evt::Event& event)
{
  INFODebug("");

  // First some checks if simulation data is there:
  if (!event.HasSimShower()) {
    ERROR("Event has no simulated data.");
    return eContinueLoop;
  }
  // is there also radio simulation data?:
  evt::ShowerSimData& simShower = event.GetSimShower();
  if (!simShower.HasRadioSimulation()) {
    ERROR("SimShower has no radio simulation.");
    return eContinueLoop;
  }

  const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
  const Vector simAxis = -simShower.GetDirection();
  const double zenith = simAxis.GetTheta(localCS);
  const double azimuth = simAxis.GetPhi(localCS);

  utl::Point posShowerMaximum;
  if (fDirectionPerStation) {

    double distanceToXmax = -1;

    if (fCalculateDistanceToXmaxOnTheFly) {
      const auto& gh = simShower.GetGHParameters();
      const double xmaxMC = gh.GetXMax();

      // initiate atm + tables (which model is used is decieded in the config)
      const atm::Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
      theAtm.InitSlantProfileModel(simShower.GetPosition(), simAxis, 5 * g / cm2);
      const atm::ProfileResult& distanceFromDepth = theAtm.EvaluateDistanceVsSlantDepth();

      // distance to shower maximum, yep, the result is negative, convention thing...
      distanceToXmax = abs(distanceFromDepth.Y(xmaxMC));
    } else {
      distanceToXmax = simShower.GetDistanceOfShowerMaximum();
      if (distanceToXmax < 0) {
        std::ostringstream err;
        err << "Distance to shower maximum (as stored in the .reas file) is negative. "
              "This might happens for mpi generated simulations. Using \"MCAxis\" in "
              " the XML config does not need the distance to the shower maximum. Abort...";
        ERROR(err);
        throw utl::NonExistentComponentException(err.str());
      }
    }
    posShowerMaximum = simShower.GetPosition() + simAxis * distanceToXmax;
  }

  // has the radio simulation data been associated with the stations?:
  if (!event.HasREvent()) {
    ERROR("Event has simulated shower but no radio event. Use RdStationAssociator!");
    return eContinueLoop;
  }

  // We need it all: rEvent, rdet and the simulated data... fetch the rest:
  REvent& rEvent = event.GetREvent();
  const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

  // we need to loop over all stations and also over all channels of each station! Here is the loop
  // over the stations which are included in the event.
  for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {

    double stationZenith = -1;
    double stationAzimuth = -1;
    if (fDirectionPerStation) {
      const Point& posStation = rDetector.GetStation(*sIt).GetPosition();
      const CoordinateSystemPtr stationCS = LocalCoordinateSystem::Create(posStation);
      const Vector stationToShowerMaximum = posShowerMaximum - posStation;
      stationZenith = stationToShowerMaximum.GetTheta(stationCS);
      stationAzimuth = stationToShowerMaximum.GetPhi(stationCS);
    } else {
      stationZenith = zenith;
      stationAzimuth = azimuth;
    }

    if (fInfoLevel >= eInfoIntermediate) {
      std::ostringstream info;
      info << "Direction used for station " << sIt->GetId() << ": "
          << "zenith = " << stationZenith / utl::deg << " deg"
          << ", azimuth = " << stationZenith / utl::deg << " deg";
      INFO(info);
    }

    // for convenience, get a reference to the source station time series
    const StationTimeSeries& efield = sIt->GetStationTimeSeries();

    // we will need to convert the station data to the shower coordinate system, for this we create a temporary FFTDataContainer
    StationFFTDataContainer StationDataSCS;
    // we have to configure the Nyquist zone of the temporary container correctly
    StationDataSCS.SetNyquistZone(sIt->GetNyquistZone());

    // at each station there should be a three component electric field vector in cartesian coordinates
    // this vector is now transfered into the shower base coordinate system with the antenna in the center.
    // Here the simulated shower incoming direction is used which seems to be a pretty fair estimate:
    StationTimeSeries& eFieldSCS = StationDataSCS.GetTimeSeries();
    ConvertToShowerCS(efield, eFieldSCS, stationZenith, stationAzimuth);

    // Now we have the eField in SCS (Shower Coordinate System) in the time domain.
    // However, we need the FFT of the three components.
    StationFrequencySpectrum& eFieldSCSSpec = StationDataSCS.GetFrequencySpectrum();

    // now that we have all prepared on the station level, we have to loop the channels to apply the transfer
    // function of the attached antenna. So loop the channels of the revent:
    for (Station::ChannelIterator cIt = sIt->ChannelsBegin(); cIt != sIt->ChannelsEnd(); ++cIt) {
      // Some preparation...
      // ...from the rdetector we will need the corresponding detector channel:
      const rdet::Channel& cDetChannel = rDetector.GetStation(*sIt).GetChannel(*cIt);

      // ...the calculated response will be saved as the channel frequency spectrum:
      ChannelFrequencySpectrum& cChannelResponseSpectrum = cIt->GetChannelFrequencySpectrum();

      // clear the channel response spectrum in case there is already data there
      cChannelResponseSpectrum.Clear();

      // set the meta data for the channel correctly, will be same as that of the original station data
      cChannelResponseSpectrum.SetBinning(sIt->GetStationFrequencySpectrum().GetBinning());
      cIt->SetNyquistZone(sIt->GetNyquistZone());

      const double designlowerfrequency = cDetChannel.GetDesignLowerFreq();
      const double designupperfrequency = cDetChannel.GetDesignUpperFreq();

      // ...this is needed to avoid multiple interpolations at the same frequency (here for technical reason)
      std::pair<std::complex<double>, std::complex<double>> cEffectiveAntennaHeight;

      // now we will loop over the efield spectrum and fold it with the transfer function:
      for (StationFrequencySpectrum::SizeType i = 0; i < eFieldSCSSpec.GetSize(); ++i) {
        // fetch the effective antenna height for the incoming direction and the frequency which belongs to the bin:

        const double frequency = StationDataSCS.GetFrequencyOfBin(i);
        if (frequency >= designlowerfrequency && frequency <= designupperfrequency) {
          cEffectiveAntennaHeight = cDetChannel.GetElectricFieldResponse(stationZenith, stationAzimuth,
                                                                         StationDataSCS.GetFrequencyOfBin(i),
                                                                         fInterpolationMode);

          // multiply the e_theta components of the effective heigths and the efield as well as the e_phi components and add them up
          // this is then the antenna response at the corresponding frequency. So push it into the channel response spectrum:
          cChannelResponseSpectrum.PushBack(
              cEffectiveAntennaHeight.first * eFieldSCSSpec[i][0]
                  + cEffectiveAntennaHeight.second * eFieldSCSSpec[i][1]);
        } else {
          cChannelResponseSpectrum.PushBack(std::complex<double>(0., 0.));
        }
      }
    }
  }

  return eSuccess;
}


void
RdAntennaStationToChannelConverter::ConvertToShowerCS(const StationTimeSeries& eFieldXYZ,
                                                      StationTimeSeries& eFieldSCS,
                                                      double sZen, double sAzi)
  const
{
  eFieldSCS.Clear();  // delete data if there is already some
  eFieldSCS.SetBinning(eFieldXYZ.GetBinning());  // have to set the binning correctly
  for (StationTimeSeries::SizeType i = 0; i < eFieldXYZ.GetSize(); ++i) {
    double tmpVector[3] = {
        (cos(sZen) * cos(sAzi) * eFieldXYZ[i][0] + cos(sZen) * sin(sAzi) * eFieldXYZ[i][1]
         - sin(sZen) * eFieldXYZ[i][2]),  // e_theta
        (-sin(sAzi) * eFieldXYZ[i][0] + cos(sAzi) * eFieldXYZ[i][1]),  // e_phi
        (sin(sZen) * cos(sAzi) * eFieldXYZ[i][0] + sin(sZen) * sin(sAzi) * eFieldXYZ[i][1]
         + cos(sZen) * eFieldXYZ[i][2])  // e_r
    };
    eFieldSCS.PushBack(Vector3D(tmpVector));
  }
}


VModule::ResultFlag
RdAntennaStationToChannelConverter::Finish()
{
  INFODebug("");
  return eSuccess;
}

}