GeometryExample.cc

Go to the documentation of this file.

#include "GeometryExample.h"

#include <utl/ErrorLogger.h>

#include <evt/Event.h>

//#include <fwk/CentralConfig.h>
#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <sdet/Station.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <utl/MathConstants.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/UTMPoint.h>
#include <utl/AxialVector.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/PhysicalConstants.h>
#include <utl/AugerUnits.h>
#include <utl/ReferenceEllipsoid.h>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_io.hpp>

using namespace det;
using namespace sdet;
using namespace fdet;
using namespace std;
//using namespace sevt;
using namespace evt;
using namespace fwk;
using namespace utl;

using namespace GeometryExampleNS;


VModule::ResultFlag
GeometryExample::Init()
{
  return eSuccess;
}


/*
  The coding example here doesn't include any of addition and subtraction
  operations and does not include most member functions (Get* /Set*). Those can
  be found in the class documentation for the geometry classes. The things
  included here are mostly utm coordinates and coordinate system pointers and an
  example using dot and cross product.
*/

VModule::ResultFlag
GeometryExample::Run(evt::Event& /*event*/)
{
  const TimeStamp times = UTCDateTime(2005,1,1,0,0).GetTimeStamp();
  Detector::GetInstance().Update(times);
  const SDetector& sDetector = Detector::GetInstance().GetSDetector();
  const FDetector& fDetector = Detector::GetInstance().GetFDetector();
  INFO("GeometryExample::Run()");

  cout << "\nCreating a Coordinate System from scratch\n"
          "-----------------------------------------\n";

  const char band = 'H';
  const string ellipsoidName("WGS84");
  const int zone = 19;
  const double northing = 6081459*m;
  const double easting = 460373*m;
  const double altitude = 1391*m;
  const ReferenceEllipsoid::EllipsoidID ellipsoid = ReferenceEllipsoid::GetEllipsoidIDFromString(ellipsoidName);

  const UTMPoint utmposition(northing, easting, altitude, zone, band, ellipsoid);
  const CoordinateSystemPtr cs1 = LocalCoordinateSystem::Create(utmposition);
  cout << "the origin of Coordinate System 1 will be: " << utmposition.GetCoordinates() << "\n";

  UTMPoint utmposition2(6081479*m, 460363*m, 1400*m, 19, 'H', ReferenceEllipsoid::GetEllipsoidIDFromString("WGS84"));
  cout << "the origin of Coordinate System 2 will be: " << utmposition2.GetCoordinates() << "\n";
  const CoordinateSystemPtr cs2 = LocalCoordinateSystem::Create(utmposition2);

  // create a Point from a UTMPoint

  const Point cs1_origin = utmposition.GetPoint();

  cout << "CS 1 origin relative to CS 2: " << cs1_origin.GetCoordinates(cs2) << "\n";

  // Create a right-handed coordinate system from a point (cs1_origin) and a direction (zenith, azimuth)
  // z-axis oposite a hypothetical shower direction,
  // y-axis in the horizontal plane,
  // x-axis pointing in the upstream direction
  const double zenith = 30*deg;
  const double azimuth = 270*deg;
  const CoordinateSystemPtr cs3(
    CoordinateSystem::RotationY(
      zenith,
      CoordinateSystem::RotationZ(
        azimuth,
        LocalCoordinateSystem::Create(cs1_origin)
      )
    )
  );

  cout << "CS 1 origin relative to CS 3: " << cs1_origin.GetCoordinates(cs3) << "\n";

  // Calculate the latitude, longitude and height on the WGS84 ellipsoid for a
  // point defined in PampaAmarilla coordinates system (defined at the approximate
  // center of the Southern site)

  const CoordinateSystemPtr pampaCS =
    CoordinateSystemRegistry::Get("PampaAmarilla");
  const ReferenceEllipsoid wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  const Point pointInPampa = Point(0., 0., 2.*km, pampaCS);
  const boost::tuple<double, double, double> latLonHeight = wgs84.PointToLatitudeLongitudeHeight(pointInPampa);
  cout << "latitude, longitude, height " << latLonHeight << endl;

  // create a vector and a point in one coordinate system

  const Point point1(0,10,20, cs1);
  const Point point3(-10,10,10, cs1);
  const Point point2(10,kPi/3,0, cs1, Point::kCylindrical);
  const Vector vector1(1,0,0, cs1);
  const Vector vector2(0,1,0, cs1);

  // read in the other using more than one basis.

  cout << "\nPoints and vectors in different coordinate systems\n"
          "--------------------------------------------------\n"
          "coordinate System 1 \t\tCoordinate System 2\n"
       << point1.GetCoordinates(cs1) << " \t\t" << point1.GetCoordinates(cs2) << "\n"
       << point2.GetCoordinates(cs1) << " \t" << point3.GetCoordinates(cs2) << "\n"
       << point3.GetCoordinates(cs1) << " \t\t\t" << point2.GetCoordinates(cs2) << "\n"
       << vector1.GetCoordinates(cs1) << " \t" << vector1.GetCoordinates(cs2) << "\n"
       << vector2.GetCoordinates(cs1) << " \t" << vector2.GetCoordinates(cs2) << "\n" ;

  cout << "NOTE that the vector's components don't depend on the coordinate system!\n";

  // The detector interface provides the coordinate systems relative to most things we need and a few default ones:

  cout << "\nUsing the Coordinate systems provided by the framework\n"
          "------------------------------------------------------\n";

  const CoordinateSystemPtr pampaCs = fwk::CoordinateSystemRegistry::Get("PampaAmarilla");

  const CoordinateSystemPtr llancaneloCs = sDetector.GetStation(221).GetLocalCoordinateSystem();

  const CoordinateSystemPtr coihuecoCs = fDetector.GetEye("Coihueco").GetEyeCoordinateSystem();
  const CoordinateSystemPtr losLeonesCs = fDetector.GetEye("Los Leones").GetEyeCoordinateSystem();
  const CoordinateSystemPtr moradosCs = fDetector.GetEye("Los Morados").GetEyeCoordinateSystem();

  const CoordinateSystemPtr siteCs = Detector::GetInstance().GetSiteCoordinateSystem();

  const Point corePosition(0,0,0, pampaCs);
  cout << "setting a 'core' position at:\n"
          " " << corePosition.GetCoordinates(pampaCs) << ", in 'Pampa Amarilla' coordinate system, \n"
          " " << corePosition.GetCoordinates(coihuecoCs) << " relative to Coihueco\n"
          " " << corePosition.GetCoordinates(siteCs) << " in the site's cs\n";

  const Vector direction(1,2*kPi/3,kPi/6, pampaCs, Vector::kSpherical);
  cout << "with direction \n  "
       << direction.GetCoordinates(pampaCs) << " in 'Pampa Amarilla' cartesian system, \n  "
       << "(" << direction.GetR(pampaCs) << ", " << direction.GetTheta(pampaCs)/deg << ", " << direction.GetPhi(pampaCs)/deg
       << ") in 'Pampa Amarilla' spherical coordinates \n  "
       << direction.GetCoordinates(coihuecoCs) << "cartesian coordinates relative to Coihueco \n";

  cout << "Note (again) how the direction's component didn't change much, as expected. "
          "The direction is a Vector, whereas the position is Point. The slight changes are due to the earth curvature.\n";

  // calculate the distance from the core to Los Leones and to llancanelo:

  const double ll_core = (fDetector.GetEye("Los Leones").GetPosition() - corePosition).GetMag();
  const double llancanelo_core = (sDetector.GetStation(221).GetPosition() - corePosition).GetMag();

  cout << "\nDistance from Los Leones to the core: " << ll_core/km << " km\n"
          "Distance from LLancanelo (Local Station #221) to the core: " << llancanelo_core/km << " km\n";

  // use cross product to calculate the SDP normal vector given the direction and core position (for Los Leones)

  AxialVector sdp = Cross(direction, corePosition - fDetector.GetEye("Los Leones").GetPosition());
  sdp.Normalize();

  cout << "SEP plane: the equation for the Shower Eye Plane is "
          "(" << sdp.GetX(losLeonesCs) << ") x + (" << sdp.GetY(losLeonesCs) << ") y + (" << sdp.GetZ(losLeonesCs) << ") z = 0\n";

  // use dot product to calculate the angle between vertical vectors at Coihueco and Los Morados.

  const Vector coihuecoVertical(0,0,1, coihuecoCs);
  const Vector moradosVertical(0,0,1, moradosCs);

  cout << "the angle between the vertical line at Coihueco and Los Morados: "
       << Angle(coihuecoVertical, moradosVertical)/deg << " degrees. The earth is curved!\n"
       << "same number, measuring the zenithal angle of the Morados vertical in Coihueco's CS: "
       << moradosVertical.GetTheta(coihuecoCs)/deg << " degrees\n";

  return eSuccess;
}


VModule::ResultFlag
GeometryExample::Finish()
{
  return eSuccess;
}