RadioGeometryUtilities.cc

Go to the documentation of this file.

#include "RadioGeometryUtilities.h"

#include <utl/Math.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/Line.h>
#include <utl/Plane.h>
#include <utl/GeometryUtilities.h>
#include <utl/AxialVector.h>

#include <utl/CoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/TransformationMatrix.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/AugerUnits.h>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <TVector3.h>

#include <vector>
#include <sstream>

using namespace utl;


RadioGeometryUtilities::RadioGeometryUtilities(const Vector& showeraxis,
                                               const CoordinateSystemPtr cs,
                                               const Vector& magneticField) :
  fcs(cs),
  fMagneticField(magneticField),
  fTransformationMatrix(3, 3),
  fInverseTransformationMatrix(3, 3)
{
  const Vector showerdirection = -showeraxis;
  const Vector vxB = Normalized(Cross(showerdirection, magneticField));
  const Vector vxvxB = Normalized(Cross(showerdirection, vxB));
  const Vector e3 = Normalized(Cross(vxB, vxvxB));

  fTransformationMatrix(0, 0) = vxB.GetX(cs);
  fTransformationMatrix(0, 1) = vxB.GetY(cs);
  fTransformationMatrix(0, 2) = vxB.GetZ(cs);
  fTransformationMatrix(1, 0) = vxvxB.GetX(cs);
  fTransformationMatrix(1, 1) = vxvxB.GetY(cs);
  fTransformationMatrix(1, 2) = vxvxB.GetZ(cs);
  fTransformationMatrix(2, 0) = e3.GetX(cs);
  fTransformationMatrix(2, 1) = e3.GetY(cs);
  fTransformationMatrix(2, 2) = e3.GetZ(cs);

  fInverseTransformationMatrix = boost::numeric::ublas::trans(fTransformationMatrix);
}


double
RadioGeometryUtilities::GetEarlyLateCorrectionFactor(const utl::Point& showerCore,
                                                     const utl::Point& stationPosition,
                                                     const utl::Point& showerMax,
                                                     const utl::Vector& showerAxis)
{
  const utl::Vector& showerDirection = -1 * showerAxis;
  const double dist = (showerMax - showerCore).GetMag();
  return 1 + (stationPosition - showerCore) * showerDirection / dist;
}


Vector
RadioGeometryUtilities::GetChargeExcessVector(const Point& core,
                                              const Vector& showeraxis,
                                              const Point& stationPosition)
{
  const Vector showeraxisUnit = Normalized(showeraxis);
  const Vector chargeExcessVector =
    core - stationPosition + showeraxisUnit * ((stationPosition - core) * showeraxisUnit);  // showeraxis has to be normalized
  return Normalized(chargeExcessVector);
}


Vector
RadioGeometryUtilities::GetLorentzVector(const Vector& showeraxis, const Vector& vMagField)
{
  return Cross(Normalized(showeraxis), Normalized(vMagField));
}


double
RadioGeometryUtilities::GetAngleToLorentzVector(const Vector& measuredEField,
                                                const Vector& showeraxis,
                                                const Vector& vMagField)
{
  return Angle(measuredEField, RadioGeometryUtilities::GetLorentzVector(showeraxis, vMagField));
}


double
RadioGeometryUtilities::GetDistanceToAxis(const Vector& ShowerAxis,
                                          const Point& CorePosition,
                                          const Point& AntennaPosition)
{
  const Line line = Line(CorePosition, ShowerAxis);
  return Distance(line, AntennaPosition);
}


double
RadioGeometryUtilities::GetSignalCorrectionFactor(const Point& core,
                                                  const Vector& showeraxis,
                                                  const Point& stationPosition,
                                                  const Vector& vMagField,
                                                  const double chargeExcessStrength)
{
  return 1 / (RadioGeometryUtilities::GetChargeExcessVector(core, showeraxis, stationPosition)
              * chargeExcessStrength
              + RadioGeometryUtilities::GetLorentzVector(showeraxis, vMagField)).GetMag();
}


Vector
RadioGeometryUtilities::GetExpectedEFieldVector(const Point& core,
                                                const Vector& showeraxis,
                                                const Point& stationPosition,
                                                const Vector& vMagField,
                                                const double chargeExcessStrength)
{
  return RadioGeometryUtilities::GetLorentzVector(showeraxis, vMagField) +
         chargeExcessStrength * RadioGeometryUtilities::GetChargeExcessVector(core, showeraxis, stationPosition);
}


