MdMuonEstimator.cc

Go to the documentation of this file.

// $Id: MdMuonEstimator.cc 34060 2021-03-11 16:43:53Z sanchez $
#include "MdMuonEstimator.h"

#include "OneBinStrategy.h"
#include "NBinStrategy.h"
#include "NBinCenteredStrategy.h"
#include "InfiniteWindowStrategy.h"

#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <mdet/MDetector.h>
#include <sdet/SDetector.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <sevt/SEvent.h>
#include <mevt/Module.h>

#include <utl/ConsecutiveEnumFactory.h>
#include <utl/Trace.h>
#include <utl/AugerUnits.h>
#include <utl/Branch.h>
#include <utl/RandomEngine.h>
#include <utl/ConfigUtils.h>
#include <utl/TimeStamp.h>
#include <utl/PhysicalConstants.h>
#include <utl/TimeInterval.h>

#include <utl/ErrorLogger.h>
#include <utl/Verbosity.h>

#include <sstream>
//#include <cassert>

#include <TH1I.h>

using namespace fwk;
using std::cout;
using std::endl;


namespace MdMuonEstimatorAG {

  const char* const MdMuonEstimator::kCountingStrategyTags[] = { "OneBin", "NBin", "NBinCentered", "InfiniteWindow" };


  VModule::ResultFlag
  MdMuonEstimator::Init() //evt::Event& event
  {
    DEBUGLOG("Initiating Muon Estimator");

    // Configuration reading.
    utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdMuonEstimator");
    fUnits.SetTimeDefault(utl::ns, "ns");
    fUnits.Configure(config);
    fLog.Configure(config);
    std::string tag;
    LoadConfig(config, "countingStrategy", tag, kCountingStrategyTags[0]);
    LoadConfig(config, "nBinWindowSize", fNBinWindowSize, 12);

    fCountingStrategy = utl::ConsecutiveEnumFactory<CountingStrategyType, eInfiniteWindow, kCountingStrategyTags>::Create(tag);

    //cout << kCountingStrategyTags[ fCountingStrategy ] << endl;
    //cout << "Muon Estimator Initiated" << endl;
    return eSuccess;
  }


