RdStationPolarizationRejector.cc

Go to the documentation of this file.

#include "RdStationPolarizationRejector.h"

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

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

#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/Vector.h>
#include <utl/CoordinateSystem.h>
#include <utl/AnalyticCoordinateTransformator.h>

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

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

#include <TMath.h>
#include <TVector3.h>

using namespace std;

namespace RdStationPolarizationRejector {

  fwk::VModule::ResultFlag
  RdStationPolarizationRejector::Init()
  {
    const auto topB = fwk::CentralConfig::GetInstance()->GetTopBranch("RdStationPolarizationRejector");
    topB.GetChild("maxBeta").GetData(fMaxBeta);
    topB.GetChild("constantChargeExcess").GetData(fConstantChargeExcess);
    topB.GetChild("maxChargeExcess").GetData(fMaxChargeExcess);
    topB.GetChild("minChargeExcess").GetData(fMinChargeExcess);
    topB.GetChild("percentile").GetData(fProbability);
    topB.GetChild("polStandardDeviations").GetData(fMaxSigmaBeta);
    topB.GetChild("InfoLevel").GetData(fInfoLevel);
    topB.GetChild("zenithCut").GetData(fZenithCut);
    topB.GetChild("noiseErrorModel").GetData(fNoiseErrorModel);
    SetProbability(GetProbability());
    INFODebug("INIT POL REJECT");
    return eSuccess;
  }

