RdTopDownStationSelector.cc

Go to the documentation of this file.

#include "RdTopDownStationSelector.h"

#include <fwk/CentralConfig.h>

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

#include <utl/Branch.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/PhysicalConstants.h>
#include <utl/Vector.h>
#include <utl/Point.h>

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

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

#include "TMath.h"


using namespace std;
using namespace fwk;


namespace RdTopDownStationSelector {

  VModule::ResultFlag
  RdTopDownStationSelector::Init()
  {
    utl::Branch topB = CentralConfig::GetInstance()->GetTopBranch("RdTopDownStationSelector");
    topB.GetChild("InfoLevel").GetData(fInfoLevel);
    topB.GetChild("maxNumberOfStations").GetData(fMaxNumberOfStations);
    topB.GetChild("minNumberOfStations").GetData(fInitialNumberOfSignalStations);
    topB.GetChild("stopAtFirstStation").GetData(fStopAtFirstStation);
    topB.GetChild("chiSquareProbabilityCutValue").GetData(fChiSquareProbabilityCutValue);
    topB.GetChild("initialSDCheck").GetData(fInitialSDCheck);
    topB.GetChild("maxIncomingDirectionDifference").GetData(fMaxIncomingDirectionDifference);
    topB.GetChild("additionalSystematicError").GetData(fAdditionalSystematicError);

    return eSuccess;
  }


  bool
  ComparePair(const std::pair<int, double>& a, const std::pair<int, double>& b)
  {
    return a.second < b.second;
  }


