RdPolarizationReconstructor.cc

Go to the documentation of this file.

#include "RdPolarizationReconstructor.h"

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

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/TimeStamp.h>
#include <utl/TimeInterval.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/Math.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/AugerCoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/UTMPoint.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSRecData.h>
#include <revt/Channel.h>

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

#include <fwk/LocalCoordinateSystem.h>

#define OUT(x) if ((x) <= fInfoLevel) cerr << "  "

using namespace revt;
using namespace fwk;
using namespace utl;
using std::string;
using std::pair;
using std::cerr;
using std::endl;


namespace RdPolarizationReconstructor {

  //Complex const I(0.0, 1.0); // To be used as shorthand

  // This is to create a coordinate system at the Rodrigo tank.
  // Unfortunately there is no survey of the earth magnetic field in every possible location
  // so we have to assume that it is constant over the complete area of the observatory.
  CoordinateSystemPtr
  GetRodrigoCS()
  {
    ReferenceEllipsoid wgs84 =
      ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
    const CoordinateSystemPtr ecef = wgs84.GetECEF();

    const UTMPoint rodrigoOrigin(6113924.83*m, 448375.63*m, 1560.0*m, 19, 'H', wgs84);
    return AugerBaseCoordinateSystem::Create(rodrigoOrigin, ecef);
  }


  typedef TraceV3C::SizeType SizeType; // Shorthand


  BoostVecRangeC
  RawStokesParameters::GetSliceX(SizeType pos, SizeType size)
  {
    return BoostVecRangeC(fSliceX, boost::numeric::ublas::range(pos, pos + size));
  }


  BoostVecRangeC
  RawStokesParameters::GetSliceY(SizeType pos, SizeType size)
  {
    return BoostVecRangeC(fSliceY, boost::numeric::ublas::range(pos, pos + size));
  }


  Complex
  RawStokesParameters::GetIndex(char channel, int index)
    const
  {
    switch (channel) {
    case 'X': return fSliceX[index];
    case 'Y': return fSliceY[index];
    case 'A': return fSliceA[index];
    case 'B': return fSliceB[index];
    case 'L': return fSliceL[index];
    case 'R': return fSliceR[index];
    }
    ERROR("The channel parameters should be one of the letters in \"XYABLR\"");
    return Complex(0,0);
  }


  // initializes the polarization data with four channels: two for the signal
  // and two more to correct for the background noise
  StokesParameters::StokesParameters(const BoostVecRangeC& signal_x, const BoostVecRangeC& signal_y,
                                     const BoostVecRangeC& noise_x, const BoostVecRangeC& noise_y,
                                     bool subtract_noise) :
    fSignal(signal_x, signal_y),
    fNoise(noise_x, noise_y),
    fStokesI(subtract_noise ? fSignal.GetRawStokesI() - fNoise.GetRawStokesI() : fSignal.GetRawStokesI()),
    fStokesQ(subtract_noise ? fSignal.GetRawStokesQ() - fNoise.GetRawStokesQ() : fSignal.GetRawStokesQ()),
    fStokesU(subtract_noise ? fSignal.GetRawStokesU() - fNoise.GetRawStokesU() : fSignal.GetRawStokesU()),
    fStokesV(subtract_noise ? fSignal.GetRawStokesV() - fNoise.GetRawStokesV() : fSignal.GetRawStokesV()),
    fStokesRatioQdivI(fStokesQ/fStokesI),
    fStokesRatioUdivI(fStokesU/fStokesI),
    fStokesRatioVdivI(fStokesV/fStokesI)
  { }


  StokesParameters::StokesParameters(const BoostVecC& signal_x, const BoostVecC& signal_y,
                                     RawStokesParameters& noise, double multiply_noise) :
    fSignal(signal_x, signal_y),
    fStokesI(fSignal.GetRawStokesI() - multiply_noise*noise.GetRawStokesI()),
    fStokesQ(fSignal.GetRawStokesQ() - multiply_noise*noise.GetRawStokesQ()),
    fStokesU(fSignal.GetRawStokesU() - multiply_noise*noise.GetRawStokesU()),
    fStokesV(fSignal.GetRawStokesV() - multiply_noise*noise.GetRawStokesV()),
    fStokesRatioQdivI(fStokesQ/fStokesI),
    fStokesRatioUdivI(fStokesU/fStokesI),
    fStokesRatioVdivI(fStokesV/fStokesI)
  { }