  fwk::VModule::ResultFlag
  RdStationPolarizationRejector::Run(evt::Event& event)
  {
    INFODebug("RUN POL REJECT");
    // nothing to do if event has no REvent
    if (!event.HasREvent()) {
      INFODebug("No REvent");
      return eContinueLoop;
    }

    // SD data is required
    if (!event.HasRecShower()) {
      if (fInfoLevel > 0) {
        WARNING("Shower has to be reconstructed before RdStationPolarizationRejector is used... aborting");
      }
      return eBreakLoop;
    }

    const det::Detector& det = det::Detector::GetInstance();
    const rdet::RDetector& rdet = det.GetRDetector();
    const utl::CoordinateSystemPtr referenceCS = det.GetReferenceCoordinateSystem();
    const utl::TimeStamp timeStamp = event.GetHeader().GetTime();

    revt::REvent& rEvent = event.GetREvent();
    const evt::ShowerRRecData& rShower = event.GetRecShower().GetRRecShower();
    const utl::Point corePosition = rShower.GetReferenceCorePosition(event);
    const utl::CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(corePosition);

    double covarianceMatrix[3];
    GetCovarianceMatrix(event, covarianceMatrix, rShower, localCS);

    const utl::Vector showerAxis = rShower.GetReferenceAxis(event);
    const double zenith = showerAxis.GetTheta(localCS);
    const double azimuth = showerAxis.GetPhi(localCS);

    if (fZenithCut) {
      if (!CheckZenith(zenith, azimuth)) {
        stringstream sstr;
        sstr << "Incoming direction zenith = " << zenith / utl::deg
             << " azimuth = " << azimuth / utl::deg
             << ". Direction does not allow for polarization rejection";
        INFODebug(sstr.str());
        return eSuccess;
      }
    }

    const fwk::MagneticFieldModel& magneticFieldModel = fwk::MagneticFieldModel::instance();
    const utl::Vector magneticField = magneticFieldModel.GetMagneticFieldVector(corePosition, timeStamp);
    const double angleToMagneticField = utl::Angle(magneticField, showerAxis);
    utl::Point stationPosition;
    const utl::RadioGeometryUtilities rdGeometryUtilities = utl::RadioGeometryUtilities(showerAxis, localCS,
                                                                                        magneticField);
    for (revt::REvent::SignalStationIterator sIt = rEvent.SignalStationsBegin(); sIt != rEvent.SignalStationsEnd();
         ++sIt) {
      revt::Station& station = *sIt;
      stationPosition = rdet.GetStation(station).GetPosition();

      if (!station.HasSignal() || station.IsSaturated()) {
        continue; // nothing to do here
      }

      revt::StationRecData& stationRecData = station.GetRecData();
      utl::Vector stationRecEfieldVector = utl::Vector(stationRecData.GetParameter(revt::eEFieldVectorEW),
                                                       stationRecData.GetParameter(revt::eEFieldVectorNS),
                                                       stationRecData.GetParameter(revt::eEFieldVectorV),
                                                       localCS);
      double angleToEfieldExpectation = rdGeometryUtilities.GetAngleToEFieldExpectation(stationRecEfieldVector,
                                                                                        corePosition, showerAxis,
                                                                                        stationPosition,
                                                                                        magneticField, 0.14);
      if (angleToEfieldExpectation > (utl::kPi / 2) * utl::radian) {
        stationRecEfieldVector = -stationRecEfieldVector;
        angleToEfieldExpectation = utl::kPi * utl::radian - angleToEfieldExpectation;
      }

      stringstream sstr;
      sstr.str("");
      sstr << "Event: " << event.GetHeader().GetId() << " Station: " << station.GetId();
      INFODebug(sstr.str());

      const double peakAmpMag = stationRecData.GetParameter(revt::ePeakAmplitudeMag);
      const double noiseRmsMag = stationRecData.GetParameter(revt::eNoiseRmsMag);
      const double snr = peakAmpMag * peakAmpMag / noiseRmsMag / noiseRmsMag;
      const double distToShowerAxis = rdGeometryUtilities.GetDistanceToAxis(showerAxis, corePosition, stationPosition);
      utl::Point semiMajorAxis, semiMinorAxis;
      GetCoreErrorEllipse(semiMajorAxis, semiMinorAxis, corePosition, covarianceMatrix, localCS);

      sstr.str("");
      sstr << "Covariance Matrix: [" << covarianceMatrix[0] << ", " << covarianceMatrix[1] << ", "
           << covarianceMatrix[2] << "]"
           << ", semi major axis: [" << semiMajorAxis.GetX(localCS) << ", " << semiMajorAxis.GetY(localCS) << "]"
           << ", semi minor axis: [" << semiMinorAxis.GetX(localCS) << ", " << semiMinorAxis.GetY(localCS) << "]";
      INFODebug(sstr.str());

      TVector3 semiMajorAxisvxBvxvxB, semiMinorAxisvxBvxvxB;
      double rotationAngle;
      GetCoreErrorEllipsevxBvxvxB(semiMajorAxisvxBvxvxB, semiMinorAxisvxBvxvxB, rotationAngle, semiMajorAxis,
                                    semiMinorAxis, corePosition, rdGeometryUtilities);

      sstr.str("");
      sstr << "semi minor axis vxb-vxvxb ["
           << semiMajorAxisvxBvxvxB.X() << ", " << semiMajorAxisvxBvxvxB.Y() << "]"
           << ", semi minor axis vxb-vxvxb ["
           << semiMinorAxisvxBvxvxB.X() << ", " << semiMinorAxisvxBvxvxB.Y() << "]";
      INFODebug(sstr.str());

      const double absSemiMajorAxisvxBvxvxB = sqrt(utl::Sqr(semiMajorAxisvxBvxvxB.X()) +
                                                   utl::Sqr(semiMajorAxisvxBvxvxB.Y()));

      const double maxChargeExcess = GetMaxChargeExcess(zenith, distToShowerAxis, absSemiMajorAxisvxBvxvxB);
      const double minChargeExcess = GetMinChargeExcess(zenith, distToShowerAxis, absSemiMajorAxisvxBvxvxB);

      const TVector3 EfieldInShowerPlane = GetEfieldInShowerPlane(stationRecEfieldVector, localCS, rdGeometryUtilities);
      const double noiseVxB = stationRecData.GetParameter(revt::eNoiseRmsVxB);
      const double noiseVxVxB = stationRecData.GetParameter(revt::eNoiseRmsVxVxB);
      const double polarizationUncertainty = GetPolAngleUncertainty(snr, EfieldInShowerPlane, noiseVxB, noiseVxVxB);

      sstr.str("");
      sstr << "Angle to EField Expectation: " << angleToEfieldExpectation / utl::degree
           << ", Angle tolerance based on uncertainty: " << fMaxSigmaBeta * polarizationUncertainty / utl::degree
           << ", sin(alpha): " << sin(angleToMagneticField)
           << ", Max charge excess: " << maxChargeExcess
           << ", Min charge excess: " << minChargeExcess;
      INFODebug(sstr.str());

      if ((angleToEfieldExpectation > fMaxBeta) && (sin(angleToMagneticField) > maxChargeExcess) &&
          (angleToEfieldExpectation > fMaxSigmaBeta * polarizationUncertainty)) {
        TVector3 stationPositionvxBvxvxB = GetStationPositionvxBvxvxB(stationPosition, corePosition,
                                                                      rdGeometryUtilities);

        sstr.str("");
        sstr << "pol reconstruction"
             << stationPositionvxBvxvxB.X() << " station x"
             << stationPositionvxBvxvxB.Y() << " station y";
        INFODebug(sstr.str());

        double localPol1, localPol2;
        GetLocalPolMaxima(localPol1, localPol2, angleToMagneticField, maxChargeExcess);

        const double measuredEfieldPol = atan2(EfieldInShowerPlane.Y(), -EfieldInShowerPlane.X());
        if (IsStationInErrorEllipse(stationPositionvxBvxvxB, semiMajorAxisvxBvxvxB, semiMinorAxisvxBvxvxB,
                                      rotationAngle)) {
          const double maxPol = std::max(localPol1, localPol2);
          const double minPol = std::min(localPol1, localPol2);

          sstr.str("");
          sstr << "station inside ellipse: measured Efield-angle: " << measuredEfieldPol / utl::degree
               << ", max Efield-angle: " << maxPol / utl::degree
               << ", min Efield-angle: " << minPol / utl::degree;
          INFODebug(sstr.str());

          if ((maxPol + fMaxSigmaBeta * polarizationUncertainty < measuredEfieldPol) ||
              (minPol - fMaxSigmaBeta * polarizationUncertainty > measuredEfieldPol)) {
            sstr.str("");
            sstr << "rejecting station: " << station.GetId();
            INFOFinal(sstr.str());
            station.SetRejectedReason(revt::ePolarizationDeselected);
          }
        } else {
          TVector3 tangentPoint1, tangentPoint2;
          GetErrorEllipseTangentPoints(tangentPoint1, tangentPoint2, stationPositionvxBvxvxB, semiMajorAxisvxBvxvxB,
                                       semiMinorAxisvxBvxvxB, rotationAngle);

          const double tangentPolMaxChargeExcess1 = GetTangentPol(stationPositionvxBvxvxB, tangentPoint1, maxChargeExcess,
                                                                  angleToMagneticField);
          const double tangentPolMinChargeExcess1 = GetTangentPol(stationPositionvxBvxvxB, tangentPoint1, minChargeExcess,
                                                                  angleToMagneticField);
          const double tangentPolMaxChargeExcess2 = GetTangentPol(stationPositionvxBvxvxB, tangentPoint2, maxChargeExcess,
                                                                  angleToMagneticField);
          const double tangentPolMinChargeExcess2 = GetTangentPol(stationPositionvxBvxvxB, tangentPoint2, minChargeExcess,
                                                                  angleToMagneticField);

          double maxPol =
            std::max(std::max(tangentPolMaxChargeExcess1, tangentPolMinChargeExcess1),
                     std::max(tangentPolMaxChargeExcess2, tangentPolMinChargeExcess2));
          double minPol =
            std::min(std::min(tangentPolMaxChargeExcess1, tangentPolMinChargeExcess1),
                     std::min(tangentPolMaxChargeExcess2, tangentPolMinChargeExcess2));

          const double ellipseTanAngle1 = atan2(tangentPoint1.Y() - stationPositionvxBvxvxB.Y(),
                                          tangentPoint1.X() - stationPositionvxBvxvxB.X());
          const double ellipseTanAngle2 = atan2(tangentPoint2.Y() - stationPositionvxBvxvxB.Y(),
                                          tangentPoint2.X() - stationPositionvxBvxvxB.X());

          sstr.str("");
          sstr << "station outside of ellipse. Tangent points: "
               << "[" << tangentPoint1.X() << ", " << tangentPoint1.Y() << "]"
               << ", [" << tangentPoint2.X() << ", " << tangentPoint2.Y() << "]"
               << ", measured Efield-angle: " << measuredEfieldPol / utl::degree
               << ", max Efield-angle:  " << maxPol / utl::degree << "\n"
               << ", min Efield-angle:  " << minPol / utl::degree << "\n";
          INFODebug(sstr.str());

          if ((std::max(ellipseTanAngle1, ellipseTanAngle2) > acos(maxChargeExcess / sin(angleToMagneticField))) &&
              (std::min(ellipseTanAngle1, ellipseTanAngle2) < acos(maxChargeExcess / sin(angleToMagneticField)))) {
            maxPol = std::max(maxPol, localPol1);
            minPol = std::min(minPol, localPol1);
          }

          if ((std::max(ellipseTanAngle1, ellipseTanAngle2) > -acos(maxChargeExcess / sin(angleToMagneticField))) &&
              (std::min(ellipseTanAngle1, ellipseTanAngle2) < -acos(maxChargeExcess / sin(angleToMagneticField)))) {
            maxPol = std::max(maxPol, localPol2);
            minPol = std::min(minPol, localPol2);
          }

          if ((maxPol + fMaxSigmaBeta * polarizationUncertainty < measuredEfieldPol) ||
              (minPol - fMaxSigmaBeta * polarizationUncertainty > measuredEfieldPol)) {
            sstr.str("");
            sstr << "rejecting station: " << station.GetId();
            INFOFinal(sstr.str());
            station.SetRejectedReason(revt::ePolarizationDeselected);
          }
        }
      }
    }

    return eSuccess;
  }

