RdAntennaChannelToStationConverter.cc

Go to the documentation of this file.

#include "RdAntennaChannelToStationConverter.h"

#include <complex>
#include <cmath>

#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerFRecData.h>
#include <revt/StationRecData.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <evt/ShowerSimData.h>

#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/StringCompare.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/MagneticFieldModel.h>

#include <det/Detector.h>

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

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

#include <fwk/CentralConfig.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>

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

namespace RdAntennaChannelToStationConverter {

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

  string usedDirectionStr = topBranch.GetChild("UsedDirection").Get<string>();

  if (StringEquivalent(usedDirectionStr, "RdDirection")) {
    fUsedDirection = eRdReconstruction;
  } else if (StringEquivalent(usedDirectionStr, "SdReconstruction")) {
    fUsedDirection = eSdReconstruction;
  } else if (StringEquivalent(usedDirectionStr, "FdReconstruction")) {
    fUsedDirection = eFdReconstruction;
  } else if (StringEquivalent(usedDirectionStr, "McTruth")) {
    fUsedDirection = eMcTruth;
  } else if (StringEquivalent(usedDirectionStr, "IndividualPerStation")) {
    fUsedDirection = eIndividualPerStation;
  } else if (StringEquivalent(usedDirectionStr, "Reference")) {
    fUsedDirection = eReference;
  } else if (StringEquivalent(usedDirectionStr, "LineOfSightStationShowerMaximum")) {
    fUsedDirection = eLineOfSightStationShowerMaximum;
  }

  fDirectionOfXmax = topBranch.GetChild("DirectionOfXmax").Get<string>();
  fXmaxEstimator = topBranch.GetChild("XmaxEstimator").Get<string>();

  topBranch.GetChild("CropFrequencies").GetData(fCropFreq);
  topBranch.GetChild("ExtendLowerLimitPercentage").GetData(fLowSmear);
  topBranch.GetChild("ExtendUpperLimitPercentage").GetData(fUpSmear);
  topBranch.GetChild("ZeroLimit").GetData(fZeroLimit);
  topBranch.GetChild("InterpolationMode").GetData(fInterpolationMode);

  topBranch.GetChild("rejectSaturatedStations").GetData(fRejectSaturatedStations);

  topBranch.GetChild("calculateWeightsFromExpectedEField").GetData(fCalculateWeightsFromExpectedEField);
  topBranch.GetChild("calculateWeightsFromAntennaPattern").GetData(fCalculateWeightsFromAntennaPattern);
  topBranch.GetChild("useUserDefinedWeight").GetData(fUseUserDefinedWeights);
  ostringstream sstr;
  if (fUseUserDefinedWeights) {
    topBranch.GetChild("weight12").GetData(fWeight12);
    topBranch.GetChild("weight13").GetData(fWeight13);
    topBranch.GetChild("weight23").GetData(fWeight23);
    sstr << "\n\tUser defined weights for the calculation of the electric field out of 3 channel"
         << " will be used, weights are:"
         << " w12: " << fWeight12 << ", w13: " << fWeight13 << ", w23: " << fWeight23;
    INFOFinal(sstr);
  }

  string whichEye = topBranch.GetChild("UsedEye").Get<string>();
  if (whichEye == "CO")
    fWhichEye = 4;
  else if (whichEye == "HE")
    fWhichEye = 5;
  else if (whichEye == "HECO")
    fWhichEye = 6;
  else
    ERROR("<UsedEye> has invalid config.");

  sstr.str("");
  sstr << "\n\tUse direction             : " << usedDirectionStr
       << "\n\tUse interpolation method  : " << fInterpolationMode
       << "\n\tReject saturated stations : " << fRejectSaturatedStations << endl;
  INFOFinal(sstr);

  if ((fCalculateWeightsFromAntennaPattern + fCalculateWeightsFromExpectedEField + fUseUserDefinedWeights) != 1) {
    ERROR("Only one option to calculate weights must be chooses.");
    return eFailure;
  }

  return eSuccess;
}