  // Part of inner loop in GetNoiseAdditionCovariance
  void
  StokesParameters::NoiseAdditionCovarianceHelper(BoostVecD& quantities, const StationRRecDataQuantities quantity_enum,
                                                  const StokesParameters& sliding_stokes, SizeType index)
  {
    switch (quantity_enum) {
    case eStokesI:
      quantities[index] = sliding_stokes.GetStokesI();
      break;
    case eStokesQ:
      quantities[index] = sliding_stokes.GetStokesQ();
      break;
    case eStokesU:
      quantities[index] = sliding_stokes.GetStokesU();
      break;
    case eStokesV:
      quantities[index] = sliding_stokes.GetStokesV();
      break;
    case eStokesRatioQdivI:
      quantities[index] = sliding_stokes.GetStokesRatioQdivI();
      break;
    case eStokesRatioUdivI:
      quantities[index] = sliding_stokes.GetStokesRatioUdivI();
      break;
    case eStokesRatioVdivI:
      quantities[index] = sliding_stokes.GetStokesRatioVdivI();
      break;
    case ePolarizationAngle:
      quantities[index] = sliding_stokes.GetPolarizationAngle();
      break;
    default:
      ERROR("Incorrect enum.");
      break;
    }
  }


  // It is assumed that the uncertainties are Gaussian. This is only approximately
  // the case for most of the values and especially for the polarization angle this can be a
  // problem when the values wrap through -0.5*pi and 0.5*pi. This wrapping leads to a problem
  // when calculating the mean when the individual values are distributed around -0.5*pi and 0.5*pi.
  // The assumption of Gaussianity becomes acceptable when the uncertainty is sufficiently small
  // but without this subroutine, errors would be made for values close to -0.5*pi and 0.5*pi.
  void
  StokesParameters::NoiseAdditionCovariancePolarizationAngleHelper(double &angle, double origval)
    const
  {
    angle -= origval;
    while (angle > kPiOnTwo)
      angle -= kPi;
    while (angle <= -kPiOnTwo)
      angle += kPi;
    angle += origval;
  }


  double
  StokesParameters::GetNoiseAdditionCovariance(const StationRRecDataQuantities quantity_enum1,
                                               const StationRRecDataQuantities quantity_enum2)
  {
    using namespace boost::numeric::ublas;
    if (fNoise.GetSize() == 0) {
      ERROR("(Co)variances can't be calculated if no noise has been provided.");
    }
    SizeType signal_size = fSignal.GetSize();
    SizeType num_sliding_windows = fNoise.GetSize() - signal_size;
    BoostVecD quantities1(num_sliding_windows);
    BoostVecD quantities2(num_sliding_windows);
    for (SizeType i = 0; i < num_sliding_windows; ++i) {
      const BoostVecRangeC sliding_noise_window_x = fNoise.GetSliceX(i, signal_size);
      const BoostVecRangeC sliding_noise_window_y = fNoise.GetSliceY(i, signal_size);
      const BoostVecRangeC signal_x = fSignal.GetSliceX(0, signal_size);
      const BoostVecRangeC signal_y = fSignal.GetSliceY(0, signal_size);
      BoostVecC double_noise_signal_x = signal_x + sliding_noise_window_x;
      BoostVecC double_noise_signal_y = signal_y + sliding_noise_window_y;
      StokesParameters sliding_stokes(double_noise_signal_x, double_noise_signal_y, fNoise, sqrt(2));
      NoiseAdditionCovarianceHelper(quantities1, quantity_enum1, sliding_stokes, i);
      NoiseAdditionCovarianceHelper(quantities2, quantity_enum2, sliding_stokes, i);
    }

    // Calculate the mean
    double mean1 = 0;
    double mean2 = 0;
    for (SizeType i = 0; i < num_sliding_windows; ++i) {
      double q1i = quantities1[i];
      double q2i = quantities2[i];
      if (quantity_enum1 == ePolarizationAngle) {
        NoiseAdditionCovariancePolarizationAngleHelper(q1i, GetPolarizationAngle());
      }
      if (quantity_enum2 == ePolarizationAngle) {
        NoiseAdditionCovariancePolarizationAngleHelper(q2i, GetPolarizationAngle());
      }
      mean1 += q1i;
      mean2 += q2i;
    }
    mean1 /= num_sliding_windows;
    mean2 /= num_sliding_windows;

    // Calculate the differences with the mean
    for (SizeType i = 0; i < num_sliding_windows; ++i) {
      quantities1[i] -= mean1;
      quantities2[i] -= mean2;
      if (quantity_enum1 == ePolarizationAngle) {
        NoiseAdditionCovariancePolarizationAngleHelper(quantities1[i], 0);
      }
      if (quantity_enum2 == ePolarizationAngle) {
        NoiseAdditionCovariancePolarizationAngleHelper(quantities2[i], 0);
      }
    }

    // Estimate the variance and return
    return inner_prod(quantities1, quantities2) / (num_sliding_windows - 1);
  }


