RdEventPostSelector.cc

Go to the documentation of this file.

#include "RdEventPostSelector.h"

#include <fwk/CentralConfig.h>

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

#include <det/Detector.h>
#include <rdet/RDetector.h>

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

#include <sevt/SEvent.h>
#include <revt/REvent.h>
#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>

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


namespace RdEventPostSelector {

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

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("MaxZenithCut").GetData(fMaxZenith);
    topBranch.GetChild("UseConverged").GetData(fUseConverged);
    topBranch.GetChild("UseSuccessfulFits").GetData(fUseSuccessfulFits);
    topBranch.GetChild("MinNumberOfStationsWithPulseFound").GetData(
        fMinNumberOfStationsWithPulseFound);
    topBranch.GetChild("UseCoincidenceCheck").GetData(fUseCoincidenceCheck);
    topBranch.GetChild("MaximumAngularDeviation").GetData(fMaximumAngularDeviation);
    topBranch.GetChild("MaximumAxisToCoreDistance").GetData(fMaximumAxisToCoreDistance);
    topBranch.GetChild("MaximumNumberOfRejectedStations").GetData(fMaximumNumberOfRejectedStations);
    topBranch.GetChild("AERA24Distance").GetData(fAERA24Distance);

    string tmp = topBranch.GetChild("UsedEye").Get<string>();
    if (tmp == "CO")
      fWhichEye = 4;
    if (tmp == "HE")
      fWhichEye = 5;
    if (tmp == "HECO")
      fWhichEye = 6;

    ostringstream info;
    info.str("");