  /* Returns maximum acceptable charge excess based on zenith and distance to shower axis
     if constantChargeExcess is set to true in steering file, the value specified there is returned */
  double
  RdStationPolarizationRejector::GetMaxChargeExcess(const double zenith,
                                                    const double distToAxis,
                                                    const double semiMajorAxisLength)
  {
    if (fConstantChargeExcess) {
      return fMaxChargeExcess;
    } else {
      int minIndex = 0;
      int maxIndex = 0;
      double minDistToAxis = distToAxis - semiMajorAxisLength;
      double maxDistToAxis = distToAxis + semiMajorAxisLength;
      while (minDistToAxis > 12.5) {
        minIndex++;
        minDistToAxis = minDistToAxis - 25.0;
      }

      while (maxDistToAxis > 12.5) {
        maxIndex++;
        maxDistToAxis = maxDistToAxis - 25.0;
      }

      double maxChargeExcess = 0;
      double currentChargeExcess;
      while (minIndex <= maxIndex) {
        if (zenith / utl::degree < 35.0) {
          currentChargeExcess = fMaxChargeExcess30[minIndex];
        } else if (zenith / utl::degree < 45.0) {
          currentChargeExcess = fMaxChargeExcess40[minIndex];
        } else if (zenith / utl::degree < 55.0) {
          currentChargeExcess = fMaxChargeExcess50[minIndex];
        } else if (zenith / utl::degree < 65.0) {
          currentChargeExcess = fMaxChargeExcess60[minIndex];
        } else {
          return 0.25;
        }

        if (currentChargeExcess > maxChargeExcess) {
          maxChargeExcess = currentChargeExcess;
        }

        minIndex++;
      }

      return maxChargeExcess;
    }
  }

