MdPileUpCorrecter.cc

Go to the documentation of this file.

#include "MdPileUpCorrecter.h"

#include "InhibitionStrategy.h"

#include <fwk/CentralConfig.h>
#include <det/VManager.h>
#include <mdet/MDetector.h>
#include <mevt/Counter.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/TimeInterval.h>
#include <utl/TabularStream.h>

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

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

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


namespace MdPileUpCorrecterAG {

  const char* const MdPileUpCorrecter::kStrategyTags[] = { "Inhibit" };

  VModule::ResultFlag
  MdPileUpCorrecter::Init()
  {
    // Configuration reading.
    utl::Branch config = fwk::CentralConfig::GetInstance()->GetTopBranch("MdPileUpCorrecter");
    fUnits.SetTimeDefault(utl::ns, "ns");
    fUnits.Configure(config);
    fLog.Configure(config);
    std::string tag;
    LoadConfig(config, "strategy", tag, kStrategyTags[0]);
    fStrategy = utl::ConsecutiveEnumFactory<Strategy, eInhibit, kStrategyTags>::Create(tag);

    // Init according to selected strategy.
    switch (fStrategy) {
    case eInhibit:
      fCorrecter.reset(new InhibitionStrategy());
      fLog("Strategy",1)(kStrategyTags[ fStrategy ])(".");
      break;
    }

    return eSuccess;
  }


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

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

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

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

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

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

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

        mevt::Module& module = *im;
        const int moduleId = module.GetId();
        const mdet::Module&  mdetModule = mdetCounter.GetModule(moduleId);

        os.str("");
        os << "Processing module " << 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 an empty (all "0s" ) module traces
        if (module.IsRejected()) {
           os.str("");
           os << "Module " << module.GetId() << " of counter " << counterId;
           os << " discarded because flagged as: " << module.GetRejectionReason();
           WARNING(os);
           continue;
        }

        double nSingleRealizationMuons = 0;
        double nAverageMuons = 0;
        // Reconstructed muons
        boost::tie(nAverageMuons, nSingleRealizationMuons) = (*fCorrecter)(module,mdetModule);

        mevt::ModuleRecData& mRecData = module.GetRecData();

        mRecData.SetNumberOfEstimatedMuons(nSingleRealizationMuons);
        mRecData.SetMeanMuons(nAverageMuons);
        // Calculate errors with the likelihood at a later step in the reconstruction, using the profile likelihood
        mRecData.SetNumberOfMuonsErrorHigh(0);
        mRecData.SetNumberOfMuonsErrorLow(0);
        // Set the lower limit to the number of muons for both saturated and unsaturated modules
        mRecData.SetNumberOfMuonsLowLimit(0);
        //Set ActiveArea
        const mdet::Scintillator& s = *(mdetModule.ScintillatorsBegin());
        const double moduleArea = s.GetArea() * im->GetRecData().GetSegmentation();
        mRecData.SetActiveArea( moduleArea );


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

        const unsigned int numberOfSamples = mdetModule.IsSiPM() ?
        mdetModule.GetFrontEndSiPM().GetBufferLength() :
        mdetModule.GetFrontEnd().GetBufferLength();
 
        std::unique_ptr<utl::TabularStream> spanTable;
        const unsigned int LogLevel = 15;
        if (fLog.GetLevel() >= LogLevel) {
          InitTable(spanTable, 5);
          *spanTable << utl::hline
                     << "Ctr"  << utl::endc //1
                     << "Mod"  << utl::endc //2
                     << "#Bin" << utl::endc //3
                     << "On"   << utl::endc //4
                     << "Inhi" << utl::endr;//5

          for ( size_t bin = 0; bin < numberOfSamples; ++bin ) {
             if (channelsOnTrace.At(bin)||channelsInhiTrace.At(bin)) 
             *spanTable                                      << utl::hline
                        << counterId                         << utl::endc
                        << moduleId                          << utl::endc
                        << bin                               << utl::endc
                        << fLog << channelsOnTrace.At(bin)   << utl::endc
                        << fLog << channelsInhiTrace.At(bin) << utl::endr 
                                                             << utl::hline;
          }
          fLog("Channels On & Inhibited span:");
          fLog(spanTable->Str());
        }

        fLog("Estimated muons in")(mRecData.IsSaturated()?"SATURATED":"")("Cter")(counterId)
            ("Mod")(moduleId)(": Single realization")(nSingleRealizationMuons)("- Average")(nAverageMuons);

      } // end module loop

    } // end counter loop

    return eSuccess;
  }


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

  void
  MdPileUpCorrecter::InitTable(std::unique_ptr<utl::TabularStream>& pt, unsigned int nCol)
    const
  {
    std::stringstream fmt;
    for(unsigned int i = 0; i < nCol; ++i)
      fmt << "|r";
    fmt << "|";
    pt.reset(new utl::TabularStream(fmt.str()));
    fLog.ApplyConfigurationOn(*pt);
  }

}