  // This method is about different types of covariances than the NoiseAdditionCovariance*() methods.
  // Here we deal with the covariances between the samples for the analytic method
  Complex
  StokesParameters::AnalyticMethodSampleCovarianceHelper(char channel1, char channel2, SizeType index)
  {
    SizeType Ms = fSignal.GetSize();
    // The map is used to cache already calculated values
    // The count() method is implicitly cast to a boolean value to determine wether the channel is in the map or not
    // Because sampleCovariances is of the type map count() can only return 0 or 1
    // The use of this count method may seem strange but specified as such in the C++ documentation
    if (!sampleCovariances.count(pair<char, char>(channel1, channel2))) {
      SizeType Mn = fNoise.GetSize();
      BoostVecC covariances(Ms);
      for (SizeType i = 0; i < Ms; ++i) {
        covariances[i] = Complex(0, 0);
        for (SizeType j = 0; j < Mn - Ms; ++j) {
          covariances[i] += fNoise.GetIndex(channel2, j) * conj(fNoise.GetIndex(channel1, j+i));
        }
        covariances[i] /= Mn - Ms; // not Mn - Ms - 1 because the baseline has already been substracted for the whole trace that's why we don't subtract the mean either
      }
      // To store the values in the map we only need to give the command below
      // The python equivalent would be something like: sampleCovariances[channel1, channel2] = covariances
      // But of course the next line is much more human readable in C++
      sampleCovariances.insert(
        pair<pair<char, char>, BoostVecC>(
          pair<char, char>(channel1, channel2),
          covariances
        )
      );
    }
    if (index >= Ms) {
      ERROR("Requested index out of range.");
    }
    return sampleCovariances[pair<char, char>(channel1, channel2)][index];
  }


  double
  StokesParameters::AnalyticMethodOneChannelVarianceHelper(char channel)
  {
    // Values are cached in the same way as in StokesParameters::AnalyticMethodSampleCovarianceHelper()
    if (!oneChannelPulseVariances.count(channel)) {
      SizeType Ms = fSignal.GetSize();
      Complex result(0.0, 0.0);
      for (SizeType i = 0; i < Ms; ++i) {
        for (SizeType j = 0; j < Ms; ++j) {
          Complex Vi = fSignal.GetIndex(channel, i);
          Complex Vjc = conj(fSignal.GetIndex(channel, j));
          Complex Covccij = i > j ?
            AnalyticMethodSampleCovarianceHelper(channel, channel, i-j) :
            conj(AnalyticMethodSampleCovarianceHelper(channel, channel, j-i));
          result += Complex(2, 0)*Vi*Vjc*Covccij;
        }
      }
      oneChannelPulseVariances.insert(pair<char, double>(channel, fabs(result.real())/(Ms*Ms)));
    }
    return oneChannelPulseVariances[channel];
  }


  // Currently this variance helper only calculates the variances in Q, U, V, QdivI, UdivI, VdivI
  // and the polarization angle. The covariances between these quantities are not implemented with the
  // analytical method.
  double
  StokesParameters::GetAnalyticMethodVariance(const StationRRecDataQuantities quantity_enum)
  {
    #define VAR(x) AnalyticMethodOneChannelVarianceHelper(x)
    #define I fStokesI
    #define Q fStokesQ
    #define U fStokesU
    #define V fStokesV
    switch (quantity_enum) {
      default:
        ERROR("This enum is not allowed.");
        return std::numeric_limits<double>::infinity();
      case eStokesI: // fall through expected: these four parameters have the same uncertainty (see GAP 2011-088)
      case eStokesQ:
      case eStokesU:
      case eStokesV:
        return VAR('X') + VAR('Y');
      case eStokesRatioQdivI:
        return (VAR('X')*pow((I - Q), 2) + VAR('Y')*pow((I + Q), 2))/pow(I, 4);
      case eStokesRatioUdivI:
        return (VAR('A')*pow((I - U), 2) + VAR('B')*pow((I + U), 2))/pow(I, 4);
      case eStokesRatioVdivI:
        return (VAR('L')*pow((I - V), 2) + VAR('R')*pow((I + V), 2))/pow(I, 4);
      case ePolarizationAngle:
        return (VAR('X') + VAR('Y'))/(4*(Q*Q)*(1 + U*U/(Q*Q)));
    }
    #undef VAR
    #undef I
    #undef Q
    #undef U
    #undef V
  }


  RdPolarizationReconstructor::RdPolarizationReconstructor() :
    fInfoLevel(0),
    fProjectOntoGroundPlane(1),
    fQuadraticNoiseSubtraction(true),
    fGeometrySource("ShowerRRecData"),
    fErrorCalculationType("NoiseAddition"),
    fSignalWindowStartOverride(0.0),
    fSignalWindowStopOverride(0.0),
    fNoiseWindowStartOverride(0.0),
    fNoiseWindowStopOverride(0.0)
  {
    // The geomagnetic field for the simulations
    // At this moment this is hardcoded here.
    // The origin of these values can be found on
    // the wiki page at
    // https://www.auger.unam.mx/AugerWiki/AERA%2C_Polarization_analysis

    // First create a UTM point at the specified location of Rodrigo
    CoordinateSystemPtr rodrigoCs(GetRodrigoCS());

  }


