FieldOfViewCalculator.cc

Go to the documentation of this file.

#include <iomanip>
#include <algorithm>
#include <vector>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/UTMPoint.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TransformationMatrix.h>
#include <utl/TimeStamp.h>
#include <utl/AugerCoordinateSystem.h>

#include <utl/Particle.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>
#include <evt/LaserData.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/Telescope.h>

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

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

#include <atm/AttenuationResult.h>
#include <atm/ScatteringResult.h>
#include <atm/Atmosphere.h>
#include <atm/ProfileResult.h>

#include "FieldOfViewCalculator.h"

#include <TH2D.h> // for debug
#include <TH1D.h> // for debug
#include <TCanvas.h>// for debug
#include <TGraph.h>// for debug
#include <TFile.h>// for debug

using namespace std;
using namespace FieldOfViewCalculatorKG;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fdet;
using namespace atm;


VModule::ResultFlag
FieldOfViewCalculator::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("FieldOfViewCalculator");

  topB.GetChild("binning").GetData(fBinning);
  topB.GetChild("maxAngleBinWidth").GetData(fMaxAngleBinWidth);

  // -------- do some output ------------
  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "         time binning: " << fBinning << " (relative to FADC sampling)" << "\n"
          "    min. angular step: "<< fMaxAngleBinWidth/deg << " deg \n ";
  INFO(info);

  return eSuccess;
}


/*
  find out the tMin and tMax for all telescopes (time from showerInitialPos to telescope)
  (-> the position along the shower track, where it enters the FOV of a tel)
  AND
  initialize the photon - distance - traces

  INFO:

  The geometry is defined relative to the "core" location:
  p = core + distance * axis"

  and not any more (as previously, before Di 7. Mär 10:14:57 CET 2017) as
  p = showerInitialPos + distance * axis
*/