  VModule::ResultFlag
  MdMuonEstimator::Run(evt::Event& event)
  {
    if (!event.HasMEvent()) {
      WARNING("\n+++++++++\nEvent without MEvent: nothing to be done. Skipping to next event.\n++++++++++\n");
      return eContinueLoop;
    }

    // Initiate counting strategy
    unsigned int traceLength = GetTraceLength(event);
    double samplingTimeTrace = GetSamplingTime(event);

    switch (fCountingStrategy) {
    case eOneBin:
      fEstimator.reset(new OneBinStrategy(traceLength, samplingTimeTrace));
      fLog("Strategy")(kCountingStrategyTags[ fCountingStrategy ])(".");
      fCountingWindowSize = 1;
      break;
    case eNBin:
      fEstimator.reset(new NBinStrategy(fNBinWindowSize, traceLength));
      fLog("#NBin Window size", 1)
        (fNBinWindowSize)("-")
        ("Strategy")(kCountingStrategyTags[ fCountingStrategy ])(".");
      fCountingWindowSize = fNBinWindowSize;
      break;
    case eNBinCentered:
      fEstimator.reset(new NBinCenteredStrategy(fNBinWindowSize, traceLength));
      fLog("#NBin Window size", 1)
        (fNBinWindowSize)("-")
        ("Strategy")(kCountingStrategyTags[ fCountingStrategy ])(".");
      fCountingWindowSize = fNBinWindowSize;
      break;
    case eInfiniteWindow:
      fEstimator.reset(new InfiniteWindowStrategy(traceLength));
      fLog("Strategy")(kCountingStrategyTags[ fCountingStrategy ])(".");
      fCountingWindowSize = traceLength; // Set to length of the trace
      break;
    default:
      ERROR("Invalid muon estimation strategy!");
      break;
    } //end switch

    mevt::MEvent& me = event.GetMEvent();
    const mdet::MDetector& theMdDetector = det::Detector::GetInstance().GetMDetector();

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

      mevt::Counter& counter = *ic;
      const int counterId = counter.GetId();

      if (counter.IsRejected()) {
         os.str("");
         os << "Skipping the rejected counter " << counterId;
         INFO(os);
         continue;
      }

      os.str("");
      os << "Processing counter " << counterId;
      INFO(os);

      const mdet::Counter& mdetCounter = theMdDetector.GetCounter(counterId);
      //int stationId = mdetCounter.GetAssociatedTankId();

      for (mevt::Counter::ModuleIterator im = counter.ModulesBegin(), em = counter.ModulesEnd(); im != em; ++im) {

        DEBUGLOG("Starting module loop");
        mevt::Module& module = *im;
        const int moduleId = module.GetId();
        //const int moduleGlobalId = counterId*100 + moduleId;

        os.str("");
        os << "Processing module " << moduleId;
        DEBUGLOG(os);

        const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

        // Check if module was flagged as rejected (by a selector) or, in the case of real data, if T1 from SD was not matched
        // this should imply empty module traces (all "0s")
        if (module.IsRejected()) {
           os.str("");
           os << "Module " << module.GetId() << " of counter " << counterId;
           os << " discarded because flagged as: " << module.GetRejectionReason();
           WARNING(os);
           continue;
        }

        if (!module.HasRecData()) {
          module.MakeRecData();
        }
        DEBUGLOG("After making mRecData");

        // Get a vector with the channel ids and bin numbers of the pattern matches
        const std::vector<std::pair<unsigned int, unsigned int>> vModulePatternMatches = GetModulePatternMatches(module);
        DEBUGLOG("After getting vModulePatternMatches");

        // Fill the channelsOn histogram
        FillChannelsOn(module, vModulePatternMatches, mdetModule);
        DEBUGLOG("After filling channelsOn");

        // The pattern matches of the channels are already filled into
        // the ChannelRecData by MdPatternFinder
        // Fill that information into the ModuleRecData
        LoadPatternMatchesInModuleRecData(module, vModulePatternMatches, mdetModule);
        DEBUGLOG("After loading pattern matches in the module reconstructed data");

        SetSegmentationAndArea(module, mdetModule); //Before EstimateNumberOfMuons!
        DEBUGLOG("After setting the segmentation and area");

        EstimateNumberOfMuons(module, vModulePatternMatches);
        DEBUGLOG("After estimating number of muons");
      } // end module loop
    } // end counter loop

    return eSuccess;
  }


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


  unsigned int
  MdMuonEstimator::GetTraceLength(evt::Event& event)
    const
  {
    unsigned int nTraceSamples = 0;

    mevt::MEvent& me = event.GetMEvent();
    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {
      mevt::Counter& counter = *ic;
      if (counter.IsRejected()) {
        continue;
      }
      const int counterId = counter.GetId();
      for (mevt::Counter::ModuleIterator im = counter.ModulesBegin(), em = counter.ModulesEnd(); im != em; ++im) {
        mevt::Module& module = *im;
        if (module.IsRejected()) {
          continue;
        }
        const int moduleId = module.GetId();

        const mdet::MDetector& theMdDetector = det::Detector::GetInstance().GetMDetector();
        const mdet::Counter& mdetCounter = theMdDetector.GetCounter(counterId);
        const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

        nTraceSamples = mdetModule.IsSiPM() ?
        mdetModule.GetFrontEndSiPM().GetBufferLength() :
        mdetModule.GetFrontEnd().GetBufferLength();

        //cout << "nTraceSamples " << nTraceSamples << endl;
        if (!nTraceSamples) {
          break;
        }
      }
      if (!nTraceSamples) {
        break;
      }
    }

    return nTraceSamples;
  }