    if (fMaxZenith < 0)
      info << " \n No zenith cut apllied";
    else
      info << " \n Maximum zenith angle: " << fMaxZenith / degree << " deg";
    if (fUseConverged)
      info << " \n Use only converged shower reconstructions";
    if (fUseSuccessfulFits)
      info << " \n Use only successfully fitted shower reconstructions";
    if (fMinNumberOfStationsWithPulseFound > 0)
      info << " \n Requiring " << fMinNumberOfStationsWithPulseFound << " stations with pulse found";
    if (fUseCoincidenceCheck == "SD" || fUseCoincidenceCheck == "FD") {
      info << " \n Checking for coincidences with SD/FD with the following parameter:";
      info << " \n Maximum angular deviation: " << fMaximumAngularDeviation / degree << " degree.";
      info << " \n Maximum axis to core distance: " << fMaximumAxisToCoreDistance / kilometer << " km.";
    }
    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdEventPostSelector::Run(Event& event)
  {
    // Check if there is an radio event at all
    if (!event.HasREvent()) {
      INFOFinal("No radio event found!");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }
    const REvent& rEvent = event.GetREvent();

    // Check for SD/FD event if coincidence check is active
    if (!event.HasSEvent() && fUseCoincidenceCheck == "SD") {
      INFOFinal("No SD event found, but needed for coincidence check. Skipping event.");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    } 

    if (!event.HasFEvent() && fUseCoincidenceCheck == "FD") {
      INFOFinal("No FD event found, but needed for coincidence check. Skipping event.");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }

    if (!event.HasRecShower()) {
      INFOFinal("Event has no RRecShower!");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }
    const ShowerRecData& recShower = event.GetRecShower();

    if (!recShower.HasRRecShower()) {
      INFOFinal("Event has no RecShower!");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }
    const ShowerRRecData& rRecShower = recShower.GetRRecShower();

    ostringstream info;

    if (fMaxZenith > 0) {
      const double thisZenith = rRecShower.GetZenith();
      if (thisZenith > fMaxZenith) {
        INFOFinal("Shower zenith is above cut value. Skipping event");
        ++fNumberOfRejectedEvents;
        ++fNumberOfRejectedEventsZenithCut;
        return eContinueLoop;
      }
    }

    if (fMaximumNumberOfRejectedStations != -1) {
      //Note: manually rejected stations (i.e. rejectionStatus=1) are not taken into account
      const int nRejectedStations = 
        count_if(rEvent.RejectedStationsBegin(), rEvent.RejectedStationsEnd(),
        [&](auto const& station){ return station.GetRejectedReason() > 1; });

      if (nRejectedStations > fMaximumNumberOfRejectedStations) {
        INFOFinal("Too many rejected stations. Skipping event.");
        ++fNumberOfRejectedEvents;
        return eContinueLoop;
      }
    }

    if (fMinNumberOfStationsWithPulseFound != 0) {
      const int numberOfStations = rEvent.GetNumberOfSignalStations();

      if (numberOfStations < fMinNumberOfStationsWithPulseFound) {
        INFOFinal("Event has not enough stations with pulse found. Skipping event.");
        info.str("");
        info << " Pulse found in " << numberOfStations << " stations, requested are "
             << fMinNumberOfStationsWithPulseFound << " stations";
        INFOIntermediate(info);
        ++fNumberOfRejectedEvents;
        ++fNumberOfRejectedEventsPulseFound;
        return eContinueLoop;
      }
    }

    // RdDirectionConvergenceChecker
    if (rRecShower.HasParameter(eFitConvergence)) {
      if (!rRecShower.GetParameter(eFitConvergence) && fUseConverged) {
        INFOFinal("Reconstruction did not converge. Skipping event.");
        ++fNumberOfRejectedEvents;
        return eContinueLoop;
      }
    } else if (!rRecShower.HasParameter(eFitConvergence) && fUseConverged) {
      INFOFinal("Parameter eFitConvergence not set. Skipping event.");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }

    // succesful wavefront fit
    if (rRecShower.HasParameter(eFitSuccess)) {
      if (!rRecShower.GetParameter(eFitSuccess) && fUseSuccessfulFits) {
        INFOFinal("Reconstruction fit was not successful. Skipping event.");
        ++fNumberOfRejectedEvents;
        return eContinueLoop;
      }
    } else if (!rRecShower.HasParameter(eFitSuccess) && fUseSuccessfulFits) {
      INFOFinal("Parameter eFitSuccess not set. Skipping event.");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }

    // coindidencxe checks
    if (fUseCoincidenceCheck == "SD" && event.GetRecShower().HasSRecShower()) {
      const ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower();
      const Vector& rShowerAxis = rRecShower.GetAxis();
      const Vector& sShowerAxis = sRecShower.GetAxis();

      if (!rRecShower.HasCorePosition()) {
        INFOFinal("no radio core reconstructed. Event is skipped.");
        ++fNumberOfRejectedEvents;
        return eContinueLoop;
      }

      const Point& rCore = rRecShower.GetCorePosition();
      const Point& sCore = sRecShower.GetCorePosition();

      const bool passes = passesCoincidenceCheck(rCore, rShowerAxis, sCore, sShowerAxis);
      if (!passes) {
        ++fNumberOfRejectedEvents;
        return eContinueLoop;
      }
    }

    if (fUseCoincidenceCheck == "FD") {
      FEvent& fdEvent = event.GetFEvent();

      for (auto eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
           eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIter) {
        if (eyeIter->GetId() != fWhichEye)  // continue until selected eye
          continue;

        const Eye& eye = *eyeIter;
        if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
          ERROR("Eye has no RecData or FRecShower but FDCore specified in the settings!");
          ++fNumberOfRejectedEvents;
          return eContinueLoop;
        }

        const ShowerFRecData& fRecShower = eye.GetRecData().GetFRecShower();
        const Vector& rShowerAxis = rRecShower.GetAxis();
        const Vector& fShowerAxis = fRecShower.GetAxis();

        if (!rRecShower.HasCorePosition()) {
          INFOFinal("no radio core reconstructed. Event is skipped");
          ++fNumberOfRejectedEvents;
          return eContinueLoop;
        }

        const Point& rCore = rRecShower.GetCorePosition();
        const Point& fCore = fRecShower.GetCorePosition();

        const bool passes = passesCoincidenceCheck(rCore, rShowerAxis, fCore, fShowerAxis);
        if (!passes) {
          ++fNumberOfRejectedEvents;
          return eContinueLoop;
        }
      }
    }

