ErrorPropagation.cc

Go to the documentation of this file.

#include <utl/config.h>
#include "ErrorPropagation.h"
#include "ConversionUtil.h"

#include <utl/NumericalErrorPropagation.h>
#include <utl/AugerException.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/ModifiedJulianDate.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/Point.h>
#include <utl/UTMPoint.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/Transformation.h>
#include <utl/Vector.h>

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

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

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>

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

#include <iostream>
#include <string>

using namespace std;
using namespace otoa;
using namespace otoa::err;
using namespace utl;
using namespace fwk;


std::vector<double>
CalculateEquatorialCoordinates::operator()(const std::vector<double>& p)
  const
{
  if (p.size() != 7) {
    ostringstream msg;
    msg << "CalculateEquatorialCoordinates needs 7 parameters, but was called with " << p.size();
    throw utl::OutOfBoundException(msg.str());
  }

  const double julianDay = p[0];
  const double gmst      = p[1];
  const double easting   = p[2];  // = x
  const double northing  = p[3];  // = y
  const double altitude  = p[4];  // = z
  double theta           = p[5];
  double phi             = p[6];

  const char band = 'H';
  const int zone = 19;
  const UTMPoint pt(northing, easting, altitude, zone, band, ReferenceEllipsoid::eWGS84);

  double latitude;
  double longitude;
  double height;
  boost::tuple<double&, double&, double&>(latitude, longitude, height) =
    pt.GetGeodeticCoordinates();

  phi += 0.5*kPi;  // phi should be -phi
                   // (Auger pointing east, and going to north, west, etc)
  phi *= -1;       // change sign of phi
  theta = 0.5*kPi - theta;  // theta

  const double sinTheta = sin(theta);
  const double cosTheta = cos(theta);
  const double cosPhi = cos(phi);
  const double sinLatitude = sin(latitude);
  const double cosLatitude = cos(latitude);
  const double dec = asin(sinLatitude*sinTheta - cosLatitude*cosTheta*cosPhi);
  const double secDec = 1 / cos(dec);
  const double haSin = cosTheta * sin(phi) * secDec;
  const double haCos = (cosLatitude*sinTheta + sinLatitude*cosTheta*cosPhi) * secDec;

  const double ha = (haSin < 0) ? 2*kPi - acos(haCos) : acos(haCos);

  // TODO lives in ConversionUtils.cc, should be in utl along with THIS code (or libnova)
  const double nutHour = otoa::CalculateNutationCorrection(julianDay);
  // local sideral time
  double lmst = fmod(gmst + nutHour, 24.)*15 + longitude/deg;

  //get the lmst between 0 and 360 degrees
  if (lmst > 0)
    lmst -= int(lmst/360)*360;
  else
    lmst += 360 - int(lmst/360)*360;

  double ra = lmst*deg - ha;

  if (ra < 0)
    ra += 2*kPi;
  if (ra > 2*kPi)
    ra -= 2*kPi;

  vector<double> result(2);
  if (-2*kPi <= ra && ra <= 2*kPi) {
    result[0] = ra;
    result[1] = dec;
  }

  return result;
}


std::vector<double>
CalculateGalacticCoordinates::operator()(const std::vector<double>& p)
  const
{
  if (p.size() != 7) {
    ostringstream msg;
    msg << "CalculateGalacticCoordinates needs 7 parameters, but was called with " << p.size();
    throw utl::OutOfBoundException(msg.str());
  }

  CalculateEquatorialCoordinates calcEqu;
  const vector<double> q = calcEqu(p);

  const double ra  = q[0];
  const double dec = q[1];

  const double raGalEq   = 282.86*deg;  // ra of the ascending node
  const double galLg2000 = 32.933*deg;  //
  const double eqGal     = 27.13*deg;   // pi/2-angle betweem eq and gal equator

  const double cosDec = cos(dec);
  const double sinDec = sin(dec);
  const double cosEqGal = cos(eqGal);
  const double galLatitude = asin(sinDec*sin(eqGal) - cosDec*sin(ra - raGalEq)*cosEqGal);
  const double secGalLatitude = 1 / cos(galLatitude);
  const double cosLon = cosDec*cos(ra - raGalEq) * secGalLatitude;
  const double sinLon = (cosDec*sin(ra - raGalEq)*sin(eqGal) + sinDec*cosEqGal) * secGalLatitude;

  double galLongitude = ((cosLon < 0 && sinLon < 0) || (cosLon > 0 && sinLon < 0)) ?
    2*kPi - acos(cosLon) + galLg2000 : acos(cosLon) + galLg2000;

  //while(GLongitude>2*kPi) GLongitude -= 2.0*kPi;
  galLongitude = fmod(galLongitude, 2*kPi);
  if (galLongitude < 0)
    galLongitude += 2*kPi;

  vector<double> result(2);
  result[0] = galLongitude;
  result[1] = galLatitude;
  return result;
}


CalculateFdCorePosition::CalculateFdCorePosition(const fevt::Eye& eye) :
  fEye(eye)
{ }


std::vector<double>
CalculateFdCorePosition::operator()(const std::vector<double>& p)
  const
{
  const double rp       = p[0];
  const double chi0     = p[1];
  const double sDPtheta = p[2];
  const double sDPphi   = p[3];

  const ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
  const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(fEye);
  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const Point eyeposition = detEye.GetPosition();

  const Vector sdp(1, sDPtheta, sDPphi, eyeCS, Vector::kSpherical);

  const Vector vertical(0, 0, 1, eyeCS);
  const Vector horizontalInSDP = Normalized(Cross(sdp, vertical)); // horizontal

  const Vector coreEyeVec = rp / sin(kPi - chi0) * horizontalInSDP;
  const Point core = eyeposition + coreEyeVec;

  // core (first of all do the utm-exception check !!!!)
  double northing = 0;
  double easting = 0;
  // altitude is fixed to local ground altitude
  try {
    const utl::UTMPoint pt(core, wgs84);
    northing = pt.GetNorthing();
    easting = pt.GetEasting();
  } catch (UTMPoint::UTMZoneException&) {
    ERROR("UTMZoneException: fd rec shower invalid");
  } catch (UTMPoint::UTMException&) {
    ERROR("UTMException: fd rec shower invalid");
  }

  std::vector<double> result(2);
  result[0] = easting;
  result[1] = northing;

  return result;
}


CalculateFdArrivalDirection::CalculateFdArrivalDirection(const fevt::Eye& eye) :
  fEye(eye)
{ }


std::vector<double>
CalculateFdArrivalDirection::operator()(const std::vector<double>& p)
  const
{
  const double chi0     = p[0];
  const double sDPtheta = p[1];
  const double sDPphi   = p[2];

  const evt::ShowerFRecData& frec = fEye.GetRecData().GetFRecShower();
  const utl::Point& core = frec.GetCorePosition();
  const fdet::Eye& detEye = det::Detector::GetInstance().GetFDetector().GetEye(fEye);
  const CoordinateSystemPtr eyeCS = detEye.GetEyeCoordinateSystem();
  const CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(core);

  const Vector sdp(1, sDPtheta, sDPphi, eyeCS, Vector::kSpherical);

  const Vector vertical(0, 0, 1, eyeCS);
  const Vector horizontalInSDP = Normalized(Cross(sdp, vertical)); // horizontal

  const Transformation rot = Transformation::Rotation(-chi0, sdp, eyeCS);

  const Vector axis = Normalized(rot * horizontalInSDP);           // apply rotation

  vector<double> result(2);
  result[0] = axis.GetTheta(localCS);
  result[1] = axis.GetPhi(localCS);

  return result;
}