double
  MdMuonEstimator::GetSamplingTime(evt::Event& event)
    const
  {
    double samplingTime = 0;

    mevt::MEvent& me = event.GetMEvent();
    for (mevt::MEvent::CounterIterator ic = me.CountersBegin(); ic != me.CountersEnd(); ++ic) {
      mevt::Counter& counter = *ic;
      if (counter.IsRejected()) {
        continue;
      }
      const int counterId = counter.GetId();
      for (mevt::Counter::ModuleIterator im = counter.ModulesBegin(), em = counter.ModulesEnd(); im != em; ++im) {
        mevt::Module& module = *im;
        if (module.IsRejected()) {
          continue;
        }
        const int moduleId = module.GetId();

        const mdet::MDetector& theMdDetector = det::Detector::GetInstance().GetMDetector();
        const mdet::Counter& mdetCounter = theMdDetector.GetCounter(counterId);
        const mdet::Module& mdetModule = mdetCounter.GetModule(moduleId);

        samplingTime = mdetModule.IsSiPM() ?
          mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
          mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();

        //cout << "samplingTime " << samplingTime << endl;
        if (!samplingTime) {
          break;
        }
      }
      if (!samplingTime) {
        break;
      }
    }

    return samplingTime;
  }


  std::vector<std::pair<unsigned int, unsigned int>>
  MdMuonEstimator::GetModulePatternMatches(mevt::Module& module)
    const
  {
    // Get pattern match bins
    std::vector<std::pair<unsigned int, unsigned int>> vMatchBins;
    for (mevt::Module::ChannelIterator channel = module.ChannelsBegin(), ech = module.ChannelsEnd();
         channel != ech; ++channel) {

      if (channel->HasRecData() && !channel->IsMasked()) {  // Use unmasked channels only
        for (mevt::ChannelRecData::PatternMatchBinIterator mut = channel->GetRecData().PatternMatchBinsBegin(),
             end = channel->GetRecData().PatternMatchBinsEnd(); mut != end; ++mut) {
          vMatchBins.push_back(std::make_pair(channel->GetId(),*mut));
        }
      }
    }

    return vMatchBins;
  }


  void
  MdMuonEstimator::FillChannelsOn(mevt::Module& module, const std::vector<std::pair<unsigned int, unsigned int>>& vModulePatternMatches, const mdet::Module& mdetModule)
    const
  {
    mevt::ModuleRecData& mRecData = module.GetRecData();

    // ChannelsOnTrace: Module-level trace that adds a 1 in every time-bin where a pattern match starts
    if ( !mRecData.HasChannelsOn())
        mRecData.MakeChannelsOn();

    utl::TraceUI& channelsOnTrace = mRecData.GetChannelsOn();

    // ChannelsInhibitedTrace: Module-level trace that adds a 1 in every time-bin of a pattern match except the first bin
    if ( !mRecData.HasChannelsInhibited() )
        mRecData.MakeChannelsInhibited();

    utl::TraceUI& channelsInhibitedTrace = mRecData.GetChannelsInhibited();

    // Get the size and binning of the module (they are the same for all channels)
    const unsigned int traceSize = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetBufferLength() :
      mdetModule.GetFrontEnd().GetBufferLength();

    // Fill auxiliary vectors to fill later the channelsOn and channelsInhibited traces
    std::vector<unsigned int> channelsOnVector(traceSize,0);
    std::vector<unsigned int> channelsInhibitedVector(traceSize,0);
    for (std::vector<std::pair<unsigned int, unsigned int>>::const_iterator it = vModulePatternMatches.begin(); it != vModulePatternMatches.end(); ++it) {
      const unsigned int bin = it->second;
      ++channelsOnVector[bin];
      for (unsigned int i = 1; i < (unsigned int) mRecData.GetWindowSize(); ++i){
        ++channelsInhibitedVector[bin+i];
      }
    }
    // Fill the ChannelsOnTrace
    for (std::vector<unsigned int>::const_iterator it = channelsOnVector.begin(); it != channelsOnVector.end(); ++it){
      channelsOnTrace.PushBack(*it);
    }
    // Fill the ChannelsInhibitedTrace
    for (std::vector<unsigned int>::const_iterator it = channelsInhibitedVector.begin(); it != channelsInhibitedVector.end(); ++it){
      channelsInhibitedTrace.PushBack(*it);
    }

    // Set the ChannelsOn binning
    const double samplingTime = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
      mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();
    const double channelsOnBinning = samplingTime * fCountingWindowSize;
    channelsOnTrace.SetBinning(channelsOnBinning);

    std::ostringstream info;
    info << "Counting channels on\n"
            "Number of bins = " << traceSize << "\n"
            "Binning (ns) = " << channelsOnBinning;
    DEBUGLOG(info);
  }

  void
  MdMuonEstimator::LoadPatternMatchesInModuleRecData(mevt::Module& module, const std::vector<std::pair<unsigned int, unsigned int>>& vModulePatternMatches, const mdet::Module& mdetModule)
    const
  {
    DEBUGLOG("Starting LoadPatternMatchesInModuleRecData");
    mevt::ModuleRecData& mRecData = module.GetRecData();

    // Fill the pattern match times vector
    const double samplingTime = mdetModule.IsSiPM() ?
      mdetModule.GetFrontEndSiPM().GetMeanSampleRatePeriod() :
      mdetModule.GetFrontEnd().GetMeanSampleRatePeriod();

    for (std::vector<std::pair<unsigned int, unsigned int>>::const_iterator it = vModulePatternMatches.begin();
         it != vModulePatternMatches.end(); ++it) {
      const double matchTime = it->second * samplingTime;
      mRecData.AddPatternMatchTime(matchTime);
      //std::cout << "Recontructed matched times in counter " << cId << ", Module " << module.GetId() << " " << matchTime << std::endl;
    }
    DEBUGLOG("After filling the pattern match times vector");
  }

