diff --git a/CMakeLists.txt b/CMakeLists.txt
index 04f3815ab447944457d6330fcdef90154e967bca..a06b7a726840ef7cf59040dbf31b22a151d9c403 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,7 +20,7 @@ if (CMAKE_BUILD_TYPE MATCHES DEBUG)
     set(CMAKE_VERBOSE_MAKEFILE ON)
 endif ()
 
-set(CMAKE_CXX_FLAGS_DEBUG "-O0 -ggdb -DDEBUG -D__DEBUG -Wall")
+set(CMAKE_CXX_FLAGS_DEBUG "-O0 -ggdb -DDEBUG -D__DEBUG -Wall -Wextra -fopenmp")
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -ftree-vectorize -ffast-math -DNODEBUG -fopenmp")
 set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELEASE} -ggdb")
 message(STATUS "Using CXX flags for ${CMAKE_BUILD_TYPE}: ${CMAKE_CXX_FLAGS_${CMAKE_BUILD_TYPE}}")
diff --git a/glauber/Fitter.cpp b/glauber/Fitter.cpp
index 0eea28ef80d65f625e62a71041fc3ee0101a1500..854fce304eb58e49429d5d322cf1c1080c77e763 100644
--- a/glauber/Fitter.cpp
+++ b/glauber/Fitter.cpp
@@ -51,8 +51,8 @@ void Glauber::Fitter::Init(int nEntries) {
 
   fNbins++;
 
-  const float min = fDataHisto.GetXaxis()->GetXmin();
-  const float max = fDataHisto.GetXaxis()->GetXmax();
+  const auto min = fDataHisto.GetXaxis()->GetXmin();
+  const auto max = fDataHisto.GetXaxis()->GetXmax();
 
   fMaxValue = min + (max - min) * fNbins / fDataHisto.GetNbinsX();
 
@@ -61,7 +61,7 @@ void Glauber::Fitter::Init(int nEntries) {
 
 }
 
-float Glauber::Fitter::Nancestors(float f) const {
+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);
@@ -70,7 +70,7 @@ float Glauber::Fitter::Nancestors(float f) const {
   return -1.;
 }
 
-float Glauber::Fitter::NancestorsMax(float f) const {
+double Glauber::Fitter::NancestorsMax(double f) const {
   const int NpartMax = fNpartHisto.GetXaxis()->GetXmax();  // some magic
   const int NcollMax = fNcollHisto.GetXaxis()->GetXmax(); //TODO
 
@@ -86,7 +86,7 @@ float Glauber::Fitter::NancestorsMax(float f) const {
  * take Glauber MC data from fSimTree
  * Populate fGlauberFitHisto with NBD x Na
  */
-void Glauber::Fitter::SetGlauberFitHisto(float f, float mu, float k, int n, Bool_t Norm2Data) {
+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));
 
@@ -97,7 +97,7 @@ void Glauber::Fitter::SetGlauberFitHisto(float f, float mu, float k, int n, Bool
     fSimTree->GetEntry(i);
     const int Na = int(Nancestors(f));
 
-    float nHits{0.};
+    double nHits{0.};
     for (int j = 0; j < Na; j++) {
       nHits += int(htemp->GetRandom());
     }
@@ -121,7 +121,7 @@ void Glauber::Fitter::NormalizeGlauberFit() {
     fDataHistoInt += fDataHisto.GetBinContent(i + 1);
   }
 
-  const float ScaleFactor = (float) fDataHistoInt / fGlauberFitHistoInt;
+  const double ScaleFactor = (double) fDataHistoInt / fGlauberFitHistoInt;
 
   //     std::cout << "Scale = " << Scale << std::endl;
   fGlauberFitHisto.Scale(ScaleFactor);
@@ -138,27 +138,27 @@ void Glauber::Fitter::NormalizeGlauberFit() {
  * @param nEvents
  * @param nIter
  */
-void Glauber::Fitter::FindMuGoldenSection(float *mu,
-                                          float *chi2,
-                                          float mu_min,
-                                          float mu_max,
-                                          float f,
-                                          float k,
+void Glauber::Fitter::FindMuGoldenSection(double *mu,
+                                          double *chi2,
+                                          double mu_min,
+                                          double mu_max,
+                                          double f,
+                                          double k,
                                           int nEvents,
                                           int nIter) {
-  const float phi{(1 + TMath::Sqrt(5)) / 2};
+  const double phi{(1. + TMath::Sqrt(5)) / 2.};
 
   /* left */
-  float mu_1 = mu_max - (mu_max - mu_min) / phi;
+  double mu_1 = mu_max - (mu_max - mu_min) / phi;
 
   /* right */
-  float mu_2 = mu_min + (mu_max - mu_min) / phi;
+  double mu_2 = mu_min + (mu_max - mu_min) / phi;
 
   SetGlauberFitHisto(f, mu_1, k, nEvents);
-  float chi2_mu1 = GetChi2();
+  double chi2_mu1 = GetChi2();
 
   SetGlauberFitHisto(f, mu_2, k, nEvents);
-  float chi2_mu2 = GetChi2();
+  double chi2_mu2 = GetChi2();
 
   for (int j = 0; j < nIter; j++) {
     if (chi2_mu1 > chi2_mu2) {
@@ -196,11 +196,11 @@ void Glauber::Fitter::FindMuGoldenSection(float *mu,
  * @param k1 upper search edge for NBD parameter
  * @param nEvents
  */
-float Glauber::Fitter::FitGlauber(float *par, Float_t f0, Int_t k0, Int_t k1, Int_t nEvents) {
-  float f_fit{-1};
-  float mu_fit{-1};
-  float k_fit{-1};
-  float Chi2Min{1e10};
+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);
 
@@ -212,7 +212,7 @@ float Glauber::Fitter::FitGlauber(float *par, Float_t f0, Int_t k0, Int_t k1, In
 
   TH1F h1("h1", "", fNbins, 0, fMaxValue);
 
-  float f, mu, k, chi2, sigma;
+  double f, mu, k, chi2, sigma;
 
   tree->Branch("histo", "TH1F", &h1);
   tree->Branch("f", &f, "f/F");
@@ -225,8 +225,8 @@ float Glauber::Fitter::FitGlauber(float *par, Float_t f0, Int_t k0, Int_t k1, In
   for (int j = k0; j <= k1; j++) {
     mu = fMaxValue / NancestorsMax(f);
     k = j;
-    const float mu_min = 0.7 * mu;
-    const float mu_max = 1.0 * mu;
+    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;
@@ -258,16 +258,16 @@ float Glauber::Fitter::FitGlauber(float *par, Float_t f0, Int_t k0, Int_t k1, In
  * Compare fGlauberFitHisto with fDataHisto
  * @return chi2 value
  */
-float Glauber::Fitter::GetChi2() const {
-  float chi2{0.0};
+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 float error2 = TMath::Power(fDataHisto.GetBinError(i), 2) + TMath::Power(fGlauberFitHisto.GetBinError(i), 2);
-    const float diff2 = TMath::Power((fGlauberFitHisto.GetBinContent(i) - fDataHisto.GetBinContent(i)), 2) / error2;
+    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;
   }
@@ -281,7 +281,7 @@ float Glauber::Fitter::GetChi2() const {
  * @param mu
  * @param k
  */
-void Glauber::Fitter::SetNBDhist(float mu, float k) {
+void Glauber::Fitter::SetNBDhist(double mu, double k) {
   // Interface for TH1F.
   const int nBins = (mu + 1.) * 3 < 10 ? 10 : (mu + 1.) * 3;
 
@@ -289,7 +289,7 @@ void Glauber::Fitter::SetNBDhist(float mu, float k) {
   fNbdHisto.SetName(Form("nbd_%f_%f", mu, k));
 
   for (int i = 0; i < nBins; ++i) {
-    const float val = NBD(i, mu, k);
+    const double val = NBD(i, mu, k);
     if (val > 1e-10) fNbdHisto.SetBinContent(i + 1, val);
     //         std::cout << "val " << val << std::endl;    
   }
@@ -302,10 +302,10 @@ void Glauber::Fitter::SetNBDhist(float mu, float k) {
  * @param k argument
  * @return NBD for a given parameters
  */
-float Glauber::Fitter::NBD(float n, float mu, float k) const {
+double Glauber::Fitter::NBD(double n, double mu, double k) const {
   // Compute NBD.
-  float F;
-  float f;
+  double F;
+  double f;
 
   if (n + k > 100.0) {
     // log method for handling large numbers
@@ -330,29 +330,25 @@ float Glauber::Fitter::NBD(float n, float mu, float k) const {
  * @param Nevents
  * @return pointer to the histogram 
  */
-std::unique_ptr<TH1F> Glauber::Fitter::GetModelHisto(const float range[2],
+std::unique_ptr<TH1F> Glauber::Fitter::GetModelHisto(const double range[2],
                                                      const TString &name,
-                                                     const float par[3],
+                                                     const double par[3],
                                                      int nEvents) {
-  const float f = par[0];
-  const float mu = par[1];
-  const float k = par[2];
+  const double f = par[0];
+  const double mu = par[1];
+  const double k = par[2];
 
-  float modelpar{-999.};
+  double modelpar{-999.};
   fSimTree->SetBranchAddress(name, &modelpar);
 
   SetNBDhist(mu, k);
-
-  //     TRandom random;  
-  //     random.SetSeed(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));
-    float nHits{0.};
+    double nHits{0.};
     for (int j = 0; j < Na; ++j) {
       nHits += (int) fNbdHisto.GetRandom();
     }
@@ -362,7 +358,7 @@ std::unique_ptr<TH1F> Glauber::Fitter::GetModelHisto(const float range[2],
     }
   }
 
-  return std::move(hModel);
+  return hModel;
 
 }
 
diff --git a/glauber/Fitter.h b/glauber/Fitter.h
index 5c64b371e0402586fece8c8a0071835567d7606a..07e02d26dc51678de254c57b0163fe3fd09942e6 100644
--- a/glauber/Fitter.h
+++ b/glauber/Fitter.h
@@ -28,28 +28,28 @@ class Fitter {
   virtual ~Fitter() = default;;
 
   void Init(int nEntries);
-  void SetGlauberFitHisto(Float_t f, Float_t mu, Float_t k, Int_t n = 10000, Bool_t Norm2Data = true);
+  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) {};
-
-  float FitGlauber(float *par, Float_t f0, Int_t k0, Int_t k1, Int_t nEvents);
-  void FindMuGoldenSection(Float_t *mu,
-                           Float_t *chi2,
-                           Float_t mu_min,
-                           Float_t mu_max,
-                           Float_t f,
-                           Float_t k,
+//  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);
 
-  Float_t GetChi2() const;
+  double GetChi2() const;
 
-  Float_t NBD(Float_t n, Float_t mu, Float_t k) const;
-  void SetNBDhist(Float_t mu, Float_t k);
-  float Nancestors(float f) const;
-  float NancestorsMax(float f) 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 Float_t range[2], const TString &name, const Float_t par[3], Int_t nEvents);
+  std::unique_ptr<TH1F> GetModelHisto(const double range[2], const TString &name, const double par[3], Int_t nEvents);
 
   //
   //         Setters
@@ -85,10 +85,10 @@ class Fitter {
   /* MC data */
   std::unique_ptr<TTree> fSimTree{nullptr};
 
-  Float_t fNpart{-1.};
-  Float_t fNcoll{-1.};
+  double fNpart{-1.};
+  double fNcoll{-1.};
 
-  Float_t fMaxValue{-1.};
+  double fMaxValue{-1.};
 
   Int_t fNbins{-1};
   Int_t fBinSize{-1};
diff --git a/glauber/main.cpp b/glauber/main.cpp
index 159951e24d1eb1c1eb0f0837a3bdfc8412c83068..0229eb2f3f8c68a8ec134b28a967f175f01470a5 100644
--- a/glauber/main.cpp
+++ b/glauber/main.cpp
@@ -16,7 +16,7 @@ int main(int argc, char *argv[]) {
     std::cout << "   ./glauber f0 k0" << std::endl;
     return -1;
   }
-  const Float_t f0 = atof(argv[1]);
+  const double_t f0 = atof(argv[1]);
   const Int_t k0 = atoi(argv[2]);
   const Int_t k1 = atoi(argv[2]);
 
@@ -61,11 +61,11 @@ int main(int argc, char *argv[]) {
   fitter.SetFitMaxBin(max_bin);
   fitter.SetOutDirName(outdir);
 
-  float par[3];
+  double par[3];
 
   auto start = std::chrono::system_clock::now();
 
-  const float chi2 = fitter.FitGlauber(par, f0, k0, k1, nevents);
+  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;
@@ -75,7 +75,7 @@ int main(int argc, char *argv[]) {
 
   Glauber::DrawHistos(fitter, true, true, true, true);
 
-  const float range[2] = {300, 350.};
+  const double range[2] = {300, 350.};
   std::unique_ptr<TH1F> hB(fitter.GetModelHisto(range, "B", par, 100000));
   hB->SaveAs("b_test.root");
 
diff --git a/src/BordersFinder.cpp b/src/BordersFinder.cpp
index 7b57873170c8f933c0a7935ff1ef3c312b07346b..c1937792ee8c5cd996db42cb5236301861ede96b 100644
--- a/src/BordersFinder.cpp
+++ b/src/BordersFinder.cpp
@@ -54,14 +54,14 @@ void BordersFinder::FindBorders() {
     
     for (int iBin=1; iBin<=histo_.GetNbinsX() && iSlice<ranges_.size() ; ++iBin)
     {
-        const float step = isSpectator_ ? ranges_.at(iSlice) : 100. - ranges_.at(iSlice);
+        const double step = isSpectator_ ? ranges_.at(iSlice) : 100. - ranges_.at(iSlice);
         const long int entriesNeeeded = step/100. * norm_;
         entriesCurrent += histo_.GetBinContent(iBin);
         
         if (entriesCurrent >= entriesNeeeded)
         {
-            const float ratio = histo_.GetBinContent(iBin)>0 ? (entriesCurrent - entriesNeeeded) / histo_.GetBinContent(iBin) : 0;
-            const float border = histo_.GetBinLowEdge(iBin) + histo_.GetBinWidth(iBin) * ratio;
+            const double ratio = histo_.GetBinContent(iBin)>0 ? (entriesCurrent - entriesNeeeded) / histo_.GetBinContent(iBin) : 0;
+            const double border = histo_.GetBinLowEdge(iBin) + histo_.GetBinWidth(iBin) * ratio;
             
             borders_.push_back(border);
             
diff --git a/src/BordersFinder.h b/src/BordersFinder.h
index c2c9647f67e0b41c3559d251b5c733c75174261c..f524605c40194af37480376791b6f35401ee0db0 100644
--- a/src/BordersFinder.h
+++ b/src/BordersFinder.h
@@ -27,8 +27,8 @@ class BordersFinder {
   void SetNormalization(long int norm) { norm_ = norm; }
   Double_t GetNormalization() const { return norm_; }
 
-  void SetRanges(const std::vector<float> &ranges) { ranges_ = ranges; }
-  void SetRanges(int n, float min, float max) {
+  void SetRanges(const std::vector<double> &ranges) { ranges_ = ranges; }
+  void SetRanges(int n, double min, double max) {
     ranges_.clear();
 //         ranges_.reserve(n+1);
 
@@ -44,7 +44,7 @@ class BordersFinder {
 
   void IsSpectator(bool is = true) { isSpectator_ = is; }
 
-  const std::vector<float> &GetRanges() const { return ranges_; }
+  const std::vector<double> &GetRanges() const { return ranges_; }
   const std::vector<double> &GetBorders() const { return borders_; }
   bool GetIsSpectator() const { return isSpectator_; }
 
@@ -53,7 +53,7 @@ class BordersFinder {
   TH1F histo_;
   Double_t norm_{-1};
 
-  std::vector<float> ranges_{};
+  std::vector<double> ranges_{};
   std::vector<double> borders_{};
 
   bool isSpectator_{false};
diff --git a/src/BordersFinder2D.cpp b/src/BordersFinder2D.cpp
index b7a8e7feae5dbf7491d8fbcbb178068866b203ec..9e682bc4f5f01393c0fd69d3b7446e02a9ff94f8 100644
--- a/src/BordersFinder2D.cpp
+++ b/src/BordersFinder2D.cpp
@@ -71,19 +71,19 @@ std::unique_ptr<TH1F> BordersFinder2D::Convert() {
  * @param norm2 seconsd line parametrization
  * @return number of entries (integral)
  */
-float BordersFinder2D::FindIntegral(const std::array<float, 2> &norm1, const std::array<float, 2> &norm2) {
-  float sum{0.};
+double BordersFinder2D::FindIntegral(const std::array<double, 2> &norm1, const std::array<double, 2> &norm2) {
+  double sum{0.};
 
   for (int iBinX = 1; iBinX <= histo2d_.GetNbinsX(); ++iBinX) {
     for (int iBinY = 1; iBinY <= histo2d_.GetNbinsY(); ++iBinY) {
-      const float entries = histo2d_.GetBinContent(iBinX, iBinY);
+      const auto entries = histo2d_.GetBinContent(iBinX, iBinY);
       if (entries == 0) continue;
 
-      const float x = histo2d_.GetXaxis()->GetBinCenter(iBinX);
-      const float y = histo2d_.GetYaxis()->GetBinCenter(iBinY);
+      const auto x = histo2d_.GetXaxis()->GetBinCenter(iBinX);
+      const auto y = histo2d_.GetYaxis()->GetBinCenter(iBinY);
 
-      const float ynorm1 = norm1[0] + norm1[1] * x;
-      const float ynorm2 = norm2[0] + norm2[1] * x;
+      const auto ynorm1 = norm1[0] + norm1[1] * x;
+      const auto ynorm2 = norm2[0] + norm2[1] * x;
 
       if (y < ynorm1 && y > ynorm2)
         sum += entries;
@@ -137,23 +137,23 @@ void BordersFinder2D::Fit2D(const TString &func) {
  * @param x argument
  * @return a0 and a1 parameters y = a0 + a1 * x
  */
-std::array<float, 2> BordersFinder2D::FindNorm(const std::vector<double> &par, float x) {
-  std::array<float, 2> ret{};
+std::array<double, 2> BordersFinder2D::FindNorm(const std::vector<double> &par, double x) {
+  std::array<double, 2> ret{};
   const auto dx = (histo2d_.GetXaxis()->GetXmax() - histo2d_.GetXaxis()->GetXmin()) / 10000.;
 
   /* left */
-  const float y1 = polN(par, x - dx);
+  const auto y1 = polN(par, x - dx);
   /* right */
-  const float y2 = polN(par, x + dx);
+  const auto y2 = polN(par, x + dx);
 
   // cx*x + cy*y + c == 0
 
   /* 1/df */
-  const float cx = 1 / (y2 - y1);
+  const auto cx = 1 / (y2 - y1);
   /* 1/2dx */
-  const float cy = 0.5 / dx;
+  const auto cy = 0.5 / dx;
 
-  const float c = -cx * x - cy * polN(par, x);
+  const auto c = -cx * x - cy * polN(par, x);
 
   ret[0] = -c / cy;
   ret[1] = -cx / cy;
diff --git a/src/BordersFinder2D.h b/src/BordersFinder2D.h
index 5dd80ca804b8e8b8b93f681f8623a1fab5a163b5..fa554590ad59cb9cca99b209235d59771874e7ff 100644
--- a/src/BordersFinder2D.h
+++ b/src/BordersFinder2D.h
@@ -30,8 +30,8 @@ class BordersFinder2D : public BordersFinder {
   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);
+  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);
 
   /**
@@ -41,9 +41,9 @@ class BordersFinder2D : public BordersFinder {
    * @param N order
    * @return
    */
-  static float polN(const std::vector<double> &par, float x) {
-    float res{0.};
-    float xn{1.};
+  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;
@@ -58,8 +58,8 @@ class BordersFinder2D : public BordersFinder {
 
   TString fitname_{""};
 
-  float xmax_{1.};
-  float ymax_{1.};
+  double xmax_{1.};
+  double ymax_{1.};
 
 //     ClassDef(BordersFinder2D, 1);
 
diff --git a/src/Container.h b/src/Container.h
index 8f300df2b9ce7d0185a577ab2ceb7580202904b0..1ed876d156c9cf330efc1b5c4237f2ae5490f09b 100644
--- a/src/Container.h
+++ b/src/Container.h
@@ -18,11 +18,11 @@ class Container {
 
   Container() = default;
 
-  void AddCentralityEstimator(uint num, float centrality) {
+  void AddCentralityEstimator(uint num, double centrality) {
     centrality_.insert(std::make_pair(num, centrality));
   }
 
-  float GetCentrality(uint num) const {
+  double GetCentrality(uint num) const {
     auto find = centrality_.find(num);
     return find != centrality_.end() ? find->second : -1;
   }
@@ -31,7 +31,7 @@ class Container {
 
  private:
 
-  std::map<uint, float> centrality_;
+  std::map<uint, double> centrality_;
 
 //     ClassDef(Container, 1);
 
diff --git a/src/Getter.cpp b/src/Getter.cpp
index 3701e93c3f04b70f0117ddb3d07cd4035d5f7079..af8fb461cd8ea8c2f73283959102df1279a632b0 100644
--- a/src/Getter.cpp
+++ b/src/Getter.cpp
@@ -5,7 +5,7 @@ ClassImp(Centrality::Getter)
 
 namespace Centrality {
 
-float Getter::GetCentrality(float value) const {
+double Getter::GetCentrality(double value) const {
   if (!isinitialized_) {
     std::cout << "Centrality::Getter is not initialized!" << std::endl;
     exit(-1);
@@ -16,12 +16,12 @@ float Getter::GetCentrality(float value) const {
 
 //     std::cout << value << " " << ranges_.at(ibin-1) << "  " << ranges_.at(ibin) << std::endl; 
 
-  const float centrality = 0.5 * (ranges_.at(ibin - 1) + ranges_.at(ibin));
+  const double centrality = 0.5 * (ranges_.at(ibin - 1) + ranges_.at(ibin));
 
   return centrality;
 }
 
-float Getter::GetCentrality(float xvalue, float yvalue) const {
+double Getter::GetCentrality(double xvalue, double yvalue) const {
   if (!isinitialized2D_) {
     std::cout << "Centrality::Getter is not initialized!" << std::endl;
     exit(-1);
@@ -31,8 +31,8 @@ float Getter::GetCentrality(float xvalue, float yvalue) const {
   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;
+    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));
@@ -42,17 +42,17 @@ float Getter::GetCentrality(float xvalue, float yvalue) const {
 }
 
 Getter *Getter::Create1DGetter(std::vector<double> borders) {
-  typedef decltype(ranges_)::value_type Floating_t;
+  typedef decltype(ranges_)::value_type doubleing_t;
 
   size_t n_borders = borders.size();
   assert(n_borders > 2);
 
-  Floating_t range_max = 100.f;
-  Floating_t range_min = 0.f;
-  Floating_t range_step = (range_max - range_min) / (n_borders - 1);
+  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 (int i = 0; i < n_borders; ++i) {
+  for (uint i = 0; i < n_borders; ++i) {
     auto rr = range_min + range_step * i;
     ranges[i] = rr;
 
diff --git a/src/Getter.h b/src/Getter.h
index 39a0a2d04d7df2259bca669e2ea5e02415a0eb5d..9864dcd291f6b4312886c49edf009afdc33a4d9b 100644
--- a/src/Getter.h
+++ b/src/Getter.h
@@ -24,8 +24,8 @@ class Getter : public TObject {
 
   Getter() = default;
 
-  float GetCentrality(float value) const;
-  float GetCentrality(float xvalue, float yvalue) const;
+  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]));
@@ -33,19 +33,19 @@ class Getter : public TObject {
   }
 
   const TAxis &GetBorders() const { return borders_; };
-  const std::vector<float> &GetRanges() const { return ranges_; };
+  const std::vector<double> &GetRanges() const { return ranges_; };
 
-  void SetRanges(const std::vector<float> &ranges) { ranges_ = ranges; }
+  void SetRanges(const std::vector<double> &ranges) { ranges_ = ranges; }
   void IsSpectator(bool is = true) { isspectator_ = is; }
 
-  void AddBorder2D(const std::array<float, 2> &border2D) {
+  void AddBorder2D(const std::array<double, 2> &border2D) {
     borders2d_.push_back(border2D);
     if (!isinitialized2D_) isinitialized2D_ = true;
   }
 
-  const std::vector<std::array<float, 2>> &GetBorders2D() const { return borders2d_; };
+  const std::vector<std::array<double, 2>> &GetBorders2D() const { return borders2d_; };
 
-  void SetMax(float x, float y) {
+  void SetMax(double x, double y) {
     xmax_ = x;
     ymax_ = y;
   }
@@ -55,11 +55,11 @@ class Getter : public TObject {
  private:
 
   TAxis borders_;
-  std::vector<std::array<float, 2>> borders2d_{};
-  std::vector<float> ranges_{};
+  std::vector<std::array<double, 2>> borders2d_{};
+  std::vector<double> ranges_{};
 
-  float xmax_{1.};
-  float ymax_{1.};
+  double xmax_{1.};
+  double ymax_{1.};
 
   bool isspectator_{false};
   bool isinitialized_{false};