src/BordersFinder2D.h

-/** @file   BordersFinder2D.h
-    @class  Centrality::BordersFinder2D
-    @author Viktor Klochkov (klochkov44@gmail.com)
-    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
-    @brief  Class to determine centrality for 2D case. Inherited from Centrality::BordersFinder.
-*/
-
-#ifndef CENTRALITY_BORDERSFINDER2D_H
-#define CENTRALITY_BORDERSFINDER2D_H
-
-#include "BordersFinder.h"
-
-#include "vector"
-#include "array"
-
-#include "TH2.h"
-#include "TF1.h"
-
-
-namespace Centrality {
-    
-class BordersFinder2D : public BordersFinder {
-
-public:
-    
-    BordersFinder2D(){}
-        
-    void SetHisto2D(TH2F&& histo2d) { histo2d_ = histo2d; }
-    TH2F&& GetHisto2D() { return std::move(histo2d_); }
-    
-    void Init();
-    std::unique_ptr<TH1F> Convert();
-    void Fit2D(const TString func);
-    std::array <float,2> FindNorm (const std::vector <double> par, float x);
-    float FindIntegral( const std::array <float,2> &norm1, const std::array <float,2> &norm2);
-    void SaveBorders2D(const std::string &filename, const std::string &getter_name);
-    
-    /**
-     *
-     * @param par parameters of polinom
-     * @param x argument
-     * @param N order
-     * @return
-     */
-    float polN (std::vector <double> par, float x)
-    {
-        float res{0.};
-        float xn{1.};
-        for (const auto ipar : par){
-            res += ipar*xn;
-            xn *= x;
-        }
-        return res;
-    }
-
-private:
-    
-    TH2F histo2d_;
-    TF1* fit_{nullptr};
-    
-    TString fitname_{""};
-    
-    float xmax_{1.};
-    float ymax_{1.};
-    
-//     ClassDef(BordersFinder2D, 1);
-    
-};
-}
-
-#endif //CENTRALITY_BORDERSFINDER2D_H

src/BordersFinder2D.hpp

+/** @file   BordersFinder2D.h
+    @class  Centrality::BordersFinder2D
+    @author Viktor Klochkov (klochkov44@gmail.com)
+    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
+    @brief  Class to determine centrality for 2D case. Inherited from Centrality::BordersFinder.
+*/
+
+#ifndef CENTRALITY_BORDERSFINDER2D_H
+#define CENTRALITY_BORDERSFINDER2D_H
+
+#include "BordersFinder.hpp"
+
+#include "array"
+#include "vector"
+
+#include "TF1.h"
+#include "TH2.h"
+
+namespace Centrality{
+
+class BordersFinder2D : public BordersFinder{
+
+ public:
+  BordersFinder2D() = default;
+
+  void SetHisto2D(TH2F&& histo2d) { histo2d_ = histo2d; }
+  TH2F&& GetHisto2D() { return std::move(histo2d_); }
+
+  void Init();
+  std::unique_ptr<TH1F> Convert();
+  void Fit2D(const TString& func);
+  std::array<double, 2> FindNorm(const std::vector<double>& par, double x);
+  double FindIntegral(const std::array<double, 2>& norm1, const std::array<double, 2>& norm2);
+  void SaveBorders2D(const std::string& filename, const std::string& getter_name);
+
+  /**
+   *
+   * @param par parameters of polinom
+   * @param x argument
+   * @param N order
+   * @return
+   */
+  static double polN(const std::vector<double>& par, double x) {
+    double res{0.};
+    double xn{1.};
+    for(const auto ipar : par) {
+      res += ipar * xn;
+      xn *= x;
+    }
+    return res;
+  }
+
+ private:
+  TH2F histo2d_;
+  TF1* fit_{nullptr};
+
+  TString fitname_{""};
+
+  double xmax_{1.};
+  double ymax_{1.};
+
+  //     ClassDef(BordersFinder2D, 1);
+};
+}// namespace Centrality
+
+#endif//CENTRALITY_BORDERSFINDER2D_H

src/BordersFinderHelper.cpp

-#include "BordersFinderHelper.h"
+#include "BordersFinderHelper.hpp"
 
 #include <iostream>
 
-#include "TRandom.h"
 #include "TCanvas.h"
 #include "TLine.h"