  Vector
  RdPolarizationReconstructor::GetArrivalDirection(const evt::Event& event)
    const
  {
    if (!event.HasRecShower() && fGeometrySource.find("Rec") != string::npos)
      ERROR("Event doesn't have a ShowerRecData object!");
    Vector axis;
    if (fGeometrySource == "ShowerRRecData") {
      if (!event.GetRecShower().HasRRecShower())
        ERROR("Event doesn't have a ShowerRRecData object!");
      axis = event.GetRecShower().GetRRecShower().GetAxis();
    } else if (fGeometrySource == "ShowerSRecData") {
      if (!event.GetRecShower().HasSRecShower())
        ERROR("Event doesn't have a ShowerSRecData object!");
      axis = event.GetRecShower().GetSRecShower().GetAxis();
    } else {
      if (!event.HasSimShower())
        ERROR("Event doesn't have a ShowerSimData object!");
      axis = event.GetSimShower().GetDirection();
    }
    return -axis / axis.GetMag();
  }


  CoordinateSystemPtr
  RdPolarizationReconstructor::GetCS(const evt::Event& event)
    const
  {
    if (!event.HasRecShower() && fGeometrySource.find("Rec") != string::npos)
      ERROR("Event doesn't have a ShowerRecData object!");
    if (fGeometrySource == "ShowerRRecData") {
      if (!event.GetRecShower().HasRRecShower())
        ERROR("Event doesn't have a ShowerRRecData object!");
      return LocalCoordinateSystem::Create(event.GetRecShower().GetRRecShower().GetCoordinateOrigin());
    } else if (fGeometrySource == "ShowerSRecData") {
      if (!event.GetRecShower().HasSRecShower())
        ERROR("Event doesn't have a ShowerSRecData object!");
      return LocalCoordinateSystem::Create(event.GetRecShower().GetSRecShower().GetBarycenter());
    } else {
      if (!event.HasSimShower())
        ERROR("Event doesn't have a ShowerSimData object!");
      return event.GetSimShower().GetLocalCoordinateSystem();
    }
  }


  void
  RdPolarizationReconstructor::GetCoreXY(const evt::Event& event,
                                         const CoordinateSystemPtr& CS, double& CoreX, double& CoreY)
    const
  {
    if (!event.HasRecShower() && fGeometrySource.find("Rec") != string::npos)
      ERROR("Event doesn't have a ShowerRecData object!");
    Point position;
    if (fGeometrySource == "ShowerRRecData") {
      if (!event.GetRecShower().HasRRecShower())
        ERROR("Event doesn't have a ShowerRRecData object!");
      position = event.GetRecShower().GetRRecShower().GetCorePosition();
    } else if (fGeometrySource == "ShowerSRecData") {
      if (!event.GetRecShower().HasSRecShower())
        ERROR("Event doesn't have a ShowerSRecData object!");
      position = event.GetRecShower().GetSRecShower().GetBarycenter();
    } else {
      if (!event.HasSimShower())
        ERROR("Event doesn't have a ShowerSimData object!");
      position = event.GetSimShower().GetPosition();
    }
    CoreX = position.GetX(CS);
    CoreY = position.GetY(CS);
  }


  Point
  RdPolarizationReconstructor::GetCore(const evt::Event& event)
    const
  {
    if (!event.HasRecShower() && fGeometrySource.find("Rec") != string::npos)
      ERROR("Event doesn't have a ShowerRecData object!");
    if (fGeometrySource == "ShowerRRecData") {
      if (!event.GetRecShower().HasRRecShower())
        ERROR("Event doesn't have a ShowerRRecData object!");
      return event.GetRecShower().GetRRecShower().GetCorePosition();
    } else if (fGeometrySource == "ShowerSRecData") {
      if (!event.GetRecShower().HasSRecShower())
        ERROR("Event doesn't have a ShowerSRecData object!");
      return event.GetRecShower().GetSRecShower().GetBarycenter();
    } else {
      if (!event.HasSimShower())
        ERROR("Event doesn't have a ShowerSimData object!");
      return event.GetSimShower().GetPosition();
    }
  }


  Vector
  RdPolarizationReconstructor::GetPerpVector(const Point& point, const Line& line)
    const
  {
    const Vector a = point - line.GetAnchor();
    const Vector& dir = line.GetDirection();
    const Vector b = (dir*a)*dir;
    const Vector perp = a - b;
    return perp;
  }