VModule::ResultFlag
RdAntennaChannelToStationConverter::Run(evt::Event& event)
{
  // check if radio event information exists
  if (!event.HasREvent()) {
    ERROR("Event has no radio event!");
    return eContinueLoop;
  }

  ostringstream msg;

  // Declare some variable which will be used inside a loop to save time
  double channelfrequency = 0.;
  Vector3C V3CZero;  // neutral 3d complex vector for initialisations

  // temporal dummy variable holding the result of the E-Field calculation.
  // Later transferred to the station frequency spectrum
  Vector3C efield;
  pair<complex<double>, complex<double> > channel1Response;
  pair<complex<double>, complex<double> > channel2Response;
  V3CZero = complex<double>(0.0), complex<double>(0.0), complex<double>(0.0);
  bool e1flag = true;
  bool e2flag = true;

  // set the direction up with the fall-back values
  const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance()
      .GetReferenceCoordinateSystem();

  double azimuth = 0.;
  double zenith = 0.;
  utl::Vector showerAxis;

  // Get direction to be used to apply the antenna pattern.
  // If the function returns false, something went wrong.
  if (!GetDirection(event, showerAxis, azimuth, zenith))
    return eFailure;

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

  // counter. counts how many times the inversion of
  // the linear equation system has failed
  unsigned int inversionfail = 0;

  for (auto& station : rEvent.StationsRange()) {

    StationFrequencySpectrum& stationspectrum = station.GetStationFrequencySpectrum();

    // initialize station frequency spectrum
    for (StationFrequencySpectrum::SizeType i = 0; i < stationspectrum.GetSize(); ++i) {
      stationspectrum[i][0] = complex<double>(0.0, 0.0);
      stationspectrum[i][1] = complex<double>(0.0, 0.0);
      stationspectrum[i][2] = complex<double>(0.0, 0.0);
    }

    // identify active channels
    vector<revt::Channel> usableChannels;
    for (auto& channel : station.ChannelsRange()){
      if (channel.IsActive()) {
        usableChannels.push_back(channel);

        if (channel.IsSaturated()) {
          station.SetSaturated();
          msg.str("");
          msg << "detected saturated channel " << channel.GetId()
              << " for station " << station.GetId() << ".";

          if(fRejectSaturatedStations) {
            msg << " Rejecting station.";
            station.SetRejectedReason(revt::eSaturated);
          }

          INFOFinal(msg);
        }
      }
    }

    // if station has no channels (i.e. no data) exclude it ( in general this should not happen)
    if (usableChannels.size() == 0) {
      msg.str("");
      msg << "Station " << station.GetId()
          << " is excluded as it does not have any active channels"
          << " (actually this should never happen!)";
      WARNING(msg);
      station.SetExcludedReason(revt::eHardwareMismatch);
      continue;
    }

    // if individual directions for each station are used, get this direction now
    if (fUsedDirection == eIndividualPerStation) {

      if (!station.HasRecData())
        ERROR("Radio station does not have RecData, but it should have by now...");
      const StationRecData& recStation = station.GetRecData();

      if (!recStation.HasParameter(eSignalArrivalDirectionX) ||
          !recStation.HasParameter(eSignalArrivalDirectionY) ||
          !recStation.HasParameter(eSignalArrivalDirectionZ)) {
        msg.str("");
        msg << "Module is configured to use individual signal arrival directions for each station, "
            << "but no direction is set for station " << station.GetId() << ".";
        ERROR(msg);
        return eFailure;
      } else {
        azimuth = recStation.GetSignalArrivalAzimuth();
        zenith = recStation.GetSignalArrivalZenith();
      }
    } else if (fUsedDirection == eLineOfSightStationShowerMaximum) {

      Point xMaxCore;
      Vector xMaxAxis;

      if (fDirectionOfXmax == "MC") {
        if (event.HasSimShower()) {
          xMaxCore = event.GetSimShower().GetPosition();
          xMaxAxis = -event.GetSimShower().GetDirection();
        } else {
          ERROR("Module is configured to use Monte Carlo truth direction, but event has no SimShower!");
          return eFailure;
        }
      } else if (event.HasRecShower()) {
        if (fDirectionOfXmax == "Sd") {
          if (!event.GetRecShower().HasSRecShower()) {
            ERROR("Module is configured to use SD direction, but event has no SRecShower!");
            return eFailure;
          } else {
            xMaxCore = event.GetRecShower().GetSRecShower().GetCorePosition();
            xMaxAxis = event.GetRecShower().GetSRecShower().GetAxis();
          }
        }
        else if (fDirectionOfXmax == "Reference") {
          evt::ShowerRRecData& rrecShower = event.GetRecShower().GetRRecShower();
          xMaxAxis = rrecShower.GetReferenceAxis(event);
          xMaxCore = rrecShower.GetReferenceCorePosition(event);
        } else {
          ERROR("The configured DirectionOfXmax was invalid.");
          return eFailure;
        }
      } else {
        ERROR("The configured DirectionOfXmax was invalid.");
        return eFailure;
      }

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

      // Estimate Xmax to estimate start value for dmax
      const double showerXmaxGuess = GetXmaxEstimator(event);
      if (isnan(showerXmaxGuess) || showerXmaxGuess / (g/cm2) <= 1) {
        msg.str("");
        msg << "Could not get xmax estimator \"" << fXmaxEstimator << "\": "
            << showerXmaxGuess / (g / cm2);
        WARNING(msg);
        return eFailure;
      }

      // distance to shower maximum,  yep, the result is negative, convention thing...
      const double xMaxDistance = abs(distanceFromDepth.Y(showerXmaxGuess));
      const Point xMaxPos = xMaxCore + xMaxAxis * xMaxDistance;

      const rdet::Station& detStation = rDetector.GetStation(station);
      const Point& stationPos = detStation.GetPosition();
      const CoordinateSystemPtr stationCS =
        fwk::LocalCoordinateSystem::Create(stationPos);

      const Vector lineXmaxStation = xMaxPos - stationPos;
      azimuth = lineXmaxStation.GetPhi(stationCS);
      zenith = lineXmaxStation.GetTheta(stationCS);

    }

    // check for invalid values (NANs) - can happen when the direction reconstruction
    // failed but the reconstruction is continued
    if (std::isnan(azimuth) || std::isnan(zenith)) {
      msg.str("");
      msg << "Illegal azimuth / zenith: " << azimuth / degree << " / " << zenith / degree
          << " deg. Check (reference) source / reconstruction.";
      ERROR(msg);
      return eContinueLoop;
    }

    if (usableChannels.size() == 2) {
      INFOIntermediate("Only two channels selected, setting weight12 to 1, the other weights to zero");
      fWeight12 = 1.;
      fWeight13 = 0;
      fWeight23 = 0;
    } else if (usableChannels.size() == 3) {

      // calculate weights out of expected efield vector
      if (fCalculateWeightsFromExpectedEField) {
        utl::Vector magneticFieldAxis =
            event.GetRecShower().GetRRecShower().GetMagneticFieldVector();
        utl::Vector efieldExpectation = utl::RadioGeometryUtilities::GetLorentzVector(
            showerAxis, magneticFieldAxis);

        msg.str("");
        msg << "Shower axis is (" << showerAxis.GetR(referenceCS) << ", "
            << showerAxis.GetTheta(referenceCS) << ", " << showerAxis.GetPhi(referenceCS) << ")";
        INFODebug(msg);

        msg.str("");
        msg << "Efield expectation is (" << efieldExpectation.GetR(referenceCS) << ", "
            << efieldExpectation.GetTheta(referenceCS) << ", "
            << efieldExpectation.GetPhi(referenceCS) << ")";
        INFODebug(msg);

        const revt::Channel& channel1 = usableChannels.at(0);
        const revt::Channel& channel2 = usableChannels.at(1);
        const revt::Channel& channel3 = usableChannels.at(2);
        const rdet::Channel& cDetChannel1 = rDetector.GetStation(station.GetId()).GetChannel(
            channel1.GetId());
        const rdet::Channel& cDetChannel2 = rDetector.GetStation(station.GetId()).GetChannel(
            channel2.GetId());
        const rdet::Channel& cDetChannel3 = rDetector.GetStation(station.GetId()).GetChannel(
            channel3.GetId());

        utl::Vector orientationChannel1 = utl::Vector(1., cDetChannel1.GetOrientationZenith(),
                                                      cDetChannel1.GetOrientationAzimuth(),
                                                      referenceCS, utl::Vector::kSpherical);
        utl::Vector orientationChannel2 = utl::Vector(1., cDetChannel2.GetOrientationZenith(),
                                                      cDetChannel2.GetOrientationAzimuth(),
                                                      referenceCS, utl::Vector::kSpherical);
        utl::Vector orientationChannel3 = utl::Vector(1., cDetChannel3.GetOrientationZenith(),
                                                      cDetChannel3.GetOrientationAzimuth(),
                                                      referenceCS, utl::Vector::kSpherical);

        double weight1 = efieldExpectation * orientationChannel1;
        double weight2 = efieldExpectation * orientationChannel2;
        double weight3 = efieldExpectation * orientationChannel3;
        msg.str("");
        msg << "weight 1 = " << weight1 << " weight2 = " << weight2 << " weight3 = " << weight3;
        INFODebug(msg);

        fWeight12 = sqrt(weight1 * weight1 + weight2 * weight2);  // weight for channel 1 and 2
        fWeight13 = sqrt(weight1 * weight1 + weight3 * weight3);  // weight for channel 1 and 3
        fWeight23 = sqrt(weight2 * weight2 + weight3 * weight3);  // weight for channel 2 and 3

        msg.str("");
        msg << "using expected e-field vector to calculate weights\n"
            << "\tweight12 = " << fWeight12 << "\n"
            << "\tweight23 = " << fWeight23 << "\n"
            << "\tweight13 = " << fWeight13 << "";
        INFOIntermediate(msg);

      } else if (!fUseUserDefinedWeights) {
        ERROR("Calculation of weights for 3 usable channels is only supported, if "
          "'CalculateWeightsFromExpectedEField' or 'fUseUserDefinedWeights' "
          "is set in XML configuration.");
        return eFailure;
      }

    } else {
      msg.str("");
      msg << "Station " << station.GetId()
          << " does not have exactly two or three active channels! Have you forgotten"
          << " to include the RdChannelSelector?";
      ERROR(msg);
      throw utl::DoesNotComputeException(
          "RdAntennaChannelToStationConverter does not support stations with"
          " number of Channels != 2 or != 3.");
    }

    for (int iPolarization = 0; iPolarization < 3; ++iPolarization) {
      // stop loop if there is nothing to do
      if (iPolarization == 0 && fWeight12 == 0) {
        continue;
      }
      if (iPolarization == 1 && fWeight13 == 0) {
        continue;
      }
      if (iPolarization == 2 && fWeight23 == 0) {
        continue;
      }

      short channelId1 = 0;
      short channelId2 = 0;
      if (iPolarization == 0) {
        channelId1 = 0;
        channelId2 = 1;
      } else if (iPolarization == 1) {
        channelId1 = 0;
        channelId2 = 2;
      } else if (iPolarization == 2) {
        channelId1 = 1;
        channelId2 = 2;
      }

      const revt::Channel& channel1 = usableChannels.at(channelId1);
      const revt::Channel& channel2 = usableChannels.at(channelId2);

      // reference to the first channel frequency spectrum for calculation
      const ChannelFrequencySpectrum& channel1spectrum =
          channel1.GetChannelFrequencySpectrum();
      // reference to the second channel frequency spectrum for calculation
      const ChannelFrequencySpectrum& channel2spectrum =
          channel2.GetChannelFrequencySpectrum();

      // check if both channel frequency spectra have the same size and binning
      // maybe we have to check CloseTo<...> for binning?
      if ((channel1spectrum.GetSize() != channel2spectrum.GetSize())
          || (channel1spectrum.GetBinning() != channel2spectrum.GetBinning())
          || (channel1.GetNyquistZone() != channel2.GetNyquistZone())) {
        if (fInfoLevel >= eInfoIntermediate) {
          msg.str("");
          msg << "The 2 channels of station " << channel1.GetStationId()
              << " have different sizes, binnings or Nyquist zones!";
          WARNING(msg);
        }
        return eFailure;
      }

      // initialize empty station frequency spectrum
      stationspectrum = StationFrequencySpectrum(channel1spectrum.GetSize(),
                                                 channel1spectrum.GetBinning(), V3CZero);

      // set the Nyquist zone of the station
      station.SetNyquistZone(channel1.GetNyquistZone());

      // up to here only the station and the station frequency spectrum are initialized

      // access to the channel detector information of channel 1
      const rdet::Channel& cDetChannel1 = rDetector.GetStation(station.GetId()).GetChannel(
          channel1.GetId());
      // access to the channel detector information of channel 2
      const rdet::Channel& cDetChannel2 = rDetector.GetStation(station.GetId()).GetChannel(
          channel2.GetId());

      // for later check if in frequency range
      const double lowerfrequency = cDetChannel1.GetDesignLowerFreq() / 100 * (100 - fLowSmear);
      // given percentage smaller than actual frequency limit is accepted
      const double upperfrequency = cDetChannel1.GetDesignUpperFreq() / 100 * (100 + fUpSmear);
      // given percentage bigger than actual frequency limit is accepted

      const bool zenithBelowHorizon = zenith < 90. * degree;

      if (!zenithBelowHorizon) {
        msg.str("");
        msg << "Station " << station.GetId()
            << ": Zenith bigger than 90 degree, will use antenna response for 89 degree.";
        WARNING(msg);
      }

      // loop over the station frequency spectrum
      // for every bin in the station frequency spectrum, the "real" e-Field is calculated
      // from the channel frequency spectra
      for (StationFrequencySpectrum::SizeType i = 0; i < stationspectrum.GetSize(); ++i) {

        // access to the frequency of the first channel at bin i
        channelfrequency = channel1.GetFrequencyOfBin(i);
        // needed for the call to the response function of the channel

        if ((!fCropFreq) || ((channelfrequency > lowerfrequency) && (channelfrequency < upperfrequency))) {

          // retrieval of the channel response
          // (pair of complex number / for each E-Field component one)
          // for channel 1
          if (zenithBelowHorizon) {
            channel1Response = cDetChannel1.GetElectricFieldResponse(zenith, azimuth,
                                                                     channelfrequency,
                                                                     fInterpolationMode);
            // for channel 2
            channel2Response = cDetChannel2.GetElectricFieldResponse(zenith, azimuth,
                                                                     channelfrequency,
                                                                     fInterpolationMode);
          } else {
            // for channel 1
            channel1Response = cDetChannel1.GetElectricFieldResponse(89. * degree, azimuth,
                                                                     channelfrequency,
                                                                     fInterpolationMode);
            // for channel 2
            channel2Response = cDetChannel2.GetElectricFieldResponse(89. * degree, azimuth,
                                                                     channelfrequency,
                                                                     fInterpolationMode);
          }

          stringstream sstr;
          sstr.str("");
          sstr << "electric field response for " << zenith / utl::deg << " , " << azimuth / utl::deg
               << ", " << channelfrequency / utl::MHz << "), channel1 = "
               << abs(channel1Response.first) << ", " << abs(channel1Response.second)
               << " channel 2 = " << abs(channel2Response.first) << ", "
               << abs(channel2Response.second);
          INFODebug(sstr.str());

          // calculate weights out of channel response
          if (fCalculateWeightsFromAntennaPattern) {
            WARNING("This option <calculateWeightsFromAntennaPattern> is not implemented yet.");
          }

          // calculation of the clear E-Field without the channel response
          // we solve the linear equation system in two ways.
          // first checking if one of the response values is close to 0.
          // in this case one of the two equations can be solved for one component
          // imediately, makeing it easier to solve the other one.
          // additionally it is checked, if a division of 0 (near 0 values of
          // the response) rises, and if that can be avoided by rearranging
          // the other function, thus preventing the division.
          // This is possible since the linear equation system can be solved
          // by either rearranging the first equation and put it into the second
          // OR by rearranging the second equation and put it into the first.
          // Both are valid methods to solve the system.
          //
          // if the inversion is not possible, it is asumed, that the
          // corresponding efield component is 0
          // Keep this in mind:
          // IT IS NOT NECESSARILY TRUE THAT THE E-FIELD COMPONENT IS 0
          // JUST BECAUSE THERE WAS NO SOLUTION FOR THE LINEAR EQUATION SYSTEM

          // Initialisation of the E-Field values to be calculated
          efield[0] = complex<double>(0.0, 0.0);
          efield[1] = complex<double>(0.0, 0.0);

          // in a first step we try to solve the linear equation system by
          // putting the rearranged first equation into the second
          // if this fails we later try to put the second equation into
          // the first one
          //
          // Flags, indicating if the first or the second efield component
          // can be determined directly by the first equation
          e1flag = true;
          e2flag = true;
          // flag indicating if the inversion of the equation system failed
          // inversionflag = false;
          // check if the first response of the first equation is near 0,
          // and if so, calculate directly the result for the second
          // efield component
          if (abs(channel1Response.first) < fZeroLimit
              && abs(channel1Response.second) > fZeroLimit) {
            efield[1] = complex<double>(channel1spectrum[i] / channel1Response.second);
            e2flag = false;
          }
          // check if the second response of the first equation is near 0,
          // and if so, calculate directly the result for the first
          // efield component
          if (abs(channel1Response.second) < fZeroLimit
              && abs(channel1Response.first) > fZeroLimit) {
            efield[0] = complex<double>(channel1spectrum[i] / channel1Response.first);
            e1flag = false;
          }

          // checking if the second efield component was calculated directly,
          // if not calculate it now
          if (e2flag) {
            // check if the divisor is near 0 to prevent a division by 0,
            // if yes then raise a response warning, since the inversion
            // is not possible
            if (abs(
                channel1Response.first * channel2Response.second
                    - channel1Response.second * channel2Response.first) < fZeroLimit) {
              efield[1] = complex<double>(0.0, 0.0);
              // count up the response warnings.
              //inversionflag = true;
              ++inversionfail;
            }
            // if no then calculate the second component of the efield
            else {
              efield[1] = complex<double>(
                  (channel1Response.first * channel2spectrum[i]
                      - channel1spectrum[i] * channel2Response.first)
                      / (channel1Response.first * channel2Response.second
                          - channel1Response.second * channel2Response.first));
            }
          }

          // checking if the first efield component was calculated directly,
          // if not calculate it now
          if (e1flag) {
            // check if the divisor is near 0 to prevent a division by 0,
            // if yes then raise a response warning, since the inversion
            // is not possible
            if (abs(channel1Response.first) < fZeroLimit) {
              efield[0] = complex<double>(0.0, 0.0);
              // count up the response warnings.
              //inversionflag = true;
              ++inversionfail;
            }
            // if no then calculate the first component of the efield
            else {
              efield[0] = complex<double>(
                  (channel1spectrum[i] / channel1Response.first)
                      - (channel1Response.second * efield[1] / channel1Response.first));
            }
          }

          // if the equation system was not solved up to now, due to a divison
          // by 0 we have to try the second option: rearranging equation two and
          // putting it into the first equation
          if (inversionfail) {

            // check if the first response of the second equation is near 0,
            // and if so, calculate directly the result for the second
            // efield component
            if (abs(channel2Response.first) < fZeroLimit
                && abs(channel2Response.second) > fZeroLimit) {
              efield[1] = complex<double>(channel2spectrum[i] / channel2Response.second);
              e2flag = false;
            }

            // check if the second response of the second equation is near 0,
            // and if so, calculate directly the result for the second
            // efield component
            if (abs(channel2Response.second) < fZeroLimit
                && abs(channel2Response.first) > fZeroLimit) {
              efield[0] = complex<double>(channel2spectrum[i] / channel2Response.first);
              e1flag = false;
            }

            // checking if the second efield component was calculated directly,
            // if not, calculate it now
            if (e2flag) {
              // check if the divisor is near 0 to prevent a division by 0,
              // if yes then raise a response warning, since the inversion
              // is not possible
              if (abs(
                  channel2Response.first * channel1Response.second
                      - channel2Response.second * channel1Response.first) < fZeroLimit) {
                efield[1] = complex<double>(0.0, 0.0);
                // count up the response warnings.
                ++inversionfail;
              }
              // if no then calculate the first component of the efield
              else {
                efield[1] = complex<double>(
                    (channel2Response.first * channel1spectrum[i]
                        - channel2spectrum[i] * channel1Response.first)
                        / (channel2Response.first * channel1Response.second
                            - channel2Response.second * channel1Response.first));
              }
            }
            // checking if the first efield component was calculated directly,
            // if not, calculate it now
            if (e1flag) {
              // check if the divisor is near 0 to prevent a division by 0,
              // if yes then raise a response warning, since the inversion
              // is not possible
              if (abs(channel2Response.first) < fZeroLimit) {
                efield[0] = complex<double>(0.0, 0.0);
                // count up the response warnings.
                ++inversionfail;
              }
              // if no then calculate the first component of the efield
              else {
                efield[0] = complex<double>(
                    (channel2spectrum[i] / channel2Response.first)
                        - (channel2Response.second * efield[1] / channel2Response.first));
              }
            }
          }

        } else {
          efield[0] = complex<double>(0.0, 0.0);
          efield[1] = complex<double>(0.0, 0.0);
        }

        // set the third component of the E-Field to 0. Just initialisation
        efield[2] = complex<double>(0.0, 0.0);

        // rotation of the coordinate system from
        // the shower coordinate system (as presented here) to the
        // global coordinate system at the station

        // rotation around y-axis
        efield = rotateYaxis(efield, zenith);
        // rotation around z-axis
        efield = rotateZaxis(efield, azimuth);

        // transfer of the calculated E-Field at bin (frequency) i
        // to the station frequency spectrum at bin i.
        double sumOfWeights = fWeight12 + fWeight13 + fWeight23;
        if (iPolarization == 0) {
          efield *= fWeight12 / sumOfWeights;
          stationspectrum[i] += efield;
        }
        if (iPolarization == 1) {
          efield *= fWeight13 / sumOfWeights;
          stationspectrum[i] += efield;
        }
        if (iPolarization == 2) {
          efield *= fWeight23 / sumOfWeights;
          stationspectrum[i] += efield;
        }
      }  // station frequency spectrum loop end
    }  //  loop 3D-Antenna ends
  }  // station loop end

  // output of the response warnings during the processing of the run
  if (inversionfail != 0) {
    if (fInfoLevel >= eInfoIntermediate) {
      msg << "Linear equation solver failed for " << inversionfail << " frequencies!\n"
          << "This is in general no problem,"
          << "if the linear equation system could be solved"
          << "for at least some frequencies." << endl;
      WARNING(msg);
    }
    msg.str("");
  }

  return eSuccess;
}

