SVDTimeGroupingModule.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
// framework
#include <framework/core/Module.h>
#include <framework/database/DBObjPtr.h>
#include <framework/datastore/StoreArray.h>
 
// svd
#include <vxd/dataobjects/VxdID.h>
#include <svd/dbobjects/SVDRecoConfiguration.h>
#include <svd/dbobjects/SVDTimeGroupingConfiguration.h>
#include <svd/dataobjects/SVDCluster.h>
 
// std
#include <string>
 
// root
#include <TH1D.h>
#include <TF1.h>
#include <TMath.h>
 
 
namespace Belle2 {
  typedef std::tuple<double, double, double> GroupInfo;
 
 
 
  class SVDTimeGroupingModule : public Module {
 
  public:
 
 
 
    SVDTimeGroupingModule();
 
 
    virtual void initialize() override;
 
    void beginRun() override;
 
    virtual void event() override;
 
 
  protected:
 
    DBObjPtr<SVDRecoConfiguration> m_recoConfig; 
    DBObjPtr<SVDTimeGroupingConfiguration> m_groupingConfig; 
    // Data members
    std::string m_svdClustersName; 
    std::string m_svdEventInfoName; 
    StoreArray<SVDCluster> m_svdClusters;
 
    bool m_forceGroupingFromDB; 
    bool m_isEnabledIn6Samples; 
    bool m_isEnabledIn3Samples; 
    bool m_useParamFromDB; 
    SVDTimeGroupingParameters m_usedParsIn6Samples;
 
    SVDTimeGroupingParameters m_usedParsIn3Samples;
 
    SVDTimeGroupingParameters m_usedPars;
 
    // helper functions
 
    void createAndFillHistorgram(TH1D& hist);
 
    void searchGausPeaksInHistogram(TH1D& hist, std::vector<GroupInfo>& groupInfoVector);
 
    void resizeToMaxSize(std::vector<GroupInfo>& groupInfoVector)
    {
      groupInfoVector.resize(m_usedPars.maxGroups, GroupInfo(0., 0., 0.));
    }
 
    void sortBackgroundGroups(std::vector<GroupInfo>& groupInfoVector);
 
    void sortSignalGroups(std::vector<GroupInfo>& groupInfoVector);
 
    void assignGroupIdsToClusters(TH1D& hist, std::vector<GroupInfo>& groupInfoVector);
 
  };
 
 
 
  inline double myGaus(const double* x, const double* par)
  {
    return par[0] * TMath::Gaus(x[0], par[1], par[2], true);
  }
 
 
  inline void addGausToHistogram(TH1D& hist,
                                 const double& integral, const double& center, const double& sigma,
                                 const double& sigmaN, const bool& isAddition = true)
  {
    int startBin = hist.FindBin(center - sigmaN * sigma);
    int   endBin = hist.FindBin(center + sigmaN * sigma);
    if (startBin < 1) startBin = 1;
    if (endBin > (hist.GetNbinsX())) endBin = hist.GetNbinsX();
 
    for (int ijx = startBin; ijx <= endBin; ijx++) {
      double tbinc       = hist.GetBinCenter(ijx);
      double tbincontent = hist.GetBinContent(ijx);
 
      if (isAddition) tbincontent += integral * TMath::Gaus(tbinc, center, sigma, true);
      else tbincontent -= integral * TMath::Gaus(tbinc, center, sigma, true);
 
      hist.SetBinContent(ijx, tbincontent);
    }
  }
 
 
  inline void subtractGausFromHistogram(TH1D& hist,
                                        const double& integral, const double& center, const double& sigma,
                                        const double& sigmaN)
  {
    addGausToHistogram(hist, integral, center, sigma, sigmaN, false);
  }
 
  inline int getSensorType(const VxdID& sensorID)
  {
    int layer  = sensorID.getLayerNumber();
    int sensor = sensorID.getSensorNumber();
    if (layer == 3)
      return 0;
    else {
      if (sensor == 1)
        return 1;
      else
        return 2;
    }
  }
 
 
} // end namespace Belle2