  VModule::ResultFlag
  RdTopDownStationSelector::Run(evt::Event& event)
  {
    if (!event.HasREvent()) {
      WARNING("No REvent");
      return eBreakLoop;
    }

    ostringstream info;
    info << "RUN in iteration " << fIteration;
    INFODebug(info);

    revt::REvent& rEvent = event.GetREvent();
    evt::ShowerRRecData& rShower = event.GetRecShower().GetRRecShower();

    const det::Detector& detector = det::Detector::GetInstance();
    const rdet::RDetector& rDetector = detector.GetRDetector();

    const utl::Vector showerAxis = rShower.GetReferenceAxis(event);
    const utl::Point showerCore = rShower.GetReferenceCorePosition(event);

    const int runId = rEvent.GetHeader().GetRunNumber();
    const int eventId = rEvent.GetHeader().GetId();
    const string headerId = event.GetHeader().GetId();

    // check if new event comes in and reset count variables
    if ((fCurrentRunId != runId) || (fCurrentEventId != eventId) || (fCurrentHeaderId != headerId)) {
      fIteration = 0;
      fCurrentRunId = runId;
      fCurrentEventId = eventId;
      fCurrentHeaderId = headerId;
      fChiSquaredNDF = 1;
      fCurrentNumberOfStations = fInitialNumberOfSignalStations;
      INFOIntermediate("A new event comes in, setting fIteration = 0.");
    }

    // check if last added station should stay in signal station list
    if (fIteration > 0) {
      if (fCurrentNumberOfStations >= 4) {
        if (!rShower.HasParameter(revt::eFitSuccess)){
          WARNING("Event has no eFitSuccess parameter. Since this is not the first iteration"
                  " (at least I think so) it should! This can happen if events with identical IDs are"
                  " reconstructed, e.g. using the Offline file format and <loop numTimes='>1'...");
        }

        if (rShower.GetParameter(revt::eFitSuccess) == 0) {
          rEvent.GetStation(fLastAddedStationId).SetRejectedReason(revt::eLargeTimeResidual);
          info.str("");
          info << "Wave fit was not successful with station " << fLastAddedStationId
               << ". Station will be rejected permanently.";
          INFOIntermediate(info);

          if (fStopAtFirstStation) {
            return eBreakLoop;
          }
        } else {
          // calculate chi2 again with additional uncertainty per station to account
          // for model uncertainties (e.g. the true shower front is not a plane wave
          double chi2 = 0;
          INFODebug("calculating chi2 manually: ");
          for (const auto& station : rEvent.SignalStationsRange()) {
            const double timeResidual = station.GetRecData().GetParameter(revt::eTimeResidual);
            const double timeResidualError = station.GetRecData().GetParameterError(revt::eTimeResidual);
            const double sigma_t2 = utl::Sqr(fAdditionalSystematicError) + utl::Sqr(timeResidualError);
            chi2 += utl::Sqr(timeResidual) / sigma_t2;

            info << "tres = " << timeResidual / utl::ns << "ns, error = "
                 << std::sqrt(sigma_t2) / utl::ns << "ns" << std::endl;
            INFODebug(info);
          }
          const double ndf = max(rShower.GetParameter(revt::eWaveFitNDF), 1.);
          const double chi2ndf = chi2 / ndf;
          info << "chi2 = " << chi2 << "/" << ndf << " = " << chi2ndf <<  std::endl;
          INFODebug(info);

          // if chi2 per ndf is considerable larger than without last added station, remove last added station
          if ((chi2ndf > 2 * fChiSquaredNDF) && (chi2ndf > 1)) {
            info.str("");
            info << "This condition is not used, but for your information: Adding Station "
                 << fLastAddedStationId << " to the fit results in a chi2/ndf = " << chi2ndf
                 << ". Previously it was " << fChiSquaredNDF
                 << ". Station will be rejected permanently.";
            INFODebug(info);
          }

          bool rejectChi2Prob = false;
          const double chi2prob = TMath::Prob(chi2, ndf);

          if (chi2prob < fChiSquareProbabilityCutValue) {
            rejectChi2Prob = true;
            info.str("");
            info << "Adding Station " << fLastAddedStationId
                 << " to the fit results in a chi2probability  = " << chi2prob
                 << ". This is smaller than the allowed probability of "
                 << fChiSquareProbabilityCutValue << ". Station will be rejected permanently.";
            INFOIntermediate(info);
          }

          // cut using preferred algorithm (currently chi2prob)
          if (rejectChi2Prob) {
            rEvent.GetStation(fLastAddedStationId).SetRejectedReason(revt::eLargeTimeResidual);
            if (fStopAtFirstStation) {
              return eBreakLoop;
            }
          } else {
            fCurrentNumberOfStations++;
            fChiSquaredNDF = chi2ndf;
            info.str("");
            info << "Adding Station " << fLastAddedStationId
                 << " does not increase chi2ndf significantly (chi2probability = " << chi2prob
                 << "). Station will be kept. In the next step " << fCurrentNumberOfStations
                 << " stations will be considered.";
            INFOIntermediate(info);
          }
        }
      } else if (fInitialSDCheck && (fCurrentNumberOfStations == 3)) {
        // to make sure that none of the closest 3 stations is a noise pulse, we check if the
        // incoming direction reconstructed from them roughly agrees with SD
        if (fCurrentNumberOfStations == 3) {
          // this variable is set to true if we find a reason to suspect that one of the stations is faulty
          bool faultyStationPresent = false;
          if (rShower.GetParameter(revt::eFitSuccess) == 0) {
            faultyStationPresent = true;
            info.str("");
            info << "Fit was not successful, suggesting that there is something wrong with one of the 3 stations"
                    " closest to the axis. Need to investigate further...";
            INFOIntermediate(info);
          } else {
            const utl::Vector referenceAxis = rShower.GetReferenceAxis(event);
            const double dotProduct = showerAxis * referenceAxis;
            const double angularDifference = acos(dotProduct / showerAxis.GetMag() / referenceAxis.GetMag());

            if (angularDifference > fMaxIncomingDirectionDifference) {
              faultyStationPresent = true;
              info.str("");
              info << "SD and RD reconstruction differ by "
                   << fMaxIncomingDirectionDifference / utl::degree
                   << " degrees. Need to investigate further...";
              INFOIntermediate(info);
            }
          }

          if (faultyStationPresent) {
            // One of the stations if faulty. We now compare measured signal time with the expectation
            // from SD to find out which one
            const utl::Point& referencePosition = rShower.GetReferenceCorePosition(event);
            const utl::Vector& referenceAxis = rShower.GetReferenceAxis(event);

            std::vector<double> expectedSignalTime;
            std::vector<double> measuredSignalTime;
            std::vector<int> stationIds;

            for (auto& station : rEvent.SignalStationsRange()) {
              measuredSignalTime.push_back(station.GetRecData().GetParameter(revt::eSignalTime));
              const int stationId = station.GetId();
              const utl::Vector stationPosition = rDetector.GetStation(stationId).GetPosition() -
                referencePosition;
              expectedSignalTime.push_back(-(referenceAxis * stationPosition / utl::kSpeedOfLight));
              stationIds.push_back(stationId);
            }

            const double avgExpectedSignalTime = TMath::Mean(expectedSignalTime.begin(),
                                                             expectedSignalTime.end());
            const double avgMeasuredSignalTime = TMath::Mean(measuredSignalTime.begin(),
                                                             measuredSignalTime.end());
            int faultyStation = 0;
            double maxTimeResidual = 0;
            for (unsigned int i = 0; i < measuredSignalTime.size(); i++) {
              const double timeResidual = abs(measuredSignalTime[i] - avgMeasuredSignalTime -
                                              expectedSignalTime[i] + avgExpectedSignalTime);
              if (timeResidual > maxTimeResidual) {
                maxTimeResidual = timeResidual;
                faultyStation = stationIds[i];
              }
            }
            rEvent.GetStation(faultyStation).SetRejectedReason(revt::eLargeTimeResidual);
          } else {
            ++fCurrentNumberOfStations;
          }
        }
      }
    }

    // break loop if the maximum number of signal stations is reached
    if (fMaxNumberOfStations < fCurrentNumberOfStations) {
      info.str("");
      info << "Maximum number of signal stations reached, breaking loop... ("
           << "current number of signal stations is " << fCurrentNumberOfStations
           << " and maximum number of stations is " << fMaxNumberOfStations << ")";
      INFOIntermediate(info);
      return eBreakLoop;
    }

    // fill distances of all suitable stations into std::vector
    // (suitable means: a signal station or a ONLY TopDown deselcted station)
    // Here all temporarily rejected stations are reseted, so vIdDistance contains all possible signal stations
    // (i.e. all signal stations that have not been rejected permanently)
    // at the end of the Main() method of the module, stations get rejected temporarily again.
    std::vector<std::pair<int, double> > vIdDistance;

    info.str("");
    info << "List of suitable stations:";
    for (auto& station : rEvent.StationsRange()) {
      // select station only if it has a signal (and therefore is not rejected)
      // or if the only rejection reason is eTopDownDeselected
      if (station.HasSignal() || station.GetRejectedReason() == revt::eTopDownDeselected){
        info << " " << station.GetId() << endl;

        // set all stations back to signal stations
        // (it just removes the eTopDownDeselected status from the station)
        station.SetSignal();
        station.ClearRejectedReason();

        const utl::Point stationPosition = rDetector.GetStation(station).GetPosition();
        const double distance = utl::RadioGeometryUtilities::GetDistanceToAxis(showerAxis, showerCore,
                                                                               stationPosition);
        vIdDistance.emplace_back(station.GetId(), distance);
      }
    }
    INFODebug(info);
    sort(vIdDistance.begin(), vIdDistance.end(), ComparePair);  // sort: increasing distances to shower axis

    // break selection process if all available stations are considered
    if (fCurrentNumberOfStations > vIdDistance.size()) {
      info.str("");
      info << "Number of available signal stations is " << vIdDistance.size()
           << " and current number of used stations is " << fCurrentNumberOfStations;
      INFOIntermediate(info);
      INFO("All available stations are considered, breaking loop.");
      return eBreakLoop;
    }

    // use only the N (fCurrentNumberOfStations) nearest stations to the shower axis
    // reject all stations temporarily but the nearest fCurrentNumberOfStations stations.
    // (This counting variable will be increased if the added station does not decrease fit quality)
    for (unsigned int i = fCurrentNumberOfStations; i < vIdDistance.size(); ++i) {
      const int stationId = vIdDistance[i].first;
      rEvent.GetStation(stationId).SetRejectedReason(revt::eTopDownDeselected);
      info.str("");
      info << "station " << stationId << " is temporarily rejected (distance to shower axis = "
           << vIdDistance[i].second / utl::m << "m)";
      INFODebug(info);
    }
    fLastAddedStationId = vIdDistance[fCurrentNumberOfStations - 1].first;
    fIteration++;

    return eSuccess;
  }


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

}