diff --git a/PWGLF/DataModel/LFStrangenessTables.h b/PWGLF/DataModel/LFStrangenessTables.h
index 057ff318f03..56322a0e47c 100644
--- a/PWGLF/DataModel/LFStrangenessTables.h
+++ b/PWGLF/DataModel/LFStrangenessTables.h
@@ -440,8 +440,10 @@ DECLARE_SOA_TABLE(StraMCCollMults_000, "AOD", "STRAMCCOLLMULTS", //! MC collisio
                   mult::MultMCFT0A, mult::MultMCFT0C, mult::MultMCNParticlesEta05, mult::MultMCNParticlesEta08, mult::MultMCNParticlesEta10, o2::soa::Marker<2>);
 DECLARE_SOA_TABLE_VERSIONED(StraMCCollMults_001, "AOD", "STRAMCCOLLMULTS", 1, //! MC collision multiplicities
                             mult::MultMCFT0A, mult::MultMCFT0C, mult::MultMCNParticlesEta05, mult::MultMCNParticlesEta08, mult::MultMCNParticlesEta10, stramccollision::TotalMultMCParticles);
+DECLARE_SOA_TABLE_VERSIONED(StraMCCollMults_002, "AOD", "STRAMCCOLLMULTS", 2, //! MC collision multiplicities
+                            mult::MultMCFT0A, mult::MultMCFT0C, mult::MultMCFV0A, mult::MultMCFDDA, mult::MultMCFDDC, mult::MultMCNParticlesEta05, mult::MultMCNParticlesEta08, mult::MultMCNParticlesEta10, stramccollision::TotalMultMCParticles);
 