    // get min distance of reference core and AERA24 stations
    double minDistance = 100 * km;
    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    for (const auto& station : rEvent.SignalStationsRange()) {
      if(station.GetId() > 124) {  // skip all not AERA24 stations
        continue;
      }

      const Point sPos = rDetector.GetStation(station).GetPosition();
      const double distance = (sPos - rRecShower.GetReferenceCorePosition(event)).GetMag();
      if(distance < minDistance) {
        minDistance = distance;
      }
    }

    if (minDistance > fAERA24Distance) {
      INFOFinal("Core position too far away from AERA24 sub array. Skipping event.");
      ++fNumberOfRejectedEvents;
      return eContinueLoop;
    }

    // if all cuts are passed, return eSuccess at the end
    return eSuccess;
  }


  VModule::ResultFlag
  RdEventPostSelector::Finish()
  {
    ostringstream sstr;
    sstr << "\n\t" << fNumberOfRejectedEventsZenithCut
         << " events rejected because of zenith angle cut (zenith was larger than "
         << fMaxZenith / degree << " deg)"
         << "\n\t" << fNumberOfRejectedEventsPulseFound
         << " events rejected because of not enough stations with pulse found"
            " (number of stations with pulse found was smaller than "
         << fMinNumberOfStationsWithPulseFound << ")"
         << "\n\t" << fNumberOfRejectedEventsAngularDeviation
         << " events rejected because of angular deviation cut"
            " (angular deviation SD/RD or FD/RD was larger than "
         << round(fMaximumAngularDeviation / degree * 10.) / 10. << " deg)"
         << "\n\t" << fNumberOfRejectedEventsAxisToCoreDistance
         << " events rejected because of axis to core distance cut"
            " (axis to core distance SD/RD or FD/RD was larger than "
         << fMaximumAxisToCoreDistance << " m)"
         << "\n\t" << fNumberOfRejectedEvents << " events rejected in total";
    INFO(sstr);

    return eSuccess;
  }


  bool
  RdEventPostSelector::passesCoincidenceCheck(const Point& rCore, const Vector& rAxis,
                                              const Point& oCore, const Vector& oAxis)
  {
    ostringstream info;

    if (fMaximumAngularDeviation >= 0) {
      const double angleDiff = acos(Normalized(rAxis) * Normalized(oAxis));
      info.str("");
      info << "Angular deviation of " << round((angleDiff / degree) * 10.) / 10
           << " degree between RD and other axis ";

      if (angleDiff / degree >= fMaximumAngularDeviation / degree) {
        info << "is larger than allowed difference of "
             << round((fMaximumAngularDeviation / degree) * 10.) / 10
             << " degree. Skipping event.";
        INFOFinal(info);
        ++fNumberOfRejectedEventsAngularDeviation;
        return false;
      } else {
        info << "is in allowed range of maximal "
             << round((fMaximumAngularDeviation / degree) * 10) / 10
             << " degree. Event is coincident candidate";
        INFOIntermediate(info);
      }
    }

    if (fMaximumAxisToCoreDistance >= 0) {
      const double axisToCoreDistance = RadioGeometryUtilities::GetDistanceToAxis(oAxis, oCore, rCore);
      info.str("");
      info << "Core distance of " << round((axisToCoreDistance / kilometer) * 100.) / 100 << " km ";

      if (axisToCoreDistance / kilometer >= fMaximumAxisToCoreDistance / kilometer) {
        info << "is larger than allowed difference of "
             << round((fMaximumAxisToCoreDistance / kilometer) * 100.) / 100
             << " km. Skipping event.";
        INFOFinal(info);
        ++fNumberOfRejectedEventsAxisToCoreDistance;
        return false;
      } else {
        info << "is in allowed range of "
             << round((fMaximumAxisToCoreDistance / kilometer) * 100) / 100
             << " km. Event is coincident candidate.";
        INFOIntermediate(info);
      }
    }

    // if comes here, event meets all coincidence criteria
    INFOIntermediate("Event meets all coincidence criteria.");
    return true;
  }

}