  // Behaves the same way as GetMaxChargeExcess but return the minimum acceptable charge excess
  double
  RdStationPolarizationRejector::GetMinChargeExcess(const double zenith,
                                                    const double distToAxis,
                                                    const double semiMajorAxisLength)
  {
    if (fConstantChargeExcess) {
      return fMinChargeExcess;
    } else {
      int minIndex = 0;
      int maxIndex = 0;
      double minDistToAxis = distToAxis - semiMajorAxisLength;
      double maxDistToAxis = distToAxis + semiMajorAxisLength;
      while (minDistToAxis > 12.5) {
        minIndex++;
        minDistToAxis = minDistToAxis - 25.0;
      }

      while (maxDistToAxis > 12.5) {
        maxIndex++;
        maxDistToAxis = maxDistToAxis - 25.0;
      }

      double minChargeExcess = 1;
      double currentChargeExcess;
      while (minIndex <= maxIndex) {
        if (zenith / utl::degree < 35.0) {
          currentChargeExcess = fMinChargeExcess30[minIndex];
        } else if (zenith / utl::degree < 45.0) {
          currentChargeExcess = fMinChargeExcess40[minIndex];
        } else if (zenith / utl::degree < 55.0) {
          currentChargeExcess = fMinChargeExcess50[minIndex];
        } else if (zenith / utl::degree < 65.0) {
          currentChargeExcess = fMinChargeExcess60[minIndex];
        } else {
          return 0;
        }

        if (currentChargeExcess < minChargeExcess) {
          minChargeExcess = currentChargeExcess;
        }

        minIndex++;
      }

      return minChargeExcess;
    }
  }

  /* Checks if zenith and azimuth are within a range where polarization rejection works relaibly.
     Returns true is pol. rejection can be performed and false otherwise */
  bool
  RdStationPolarizationRejector::CheckZenith(const double zenith, double azimuth)
  {
    int azimuth_index = 0;
    if (azimuth < 0) {
      azimuth = azimuth + 360. * utl::degree;
    }

    while (azimuth > 30 * utl::degree) {
      azimuth_index++;
      azimuth -= 30. * utl::degree;
    }

    if (fMaxZenith[azimuth_index] > zenith) {
      return true;
    } else {
      return false;
    }
  }

