InhibitionStrategy.cc

Go to the documentation of this file.

#include "InhibitionStrategy.h"
//
#include <mevt/Module.h>
#include <mdet/Module.h>
//
#include <utl/Trace.h>
#include <utl/TabularStream.h>

namespace MdPileUpCorrecterAG {

  InhibitionStrategy::InhibitionStrategy()
  {
  }

  boost::tuple<double, double>
  InhibitionStrategy::operator()(mevt::Module& module, const mdet::Module& mdetModule )
    const
  {
    double nSingleRealizationMuons = 0;
    double nAverageMuons = 0;

    const unsigned int numberOfSamples = mdetModule.IsSiPM() ?
    mdetModule.GetFrontEndSiPM().GetBufferLength() :
    mdetModule.GetFrontEnd().GetBufferLength();

    mevt::ModuleRecData& mRecData = module.GetRecData();
    unsigned int ns = mRecData.GetSegmentation();

    //Create the Channels On trace
    if (!mRecData.HasChannelsOn())
      mRecData.MakeChannelsOn();

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

    //Create the Channels Inhibited trace
     if (!mRecData.HasChannelsInhibited())
      mRecData.MakeChannelsInhibited();

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

    for ( size_t bin = 0; bin < numberOfSamples; ++bin ) {

      size_t chanOn = 0;
      size_t chanIn = 0;
      for (mevt::Module::ChannelIterator ch = module.ChannelsBegin(), ech = module.ChannelsEnd(); ch != ech; ++ch) {

        mevt::Channel& channel = *ch;

        if (channel.HasRecData() && !channel.IsMasked() && channel.GetRecData().HasTrace() ) {//use unmasked channels only with reconstructed trace

          const utl::TraceUI traceRec = channel.GetRecData().GetTrace();

          if ( traceRec.GetSize() ){

          if ( traceRec.At(bin) == 1 )
            chanOn++;
          else
          if ( traceRec.At(bin) == 2 )
            chanIn++;

         }
        }
      }//end loop on channels at fixed bin

      //Set values of channels On and Inhibited per bin
      channelsOnTrace.PushBack(chanOn);
      channelsInhiTrace.PushBack(chanIn);
      // If saturated (this is, even a single bin in the channels ON trace
      // with number of channels on equal to
      // segmentation) not possible to estimate the number of muons
      if ( chanOn + chanIn == ns )
        mRecData.SetSaturated(true);
     }//end loop on bins

    if ( mRecData.IsSaturated() )
      return boost::tuple<double, double>(-1, -1);

    //Create the muons Vs time trace
    if (!mRecData.HasNumberOfMuonsVsTime())
      mRecData.MakeNumberOfMuonsVsTime();

    utl::Trace<double>& muSTrace = mRecData.GetNumberOfMuonsVsTime();
    muSTrace.Clear();

    if (!mRecData.HasMeanMuonsVsTime())
      mRecData.MakeMeanMuonsVsTime();

    utl::Trace<double>& muATrace = mRecData.GetMeanMuonsVsTime();
    muATrace.Clear();

    //Apply pile-up correction in bin-by-bin basis
    for ( size_t bin = 0; bin < numberOfSamples; ++bin ) {

        unsigned int k  = channelsOnTrace.At(bin);

        if ( !k ) continue;//Nothing to be done

        unsigned int ni = channelsInhiTrace.At(bin);

        //NOTE: ni==ns if k==0
        unsigned int bin_segmentation = ns-ni;

        double S = double(ns)/bin_segmentation * std::log(1 - double(k)/bin_segmentation)/std::log(1 - 1.0/bin_segmentation);
        double A = ns * std::log(1 - double(k)/bin_segmentation);

        nSingleRealizationMuons += S;
        nAverageMuons           -= A;

        muSTrace.PushBack(S);
        muATrace.PushBack(A);

    }
    return boost::tuple<double, double>(nAverageMuons, nSingleRealizationMuons);
  }


}