NBinStrategy.cc

Go to the documentation of this file.

#include "NBinStrategy.h"
#include <TH1I.h>

namespace MdMuonEstimatorAG {

  NBinStrategy::NBinStrategy(unsigned int nBinWindowSize, unsigned int nTraceLength) :
    fNBinWindowSize(nBinWindowSize),
    fTraceSize(nTraceLength)
  {
  }
  
  double
  NBinStrategy::operator()(std::vector<std::pair<unsigned int, unsigned int>>& vModulePatternMatches, mevt::ModuleRecData& mRecData)
    const
  {
    //std::cout << "In N bin strategy" << std::endl;
    //std::cout << "Size " << vModulePatternMatches.size() << std::endl;

    double mu = 0;
    double nmu = 0;
    bool isSat = 0;

    if (vModulePatternMatches.size()>0) {
      //std::cout << "Channels with signal ";
      //for (std::vector<std::pair<unsigned int, unsigned int>>::const_iterator it = vModulePatternMatches.begin(); it != vModulePatternMatches.end(); ++it) {
        //std::cout << it->first << " ";
      //}
      //std::cout << std::endl;

      // Build a histogram of pattern matches
      TH1I hChannelsOn("h1", "Pattern matches", fTraceSize, 1, fTraceSize);
      for (std::vector<std::pair<unsigned int, unsigned int>>::const_iterator it = vModulePatternMatches.begin(); it != vModulePatternMatches.end(); ++it) {
        const unsigned int tbin = it->second + 1; // beware shift between vector and TH1 indexes
        hChannelsOn.AddBinContent(tbin,1);  // load the histogram
      }
      // Rebin the histogram in windows of size fNBinWindowSize
      if (fNBinWindowSize <= fTraceSize) {
        hChannelsOn.Rebin(fNBinWindowSize);  // group bins according to the window size
      } else {
        hChannelsOn.Rebin(fTraceSize);
      }

      // Count
      unsigned int nActiveCh = mRecData.GetSegmentation(); // = number of active channels
      for (int i = 0; i <= hChannelsOn.GetNbinsX(); ++i) {
        unsigned int ki = hChannelsOn.GetBinContent(i);
        if (ki > 0) {
          mu -= (double) nActiveCh * log(1. - (double) ki/(double) nActiveCh);
          nmu += log(1. - (double) ki/(double) nActiveCh) / log(1.-1./(double) nActiveCh);
        }
        if (ki >= nActiveCh){
          isSat = 1;
          //std::cout << "Module is saturated" << std::endl;
        }
      }
    }
    //std::cout << "Mu " << mu << std::endl;
    //std::cout << "Nmu " << nmu << std::endl;
    //std::cout << "IsSat " << isSat << std::endl;
    
    mRecData.SetMeanMuons(mu);
    mRecData.SetNumberOfEstimatedMuons(nmu);
    mRecData.SetSaturated(isSat);
    
    return nmu;
  }  
}