diff --git a/Common/Tools/Multiplicity/README.md b/Common/Tools/Multiplicity/README.md
new file mode 100644
index 00000000000..c38f81a4f6e
--- /dev/null
+++ b/Common/Tools/Multiplicity/README.md
@@ -0,0 +1,17 @@
+# Multiplicity/centrality tools in O2/O2Physics
+
+This directory serves to aggregate all files necessary for multiplicity
+and centrality calibration going from pp to Pb-Pb. It offers functionality
+to perform simple slicing in percentiles, such as what is done in
+proton-proton collisions, as well as the Glauber + analytical NBD fitter
+used to perform anchoring and hadronic event distribution estimates in
+nucleus-nucleus collisions.
+
+## Files
+* `README.md` this readme
+* `CMakeLists.txt` definition of source files that need compiling
+* `multCalibrator.cxx/h` a class to do percentile slicing of a given histogram. Used for all systems.
+* `multMCCalibrator.cxx/h` a class to perform data-to-mc matching of average Nch.
+* `multGlauberNBDFitter.cxx/h` a class to do glauber fits.
+* `multModule.h` a class to perform calculations of multiplicity and centrality tables for analysis. Meant to be used inside the main core service wagon 'multcenttable'.
+
diff --git a/Common/Tools/Multiplicity/macros/README.md b/Common/Tools/Multiplicity/macros/README.md
new file mode 100644
index 00000000000..1044342cd8b
--- /dev/null
+++ b/Common/Tools/Multiplicity/macros/README.md
@@ -0,0 +1,69 @@
+# Example multiplicity calibration macros
+
+You will find some example macros in this directory that will allow for the calculation of a glauber fit and the corresponding calibration histograms.
+
+A simplified description of the procedure necessary to generate a
+Glauber + NBD fit to a certain histogram is described below.
+
+## Performing a Glauber + NBD fit
+
+### First step: calculation of Glauber MC sample
+
+The Glauber + NBD model assumes that the multiplicity / amplitude
+distribution that one intends to fit can be described as coming
+from a certain number N_{ancestor} of particle-emitting sources called 'ancestors'. Each ancestor is assumed to produce particles
+according to a negative binominal distribution. Further,
+the number of ancestors is assumed to be related to the basic
+glauber quantities N_{part} and N_{coll} via the formula:
+
+N_{ancestor} = f * N_{part}  + (1-f) N_{coll}
+
+Usually, the value f is fixed at 0.8, and thus the range of
+N_{ancestors} is typically 0-700 in Pb-Pb collisions.
+
+In order to allow for Glauber + NBD fitting, the correlation of
+(N_{part}, N_{coll}) needs to be known. For that purpose,
+a tree of Glauber MC needs to be generated using TGlauberMC using the relevant nuclei, which also involves choosing an appropriate nuclear profile.
+
+Once TGlauberMC has been used to produce a tree with N_{part} and
+N_{coll} values, their correlation needs to be saved to a 2D histogram that serves as input to the Glauber + NBD fitter used in
+O2/O2Physics. This is done using the macro called `saveCorrelation.C` contained in this directory, which produces a file named `baseHistos.root`. The file `saveCorrelation.C` serves as
+an example for the Pb-Pb case and minor adaptation may be needed
+in case other nuclei are to be used.
+
+### Second step: execution of Glauber + NBD fitter
+
+The fitting procedure is done via the macro ``. Because
+the numerical fitting utilizes a convolution of a N_{ancestor}
+distribution and NBD distributions, it will not be as fast
+as typical one-dimensional fits: typical times of 10-100 seconds
+are not unusual. The macro will produce an output file
+that contains:
+
+* the original fitted histogram
+* the actual glauber fit distribution, plotted all the way
+to zero multiplicity
+
+This output can then be used to extract percentiles.
+
+### Third step: calculation of percentile boundaries
+
+Once both the data and glauber fit distributions are known,
+the next step is to create a 'stitched' data/glauber distribution in
+which the distribution at low multiplicities follows
+the glauber fit and the distribution at higher multiplicities
+follows the data. The multiplicity value in which the switch from
+glauber to data takes place is called 'anchor point'.
+
+Because of the fact that this 'stitching' procedure may need to be tuned and the actual glauber fit is slow, the stitching is done
+in a separate macro called `runCalibration.C`. It is provided
+in this directory as well and it is the third and last step in calculating percentile boundaries. The macro will printout some
+key boundaries as well as save an output file with the calibration.
+
+*Bonus*: at the stage in which the glauber fit has been done
+and all information is available, a de-convolution process
+can be used to calculate the average N_{part} and N_{coll}
+in centrality slices. This functionality is also provided
+in O2Physics as part of the `multGlauberNBDFitter` and the
+`runCalibration.C` macro can optionally also perform that
+deconvolution. *Warning*: this procedure might take a mooment.
\ No newline at end of file
diff --git a/Common/Tools/Multiplicity/macros/runCalibration.C b/Common/Tools/Multiplicity/macros/runCalibration.C
new file mode 100644
index 00000000000..7a5e76c12d0
--- /dev/null
+++ b/Common/Tools/Multiplicity/macros/runCalibration.C
@@ -0,0 +1,172 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+//
+/// \file saveCorrelation.C
+/// \brief
+/// \author ALICE
+
+#include <iostream>
+
+/// @brief function to calibrate centrality
+/// @param lInputFileName name of input file.
+/// @param anchorPointPercentage anchor point percentage to use
+/// @param matchRange width of region in which data/glauber matching is to be done in rolling anchoring test
+/// @param doNpartNcoll wether or not to attempt calculating Npart, Ncoll in centrality bins
+void runCalibration(TString lInputFileName = "results/AR_544122_glauberNBD_ancestorMode2_hFT0C_BCs.root", double anchorPointPercentage = 90.0, double matchRange = 200.0, bool doNpartNcoll = false)
+{
+  TFile* file = new TFile(lInputFileName.Data(), "READ");
+  file->ls();
+
+  TH1F* hData = (TH1F*)file->Get("hV0MUltraFine");
+  TH1F* hGlauberParameters = (TH1F*)file->Get("hGlauberParameters");
+  TH1F* hGlauberFitRange = (TH1F*)file->Get("hGlauberFitRange");
+  hData->SetName("hData");
+  TH1F* hStitched = (TH1F*)hData->Clone("hStitched");
+  TH1F* hFit = (TH1F*)file->Get("hGlauber");
+
+  TCanvas* c1 = new TCanvas("c1", "", 800, 600);
+  c1->SetLeftMargin(0.17);
+  c1->SetBottomMargin(0.17);
+  c1->SetRightMargin(0.15);
+  c1->SetTopMargin(0.05);
+  c1->SetTicks(1, 1);
+  c1->SetLogz();
+  c1->SetFrameFillStyle(0);
+  c1->SetFillStyle(0);
+
+  cout << "Data bin width: " << hData->GetBinWidth(1) << endl;
+  cout << "Fit bin width: " << hFit->GetBinWidth(1) << endl;
+  cout << "Match range to use: " << matchRange << endl;
+
+  //____________________________________________
+  double anchorPointFraction = anchorPointPercentage / 100.f;
+  double anchorPoint = -1; // the anchor point value in raw
+
+  //____________________________________________
+  // doing partial integration up to certain point for finding anchor point bin
+  for (int ii = 1; ii < hData->GetNbinsX() + 1; ii++) {
+    // renormalize data curve
+    int bin1 = ii + 1;
+    int bin2 = hData->FindBin(hData->GetBinLowEdge(ii + 1) + matchRange + 1e-3);
+    double matchRangeData = hData->Integral(bin1, bin2);
+    double matchRangeFit = hFit->Integral(bin1, bin2);
+
+    // rescale fit to match in the vicinity of the region we're at
+    hFit->Scale(matchRangeData / matchRangeFit);
+
+    double integralFit = hFit->Integral(1, ii);
+    double integralData = hData->Integral(ii + 1, hData->GetNbinsX() + 1);
+    double integralAll = integralFit + integralData;
+
+    cout << "at bin #" << ii << ", integrated up to " << hData->GetBinLowEdge(ii + 1) << " fraction above this value is: " << integralData / integralAll << endl;
+    anchorPoint = hData->GetBinLowEdge(ii + 1);
+
+    if (integralData / integralAll < anchorPointFraction)
+      break;
+  }
+
+  //____________________________________________
+  for (int ii = 1; ii < hData->GetNbinsX() + 1; ii++) {
+    // renormalize data curve
+    if (hData->GetBinCenter(ii) < anchorPoint)
+      hStitched->SetBinContent(ii, hFit->GetBinContent(ii));
+  }
+
+  cout << "Anchor point determined to be: " << anchorPoint << endl;
+  cout << "Preparing stitched histogram ... " << endl;
+
+  hFit->SetLineColor(kRed);
+  hStitched->SetLineColor(kBlue);
+
+  hData->GetYaxis()->SetTitleSize(0.055);
+  hData->GetXaxis()->SetTitleSize(0.055);
+  hData->GetYaxis()->SetLabelSize(0.04);
+  hData->GetXaxis()->SetLabelSize(0.04);
+  hData->SetTitle("");
+  hData->Draw("hist");
+  hFit->Draw("hist same");
+  hStitched->Draw("hist same");
+
+  // All fine, let's try the calibrator
+  multCalibrator* lCalib = new multCalibrator("lCalib");
+  lCalib->SetAnchorPointPercentage(100.0f);
+  lCalib->SetAnchorPointRaw(-1e-6);
+
+  // Set standard Pb-Pb boundaries
+  lCalib->SetStandardOnePercentBoundaries();
+
+  TString calibFileName = lInputFileName.Data();
+  calibFileName.ReplaceAll("glauberNBD", "calibration");
+  TFile* fileCalib = new TFile(calibFileName.Data(), "RECREATE");
+
+  TH1F* hCalib = lCalib->GetCalibrationHistogram(hStitched, "hCalib");
+
+  TCanvas* c2 = new TCanvas("c2", "", 800, 600);
+  c2->SetLeftMargin(0.17);
+  c2->SetBottomMargin(0.17);
+  c2->SetRightMargin(0.15);
+  c2->SetTopMargin(0.05);
+  c2->SetTicks(1, 1);
+  //  c2->SetLogz();
+  c2->SetFrameFillStyle(0);
+  c2->SetFillStyle(0);
+
+  hCalib->GetYaxis()->SetTitleSize(0.055);
+  hCalib->GetXaxis()->SetTitleSize(0.055);
+  hCalib->GetYaxis()->SetLabelSize(0.04);
+  hCalib->GetXaxis()->SetLabelSize(0.04);
+  hCalib->SetTitle("");
+  hCalib->Draw();
+
+  fileCalib->cd();
+
+  hData->Write();
+  hCalib->Write();
+  hStitched->Write();
+  hFit->Write();
+
+  if (doNpartNcoll) {
+    cout << "Will now attempt to calculate % -> Np, Nc map..." << endl;
+
+    TProfile* hProfileNpart = new TProfile("hProfileNpart", "", 100, 0, 100);
+    TProfile* hProfileNcoll = new TProfile("hProfileNcoll", "", 100, 0, 100);
+    TH2F* h2dNpart = new TH2F("h2dNpart", "", 100, 0, 100, 500, -0.5f, 499.5f);
+    TH2F* h2dNcoll = new TH2F("h2dNcoll", "", 100, 0, 100, 3000, -0.5f, 2999.5);
+
+    // Replay
+    multGlauberNBDFitter* g = new multGlauberNBDFitter("lglau");
+    TF1* fitfunc = g->GetGlauberNBD();
+
+    // Step 1: open the (Npart, Ncoll) pair information, provide
+    TFile* fbasefile = new TFile("basehistos.root", "READ");
+    TH2D* hNpNc = (TH2D*)fbasefile->Get("hNpNc");
+    g->SetNpartNcollCorrelation(hNpNc);
+    g->InitializeNpNc();
+
+    fitfunc->SetParameter(0, hGlauberParameters->GetBinContent(1));
+    fitfunc->SetParameter(1, hGlauberParameters->GetBinContent(2));
+    fitfunc->SetParameter(2, hGlauberParameters->GetBinContent(3));
+    fitfunc->SetParameter(3, hGlauberParameters->GetBinContent(4));
+    fitfunc->SetParameter(4, hGlauberParameters->GetBinContent(5));
+
+    Double_t lMax = hData->GetBinLowEdge(hData->GetNbinsX() + 1);
+
+    // uncomment if Np Nc needed -> Warning, slow!
+    g->CalculateAvNpNc(hProfileNpart, hProfileNcoll, h2dNpart, h2dNcoll, hCalib, 0, lMax);
+
+    hProfileNpart->Write();
+    hProfileNcoll->Write();
+    h2dNpart->Write();
+    h2dNcoll->Write();
+  }
+
+  fileCalib->Write();
+}
diff --git a/Common/Tools/Multiplicity/macros/runGlauberFit.C b/Common/Tools/Multiplicity/macros/runGlauberFit.C
new file mode 100755
index 00000000000..fe33c209308
--- /dev/null
+++ b/Common/Tools/Multiplicity/macros/runGlauberFit.C
@@ -0,0 +1,386 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+//
+/// \file saveCorrelation.C
+/// \brief
+/// \author ALICE
+
+#include "multCalibrator.h"
+#include "multGlauberNBDFitter.h"
+
+#include "TCanvas.h"
+#include "TDirectory.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TGraph.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TLatex.h"
+#include "TLegend.h"
+#include "TLine.h"
+#include "TStopwatch.h"
+#include "TStyle.h"
+#include "TSystem.h"
+#include "TTree.h"
+
+#include <iostream>
+
+//________________________________________________________________
+Double_t FastIntegrate(TF1* f1, Double_t a, Double_t b, Int_t n = 5)
+{
+  // Do fast integration with N sampling points
+  const Int_t nc = n;
+  Double_t x[nc], y[nc];
+  Double_t lWidth = (b - a) / ((double)(n - 1));
+  for (Int_t ii = 0; ii < n; ii++) {
+    x[ii] = a + ((double)(ii)) * lWidth;
+    y[ii] = f1->Eval(x[ii]);
+  }
+  // Now go via trapezoids, please (this probably has a name)
+  Double_t lIntegral = 0;
+  for (Int_t ii = 0; ii < n - 1; ii++) {
+    lIntegral += 0.5 * lWidth * (y[ii] + y[ii + 1]);
+  }
+  return lIntegral / (b - a);
+}
+
+Double_t GetBoundaryForPercentile(TH1* histo, Double_t lPercentileRequested)
+{
+  // This function returns the boundary for a specific percentile.
+  Double_t lReturnValue = 0.0;
+  Double_t lPercentile = 100.0 - lPercentileRequested;
+
+  const Long_t lNBins = histo->GetNbinsX();
+  Double_t lCountDesired = lPercentile * histo->GetEntries() / 100;
+  Long_t lCount = 0;
+  for (Long_t ibin = 1; ibin < lNBins; ibin++) {
+    lCount += histo->GetBinContent(ibin);
+    if (lCount >= lCountDesired) {
+      // Found bin I am looking for!
+      Double_t lWidth = histo->GetBinWidth(ibin);
+      Double_t lLeftPercentile = 100. * (lCount - histo->GetBinContent(ibin)) / histo->GetEntries();
+      Double_t lRightPercentile = 100. * lCount / histo->GetEntries();
+
+      Double_t lProportion = (lPercentile - lLeftPercentile) / (lRightPercentile - lLeftPercentile);
+
+      lReturnValue = histo->GetBinLowEdge(ibin) + lProportion * lWidth;
+      break;
+    }
+  }
+  return lReturnValue;
+}
+
+/// @brief master glauber fit function
+/// @param lInputFileName input file name (from the grid, typically centrality-studies task)
+/// @param histogramName histogram name: histogram to use within the input file
+/// @param ancestorMode ancestor mode: 0: truncation, 1: rounding, 2: effective / non-integer (default: 2)
+/// @param lFreek free k: keep k value free (default Pb-Pb: fixed at 1.5)
+/// @param use_dMu_dNanc use dMu/dNanc: allow for a varying production of Nch vs ancestor if Nancestor is large. Models detector saturation in effective manner.
+/// @param lFreef free f: keep f value free (default Pb-Pb: fixed at 0.8)
+/// @param lfvalue f value: the value to use for fixed f
+/// @param outputFile name of output file
+int runGlauberFit(TString lInputFileName = "AnalysisResultsLHC24ar.root", TString histogramName = "hFT0C_BCs", int ancestorMode = 2, Bool_t lFreek = kFALSE, Bool_t use_dMu_dNanc = kFALSE, Bool_t lFreef = kFALSE, Float_t lfvalue = 0.800, TString outputFile = "output.root")
+{
+  gStyle->SetLineScalePS(1);
+  gStyle->SetOptStat(0);
+
+  cout << "Starting!" << endl;
+  TFile* file = new TFile(lInputFileName.Data(), "READ");
+  if (!file)
+    cout << "Problem with file!" << endl;
+  TH1F* hV0Mfine = 0x0;
+
+  hV0Mfine = (TH1F*)file->Get(Form("centrality-study/%s", histogramName.Data()));
+
+  // disregard bin zero
+  cout << "Received bin zero content: " << hV0Mfine->GetBinContent(0) << ", will set to zero..." << endl;
+  hV0Mfine->SetBinContent(0, 0);
+
+  if (!hV0Mfine)
+    cout << "Problem with histogram!" << endl;
+
+  cout << "Input histogram has been received successfully! Information: " << endl;
+
+  cout << "Counts: " << hV0Mfine->GetEntries() << endl;
+  cout << "NbinsX: " << hV0Mfine->GetNbinsX() << endl;
+  cout << "MaxX: " << hV0Mfine->GetBinLowEdge(hV0Mfine->GetNbinsX() + 1) << endl;
+
+  cout << "Creating output file..." << endl;
+  TString lProcessedFileName = lInputFileName.Data();
+  TString lkMode = "fixedK";
+  if (lFreek)
+    lkMode = "freeK";
+  TString lSaturationMode = "fixedMu";
+  if (use_dMu_dNanc)
+    lSaturationMode = "freeMu";
+
+  TFile* fOutput = new TFile(outputFile.Data(), "RECREATE");
+
+  TH1F* hV0M = (TH1F*)hV0Mfine->Clone("hV0M");
+  TH1F* hV0MUltraFine = (TH1F*)hV0Mfine->Clone("hV0MUltraFine");
+  hV0M->SetName("hV0M");
+  hV0M->SetTitle("");
+
+  //____________________________________________
+  // maximum fit range estimate (avoid tails)
+  // may need adjusting
+  Double_t lFitRangeMax = GetBoundaryForPercentile(hV0Mfine, 0.008);
+  cout << "Fit range max estimated from histogram: " << lFitRangeMax << endl;
+
+  // minimum fit range estimate (guess region that may be unfittable)
+  // may need adjusting
+  Double_t maxPercent = 0.01;
+  Double_t fractionOfMax = 0.012;
+  Double_t lFitRange = fractionOfMax * GetBoundaryForPercentile(hV0Mfine, maxPercent);
+
+  // adjust if low mult (Ntracks, typically)
+  Double_t maxRangeForTracks = 10000;
+  Double_t fractionOfMaxBroader = 0.02;
+  if (lFitRangeMax < maxRangeForTracks)
+    lFitRange = fractionOfMaxBroader * GetBoundaryForPercentile(hV0Mfine, maxPercent);
+
+  cout << "Fit range min estimated from histogram: " << lFitRange << endl;
+
+  //____________________________________________
+  // rebinning matters
+  int rebinFactor = 20;
+  if (lFitRangeMax < maxRangeForTracks)
+    rebinFactor = 1;
+
+  cout << "Creating rebinned histogram with rebin factor: " << rebinFactor << endl;
+  hV0M->Rebin(rebinFactor);
+
+  //____________________________________________
+  // simple plots for inspection
+  TCanvas* c1 = new TCanvas("c1", "", 1300, 900);
+  c1->SetFrameFillStyle(0);
+  c1->SetFillStyle(0);
+  c1->Divide(1, 2);
+  c1->cd(1)->SetFrameFillStyle(0);
+  c1->cd(1)->SetFillStyle(0);
+  c1->cd(2)->SetFrameFillStyle(0);
+  c1->cd(2)->SetFillStyle(0);
+
+  c1->cd(1);
+  c1->cd(1)->SetLogy();
+  c1->cd(1)->SetTicks(1, 1);
+  c1->cd(1)->SetPad(0, 0.5, 1, 1);
+  c1->cd(2)->SetPad(0, 0.0, 1, .5);
+
+  c1->cd(1)->SetBottomMargin(0.001);
+  c1->cd(1)->SetRightMargin(0.25);
+  c1->cd(1)->SetTopMargin(0.02);
+  c1->cd(1)->SetLeftMargin(0.07);
+
+  c1->cd(2)->SetBottomMargin(0.14);
+  c1->cd(2)->SetRightMargin(0.25);
+  c1->cd(2)->SetTopMargin(0.001);
+  c1->cd(2)->SetLeftMargin(0.07);
+  c1->cd(2)->SetTicks(1, 1);
+  c1->cd(1);
+
+  hV0M->GetXaxis()->SetRangeUser(0, lFitRangeMax);
+  hV0M->GetYaxis()->SetRangeUser(0.25, hV0M->GetMaximum() * 3);
+  hV0M->SetLineColor(kBlack);
+  hV0M->SetMarkerStyle(20);
+  hV0M->SetMarkerColor(kBlack);
+  hV0M->SetMarkerSize(0.5);
+  hV0M->GetYaxis()->SetTitleSize(0.07);
+  hV0M->GetYaxis()->SetLabelSize(0.05);
+  hV0M->GetYaxis()->SetTitle("Count");
+  hV0M->GetYaxis()->SetTitleOffset(0.5);
+  hV0M->GetXaxis()->SetLabelSize(0.05);
+  hV0M->GetXaxis()->SetTitleSize(0.06);
+  hV0M->GetXaxis()->SetTitle("FT0A+C Amplitude");
+  hV0M->GetYaxis()->SetTickLength(0.015);
+  hV0M->SetStats(0);
+  hV0M->Draw("E");
+
+  // Stand back! Imma gonna do GLAUBER FITTIN'
+  multGlauberNBDFitter* g = new multGlauberNBDFitter("lglau");
+  g->SetAncestorMode(ancestorMode);
+
+  // Step 1: open the (Npart, Ncoll) pair information, provide
+  TFile* fbasefile = new TFile("basehistos.root", "READ");
+  TH2D* hNpNc = (TH2D*)fbasefile->Get("hNpNc");
+  // return to proper scope
+  fOutput->cd();
+  g->SetNpartNcollCorrelation(hNpNc);
+  g->SetInputV0M(hV0M);
+  g->SetFitRange(lFitRange, lFitRangeMax);
+  // Step 3: go for it ...
+  g->SetNorm(1.53527e+08);
+  TString lString = "REM0";
+  g->SetFitOptions(lString.Data());
+  g->SetFitNpx(100000);
+
+  TF1* fitfunc = g->GetGlauberNBD();
+
+  //____________________________________________
+  //
+  // set initial fit parameters here
+  // may require manual tuning depending on data!
+
+  Double_t guessedMu = lFitRangeMax / 53968.4 * 0.175 * 3.53971e+02;
+  cout << "Guessed GlauberNBD mu value: " << guessedMu << endl;
+
+  fitfunc->SetParameter(0, guessedMu); // mu value
+  fitfunc->SetParLimits(0, 0.25 * guessedMu, guessedMu * 2);
+
+  if (!lFreek) {
+    fitfunc->FixParameter(1, 1.5); // k value
+  } else {
+    fitfunc->SetParLimits(1, 0.01, 35);
+    fitfunc->SetParameter(1, 1.5);
+  }
+  if (!lFreef) {
+    fitfunc->FixParameter(2, lfvalue); // f value
+  } else {
+    fitfunc->SetParLimits(2, 0.55, 0.97);
+    fitfunc->SetParameter(2, 0.800);
+  }
+  fitfunc->SetParLimits(3, 0.1e+5, 800e+10); // normalization
+  fitfunc->SetParameter(3, 0.80832e+08);
+
+  // dMu/dNanc
+  fitfunc->SetParameter(4, 0);
+  if (!use_dMu_dNanc) {
+    fitfunc->FixParameter(4, 0);
+  }
+
+  // fitfunc->SetParameter(4,-1.15443e-01);
+  // fitfunc->SetParameter(4,-4.55957e-02);
+
+  // dk/dNanc
+  fitfunc->FixParameter(5, 0);
+  // fitfunc->SetParameter(5,1.63590e-03);
+
+  // d2Mu/dNanc2
+  fitfunc->FixParameter(6, 0.0);
+  // fitfunc->SetParameter(6,4.02271e-05);
+
+  fitfunc->FixParameter(7, 0.0);
+  // fitfunc->SetParameter(7,-1.24349e-06);
+
+  //____________________________________________
+  // handle internals for fitting: needs to
+  // be done before the fit is attempted!
+  g->InitializeNpNc();
+  g->InitAncestor();
+
+  cout << "WILL NOW ATTEMPT GLAUBER FIT" << endl;
+  cout << "This will take a while. Please wait..." << endl;
+  Int_t lFitStatus = 0;
+  lFitStatus = g->DoFit();
+  Int_t lAttempts = 1;
+  Int_t lMaxAttempts = 10;
+  while (lAttempts < lMaxAttempts && lFitStatus == 0) {
+    // insist on fitting until it works
+    cout << "Attempting fit again (" << lAttempts << " attempt)..." << endl;
+    lFitStatus = g->DoFit();
+  }
+  cout << "Final fit status: " << lFitStatus << endl;
+
+  gStyle->SetOptStat(0);
+  // Do a ratio plot
+  TH1D* hGlauber = (TH1D*)hV0MUltraFine->Clone("hGlauber");
+  TH1D* hRatio = (TH1D*)hV0MUltraFine->Clone("hRatio");
+  hGlauber->Reset();
+
+  c1->cd(1);
+  fitfunc->SetLineColor(kRed);
+  fitfunc->SetLineWidth(2);
+  fitfunc->Draw("same");
+
+  cout << "Calculating glauber function histogram with the same binning as data input... please wait..." << endl;
+  for (Int_t ii = 1; ii < hGlauber->GetNbinsX() + 1; ii++) {
+    Double_t lFuncVal = FastIntegrate(fitfunc, hGlauber->GetBinLowEdge(ii), hGlauber->GetBinLowEdge(ii + 1), 4);
+    hGlauber->SetBinContent(ii, lFuncVal);
+    Int_t printEveryThisManyBins = 500;
+    if (ii % printEveryThisManyBins == 0) {
+      cout << "At integration #" << ii << "/" << hGlauber->GetNbinsX() + 1 << "..." << endl;
+    }
+  }
+  cout << "Glauber function evaluated. Should go quickly now." << endl;
+
+  c1->cd(2);
+  Float_t lLoRangeRatio = 0.35;
+  Float_t lHiRangeRatio = 1.65;
+  hRatio->Divide(hGlauber);
+  hRatio->GetYaxis()->SetTitle("Data/Fit");
+  hRatio->GetXaxis()->SetTitle("FT0C amplitude");
+  hRatio->GetYaxis()->SetTitleSize(0.055);
+  hRatio->GetYaxis()->SetTitleOffset(0.7);
+  hRatio->GetXaxis()->SetTitleSize(0.055);
+  hRatio->GetYaxis()->SetLabelSize(0.045);
+  hRatio->GetXaxis()->SetLabelSize(0.045);
+  hRatio->GetYaxis()->SetRangeUser(lLoRangeRatio, lHiRangeRatio);
+  hRatio->GetXaxis()->SetRangeUser(0, lFitRangeMax);
+  hRatio->SetMarkerStyle(20);
+  hRatio->SetMarkerColor(kGray + 2);
+  hRatio->SetLineColor(kGray + 2);
+  // hRatio->SetMarkerSize(1.0);
+  hRatio->SetMarkerSize(.7);
+  hRatio->SetStats(0);
+  // hRatioWide->SetStats(0);
+
+  hRatio->Draw("hist");
+
+  TLine* line = new TLine(0, 1, lFitRangeMax, 1);
+  line->SetLineStyle(7);
+  line->SetLineColor(kGray + 1);
+  line->Draw();
+
+  TLine* lFitRangeLine = new TLine(lFitRange, lLoRangeRatio, lFitRange, 0.9);
+  lFitRangeLine->SetLineColor(kBlue);
+  lFitRangeLine->SetLineWidth(1);
+  lFitRangeLine->SetLineStyle(2);
+  lFitRangeLine->Draw();
+
+  TH1D* hRatioGrayed = (TH1D*)hRatio->Clone("hRatioGrayed");
+  hRatioGrayed->SetMarkerColor(kGray + 2);
+  hRatioGrayed->SetLineColor(kGray + 2);
+  hRatioGrayed->Draw("same");
+
+  hRatio->SetLineWidth(1);
+  hRatio->Draw("same hist");
+  // hRatioWide->Draw("same");
+
+  c1->cd(1);
+  TLatex* lat = new TLatex();
+  lat->SetNDC();
+  Float_t lPosText = 0.76;
+  Float_t lYShift = 0.25;
+  lat->SetTextSize(0.042);
+
+  // save the glauber parameters explicitly
+  TH1D* hGlauberParameters = new TH1D("hGlauberParameters", "", 10, 0, 10);
+  TH1D* hGlauberFitRange = new TH1D("hGlauberFitRange", "", 10, 0, 10);
+
+  // fitfunc
+  hGlauberParameters->SetBinContent(1, fitfunc->GetParameter(0));
+  hGlauberParameters->SetBinContent(2, fitfunc->GetParameter(1));
+  hGlauberParameters->SetBinContent(3, fitfunc->GetParameter(2));
+  hGlauberParameters->SetBinContent(4, fitfunc->GetParameter(3));
+  hGlauberParameters->SetBinContent(5, fitfunc->GetParameter(4));
+  hGlauberParameters->Write();
+
+  Double_t lLoRangeGlauber, lHiRangeGlauber;
+  fitfunc->GetRange(lLoRangeGlauber, lHiRangeGlauber);
+  hGlauberFitRange->SetBinContent(1, lLoRangeGlauber);
+  hGlauberFitRange->SetBinContent(2, lHiRangeGlauber);
+  hGlauberFitRange->Write();
+
+  hRatio->Write();
+  fOutput->Write();
+
+  return 0;
+}
diff --git a/Common/Tools/Multiplicity/macros/saveCorrelation.C b/Common/Tools/Multiplicity/macros/saveCorrelation.C
new file mode 100644
index 00000000000..b23415639ba
--- /dev/null
+++ b/Common/Tools/Multiplicity/macros/saveCorrelation.C
@@ -0,0 +1,64 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+//
+/// \file saveCorrelation.C
+/// \brief This file provides a simple macro to convert TGlauberMC output tuples
+/// into the 2D correlation histogram necessary for the ALICE machinery
+/// that performs Glauber + NBD fits.
+/// \author ALICE
+
+#include <iostream>
+
+/// @brief function to save Npart x Ncoll correlation to file for glauber fits
+/// @param filename input TGlauberMC ntuple file
+/// @param outputFile output file for Npart x Ncoll correlation TH2D
+void saveCorrelation(TString filename = "gmc-PbPb-snn68.21-md0.40-nd-1.0-rc1-smax99.0.root", TString outputFile = "basehistos.root")
+{
+  TFile* fin = new TFile(filename.Data(), "READ");
+  TNtuple* ntup = (TNtuple*)fin->Get("nt_Pb_Pb");
+
+  // try other Pb nuclear profiles in case "Pb" - "Pb" not found
+  if (!ntup) {
+    ntup = (TNtuple*)fin->Get("nt_Pbpn_Pbpn");
+  }
+  if (!ntup) {
+    ntup = (TNtuple*)fin->Get("nt_PbpnVar1_PbpnVar1");
+  }
+  if (!ntup) {
+    ntup = (TNtuple*)fin->Get("nt_PbpnVar2_PbpnVar2");
+  }
+  if (!ntup) {
+    ntup = (TNtuple*)fin->Get("nt_PbpnVar3_PbpnVar3");
+  }
+  if (!ntup) {
+    ntup = (TNtuple*)fin->Get("nt_PbpnVar4_PbpnVar4");
+  }
+
+  if (!ntup) {
+    cout << "No tree found!" << endl;
+    return;
+  }
+
+  cout << "Glauber tree entries: " << ntup->GetEntries() << endl;
+
+  TFile* fout = new TFile(outputFile.Data(), "RECREATE");
+
+  // 2D correlation plot necessary for Glauber + NBD fitting
+  // The provided range should be enough for Pb-Pb
+  TH2D* hNpNc = new TH2D("hNpNc", "", 500, -0.5, 499.5, 2500, -0.5, 2499.5);
+
+  // let's draw this on screen for inspection
+  TCanvas* c1 = new TCanvas("c1", "", 800, 600);
+  c1->SetTicks(1, 1);
+  c1->SetLogz();
+  ntup->Draw("Ncoll:Npart>>hNpNc", "", "colz");
+  fout->Write();
+}