MdEventSelector.cc

Go to the documentation of this file.

#include "MdEventSelector.h"
#include <utl/ErrorLogger.h>
//
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
//
#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerMRecData.h>
#include <sevt/Station.h>
//
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/ConfigUtils.h>

using namespace std;
using namespace utl;
using namespace fwk;


namespace MdEventSelectorAG {

  MdEventSelector::MdEventSelector()
  { }


  MdEventSelector::~MdEventSelector()
  { }


  VModule::ResultFlag
  MdEventSelector::Init()
  {
    // Configuration reading
    utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdEventSelector");

    //Minimum Number of good stations to try the MLDF fit
    LoadConfig(config, "MinNumberOfCountersAtBeginning", fMinNumberOfCounters, 3);
    LoadConfig(config, "silentRange", fSilentRange, 5000);
    LoadConfig(config, "selectEventByTheta", fSelectEventByTheta, false);
    LoadConfig(config, "maxTheta", fMaxTheta, 90);

    LoadConfig(config, "T5NumberOfActiveStations", fT5NumberOfActiveStations, 6);
    LoadConfig(config, "RejectNonT5Events", fRejectNonT5Events, 0);
    LoadConfig(config, "FirstCrownMinDistance", fFirstCrownMinDistance, 600);
    LoadConfig(config, "FirstCrownMaxDistance", fFirstCrownMaxDistance, 1000);

    LoadConfig(config, "RejectSaturatedEvents", fRejectSaturatedEvents, 0);
    LoadConfig(config, "RejectDenseStations", fRejectDenseStations, 0);

    LoadConfig(config, "RejectTimeOutlierCounters", fRejectTimeOutlierCounters, 0);
    LoadConfig(config, "timeOutlierMuons", fTimeOutlierMuons, 10);
    LoadConfig(config, "timeOutlierDistance", fTimeOutlierDistance, 0);
    LoadConfig(config, "timeOutlierConstant", fTimeOutlierConstant, 0);
    LoadConfig(config, "timeOutlierSlope", fTimeOutlierSlope, 0);

    theMDetector = &(det::Detector::GetInstance().GetMDetector());
    return eSuccess;
  }


  VModule::ResultFlag
  MdEventSelector::Run(evt::Event& event)
  {
    //INFO("");

    ostringstream os;

    if (!event.HasMEvent()) {
      WARNING("Event without muon detector event: nothing to be done. Skipping to next event.");
      return eContinueLoop;
    }

    mevt::MEvent& me = event.GetMEvent();
    sevt::SEvent& se = event.GetSEvent();

    if (fRejectSaturatedEvents && me.IsSaturated()) {
      os << "Saturated event " << event.GetHeader().GetId() << " will not be reconstructed." << std::endl;
      WARNING(os);
      return eContinueLoop;
    }

    RejectFarSilents(me);

    if (fSelectEventByTheta && IsHighTheta(event)) {
      os.str("");
      os << "Skipping to next event. Zenith angle is greater than " << fMaxTheta/utl::deg << " degrees";
      INFO(os);
      return eContinueLoop;
    }

    RejectCloseSimDetectors(me);

    if (fRejectDenseStations)
      RejectDenseDetectors(me);

    if (fRejectTimeOutlierCounters) {
      RejectTimeOutliers(me, se);
    }

    if (!HasEnoughDetectors(me)) {
      INFO("Not enough detectors to reconstruct event, skipping...");
      return eContinueLoop;
    }

    // T5 processing
    const bool t5 = IsMEventT5(me);
    SetT5(event, t5);

    if (t5) {
      ++RunController::GetInstance().GetRunData().GetNamedCounters()["MdEventSelector/T5"];
    }

    if (fRejectNonT5Events && !t5) {
      INFO("Rejecting non T5 event");
      return eContinueLoop;
    }

    return eSuccess;
  }


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


  void
  MdEventSelector::SetT5(evt::Event& event, const bool t5)
    const
  {
    if (!event.HasRecShower())
      event.MakeRecShower();

    evt::ShowerRecData& shRec = event.GetRecShower();
    if (!shRec.HasMRecShower())
      shRec.MakeMRecShower();

    evt::ShowerMRecData& shMRec = shRec.GetMRecShower();
    shMRec.SetT5Trigger(t5);
  }