// comment by TH: It could make sense to move the rotation functionality into the SVector class
//                (for the specialized case of 3 components) and also use it from inside the
//                RdAntennaStationToChannelConverter module.

// function to invert the transformation to the shower coordinate system.
// first part, rotation around y-axis.
// separated into two functions for better transparency and flexibility
utl::Vector3C RdAntennaChannelToStationConverter::rotateYaxis(utl::Vector3C vecC, double angle)
{
  utl::Vector3C resC;

  // transformation matrix used with input vector
  // canonical rotation matrix around the Y axis
  resC[0] = complex<double>(vecC[0] * cos(angle) + vecC[2] * sin(angle));
  resC[1] = complex<double>(vecC[1]);
  resC[2] = complex<double>(vecC[2] * cos(angle) - vecC[0] * sin(angle));

  return resC;
}

// function to invert the transformation to the shower coordinate system.
// second part, rotation around z-axis.
// separated into two functions for better transparency and flexibility
utl::Vector3C RdAntennaChannelToStationConverter::rotateZaxis(utl::Vector3C vecC, double angle)
{
  utl::Vector3C resC;

  // transformation matrix used with input vector
  // canonical rotation matrix around the X axis
  resC[0] = complex<double>(vecC[0] * cos(angle) - vecC[1] * sin(angle));
  resC[1] = complex<double>(vecC[0] * sin(angle) + vecC[1] * cos(angle));
  resC[2] = complex<double>(vecC[2]);

  return resC;

}

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