double
RadioGeometryUtilities::GetAngleToEFieldExpectation(const Vector& measuredEField,
                                                    const Point& core,
                                                    const Vector& showeraxis,
                                                    const Point& stationPosition,
                                                    const Vector& vMagField,
                                                    const double chargeExcessStrength)
{
  return
    Angle(
      measuredEField,
      RadioGeometryUtilities::GetExpectedEFieldVector(core, showeraxis, stationPosition, vMagField, chargeExcessStrength)
    );
}


const Point
RadioGeometryUtilities::GetCoreAtObservationLevel(const Point& core,
                                                  const Vector& showerAxis,
                                                  const CoordinateSystemPtr cs,
                                                  const double altitudeObsLvl)
{
  const std::vector<Point> intersections =
    Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84), altitudeObsLvl, Line(core, showerAxis));

  if (intersections.empty())
    WARNING("Shower does not intersect observation level.");

  // pick closest point (there are typically two solutions)
  Vector toPlane;
  double closest2 = 0;
  bool first = true;

  for (const auto& inter : intersections) {
    const Vector coreToPoint = inter - core;
    const double dist2 = coreToPoint.GetMag2();
    if (first || dist2 < closest2) {
      closest2 = dist2;
      toPlane = coreToPoint;
      first = false;
    }
  }

  const Point newCore = Point(0, 0, 0, CoordinateSystem::Translation(toPlane, cs));

  std::ostringstream inf;
  inf << "Get new core displaced by " << (toPlane * showerAxis) / m << "m "
         "along shower axis for an observation level of " << altitudeObsLvl / m << "m";
  INFO(inf);

  return newCore;
}


void
RadioGeometryUtilities::GetVectorInShowerPlaneVxB(double& x, double& y, double& z, const Point& point)
  const
{
  boost::numeric::ublas::vector<double> v(3);
  v(0) = point.GetX(fcs);
  v(1) = point.GetY(fcs);
  v(2) = point.GetZ(fcs);
  const boost::numeric::ublas::vector<double> w = boost::numeric::ublas::prod(fTransformationMatrix, v);
  x = w(0);
  y = w(1);
  z = w(2);
}


void
RadioGeometryUtilities::GetVectorInShowerPlaneVxB(double& x, double& y, double& z, const TVector3& vector)
  const
{
  boost::numeric::ublas::vector<double> v(3);
  v(0) = vector[0];
  v(1) = vector[1];
  v(2) = vector[2];
  const boost::numeric::ublas::vector<double> w = boost::numeric::ublas::prod(fTransformationMatrix, v);
  x = w(0);
  y = w(1);
  z = w(2);
}


TraceV3D
RadioGeometryUtilities::GetTraceInShowerPlaneVxB(const TraceV3D& trace)
  const
{
  TraceV3D rotatedTrace;
  rotatedTrace.SetBinning(trace.GetBinning());
  //double max1 = 0;
  //double max2 = 0;
  for (unsigned int i = 0; i < trace.GetSize(); ++i) {
    boost::numeric::ublas::vector<double> v(3);
    v(0) = trace[i][0];
    v(1) = trace[i][1];
    v(2) = trace[i][2];
    const boost::numeric::ublas::vector<double> w = boost::numeric::ublas::prod(fTransformationMatrix, v);
    /*const double max_tmp1 = sqrt(v(0)*v(0) + v(1)*v(1) + v(2)*v(2));
    const double max_tmp2 = sqrt(w(0)*w(0) + w(1)*w(1) + w(2)*w(2));
    if (max_tmp1 > max1)
      max1 = max_tmp1;
    if (max_tmp2 > max2)
      max2 = max_tmp2;*/
    Vector3D ww;
    ww = w(0), w(1), w(2);
    rotatedTrace.PushBack(ww);
  }
  return rotatedTrace;
}


TraceV3D
RadioGeometryUtilities::GetTraceFromShowerPlaneVxB(const TraceV3D& trace)
  const
{
  TraceV3D rotatedTrace;
  rotatedTrace.SetBinning(trace.GetBinning());
  for (unsigned int i = 0; i < trace.GetSize(); ++i) {
    boost::numeric::ublas::vector<double> v(3);
    v(0) = trace[i][0];
    v(1) = trace[i][1];
    v(2) = trace[i][2];
    const boost::numeric::ublas::vector<double> w = boost::numeric::ublas::prod(fInverseTransformationMatrix, v);
    Vector3D ww;
    ww = w(0), w(1), w(2);
    rotatedTrace.PushBack(ww);
  }
  return rotatedTrace;
}