  /* Sets the covariance matrix of the SD core position uncertainty. However,
     since the two non-diagonal entries are equal, only one of them is saved */
  void
  RdStationPolarizationRejector::GetCovarianceMatrix(const evt::Event& event,
                                                     double covarianceMatrix[3],
                                                     const evt::ShowerRRecData& rShower,
                                                     const utl::CoordinateSystemPtr& cs)
  {
    const utl::Vector coreError = rShower.GetReferenceCoreError(event);
    const double correlation = rShower.GetReferenceCoreErrorCorrelationXY(event);
    const double sigmaX = coreError.GetX(cs);
    const double sigmaY = coreError.GetY(cs);
    covarianceMatrix[0] = sigmaX * sigmaX;
    covarianceMatrix[1] = sigmaX * sigmaY * correlation;
    covarianceMatrix[2] = sigmaY * sigmaY;
  }

  void
  RdStationPolarizationRejector::SetProbability(const double prob)
  {
    fProbability = prob;
    fChiSquared = TMath::ChisquareQuantile(prob, 2);
  }

  /* Returns the end points of the semi major and semi minor axis of the error ellipse of the SD core position
     Both axes are calculated relative to the a priori SD core position */
  void
  RdStationPolarizationRejector::GetCoreErrorEllipse(utl::Point& semiMajorAxis,
                                                     utl::Point& semiMinorAxis,
                                                     const utl::Point& corePosition,
                                                     const double covarianceMatrix[3],
                                                     const utl::CoordinateSystemPtr& cs)
  {
    double coreErrorEllipseSemiMajorAxis[2];
    double coreErrorEllipseSemiMinorAxis[2];
    if (covarianceMatrix[1] == 0) {
      if (abs(covarianceMatrix[0]) > abs(covarianceMatrix[2])) {
        coreErrorEllipseSemiMajorAxis[0] = sqrt(covarianceMatrix[0] * fChiSquared);
        coreErrorEllipseSemiMajorAxis[1] = 0;
        coreErrorEllipseSemiMinorAxis[0] = 0;
        coreErrorEllipseSemiMinorAxis[1] = sqrt(covarianceMatrix[2] * fChiSquared);
      } else {
        coreErrorEllipseSemiMinorAxis[0] = sqrt(covarianceMatrix[0] * fChiSquared);
        coreErrorEllipseSemiMinorAxis[1] = 0;
        coreErrorEllipseSemiMajorAxis[0] = 0;
        coreErrorEllipseSemiMajorAxis[1] = sqrt(covarianceMatrix[2] * fChiSquared);
      }
    } else {
      double trace = covarianceMatrix[0] + covarianceMatrix[2];
      double determinant = covarianceMatrix[0] * covarianceMatrix[2] - covarianceMatrix[1] * covarianceMatrix[1];
      double eigenValue1 = trace / 2 + sqrt(trace * trace / 4 - determinant);
      double eigenValue2 = trace / 2 - sqrt(trace * trace / 4 - determinant);
      if (abs(eigenValue1) > abs(eigenValue2)) {
        double eigenVector1Norm = sqrt(
          (eigenValue1 - covarianceMatrix[2]) * (eigenValue1 - covarianceMatrix[2]) + covarianceMatrix[1] *
          covarianceMatrix[1]);
        coreErrorEllipseSemiMajorAxis[0] = sqrt(eigenValue1 * fChiSquared) * (eigenValue1 - covarianceMatrix[2]) /
                                           eigenVector1Norm;
        coreErrorEllipseSemiMajorAxis[1] = sqrt(eigenValue1 * fChiSquared) * covarianceMatrix[1] / eigenVector1Norm;
        double eigenVector2Norm = sqrt(
          (eigenValue2 - covarianceMatrix[2]) * (eigenValue2 - covarianceMatrix[2]) + covarianceMatrix[1] *
          covarianceMatrix[1]);
        coreErrorEllipseSemiMinorAxis[0] = sqrt(eigenValue2 * fChiSquared) * (eigenValue2 - covarianceMatrix[2]) /
                                           eigenVector2Norm;
        coreErrorEllipseSemiMinorAxis[1] = sqrt(eigenValue2 * fChiSquared) * covarianceMatrix[1] / eigenVector2Norm;
      } else {
        double eigenVector1Norm = sqrt(
          (eigenValue1 - covarianceMatrix[2]) * (eigenValue1 - covarianceMatrix[2]) + covarianceMatrix[1] *
          covarianceMatrix[1]);
        coreErrorEllipseSemiMinorAxis[0] = sqrt(eigenValue1 * fChiSquared) * (eigenValue1 - covarianceMatrix[2]) /
                                           eigenVector1Norm;
        coreErrorEllipseSemiMinorAxis[1] = sqrt(eigenValue1 * fChiSquared) * covarianceMatrix[1] / eigenVector1Norm;
        double eigenVector2Norm = sqrt(
          (eigenValue2 - covarianceMatrix[2]) * (eigenValue2 - covarianceMatrix[2]) + covarianceMatrix[1] *
          covarianceMatrix[1]);
        coreErrorEllipseSemiMajorAxis[0] = sqrt(eigenValue2 * fChiSquared) * (eigenValue2 - covarianceMatrix[2]) /
                                           eigenVector2Norm;
        coreErrorEllipseSemiMajorAxis[1] = sqrt(eigenValue2 * fChiSquared) * covarianceMatrix[1] / eigenVector2Norm;
      }
    }

    semiMajorAxis = utl::Point(coreErrorEllipseSemiMajorAxis[0] + corePosition.GetX(
                                 cs), coreErrorEllipseSemiMajorAxis[1] + corePosition.GetY(cs), corePosition.GetZ(
                                 cs), cs);
    semiMinorAxis = utl::Point(coreErrorEllipseSemiMinorAxis[0] + corePosition.GetX(
                                 cs), coreErrorEllipseSemiMinorAxis[1] + corePosition.GetY(cs), corePosition.GetZ(
                                 cs), cs);
  }

