NBinCenteredStrategy.cc

Go to the documentation of this file.

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

namespace MdMuonEstimatorAG {

  NBinCenteredStrategy::NBinCenteredStrategy(unsigned int nBinWindowSize, unsigned int nTraceLength) :
    fNBinWindowSize(nBinWindowSize),
    fTraceSize(nTraceLength)
  {
  }
  
  double
  NBinCenteredStrategy::operator()(std::vector<std::pair<unsigned int, unsigned int>>& vModulePatternMatches, mevt::ModuleRecData& mRecData)
    const
  {
    //std::cout << "In N bin centered 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;
         
      // Create a vector that contains the sum of pattern matches (any of the bins of the pattern contributes) in each time bin of the trace
      unsigned int vSummedPatternMatches[fTraceSize] = {0};
      for (std::vector<std::pair<unsigned int, unsigned int>>::const_iterator it = vModulePatternMatches.begin(); it != vModulePatternMatches.end(); ++it) {
        for (unsigned int i = 0 ; i < mRecData.GetWindowSize() ; ++i) { //mRecData.GetWindowSize() = size of the pattern match
          if ((int) it->second - (int) i >= 0) {
            ++vSummedPatternMatches[(int) it->second - (int) i];
          }
        }
      }
      //std::cout << "SUMMED VECTOR" << std::endl;
      //for (unsigned int i = 0; i < nTraceSamples; ++i) { 
      //    std::cout << vSummedPatternMatches[i] << " ";
      //}
      //std::cout << std::endl;
      
      //Find the (first) bin where the sum is maximum
      unsigned int imax = 0; // imax is the bin with maximum counts
      unsigned int maxcounts = 0; // maxcounts is the number of pattern matches (any of the bins of the pattern contributes) in the bin imax
      for (unsigned int i = 0; i < fTraceSize; ++i) { 
        if (vSummedPatternMatches[i] > maxcounts) {
          maxcounts = vSummedPatternMatches[i];
          imax = i;
        }
      }
      
      // Define the start time and size of the histogram such that there is a window of size fNBinWindowSize opening at time-bin imax
      unsigned int histogramStartTime = imax % fNBinWindowSize + 1; // beware shift between vector and TH1 indexes
      unsigned int histogramSize = fTraceSize - histogramStartTime + 1;
      //std::cout << " histogramStartTime " << histogramStartTime << " histogramSize " << histogramSize << endl;
      
      // Build a histogram of pattern matches
      TH1I hChannelsOn("h", "Pattern matches", histogramSize, histogramStartTime, histogramSize); //XXX Check!!!!
      for (std::vector<std::pair<unsigned int, unsigned int>>::const_iterator it = vModulePatternMatches.begin(); it != vModulePatternMatches.end(); ++it) {
        const unsigned int bin = it->second + 1; // beware shift between vector and TH1 indexes
        hChannelsOn.Fill(bin); // 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;
  } 
}