RdDirectionConvergenceChecker.cc

Go to the documentation of this file.

#include "RdDirectionConvergenceChecker.h"

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>

#include <revt/REvent.h>
#include <revt/Header.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>


using namespace std;
using namespace fwk;
using namespace utl;
using namespace revt;
using namespace evt;


namespace RdDirectionConvergenceChecker {


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

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("AngleEpsilon").GetData(fAngleEpsilon);
    topBranch.GetChild("MaxIterations").GetData(fMaxIterations);
    topBranch.GetChild("CheckMirrorZenith").GetData(fCheckMirrorZenith);

    return eSuccess;
  }


  VModule::ResultFlag
  RdDirectionConvergenceChecker::Run(evt::Event& event)
  {
    // nothing to do if there is no REvent
    if (!event.HasREvent()) {
      INFOFinal("eContinueLoop: No REvent, so no reconstruction");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();
    ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

    const int eventId = rEvent.GetHeader().GetId();
    ostringstream info;

    // Check if we have a new event and set the iteration
    if (!showerrrec.HasParameter(eNumberOfIterationsInDirectionFit)) {
      info.str("");
      info << "New event with ID: " << eventId;
      INFODebug(info);

      // Reset the vectors
      fAllAzimuths.clear();
      fAllZeniths.clear();

      showerrrec.SetParameter(eNumberOfIterationsInDirectionFit, 0, false);
    }

    // read out, increase and write back iteration counter
    unsigned int currentIteration = showerrrec.GetParameter(eNumberOfIterationsInDirectionFit) + 1;
    showerrrec.SetParameter(eNumberOfIterationsInDirectionFit, currentIteration, false);

    if (currentIteration == 1) {
      /* There is no previous direction to compare against yet.
         If we have no reconstructed radio incident direction after one iteration, we break the loop,
         since then, the fit seems to fail allways. otherwise we just continue with the next iteration. */
      if (!showerrrec.HasAxis()) {
        INFO("No reconstructed shower axis was found in first iteration. Breaking the loop without axis!");
        showerrrec.SetParameter(eFitConvergence, false);
        return eBreakLoop;
      }

      // Save the current angle to the vector containing all reconstructed directions.
      fAllAzimuths.push_back(showerrrec.GetAzimuth());
      fAllZeniths.push_back(showerrrec.GetZenith());

      return eSuccess;
    }

    info.str("");
    info << "Currently in " << currentIteration << " iteration for event " << eventId;
    INFOIntermediate(info);

    // Check if the maximum number of iterations is reached
    if (currentIteration > fMaxIterations) {
      // if we tried mirrored zenith but did not find convergence, go back to not mirrored zenith
      if (showerrrec.HasParameter(eMirroredZenith)) {
        showerrrec.SetParameter(eShowerAxisX, ePreFitShowerAxisX, false);
        showerrrec.SetParameter(eShowerAxisY, ePreFitShowerAxisY, false);
        showerrrec.SetParameter(eShowerAxisZ, ePreFitShowerAxisZ, false);
        showerrrec.SetParameter(eNumberOfIterationsInDirectionFit, 1, false);

        fAllAzimuths.clear();
        fAllZeniths.clear();

        fAllAzimuths.push_back(showerrrec.GetAzimuth());
        fAllZeniths.push_back(showerrrec.GetZenith());

        INFOFinal("Iteration with mirrowed zenith did not converge. "
                  "Restart iteration with original axis of the pre-fit");
        return eSuccess;
      } else {
        INFO("Iteration limit reached without convergence! Killing Loop!");
        showerrrec.SetParameter(eFitConvergence, false);
        return eBreakLoop;
      }
    }

    // check if a previous iteration has failed
    if (showerrrec.HasParameter(eFitSuccess) && showerrrec.GetParameter(eFitSuccess) == false) {
      INFO("Previous fit failed, killing loop!");
      showerrrec.SetParameter(eFitConvergence, false);
      return eBreakLoop;
    }

    // fit worked, get the reconstructed direction and check for convergence
    const double curZen = showerrrec.GetZenith();
    const double curAzi = showerrrec.GetAzimuth();

    // Calculate difference between current and previous azimuth and zenith angles
    const double angleDiff = GetAngularDifference(curZen, curAzi, fAllZeniths.back(), fAllAzimuths.back());
    info.str("");
    info << "angular difference: " << angleDiff << " rad, " << angleDiff / utl::deg << " deg";
    INFOIntermediate(info);

    if (fInfoLevel >= eInfoDebug) {
      info.str("");
      info << "Current iteration: " << currentIteration << ", "
           << "last azimuth: " << fAllAzimuths.back() / deg << ", current one: " << curAzi / deg << ", "
           << "last zenith: " << fAllZeniths.back() / deg << ", curZen one: " << curZen / deg << ", "
           << "angular difference: " << angleDiff / deg;
      INFODebug(info);
    }

    // Check if the space angle difference is below the convergence limit
    if (abs(angleDiff) <= fAngleEpsilon){
      /* check if we would also have convergence for a down-going shower, far more probable
         but fitting can run into minimum for a up-going zenith (especially noisy events with few stations).
         Save information in eMirroredZenith. */
      if (fCheckMirrorZenith && curZen / deg > 90 && !showerrrec.HasParameter(eMirroredZenith)) {
        const auto oldZAxis = showerrrec.GetParameter(eShowerAxisZ);
        showerrrec.SetParameter(eShowerAxisZ, -oldZAxis, false);
        showerrrec.SetParameter(eMirroredZenith, true, true);
        showerrrec.SetParameter(eNumberOfIterationsInDirectionFit, 1, false);

        INFOFinal("Upwarding going shower detected, restart iteration with mirrowed zenith angle");
        return eSuccess;
      }

      INFOFinal("Convergence reached!");
      showerrrec.SetParameter(eFitConvergence, true);
      return eBreakLoop;
    }

    /* If we have more than 4 iterations, we do a check for alternating reconstructed directions.
       if the plane fit is triggering on two different peaks, reconstructing different but alternating directions
       in each iteration, we want to stop the process, since this is a signal for a noise dominated
       (or otherwise unphysical) reconstruction */
    if (fAllAzimuths.size() > 4 && fAllZeniths.size() > 4) {
      // copy last 4 elements to new vector and reverse order
      vector<double> latestAzi(fAllAzimuths.end() - 4, fAllAzimuths.end());
      vector<double> latestZen(fAllZeniths.end() - 4, fAllZeniths.end());

      reverse(latestAzi.begin() , latestAzi.end());
      reverse(latestZen.begin() , latestZen.end());

      info.str("");
      info << "Azimuth 0 - 1: " << (latestAzi[0] - latestAzi[1]) / degree << "\n"
           << "Azimuth 1 - 2: " << (latestAzi[1] - latestAzi[2]) / degree << "\n"
           << "Azimuth 2 - 3: " << (latestAzi[2] - latestAzi[3]) / degree << "\n"
           << "Zenith  0 - 1: " << (latestZen[0] - latestZen[1]) / degree << "\n"
           << "Zenith  1 - 2: " << (latestZen[1] - latestZen[2]) / degree << "\n"
           << "Zenith  2 - 3: " << (latestZen[2] - latestZen[3]) / degree << "\n";
      INFODebug(info);

      const double sa1 = GetAngularDifference(latestZen[0], latestAzi[0], latestZen[1], latestAzi[1]);
      const double sa2 = GetAngularDifference(latestZen[1], latestAzi[1], latestZen[2], latestAzi[2]);
      const double sa3 = GetAngularDifference(latestZen[2], latestAzi[2], latestZen[3], latestAzi[3]);

      if ((abs(abs(sa1) - abs(sa2)) < 0.001*degree) && (abs(abs(sa2) - abs(sa3)) < 0.001*degree ) ) {
        INFO("Alternation detected, breaking loop ...");
        showerrrec.SetParameter(eFitConvergence, false);
        return eBreakLoop;
      }
    }

    // Save the current angle to the vector containing all reconstructed directions.
    // start all over again
    fAllAzimuths.push_back(curAzi);
    fAllZeniths.push_back(curZen);

    return eSuccess;
  }


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


  double
  RdDirectionConvergenceChecker::GetAngularDifference(const double zen1, const double azim1,
                                                      const double zen2, const double azim2)
  {
    Vector3D curDirection;
    Vector3D lastDirection;

    curDirection = sin(zen1)*cos(azim1), sin(zen1)*sin(azim1), cos(zen1);
    lastDirection = sin(zen2)*cos(azim2), sin(zen2)*sin(azim2), cos(zen2);

    return Angle(curDirection, lastDirection);
  }

}