  VModule::ResultFlag
  RdPolarizationReconstructor::Init()
  {
    INFO("RdPolarizationReconstructor::Init()");

    // Read in the configurations of the xml file
    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdPolarizationReconstructor");
    topBranch.GetChild("QuadraticNoiseSubtraction").GetData(fQuadraticNoiseSubtraction);
    topBranch.GetChild("ProjectOntoGroundPlane").GetData(fProjectOntoGroundPlane);
    topBranch.GetChild("AssumeGroundPlaneIsHorizontalAndFlat").GetData(fAssumeGroundPlaneIsHorizontalAndFlat);
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("GeometrySource").GetData(fGeometrySource);
    topBranch.GetChild("ErrorCalculationType").GetData(fErrorCalculationType);
    topBranch.GetChild("SignalWindowStartOverride").GetData(fSignalWindowStartOverride);
    topBranch.GetChild("SignalWindowStopOverride").GetData(fSignalWindowStopOverride);
    topBranch.GetChild("NoiseWindowStartOverride").GetData(fNoiseWindowStartOverride);
    topBranch.GetChild("NoiseWindowStopOverride").GetData(fNoiseWindowStopOverride);

    return eSuccess;
  }


  VModule::ResultFlag
  RdPolarizationReconstructor::Run(evt::Event& event)
  {
    INFO("RdPolarizationReconstructor::Run()");

    // Check wether there are events at all
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    // The rdet::RDetector object is needed to retrieve the coordinate system that corresponds
    // to the StationTimeSeries instance. These lines have to be removed when the appropriate
    // functionality has been incorporated in the StationTimeSeries.
    REvent& rEvent = event.GetREvent();
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    for (REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin();
         sIt != rEvent.SignalStationsEnd(); ++sIt) {

      // Creating the analytic signal
      StationFFTDataContainer eFieldDataImag = sIt->GetFFTDataContainer();
      FFTDataContainerAlgorithm::HilbertTransform(eFieldDataImag);
      const StationTimeSeries &TimeSeriesReal = sIt->GetFFTDataContainer().GetConstTimeSeries();
      const StationTimeSeries &TimeSeriesImag = eFieldDataImag.GetConstTimeSeries();
      TraceV3C analyticSignal;
      for (unsigned int i = 0; i < TimeSeriesReal.GetSize(); ++i) {
        Vector3C c;
        c[0] = Complex(TimeSeriesReal[i][0], TimeSeriesImag[i][0]);
        c[1] = Complex(TimeSeriesReal[i][1], TimeSeriesImag[i][1]);
        c[2] = Complex(TimeSeriesReal[i][2], TimeSeriesImag[i][2]);
        analyticSignal.PushBack(c);
      }

      SizeType size = analyticSignal.GetSize();
      // the 3D time series need to be translated to a 2D signal that is perpendicular
      // to the Poynting-vector (= arrival-direction) of the pulse
      BoostVecC traceX(size);
      BoostVecC traceY(size);
      double observerAngle;
      double alpha;
      double chxfraction;
      // traceX and traceY will be filled with the rotated complex electric field
      GetRotatedPolarizationTraces(event, *sIt, rDetector, analyticSignal, traceX, traceY,
                                    observerAngle, alpha);
      // stokes parameters are calculated from traceX and traceY using the RecData signal window and noise window
      CalculateAndFillStokesParameters(*sIt, traceX, traceY, observerAngle, alpha);
      // Charge excess relative contribution "a"
      CalculateAndFillChargeExcessFraction(*sIt, alpha, chxfraction);

    }
    return eSuccess;
  }