  /* Transforms the SD core position error ellipse into the vxB-vxvxB-system and projects it onto the shower plane
     Returns the semi major and semi minor axis of the resulting ellipse */
  void
  RdStationPolarizationRejector::GetCoreErrorEllipsevxBvxvxB(TVector3& semiMajorAxisvxBvxvxB,
                                                             TVector3& semiMinorAxisvxBvxvxB,
                                                             double& rotationAngleInVxB,
                                                             const utl::Point& semiMajorAxis,
                                                             const utl::Point& semiMinorAxis,
                                                             const utl::Point& corePosition,
                                                             const utl::RadioGeometryUtilities& rdGeometryUtilities)
  {
    double ax, ay, az;
    double bx, by, bz;
    double coreX, coreY, coreZ;

    rdGeometryUtilities.GetVectorInShowerPlaneVxB(ax, ay, az, semiMajorAxis);
    rdGeometryUtilities.GetVectorInShowerPlaneVxB(bx, by, bz, semiMinorAxis);
    rdGeometryUtilities.GetVectorInShowerPlaneVxB(coreX, coreY, coreZ, corePosition);

    ax = ax - coreX;
    ay = ay - coreY;
    az = az - coreZ;
    bx = bx - coreX;
    by = by - coreY;
    bz = bz - coreZ;

    double aDotB = ax * bx + ay * by;
    if (aDotB == 0) {
      if (ax * ax + ay * ay > bx * bx + by * by) {
        semiMajorAxisvxBvxvxB = TVector3(ax, ay, 0);
        semiMinorAxisvxBvxvxB = TVector3(bx, by, 0);
      } else {
        semiMinorAxisvxBvxvxB = TVector3(ax, ay, 0);
        semiMajorAxisvxBvxvxB = TVector3(bx, by, 0);
      }
    } else {
      const double aSquared = ax * ax + ay * ay;
      const double bSquared = bx * bx + by * by;
      const double tanPhi1 = -0.5 * (aSquared - bSquared) / aDotB +
                             sqrt(0.25 * ((aSquared - bSquared) * (aSquared - bSquared) / aDotB / aDotB) + 1);
      const double tanPhi2 = -0.5 * (aSquared - bSquared) / aDotB -
                             sqrt(0.25 * ((aSquared - bSquared) * (aSquared - bSquared) / aDotB / aDotB) + 1);
      const double phi1 = atan(tanPhi1);
      const double phi2 = atan(tanPhi2);
      const TVector3 axis1 = TVector3(ax * cos(phi1) + bx * sin(phi1), ay * cos(phi1) + by * sin(phi1), 0);
      const TVector3 axis2 = TVector3(ax * cos(phi2) + bx * sin(phi2), ay * cos(phi2) + by * sin(phi2), 0);
      if (axis1.X() * axis1.X() + axis1.Y() * axis1.Y() > axis2.X() * axis2.X() + axis2.Y() * axis2.Y()) {
        semiMajorAxisvxBvxvxB = axis1;
        semiMinorAxisvxBvxvxB = axis2;
        rotationAngleInVxB = atan2(axis1.Y(), axis1.X());
      } else {
        semiMajorAxisvxBvxvxB = axis2;
        semiMinorAxisvxBvxvxB = axis1;
        rotationAngleInVxB = atan2(axis2.Y(), axis2.X());
      }

      if (rotationAngleInVxB > utl::kPi / 2.0) {
        rotationAngleInVxB = rotationAngleInVxB - utl::kPi;
      } else if (rotationAngleInVxB < -utl::kPi / 2.0) {
        rotationAngleInVxB = rotationAngleInVxB + utl::kPi;
      }
    }
  }