  bool
  MdEventSelector::HasEnoughDetectors(const mevt::MEvent& me)
    const
  {
    const unsigned int nCandidateCounters = me.GetNumberOfCandidateCounters();

    if (nCandidateCounters < fMinNumberOfCounters)
      return false;

    // Check for counters with muons
    unsigned int nCountersWithMuons = 0;
    for (mevt::MEvent::CounterConstIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {
      if (ic->IsCandidate() && !ic->IsEmpty()) {
        ++nCountersWithMuons;
      }
    }

    return nCountersWithMuons >= fMinNumberOfCounters;
  }


  void
  MdEventSelector::RejectFarSilents(mevt::MEvent& me)
    const
  {
    //INFO("");

    if (!me.GetNumberOfCandidateCounters())
      return;

    const mevt::Counter* hottestCounter = FindHottestCounter(me);

    const unsigned int hottestId = hottestCounter->GetId();
    utl::Point hottestPosition = theMDetector->GetCounter(hottestId).GetPosition();

    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {

      const unsigned int counterId = ic->GetId();
      const utl::Point counterPosition = theMDetector->GetCounter(counterId).GetPosition();

      const utl::Vector moduleToHottest = counterPosition-hottestPosition;
      const double distance = moduleToHottest.GetMag();

      if (ic->IsSilent() && distance > fSilentRange)
        ic->SetRejected("Far silent");

    } // end counters loop
  }


  const mevt::Counter*
  MdEventSelector::FindHottestCounter(const mevt::MEvent& me)
    const
  {
    return me.IsSaturated() ? FindHottestCounterSaturated(me) : FindHottestCounterUnsaturated(me);
  }


  const mevt::Counter*
  MdEventSelector::FindHottestCounterUnsaturated(const mevt::MEvent& me)
    const
  {
    //INFO("");

    double maxDensity = -1;
    const mevt::Counter* hottestCounter = nullptr;

    for (mevt::MEvent::CounterConstIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {

      //cout << "Counter " << ic->GetId() << " IsCandidate? " << std::boolalpha << ic->IsCandidate();<< endl;

      if (ic->IsCandidate() && !ic->IsSaturated()) {

        const double density = ic->GetMuonDensity();

        //cout << "Counter " << ic->GetId() << ", density = " << density << endl;

        if (density > maxDensity) {
          hottestCounter = &(*ic);
          maxDensity = density;
        }

      }

    }

    //cout << "Hottest counter " <<  hottestCounter->GetId() << endl;

    return hottestCounter;
  }


  const mevt::Counter*
  MdEventSelector::FindHottestCounterSaturated(const mevt::MEvent& me)
    const
  {
    //INFO("");

    unsigned int maxWindowsOn = 0;
    const mevt::Counter* hottestCounter = nullptr;

    for (mevt::MEvent::CounterConstIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {

      //cout << "Counter " << ic->GetId() << " IsCandidate? " << std::boolalpha << ic->IsCandidate() << endl;

      if (ic->IsCandidate() && ic->IsSaturated()) {

        const double windowsOn = ic->GetNumberOfChannelsOn();

        //cout << "Counter " << ic->GetId() << ", windows on = " << windowsOn << endl;

        if (windowsOn>maxWindowsOn) {
          hottestCounter = &(*ic);
          maxWindowsOn = windowsOn;
        }

      }

    }

    //cout << "Hottest counter " <<  hottestCounter->GetId() << endl;

    return hottestCounter;
  }


  bool
  MdEventSelector::IsHighTheta(const evt::Event& event)
    const
  {
    // nothing to do if there is no RecShower
    if (!(event.HasRecShower() && event.GetRecShower().HasSRecShower())) {
      ERROR("Current event does not have a recontructed shower.");
      return false;
    }

    const evt::ShowerSRecData& sRecShower = event.GetRecShower().GetSRecShower();

    // all these lines are required to find the zenith angle from the SD reconstruction
    const utl::Point& rCore = sRecShower.GetCorePosition();
    const utl::Vector& rAxis = sRecShower.GetAxis();
    const utl::CoordinateSystemPtr coreCS = fwk::LocalCoordinateSystem::Create(rCore);
    const double theta = rAxis.GetTheta(coreCS)/utl::deg;

    return theta > fMaxTheta/utl::deg;
  }


  void
  MdEventSelector::RejectCloseSimDetectors(mevt::MEvent& me)
    const
  {
    if (!me.GetNumberOfCandidateCounters())
      return;

    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {
      if (ic->HasSimData() && ic->GetSimData().IsInsideMinRadius())
        ic->SetRejected("Inside minimal radius");
    }
  }


  void
  MdEventSelector::RejectDenseDetectors(mevt::MEvent& me)
    const
  {
    if (!me.GetNumberOfCandidateCounters())
      return;

    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {
      const mdet::Counter& counter = theMDetector->GetCounter(ic->GetId());
      if (counter.IsDense())
        ic->SetRejected("In dense array");
    }
  }


  void
  MdEventSelector::RejectTimeOutliers(mevt::MEvent& me, const sevt::SEvent& se)
    const
  {
    ostringstream os;

    if (!me.GetNumberOfCandidateCounters())
      return;

    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {

      if (ic->IsRejected() || !ic->HasT50())
        continue;

      const utl::TimeStamp mdT50 = ic->GetSignalT50();

      // Get the signal start time of the partner sd station
      const mdet::Counter& dCounter = theMDetector->GetCounter(ic->GetId());
      const int sdStationId = dCounter.GetAssociatedTankId();

      if (!se.HasStation(sdStationId))
        continue;

      const sevt::Station& sdStation = se.GetStation(sdStationId);

      if (!sdStation.HasRecData())
        continue;

      const utl::TimeStamp sdSignalStartTime = sdStation.GetRecData().GetSignalStartTime();

      const utl::TimeInterval md2sdTimeInterval = mdT50 - sdSignalStartTime;

      const double distance = sdStation.GetRecData().GetSPDistance();  // distance to axis from the SD reconstruction

      // Reject counters with few muons and the t50 too late with respecto to the SD t0
      if (int(ic->GetNumberOfChannelsOn()) < fTimeOutlierMuons && IsTimeOutlier(md2sdTimeInterval, distance)) {
        ic->SetRejected("Time outlier");
        os.str("");
        os << "Counter " << dCounter.GetId() << " rejected because its t50 is out of time with the SD t0.";
        INFO(os);
      }

    }
  }


  bool
  MdEventSelector::IsTimeOutlier(const double md2sdTimeInterval, const double distance)
    const
  {
    // broken function used to cut time outliers: constant below and linear above fTimeOutlierDistance
    if (distance < fTimeOutlierDistance)
      return md2sdTimeInterval > fTimeOutlierConstant;
    else
      return md2sdTimeInterval > fTimeOutlierConstant + fTimeOutlierSlope * (distance - fTimeOutlierDistance);
  }


  bool
  MdEventSelector::IsMEventT5(const mevt::MEvent& me)
    const
  {
    //INFO("");

    const mevt::Counter* const hottestCounter = FindHottestCounter(me);

    const unsigned int hottestId = hottestCounter->GetId();
    utl::Point hottestPosition = theMDetector->GetCounter(hottestId).GetPosition();

    int functioningNeighbors = 0;

    for (mevt::MEvent::CounterConstIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {

      if (ic->IsCandidate() || ic->IsSilent()) {

        const unsigned int counterId = ic->GetId();
        const utl::Point counterPosition = theMDetector->GetCounter(counterId).GetPosition();

        const utl::Vector moduleToHottest = counterPosition-hottestPosition;
        const double distance = moduleToHottest.GetMag();

        //cout << "Counter " << ic->GetId() << ", distance = " << distance << endl;

        if (fFirstCrownMinDistance < distance && distance < fFirstCrownMaxDistance) {
          ++functioningNeighbors;
        }

      }

    }

    if (functioningNeighbors >= fT5NumberOfActiveStations) {
      INFO("Event is MD T5");
      return true;
    } else {
      ostringstream info;
      info << "Event is not MD T5: only " << functioningNeighbors << " "
              "(out of " << fT5NumberOfActiveStations << ") "
              "functioning neighbours of the hottest counter "
           << hottestCounter->GetId() << " found";
      INFO(info);

      return false;
    }

  }

}