  // This method generates two traces that are rotated to the appropriate coordinate system
  // for further polarization analysis. An explicit coordinate transformation is necessary
  // because we need to go to a two dimensional system that is perpendicular to the Poynting vector
  // (The Poynting vector is approximated by the arrival direction)
  // NOTE: If in the future something more elaborate than a plane fit is done then this code
  // needs to be modified with a better approximation of the Poynting vector
  void
  RdPolarizationReconstructor::GetRotatedPolarizationTraces(const evt::Event& event, const Station& station,
                                                            const rdet::RDetector& rDetector, const TraceV3C& analyticSignal,
                                                            BoostVecC& traceX, BoostVecC& traceY, double& observerAngle,
                                                            double& alpha)
    const
  {
    CoordinateSystemPtr CS = GetCS(event);
    // The components of the electric field are now in a vector that has a coordinate frame
    // associated to it

    // Get the arrival direction and normalize to make sure it has length one
    // (which should be the case at this moment but with new definitions this
    // may not be the case any more)
    Vector arrivalDirection = GetArrivalDirection(event);

    OUT(3) << "Event: " << event.GetREvent().GetHeader().GetId() << " Station: " << station.GetId() << endl;
    OUT(3) << "arrivalDirection = " << arrivalDirection.GetX(CS) << " " << arrivalDirection.GetY(CS) << " " << arrivalDirection.GetZ(CS) << endl;

    arrivalDirection /= arrivalDirection.GetMag();
    //magnetic field vector
    utl::Vector bfield = event.GetRecShower().GetRRecShower().GetMagneticFieldVector();
    // Calculate the basis vectors for the vxB frame
    // So the X direction is defined as -1*-vxB
    AxialVector baseVectorX = Cross(arrivalDirection, bfield/bfield.GetMag());

    OUT(3) << "bfield = " << bfield.GetX(CS) << " " << bfield.GetY(CS) << " " << bfield.GetZ(CS) << endl;
    OUT(3) << "-vxB (normalized v and B) = " << -baseVectorX.GetX(CS) << " " << -baseVectorX.GetY(CS) << " " << -baseVectorX.GetZ(CS) << endl;

    baseVectorX /= baseVectorX.GetMag(); // baseVectorX now holds vxB (normalized)
    AxialVector baseVectorY = Cross(arrivalDirection, baseVectorX);
    baseVectorY /= baseVectorY.GetMag(); // baseVectorY is chosen perpendicular to she shower axis
    // baseVectorX, baseVectorY and a Z vector along the arrival direction now constitute the
    // vxB coordinate frame

    // If ProjectOntoGroundPlane == 1 then Calculate the pojection of vxB
    // onto the ground plane and produce a different coordinate system
    // where baseVectorX and baseVectorY lie on the ground plane
    if (fProjectOntoGroundPlane) {
      AxialVector projectedVectorX(baseVectorX.GetX(CS), baseVectorX.GetY(CS), 0, CS);
      AxialVector projectedVectorY(-baseVectorX.GetY(CS), baseVectorX.GetX(CS), 0, CS);
      projectedVectorX /= projectedVectorX.GetMag();
      projectedVectorY /= projectedVectorY.GetMag();
      baseVectorX = projectedVectorX;
      baseVectorY = projectedVectorY;
      OUT(3) << "projected -vxB (normalized v and B) = " << -baseVectorX.GetX(CS) << " " << -baseVectorX.GetY(CS) << " " << -baseVectorX.GetZ(CS) << endl;
    }

    // Get the core position
    Point core = GetCore(event);

    // Get the station position
    Point stationPos = rDetector.GetStation(station).GetPosition();

    // Set the station Z position to the same value as the core Z position
    // if AssumeGroundPlaneIsHorizontalAndFlat is 1 in the settings
    if (fAssumeGroundPlaneIsHorizontalAndFlat) {
      stationPos = Point(stationPos.GetX(CS), stationPos.GetY(CS), core.GetZ(CS), CS);
    }

    // Calculate the vector with length 1 that is perpendicular to the shower axis
    // and points to the station
    Vector vectorPointingToStation;
    if (!fProjectOntoGroundPlane) {
      // Most general behavior
      vectorPointingToStation = GetPerpVector(stationPos, Line(core, arrivalDirection));
    } else {
      // Calculations for the projected case
      // In this case vectorPointingToStation is not really a perpendictular vector.
      // It is the vector that points from the shower core to the station
      Vector coreProj(core.GetX(CS), core.GetY(CS), 0, CS);
      Vector stationProj(stationPos.GetX(CS), stationPos.GetY(CS), 0, CS);
      vectorPointingToStation = stationProj - coreProj;
    }
    vectorPointingToStation /= vectorPointingToStation.GetMag();

    // Determine the opening angle between the shower axis and the geomagnetic field
    alpha = Angle(arrivalDirection, bfield/bfield.GetMag());
    // // Use the full coordinate frame to allow the angle to run from 0 to 2pi
    // if (Cross(arrivalDirection, bfield).GetZ(CS) < 0) {
    //   alpha = 2*kPi - alpha;
    // }

    // Determine the observer angle in the vxB frame
    observerAngle = Angle(baseVectorX, vectorPointingToStation);

    OUT(3) << "" << endl;
    OUT(3) << "baseVectorX = " << baseVectorX.GetX(CS) << " " << baseVectorX.GetY(CS) << " " << baseVectorX.GetZ(CS) << endl;
    OUT(3) << "baseVectorY = " << baseVectorY.GetX(CS) << " " << baseVectorY.GetY(CS) << " " << baseVectorY.GetZ(CS) << endl;
    OUT(3) << "vectorPointingToStation = " << vectorPointingToStation.GetX(CS) << " " << vectorPointingToStation.GetY(CS) << " " << vectorPointingToStation.GetZ(CS) << endl;
    OUT(3) << "observerAngle = " << observerAngle << endl;

    // Use the full coordinate frame to allow the angle to run from 0 to 2pi
    if (baseVectorY * vectorPointingToStation < 0) {
      observerAngle = 2*kPi - observerAngle;
    }

    OUT(3) << "observerAngle in [0, 2*Pi) = " << observerAngle << endl;

    for (SizeType i = 0; i != analyticSignal.GetSize(); ++i) {
      Vector3C electricField3C = analyticSignal[i];

      Vector electricFieldVectorRe(electricField3C[0].real(), electricField3C[1].real(), electricField3C[2].real(), CS);
      Vector electricFieldVectorIm(electricField3C[0].imag(), electricField3C[1].imag(), electricField3C[2].imag(), CS);

      double electricFieldXRe = electricFieldVectorRe*baseVectorX;
      double electricFieldYRe = electricFieldVectorRe*baseVectorY;
      double electricFieldXIm = electricFieldVectorIm*baseVectorX;
      double electricFieldYIm = electricFieldVectorIm*baseVectorY;
      // electricFieldX and electricFieldY are now the components perpendicular to the
      // Poynting vector which is approximated by the arrival direction. electricFieldX
      // is parallel to vxB. electricFieldY is perpendicular to vxB and the arrival direction.
      traceX[i] = Complex(electricFieldXRe, electricFieldXIm);
      traceY[i] = Complex(electricFieldYRe, electricFieldYIm);
    }
  }