  TVector3
  RdStationPolarizationRejector::GetStationPositionvxBvxvxB(const utl::Point stationPosition,
                                                            const utl::Point& corePosition,
                                                            const utl::RadioGeometryUtilities& rdGeometryUtilities)
  {
    double x, y, z;
    double coreX, coreY, coreZ;
    rdGeometryUtilities.GetVectorInShowerPlaneVxB(x, y, z, stationPosition);
    rdGeometryUtilities.GetVectorInShowerPlaneVxB(coreX, coreY, coreZ, corePosition);
    x = x - coreX;
    y = y - coreY;
    z = z - coreZ;
    return TVector3(x, y, z);
  }

  bool
  RdStationPolarizationRejector::IsStationInErrorEllipse(const TVector3& stationPositionvxBvxvxB,
                                                         const TVector3& semiMajorAxisVxB,
                                                         const TVector3& semiMinorAxisVxB,
                                                         const double rotationAngle)
  {
    TVector3 rotatedPosition = stationPositionvxBvxvxB;
    rotatedPosition.RotateZ(-rotationAngle);
    const double xStation = rotatedPosition.X();
    const double yStation = rotatedPosition.Y();
    if (xStation * xStation /
        (semiMajorAxisVxB.X() * semiMajorAxisVxB.X() + semiMajorAxisVxB.Y() * semiMajorAxisVxB.Y()) + yStation *
        yStation /
        (semiMinorAxisVxB.X() * semiMinorAxisVxB.X() + semiMinorAxisVxB.Y() * semiMinorAxisVxB.Y()) > 1) {
      return false;
    } else {
      return true;
    }
  }

  /* Returns the coordinates of the points where tangents of the ellipse that go through the point
     statuibPositionvxBvxvxB touch the ellipse */
  void
  RdStationPolarizationRejector::GetErrorEllipseTangentPoints(TVector3& ellipseTangentPoint1,
                                                              TVector3& ellipseTangentPoint2,
                                                              const TVector3& stationPositionvxBvxvxB,
                                                              const TVector3& semiMajorAxisVxB,
                                                              const TVector3& semiMinorAxisVxB,
                                                              const double rotationAngle)
  {
    TVector3 rotatedStationPosition = stationPositionvxBvxvxB;
    rotatedStationPosition.RotateZ(-rotationAngle);
    const double semiMajorAxisLength2 = semiMajorAxisVxB.X() * semiMajorAxisVxB.X() + semiMajorAxisVxB.Y() *
                                        semiMajorAxisVxB.Y();
    const double semiMinorAxisLength2 = semiMinorAxisVxB.X() * semiMinorAxisVxB.X() + semiMinorAxisVxB.Y() *
                                        semiMinorAxisVxB.Y();
    const double divisor = rotatedStationPosition.Y() * rotatedStationPosition.Y() * semiMajorAxisLength2 +
                           rotatedStationPosition.X() * rotatedStationPosition.X() * semiMinorAxisLength2;
    const double squareRootTerm = pow(
      semiMinorAxisLength2 * semiMajorAxisLength2 * rotatedStationPosition.X() / divisor,
      2) +
                                  (pow(rotatedStationPosition.Y() * semiMajorAxisLength2,
                                       2) - pow(semiMajorAxisLength2, 2) * semiMinorAxisLength2) / divisor;
    const double x1 = semiMinorAxisLength2 * rotatedStationPosition.X() / divisor + sqrt(squareRootTerm);
    const double x2 = semiMinorAxisLength2 * rotatedStationPosition.X() / divisor - sqrt(squareRootTerm);
    const double y1 = semiMinorAxisLength2 / rotatedStationPosition.Y() - semiMinorAxisLength2 / semiMajorAxisLength2 *
                      x1 * rotatedStationPosition.X() / rotatedStationPosition.Y();
    const double y2 = semiMinorAxisLength2 / rotatedStationPosition.Y() - semiMinorAxisLength2 / semiMajorAxisLength2 *
                      x2 * rotatedStationPosition.X() / rotatedStationPosition.Y();
    ellipseTangentPoint1 = TVector3(x1, y1, 0);
    ellipseTangentPoint1.RotateZ(rotationAngle);
    ellipseTangentPoint2 = TVector3(x2, y2, 0);
    ellipseTangentPoint2.RotateZ(rotationAngle);
  }