bool RdAntennaChannelToStationConverter::GetDirection(
    const evt::Event& event, utl::Vector& showerAxis, double& azimuth, double& zenith)
  const
{
  CoordinateSystemPtr coordinateOriginCS;
  /* Yes... for the options "eReference", "eSdReconstruction",
     "eFdReconstruction" the coordinate origin stored in the RRecShower
     is used to calculate azimuth and zenith angle. If "eMcTruth" take the MC core.
     If "eRdReconstruction" take azimuth and zenith stored in RRecShower. */
  if (fUsedDirection == eReference || fUsedDirection == eSdReconstruction ||
      fUsedDirection == eFdReconstruction || fUsedDirection == eLineOfSightStationShowerMaximum) {
    coordinateOriginCS = fwk::LocalCoordinateSystem::Create(
      event.GetRecShower().GetRRecShower().GetCoordinateOrigin());
  }

  // check which direction is to be used (for eIndividualPerStation, look further down)
  if (fUsedDirection == eIndividualPerStation) {
    /* If individual stations are used, then take the shower axis from the radio reconstruction:
       The shower axis will be needed to calculate the magnetic field expectation (Lorentz angle etc.) */
    if (event.GetRecShower().HasRRecShower()) {
      INFOIntermediate("Using shower axis from radio reconstruction, "
                       "but individual arrival directions for each station");
      showerAxis = event.GetRecShower().GetRRecShower().GetAxis();
    } else {
      ERROR("Module is configured to use individual direction, which still take the overall "
        " shower axis as RD direction, but event has no RRecShower!");
      // This may not cause a eFailure. The shower Axis is only needed if number of channels != 2.
    }
  }
  else if (fUsedDirection == eReference || fUsedDirection == eLineOfSightStationShowerMaximum) {
    // See comments in xml.in
    showerAxis = event.GetRecShower().GetRRecShower().GetReferenceAxis(event);
    azimuth = showerAxis.GetPhi(coordinateOriginCS);
    zenith = showerAxis.GetTheta(coordinateOriginCS);
  }
  else if (fUsedDirection == eMcTruth) {
    if (event.HasSimShower()) {
      const CoordinateSystemPtr localCS = event.GetSimShower().GetLocalCoordinateSystem();
      showerAxis = -event.GetSimShower().GetDirection();
      zenith = showerAxis.GetTheta(localCS);
      azimuth = showerAxis.GetPhi(localCS);
    } else {
      ERROR("Module is configured to use Monte Carlo truth direction, but event has no SimShower!");
      return false;
    }
  } else if (event.HasRecShower()) {
    // which axis has to be used?
    if (fUsedDirection == eSdReconstruction) {
      INFOIntermediate("Using direction obtained from SD reconstruction");

      // checking if it has a SRecShower
      if (!event.GetRecShower().HasSRecShower()) {
        ERROR("Module is configured to use SD direction, but event has no SRecShower!");
        return false;
      } else {  // get angles of the SD axis in the coordinate system valid at the RD Coordinate Origin
        showerAxis = event.GetRecShower().GetSRecShower().GetAxis();
        azimuth = showerAxis.GetPhi(coordinateOriginCS);
        zenith = showerAxis.GetTheta(coordinateOriginCS);
      }
    }
    else if (fUsedDirection == eRdReconstruction) {
      if (event.GetRecShower().HasRRecShower()) {
        INFOIntermediate("Using angles from radio reconstruction");
        zenith = event.GetRecShower().GetRRecShower().GetZenith();
        azimuth = event.GetRecShower().GetRRecShower().GetAzimuth();
        showerAxis = event.GetRecShower().GetRRecShower().GetAxis();
      } else {
        ERROR("Module is configured to use RD direction, but event has no RRecShower!");
        return false;
      }
    }
    else if (fUsedDirection == eFdReconstruction) {
      INFOIntermediate("Using direction obtained from FD reconstruction");

      //checking if it has a FRecShower
      if (!event.HasFEvent()) {
        ERROR("Module is configured to use FD direction, but event has no FEvent!");
        return false;
      } else {
        const FEvent& fdEvent = event.GetFEvent();
        for (FEvent::ConstEyeIterator eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
            eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIter) {
          if (eyeIter->GetId() != fWhichEye) continue; // continue until selected eye
          const fevt::Eye& eye = *eyeIter;
          if(!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
            ERROR("Eye has no RecData or FRecShower but FD direction specified in the settings!");
            return false;
          } else {
            // get angles of the FD axis in the coordinate system valid at the RD Coordinate Origin
            showerAxis = eye.GetRecData().GetFRecShower().GetAxis();
            azimuth = showerAxis.GetPhi(coordinateOriginCS);
            zenith = showerAxis.GetTheta(coordinateOriginCS);
          }
        }
      }
    }
  } else {
    ERROR("No reconstructed shower found, but shower direction is required to apply antenna model! "
      "Alternatively, individual signal arrival direction can be defined for each station, "
      "when this feature is activated in the XML file.");
    return false;
  }

  return true;
}

