OneBinStrategy.cc

Go to the documentation of this file.

#include "OneBinStrategy.h"
#include <TH1I.h>
#include <utl/Trace.h>

namespace MdMuonEstimatorAG {

  OneBinStrategy::OneBinStrategy(unsigned int nTraceLength, double sampTime) :
    fTraceSize(nTraceLength),
    fSamplingTime(sampTime)
  {
  }

  double
  OneBinStrategy::operator()(std::vector<std::pair<unsigned int, unsigned int>>& vModulePatternMatches, mevt::ModuleRecData& mRecData)
    const
  {
    //std::cout << "In one bin strategy" << std::endl;
    //std::cout << "Size " << vModulePatternMatches.size() << std::endl;

    double mu = 0;
    double nmu = 0;
    bool isSat = 0;
    std::vector<double> vMu(fTraceSize, 0); //vector of length fTraceSize initialized to 0
    std::vector<double> vNmu(fTraceSize, 0); //vector of length fTraceSize initialized to 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
      }

      //std::cout << "PatternMatch histogram" << std::endl;
      //for (int i = 0; i <= hChannelsOn.GetNbinsX(); ++i) {
        //if (hChannelsOn.GetBinContent(i)>0)
          //std::cout << i << " ";
      //}
      //std::cout << std::endl;
      //for (int i = 0; i <= hChannelsOn.GetNbinsX(); ++i) {
        //if (hChannelsOn.GetBinContent(i)>0)
          //std::cout << hChannelsOn.GetBinContent(i) << " ";
      //}
      //std::cout << std::endl;


      // 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) {
          // Count inhibited channels for this time bin
          // These are the number of pattern matches beginning at most mRecData.GetWindowSize()-1 (size of the pattern match) bins before i
          unsigned int numInhibitedCh = 0;
          for (int deltai = 1; ((deltai <= i) & (deltai < (int) mRecData.GetWindowSize())); ++deltai) {
            numInhibitedCh += hChannelsOn.GetBinContent(i-deltai);
            //std::cout << deltai << " ";
          }
          //std::cout << " , inhibChannels: " << numInhibitedCh << endl;
          if (numInhibitedCh < nActiveCh) { // if not all channels are inhibited = if the module is not saturated
            double mui = - (double) nActiveCh * log(1 - (double) ki/((double) nActiveCh- (double) numInhibitedCh));
            double nmui = ((double) nActiveCh/((double) nActiveCh- (double) numInhibitedCh))* log(1 - (double) ki/((double) nActiveCh- (double) numInhibitedCh))/log(1.-1./((double) nActiveCh- (double) numInhibitedCh));
            mu += mui;
            nmu += nmui;
            // Fill the reconstructed muon/number of muons trace
            vMu[i] = mui;
            vNmu[i] = nmui;
          }
          if (ki + numInhibitedCh >= nActiveCh) {
            isSat = 1;
            //std::cout << "Module is saturated" << std::endl;
          }
        }
      }
    }

    // Fill the reconstructed mu and nmu traces

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

    if (!mRecData.HasNumberOfMuonsVsTime())
      mRecData.MakeNumberOfMuonsVsTime();

    utl::Trace<double>& muTrace = mRecData.GetMeanMuonsVsTime();
    utl::Trace<double>& nmuTrace = mRecData.GetNumberOfMuonsVsTime();
    muTrace.Clear();
    nmuTrace.Clear();

    for (unsigned int i = 0; i < fTraceSize; ++i) {
      const double vMui = vMu[i];
      muTrace.PushBack(vMui);
      const double vNmui = vNmu[i];
      nmuTrace.PushBack(vNmui);
    }
    // Set the binning
    muTrace.SetBinning(fSamplingTime);
    nmuTrace.SetBinning(fSamplingTime);

    // Print out the traces
    //for (utl::TraceD::ConstIterator it = muTrace.Begin(); it != muTrace.End(); ++it) {
      //if (*it != 0)
        //std::cout << *it << " ";
    //}
    //std::cout << std::endl;

    //for (utl::TraceD::ConstIterator it = nmuTrace.Begin(); it != nmuTrace.End(); ++it) {
      //if (*it != 0)
        //std::cout << *it << " ";
    //}
    //std::cout << std::endl;


    // Fill the reconstructed mu, nmu, and isSat
    //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;
  }
}