  /* Returns the polar angles the most extreme polarization directions possible
     if the station is inside the SD core error ellipse */
  void
  RdStationPolarizationRejector::GetLocalPolMaxima(double& pol1,
                                                   double& pol2,
                                                   const double angleToMagneticField,
                                                   double maxChargeExcess)
  {
    const double phi = acos(maxChargeExcess / sin(angleToMagneticField));
    pol1 = atan2(maxChargeExcess * sin(phi), -maxChargeExcess * cos(phi) + sin(angleToMagneticField));
    pol2 = atan2(-maxChargeExcess * sin(phi), -maxChargeExcess * cos(phi) + sin(angleToMagneticField));
  }

  /* Returns the polar angle that the station at stationPosition vxBvxvxB would
     measure if the  SD core was at tangent point */
  double
  RdStationPolarizationRejector::GetTangentPol(const TVector3& stationPositionvxBvxvxB,
                                               const TVector3& tangentPoint,
                                               const double chargeExcess,
                                               const double angleToMagneticField)
  {
    const TVector3 diff = TVector3(stationPositionvxBvxvxB.X() - tangentPoint.X(),
                                   stationPositionvxBvxvxB.Y() - tangentPoint.Y(), 0);
    const double absDiff = sqrt(diff.X() * diff.X() + diff.Y() * diff.Y());
    TVector3 EFieldGeomagnetic = TVector3(-sin(angleToMagneticField), 0, 0);
    TVector3 EFieldChargeExcess = TVector3(-chargeExcess * diff.X() / absDiff, -chargeExcess * diff.Y() / absDiff, 0);
    TVector3 EField = EFieldGeomagnetic + EFieldChargeExcess;
    return atan2(EField.Y(), -EField.X());
  }

  TVector3
  RdStationPolarizationRejector::GetEfieldInShowerPlane(const utl::Vector& efield,
                                                        const utl::CoordinateSystemPtr& cs,
                                                        const utl::RadioGeometryUtilities& rdGeometryUtilities)
  {
    const utl::Point stationRecEfieldVector = utl::Point(efield.GetX(cs), efield.GetY(cs), efield.GetZ(cs), cs);
    double x;
    double y;
    double z;
    rdGeometryUtilities.GetVectorInShowerPlaneVxB(x, y, z, stationRecEfieldVector);
    return TVector3(x, y, z);
  }

  // Returns the uncertainty of the measured polarization angle based on the signal-to-noise-ratio
  double
  RdStationPolarizationRejector::GetPolAngleUncertainty(const double signalToNoise,
                                                        const TVector3 efield,
                                                        const double noiseVxB,
                                                        const double noiseVxVxB)
  {
    if (fNoiseErrorModel) {
      std::stringstream sstr;
      sstr << "EField is [" << efield[0] << ", " << efield[1] << "]. vxB noise RMS is " << noiseVxB
           << ", vxvxB noise RMS is " << noiseVxVxB << ".";
      INFODebug(sstr.str());
      double vxBError = (13.89 / pow(10., 6) + .594 * noiseVxB) * efield[1] / efield.Mag() / efield.Mag();
      double vxvxBError = (21.14 / pow(10., 6) + .487 * noiseVxVxB) * efield[0] / efield.Mag() / efield.Mag();
      return std::sqrt(vxBError * vxBError + vxvxBError * vxvxBError + (1. * utl::degree) * (1. * utl::degree));
    } else {
      const double sigma1 = (23.6 / sqrt(signalToNoise) - 4.3 / signalToNoise + 41.8 / pow(signalToNoise, 1.5)) * utl::degree;
      const double sigma2 = 1 * utl::deg;    // 1 degree alignment uncertainty
      return sqrt(sigma1 * sigma1 + sigma2 * sigma2);
    }
  }

  fwk::VModule::ResultFlag
  RdStationPolarizationRejector::Finish()
  {
    return eSuccess;
  }

}