+#include "TRandom.h"
+
+#include "TStyle.h"
 
 namespace Centrality {
 
-void BordersFinderHelper::QA(const Getter& getter, const TH1F& histo) const
-{
-    auto ranges = getter.GetRanges();
-    std::sort (ranges.begin(), ranges.end());
-           
-    TH1D hCentr("hCentr", "", ranges.size()-1, &(ranges.front()) );
-    
-    std::unique_ptr <TRandom> r {new TRandom};
-        
-    for (Int_t iBin = 0; iBin<histo.GetNbinsX(); ++iBin )
-    {
-        const Float_t Mult = histo.GetBinCenter(iBin+1);
-        
-        for (Int_t j = 0; j<histo.GetBinContent(iBin+1); ++j )
-        {
-            hCentr.Fill( getter.GetCentrality(Mult + (r->Rndm()-0.5)) );
-        }
-    }
-    hCentr.Write(name_ + "_qa");
-    if (ispdf_)
-    {
-        std::unique_ptr <TCanvas> c {new TCanvas("c", "", 1200, 800)};
-        hCentr.Draw();
-        c->SaveAs(name_ + "_qa" + ".pdf");
+void BordersFinderHelper::QA(const Getter& getter, const TH1F& histo) const {
+  auto ranges = getter.GetRanges();
+  std::sort(ranges.begin(), ranges.end());
+
+  TH1D hCentr("hCentr", "", ranges.size() - 1, &(ranges.front()));
+
+  std::unique_ptr<TRandom> r{new TRandom};
+
+  for(Int_t iBin = 0; iBin < histo.GetNbinsX(); ++iBin) {
+    const Float_t Mult = histo.GetBinCenter(iBin + 1);
+
+    for (Int_t j = 0; j < histo.GetBinContent(iBin + 1); ++j) {
+      hCentr.Fill(getter.GetCentrality(Mult + (r->Rndm() - 0.5)));
     }
+  }
+  hCentr.Write(name_ + "_qa");
+  if (ispdf_) {
+    std::unique_ptr<TCanvas> c{new TCanvas("c", "", 1200, 800)};
+    hCentr.Draw();
+    c->SaveAs(name_ + "_qa" + ".pdf");
+  }
 }
 
-void BordersFinderHelper::PlotHisto(const Getter& getter, TH1F& histo) const
-{
-    std::unique_ptr <TCanvas> c {new TCanvas("c", "", 1200, 800)};
-    histo.Draw();
-    
-    auto borders = getter.GetBorders();
-    TLine *line;
-
-    for (int i=0; i<borders.GetNbins(); ++i)
-    {
-        const float border = borders.GetBinLowEdge(i+1);
-        const int height = histo.GetBinContent( histo.FindBin(border) );
-        
-        line = new TLine(border, 0, border ,height); 
-        line->Draw("same");
-    }
-    c->SetLogy();
-    c->Write(name_ + "_histo_1d");
+void BordersFinderHelper::PlotHisto(const Getter& getter, TH1F& histo) const {
+  std::unique_ptr<TCanvas> c{new TCanvas("c", "", 1200, 800)};
+  histo.Draw();
+
+  const auto& borders = getter.GetBorders();
+  TLine* line;
+
+  for(int i = 0; i < borders.GetNbins(); ++i) {
+    const float border = borders.GetBinLowEdge(i + 1);
+    const int height = histo.GetBinContent(histo.FindBin(border));
 
-    if (ispdf_)
-        c->SaveAs(name_ + "_histo_1d" + ".pdf");
+    line = new TLine(border, 0, border, height);
+    line->Draw("same");
+  }
+  c->SetLogy();
+  c->Write(name_ + "_histo_1d");
 
-    delete line;
+  if (ispdf_)
+    c->SaveAs(name_ + "_histo_1d" + ".pdf");
+
+  delete line;
 }
 
-void BordersFinderHelper::PlotHisto2D(const Getter& getter, TH2F& histo, TF1& func) const
-{
-    std::unique_ptr <TCanvas> c {new TCanvas("c", "", 1200, 800)};
-    histo.Draw("colz");
-    func.Draw("same");
-    
-    const auto& borders = getter.GetBorders2D();
-    TLine *line{nullptr};
-
-    for (uint i=0; i<borders.size(); ++i)
-    {
-        const float x1 = 0.;
-        const float x2 = 1.;
-        
-        const float y1 = borders.at(i)[0] + borders.at(i)[1] * x1;
-        const float y2 = borders.at(i)[0] + borders.at(i)[1] * x2;
-        
-        line = new TLine(x1, y1, x2 ,y2); 
-        line->Draw("same");
+void BordersFinderHelper::PlotHisto2D(const Getter& getter, TH2F& histo, TF1& func) const {
+  gStyle->SetOptStat(0000);
+
+  std::unique_ptr<TCanvas> c{new TCanvas("c", "", 1000, 1000)};
+  histo.GetXaxis()->SetRangeUser(0, 1.1);
+  histo.GetYaxis()->SetRangeUser(0, 1.1);
+  histo.Draw("colz");
+  func.SetLineColor(kBlack);
+  func.Draw("same");
+
+  const auto& borders = getter.GetBorders2D();
+  TLine* line{nullptr};
+
+  for(uint i = 0; i < borders.size() - 1; ++i) {
+
+    float x1{0.};
+    float x2{1.};
+    float x{0.5};
+
+    for(int iter = 0; iter < 10; ++iter) {
+      x = (x1 + x2) / 2;
+
+      if((func.Eval(x1) - borders.at(i)[0] - borders.at(i)[1] * x1)
+        * (func.Eval(x) - borders.at(i)[0] - borders.at(i)[1] * x)
+        < 0) {
+        x2 = x;
+      } else {
+        x1 = x;
+      }
     }
+    x1 = x - 0.06 * borders.at(i)[1];
+    x2 = x + 0.06 * borders.at(i)[1];
+    if (x1 < 0) x1 = 0;
 
-    c->Write(name_ + "_histo_2d");
-    
-    if (ispdf_)
-        c->SaveAs(name_ + "_histo_2d" + ".pdf");
+    float y1 = borders.at(i)[0] + borders.at(i)[1] * x1;
+    float y2 = borders.at(i)[0] + borders.at(i)[1] * x2;
 
-    delete line;
-}
+    //       if (y1 < 0) y1=0;
+
+    line = new TLine(x1, y1, x2, y2);
+    line->SetLineWidth(2);
+    line->SetLineColor(kRed);
+    line->Draw("same");
+  }
 
-    
+  c->SetLogz(true);
+  c->Write(name_ + "_histo_2d");
 
+  if(ispdf_) {
+    c->SaveAs(name_ + "_histo_2d" + ".pdf");
+  }
+
+  delete line;
 }
+
+}// namespace Centrality

src/BordersFinderHelper.h → src/BordersFinderHelper.hpp

@@ -8,31 +8,30 @@
 #ifndef CENTRALITY_BORDERSFINDERHELPER_H
 #define CENTRALITY_BORDERSFINDERHELPER_H
 
+#include "TF1.h"
 #include "TH1.h"
 #include "TH2.h"
-#include "TF1.h"
 
-#include "Getter.h"
-
-namespace Centrality {
-    
-class BordersFinderHelper {
-
-public:
-    
-    BordersFinderHelper(){}
-    void QA(const Getter& getter, const TH1F &histo) const;
-    void PlotHisto(const Getter& getter, TH1F& histo) const;
-    void PlotHisto2D(const Getter& getter, TH2F& histo, TF1& func) const;
-
-    void SetName(const TString name) { name_ = name; }
-    void SetIsPdf(bool is = true) { ispdf_ = is; }
-private:
-    
-    TString name_{"test"};
-    bool ispdf_{false};
+#include "Getter.hpp"
+
+namespace Centrality{
+
+class BordersFinderHelper{
+
+ public:
+  BordersFinderHelper() = default;
+  void QA(const Getter& getter, const TH1F& histo) const;
+  void PlotHisto(const Getter& getter, TH1F& histo) const;
+  void PlotHisto2D(const Getter& getter, TH2F& histo, TF1& func) const;
+
+  void SetName(const TString& name) { name_ = name; }
+  void SetIsPdf(bool is = true) { ispdf_ = is; }
+
+ private:
+  TString name_{"test"};
+  bool ispdf_{false};
 };
 
-}
+}// namespace Centrality
 
-#endif //CENTRALITY_BORDERSFINDERHELPER_H
+#endif//CENTRALITY_BORDERSFINDERHELPER_H

src/CMakeLists.txt

+set(SOURCES
+        BordersFinder.cpp
+        BordersFinderHelper.cpp
+        Getter.cpp
+        BordersFinder2D.cpp
+        Fitter.cpp
+        FitterHelper.cpp
+        )
+
+string(REPLACE ".cpp" ".hpp" HEADERS "${SOURCES}")
+
+set(DICT_FILE_NAME G__${PROJECT_NAME})
+set(PCM_FILE_NAME lib${PROJECT_NAME})
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+
+add_library(Centrality SHARED ${SOURCES} ${DICT_FILE_NAME})
+ROOT_GENERATE_DICTIONARY(${DICT_FILE_NAME} ${HEADERS} LINKDEF CentralityLinkDef.hpp)
+target_link_libraries(Centrality ${ROOT_LIBRARIES})
+
+install(TARGETS Centrality EXPORT CentralityTargets
+        LIBRARY DESTINATION lib
+        ARCHIVE DESTINATION lib
+        RUNTIME DESTINATION bin
+        INCLUDES DESTINATION include
+        )
+
+install(
+        FILES
+        ${HEADERS}
+        DESTINATION
+        include/Centrality
+        COMPONENT
+        Devel
+)
+
+install(
+        FILES
+        "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}_rdict.pcm"
+        DESTINATION
+        lib
+        OPTIONAL
+)
+
+install(
+        FILES
+        "${CMAKE_CURRENT_BINARY_DIR}/${PCM_FILE_NAME}.rootmap"
+        DESTINATION
+        lib
+        OPTIONAL
+)

src/CentralityLinkDef.hpp

+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class Glauber::Fitter + ;
+#pragma link C++ class Centrality::Getter + ;
+
+#endif

src/Container.cpp

-#include "Container.h"
-// ClassImp(Centrality::Container)
-
-namespace Centrality {
-
-}

src/Container.h

-/** @file   Container.h
-    @class  Centrality::Container
-    @author Viktor Klochkov (klochkov44@gmail.com)
-    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
-    @brief  Class to store centrality information
-*/
-
-#ifndef CENTRALITY_CONTAINER_H
-#define CENTRALITY_CONTAINER_H
-
-#include <map>
-
-namespace Centrality {
-
-class Container{
-
-public:
-    
-    Container(){}
-    
-    void AddCentralityEstimator(uint num, float centrality)
-    {
-        centrality_.insert(std::make_pair(num, centrality));
-    }
-    
-    float GetCentrality(uint num) const 
-    { 
-        auto find = centrality_.find(num);        
-        return find != centrality_.end() ? find->second : -1;
-    }
-    
-    void Clear() { centrality_.clear(); }
-    
-private:
-    
-    std::map <uint, float> centrality_;
-    
-//     ClassDef(Container, 1);
-    
-};
-}
-
-#endif //CENTRALITY_CONTAINER_H

src/Fitter.cpp

+#include "Fitter.hpp"
+
+#include "TF1.h"
+#include "TFile.h"
+#include "TMath.h"
+#include "TRandom.h"
+#include "TTree.h"
+
+ClassImp(Glauber::Fitter)
+
+// -----   Default constructor   -------------------------------------------
+Glauber::Fitter::Fitter(std::unique_ptr<TTree> tree) {
+  fSimTree = std::move(tree);
+  std::cout << fSimTree->GetEntries() << std::endl;
+
+  if(!fSimTree) {
+    std::cout << "SetSimHistos: *** Error - " << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  fSimTree->SetBranchAddress("Npart", &fNpart);
+  fSimTree->SetBranchAddress("Ncoll", &fNcoll);
+}
+
+void Glauber::Fitter::Init(int nEntries) {
+  if (nEntries < 0 || nEntries > fSimTree->GetEntries()) {
+    std::cout << "Init: *** ERROR - number of entries < 0 or less that number of entries in input tree" << std::endl;
+    std::cout << "Init: *** number of entries in input tree = " << fSimTree->GetEntries() << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  const int NpartMax = int(fSimTree->GetMaximum("Npart"));
+  const int NcollMax = int(fSimTree->GetMaximum("Ncoll"));
+
+  fNpartHisto = TH1F("fNpartHisto", "Npart", NpartMax / fBinSize, 0, NpartMax);
+  fNcollHisto = TH1F("fNcollHisto", "Ncoll", NcollMax / fBinSize, 0, NcollMax);
+
+  for (int i = 0; i < nEntries; i++) {
+    fSimTree->GetEntry(i);
+    fNcollHisto.Fill(fNcoll);
+    fNpartHisto.Fill(fNpart);
+  }
+  std::cout << fSimTree->GetEntries() << std::endl;
+
+  fNbins = fDataHisto.GetNbinsX();
+
+  while (fDataHisto.GetBinContent(fNbins - 1) == 0)
+    fNbins--;
+
+  fNbins++;
+
+  const auto min = fDataHisto.GetXaxis()->GetXmin();
+  const auto max = fDataHisto.GetXaxis()->GetXmax();
+
+  fMaxValue = min + (max - min) * fNbins / fDataHisto.GetNbinsX();
+
+  std::cout << "fNbins = " << fNbins << std::endl;
+  std::cout << "fMaxValue = " << fMaxValue << std::endl;
+}
+
+double Glauber::Fitter::Nancestors(double f) const {
+  if(fMode == "Default") { return f * fNpart + (1 - f) * fNcoll; }
+  else if(fMode == "PSD") {
+    return f - fNpart;
+  } else if(fMode == "Npart") {
+    return TMath::Power(fNpart, f);
+  } else if(fMode == "Ncoll") {
+    return TMath::Power(fNcoll, f);
+  }
+
+  return -1.;
+}
+
+double Glauber::Fitter::NancestorsMax(double f) const {
+  const int NpartMax = fNpartHisto.GetXaxis()->GetXmax();// some magic
+  const int NcollMax = fNcollHisto.GetXaxis()->GetXmax();//TODO
+
+  if(fMode == "Default") { return f * NpartMax + (1 - f) * NcollMax; }
+  else if(fMode == "PSD") {
+    return NpartMax;
+  } else if(fMode == "Npart") {
+    return TMath::Power(NpartMax, f);
+  } else if(fMode == "Ncoll") {
+    return TMath::Power(NcollMax, f);
+  }
+
+  return -1.;
+}
+
+/*
+ * take Glauber MC data from fSimTree
+ * Populate fGlauberFitHisto with NBD x Na
+ */
+void Glauber::Fitter::SetGlauberFitHisto(double f, double mu, double k, int n, Bool_t Norm2Data) {
+  fGlauberFitHisto = TH1F("glaub", "", fNbins * 1.3, 0, 1.3 * fMaxValue);
+  fGlauberFitHisto.SetName(Form("glaub_%4.2f_%6.4f_%4.2f_%d", f, mu, k, n));
+
+  SetNBDhist(mu, k);
+
+  std::unique_ptr<TH1F> htemp{(TH1F*) fNbdHisto.Clone("htemp")};// WTF??? Not working without pointer
+  for(int i = 0; i < n; i++) {
+    fSimTree->GetEntry(i);
+    const int Na = int(Nancestors(f));
+
+    double nHits{0.};
+    for(int j = 0; j < Na; j++) {
+      nHits += int(htemp->GetRandom());
+    }
+    fGlauberFitHisto.Fill(nHits);
+  }
+
+  if (Norm2Data)
+    NormalizeGlauberFit();
+}
+
+void Glauber::Fitter::NormalizeGlauberFit() {
+
+  int fGlauberFitHistoInt{0};
+  int fDataHistoInt{0};
+
+  const int lowchibin = fFitMinBin;
+  const int highchibin = fFitMaxBin < fNbins ? fFitMaxBin : fNbins;
+
+  for (int i = lowchibin; i < highchibin; i++) {
+    fGlauberFitHistoInt += fGlauberFitHisto.GetBinContent(i + 1);
+    fDataHistoInt += fDataHisto.GetBinContent(i + 1);
+  }
+
+  const double ScaleFactor = (double) fDataHistoInt / fGlauberFitHistoInt;
+
+  //     std::cout << "Scale = " << Scale << std::endl;
+  fGlauberFitHisto.Scale(ScaleFactor);
+}
+
+/**
+ * 
+ * @param mu mean value of negative binominal distribution (we are looking for it)
+ * @param chi2 return value (indicates good match)
+ * @param mu_min lower search edge for mean value NBD
+ * @param mu_max upper search edge for mean value NBD
+ * @param f parameter of Na
+ * @param k NBD parameter
+ * @param nEvents
+ * @param nIter
+ */
+void Glauber::Fitter::FindMuGoldenSection(double* mu,
+                                          double* chi2,
+                                          double mu_min,
+                                          double mu_max,
+                                          double f,
+                                          double k,
+                                          int nEvents,
+                                          int nIter) {
+  const double phi{(1. + TMath::Sqrt(5)) / 2.};
+
+  /* left */
+  double mu_1 = mu_max - (mu_max - mu_min) / phi;
+
+  /* right */
+  double mu_2 = mu_min + (mu_max - mu_min) / phi;
+
+  SetGlauberFitHisto(f, mu_1, k, nEvents);
+  double chi2_mu1 = GetChi2();
+
+  SetGlauberFitHisto(f, mu_2, k, nEvents);
+  double chi2_mu2 = GetChi2();
+
+  for (int j = 0; j < nIter; j++) {
+    if (chi2_mu1 > chi2_mu2) {
+      mu_min = mu_1;
+      mu_1 = mu_2;
+      mu_2 = mu_min + (mu_max - mu_min) / phi;
+      chi2_mu1 = chi2_mu2;
+      SetGlauberFitHisto(f, mu_2, k, nEvents);
+      chi2_mu2 = GetChi2();
+    } else {
+      mu_max = mu_2;
+      mu_2 = mu_1;
+      mu_1 = mu_max - (mu_max - mu_min) / phi;
+      chi2_mu2 = chi2_mu1;
+      SetGlauberFitHisto(f, mu_1, k, nEvents);
+      chi2_mu1 = GetChi2();
+    }
+
+    std::cout << "mu1 = " << mu_1 << " mu2 = " << mu_2 << " chi2_mu1 = " << chi2_mu1 << " chi2_mu2 = " << chi2_mu2
+              << std::endl;
+  }
+
+  /* take min(mu) */
+  *mu = (chi2_mu1 < chi2_mu2) ? mu_1 : mu_2;
+  /* take min(chi2) */
+  *chi2 = (chi2_mu1 < chi2_mu2) ? chi2_mu1 : chi2_mu2;
+}
+
+/**
+ * Find the best match
+ *
+ * @param return value of best fit parameters
+ * @param f0 parameter of Na, for which chi2 will be calculated
+ * @param k0 lower search edge for NBD parameter
+ * @param k1 upper search edge for NBD parameter
+ * @param nEvents
+ */
+double Glauber::Fitter::FitGlauber(double* par, double f0, Int_t k0, Int_t k1, Int_t nEvents) {
+  double f_fit{-1};
+  double mu_fit{-1};
+  double k_fit{-1};
+  double Chi2Min{1e10};
+
+  const TString filename = Form("%s/fit_%4.2f_%d_%d_%d.root", fOutDirName.Data(), f0, k0, k1, fFitMinBin);
+
+  //     std::unique_ptr<TFile> file {TFile::Open(filename, "recreate")};
+  //     std::unique_ptr<TTree> tree {new TTree("test_tree", "tree" )};
+
+  TFile* file{TFile::Open(filename, "recreate")};
+  TTree* tree{new TTree("test_tree", "tree")};
+
+  TH1F h1("h1", "", fNbins, 0, fMaxValue);
+
+  double f, mu, k, chi2, sigma;
+
+  tree->Branch("histo", "TH1F", &h1);
+  tree->Branch("f", &f, "f/F");
+  tree->Branch("mu", &mu, "mu/F");
+  tree->Branch("k", &k, "k/F");
+  tree->Branch("chi2", &chi2, "chi2/F");
+  tree->Branch("sigma", &sigma, "sigma/F");
+
+  f = f0;
+  for (int j = k0; j <= k1; j++) {
+    mu = fMaxValue / NancestorsMax(f);
+    k = j;
+    const double mu_min = 0.7 * mu;
+    const double mu_max = 1.0 * mu;
+
+    FindMuGoldenSection(&mu, &chi2, mu_min, mu_max, f, k, nEvents, 10);
+    sigma = (mu / k + 1) * mu;
+    h1 = fGlauberFitHisto;
+
+    tree->Fill();
+
+    if (chi2 < Chi2Min) {
+      f_fit = f;
+      mu_fit = mu;
+      k_fit = k;
+      Chi2Min = chi2;
+      fBestFitHisto = fGlauberFitHisto;
+    }
+  }
+  tree->Write();
+  file->Write();
+  file->Close();
+
+  par[0] = f_fit;
+  par[1] = mu_fit;
+  par[2] = k_fit;
+
+  return Chi2Min;
+}
+
+/**
+ * Compare fGlauberFitHisto with fDataHisto
+ * @return chi2 value
+ */
+double Glauber::Fitter::GetChi2() const {
+  double chi2{0.0};
+
+  const int lowchibin = fFitMinBin;
+  const int highchibin = fFitMaxBin < fNbins ? fFitMaxBin : fNbins;
+
+  for (int i = lowchibin; i <= highchibin; ++i) {
+    if (fDataHisto.GetBinContent(i) < 1.0) continue;
+    const double error2 = TMath::Power(fDataHisto.GetBinError(i), 2) + TMath::Power(fGlauberFitHisto.GetBinError(i), 2);
+    const double diff2 = TMath::Power((fGlauberFitHisto.GetBinContent(i) - fDataHisto.GetBinContent(i)), 2) / error2;
+
+    chi2 += diff2;
+  }
+
+  chi2 = chi2 / (highchibin - lowchibin + 1);
+  return chi2;
+}
+
+/**
+ * Popuates histogram nbd_<mean>_<k> with values of NBD
+ * @param mu
+ * @param k
+ */
+void Glauber::Fitter::SetNBDhist(double mu, double k) {
+  // Interface for TH1F.
+  const int nBins = (mu + 1.) * 3 < 10 ? 10 : (mu + 1.) * 3;
+
+  fNbdHisto = TH1F("fNbdHisto", "", nBins, 0, nBins);
+  fNbdHisto.SetName(Form("nbd_%f_%f", mu, k));
+
+  for (int i = 0; i < nBins; ++i) {
+    const double val = NBD(i, mu, k);
+    if (val > 1e-10) fNbdHisto.SetBinContent(i + 1, val);
+    //         std::cout << "val " << val << std::endl;
+  }
+}
+
+/**
+ * Negative Binomial Distribution (by definition)
+ * @param n argument
+ * @param mu mean
+ * @param k argument
+ * @return NBD for a given parameters
+ */
+double Glauber::Fitter::NBD(double n, double mu, double k) const {
+  // Compute NBD.
+  double F;
+  double f;
+
+  if (n + k > 100.0) {
+    // log method for handling large numbers
+    F = TMath::LnGamma(n + k) - TMath::LnGamma(n + 1.) - TMath::LnGamma(k);
+    f = n * TMath::Log(mu / k) - (n + k) * TMath::Log(1.0 + mu / k);
+    F += f;
+    F = TMath::Exp(F);
+  } else {
+    F = TMath::Gamma(n + k) / (TMath::Gamma(n + 1.) * TMath::Gamma(k));
+    f = n * TMath::Log(mu / k) - (n + k) * TMath::Log(1.0 + mu / k);
+    f = TMath::Exp(f);
+    F *= f;
+  }
+
+  return F;
+}
+/**
+ * Creates histo with a given model parameter distribution
+ * @param range observable range
+ * @param name name of the MC-Glauber model parameter 
+ * @param par array with fit parameters
+ * @param Nevents
+ * @return pointer to the histogram 
+ */
+std::unique_ptr<TH1F> Glauber::Fitter::GetModelHisto(const double range[2],
+                                                     const TString& name,
+                                                     const double par[3],
+                                                     int nEvents) {
+  const double f = par[0];
+  const double mu = par[1];
+  const double k = par[2];
+
+  double modelpar{-999.};
+  fSimTree->SetBranchAddress(name, &modelpar);
+
+  SetNBDhist(mu, k);
+  std::unique_ptr<TH1F> hModel(new TH1F("hModel", "name", 100, fSimTree->GetMinimum(name), fSimTree->GetMaximum(name)));
+
+#pragma omp parallel for
+  for (int i = 0; i < nEvents; i++) {
+    fSimTree->GetEntry(i);
+    const int Na = int(Nancestors(f));
+    double nHits{0.};
+    for (int j = 0; j < Na; ++j) {
+      nHits += (int) fNbdHisto.GetRandom();
+    }
+
+    if (nHits > range[0] && nHits < range[1]) {
+      hModel->Fill(modelpar);
+    }
+  }
+
+  return hModel;
+}

src/Fitter.hpp

+/** @file   Fitter.h
+ *    @class  Glauber::Fitter
+ *    @author Viktor Klochkov (klochkov44@gmail.com)
+ *    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
+ *    @brief  Class to fit histo with Glauber based function
+ */
+
+#ifndef GlauberFitter_H
+#define GlauberFitter_H 1
+
+#include "TH1F.h"
+#include "TNamed.h"
+#include "TString.h"
+#include "TTree.h"
+
+#include <vector>
+// #include "TMinuit.h"
+
+namespace Glauber{
+class Fitter{
+
+ public:
+  /**   Default constructor   **/
+  Fitter() = default;;
+  explicit Fitter(std::unique_ptr<TTree> tree);
+  /**   Destructor   **/
+  virtual ~Fitter() = default;;
+
+  void Init(int nEntries);
+  void SetGlauberFitHisto(double f, double mu, double k, Int_t n = 10000, Bool_t Norm2Data = true);
+  void NormalizeGlauberFit();
+  //  void DrawHistos(Bool_t isSim = true, Bool_t isData = true, Bool_t isGlauber = false, Bool_t isNBD = false) {};
+
+  double FitGlauber(double* par, double f0, Int_t k0, Int_t k1, Int_t nEvents);
+  void FindMuGoldenSection(double* mu,
+                           double* chi2,
+                           double mu_min,
+                           double mu_max,
+                           double f,
+                           double k,
+                           Int_t nEvents = 10000,
+                           Int_t nIter = 5);
+
+  double GetChi2() const;
+
+  double NBD(double n, double mu, double k) const;
+  void SetNBDhist(double mu, double k);
+  double Nancestors(double f) const;
+  double NancestorsMax(double f) const;
+
+  std::unique_ptr<TH1F> GetModelHisto(const double range[2], const TString& name, const double par[3], Int_t nEvents);
+
+  //
+  //         Setters
+  //
+  void SetInputHisto(const TH1F& h) { fDataHisto = h; }
+  void SetFitMinBin(Int_t min) { fFitMinBin = min; }
+  void SetFitMaxBin(Int_t min) { fFitMaxBin = min; }
+  void SetNormMinBin(Int_t min) { fNormMinBin = min; }
+  void SetBinSize(Int_t size) { fBinSize = size; }
+  void SetOutDirName(const TString& name) { fOutDirName = name; }
+  void SetMode(const TString& mode) { fMode = mode; }
+
+  //
+  //         Getters
+  //
+  TH1F GetGlauberFitHisto() const { return fGlauberFitHisto; }
+  TH1F GetDataHisto() const { return fDataHisto; }
+  TH1F GetNBDHisto() const { return fNbdHisto; }
+  TH1F GetNpartHisto() const { return fNpartHisto; }
+  TH1F GetNcollHisto() const { return fNcollHisto; }
+  TH1F GetBestFiHisto() const { return fBestFitHisto; }
+
+ private:
+  /**   Data members  **/
+  TH1F fNpartHisto;
+  TH1F fNcollHisto;
+  TH1F fDataHisto;
+  TH1F fNbdHisto;
+  TH1F fGlauberFitHisto;
+  TH1F fBestFitHisto;
+
+  /* MC data */
+  std::unique_ptr<TTree> fSimTree{nullptr};
+
+  double fNpart{-1.};
+  double fNcoll{-1.};
+
+  double fMaxValue{-1.};
+
+  Int_t fNbins{-1};
+  Int_t fBinSize{-1};
+
+  Int_t fFitMinBin{-1};
+  Int_t fFitMaxBin{-1};
+
+  Int_t fNormMinBin{-1};
+
+  TString fMode{"Default"};
+
+  TString fOutDirName{""};
+ ClassDef(Fitter, 2);
+};
+}// namespace Glauber
+
+#endif

src/FitterHelper.cpp

+#include "FitterHelper.hpp"

src/FitterHelper.hpp

+/** @file   FitterHelper.h
+ *    @author Viktor Klochkov (klochkov44@gmail.com)
+ *    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
+ *    @brief  Methods for fit QA
+ */
+
+#ifndef GlauberFitterHelper_H
+#define GlauberFitterHelper_H 1
+
+#include "TCanvas.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "TLegend.h"
+#include "TPad.h"
+
+#include "Fitter.hpp"
+
+namespace Glauber{
+inline void DrawHistos(const Fitter& fit, Bool_t isSim, Bool_t isData, Bool_t isGlauber, Bool_t isNBD) {
+  std::unique_ptr<TCanvas> c1{new TCanvas("c1", "canvas", 1500, 900)};
+
+  c1->Divide(2, 2);
+
+  std::unique_ptr<TPad> c1_1{(TPad*) c1->GetListOfPrimitives()->FindObject("c1_1")};
+  std::unique_ptr<TPad> c1_2{(TPad*) c1->GetListOfPrimitives()->FindObject("c1_2")};
+  std::unique_ptr<TPad> c1_4{(TPad*) c1->GetListOfPrimitives()->FindObject("c1_4")};
+
+  c1_1->SetLogy(1);
+  c1_2->SetLogy(1);
+  c1_4->SetLogy(1);
+
+  /*const*/ TH1F hGlaub = fit.GetGlauberFitHisto();
+  /*const*/ TH1F hData = fit.GetDataHisto();
+  /*const*/ TH1F hNBD = fit.GetNBDHisto();
+  /*const*/ TH1F hNcoll = fit.GetNcollHisto();
+  /*const*/ TH1F hNpart = fit.GetNpartHisto();
+  /*const*/ TH1F hBestFit = fit.GetBestFiHisto();
+
+  std::unique_ptr<TFile> fOut{TFile::Open("glauber_qa.root", "recreate")};
+
+  if(isSim) {
+    c1->cd(1);
+    hNcoll.SetLineColor(2);
+
+    hNcoll.Draw();
+    hNpart.Draw("same");
+
+    std::unique_ptr<TLegend> legSim{new TLegend(0.6, 0.75, 0.75, 0.83)};
+    legSim->AddEntry(&hNpart, "Npart", "l");
+    legSim->AddEntry(&hNcoll, "hNcoll", "l");
+    legSim->Draw("same");
+
+    hNcoll.Write();
+    hNpart.Write();
+  }
+
+  if(isData) {
+    c1->cd(2);
+    hData.Draw();
+    hData.Write();
+    if(isGlauber) {
+      hBestFit.SetLineColor(kRed);
+      hBestFit.Draw("same");
+
+      std::unique_ptr<TLegend> legData{new TLegend(0.6, 0.75, 0.75, 0.83)};
+      legData->AddEntry(&hBestFit, "Fit", "l");
+      legData->AddEntry(&hData, "Data", "l");
+      legData->Draw("same");
+      hBestFit.Write();
+    }
+  }
+
+  if(isNBD) {
+    c1->cd(3);
+    hNBD.Draw();
+    hNBD.Write();
+  }
+
+  if(isGlauber) {
+    c1->cd(4);
+    hBestFit.Draw();
+  }
+
+  c1->Write();
+  c1->SaveAs("glauber.pdf");
+  fOut->Close();
+}
+
+}// namespace Glauber
+
+#endif

src/Getter.cpp

-#include "Getter.h"
+#include "Getter.hpp"
+
 #include <iostream>
 
 ClassImp(Centrality::Getter)
 
-namespace Centrality {
+namespace Centrality{
 
+double Getter::GetCentrality(double value) const {
+  if(!isinitialized_) {
+    std::cout << "Centrality::Getter is not initialized!" << std::endl;
+    exit(-1);
+  }
+  const int ibin = borders_.FindBin(value);
+  if(ibin == 0 || ibin > borders_.GetNbins()) {
+    return -1;
+  }
 
-float Getter::GetCentrality(float value) const 
-{
-    if (!isinitialized_) 
-    {
-        std::cout << "Centrality::Getter is not initialized!" << std::endl;
-        exit(-1);
-    }
-    const int ibin = borders_.FindBin(value);    
-    if ( ibin == 0 || ibin > borders_.GetNbins() )
-        return -1;
-    
-//     std::cout << value << " " << ranges_.at(ibin-1) << "  " << ranges_.at(ibin) << std::endl; 
-    
-    const float centrality = 0.5 * ( ranges_.at(ibin-1) + ranges_.at(ibin) );
-
-    return centrality;
+  //     std::cout << value << " " << ranges_.at(ibin-1) << "  " << ranges_.at(ibin) << std::endl;
+
+  const double centrality = 0.5 * (ranges_.at(ibin - 1) + ranges_.at(ibin));
+
+  return centrality;
 }
 
-float Getter::GetCentrality(float xvalue, float yvalue) const 
-{
-    if (!isinitialized2D_) 
-    {
-        std::cout << "Centrality::Getter is not initialized!" << std::endl;
-        exit(-1);
-    }
-    
-    xvalue /= xmax_;
-    yvalue /= ymax_;
-    
-    for (uint iborder=0; iborder<borders2d_.size()-1; ++iborder )
-    {
-        const float y1 = borders2d_.at(iborder)[0] + borders2d_.at(iborder)[1]*xvalue ;
-        const float y2 = borders2d_.at(iborder+1)[0] + borders2d_.at(iborder+1)[1]*xvalue ;
-        
-        if ( yvalue < y1 && yvalue > y2)
-            return 0.5 * ( ranges_.at(iborder-1) + ranges_.at(iborder) );
+double Getter::GetCentrality(double xvalue, double yvalue) const {
+  if(!isinitialized2D_) {
+    std::cout << "Centrality::Getter is not initialized!" << std::endl;
+    exit(-1);
+  }
+
+  xvalue /= xmax_;
+  yvalue /= ymax_;
+
+  for(uint iborder = 0; iborder < borders2d_.size() - 1; ++iborder) {
+    const double y1 = borders2d_.at(iborder)[0] + borders2d_.at(iborder)[1] * xvalue;
+    const double y2 = borders2d_.at(iborder + 1)[0] + borders2d_.at(iborder + 1)[1] * xvalue;
+
+    if(yvalue < y1 && yvalue > y2) {
+      return 0.5 * (ranges_.at(iborder - 1) + ranges_.at(iborder));
     }
-    
-    return -1.;
+  }
+
+  return -1.;
 }
 
+Getter* Getter::Create1DGetter(std::vector<double> borders) {
+  typedef decltype(ranges_)::value_type doubleing_t;
+
+  size_t n_borders = borders.size();
+  assert(n_borders > 2);
+
+  doubleing_t range_max = 100.f;
+  doubleing_t range_min = 0.f;
+  doubleing_t range_step = (range_max - range_min) / (n_borders - 1);
+
+  decltype(ranges_) ranges(n_borders);
+  for(uint i = 0; i < n_borders; ++i) {
+    auto rr = range_min + range_step * i;
+    ranges[i] = rr;
+
+    std::cout << rr << "%: " << borders[i] << std::endl;
+  }
+
+  TAxis ax_borders(static_cast<Int_t>(n_borders - 1), &(borders[0]));
+
+  auto getter = new Getter;
+  getter->ranges_ = std::move(ranges);
+  getter->borders_ = ax_borders;
+  getter->isinitialized_ = true;
+
+  return getter;
+}
 }

src/Getter.h

-/** @file   Getter.h
-    @class  Centrality::Getter
-    @author Viktor Klochkov (klochkov44@gmail.com)
-    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
-    @brief  Class to calculate centrality class
-*/
-
-#ifndef CENTRALITY_GETTER_H
-#define CENTRALITY_GETTER_H
-
-#include "vector"
-#include "array"
-
-#include "TObject.h"
-#include "TAxis.h"
-#include "TRandom.h"
-
-namespace Centrality {
-    
-class Getter : public TObject {
-
-public:
-    
-    Getter(){}
-    
-    float GetCentrality(float value) const;
-    float GetCentrality(float xvalue, float yvalue) const;
-    
-    void SetBorders(const std::vector<double> &borders)
-    {
-        borders_ = TAxis( borders.size()-1, &(borders[0]) );
-        isinitialized_ = true;
-    }
-    
-    const TAxis& GetBorders() const { return  borders_; };
-    const std::vector<float>& GetRanges() const { return  ranges_; };
-    
-    void SetRanges(const std::vector<float> &ranges) { ranges_ = ranges; }
-    void IsSpectator(bool is=true) { isspectator_ = is; }
-    
-    void AddBorder2D( const std::array<float,2> &border2D ) 
-    {       
-        borders2d_.push_back(border2D);
-        if (!isinitialized2D_) isinitialized2D_ = true;
-    }
-    
-    const std::vector <std::array<float,2>>& GetBorders2D() const { return  borders2d_; };
-    
-    void SetMax(float x, float y) { xmax_=x; ymax_=y; }
-    
-private:
-    
-    TAxis borders_;
-    std::vector <std::array<float,2>> borders2d_{};     
-    std::vector<float> ranges_{};
-
-    float xmax_{1.};
-    float ymax_{1.};    
-    
-    bool isspectator_{false};
-    bool isinitialized_{false};
-    bool isinitialized2D_{false};
-    
-    ClassDef(Getter, 2);
-    
-};
-}
-
-#endif //CENTRALITY_GETTER_H

src/Getter.hpp

+/** @file   Getter.h
+    @class  Centrality::Getter
+    @author Viktor Klochkov (klochkov44@gmail.com)
+    @author Ilya Selyuzhenkov (ilya.selyuzhenkov@gmail.com)
+    @brief  Class to calculate centrality class
+*/
+
+#ifndef CENTRALITY_GETTER_H
+#define CENTRALITY_GETTER_H
+
+#include "array"
+#include "vector"
+#include <cassert>
+
+#include "TAxis.h"
+#include "TObject.h"
+#include "TRandom.h"
+
+namespace Centrality{
+
+class Getter : public TObject{
+
+ public:
+  Getter() = default;
+
+  double GetCentrality(double value) const;
+  double GetCentrality(double xvalue, double yvalue) const;
+
+  void SetBorders(const std::vector<double>& borders) {
+    borders_ = TAxis(borders.size() - 1, &(borders[0]));
+    isinitialized_ = true;
+  }
+
+  const TAxis& GetBorders() const { return borders_; };
+  const std::vector<double>& GetRanges() const { return ranges_; };
+
+  void SetRanges(const std::vector<double>& ranges) { ranges_ = ranges; }
+  void IsSpectator(bool is = true) { isspectator_ = is; }
+
+  void AddBorder2D(const std::array<double, 2>& border2D) {
+    borders2d_.push_back(border2D);
+    if(!isinitialized2D_) { isinitialized2D_ = true; }
+  }
+
+  const std::vector<std::array<double, 2>>& GetBorders2D() const { return borders2d_; };
+
+  void SetMax(double x, double y) {
+    xmax_ = x;
+    ymax_ = y;
+  }
+
+  static Getter* Create1DGetter(std::vector<double> borders);
+
+ private:
+  TAxis borders_;
+  std::vector<std::array<double, 2>> borders2d_{};
+  std::vector<double> ranges_{};
+
+  double xmax_{1.};
+  double ymax_{1.};
+
+  bool isspectator_{false};
+  bool isinitialized_{false};
+  bool isinitialized2D_{false};
+
+ ClassDef(Getter, 2);
+};
+}// namespace Centrality
+
+#endif//CENTRALITY_GETTER_H

tasks/CMakeLists.txt

+# Create a main program using the library
+add_executable(main main.cpp)
+target_link_libraries(main Centrality ${ROOT_LIBRARIES})
+target_include_directories(main PUBLIC ${CMAKE_SOURCE_DIR}/src)
+
+add_executable(glauber glauber.cpp)
+target_link_libraries(glauber ${ROOT_LIBRARIES} Centrality)
+target_include_directories(glauber PUBLIC ${CMAKE_SOURCE_DIR}/src)
+
+if (1)
+    #    include(AnalysisTree)
+
+    add_executable(fill_centrality fill_centrality.cpp)
+    target_link_libraries(fill_centrality Centrality ATCentrality AnalysisTreeBase AnalysisTreeInfra ${ROOT_LIBRARIES})
+    target_include_directories(fill_centrality PUBLIC ${CMAKE_SOURCE_DIR}/src ${CMAKE_SOURCE_DIR}/at_interface ${PROJECT_INCLUDE_DIRECTORIES})
+
+    install(TARGETS fill_centrality RUNTIME DESTINATION bin)
+endif ()
+
+
+install(TARGETS glauber RUNTIME DESTINATION bin)
+install(TARGETS main RUNTIME DESTINATION bin)
\ No newline at end of file

tasks/fill_centrality.cpp

+#include "CentralityFiller.hpp"
+#include "AnalysisTree/TaskManager.hpp"
+
+using namespace AnalysisTree;
+
+void fill_centrality(const std::string& filelist, const std::string& centrality_file, const std::string& output) {
+
+  auto* man = TaskManager::GetInstance();
+  man->SetOutputName(output, "aTree");
+  man->SetOutputTreeConfig(OutputTreeConfig(eBranchWriteMode::kCopyTree));
+
+  auto* task = new CentralityFiller(centrality_file, "centr_getter_1d");
+  task->SetInput("RecEventHeader", "M");
+  task->SetOutput("AnaEventHeader", "centrality_tracks");
+
+  man->AddTask(task);
+
+  man->Init({filelist}, {"rTree"});
+  man->Run(-1);
+  man->Finish();
+}
+
+int main(int argc, char** argv) {
+  if(argc <= 2) {
+    std::cout << "Not enough arguments! Please use:" << std::endl;
+    std::cout << "   ./fill_centrality filelist centrality_file" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  const std::string filelist = argv[1];
+  const std::string centrality_file = argv[2];
+  const std::string output_file = "centrality.analysistree.root";
+
+  fill_centrality(filelist, centrality_file, output_file);
+  return EXIT_SUCCESS;
+}

tasks/glauber.cpp

+#include <chrono>
+#include <iostream>
+
+#include "Fitter.hpp"
+#include "FitterHelper.hpp"
+
+#include "TFile.h"
+#include "TH1.h"
+#include "TLegend.h"
+
+using namespace Glauber;
+
+int main(int argc, char* argv[]) {
+  if(argc < 2) {
+    std::cout << "Wrong number of parameters! Executable usage:" << std::endl;
+    std::cout << "   ./glauber f0 k0" << std::endl;
+    return -1;
+  }
+  const double_t f0 = atof(argv[1]);
+  const Int_t k0 = atoi(argv[2]);
+  const Int_t k1 = atoi(argv[2]);
+
+  // *****************************************
+  // Modify this part according to your needs
+  // *****************************************
+
+  ///  |   mode    |   function for Na     |
+  ///  |   Default |  Npart + (1-f)*Ncoll  |
+  ///  |   PSD     |     f - Npart         |
+  ///  |   Npart   |     Npart^f           |
+  ///  |   Ncoll   |     Ncoll^f           |
+  const TString mode = "Default";
+
+  const TString glauber_filename = "../input/glauber_auau_sigma_30_100k.root";// input files
+  const TString glauber_treename = "nt_Au_Au";
+  const TString in_filename = "../input/test_input.root";
+  const TString histoname = "hMreco";
+
+  const Int_t min_bin = 50;   // not fitting low multiplicity region due to trigger bias, etc
+  const Int_t max_bin = 10000;// very large number to fit the whole histo
+  const Int_t nevents = 100000;
+
+  const TString outdir = ".";
+  // *****************************************
+  // *****************************************
+
+  std::unique_ptr<TFile> glauber_file{TFile::Open(glauber_filename, "read")};
+  std::unique_ptr<TTree> glauber_tree{(TTree*) glauber_file->Get(glauber_treename)};
+
+  std::unique_ptr<TFile> f{TFile::Open(in_filename)};
+  TH1F* hdata = (TH1F*) f->Get(histoname);
+
+  Fitter fitter(std::move(glauber_tree));
+
+  fitter.SetMode(mode);
+  fitter.SetInputHisto(*hdata);
+  fitter.SetBinSize(1);
+  fitter.Init(nevents);
+
+  fitter.SetFitMinBin(min_bin);
+  fitter.SetFitMaxBin(max_bin);
+  fitter.SetOutDirName(outdir);
+
+  double par[3];
+
+  auto start = std::chrono::system_clock::now();
+
+  const double chi2 = fitter.FitGlauber(par, f0, k0, k1, nevents);
+
+  auto end = std::chrono::system_clock::now();
+  std::chrono::duration<double> elapsed_seconds = end - start;
+  std::cout << "elapsed time: " << elapsed_seconds.count() << " s\n";
+
+  std::cout << "f = " << par[0] << "    mu = " << par[1] << "    k = " << par[2] << "    chi2 = " << chi2 << std::endl;
+
+  Glauber::DrawHistos(fitter, true, true, true, true);
+
+  const double range[2] = {300, 350.};
+  std::unique_ptr<TH1F> hB(fitter.GetModelHisto(range, "B", par, 100000));
+  hB->SaveAs("b_test.root");
+
+  std::cout << "END!" << std::endl;
+
+  return 0;
+}

tasks/main.cpp

+#include "BordersFinder.hpp"
+#include "BordersFinder2D.hpp"
+#include "Getter.hpp"
+
+#include <array>
+#include <chrono>
+#include <iomanip>
+#include <iostream>
+#include <random>
+
+#include "TFile.h"
+#include "TH2.h"
+#include "TTree.h"
+#include <TROOT.h>
+#include <TSystem.h>
+
+int main(int argc, char** argv) {
+
+  auto start = std::chrono::system_clock::now();
+  ROOT::EnableImplicitMT(2);
+
+  if(argc <= 4) {
+    std::cout << "Not enough arguments! Please use:" << std::endl;
+    std::cout << "   ./main filename histoname is_spectator is2d" << std::endl;
+    return -1;
+  }
+
+  const char* filename{argv[1]};
+  const char* histoname{argv[2]};
+  const bool is_spectator = strcmp(argv[3], "true") == 0;
+  const bool is_2d = strcmp(argv[4], "true") == 0;
+
+  std::unique_ptr<TFile> fIn{TFile::Open(filename, "read")};
+  std::string outfilename = "test.root";
+  if(!is_2d) {
+    std::unique_ptr<TH1F> histo{(TH1F*) (fIn->Get(histoname))};
+    Centrality::BordersFinder bf;
+    bf.SetHisto(*histo);
+    bf.SetRanges(20, 0, 100);// number of bins, min, max value
+    //   bf.SetRanges( {0,10,30,60,100} );  // centrality bins borders with array
+    bf.IsSpectator(is_spectator);// true if impact parameter b correlated with estimator (spectators eneggy),
+    // false - anticorrelated (multiplicity of produced particles)
+
+    bf.FindBorders();
+    bf.SaveBorders(outfilename, "centr_getter_1d");
+  } else {
+    std::unique_ptr<TH2F> histo2d{(TH2F*) (fIn->Get(argv[2]))};
+    Centrality::BordersFinder2D bf;
+    bf.SetHisto2D(std::move(*histo2d));
+    std::unique_ptr<TH1F> h1d = bf.Convert();
+
+    bf.SetHisto(*h1d);
+    bf.SetRanges(20, 0, 100);// number of bins, min, max value
+    //   bf.SetRanges( {0,10,30,60,100} );  // centrality bins borders with array
+    bf.IsSpectator(is_spectator);// true if impact parameter b correlated with estimator (spectators energy),
+    // false - anticorrelated (multiplicity of produced particles)
+
+    bf.FindBorders();
+    bf.SaveBorders2D(outfilename, "centr_getter_2d");
+  }
+
+  auto end = std::chrono::system_clock::now();
+  std::chrono::duration<double> elapsed_seconds = end - start;
+  std::cout << "elapsed time: " << elapsed_seconds.count() << " s\n";
+  return 0;
+}