VModule::ResultFlag
FieldOfViewCalculator::Run(evt::Event& event)
{
  if (!event.HasFEvent()) {
    ERROR("Missing FEvent. Check your module configurations and ModuleSequence.");
    return eFailure;
  }
  fevt::FEvent& FDevent = event.GetFEvent();

  if (!event.HasSimShower()) {
    ERROR("Event has no simulated shower.");
    return eFailure;
  }
  evt::ShowerSimData& simShower = event.GetSimShower();

  // Shower initial and final positions
  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  const Point core(0, 0, 0, showerCS);
  const Vector showerAxis(0, 0, -1, showerCS);

  //const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::GetWGS84();
  //const double coreZ = wgs84.PointToLatitudeLongitudeHeight(core).get<2>();

  Detector& detector = Detector::GetInstance();
  const Atmosphere& atmo = detector.GetAtmosphere();

  // --------------------------------------------------------------------------
  // find the part of the shower, that was emitting light

  detector.GetAtmosphere().InitSlantProfileModel(core, -showerAxis, 10.*g/cm2);

  const ProfileResult& distanceVsSlantDepth = atmo.EvaluateDistanceVsSlantDepth();
  double atmDepthMin = distanceVsSlantDepth.MinX();
  double atmDepthMax = distanceVsSlantDepth.MaxX();
  double distStart = distanceVsSlantDepth.Y(atmDepthMin);
  double distStop = distanceVsSlantDepth.Y(atmDepthMax);

  // before injection there is no dE/dX profile
  double profStartX = simShower.GetXInject(); 
  double profStartDist = distanceVsSlantDepth.Y(profStartX);
  if (profStartDist > distStart)
    distStart = profStartDist;

#if 1
  {
    // debug
    const double distInitFinal = distStop - distStart;
    ostringstream info;
    info << "Shower geometry parameters: distance_start=" << distStart/km << " km, "
         << " distance_stop=" << distStop/km << " km, "
         << " length=" << distInitFinal/km << " km";
    INFO(info);
  }
#endif

  TimeStamp timeAtCore = simShower.GetTimeStamp();
  const FDetector& detFD = detector.GetFDetector();
  unsigned int sumAllBins = 0;

  for (FDetector::EyeIterator iEye = detFD.EyesBegin();
       iEye != detFD.EyesEnd(); ++iEye) {

    const unsigned int eyeId = iEye->GetId();

    if (FDevent.HasEye(eyeId, fevt::ComponentSelector::eExists)) {
      ERROR("EYE already there!!!!");
    }
    FDevent.MakeEye(eyeId, fevt::ComponentSelector::eExists);
    fevt::Eye& eyeEvent = FDevent.GetEye(eyeId, fevt::ComponentSelector::eExists);

    for (fdet::Eye::TelescopeIterator iTel = iEye->TelescopesBegin();
         iTel != iEye->TelescopesEnd(); ++iTel) {

      const unsigned int telId = iTel->GetId();

      if (eyeEvent.HasTelescope(telId, fevt::ComponentSelector::eExists)) {
        ERROR("TELESCOPE already there!!!!");
      }
      eyeEvent.MakeTelescope(telId, fevt::ComponentSelector::eExists);
      fevt::Telescope& telEvent = eyeEvent.GetTelescope(telId, fevt::ComponentSelector::eExists);

      /*
  for (fevt::FEvent::EyeIterator iEye = fevent.EyesBegin(fevt::ComponentSelector::eInDAQ);
       iEye != fevent.EyesEnd(fevt::ComponentSelector::eInDAQ) ; ++iEye) {

    if (iEye->GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;

    for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

      fevt::Telescope& telEvent = *iTel;
      unsigned int telId = telEvent.GetId();

      const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);
      */

      const fdet::Telescope& detTel = *iTel;
      const Point telPosition = detTel.GetPosition();
      const Vector d_initial = core - telPosition;

      ostringstream info;
      info << "Eye=" << eyeId << " Tel=" << telId;

      if (!telEvent.HasSimData())
        telEvent.MakeSimData();
      fevt::TelescopeSimData& telSim = telEvent.GetSimData();

      const Vector axis = detTel.GetAxis();

      /*
        Calculate tmin/tmax of shower in FOV of telescope
      */

      const double b1 = std::pow(d_initial * axis, 2);
      const double b2 = 2 * (d_initial * axis) * (axis * showerAxis);
      const double b3 = std::pow(axis * showerAxis, 2);
      const double a1 = d_initial * d_initial;
      const double a2 = 2 * (d_initial * showerAxis);
      const double a3 = 1; // == showerAxis * showerAxis;

      const double cosFOV = cos(detTel.GetCamera().GetFieldOfView());
      const double cosFOV2 = cosFOV * cosFOV;

      const double A = cosFOV2 * a3 - b3;
      const double B = cosFOV2 * a2 - b2;
      const double C = cosFOV2 * a1 - b1;

      // quadratic equation gets complex ?
      if (B*B - 4*A*C < 0) {
        //info << " - not in field of view";
        //INFO(info);
        continue;
      } // if quadratic equation is not solvable

      double distMin = (-B - sqrt(B*B - 4*A*C)) / (2*A);
      double distMax = (-B + sqrt(B*B - 4*A*C)) / (2*A);

      /* **
         cerr << "1: distMin " << setw (10) << distMin
         << " distMax " << distMax
         << " distGnd: " << distGround
         << endl;
         *** */

      // check order in time / distance
      if (distMin > distMax)
        swap(distMin, distMax);

      /* **
         cerr << "2: distMin " << setw (10) << distMin
         << " distMax " << distMax
         << " distGnd: " << distGround
         << endl;
         *** */

      //const Point p1 = showerInitialPos + distMin*showerAxis;
      //const Point p2 = showerInitialPos + distMax*showerAxis;
      const Point pMin = core + distMin*showerAxis;
      const Point pMax = core + distMax*showerAxis;

      // check for non-FOV shower parts
      const bool invalidPointMin = ((pMin - telPosition)*axis < 0);
      const bool invalidPointMax = ((pMax - telPosition)*axis < 0);

      if (invalidPointMin && invalidPointMax) {  // --> (case 1)

        //info << " - not in field of view" ;
        //INFO(info);
        continue;

      } else if (invalidPointMax) {         // --> (case 2)

        if (distMin > distStart) {               // --> (case 2a)

          distMax = distMin;
          distMin = distStart;

        } else {                       // --> (case 2b)

          //info << " - not in field of view";
          //INFO(info);
          continue;

        }

      } else if (invalidPointMin) {         // --> (case 3)

        if (distMax<distStop) {      // --> (case 3a)

          distMin = distMax;
          distMax = distStop;

        } else {                       // --> (case 3b)

          //info << " - not in field of view";
          //INFO(info);
          continue;

        }
      }

      /* **
         cerr << "3: distMin " << setw (10) << distMin
         << " distMax " << distMax
         << " distGnd: " << distGround
         << endl;
         *** */

      // check for points before profile start or below ground level
      distMin = max(distMin, distStart);
      distMax = min(distMax, distStop);

      /* **
         cerr << "4: distMin " << setw (10) << distMin
         << " distMax " << distMax
         << " distGnd: " << distGround
         << endl;
         *** */

      if (distMin >= distMax) {
        //info << " - not in field of view ";
        //INFO(info);
        continue;
      }

      // here we finally found a part of the shower inside FOV !!!!!

      //const Point pointMin = showerInitialPos + distMin*showerAxis;
      //const Point pointMax = showerInitialPos + distMax*showerAxis;
      const Point pointMin = core + distMin*showerAxis;
      const Point pointMax = core + distMax*showerAxis;

      //tMinAtDia = ((showerInitialPos - pointMin).GetMag() +
      double tMinAtDia = (distMin + (pointMin - telPosition).GetMag()) / kSpeedOfLight;

      /* **
         cerr << "tMinAtDia " << tMinAtDia
         << " dt " << (tMinAtDia - tMinAtDia_) / ns
         << " ds " << (tMinAtDia - tMinAtDia_) * kSpeedOfLight / m
         << " distMin " << distMin
         << endl;
         *** */

      //tMaxAtDia = ((showerInitialPos - pointMax).GetMag() +
      double tMaxAtDia = (distMax + (pointMax - telPosition).GetMag()) / kSpeedOfLight;

      /* **
         cerr << "tMaxAtDia " <<  tMaxAtDia
         << " dt " << (tMaxAtDia - tMaxAtDia_) / ns
         << " ds " << (tMaxAtDia - tMaxAtDia_) * kSpeedOfLight / m
         << " distMax " << distMax
         << endl;
         *** */

      // time must move FORWARD
      if (tMinAtDia > tMaxAtDia)
        swap(tMinAtDia, tMaxAtDia);

      /*
        Found tMinAtDia/tMaxAtDia -> initialize simteldata
      */

      const double deltaT = tMaxAtDia - tMinAtDia;
      unsigned int numbins = deltaT / (detTel.GetCamera().GetFADCBinSize() / fBinning);

      if (numbins < 2) {
        info << " - track too short: deltaT/ns=" << deltaT/ns;
        INFO(info);
        continue;
      }

      const double minimalDist = min((telPosition - pointMin).GetMag(),
                                     (telPosition - pointMax).GetMag());
      double deltaD = fabs(distMax - distMin) / numbins;
      double deltaPhi = atan(deltaD / minimalDist);

      // RU Mi 8. Mär 11:56:23 CET 2017 It is actually only the
      // angular bin-width that matters for FD telescopes. The
      // previous logic was producing either to fine (bad
      // performance), or to sparse binning (bad physics).

      if (deltaPhi>fMaxAngleBinWidth) { // FOR VERY CLOSE SHOWERS
        deltaPhi = fMaxAngleBinWidth;
        deltaD = tan(deltaPhi) * minimalDist;
        numbins = fabs(distMax-distMin) / deltaD;
      }

      const double tracebin = deltaT / numbins;

      info << resetiosflags(std::ios::scientific) << setiosflags(std::ios::fixed);
      info //<< " : tMinAtDia=" << tMinAtDia << " tMaxAtDia=" << tMaxAtDia << ", "
        << ", nBins=" << setw(4) << numbins << ", dT=" << setw(5) << setprecision(1) << tracebin/ns << " ns, "
        //<< " dMin=" << distMin/m << " dMax=" << distMax/m << ", "
        << "dDist=" << setw(7) << setprecision(2) << deltaD/m << " m, "
        << "dAngle=" << setw(5) << setprecision(3) << deltaPhi/deg << " deg";
      //info << "; dmin=" << distMin << " dmax=" << distMax << " tmin=" << tMinAtDia << " tmax=" << tMaxAtDia << " mdist=" << minimalDist << " tMinAtDia=" << tMinAtDia;
      INFO(info);

      //TimeStamp traceAtDiaTime = timeAtStart + TimeInterval(tMinAtDia);
      TimeStamp traceAtDiaTime = timeAtCore + TimeInterval(tMinAtDia);
      telSim.SetPhotonsStartTime(traceAtDiaTime);

      // create event-structure
      // TRACE POSITION --------------------------------------------
      if (!telSim.HasDistanceTrace())
        telSim.MakeDistanceTrace(numbins, tracebin);

      sumAllBins += numbins;

    } // tel loop

  } // eye loop

  ostringstream info;
  info << "Looping over " << sumAllBins << " bins in all telescopes combined";
  INFO(info);

  return eSuccess;
}


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