double
RdAntennaChannelToStationConverter::GetXmaxEstimator(const evt::Event& event)
  const
{
  double xmax = 0;

  if (fXmaxEstimator == "MC") {
    if (!event.HasSimShower()) {
      ERROR("Event has no sim shower but was requested.");
      return nan("");
    }
    const evt::ShowerSimData& simShower = event.GetSimShower();
    const auto& gh = simShower.GetGHParameters();
    xmax = gh.GetXMax();
  }
  else {
    if (!event.GetRecShower().HasSRecShower()) {
      ERROR("Event has no SD rec shower. Can not estimate Xmax with utl::XmaxParam::Mean.");
      return nan("");
    }
    const evt::ShowerSRecData& sdShower = event.GetRecShower().GetSRecShower();
    const double crEnergy = sdShower.GetEnergy();
    utl::HadronicInteractionModel hadModel;

    if (fXmaxEstimator == "ParamWithSDEnergy-Sibyll2.3d")
      hadModel = utl::HadronicInteractionModel::eSibyll23d;
    else if (fXmaxEstimator == "ParamWithSDEnergy-EPOS-LHC")
      hadModel = utl::HadronicInteractionModel::eEPOSLHC;
    else if (fXmaxEstimator == "ParamWithSDEnergy-QGSJETII-04")
      hadModel = utl::HadronicInteractionModel::eQGSJETII04;
    else {
      ERROR("Wrong xml configuration!");
      return nan("");
    }

    xmax = 0.5 * (XmaxParam::Mean(crEnergy, 1, hadModel) + XmaxParam::Mean(crEnergy, 56, hadModel));
  }

  return xmax;
}

}