Point
RadioGeometryUtilities::GetVectorFromShowerPlaneVxB(const double x, const double y, const double z,
                                                    const bool verticalZero)
  const
{
  boost::numeric::ublas::vector<double> v(3);
  v(0) = x;
  v(1) = y;
  v(2) = z;
  boost::numeric::ublas::vector<double> w = boost::numeric::ublas::prod(fInverseTransformationMatrix, v);
  if (verticalZero)
    w(2) = 0;

  return Point(w(0), w(1), w(2), fcs);
}


Point
RadioGeometryUtilities::GetVectorFromShowerPlaneVxB(const double x, const double y)
  const
{
  boost::numeric::ublas::vector<double> v(3);
  v(0) = x;
  v(1) = y;
  v(2) = GetHeightInShowerPlane(x, y);
  const boost::numeric::ublas::vector<double> w = boost::numeric::ublas::prod(fInverseTransformationMatrix, v);

  return Point(w(0), w(1), w(2), fcs);
}


double
RadioGeometryUtilities::GetHeightInShowerPlane(const double x, const double y)
  const
{
  return -1. * (fTransformationMatrix(0, 2) * x + fTransformationMatrix(1, 2) * y) / fTransformationMatrix(2, 2);
}


double
RadioGeometryUtilities::GetAngleToEFieldExpectation2D(const Vector& measuredEField,
                                                      const Point& core, const Vector& showeraxis,
                                                      const Point& stationPosition, const Vector& vMagField,
                                                      const double chargeExcessStrength, const CoordinateSystemPtr localCS)
{
  const Vector vxB = Normalized(Cross(showeraxis, -vMagField));
  const Vector vxvxB = Normalized(Cross(showeraxis, vxB));
  const Vector e3 = Normalized(Cross(vxB, vxvxB));

  double transfomationMatix[3][3];
  transfomationMatix[0][0] = vxB.GetX(localCS);
  transfomationMatix[1][0] = vxB.GetY(localCS);
  transfomationMatix[2][0] = vxB.GetZ(localCS);
  transfomationMatix[0][1] = vxvxB.GetX(localCS);
  transfomationMatix[1][1] = vxvxB.GetY(localCS);
  transfomationMatix[2][1] = vxvxB.GetZ(localCS);
  transfomationMatix[0][2] = e3.GetX(localCS);
  transfomationMatix[1][2] = e3.GetY(localCS);
  transfomationMatix[2][2] = e3.GetZ(localCS);

  //TransformationMatrix transM =
  //    TransformationMatrix::TransformToBasis(vxB.GetX(localCS), vxvxB.GetX(localCS), e3.GetX(localCS),
  //                                                vxB.GetY(localCS), vxvxB.GetY(localCS), e3.GetY(localCS),
  //                                                vxB.GetZ(localCS), vxvxB.GetZ(localCS), e3.GetZ(localCS));
  const Vector exp_efield =
    RadioGeometryUtilities::GetExpectedEFieldVector(core, showeraxis, stationPosition, vMagField, chargeExcessStrength);
  const double v_exp_efield[3] =
    { exp_efield.GetX(localCS), exp_efield.GetY(localCS), exp_efield.GetZ(localCS) };
  double v[3] = { 0 };
  for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
      v[i] += transfomationMatix[j][i] * v_exp_efield[j];
  const TVector3 exp_efield_transformed(v);

  const double v_measured_efield[3] =
    { measuredEField.GetX(localCS), measuredEField.GetY(localCS), measuredEField.GetZ(localCS) };
  for (int i = 0; i < 3; ++i) {
    v[i] = 0;
    for (int j = 0; j < 3; ++j) {
      v[i] += transfomationMatix[j][i] * v_measured_efield[j];
    }
  }
  TVector3 measured_efield_transformed(v);
  if (measured_efield_transformed.x() > 0)
    measured_efield_transformed *= -1;

  const double diff = NormalizeAngleMinusPiPi(exp_efield_transformed.Phi() - measured_efield_transformed.Phi());
  return diff;
}