void
  MdMuonEstimator::SetSegmentationAndArea(mevt::Module& module, const mdet::Module& mdetModule)
    const
  {
    mevt::ModuleRecData& mRecData = module.GetRecData();

    const size_t segmentation = module.GetNumberOfActiveChannels();
    mRecData.SetSegmentation(segmentation);
    DEBUGLOG("After setting the segmentation");

    const mdet::Scintillator& scintillator = *(mdetModule.ScintillatorsBegin());
    const double scintillatorArea = scintillator.GetArea();
    const double moduleArea = scintillatorArea * segmentation;
    mRecData.SetActiveArea(moduleArea);
    DEBUGLOG("After setting the active area");
  }

void
  MdMuonEstimator::EstimateNumberOfMuons(mevt::Module& module, std::vector<std::pair<unsigned int, unsigned int>> vModulePatternMatches)
    const
  {
    // TODO: Add error estimation
    mevt::ModuleRecData& mRecData = module.GetRecData();
    std::ostringstream os;

    // Set the number of estimated muons
    const double numberOfMuons = (*fEstimator)(vModulePatternMatches, mRecData);
    mRecData.SetNumberOfEstimatedMuons(numberOfMuons);
    DEBUGLOG("After setting the number of muons");

    // Calculate errors with the likelihood at a later step in the reconstruction, using the profile likelihood
    mRecData.SetNumberOfMuonsErrorHigh(0);
    mRecData.SetNumberOfMuonsErrorLow(0);

    if (!mRecData.IsSaturated()) {
      os.str("");
      os << "module " << module.GetId() << ", "
          << std::fixed << std::setprecision(1)
          << "area = " << mRecData.GetActiveArea() << " m2, "
          << "bars = " << mRecData.GetSegmentation() << ", "
          << "muons = " << numberOfMuons; // << " +" << errorHigh << " -" << errorLow;
      INFO(os);
    } else {
      os.str("");
      os << "module " << module.GetId() << " is saturated";
      INFO(os);
    }
    // Set the lower limit to the number of muons for both saturated and unsaturated modules
    mRecData.SetNumberOfMuonsLowLimit(0);

    DEBUGLOG("Finishing EstimateNumberOfMuons");
  }

}