CorsikaShowerFileGeometryProducer.cc

Go to the documentation of this file.

#include <io/CorsikaShowerFileGeometryProducer.h>

#include <io/RawCorsikaFile.h>
#include <io/CorsikaBlock.h>
#include <io/CorsikaIOException.h>
#include <io/CorsikaUtilities.h>

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/CoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/UTMPoint.h>
#include <utl/Point.h>
#include <utl/Math.h>
#include <utl/Vector.h>
#include <utl/GeometryUtilities.h>

#include <det/Detector.h>

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

#include <vector>
#include <sstream>

using namespace io;
using namespace utl;
using namespace std;


CorsikaShowerFileGeometryProducer::
CorsikaShowerFileGeometryProducer(const double corsikaZenith,
                                  const double corsikaAzimuth,
                                  const double obsLevelHeight,
                                  const bool multiCore,
                                  const utl::CoordinateSystemPtr refCS,
                                  const double corsikaRotation,
                                  const bool correctForRealisticMagneticFieldDeclination) :
  fCORSIKAZenith(corsikaZenith),
  fCORSIKAAzimuth(corsikaAzimuth),
  fObsLevelHeight(obsLevelHeight),
  fHasMultiCore(multiCore),
  fRefCS(refCS),
  fCoordinateRotation(corsikaRotation),
  fCorrectForRealisticMagneticFieldDeclination(correctForRealisticMagneticFieldDeclination)
{ }


/*
  This crazy function returns the coordinate system corresponding to
  the surface where CORSIKA saved particles coordinates. CORSIKA can
  save flat, as well as curved surfaces. The "shower core" in CORSIKA
  is always at 0,0 and this is the center of the CS.

  Warning: if you use an option in CORSIKA like CSCAT, TELESCOPE,
  AUGERHIT, etc. that makes assumptions on the location of specific
  "detectors", "volumes" or "positions", this must be 100% reflected
  here. The consquence can be that CS z-axis does not strictly point
  upwards, but is deflected to some extend, etc. Be careful. Don't
  make any implicit assumptions! Use points, vectors and CSs
  consitently and you will be fine.
*/

utl::CoordinateSystemPtr
CorsikaShowerFileGeometryProducer::MakeGroundParticleCoordinateSystem(const utl::Point& anyPointOnShowerAxis)
{
  /* There are several options that have to be considered.

   1) Normal case: CORSIKA observation level height must correspond
      to Offline altitude. CS origin is shifted to height altitude.

   2) Case where it is important to retain the alignment of CORSIKA
      output with other external positions (TELESCOPS, AUGERHIT,
      Radio-Antennas, COREAS...), which are given in a particular _fixed_
      reference frame within CORSIKA. This requires to put the CS
      into this particular CORSIKA reference frame. For example if
      site-CS is used for "TELESCOPES" then CS must be on altitute
      where the z-coordinate in "site-CS" corresponds to the CORSIKA
      observation level. (this is tested for flat CORSIKA. need to confirm
      for curved and curved-output, too). External(core)-CS will
      then have a z-axis parallel to site-CS, not strictly pointing
      upward. Zenith, azimuth are relative to this particular CS,
      too!

   PROBLEM TO SOLVE: Which reference frame should we select for
                     Auger in CORSIKA as a convention ??? Need to decide !!!
  */

  static bool showThisMessageOnlyOnce = true;
  const bool useAugerCS = !fHasMultiCore;

  const utl::CoordinateSystemPtr cs = fwk::LocalCoordinateSystem::Create(anyPointOnShowerAxis);
  utl::CoordinateSystemPtr coreCS;

  double north = 0;
  if (fCorrectForRealisticMagneticFieldDeclination) {
    // magnetic field declination
    north = det::Detector::GetInstance().GetGeoMagneticField().GetDeclination(anyPointOnShowerAxis);
    ostringstream info;
    info << "Magnetic field declination: " << north/utl::degree << " deg";
    INFO(info);
  }

  if (useAugerCS) {

    // Option 1
    // second: move CORSIKA CS along shower axis so that Auger altitude
    // is identical to CORSIKA observation level height create vector

    // pointing along the momvement of the shower in CORSIKA frame (see
    // Fig 1 of CORSIKA_GUIDE7.5600.pdf). This is different from normal
    // Auger convention!

    const double sinZenith = sin(fCORSIKAZenith);
    const double cosZenith = cos(fCORSIKAZenith);

    // This is for Auger azimuth, corsika phi always relative to magnetic north
    // CORSIKA direction is direction of motion of CR, not direction from where CR is coming
    const double augerAzimuth = Corsika::CorsikaAzimuthToAuger(fCORSIKAAzimuth + kPi, north);
    const Vector direction(-cos(augerAzimuth)*sinZenith,
                           -sin(augerAzimuth)*sinZenith,
                           -cosZenith, cs);

    const vector<Point> intersections =
      utl::Intersection(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84),
                        fObsLevelHeight,
                        Line(anyPointOnShowerAxis, direction));

    {
      const UTMPoint theUTMPoint(anyPointOnShowerAxis, ReferenceEllipsoid::eWGS84);
      const double altitude = theUTMPoint.GetHeight();
      ostringstream inf;
      inf << "Observation level height: " << fObsLevelHeight/m << " m, Offline core height: " << altitude/m << " m";
      INFO(inf);
    }

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

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

    for (const auto& inter : intersections) {
      const Vector coreToPoint = inter - anyPointOnShowerAxis;
      const double dist = coreToPoint.GetMag();
      if (first) {
        closest = dist;
        toPlane = coreToPoint;
        first = false;
      } else if (dist < closest) {
        closest = dist;
        toPlane = coreToPoint;
      }
    }

    if (showThisMessageOnlyOnce) {
      ostringstream inf;
      inf << "CORSIKA particle ground plane is displaced by " << (toPlane*direction)/m << "m "
             "along shower axis with respect to Offline ground level.";
      INFO(inf);
      showThisMessageOnlyOnce = false;
    }

    coreCS = fwk::LocalCoordinateSystem::Create(Point(0, 0, 0, CoordinateSystem::Translation(toPlane, cs)));

  // *******************************************************************
  } else { //  CS is not an AugerCoordinateSystem

    // Option 2

    // corsika z-axis might not be pointing straight up !

    const Vector shift = anyPointOnShowerAxis - Point(0, 0, 0, fRefCS);

    // now move refCS origin to core
    coreCS = CoordinateSystem::Translation(shift, fRefCS);

    if (showThisMessageOnlyOnce) {
      ostringstream inf;
      inf << "CORSIKA particle ground plane is rotated by "
          << Angle(Vector(0, 0, 1, coreCS), Vector(0, 0, 1, cs))/deg << " deg "
             "in zenith with respect to Offline ground level.";
      INFO(inf);
      showThisMessageOnlyOnce = false;
    }
  }

  /*
     Final step:

     rotate Auger CS into CORSIKA CS

     in CORSIKA x axis points to magnetic north. In Offline x-axis points to geographic east.
  */
  const double angle = Corsika::CorsikaAzimuthToAuger(fCoordinateRotation, north);
  ostringstream inf;
  inf << "Rotation (around z axis) between Auger CS and Corsika CS: " << angle / deg << " deg.";
  INFO(inf);
  const CoordinateSystemPtr particlePlaneCS = CoordinateSystem::RotationZ(angle, coreCS);
  return particlePlaneCS;
}