-using StraMCCollMults = StraMCCollMults_001;
+using StraMCCollMults = StraMCCollMults_002;
 
 namespace dautrack
 {
diff --git a/PWGLF/DataModel/mcCentrality.h b/PWGLF/DataModel/mcCentrality.h
index 94544949017..2bb14880953 100644
--- a/PWGLF/DataModel/mcCentrality.h
+++ b/PWGLF/DataModel/mcCentrality.h
@@ -20,6 +20,7 @@
 #ifndef PWGLF_DATAMODEL_MCCENTRALITY_H_
 #define PWGLF_DATAMODEL_MCCENTRALITY_H_
 
+// O2 includes
 #include "Common/DataModel/Centrality.h"
 
 #include <Framework/ASoA.h>
@@ -38,8 +39,12 @@ DECLARE_SOA_TABLE(McCentFV0As, "AOD", "MCCENTFV0A", o2::soa::Marker<1>, cent::Ce
 DECLARE_SOA_TABLE(McCentFT0Ms, "AOD", "MCCENTFT0M", o2::soa::Marker<2>, cent::CentFT0M);
 DECLARE_SOA_TABLE(McCentFT0As, "AOD", "MCCENTFT0A", o2::soa::Marker<3>, cent::CentFT0A);
 DECLARE_SOA_TABLE(McCentFT0Cs, "AOD", "MCCENTFT0C", o2::soa::Marker<4>, cent::CentFT0C);
-DECLARE_SOA_TABLE(McCentFDDMs, "AOD", "MCCENTFDDM", o2::soa::Marker<5>, cent::CentFDDM);
-DECLARE_SOA_TABLE(McCentNTPVs, "AOD", "MCCENTNTPV", o2::soa::Marker<6>, cent::CentNTPV);
+DECLARE_SOA_TABLE(McCentFT0CVariant1s, "AOD", "MCCENTFT0CVAR1", o2::soa::Marker<5>, cent::CentFT0CVariant1);
+DECLARE_SOA_TABLE(McCentFT0CVariant2s, "AOD", "MCCENTFT0CVAR2", o2::soa::Marker<6>, cent::CentFT0CVariant2);
+DECLARE_SOA_TABLE(McCentFDDMs, "AOD", "MCCENTFDDM", o2::soa::Marker<7>, cent::CentFDDM);
+DECLARE_SOA_TABLE(McCentNTPVs, "AOD", "MCCENTNTPV", o2::soa::Marker<8>, cent::CentNTPV);
+DECLARE_SOA_TABLE(McCentNGlobals, "AOD", "MCCENTNGLOBAL", o2::soa::Marker<9>, cent::CentNGlobal);
+DECLARE_SOA_TABLE(McCentMFTs, "AOD", "MCCENTMFT", o2::soa::Marker<10>, cent::CentMFT);
 
 } // namespace o2::aod
 
diff --git a/PWGLF/TableProducer/Common/mcCentrality.cxx b/PWGLF/TableProducer/Common/mcCentrality.cxx
index 3e017c5808b..2a73fa6f785 100644
--- a/PWGLF/TableProducer/Common/mcCentrality.cxx
+++ b/PWGLF/TableProducer/Common/mcCentrality.cxx
@@ -17,28 +17,20 @@
 /// \brief  Task to produce the table for the equalized multiplicity into centrality bins
 ///
 
-#include "PWGLF/DataModel/mcCentrality.h"
+#include "PWGLF/Utils/mcCentralityModule.h"
 
-#include "PWGLF/Utils/inelGt.h"
+#include "Common/DataModel/Multiplicity.h"
 
 #include <CCDB/BasicCCDBManager.h>
-#include <CCDB/CcdbApi.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
 #include <Framework/Configurable.h>
 #include <Framework/HistogramRegistry.h>
-#include <Framework/HistogramSpec.h>
 #include <Framework/InitContext.h>
-#include <Framework/O2DatabasePDGPlugin.h>
 #include <Framework/OutputObjHeader.h>
 #include <Framework/runDataProcessing.h>
 
-#include <TFile.h>
-#include <TH1.h>
-
-#include <chrono>
-#include <cstdint>
 #include <string>
 
 using namespace o2;
@@ -47,146 +39,37 @@ using namespace o2::framework::expressions;
 
 /// Task to produce the response table
 struct McCentrality {
-
-  // Tables to produce
-  Produces<aod::McCentFV0As> centFV0A;
-  Produces<aod::McCentFT0Ms> centFT0M;
-  Produces<aod::McCentFT0As> centFT0A;
-  Produces<aod::McCentFT0Cs> centFT0C;
-  Produces<aod::McCentFDDMs> centFDDM;
-  Produces<aod::McCentNTPVs> centNTPV;
-
   // Input parameters
+  o2::framework::Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
   Service<o2::ccdb::BasicCCDBManager> ccdb;
-  Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
-  Configurable<int64_t> ccdbTimestamp{"ccdbTimestamp", -1, "timestamp of the object used to query in CCDB the detector response. If 0 the object corresponding to the run number is used, if < 0 the latest object is used"};
-  Configurable<std::string> path{"path", "/tmp/InputCalibMC.root", "path to calib file or ccdb path if begins with ccdb://"};
-  Configurable<bool> selectPrimaries{"selectPrimaries", true, "Select only primary particles"};
-  Service<o2::framework::O2DatabasePDG> pdgDB;
-  ConfigurableAxis binsPercentile{"binsPercentile", {VARIABLE_WIDTH, 0, 0.001, 0.01, 1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}, "Binning of the percentile axis"};
-  ConfigurableAxis binsMultiplicity{"binsMultiplicity", {1000, 0, 5000}, "Binning of the multiplicity axis"};
-  Configurable<bool> fillFt0M{"fillFt0M", true, "Fills the FT0M histogram"};
-  Configurable<bool> fillFt0A{"fillFt0A", false, "Fills the FT0A histogram"};
-  Configurable<bool> fillFt0C{"fillFt0C", false, "Fills the FT0C histogram"};
-  Configurable<bool> doNotCrashOnNull{"doNotCrashOnNull", false, "If ccdb object does not exist, fill with dummy values"};
-
-  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
-  TH1F* h1dFT0M = nullptr;
-  TH1F* h1dFT0A = nullptr;
-  TH1F* h1dFT0C = nullptr;
-  // TH1F* h1dFDD;
-  // TH1F* h1dNTP;
+  o2::pwglf::mccentrality::products products;
+  o2::pwglf::mccentrality::coreConfigurables baseOpts;
+  o2::pwglf::mccentrality::BuilderModule mcCentralityBuilderModule;
 
-  o2::pwglf::ParticleCounter<o2::framework::O2DatabasePDG> mCounter;
+  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
-  void init(o2::framework::InitContext& /*initContext*/)
+  void init(o2::framework::InitContext& initContext)
   {
     // Set up the CCDB
     ccdb->setURL(ccdbUrl.value);
     ccdb->setCaching(true);
     ccdb->setLocalObjectValidityChecking();
-    ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
     ccdb->setFatalWhenNull(false);
 
-    mCounter.mPdgDatabase = pdgDB.service;
-    mCounter.mSelectPrimaries = selectPrimaries.value;
-
-    if (fillFt0M) {
-      histos.add("FT0M/percentile", "FT0M percentile.", HistType::kTH1D, {{binsPercentile, "FT0M percentile"}});
-      histos.add("FT0M/percentilevsMult", "FT0M percentile.", HistType::kTH2D, {{binsPercentile, "FT0M percentile"}, {binsMultiplicity, "FT0M mult."}});
-    }
-    if (fillFt0A) {
-      histos.add("FT0A/percentile", "FT0A percentile.", HistType::kTH1D, {{binsPercentile, "FT0A percentile"}});
-      histos.add("FT0A/percentilevsMult", "FT0A percentile.", HistType::kTH2D, {{binsPercentile, "FT0A percentile"}, {binsMultiplicity, "FT0A mult."}});
-    }
-    if (fillFt0C) {
-      histos.add("FT0C/percentile", "FT0C percentile.", HistType::kTH1D, {{binsPercentile, "FT0C percentile"}});
-      histos.add("FT0C/percentilevsMult", "FT0C percentile.", HistType::kTH2D, {{binsPercentile, "FT0C percentile"}, {binsMultiplicity, "FT0C mult."}});
-    }
-
-    TList* lOfInput = nullptr;
-    if (path.value.rfind("ccdb://", 0) == 0) { // Getting post calib. from CCDB
-      path.value.replace(0, 7, "");
-      lOfInput = ccdb->get<TList>(path);
-      if (!lOfInput) {
-        if (doNotCrashOnNull) {
-          LOG(info) << "Could not find the calibration TList from CCDB in path " << path << ", will fill tables with dummy values";
-        } else {
-          LOG(fatal) << "Could not find the calibration TList from CCDB in path " << path;
-          return;
-        }
-      }
-    } else { // Getting post calib. from file
-      TFile* f = TFile::Open(path.value.c_str(), "READ");
-      if (!f) {
-        LOG(fatal) << "The input file " << path << " is not valid";
-      }
-      if (!f->IsOpen()) {
-        LOG(fatal) << "The input file " << f->GetName() << " is not open";
-      }
-      lOfInput = static_cast<TList*>(f->Get("ccdb_object"));
-      if (!lOfInput) {
-        f->ls();
-        LOG(fatal) << "The input file " << path.value << " does not contain the TList ccdb_object";
-      }
-    }
-    auto getHist = [this, lOfInput](const char* name) -> TH1F* {
-      if (!lOfInput) {
-        return nullptr;
-      }
-      auto hist = static_cast<TH1F*>(lOfInput->FindObject(name));
-      if (!hist) {
-        lOfInput->ls();
-        if (this->doNotCrashOnNull) {
-          LOG(info) << "Could not open histogram " << name << " from TList, will fill tables with dummy values";
-        } else {
-          LOG(fatal) << "Could not open histogram " << name << " from TList";
-        }
-      }
-      return hist;
-    };
-    if (fillFt0M) {
-      h1dFT0M = getHist("h1dFT0M");
-    }
-    if (fillFt0A) {
-      h1dFT0A = getHist("h1dFT0A");
-    }
-    if (fillFt0C) {
-      h1dFT0C = getHist("h1dFT0C");
-    }
+    mcCentralityBuilderModule.init(baseOpts, histos, initContext);
   }
 
   // Full tables (independent on central calibrations)
-  void process(aod::McCollision const& /*mcCollision*/,
-               aod::McParticles const& mcParticles)
+  void process(soa::Join<aod::McCollisions, aod::MultMCExtras> const& mcCollisions,
+               aod::BCsWithTimestamps const& bcs)
   {
-    const float nFT0A = mCounter.countFT0A(mcParticles);
-    const float nFT0C = mCounter.countFT0C(mcParticles);
-    const float nFT0M = nFT0A + nFT0C;
-    // const float nFV0A = mCounter.countFV0A(mcParticles);
-
-    if (fillFt0M) {
-      const float valueCentFT0M = h1dFT0M ? h1dFT0M->GetBinContent(h1dFT0M->FindBin(nFT0M)) : 105.0f;
-      centFT0M(valueCentFT0M);
-      histos.fill(HIST("FT0M/percentile"), valueCentFT0M);
-      histos.fill(HIST("FT0M/percentilevsMult"), valueCentFT0M, nFT0M);
-    }
-    if (fillFt0A) {
-      const float valueCentFT0A = h1dFT0A ? h1dFT0A->GetBinContent(h1dFT0A->FindBin(nFT0A)) : 105.0f;
-      centFT0A(valueCentFT0A);
-      histos.fill(HIST("FT0A/percentile"), valueCentFT0A);
-      histos.fill(HIST("FT0A/percentilevsMult"), valueCentFT0A, nFT0A);
-    }
-    if (fillFt0C) {
-      const float valueCentFT0C = h1dFT0C ? h1dFT0C->GetBinContent(h1dFT0C->FindBin(nFT0C)) : 105.0f;
-      centFT0C(valueCentFT0C);
-      histos.fill(HIST("FT0C/percentile"), valueCentFT0C);
-      histos.fill(HIST("FT0C/percentilevsMult"), valueCentFT0C, nFT0C);
-    }
-    // const float valueCentFV0A = h1dFT0M->GetBinContent(h1dFT0M->FindBin(nFV0A));
-    // centFV0A(valueCentFV0A);
+    mcCentralityBuilderModule.dataProcess(ccdb, histos, bcs, mcCollisions, products);
   }
 };
 
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<McCentrality>(cfgc)}; }
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<McCentrality>(cfgc)};
+}
diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt
index 83b05e9a1bc..4ec8dc02d4c 100644
--- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt
+++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt
@@ -188,3 +188,8 @@ o2physics_add_dpl_workflow(phi-strange-correlator
     SOURCES phiStrangeCorrelator.cxx
     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(strangeness-mc-centrality
+    SOURCES strangenessMcCentrality.cxx
+    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+    COMPONENT_NAME Analysis)
diff --git a/PWGLF/TableProducer/Strangeness/Converters/stramccollmultconverter2.cxx b/PWGLF/TableProducer/Strangeness/Converters/stramccollmultconverter2.cxx
new file mode 100644
index 00000000000..cc56470aacd
--- /dev/null
+++ b/PWGLF/TableProducer/Strangeness/Converters/stramccollmultconverter2.cxx
@@ -0,0 +1,46 @@
+// 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.
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/runDataProcessing.h>
+
+using namespace o2;
+using namespace o2::framework;
+
+// Converts V0 version 001 to 002
+struct stramccollmultconverter2 {
+  Produces<aod::StraMCCollMults_002> straMCCollMults_002;
+
+  void process(aod::StraMCCollMults_001 const& straMCcolls)
+  {
+    straMCCollMults_002.reserve(straMCcolls.size());
+    for (auto& straMCcoll : straMCcolls) {
+      straMCCollMults_002(straMCcoll.multMCFT0A(),
+                          straMCcoll.multMCFT0C(),
+                          -1., // dummy value multMCFV0A
+                          -1., // dummy value multMCFDDA
+                          -1., // dummy value multMCFDDC
+                          straMCcoll.multMCNParticlesEta05(),
+                          straMCcoll.multMCNParticlesEta08(),
+                          straMCcoll.multMCNParticlesEta10(),
+                          straMCcoll.totalMultMCParticles());
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<stramccollmultconverter2>(cfgc)};
+}
diff --git a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx
index 94cf60e61d6..0c07172d25d 100644
--- a/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx
+++ b/PWGLF/TableProducer/Strangeness/strangederivedbuilder.cxx
@@ -264,7 +264,7 @@ struct strangederivedbuilder {
   std::vector<uint32_t> genOmegaPlus;
 
   // create collision indices beforehand
-  std::vector<uint64_t> TrackCollIndices; // index -1: no collision
+  std::vector<uint64_t> TrackGlobalBc;
   std::vector<int> V0CollIndices;         // index -1: no collision
   std::vector<int> CascadeCollIndices;    // index -1: no collision
   std::vector<int> KFCascadeCollIndices;  // index -1: no collision
@@ -280,7 +280,7 @@ struct strangederivedbuilder {
   }
 
   template <typename VType>
-  void getCfg(o2::framework::InitContext& initContext, const std::string& name, VType& v, const std::string& task)
+  void getCfg(o2::framework::InitContext& initContext, const std::string& name, VType& v, const std::string task)
   {
     if (!o2::common::core::getTaskOptionValue(initContext, task, name, v, inheritEvtSelFromTaskVerbose)) {
       LOG(fatal) << "Could not get " << name << " from " << task << " task";
@@ -288,7 +288,7 @@ struct strangederivedbuilder {
   }
 
   template <typename TCollision>
-  bool isCollisionAccepted(TCollision collision, std::array<int, o2::aod::straselections::kNsel>& nSelected)
+  bool isCollisionAccepted(TCollision const& collision, std::array<int, o2::aod::straselections::kNsel>& nSelected)
   // check whether the collision passes our collision selections
   {
     if (requireTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) {
@@ -552,13 +552,13 @@ struct strangederivedbuilder {
   void populateCollisionTables(coll const& collisions, udcoll const& udCollisions, tracks const& Tracks, v0d const& V0s, cad const& Cascades, kfcad const& KFCascades, tracad const& TraCascades, bcType const& /*bcs*/)
   {
     // create collision indices beforehand
-    TrackCollIndices.clear();
+    TrackGlobalBc.clear();
     V0CollIndices.clear();
     CascadeCollIndices.clear();
     KFCascadeCollIndices.clear();
     TraCascadeCollIndices.clear();
 
-    TrackCollIndices.resize(Tracks.size(), 0);            // index -1: no collision
+    TrackGlobalBc.resize(Tracks.size(), 0);
     V0CollIndices.resize(V0s.size(), -1);                 // index -1: no collision
     CascadeCollIndices.resize(Cascades.size(), -1);       // index -1: no collision
     KFCascadeCollIndices.resize(KFCascades.size(), -1);   // index -1: no collision
@@ -645,7 +645,7 @@ struct strangederivedbuilder {
       }
 
       for (const auto& track : TrackTable_thisColl)
-        TrackCollIndices[track.globalIndex()] = bc.globalBC();
+        TrackGlobalBc[track.globalIndex()] = bc.globalBC();
 
       if (fillOnlySelectedCollisions && !isCollisionAccepted(collision, totalNbrCollisionsPerSelection)) {
         continue;
@@ -854,6 +854,7 @@ struct strangederivedbuilder {
       products.strangeMCColl(mccollision.posX(), mccollision.posY(), mccollision.posZ(),
                              mccollision.impactParameter(), mccollision.eventPlaneAngle(), mccollision.generatorsID());
       products.strangeMCMults(mccollision.multMCFT0A(), mccollision.multMCFT0C(),
+                              mccollision.multMCFV0A(), mccollision.multMCFDDA(), mccollision.multMCFDDC(),
                               mccollision.multMCNParticlesEta05(),
                               mccollision.multMCNParticlesEta08(),
                               mccollision.multMCNParticlesEta10(),
@@ -1047,7 +1048,7 @@ struct strangederivedbuilder {
                                    aod::dautrack::packing::packInInt8(tr.tpcNSigmaPr()));
           // populate daughter-level TOF information
           if (tr.hasTOF()) {
-            products.dauTrackTOFPIDs(TrackCollIndices[tr.globalIndex()], products.dauTrackExtras.lastIndex(), tr.tofSignal(), tr.tofEvTime(), tr.tofEvTimeErr(), tr.length(), tr.tofExpMom());
+            products.dauTrackTOFPIDs(TrackGlobalBc[tr.globalIndex()], products.dauTrackExtras.lastIndex(), tr.tofSignal(), tr.tofEvTime(), tr.tofEvTimeErr(), tr.length(), tr.tofExpMom());
           }
         } else {
           // populate with empty fully-compatible Nsigmas if no corresponding table available
diff --git a/PWGLF/TableProducer/Strangeness/strangenessMcCentrality.cxx b/PWGLF/TableProducer/Strangeness/strangenessMcCentrality.cxx
new file mode 100644
index 00000000000..31c969b5303
--- /dev/null
+++ b/PWGLF/TableProducer/Strangeness/strangenessMcCentrality.cxx
@@ -0,0 +1,71 @@
+// 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   mcCentrality.cxx
+/// \author Romain Schotter romain.schotter@cern.ch
+/// \brief  Task to produce the MC centrality table for strangeness derived data
+///
+
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+#include "PWGLF/Utils/mcCentralityModule.h"
+
+#include <CCDB/BasicCCDBManager.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/Configurable.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/InitContext.h>
+#include <Framework/OutputObjHeader.h>
+#include <Framework/runDataProcessing.h>
+
+#include <string>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+/// Task to produce the response table
+struct StrangenessMcCentrality {
+  // Input parameters
+  o2::framework::Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+  Service<o2::ccdb::BasicCCDBManager> ccdb;
+
+  o2::pwglf::mccentrality::products products;
+  o2::pwglf::mccentrality::coreConfigurables baseOpts;
+  o2::pwglf::mccentrality::BuilderModule mcCentralityBuilderModule;
+
+  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+  void init(o2::framework::InitContext& initContext)
+  {
+    // Set up the CCDB
+    ccdb->setURL(ccdbUrl.value);
+    ccdb->setCaching(true);
+    ccdb->setLocalObjectValidityChecking();
+    ccdb->setFatalWhenNull(false);
+
+    mcCentralityBuilderModule.init(baseOpts, histos, initContext);
+  }
+
+  // Full tables (independent on central calibrations)
+  void process(aod::StraMCCollMults const& mcCollisions,
+               aod::StraStamps const& bcs)
+  {
+    mcCentralityBuilderModule.dataProcess(ccdb, histos, bcs, mcCollisions, products);
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<StrangenessMcCentrality>(cfgc)};
+}
diff --git a/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx b/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx
index 2a10605f42d..aa35412c7f4 100644
--- a/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx
+++ b/PWGLF/Tasks/Strangeness/derivedlambdakzeroanalysis.cxx
@@ -31,6 +31,7 @@
 #include "PWGLF/DataModel/LFStrangenessMLTables.h"
 #include "PWGLF/DataModel/LFStrangenessPIDTables.h"
 #include "PWGLF/DataModel/LFStrangenessTables.h"
+#include "PWGLF/DataModel/mcCentrality.h"
 #include "PWGUD/Core/SGSelector.h"
 
 #include "Common/CCDB/EventSelectionParams.h"
@@ -102,7 +103,7 @@ enum CentEstimator {
 struct derivedlambdakzeroanalysis {
   HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
 
-  bool isRun3;
+  bool isRun3 = true;
 
   // master analysis switches
   Configurable<bool> analyseK0Short{"analyseK0Short", true, "process K0Short-like candidates"};
@@ -113,6 +114,7 @@ struct derivedlambdakzeroanalysis {
   Configurable<std::string> irSource{"irSource", "", "Estimator of the interaction rate (Recommended: pp --> T0VTX, Pb-Pb --> ZNC hadronic)"};
   Configurable<int> centralityEstimator{"centralityEstimator", kCentFT0C, "Run 3 centrality estimator (0:CentFT0C, 1:CentFT0M, 2:CentFT0CVariant1, 3:CentMFT, 4:CentNGlobal, 5:CentFV0A)"};
   Configurable<bool> doUPCanalysis{"doUPCanalysis", true, "Study V0s in hadronic and UPC collisions"};
+  Configurable<bool> useMcCentrality{"useMcCentrality", false, "Use MC centrality?"};
 
   Configurable<bool> doEventQA{"doEventQA", false, "do event QA histograms"};
   Configurable<bool> doCompleteTopoQA{"doCompleteTopoQA", false, "do topological variable QA histograms"};
@@ -311,8 +313,8 @@ struct derivedlambdakzeroanalysis {
   o2::ccdb::CcdbApi ccdbApi;
   Service<o2::ccdb::BasicCCDBManager> ccdb;
   ctpRateFetcher rateFetcher;
-  int mRunNumber;
-  float magField;
+  int mRunNumber = 0.;
+  float magField = 0.;
   std::map<std::string, std::string> metadata;
   o2::parameters::GRPMagField* grpmag = nullptr;
 
@@ -437,24 +439,24 @@ struct derivedlambdakzeroanalysis {
                               selPhysPrimAntiLambda, // for mc tagging
   };
 
-  uint64_t maskTopological;
-  uint64_t maskTopoNoV0Radius;
-  uint64_t maskTopoNoDCANegToPV;
-  uint64_t maskTopoNoDCAPosToPV;
-  uint64_t maskTopoNoCosPA;
-  uint64_t maskTopoNoDCAV0Dau;
-  uint64_t maskTrackProperties;
+  uint64_t maskTopological = 0;
+  uint64_t maskTopoNoV0Radius = 0;
+  uint64_t maskTopoNoDCANegToPV = 0;
+  uint64_t maskTopoNoDCAPosToPV = 0;
+  uint64_t maskTopoNoCosPA = 0;
+  uint64_t maskTopoNoDCAV0Dau = 0;
+  uint64_t maskTrackProperties = 0;
 
-  uint64_t maskK0ShortSpecific;
-  uint64_t maskLambdaSpecific;
-  uint64_t maskAntiLambdaSpecific;
+  uint64_t maskK0ShortSpecific = 0;
+  uint64_t maskLambdaSpecific = 0;
+  uint64_t maskAntiLambdaSpecific = 0;
 
-  uint64_t maskSelectionK0Short;
-  uint64_t maskSelectionLambda;
-  uint64_t maskSelectionAntiLambda;
+  uint64_t maskSelectionK0Short = 0;
+  uint64_t maskSelectionLambda = 0;
+  uint64_t maskSelectionAntiLambda = 0;
 
-  uint64_t secondaryMaskSelectionLambda;
-  uint64_t secondaryMaskSelectionAntiLambda;
+  uint64_t secondaryMaskSelectionLambda = 0;
+  uint64_t secondaryMaskSelectionAntiLambda = 0;
 
   void init(InitContext const&)
   {
@@ -733,6 +735,9 @@ struct derivedlambdakzeroanalysis {
 
     histos.add("hEventCentrality", "hEventCentrality", kTH1D, {axisConfigurations.axisCentralityFine});
     histos.add("hCentralityVsNch", "hCentralityVsNch", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisNch});
+    if (useMcCentrality) {
+      histos.add("hEventVsGenCentrality", "hEventVsGenCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisCentralityFine});
+    }
     if (doEventQA) {
       if (isRun3) {
         histos.add("hEventSelectionVsCentrality", "hEventSelectionVsCentrality", kTH2D, {{23, -0.5f, +22.5f}, axisConfigurations.axisCentralityFine});
@@ -1190,10 +1195,16 @@ struct derivedlambdakzeroanalysis {
 
     // Creation of histograms: MC generated
     if ((doprocessGeneratedRun3 || doprocessGeneratedRun2)) {
-      histos.add("hGenEvents", "hGenEvents", kTH2D, {{axisConfigurations.axisNch}, {2, -0.5f, +1.5f}});
-      histos.get<TH2>(HIST("hGenEvents"))->GetYaxis()->SetBinLabel(1, "All gen. events");
-      histos.get<TH2>(HIST("hGenEvents"))->GetYaxis()->SetBinLabel(2, "Gen. with at least 1 rec. events");
-      histos.add("hGenEventCentrality", "hGenEventCentrality", kTH1D, {{101, 0.0f, 101.0f}});
+      if (useMcCentrality) {
+        histos.add("hGenEventsVsMcCentrality", "hGenEventsVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, {2, -0.5f, +1.5f}});
+        histos.get<TH2>(HIST("hGenEventsVsMcCentrality"))->GetYaxis()->SetBinLabel(1, "All gen. events");
+        histos.get<TH2>(HIST("hGenEventsVsMcCentrality"))->GetYaxis()->SetBinLabel(2, "Gen. with at least 1 rec. events");
+      } else {
+        histos.add("hGenEvents", "hGenEvents", kTH2D, {{axisConfigurations.axisNch}, {2, -0.5f, +1.5f}});
+        histos.get<TH2>(HIST("hGenEvents"))->GetYaxis()->SetBinLabel(1, "All gen. events");
+        histos.get<TH2>(HIST("hGenEvents"))->GetYaxis()->SetBinLabel(2, "Gen. with at least 1 rec. events");
+      }
+      histos.add("hGenEventCentrality", "hGenEventCentrality", kTH1D, {axisConfigurations.axisCentralityFine});
 
       histos.add("hCentralityVsNcoll_beforeEvSel", "hCentralityVsNcoll_beforeEvSel", kTH2D, {axisConfigurations.axisCentrality, {50, -0.5f, 49.5f}});
       histos.add("hCentralityVsNcoll_afterEvSel", "hCentralityVsNcoll_afterEvSel", kTH2D, {axisConfigurations.axisCentrality, {50, -0.5f, 49.5f}});
@@ -1208,21 +1219,39 @@ struct derivedlambdakzeroanalysis {
       histos.add("h2dGenOmegaMinus", "h2dGenOmegaMinus", kTH2D, {axisConfigurations.axisCentrality, axisConfigurations.axisPt});
       histos.add("h2dGenOmegaPlus", "h2dGenOmegaPlus", kTH2D, {axisConfigurations.axisCentrality, axisConfigurations.axisPt});
 
-      histos.add("h2dGenK0ShortVsMultMC_RecoedEvt", "h2dGenK0ShortVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenLambdaVsMultMC_RecoedEvt", "h2dGenLambdaVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenAntiLambdaVsMultMC_RecoedEvt", "h2dGenAntiLambdaVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenXiMinusVsMultMC_RecoedEvt", "h2dGenXiMinusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenXiPlusVsMultMC_RecoedEvt", "h2dGenXiPlusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenOmegaMinusVsMultMC_RecoedEvt", "h2dGenOmegaMinusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenOmegaPlusVsMultMC_RecoedEvt", "h2dGenOmegaPlusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-
-      histos.add("h2dGenK0ShortVsMultMC", "h2dGenK0ShortVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenLambdaVsMultMC", "h2dGenLambdaVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenAntiLambdaVsMultMC", "h2dGenAntiLambdaVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenXiMinusVsMultMC", "h2dGenXiMinusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenXiPlusVsMultMC", "h2dGenXiPlusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenOmegaMinusVsMultMC", "h2dGenOmegaMinusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
-      histos.add("h2dGenOmegaPlusVsMultMC", "h2dGenOmegaPlusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+      if (useMcCentrality) {
+        histos.add("h2dGenK0ShortVsMcCentrality_RecoedEvt", "h2dGenK0ShortVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenLambdaVsMcCentrality_RecoedEvt", "h2dGenLambdaVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenAntiLambdaVsMcCentrality_RecoedEvt", "h2dGenAntiLambdaVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenXiMinusVsMcCentrality_RecoedEvt", "h2dGenXiMinusVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenXiPlusVsMcCentrality_RecoedEvt", "h2dGenXiPlusVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaMinusVsMcCentrality_RecoedEvt", "h2dGenOmegaMinusVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaPlusVsMcCentrality_RecoedEvt", "h2dGenOmegaPlusVsMcCentrality_RecoedEvt", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+
+        histos.add("h2dGenK0ShortVsMcCentrality", "h2dGenK0ShortVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenLambdaVsMcCentrality", "h2dGenLambdaVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenAntiLambdaVsMcCentrality", "h2dGenAntiLambdaVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenXiMinusVsMcCentrality", "h2dGenXiMinusVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenXiPlusVsMcCentrality", "h2dGenXiPlusVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaMinusVsMcCentrality", "h2dGenOmegaMinusVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaPlusVsMcCentrality", "h2dGenOmegaPlusVsMcCentrality", kTH2D, {axisConfigurations.axisCentralityFine, axisConfigurations.axisPt});
+      } else {
+        histos.add("h2dGenK0ShortVsMultMC_RecoedEvt", "h2dGenK0ShortVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenLambdaVsMultMC_RecoedEvt", "h2dGenLambdaVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenAntiLambdaVsMultMC_RecoedEvt", "h2dGenAntiLambdaVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenXiMinusVsMultMC_RecoedEvt", "h2dGenXiMinusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenXiPlusVsMultMC_RecoedEvt", "h2dGenXiPlusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaMinusVsMultMC_RecoedEvt", "h2dGenOmegaMinusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaPlusVsMultMC_RecoedEvt", "h2dGenOmegaPlusVsMultMC_RecoedEvt", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+
+        histos.add("h2dGenK0ShortVsMultMC", "h2dGenK0ShortVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenLambdaVsMultMC", "h2dGenLambdaVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenAntiLambdaVsMultMC", "h2dGenAntiLambdaVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenXiMinusVsMultMC", "h2dGenXiMinusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenXiPlusVsMultMC", "h2dGenXiPlusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaMinusVsMultMC", "h2dGenOmegaMinusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+        histos.add("h2dGenOmegaPlusVsMultMC", "h2dGenOmegaPlusVsMultMC", kTH2D, {axisConfigurations.axisNch, axisConfigurations.axisPt});
+      }
 
       if (doSecondaryV0s) {
         histos.add("h2dGenSecLambda", "h2dGenSecLambda", kTH2D, {axisConfigurations.axisCentrality, axisConfigurations.axisPt});
@@ -1266,22 +1295,42 @@ struct derivedlambdakzeroanalysis {
   // If takeMcCentrality is enabled, the centrality is taken from the MC collision; otherwise it is taken
   // from the reconstructed collision. Returns -1 if no corresponding centrality estimator is found or if no MC collision is associated to the recoed collision.
   template <typename TCollision>
-  auto getCentralityRun3(TCollision const& collision)
+  auto getCentralityRun3(TCollision const& collision, bool takeMcCentrality = false)
   {
-    if (centralityEstimator == kCentFT0C)
-      return collision.centFT0C();
-    else if (centralityEstimator == kCentFT0M)
-      return collision.centFT0M();
-    else if (centralityEstimator == kCentFT0CVariant1)
-      return collision.centFT0CVariant1();
-    else if (centralityEstimator == kCentMFT)
-      return collision.centMFT();
-    else if (centralityEstimator == kCentNGlobal)
-      return collision.centNGlobal();
-    else if (centralityEstimator == kCentFV0A)
-      return collision.centFV0A();
-
-    return -1.f;
+    // Helper lambda to extract centrality from any object exposing the cent* columns
+    auto extractCentrality = [this](auto const& coll) -> float {
+      switch (centralityEstimator) {
+        case kCentFT0C:
+          return coll.centFT0C();
+        case kCentFT0M:
+          return coll.centFT0M();
+        case kCentFT0CVariant1:
+          return coll.centFT0CVariant1();
+        // case kCentMFT:          return coll.centMFT();
+        case kCentNGlobal:
+          return coll.centNGlobal();
+        case kCentFV0A:
+          return coll.centFV0A();
+        default:
+          return -1.f;
+      }
+    };
+
+    // check if we are in MC
+    if constexpr (requires { collision.straMCCollisionId(); }) {
+      if (takeMcCentrality) { // check if we need to extract centrality from recoed or generated collision
+        // if recoed collision has a MC collision, de-reference it and extract centrality from it
+        // otherwise return -1.f
+        if (collision.has_straMCCollision()) {
+          auto mcCollision = collision.template straMCCollision_as<soa::Join<aod::StraMCCollisions, aod::StraMCCollMults, aod::McCentFV0As, aod::McCentFT0Ms, aod::McCentFT0Cs, aod::McCentFT0CVariant1s, aod::McCentNGlobals>>();
+          return extractCentrality(mcCollision);
+        } else {
+          return -1.f;
+        }
+      }
+    }
+
+    return extractCentrality(collision);
   }
 
   // ______________________________________________________
@@ -1297,7 +1346,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TCollision>
-  void initCCDB(TCollision collision)
+  void initCCDB(TCollision const& collision)
   {
     if (mRunNumber == collision.runNumber()) {
       return;
@@ -1358,7 +1407,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TV0, typename TCollision>
-  uint64_t computeReconstructionBitmap(TV0 v0, TCollision collision, float rapidityLambda, float rapidityK0Short, float /*pT*/)
+  uint64_t computeReconstructionBitmap(TV0 const& v0, TCollision const& collision, float rapidityLambda, float rapidityK0Short, float /*pT*/)
   // precalculate this information so that a check is one mask operation, not many
   {
     uint64_t bitMap = 0;
@@ -1621,7 +1670,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TV0>
-  uint64_t computeMCAssociation(TV0 v0)
+  uint64_t computeMCAssociation(TV0 const& v0)
   // precalculate this information so that a check is one mask operation, not many
   {
     uint64_t bitMap = 0;
@@ -1804,7 +1853,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TV0>
-  void analyseCandidate(TV0 v0, float pt, float centrality, uint64_t selMap, uint8_t gapSide, int& nK0Shorts, int& nLambdas, int& nAntiLambdas)
+  void analyseCandidate(TV0 const& v0, float pt, float centrality, uint64_t selMap, uint8_t gapSide, int& nK0Shorts, int& nLambdas, int& nAntiLambdas)
   // precalculate this information so that a check is one mask operation, not many
   {
     bool passK0ShortSelections = false;
@@ -2207,7 +2256,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TV0>
-  void analyseCollisionAssociation(TV0 /*v0*/, float pt, int mcNch, bool correctAssociation, uint64_t selMap)
+  void analyseCollisionAssociation(TV0 const& /*v0*/, float pt, int mcNch, bool correctAssociation, uint64_t selMap)
   // analyse collision association
   {
     // __________________________________________
@@ -2230,7 +2279,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TV0>
-  void fillFeeddownMatrix(TV0 v0, float pt, float centrality, uint64_t selMap)
+  void fillFeeddownMatrix(TV0 const& v0, float pt, float centrality, uint64_t selMap)
   // fill feeddown matrix for Lambdas or AntiLambdas
   // fixme: a potential improvement would be to consider mass windows for the l/al
   {
@@ -2284,7 +2333,7 @@ struct derivedlambdakzeroanalysis {
   }
 
   template <typename TCollision>
-  bool isEventAccepted(TCollision collision, bool fillHists)
+  bool isEventAccepted(TCollision const& collision, bool fillHists)
   // check whether the collision passes our collision selections
   {
     float centrality = -1.0f;
@@ -2292,7 +2341,7 @@ struct derivedlambdakzeroanalysis {
       histos.fill(HIST("hEventSelection"), 0. /* all collisions */);
       if (doEventQA) {
         if constexpr (requires { collision.centFT0C(); }) { // check if we are in Run 3
-          centrality = getCentralityRun3(collision);
+          centrality = getCentralityRun3(collision, useMcCentrality);
         }
         histos.fill(HIST("hEventSelectionVsCentrality"), 0. /* all collisions */, centrality);
       }
@@ -2672,7 +2721,7 @@ struct derivedlambdakzeroanalysis {
   void fillReconstructedEventProperties(TCollision const& collision, float& centrality, float& collisionOccupancy, double& interactionRate, int& gapSide, int& selGapSide)
   {
     if constexpr (requires { collision.centFT0C(); }) { // check if we are in Run 3
-      centrality = getCentralityRun3(collision);
+      centrality = getCentralityRun3(collision, useMcCentrality);
       collisionOccupancy = eventSelections.useFT0CbasedOccupancy ? collision.ft0cOccupancyInTimeRange() : collision.trackOccupancyInTimeRange();
       // Fetch interaction rate only if required (in order to limit ccdb calls)
       interactionRate = !irSource.value.empty() ? rateFetcher.fetch(ccdb.service, collision.timestamp(), collision.runNumber(), irSource) * 1.e-3 : -1;
@@ -2682,6 +2731,10 @@ struct derivedlambdakzeroanalysis {
         centrality = hRawCentrality->GetBinContent(hRawCentrality->FindBin(centralityEstimator == kCentFT0M ? collision.multFT0A() + collision.multFT0C() : collision.multFT0C()));
       }
 
+      if (useMcCentrality) {
+        histos.fill(HIST("hEventVsGenCentrality"), getCentralityRun3(collision, false), getCentralityRun3(collision, true));
+      }
+
       if (doUPCanalysis) {
         // gap side
         gapSide = collision.gapSide();
@@ -2754,7 +2807,15 @@ struct derivedlambdakzeroanalysis {
         continue;
       }
 
-      histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/);
+      if constexpr (run3) {
+        if (useMcCentrality) {
+          histos.fill(HIST("hGenEventsVsMcCentrality"), getCentralityRun3(mcCollision), 0 /* all gen. events*/);
+        } else {
+          histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/);
+        }
+      } else {
+        histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/);
+      }
 
       auto groupedCollisions = getGroupedCollisions<run3>(collisions, mcCollision.globalIndex());
       // Check if there is at least one of the reconstructed collisions associated to this MC collision
@@ -2772,7 +2833,7 @@ struct derivedlambdakzeroanalysis {
         if constexpr (run3) { // check if we are in Run 3
           if (biggestNContribs < collision.multPVTotalContributors()) {
             biggestNContribs = collision.multPVTotalContributors();
-            centrality = getCentralityRun3(collision);
+            centrality = getCentralityRun3(collision, useMcCentrality);
           }
         } else { // we are in Run 2: there should be only one collision in groupedCollisions
           centrality = eventSelections.useSPDTrackletsCent ? collision.centRun2SPDTracklets() : collision.centRun2V0M();
@@ -2790,7 +2851,15 @@ struct derivedlambdakzeroanalysis {
       histos.fill(HIST("hEventPVzMC"), mcCollision.posZ());
 
       if (atLeastOne) {
-        histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 1 /* at least 1 rec. event*/);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("hGenEventsVsMcCentrality"), getCentralityRun3(mcCollision), 1 /* at least 1 rec. event*/);
+          } else {
+            histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 1 /* at least 1 rec. event*/);
+          }
+        } else {
+          histos.fill(HIST("hGenEvents"), mcCollision.multMCNParticlesEta05(), 1 /* at least 1 rec. event*/);
+        }
 
         histos.fill(HIST("hGenEventCentrality"), centrality);
       }
@@ -2869,8 +2938,8 @@ struct derivedlambdakzeroanalysis {
 
   // ______________________________________________________
   // Simulated processing (subscribes to MC information too)
-  template <typename TCollision, typename TV0s>
-  void analyzeRecoedV0sInMonteCarlo(TCollision const& collision, TV0s const& fullV0s)
+  template <typename TCollision, typename TV0s, typename TMCCollisions>
+  void analyzeRecoedV0sInMonteCarlo(TCollision const& collision, TV0s const& fullV0s, TMCCollisions const& /*mcCollisions*/)
   {
     // Fire up CCDB
     if ((mlConfigurations.useK0ShortScores && mlConfigurations.calculateK0ShortScores) ||
@@ -2947,7 +3016,7 @@ struct derivedlambdakzeroanalysis {
         bool correctCollision = false;
         int mcNch = -1;
         if (collision.has_straMCCollision()) {
-          auto mcCollision = collision.template straMCCollision_as<soa::Join<aod::StraMCCollisions, aod::StraMCCollMults>>();
+          auto mcCollision = collision.template straMCCollision_as<TMCCollisions>();
           mcNch = mcCollision.multMCNParticlesEta05();
           correctCollision = (v0MC.straMCCollisionId() == mcCollision.globalIndex());
         }
@@ -2993,7 +3062,7 @@ struct derivedlambdakzeroanalysis {
           ymc > v0Selections.rapidityMaxCut)
         continue;
 
-      auto mcCollision = v0MC.template straMCCollision_as<soa::Join<aod::StraMCCollisions, aod::StraMCCollMults>>();
+      auto mcCollision = v0MC.template straMCCollision_as<TMCCollisions>();
       if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) {
         continue;
       }
@@ -3009,33 +3078,81 @@ struct derivedlambdakzeroanalysis {
       if (listBestCollisionIdx[mcCollision.globalIndex()] > -1) {
         auto collision = collisions.iteratorAt(listBestCollisionIdx[mcCollision.globalIndex()]);
         if constexpr (requires { collision.centFT0C(); }) { // check if we are in Run 3
-          centrality = getCentralityRun3(collision);
+          centrality = getCentralityRun3(collision, useMcCentrality);
         } else { // no, we are in Run 2
           centrality = eventSelections.useSPDTrackletsCent ? collision.centRun2SPDTracklets() : collision.centRun2V0M();
         }
 
         if (v0MC.pdgCode() == PDG_t::kK0Short) {
-          histos.fill(HIST("h2dGenK0ShortVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenK0ShortVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenK0ShortVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenK0ShortVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
         if (v0MC.pdgCode() == PDG_t::kLambda0) {
-          histos.fill(HIST("h2dGenLambdaVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenLambdaVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenLambdaVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenLambdaVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
         if (v0MC.pdgCode() == PDG_t::kLambda0Bar) {
-          histos.fill(HIST("h2dGenAntiLambdaVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenAntiLambdaVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenAntiLambdaVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenAntiLambdaVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
       }
 
       if (v0MC.pdgCode() == PDG_t::kK0Short) {
         histos.fill(HIST("h2dGenK0Short"), centrality, ptmc);
-        histos.fill(HIST("h2dGenK0ShortVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenK0ShortVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenK0ShortVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenK0ShortVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
       if (v0MC.pdgCode() == PDG_t::kLambda0) {
         histos.fill(HIST("h2dGenLambda"), centrality, ptmc);
-        histos.fill(HIST("h2dGenLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenLambdaVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
       if (v0MC.pdgCode() == PDG_t::kLambda0Bar) {
         histos.fill(HIST("h2dGenAntiLambda"), centrality, ptmc);
-        histos.fill(HIST("h2dGenAntiLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenAntiLambdaVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenAntiLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenAntiLambdaVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
     }
 
@@ -3057,7 +3174,7 @@ struct derivedlambdakzeroanalysis {
           ymc > v0Selections.rapidityMaxCut)
         continue;
 
-      auto mcCollision = cascMC.template straMCCollision_as<soa::Join<aod::StraMCCollisions, aod::StraMCCollMults>>();
+      auto mcCollision = cascMC.template straMCCollision_as<TMCCollisions>();
       if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) {
         continue;
       }
@@ -3073,22 +3190,54 @@ struct derivedlambdakzeroanalysis {
       if (listBestCollisionIdx[mcCollision.globalIndex()] > -1) {
         auto collision = collisions.iteratorAt(listBestCollisionIdx[mcCollision.globalIndex()]);
         if constexpr (requires { collision.centFT0C(); }) { // check if we are in Run 3
-          centrality = getCentralityRun3(collision);
+          centrality = getCentralityRun3(collision, useMcCentrality);
         } else { // no, we are in Run 2
           centrality = eventSelections.useSPDTrackletsCent ? collision.centRun2SPDTracklets() : collision.centRun2V0M();
         }
 
         if (cascMC.pdgCode() == PDG_t::kXiMinus) {
-          histos.fill(HIST("h2dGenXiMinusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenXiMinusVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenXiMinusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenXiMinusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
         if (cascMC.pdgCode() == PDG_t::kXiPlusBar) {
-          histos.fill(HIST("h2dGenXiPlusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenXiPlusVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenXiPlusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenXiPlusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
         if (cascMC.pdgCode() == PDG_t::kOmegaMinus) {
-          histos.fill(HIST("h2dGenOmegaMinusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenOmegaMinusVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenOmegaMinusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenOmegaMinusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
         if (cascMC.pdgCode() == PDG_t::kOmegaPlusBar) {
-          histos.fill(HIST("h2dGenOmegaPlusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          if constexpr (run3) {
+            if (useMcCentrality) {
+              histos.fill(HIST("h2dGenOmegaPlusVsMcCentrality_RecoedEvt"), getCentralityRun3(mcCollision), ptmc);
+            } else {
+              histos.fill(HIST("h2dGenOmegaPlusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+            }
+          } else {
+            histos.fill(HIST("h2dGenOmegaPlusVsMultMC_RecoedEvt"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
         }
 
         if (doSecondaryV0s && std::abs(cascMC.pdgCodeV0()) == kLambda0) {
@@ -3116,19 +3265,51 @@ struct derivedlambdakzeroanalysis {
 
       if (cascMC.pdgCode() == PDG_t::kXiMinus) {
         histos.fill(HIST("h2dGenXiMinus"), centrality, ptmc);
-        histos.fill(HIST("h2dGenXiMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenXiMinusVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenXiMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenXiMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
       if (cascMC.pdgCode() == PDG_t::kXiPlusBar) {
         histos.fill(HIST("h2dGenXiPlus"), centrality, ptmc);
-        histos.fill(HIST("h2dGenXiPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenXiPlusVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenXiPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenXiPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
       if (cascMC.pdgCode() == PDG_t::kOmegaMinus) {
         histos.fill(HIST("h2dGenOmegaMinus"), centrality, ptmc);
-        histos.fill(HIST("h2dGenOmegaMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenOmegaMinusVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenOmegaMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenOmegaMinusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
       if (cascMC.pdgCode() == PDG_t::kOmegaPlusBar) {
         histos.fill(HIST("h2dGenOmegaPlus"), centrality, ptmc);
-        histos.fill(HIST("h2dGenOmegaPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        if constexpr (run3) {
+          if (useMcCentrality) {
+            histos.fill(HIST("h2dGenOmegaPlusVsMcCentrality"), getCentralityRun3(mcCollision), ptmc);
+          } else {
+            histos.fill(HIST("h2dGenOmegaPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+          }
+        } else {
+          histos.fill(HIST("h2dGenOmegaPlusVsMultMC"), mcCollision.multMCNParticlesEta05(), ptmc);
+        }
       }
 
       if (doSecondaryV0s && std::abs(cascMC.pdgCodeV0()) == kLambda0) {
@@ -3194,21 +3375,21 @@ struct derivedlambdakzeroanalysis {
 
   // ______________________________________________________
   // Simulated processing in Run 3 (subscribes to MC information too)
-  void processMonteCarloRun3(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraEvSelExtras, aod::StraStamps, aod::StraEvTimes, aod::StraCollLabels>::iterator const& collision, V0McCandidates const& fullV0s, DauTracks const&, aod::MotherMCParts const&, soa::Join<aod::StraMCCollisions, aod::StraMCCollMults> const& /*mccollisions*/, soa::Join<aod::V0MCCores, aod::V0MCCollRefs> const&)
+  void processMonteCarloRun3(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraEvSelExtras, aod::StraStamps, aod::StraEvTimes, aod::StraCollLabels>::iterator const& collision, V0McCandidates const& fullV0s, DauTracks const&, aod::MotherMCParts const&, soa::Join<aod::StraMCCollisions, aod::StraMCCollMults, aod::McCentFV0As, aod::McCentFT0Ms, aod::McCentFT0Cs, aod::McCentFT0CVariant1s, aod::McCentNGlobals> const& mccollisions, soa::Join<aod::V0MCCores, aod::V0MCCollRefs> const&)
   {
-    analyzeRecoedV0sInMonteCarlo(collision, fullV0s);
+    analyzeRecoedV0sInMonteCarlo(collision, fullV0s, mccollisions);
   }
 
   // ______________________________________________________
   // Simulated processing in Run 2 (subscribes to MC information too)
-  void processMonteCarloRun2(soa::Join<aod::StraCollisions, aod::StraCentsRun2, aod::StraEvSelsRun2, aod::StraStamps, aod::StraEvTimes, aod::StraCollLabels>::iterator const& collision, V0McCandidates const& fullV0s, DauTracks const&, aod::MotherMCParts const&, soa::Join<aod::StraMCCollisions, aod::StraMCCollMults> const& /*mccollisions*/, soa::Join<aod::V0MCCores, aod::V0MCCollRefs> const&)
+  void processMonteCarloRun2(soa::Join<aod::StraCollisions, aod::StraCentsRun2, aod::StraEvSelsRun2, aod::StraStamps, aod::StraEvTimes, aod::StraCollLabels>::iterator const& collision, V0McCandidates const& fullV0s, DauTracks const&, aod::MotherMCParts const&, soa::Join<aod::StraMCCollisions, aod::StraMCCollMults> const& mccollisions, soa::Join<aod::V0MCCores, aod::V0MCCollRefs> const&)
   {
-    analyzeRecoedV0sInMonteCarlo(collision, fullV0s);
+    analyzeRecoedV0sInMonteCarlo(collision, fullV0s, mccollisions);
   }
 
   // ______________________________________________________
   // Simulated processing in Run 3 (subscribes to MC information too)
-  void processGeneratedRun3(soa::Join<aod::StraMCCollisions, aod::StraMCCollMults> const& mcCollisions, soa::Join<aod::V0MCCores, aod::V0MCCollRefs> const& V0MCCores, soa::Join<aod::CascMCCores, aod::CascMCCollRefs> const& CascMCCores, soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraEvSelExtras, aod::StraStamps, aod::StraEvTimes, aod::StraCollLabels> const& collisions)
+  void processGeneratedRun3(soa::Join<aod::StraMCCollisions, aod::StraMCCollMults, aod::McCentFV0As, aod::McCentFT0Ms, aod::McCentFT0Cs, aod::McCentFT0CVariant1s, aod::McCentNGlobals> const& mcCollisions, soa::Join<aod::V0MCCores, aod::V0MCCollRefs> const& V0MCCores, soa::Join<aod::CascMCCores, aod::CascMCCollRefs> const& CascMCCores, soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraEvSelExtras, aod::StraStamps, aod::StraEvTimes, aod::StraCollLabels> const& collisions)
   {
     analyzeGeneratedV0s<true>(mcCollisions, V0MCCores, CascMCCores, collisions);
   }
diff --git a/PWGLF/Utils/mcCentralityModule.h b/PWGLF/Utils/mcCentralityModule.h
new file mode 100644
index 00000000000..cf2b9de0a3a
--- /dev/null
+++ b/PWGLF/Utils/mcCentralityModule.h
@@ -0,0 +1,737 @@
+// 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   mcCentralityModule.cxx
+/// \author Romain Schotter romain.schotter@cern.ch
+/// \brief  Module to produce MC centrality table based on ANY table containing MC multiplicity
+///
+
+#ifndef PWGLF_UTILS_MCCENTRALITYMODULE_H_
+#define PWGLF_UTILS_MCCENTRALITYMODULE_H_
+
+#include "PWGLF/DataModel/mcCentrality.h"
+
+#include <CCDB/BasicCCDBManager.h>
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/Array2D.h>
+#include <Framework/Configurable.h>
+#include <Framework/DataSpecUtils.h>
+#include <Framework/DeviceSpec.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/HistogramSpec.h>
+#include <Framework/RunningWorkflowInfo.h>
+#include <Framework/runDataProcessing.h>
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+#include <TRandom3.h>
+#include <TString.h>
+
+#include <RtypesCore.h>
+
+#include <concepts>
+#include <cstddef>
+#include <cstdint>
+#include <cstdio>
+#include <cstdlib>
+#include <memory>
+#include <string>
+#include <type_traits>
+#include <vector>
+
+//__________________________________________
+// strangeness builder module
+
+namespace o2
+{
+namespace pwglf
+{
+namespace mccentrality // avoid polluting other namespaces
+{
+
+// statics necessary for the configurables in this namespace
+static constexpr int nParameters = 1;
+static const std::vector<std::string> tableNames{
+  "McCentFV0As",
+  "McCentFT0Ms",
+  "McCentFT0As",
+  "McCentFT0Cs",
+  "McCentFT0CVariant1s",
+  "McCentFT0CVariant2s",
+  "McCentFDDMs",
+  "McCentNTPVs",
+  "McCentNGlobals",
+  "McCentMFTs"};
+
+static constexpr int nTablesConst = 10;
+static const std::vector<std::string> parameterNames{"enable"};
+static const int defaultParameters[nTablesConst][nParameters]{
+  {-1},
+  {-1},
+  {-1},
+  {-1},
+  {-1},
+  {-1},
+  {-1},
+  {-1},
+  {-1},
+  {-1}};
+
+// table index : match order above
+enum tableIndex { kFV0A = 0,
+                  kFT0M,
+                  kFT0A,
+                  kFT0C,
+                  kFT0CVariant1,
+                  kFT0CVariant2,
+                  kFDDM,
+                  kNTPV,
+                  kNGlobal,
+                  kMFT,
+                  kNestimators };
+
+static constexpr const char* DirList[] = {
+  "FV0A",
+  "FT0M",
+  "FT0A",
+  "FT0C",
+  "FT0CVariant1",
+  "FT0CVariant2",
+  "FDDM",
+  "NTPV",
+  "NGlobal",
+  "MFT"};
+
+// mcCentralityModule: 1st-order configurables
+struct coreConfigurables : o2::framework::ConfigurableGroup {
+  o2::framework::Configurable<o2::framework::LabeledArray<int>> enabledTables{"enabledTables",
+                                                                              {defaultParameters[0], nTablesConst, nParameters, tableNames, parameterNames},
+                                                                              "Produce this table: -1 for autodetect; otherwise, 0/1 is false/true"};
+  std::vector<int> mEnabledTables; // Vector of enabled tables
+
+  o2::framework::Configurable<bool> recalibrateCentrality{"recalibrateCentrality", false, "If true, re-calibrate the MC centrality for the binning in binPercentile."};
+  o2::framework::Configurable<int> recalibrateMode{"recalibrateMode", 1, "Strategy to calibrate MC centrality? 0: from low to high mult.; 1: from high to low mult."};
+  o2::framework::Configurable<int> minEntries{"minEntries", 100, "Minimum number of entries for estimating the mean value for the recalibration"};
+  o2::framework::Configurable<bool> doNotCrashOnNull{"doNotCrashOnNull", false, "If ccdb object does not exist, fill with dummy values"};
+
+  o2::framework::Configurable<bool> assignCentralityPerCandidate{"assignCentralityPerCandidate", false, "If true, assign centrality by sampling P(centrality|generated multiplicity) per MC collision instead of the class-averaged (statistical) value"};
+  o2::framework::Configurable<ULong64_t> centralitySamplingSeed{"centralitySamplingSeed", 137, "Base seed of the RNG used to sample centrality in candidate-by-candidate mode (assignCentralityPerCandidate)"};
+  o2::framework::ConfigurableAxis binsPercentile{"binsPercentile", {o2::framework::VARIABLE_WIDTH, 0, 1.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}, "Binning of the percentile axis"};
+  o2::framework::ConfigurableAxis binsPercentileFine{"binsPercentileFine", {o2::framework::VARIABLE_WIDTH, 0, 0.001, 0.01, 1.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}, "Binning of the percentile axis"};
+  o2::framework::ConfigurableAxis binsMultiplicity{"binsMultiplicity", {1000, 0, 5000}, "Binning of the multiplicity axis"};
+
+  // ccdb configurables
+  o2::framework::Configurable<std::string> path{"path", "/tmp/InputCalibMC.root", "path to calib file or ccdb path if begins with ccdb://"};
+
+  // debug option
+  o2::framework::Configurable<bool> verbose{"verbose", false, "If true, display more messages"};
+};
+
+struct products : o2::framework::ProducesGroup {
+  // Tables to produce
+  o2::framework::Produces<aod::McCentFV0As> centFV0A;
+  o2::framework::Produces<aod::McCentFT0Ms> centFT0M;
+  o2::framework::Produces<aod::McCentFT0As> centFT0A;
+  o2::framework::Produces<aod::McCentFT0Cs> centFT0C;
+  o2::framework::Produces<aod::McCentFT0CVariant1s> centFT0CVariant1;
+  o2::framework::Produces<aod::McCentFT0CVariant2s> centFT0CVariant2;
+  o2::framework::Produces<aod::McCentFDDMs> centFDDM;
+  o2::framework::Produces<aod::McCentNTPVs> centNTPV;
+  o2::framework::Produces<aod::McCentNGlobals> centNGlobal;
+  o2::framework::Produces<aod::McCentMFTs> centMFT;
+};
+
+template <typename T>
+concept HasMcMults = requires(typename T::iterator a) {
+  { a.multMCFT0A() } -> std::convertible_to<int>;
+  { a.multMCFT0C() } -> std::convertible_to<int>;
+  { a.multMCFV0A() } -> std::convertible_to<int>;
+  { a.multMCFDDA() } -> std::convertible_to<int>;
+  { a.multMCFDDC() } -> std::convertible_to<int>;
+  { a.multMCNParticlesEta08() } -> std::convertible_to<int>;
+  { a.multMCNParticlesEta05() } -> std::convertible_to<int>;
+};
+
+/// Task to produce the response table
+struct BuilderModule {
+  // Input parameters
+  o2::framework::Service<o2::ccdb::BasicCCDBManager> ccdb;
+
+  // declaration of structs here
+  // (N.B.: will be invisible to the outside, create your own copies)
+  coreConfigurables baseOpts;
+
+  TList* MCCentralityCalibObjects = nullptr;
+  TH1D* h1dFV0A = nullptr;
+  TH1D* h1dFT0M = nullptr;
+  TH1D* h1dFT0A = nullptr;
+  TH1D* h1dFT0C = nullptr;
+  TH1D* h1dFT0CVariant1 = nullptr;
+  TH1D* h1dFT0CVariant2 = nullptr;
+  TH1D* h1dFDDM = nullptr;
+  TH1D* h1dNTPV = nullptr;
+  TH1D* h1dNGlobal = nullptr;
+  TH1D* h1dMFT = nullptr;
+
+  // Calibration objects for candidate-by-candidate sampling (assignCentralityPerCandidate):
+  // x = centrality (%) of the matched reconstructed MC collision, y = generated multiplicity (|eta| < 0.5)
+  std::vector<TH2D*> h2dCentVsGenMult;
+  TRandom3 fRandomSampler;
+
+  // QA of the on-the-fly calibration (recalibrateCentrality): mean <N_PV>_{|eta|<0.5} in data and MC
+  // (reco), one entry per centrality class -- filled in extractCentralityCalibration().
+  // Must be TProfile, not TH1D+SetBinContent/SetBinError: every grid job redundantly derives the exact
+  // same calibration from the same CCDB input, and grid-merged output histograms are combined via a
+  // plain bin-content sum (TH1::Merge), which would turn a mean into N_jobs*mean. TProfile stores the
+  // underlying sum(w)/sum(wy)/sum(wy^2) per bin and merges those correctly, so the mean/error recomputed
+  // from the merged output stay correct however many jobs contributed. For the same reason there is no
+  // ratio histogram here: a ratio of two already-merged means can't be reconstructed after the fact, so
+  // compute PVMC/PVData downstream from the merged output if you need it, not inside the task.
+  std::vector<std::shared_ptr<TProfile>> hCalibPVData;
+  std::vector<std::shared_ptr<TProfile>> hCalibPVMC;
+
+  int nEnabledTables = 0;
+  int mRunNumber;
+
+  // TAxis
+  int nCentBins;
+  std::vector<double> centralityBins;
+
+  // Registers the on-the-fly-calibration QA histograms for one estimator (no-op unless recalibrateCentrality
+  // is on): mean <N_PV>_{|eta|<0.5} in data and in MC (reco), one entry per centrality class, as TProfile
+  // (see the comment on hCalibPVData/hCalibPVMC for why TH1D+SetBinContent doesn't survive grid merging).
+  // Uses the runtime add() overload (returning the type-erased HistPtr) since the estimator name/directory
+  // is only known at runtime here, unlike fillHistograms() which relies on a compile-time HIST() lookup.
+  template <typename THistoRegistry>
+  void registerCalibQAHistos(THistoRegistry& histos, int idx, const char* name)
+  {
+    if (!baseOpts.recalibrateCentrality) {
+      return;
+    }
+    o2::framework::AxisSpec axisCentClass{baseOpts.binsPercentile, Form("%s percentile (%%)", name)};
+    hCalibPVData[idx] = std::get<std::shared_ptr<TProfile>>(histos.add(Form("%s/calibPVData", name), Form("#LT #it{N}_{PV}^{Data} #GT vs %s percentile;%s percentile (%%);#LT #it{N}_{PV}^{Data} #GT_{|#it{#eta}|<0.5}", name, name), o2::framework::HistType::kTProfile, {axisCentClass}));
+    hCalibPVMC[idx] = std::get<std::shared_ptr<TProfile>>(histos.add(Form("%s/calibPVMC", name), Form("#LT #it{N}_{PV}^{MC} #GT (reco) vs %s percentile;%s percentile (%%);#LT #it{N}_{PV}^{MC} #GT_{|#it{#eta}|<0.5} (reco)", name, name), o2::framework::HistType::kTProfile, {axisCentClass}));
+  }
+
+  // Registers the candidate-by-candidate sampling QA (no-op unless assignCentralityPerCandidate is on):
+  // sampled percentile vs the generated multiplicity (|eta| < 0.5) actually used to key the sampling, so
+  // it can be visually compared against the input calibration (hGenMultEta05VsCentrality<estimator>) to
+  // check that the sampled distribution reproduces the expected P(centrality|genMult). This is a plain
+  // TH2D (not a TProfile): unlike the calibration QA above, it stores a distribution to look at, not a
+  // mean to merge, so ordinary bin-content-sum merging across grid jobs is exactly what's wanted here.
+  template <typename THistoRegistry>
+  void registerCandidateQAHistos(THistoRegistry& histos, const char* name)
+  {
+    if (!baseOpts.assignCentralityPerCandidate) {
+      return;
+    }
+    histos.add(Form("%s/percentileVsGenMultEta05", name), Form("Sampled %s percentile vs generated mult (|#it{#eta}|<0.5);%s percentile (%%);generated mult (|#it{#eta}|<0.5)", name, name), o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, Form("%s percentile", name)}, {baseOpts.binsMultiplicity, "generated mult (|#eta|<0.5)"}});
+  }
+
+  template <typename TBaseConfigurables, typename THistoRegistry, typename TInitContext>
+  void init(TBaseConfigurables const& inputBaseOpts, THistoRegistry& histos, TInitContext& context)
+  {
+    // read in configurations from the task where it's used
+    // could be grouped even further, but should work
+    baseOpts = inputBaseOpts;
+
+    baseOpts.mEnabledTables.resize(nTablesConst, 0);
+    h2dCentVsGenMult.resize(nTablesConst, nullptr);
+    hCalibPVData.resize(nTablesConst);
+    hCalibPVMC.resize(nTablesConst);
+    fRandomSampler.SetSeed(baseOpts.centralitySamplingSeed);
+
+    LOGF(info, "Checking if MC centrality is required");
+    auto& workflows = context.services().template get<o2::framework::RunningWorkflowInfo const>();
+
+    nEnabledTables = 0;
+
+    TString listOfRequestors[nTablesConst];
+    for (int i = 0; i < nTablesConst; i++) {
+      int f = baseOpts.enabledTables->get(tableNames[i].c_str(), "enable");
+      if (f == 1) {
+        baseOpts.mEnabledTables[i] = 1;
+        listOfRequestors[i] = "manual enabling";
+        nEnabledTables++;
+      }
+      if (f == -1) {
+        // autodetect this table in other devices
+        for (o2::framework::DeviceSpec const& device : workflows.devices) {
+          // Step 1: check if this device subscribed to the V0data table
+          for (auto const& input : device.inputs) {
+            if (o2::framework::DataSpecUtils::partialMatch(input.matcher, o2::header::DataOrigin("AOD"))) {
+              auto&& [origin, description, version] = o2::framework::DataSpecUtils::asConcreteDataMatcher(input.matcher);
+              std::string tableNameWithVersion = tableNames[i];
+              if (version > 0) {
+                tableNameWithVersion += Form("_%03d", version);
+              }
+              if (input.matcher.binding == tableNameWithVersion) {
+                LOGF(info, "Device %s has subscribed to %s (version %i)", device.name, tableNames[i], version);
+                listOfRequestors[i].Append(Form("%s ", device.name.c_str()));
+                baseOpts.mEnabledTables[i] = 1;
+                nEnabledTables++;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    if (nEnabledTables == 0) {
+      LOGF(info, "MC centrality not required. Will suppress all functionality, including logs, from this point forward.");
+      return;
+    }
+    mRunNumber = 0;
+
+    // TAxis
+    o2::framework::AxisSpec axisBinsPercentile{baseOpts.binsPercentile, "axisBinsPercentile"};
+    nCentBins = axisBinsPercentile.binEdges.size() - 1;
+    for (std::size_t iCent = 0; iCent < axisBinsPercentile.binEdges.size(); iCent++) {
+      centralityBins.push_back(axisBinsPercentile.binEdges[iCent]);
+    }
+
+    if (baseOpts.mEnabledTables[kFV0A]) {
+      histos.add("FV0A/percentile", "FV0A percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FV0A percentile"}});
+      histos.add("FV0A/percentilevsMult", "FV0A percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FV0A percentile"}, {baseOpts.binsMultiplicity, "FV0A mult."}});
+      registerCalibQAHistos(histos, kFV0A, "FV0A");
+      registerCandidateQAHistos(histos, "FV0A");
+    }
+    if (baseOpts.mEnabledTables[kFT0M]) {
+      histos.add("FT0M/percentile", "FT0M percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FT0M percentile"}});
+      histos.add("FT0M/percentilevsMult", "FT0M percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FT0M percentile"}, {baseOpts.binsMultiplicity, "FT0M mult."}});
+      registerCalibQAHistos(histos, kFT0M, "FT0M");
+      registerCandidateQAHistos(histos, "FT0M");
+    }
+    if (baseOpts.mEnabledTables[kFT0A]) {
+      histos.add("FT0A/percentile", "FT0A percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FT0A percentile"}});
+      histos.add("FT0A/percentilevsMult", "FT0A percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FT0A percentile"}, {baseOpts.binsMultiplicity, "FT0A mult."}});
+      registerCalibQAHistos(histos, kFT0A, "FT0A");
+      registerCandidateQAHistos(histos, "FT0A");
+    }
+    if (baseOpts.mEnabledTables[kFT0C]) {
+      histos.add("FT0C/percentile", "FT0C percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FT0C percentile"}});
+      histos.add("FT0C/percentilevsMult", "FT0C percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FT0C percentile"}, {baseOpts.binsMultiplicity, "FT0C mult."}});
+      registerCalibQAHistos(histos, kFT0C, "FT0C");
+      registerCandidateQAHistos(histos, "FT0C");
+    }
+    if (baseOpts.mEnabledTables[kFT0CVariant1]) {
+      histos.add("FT0CVariant1/percentile", "FT0CVariant1 percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FT0CVariant1 percentile"}});
+      histos.add("FT0CVariant1/percentilevsMult", "FT0CVariant1 percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FT0CVariant1 percentile"}, {baseOpts.binsMultiplicity, "FT0C mult."}});
+      registerCalibQAHistos(histos, kFT0CVariant1, "FT0CVariant1");
+      registerCandidateQAHistos(histos, "FT0CVariant1");
+    }
+    if (baseOpts.mEnabledTables[kFT0CVariant2]) {
+      histos.add("FT0CVariant2/percentile", "FT0CVariant2 percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FT0CVariant2 percentile"}});
+      histos.add("FT0CVariant2/percentilevsMult", "FT0CVariant2 percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FT0CVariant2 percentile"}, {baseOpts.binsMultiplicity, "FT0C mult."}});
+      registerCalibQAHistos(histos, kFT0CVariant2, "FT0CVariant2");
+      registerCandidateQAHistos(histos, "FT0CVariant2");
+    }
+    if (baseOpts.mEnabledTables[kFDDM]) {
+      histos.add("FDDM/percentile", "FDDM percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "FDDM percentile"}});
+      histos.add("FDDM/percentilevsMult", "FDDM percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "FDDM percentile"}, {baseOpts.binsMultiplicity, "FDDM mult."}});
+      registerCalibQAHistos(histos, kFDDM, "FDDM");
+      registerCandidateQAHistos(histos, "FDDM");
+    }
+    if (baseOpts.mEnabledTables[kNTPV]) {
+      histos.add("NTPV/percentile", "NTPV percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "NTPV percentile"}});
+      histos.add("NTPV/percentilevsMult", "NTPV percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "NTPV percentile"}, {baseOpts.binsMultiplicity, "NTPV mult."}});
+      registerCalibQAHistos(histos, kNTPV, "NTPV");
+      registerCandidateQAHistos(histos, "NTPV");
+    }
+    if (baseOpts.mEnabledTables[kNGlobal]) {
+      histos.add("NGlobal/percentile", "NGlobal percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "NGlobal percentile"}});
+      histos.add("NGlobal/percentilevsMult", "NGlobal percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "NGlobal percentile"}, {baseOpts.binsMultiplicity, "NGlobal mult."}});
+      registerCalibQAHistos(histos, kNGlobal, "NGlobal");
+      registerCandidateQAHistos(histos, "NGlobal");
+    }
+    if (baseOpts.mEnabledTables[kMFT]) {
+      histos.add("MFT/percentile", "MFT percentile.", o2::framework::HistType::kTH1D, {{baseOpts.binsPercentileFine, "MFT percentile"}});
+      histos.add("MFT/percentilevsMult", "MFT percentile.", o2::framework::HistType::kTH2D, {{baseOpts.binsPercentileFine, "MFT percentile"}, {baseOpts.binsMultiplicity, "MFT mult."}});
+    }
+
+    LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*");
+    LOGF(info, " MC centrality: basic configuration listing");
+    LOGF(info, "*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*");
+    // list enabled tables
+    for (int i = 0; i < nTablesConst; i++) {
+      // printout to be improved in the future
+      if (baseOpts.mEnabledTables[i]) {
+        LOGF(info, " -~> Table enabled: %s, requested by %s", tableNames[i], listOfRequestors[i].Data());
+      }
+    }
+  }
+
+  template <typename THist>
+  THist* getHist(const char* name)
+  {
+    if (!this->MCCentralityCalibObjects) {
+      return (THist*)0x0;
+    }
+    auto hist = reinterpret_cast<THist*>(this->MCCentralityCalibObjects->FindObject(name));
+    if (!hist) {
+      if (this->baseOpts.verbose) {
+        this->MCCentralityCalibObjects->ls();
+      }
+      if (this->baseOpts.doNotCrashOnNull) {
+        LOG(info) << "Could not open histogram " << name << " from TList, will fill tables with dummy values";
+      } else {
+        LOG(fatal) << "Could not open histogram " << name << " from TList";
+      }
+    }
+    return hist;
+  }
+
+  TH1D* extractCentralityCalibration(int idx, const char* name, bool reverse = false)
+  {
+
+    auto CalibMC = [&](TString const& estimator) {
+      LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Starting...\n", estimator.Data());
+
+      std::vector<double> percentile_center(nCentBins);
+      std::vector<double> epercentile_center(nCentBins);
+
+      // Histograms
+      TH2D* h2dMultVsCent_Data = getHist<TH2D>(Form("h2dMultVsCent%s_Data", estimator.Data()));
+      TH2D* h2dMultRecoVsMultGen_MC = getHist<TH2D>(Form("hMultEta05VsGenMult%s", estimator.Data()));
+      if (!h2dMultVsCent_Data || !h2dMultRecoVsMultGen_MC) {
+        return (TH1D*)0x0;
+      }
+
+      // QA of this on-the-fly calibration: mean <N_PV>_{|eta|<0.5} in data/MC per centrality class,
+      // registered by registerCalibQAHistos() (only non-null when recalibrateCentrality is on). Filled
+      // bin-by-bin from projData/projMC below (TProfile::Fill(x, y, weight)), not via SetBinContent, so
+      // that the mean/error reconstructed from the merged grid output stay statistically correct.
+      TProfile* hPVData = hCalibPVData[idx].get();
+      TProfile* hPVMC = hCalibPVMC[idx].get();
+
+      TH1D* h1dCalib = h2dMultRecoVsMultGen_MC->ProjectionX(Form("h1d%s", estimator.Data()), 1, h2dMultRecoVsMultGen_MC->GetNbinsX());
+      h1dCalib->Reset();
+      h1dCalib->SetTitle(Form("%s calibration object", estimator.Data()));
+      h1dCalib->GetXaxis()->SetTitle(Form("#it{N}_{%s, gen.}", estimator.Data()));
+      h1dCalib->GetYaxis()->SetTitle(Form("%s percentile (%%)", estimator.Data()));
+
+      // NOTE: candidate-by-candidate assignment (assignCentralityPerCandidate) does not go through this
+      // mean-matching path at all. It samples directly from the "hGenMultEta05VsCentrality<estimator>"
+      // joint histogram (x = reco centrality of the matched MC collision, y = generated mult |eta|<0.5)
+      // retrieved in initCCDB() below, see h2dCentVsGenMult and dataProcess().
+
+      if (reverse) {
+        for (int i = 0; i < nCentBins; i++) {
+          int irev = i;
+          percentile_center[i] = (centralityBins[irev] + centralityBins[irev + 1]) / 2;
+          epercentile_center[i] = (centralityBins[irev] - centralityBins[irev + 1]) / 2;
+        }
+
+        int startBinMc = h1dCalib->GetNbinsX();
+        for (int i = 0; i < nCentBins; i++) { // Loop over centrality bins
+          // start from the end (from the high multiplicity collisions)
+          int irev = i;
+          TH1D* projData = h2dMultVsCent_Data->ProjectionY(Form("projData_%d", i), h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev] + 1e-5),
+                                                           h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev + 1] - 1e-5));
+
+          double meanMult_Data = projData->GetMean();
+
+          double meanMult_MC = -1;
+          double diffMultDataVsMC = 1e+09;
+          int endBinMc = startBinMc;
+          for (int j = startBinMc; j >= 1; j--) {
+            // Loop over MC bins
+            TH1D* projMC = h2dMultRecoVsMultGen_MC->ProjectionY("", j, startBinMc);
+            int nEntries = projMC->Integral();
+            double meanMC = projMC->GetMean();
+            double ldiff = std::abs(meanMC - meanMult_Data);
+            // std::cout << ldiff << std::endl;
+
+            // if less than 100 entries for estimating the mean, do not consider it
+            if (nEntries < baseOpts.minEntries)
+              continue;
+
+            // Find the minimum difference and corresponding MC bin
+            if (ldiff < diffMultDataVsMC) {
+              diffMultDataVsMC = ldiff;
+              meanMult_MC = meanMC;
+              endBinMc = j;
+            }
+          }
+          if (i == nCentBins - 1) {
+            endBinMc = 1;
+          }
+
+          TH1D* projMC = h2dMultRecoVsMultGen_MC->ProjectionY(Form("projMC_%d", i), endBinMc, startBinMc);
+
+          // Refill the per-class mean from the underlying (data/MC) multiplicity distributions, bin by
+          // bin, so the TProfile accumulates real sum(w)/sum(wy)/sum(wy^2) rather than a precomputed mean.
+          if (hPVData) {
+            for (int ybin = 1; ybin <= projData->GetNbinsX(); ybin++) {
+              double w = projData->GetBinContent(ybin);
+              if (w > 0) {
+                hPVData->Fill(percentile_center[i], projData->GetBinCenter(ybin), w);
+              }
+            }
+          }
+          if (hPVMC) {
+            for (int ybin = 1; ybin <= projMC->GetNbinsX(); ybin++) {
+              double w = projMC->GetBinContent(ybin);
+              if (w > 0) {
+                hPVMC->Fill(percentile_center[i], projMC->GetBinCenter(ybin), w);
+              }
+            }
+          }
+
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Data centrality bin \e[1;31m%g-%g%%\e[0;00m\n", estimator.Data(), centralityBins[irev], centralityBins[irev + 1]);
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> MC multiplicity range %g-%g\n", estimator.Data(), projMC->GetBinLowEdge(endBinMc), projMC->GetBinLowEdge(startBinMc));
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> <PV> data = %.4f Vs <PV> MC = %.4f (MC/Data = %.4f%%)\n", estimator.Data(), meanMult_Data, meanMult_MC, (meanMult_MC - meanMult_Data) * 100 / meanMult_Data);
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> N entries Data = %g Vs N entries MC = %g\n", estimator.Data(), projData->Integral(), projMC->Integral());
+          LOGF(info, "\n");
+          for (int ibin = 1; ibin <= h1dCalib->GetNbinsX(); ibin++) {
+            if (ibin <= startBinMc && ibin >= endBinMc) {
+              h1dCalib->SetBinContent(ibin, percentile_center[i]);
+            }
+          }
+          startBinMc = endBinMc;
+        } // End loop over centrality bins
+      } else {
+        for (int i = 0; i < nCentBins; i++) {
+          int irev = nCentBins - i;
+          percentile_center[i] = (centralityBins[irev - 1] + centralityBins[irev]) / 2;
+          epercentile_center[i] = (centralityBins[irev] - centralityBins[irev - 1]) / 2;
+        }
+        int startBinMc = 1;
+        for (int i = 0; i < nCentBins; i++) { // Loop over centrality bins
+          // start from the end (from the low multiplicity collisions)
+          int irev = nCentBins - i;
+          TH1D* projData = h2dMultVsCent_Data->ProjectionY(Form("projData_%d", i), h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev - 1] + 1e-5),
+                                                           h2dMultVsCent_Data->GetXaxis()->FindBin(centralityBins[irev] - 1e-5));
+
+          double meanMult_Data = projData->GetMean();
+
+          double meanMult_MC = -1;
+          double diffMultDataVsMC = 1e+09;
+          int endBinMc = h1dCalib->GetNbinsX();
+          for (int j = startBinMc; j <= h1dCalib->GetNbinsX(); j++) {
+            // Loop over MC bins
+            TH1D* projMC = h2dMultRecoVsMultGen_MC->ProjectionY("", startBinMc, endBinMc);
+            int nEntries = projMC->Integral();
+            double meanMC = projMC->GetMean();
+            double ldiff = std::abs(meanMC - meanMult_Data);
+            // std::cout << ldiff << std::endl;
+
+            // if less than 100 entries for estimating the mean, do not consider it
+            if (nEntries < baseOpts.minEntries)
+              continue;
+
+            // Find the minimum difference and corresponding MC bin
+            if (ldiff < diffMultDataVsMC) {
+              diffMultDataVsMC = ldiff;
+              meanMult_MC = meanMC;
+              endBinMc = j;
+            }
+          }
+          if (i == nCentBins - 1) {
+            endBinMc = h2dMultRecoVsMultGen_MC->GetNbinsX();
+          }
+
+          TH1D* projMC = h2dMultRecoVsMultGen_MC->ProjectionY(Form("projMC_%d", i), startBinMc, endBinMc);
+
+          // Refill the per-class mean from the underlying (data/MC) multiplicity distributions, bin by
+          // bin, so the TProfile accumulates real sum(w)/sum(wy)/sum(wy^2) rather than a precomputed mean.
+          if (hPVData) {
+            for (int ybin = 1; ybin <= projData->GetNbinsX(); ybin++) {
+              double w = projData->GetBinContent(ybin);
+              if (w > 0) {
+                hPVData->Fill(percentile_center[i], projData->GetBinCenter(ybin), w);
+              }
+            }
+          }
+          if (hPVMC) {
+            for (int ybin = 1; ybin <= projMC->GetNbinsX(); ybin++) {
+              double w = projMC->GetBinContent(ybin);
+              if (w > 0) {
+                hPVMC->Fill(percentile_center[i], projMC->GetBinCenter(ybin), w);
+              }
+            }
+          }
+
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Data centrality bin \e[1;31m%g-%g%%\e[0;00m\n", estimator.Data(), centralityBins[irev - 1], centralityBins[irev]);
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> MC multiplicity range %g-%g\n", estimator.Data(), projMC->GetBinLowEdge(startBinMc), projMC->GetBinLowEdge(endBinMc + 1));
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> <PV> data = %.4f Vs <PV> MC = %.4f\n", estimator.Data(), meanMult_Data, meanMult_MC);
+          LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> N entries Data = %g Vs N entries MC = %g\n", estimator.Data(), projData->Integral(), projMC->Integral());
+          for (int ibin = 1; ibin <= h1dCalib->GetNbinsX(); ibin++) {
+            if (ibin <= endBinMc && ibin >= startBinMc) {
+              h1dCalib->SetBinContent(ibin, percentile_center[i]);
+            }
+          }
+
+          startBinMc = endBinMc + 1;
+        } // End loop over centrality bins
+      }
+
+      LOGF(info, "\e[1;31mCalibration for %s estimator\e[0;00m -> Done!\n", estimator.Data());
+
+      return h1dCalib;
+    };
+    return CalibMC(name);
+  }
+
+  // Loads the calibration object(s) for a single estimator: either the class-averaged TH1D (statistical
+  // assignment) or, in candidate-by-candidate mode, the raw joint TH2D used to sample P(centrality|genMult).
+  void loadEstimatorCalibration(int idx, const char* name, TH1D*& h1dOut)
+  {
+    if (!baseOpts.mEnabledTables[idx]) {
+      return;
+    }
+    if (baseOpts.assignCentralityPerCandidate) {
+      // x = reco centrality of the matched MC collision, y = generated multiplicity (|eta| < 0.5)
+      // filled directly by the calibration-producing task (e.g. centralityQa.cxx); used as-is, no
+      // mean-matching required since it is already a genuine joint distribution.
+      h2dCentVsGenMult[idx] = getHist<TH2D>(Form("hGenMultEta05VsCentrality%s", name));
+      return;
+    }
+    h1dOut = baseOpts.recalibrateCentrality ? extractCentralityCalibration(idx, name, baseOpts.recalibrateMode) : getHist<TH1D>(Form("h1d%s", name));
+  }
+
+  template <typename TCCDB, typename TBCs>
+  bool initCCDB(TCCDB& ccdb, TBCs const& bcs)
+  {
+    if (!bcs.size()) {
+      LOGF(warn, "No BC found, skipping this DF.");
+      return false; // signal to skip this DF
+    }
+
+    if (mRunNumber == bcs.iteratorAt(0).runNumber()) {
+      return true;
+    }
+    // mark this run as configured
+    mRunNumber = bcs.iteratorAt(0).runNumber();
+
+    MCCentralityCalibObjects = ccdb->template getForRun<TList>(baseOpts.path, mRunNumber);
+    if (MCCentralityCalibObjects) {
+      LOGF(info, "loaded TList with this many objects: %i", MCCentralityCalibObjects->GetEntries());
+    }
+
+    loadEstimatorCalibration(kFV0A, "FV0A", h1dFV0A);
+    loadEstimatorCalibration(kFT0M, "FT0M", h1dFT0M);
+    loadEstimatorCalibration(kFT0A, "FT0A", h1dFT0A);
+    loadEstimatorCalibration(kFT0C, "FT0C", h1dFT0C);
+    loadEstimatorCalibration(kFT0CVariant1, "FT0CVariant1", h1dFT0CVariant1);
+    loadEstimatorCalibration(kFT0CVariant2, "FT0CVariant2", h1dFT0CVariant2);
+    loadEstimatorCalibration(kFDDM, "FDDM", h1dFDDM);
+    loadEstimatorCalibration(kNTPV, "NTPV", h1dNTPV);
+    loadEstimatorCalibration(kNGlobal, "NGlobal", h1dNGlobal);
+    if (baseOpts.mEnabledTables[kMFT]) {
+      // to be added later (multMCMFT not yet available), kept as a dedicated block on purpose
+    }
+    LOG(info) << "Fully configured for run: " << mRunNumber;
+
+    return true;
+  }
+
+  template <int tableIndex>
+  void fillHistograms(o2::framework::HistogramRegistry& histos, double percentile, double multiplicity, double genMultEta05)
+  {
+    histos.fill(HIST(DirList[tableIndex]) + HIST("/percentile"), percentile);
+    histos.fill(HIST(DirList[tableIndex]) + HIST("/percentilevsMult"), percentile, multiplicity);
+    if (baseOpts.assignCentralityPerCandidate) {
+      // QA of the sampling itself: sampled percentile vs the generated multiplicity it was actually keyed
+      // on (registered by registerCandidateQAHistos()), so it can be compared to the input calibration.
+      histos.fill(HIST(DirList[tableIndex]) + HIST("/percentileVsGenMultEta05"), percentile, genMultEta05);
+    }
+  }
+
+  // Combines the run number, a per-DF collision sequence number and the estimator index into a single RNG
+  // seed, so that a given generated collision always draws the same sampled centrality for a given
+  // estimator, regardless of processing order (note: the sequence number is only unique within one DF,
+  // it is not a persistent table index).
+  // This is a simple hash-combine, not a cryptographic hash: the strides just keep the three fields from
+  // landing on the same seed value for the collision counts typically seen in one run, they don't
+  // guarantee it. A rare accidental collision only means two unrelated collisions share a random draw,
+  // which has no systematic effect on the sampled distribution.
+  ULong64_t computeSamplingSeed(int64_t collisionIndex, int tableIdx) const
+  {
+    static constexpr ULong64_t kCollisionStride = 131ull; // > kNestimators, so tableIdx can't alias into the collision term
+    static constexpr ULong64_t kRunStride = 1000003ull;   // prime, well above kCollisionStride * (typical collisions per DF)
+    return baseOpts.centralitySamplingSeed.value + kRunStride * static_cast<ULong64_t>(mRunNumber) +
+           kCollisionStride * static_cast<ULong64_t>(collisionIndex) + static_cast<ULong64_t>(tableIdx);
+  }
+
+  //__________________________________________________
+  template <typename TCCDB, typename THistoRegistry, typename TBCs, typename TProducts>
+  void dataProcess(TCCDB& ccdb, THistoRegistry& histos, TBCs const& bcs, const HasMcMults auto& mcCollisions, TProducts& products)
+  {
+    if (nEnabledTables == 0) {
+      return; // fully suppressed
+    }
+
+    if (!initCCDB(ccdb, bcs))
+      return;
+
+    // Local per-DF sequence number used only to seed the candidate-by-candidate sampling below: not every
+    // joined mcCollisions type carries an index column (globalIndex()), so we cannot rely on that instead.
+    int64_t mcCollisionCounter = 0;
+    for (const auto& mcCollision : mcCollisions) {
+      const double nFV0A = mcCollision.multMCFV0A();
+      const double nFT0A = mcCollision.multMCFT0A();
+      const double nFT0C = mcCollision.multMCFT0C();
+      const double nFDDA = mcCollision.multMCFDDA();
+      const double nFDDC = mcCollision.multMCFDDC();
+      const double nGlobal = mcCollision.multMCNParticlesEta08();
+      const double nGenMultEta05 = mcCollision.multMCNParticlesEta05();
+      // const double nMFT = mcCollision.multMCMFT(); // to be added later
+      const double nFT0M = nFT0A + nFT0C;
+      const double nFDDM = nFDDA + nFDDC;
+
+      auto populateTable = [&](auto& table, TH1D* h1dMultCalib, TH2D* h2dCentCalib, double multiplicity, auto tableIndex) {
+        double percentile = 105.0f;
+        if (baseOpts.mEnabledTables[tableIndex]) {
+          if (baseOpts.assignCentralityPerCandidate && h2dCentCalib) {
+            // Sample P(centrality | genMult) for this generated collision: slice h2dCentCalib at the
+            // y-bin (generated mult |eta|<0.5) matching this collision and draw from that 1-D distribution.
+            fRandomSampler.SetSeed(computeSamplingSeed(mcCollisionCounter, tableIndex.value));
+            int yBin = h2dCentCalib->GetYaxis()->FindBin(nGenMultEta05);
+            std::unique_ptr<TH1D> pdf{h2dCentCalib->ProjectionX(Form("pdfCent_%d", tableIndex.value), yBin, yBin)};
+            if (pdf->Integral() > 0) {
+              percentile = pdf->GetRandom(&fRandomSampler);
+            }
+            fillHistograms<tableIndex.value>(histos, percentile, multiplicity, nGenMultEta05);
+          } else if (h1dMultCalib) {
+            int bin = h1dMultCalib->FindBin(multiplicity);
+            percentile = h1dMultCalib->GetBinContent(bin);
+            fillHistograms<tableIndex.value>(histos, percentile, multiplicity, nGenMultEta05);
+          }
+          table(percentile);
+        }
+        return percentile;
+      };
+
+      populateTable(products.centFV0A, h1dFV0A, h2dCentVsGenMult[kFV0A], nFV0A, std::integral_constant<int, kFV0A>{});
+      populateTable(products.centFT0M, h1dFT0M, h2dCentVsGenMult[kFT0M], nFT0M, std::integral_constant<int, kFT0M>{});
+      populateTable(products.centFT0A, h1dFT0A, h2dCentVsGenMult[kFT0A], nFT0A, std::integral_constant<int, kFT0A>{});
+      populateTable(products.centFT0C, h1dFT0C, h2dCentVsGenMult[kFT0C], nFT0C, std::integral_constant<int, kFT0C>{});
+      populateTable(products.centFT0CVariant1, h1dFT0CVariant1, h2dCentVsGenMult[kFT0CVariant1], nFT0C, std::integral_constant<int, kFT0CVariant1>{});
+      populateTable(products.centFT0CVariant2, h1dFT0CVariant2, h2dCentVsGenMult[kFT0CVariant2], nFT0C, std::integral_constant<int, kFT0CVariant2>{});
+      populateTable(products.centFDDM, h1dFDDM, h2dCentVsGenMult[kFDDM], nFDDM, std::integral_constant<int, kFDDM>{});
+      populateTable(products.centNTPV, h1dNTPV, h2dCentVsGenMult[kNTPV], nGlobal, std::integral_constant<int, kNTPV>{});
+      populateTable(products.centNGlobal, h1dNGlobal, h2dCentVsGenMult[kNGlobal], nGlobal, std::integral_constant<int, kNGlobal>{});
+      // populateTable(products.centMFT, h1dMFT, h2dCentVsGenMult[kMFT], nMFT, std::integral_constant<int, kMFT>{}); // to be added later
+
+      mcCollisionCounter++;
+    }
+  }
+}; // end mcCentralityModule
+
+} // namespace mccentrality
+} // namespace pwglf
+} // namespace o2
+
+#endif // PWGLF_UTILS_MCCENTRALITYMODULE_