  void
  RdPolarizationReconstructor::CheckSignalAndNoiseRanges(SizeType signal_start_bin, SizeType signal_stop_bin,
                                                         SizeType noise_start_bin, SizeType noise_stop_bin, SizeType size)
    const
  {
    OUT(2) << "Ranges for extracting polarization information:" << endl;
    OUT(2) << "    - signal range = [" << signal_start_bin << ", " << signal_stop_bin << ") samples" << endl;
    OUT(2) << "    - noise range  = [" << noise_start_bin << ", " << noise_stop_bin << ") samples" << endl;

    if (signal_start_bin > size)
      ERROR("signal_start_bin out of bounds.");
    if (signal_stop_bin > size)
      ERROR("signal_stop_bin out of bounds.");
    if (noise_start_bin > size)
      ERROR("noise_start_bin out of bounds.");
    if (noise_stop_bin > size)
      ERROR("noise_stop_bin out of bounds.");
    if (signal_start_bin >= signal_stop_bin)
      ERROR("invalid values: can't allow: signal_start_bin >= signal_stop_bin");
    if (noise_start_bin >= noise_stop_bin)
      ERROR("invalid values: can't allow: noise_start_bin >= noise_top_bin");
    if (signal_stop_bin - signal_start_bin > noise_stop_bin - noise_start_bin)
      ERROR("the signal region should not be larger than the noise region");
    if (((noise_start_bin >= signal_start_bin) && (noise_start_bin < signal_stop_bin)) ||
        ((noise_stop_bin >= signal_start_bin) && (noise_stop_bin < signal_stop_bin)) ||
        ((noise_start_bin <= signal_start_bin) && (noise_stop_bin > signal_stop_bin)))
      ERROR("the signal and noise regions may not overlap");
    if ((signal_stop_bin - signal_start_bin)*5 > noise_stop_bin - noise_start_bin)
      WARNING("for accurate error determination the signal region should be much smaller than the noise region");
  }


