RdStationSimulationRejector.cc

Go to the documentation of this file.

#include "RdStationSimulationRejector.h"

#include <fwk/CentralConfig.h>

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

#include <utl/Branch.h>
#include <utl/ErrorLogger.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystem.h>
#include <utl/Point.h>

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

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

#include <sevt/SEvent.h>
#include <sevt/StationTriggerData.h>


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


namespace RdStationSimulationRejector {

  VModule::ResultFlag
  RdStationSimulationRejector::Init()
  {
    const Branch topB = CentralConfig::GetInstance()->GetTopBranch("RdStationSimulationRejector");
    topB.GetChild("InfoLevel").GetData(fInfoLevel);
    topB.GetChild("RejectStationAffectedByThinning").GetData(fRejectStationAffectedByThinning);
    topB.GetChild("RemoveNonTriggeredRdStationsFromEvent").GetData(fRemoveNonTriggeredRdStationsFromEvent);

    ostringstream info;
    info << "Module config:"
         << "\n\tReject station affected by thinning: " << fRejectStationAffectedByThinning
         << "\n\tReject RD station when WCD is not triggered: " << fRemoveNonTriggeredRdStationsFromEvent;
    INFOFinal(info);

    return eSuccess;
  }


  VModule::ResultFlag
  RdStationSimulationRejector::Run(Event& event)
  {
    if (!event.HasREvent())
      return eSuccess;

    REvent& rEvent = event.GetREvent();
    const Detector& detector = Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();

    unsigned int numberOfNotTriggeredRDStations = 0;
    unsigned int numberOfAcceptedStations = 0;
    ostringstream info;

    for (auto& station : rEvent.StationsRange()) {
      const int rStationID = station.GetId();

      // For the Rd we will just get data for stations with a trigger in the WCD.
      // If corresponding WCD did not trigger do not add Rd station to REvent.
      if (rStationID > rDetector.GetRdSdStationIdLink() && fRemoveNonTriggeredRdStationsFromEvent) {
        const int sStationID = rStationID - rDetector.GetRdSdStationIdLink();

        const sevt::SEvent& sEvent = event.GetSEvent();
        if (!sEvent.HasStation(sStationID) || !sEvent.GetStation(sStationID).HasTriggerData()) {
          info.str("");
          info << "Not add RD station " << rStationID
              << ". No corresponding SD station in has trigger data.";
          INFODebug(info);
          ++numberOfNotTriggeredRDStations;
          station.SetExcludedReason(eStationWithoutTrigger);
          continue;
        }

        const sevt::StationTriggerData& trig = sEvent.GetStation(sStationID).GetTriggerData();
        if (!trig.IsT2() && !trig.IsT1()) {  // This is not equal to station.IsSilent() but this is desired here
          info.str("");
          info << "Not add RD station " << rStationID
              << ". No corresponding SD station in has trigger data.";
          INFODebug(info);
          ++numberOfNotTriggeredRDStations;
          station.SetExcludedReason(eStationWithoutTrigger);
          continue;
        }
      }

      if (fRejectStationAffectedByThinning) {
        const ShowerSimData& simshow = event.GetSimShower();
        const RadioSimulation& rsim = simshow.GetRadioSimulation();
        const CoordinateSystemPtr corsys = simshow.GetLocalCoordinateSystem();

        // exclude (most) simulated pulses (stations) which are affected by thinning
        const Vector refShowerAxis = -1 * simshow.GetDirection();  // sim shower direction defined in particle direction
        const double zenith = refShowerAxis.GetTheta(corsys);

        // not optimized for below
        if (zenith > 65*deg) {
          const Point& simCore = simshow.GetPosition();
          const double maximumAllowedDistance = RdStationSimulationRejector::thinning_distance(zenith);

          const Point Position = rDetector.GetStation(station).GetPosition();
          const double axisDistance = RadioGeometryUtilities::GetDistanceToAxis(refShowerAxis, simCore, Position);

          // not quite correct since the early late correction should be included
          if (axisDistance > maximumAllowedDistance) {
            info.str("");
            info << "Station " << rStationID << " excluded because of thinning:";
            INFOFinal(info);
            station.SetExcludedReason(eThinningArtefacts);
            continue;
          } else {
            numberOfAcceptedStations++;
          }
        }

        if (numberOfAcceptedStations / rsim.GetNumPulses() < 0.3) {
          WARNING("A lot of stations have been excluded due to the thinning cut.\n"
                  "Please check if the core definition (e.g. SimShowerData or RadioSimulation core) is correct!");
        }
      }
    }

    if (numberOfNotTriggeredRDStations) {
      info.str("");
      info << numberOfNotTriggeredRDStations << " RD stations have not been triggered.";
      INFO(info);
    }

    return eSuccess;
  }


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

}