  void
  RdPolarizationReconstructor::CalculateAndFillStokesParameters(Station& station, BoostVecC& traceX,
                                                                BoostVecC& traceY, double observerAngle, double alpha)
  {
    using namespace boost::numeric::ublas;
    double binning = station.GetStationTimeSeries().GetBinning();
    SizeType total_size = station.GetStationTimeSeries().GetSize();
    // Use start and stop times from StationRecData

    double signal_start_time = fSignalWindowStartOverride == 0.0 ?
      station.GetRecData().GetParameter(eSignalWindowStart) : fSignalWindowStartOverride;
    double signal_stop_time =  fSignalWindowStopOverride  == 0.0 ?
      station.GetRecData().GetParameter(eSignalWindowStop) : fSignalWindowStopOverride;
    double noise_start_time =  fNoiseWindowStartOverride  == 0.0 ?
      station.GetRecData().GetParameter(eNoiseWindowStart) : fNoiseWindowStartOverride;
    double noise_stop_time =   fNoiseWindowStopOverride   == 0.0 ?
      station.GetRecData().GetParameter(eNoiseWindowStop) : fNoiseWindowStopOverride;
    // translate the times to bins

    OUT(3) << "signal_start_time   = " << signal_start_time/ns    << " ns" << endl;
    OUT(3) << "signal_stop_time    = " << signal_stop_time/ns     << " ns" << endl;
    OUT(3) << "noise_start_time    = " << noise_start_time/ns     << " ns" << endl;
    OUT(3) << "noise_stop_time     = " << noise_stop_time/ns      << " ns" << endl;
    OUT(3) << "binning             = " << binning/ns              << " ns" << endl;

    SizeType signal_start_bin = round(signal_start_time / binning);
    SizeType signal_stop_bin = round(signal_stop_time / binning);
    SizeType noise_start_bin = round(noise_start_time / binning);
    SizeType noise_stop_bin = round(noise_stop_time / binning);
    // test wether the start and stop times contain decent values
    CheckSignalAndNoiseRanges(signal_start_bin, signal_stop_bin, noise_start_bin, noise_stop_bin, total_size);
    // cut out the signal and the noise region from the trace
    BoostVecRangeC signal_x(traceX, range(signal_start_bin, signal_stop_bin));
    BoostVecRangeC signal_y(traceY, range(signal_start_bin, signal_stop_bin));
    BoostVecRangeC noise_x(traceX, range(noise_start_bin, noise_stop_bin));
    BoostVecRangeC noise_y(traceY, range(noise_start_bin, noise_stop_bin));

    // here the stokes parameters are calculated
    StokesParameters stokes(signal_x, signal_y, noise_x, noise_y, fQuadraticNoiseSubtraction);
    // the stokes parameters are filled into the reconstruction data structure here
    StationRecData &recData(station.GetRecData());
    recData.SetParameter(eStokesI, stokes.GetStokesI());
    recData.SetParameter(eStokesQ, stokes.GetStokesQ());
    recData.SetParameter(eStokesU, stokes.GetStokesU());
    recData.SetParameter(eStokesV, stokes.GetStokesV());
    recData.SetParameter(eStokesRatioQdivI, stokes.GetStokesRatioQdivI());
    recData.SetParameter(eStokesRatioUdivI, stokes.GetStokesRatioUdivI());
    recData.SetParameter(eStokesRatioVdivI, stokes.GetStokesRatioVdivI());
    recData.SetParameter(ePolarizationAngle, stokes.GetPolarizationAngle());

    recData.SetParameter(eObserverAngle, observerAngle);
    recData.SetParameter(eAlpha, alpha);

    // the errors, i.e. variances and covariances for the stokes parameters are filled here
    // if requested by the configuration
    StationRRecDataQuantities quantities1[8] = {
      eStokesI,
      eStokesQ,
      eStokesU,
      eStokesV,
      eStokesRatioQdivI,
      eStokesRatioUdivI,
      eStokesRatioVdivI,
      ePolarizationAngle
    };
    if (fErrorCalculationType == "NoiseAddition") {
      for (int i = 0; i < 8; ++i) {
        for (int j = i; j < 8; ++j) {
          StationRRecDataQuantities ei = quantities1[i];
          StationRRecDataQuantities ej = quantities1[j];
          recData.SetParameterCovariance(ei, ej, stokes.GetNoiseAdditionCovariance(ei, ej));
        }
      }
    }
    // The analytical calculation only sets the variances, no covariances are calculated
    if (fErrorCalculationType == "Analytical") {
      for (int i = 0; i < 8; ++i) {
        StationRRecDataQuantities ei = quantities1[i];
        recData.SetParameterCovariance(ei, ei, stokes.GetAnalyticMethodVariance(ei));
      }
    }
  }


  void
  RdPolarizationReconstructor::CalculateAndFillChargeExcessFraction(Station& station, double& alpha, double& chxfraction)
  {
    // Get the polarization Angle and Observer Angle
    double polAngle = station.GetRecData().GetParameter(ePolarizationAngle);
    double obsAngle = station.GetRecData().GetParameter(eObserverAngle);

    // Calculate "a" using polarization angle as function of the observer angle
    chxfraction = (sin(alpha) * tan(polAngle)) / (sin(obsAngle) - cos(obsAngle)*tan(polAngle));

    // Estimate lower-limit uncertainty, approximation but agrees with MC though
    double polAngErr = station.GetRecData().GetParameterError(ePolarizationAngle);
    double chxFracError = abs((sin(alpha)*tan(polAngle + polAngErr))/(sin(obsAngle) - cos(obsAngle)*tan(polAngle + polAngErr)) - (sin(alpha)*tan(polAngle - polAngErr))/(sin(obsAngle) - cos(obsAngle)*tan(polAngle - polAngErr)))/2;

    // Set Parameters
    StationRecData &recData(station.GetRecData());
    recData.SetParameter(eCxxFraction, chxfraction);
    recData.SetParameterError(eCxxFraction, chxFracError);
  }


  VModule::ResultFlag
  RdPolarizationReconstructor::Finish()
  {
    INFO("RdPolarizationReconstructor::Finish()");

    return eSuccess;
  }

}