diff --git a/PWGHF/HFC/TableProducer/CMakeLists.txt b/PWGHF/HFC/TableProducer/CMakeLists.txt
index d0c1fb800ff..299981d8c1f 100644
--- a/PWGHF/HFC/TableProducer/CMakeLists.txt
+++ b/PWGHF/HFC/TableProducer/CMakeLists.txt
@@ -98,3 +98,8 @@ o2physics_add_dpl_workflow(producer-charm-hadrons-v0-femto-dream
                     SOURCES producerCharmHadronsV0FemtoDream.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore O2Physics::EventFilteringUtils
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(producer-charm-hadrons-charm-femto-dream
+                    SOURCES producerCharmHadronsCharmFemtoDream.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore O2Physics::EventFilteringUtils
+                    COMPONENT_NAME Analysis)
diff --git a/PWGHF/HFC/TableProducer/producerCharmHadronsCharmFemtoDream.cxx b/PWGHF/HFC/TableProducer/producerCharmHadronsCharmFemtoDream.cxx
new file mode 100644
index 00000000000..3a58131ad6a
--- /dev/null
+++ b/PWGHF/HFC/TableProducer/producerCharmHadronsCharmFemtoDream.cxx
@@ -0,0 +1,405 @@
+// 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 producerCharmHadronsCharmFemtoDream.cxx
+/// \brief Produce reduced D-meson candidate tables for charm-charm femtoscopy
+/// \author Biao Zhang, Heidelberg University, biao.zhang@cern.ch
+
+#include "PWGCF/DataModel/FemtoDerived.h"
+#include "PWGHF/Core/CentralityEstimation.h"
+#include "PWGHF/Core/HfMlResponseD0ToKPi.h"
+#include "PWGHF/Core/HfMlResponseDstarToD0Pi.h"
+#include "PWGHF/DataModel/CandidateReconstructionTables.h"
+#include "PWGHF/DataModel/CandidateSelectionTables.h"
+#include "PWGHF/Utils/utilsBfieldCCDB.h"
+#include "PWGHF/Utils/utilsEvSelHf.h"
+
+#include "Common/Core/ZorroSummary.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/Multiplicity.h"
+#include "Tools/ML/MlResponse.h"
+
+#include <CCDB/BasicCCDBManager.h>
+#include <CCDB/CcdbApi.h>
+#include <CommonConstants/PhysicsConstants.h>
+#include <DetectorsBase/Propagator.h>
+#include <Framework/ASoA.h>
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/Array2D.h>
+#include <Framework/Configurable.h>
+#include <Framework/Expressions.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/HistogramSpec.h>
+#include <Framework/InitContext.h>
+#include <Framework/runDataProcessing.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdint>
+#include <string>
+#include <vector>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+struct HfProducerCharmHadronsCharmFemtoDream {
+  enum MlMode : uint8_t {
+    NoMl = 0,
+    FillMlFromSelector = 1,
+    FillMlFromNewBDT = 2
+  };
+  static constexpr std::size_t NMlClasses = 3;
+  static constexpr std::size_t NMlPtEdgesMin = 2;
+  static constexpr int ND0Hypotheses = 2;
+  // Each species needs its own model, feature order and pT-dependent cuts.
+  struct HfMlConfig : ConfigurableGroup {
+    std::string prefix;
+    static inline const std::array<double, 3> defaultCuts{1., 0., 0.};
+    Configurable<int> applyMlMode{"applyMlMode", FillMlFromSelector, "0: no ML, 1: selector scores, 2: new BDT after selector"};
+    Configurable<std::vector<double>> binsPtMl{"binsPtMl", std::vector<double>{0., 36.}, "pT bin limits for new BDT"};
+    Configurable<LabeledArray<double>> cutsMl{"cutsMl", {defaultCuts.data(), 1, 3}, "New BDT cuts per pT bin: background, prompt, nonprompt"};
+    Configurable<std::vector<int>> cutDirMl{"cutDirMl", std::vector<int>{0, 1, 1}, "Reject scores above (0), below (1), or do not cut (2)"};
+    Configurable<int> nClassesMl{"nClassesMl", 3, "Three output classes: background, prompt, nonprompt"};
+    Configurable<std::vector<std::string>> namesInputFeatures{"namesInputFeatures", std::vector<std::string>{}, "Ordered input feature names for new BDT"};
+    Configurable<std::vector<std::string>> onnxFileNames{"onnxFileNames", std::vector<std::string>{}, "Model files, one per pT bin"};
+    Configurable<std::vector<std::string>> modelPathsCCDB{"modelPathsCCDB", std::vector<std::string>{}, "CCDB model paths, one per pT bin"};
+    Configurable<int64_t> timestampCCDB{"timestampCCDB", -1, "Timestamp used to retrieve models"};
+    Configurable<bool> loadModelsFromCCDB{"loadModelsFromCCDB", false, "Load new BDT from CCDB instead of local files"};
+  };
+
+  Produces<aod::FDCollisions> collisions;
+  Produces<aod::FDColMasks> rowMasks;
+  Produces<aod::FDHfCand2Prong> d0Rows;
+  Produces<aod::FDHfCandDstar> dstarRows;
+
+  Configurable<int> selectionFlagD0{"selectionFlagD0", 1, "Minimum D0 selector decision"};
+  Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "CCDB URL"};
+  Configurable<std::string> ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "Run 3 magnetic field"};
+  HfMlConfig mlD0{.prefix = "mlD0"}, mlDstar{.prefix = "mlDstar"};
+
+  Service<o2::ccdb::BasicCCDBManager> ccdb{};
+  o2::ccdb::CcdbApi ccdbApi;
+  o2::hf_evsel::HfEventSelection hfEvSel;
+  o2::analysis::HfMlResponseD0ToKPi<float> hfMlResponseD0;
+  o2::analysis::HfMlResponseDstarToD0Pi<float> hfMlResponseDstar;
+
+  using Collisions = soa::Join<aod::Collisions, aod::EvSels, aod::Mults>;
+  using CollisionsWithFT0C = soa::Join<Collisions, aod::CentFT0Cs>;
+  using D0s = soa::Join<aod::HfCand2Prong, aod::HfSelD0>;
+  using Dstars = soa::Join<aod::HfCandDstars, aod::HfD0FromDstar, aod::HfSelDstarToD0Pi>;
+
+  using D0sMl = soa::Join<aod::HfCand2ProngWPid, aod::HfSelD0, aod::HfMlD0>;
+  using DstarsMl = soa::Join<aod::HfCandDstarsWPid, aod::HfD0FromDstar, aod::HfSelDstarToD0Pi, aod::HfMlDstarToD0Pi>;
+
+  Filter filterSelectCandidateD0 = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0);
+  Filter filterSelectCandidateDstar = aod::hf_sel_candidate_dstar::isSelDstarToD0Pi == true;
+
+  HistogramRegistry registry{"registry"};
+  OutputObj<ZorroSummary> zorroSummary{"zorroSummary"};
+
+  bool hasD0 = false;
+  bool hasDstar = false;
+  int runNumber = -1;
+
+  template <typename Response>
+  void initMl(HfMlConfig const& cfg, Response& response, bool withMl)
+  {
+    if (cfg.applyMlMode.value < NoMl || cfg.applyMlMode.value > FillMlFromNewBDT) {
+      LOGP(fatal, "{}: invalid applyMlMode", cfg.prefix);
+    }
+    if (cfg.applyMlMode.value != FillMlFromNewBDT) {
+      return;
+    }
+    if (!withMl) {
+      LOGP(fatal, "{}: new BDT requires an Ml producer process", cfg.prefix);
+    }
+    auto const& edges = cfg.binsPtMl.value;
+    const bool areEdgesFinite = std::all_of(edges.begin(), edges.end(), [](double x) { return std::isfinite(x); });
+    const bool areEdgesIncreasing = std::adjacent_find(edges.begin(), edges.end(), [](double a, double b) { return a >= b; }) == edges.end();
+    if (edges.size() < NMlPtEdgesMin || !areEdgesFinite || !areEdgesIncreasing) {
+      LOGP(fatal, "{}: binsPtMl must be finite and strictly increasing", cfg.prefix);
+    }
+    const auto nBins = edges.size() - 1;
+    if (cfg.nClassesMl.value != NMlClasses || cfg.cutDirMl.value.size() != NMlClasses || cfg.cutsMl.value.rows() != nBins || cfg.cutsMl.value.cols() != NMlClasses || cfg.onnxFileNames.value.size() != nBins || cfg.namesInputFeatures.value.empty() || (cfg.loadModelsFromCCDB.value && cfg.modelPathsCCDB.value.size() != nBins)) {
+      LOGP(fatal, "{}: provide three classes, cuts/models for every pT bin, and input features", cfg.prefix);
+    }
+    for (const auto& direction : cfg.cutDirMl.value) {
+      if (direction < o2::cuts_ml::CutGreater || direction > o2::cuts_ml::CutNot) {
+        LOGP(fatal, "{}: invalid cutDirMl", cfg.prefix);
+      }
+    }
+    for (unsigned int bin = 0; bin < nBins; ++bin) {
+      for (unsigned int score = 0; score < NMlClasses; ++score) {
+        if (!std::isfinite(cfg.cutsMl.value.get(bin, score))) {
+          LOGP(fatal, "{}: cutsMl must be finite", cfg.prefix);
+        }
+      }
+    }
+    response.configure(edges, cfg.cutsMl.value, cfg.cutDirMl.value, cfg.nClassesMl.value);
+    response.cacheInputFeaturesIndices(cfg.namesInputFeatures.value);
+    if (cfg.loadModelsFromCCDB.value) {
+      ccdbApi.init(ccdbUrl.value);
+      response.setModelPathsCCDB(cfg.onnxFileNames.value, ccdbApi, cfg.modelPathsCCDB.value, cfg.timestampCCDB.value);
+    } else {
+      response.setModelPathsLocal(cfg.onnxFileNames.value);
+    }
+    response.init();
+  }
+
+  void init(InitContext const&)
+  {
+    const std::array<bool, 8> processSwitches{doprocessD0D0, doprocessD0Dstar, doprocessD0D0Ml, doprocessD0DstarMl, doprocessD0D0WithFT0C, doprocessD0DstarWithFT0C, doprocessD0D0MlWithFT0C, doprocessD0DstarMlWithFT0C};
+    if (std::count(processSwitches.begin(), processSwitches.end(), true) != 1) {
+      LOGP(fatal, "Enable exactly one charm-charm producer process");
+    }
+    if (selectionFlagD0 < 1) {
+      LOGP(fatal, "selectionFlagD0 must be positive");
+    }
+    ccdb->setURL(ccdbUrl);
+    ccdb->setCaching(true);
+    ccdb->setLocalObjectValidityChecking();
+    const bool withMl = doprocessD0D0Ml || doprocessD0DstarMl || doprocessD0D0MlWithFT0C || doprocessD0DstarMlWithFT0C;
+    initMl(mlD0, hfMlResponseD0, withMl);
+    if (doprocessD0Dstar || doprocessD0DstarMl || doprocessD0DstarWithFT0C || doprocessD0DstarMlWithFT0C) {
+      initMl(mlDstar, hfMlResponseDstar, withMl);
+    }
+    hfEvSel.init(registry, &zorroSummary);
+    const AxisSpec eventStage{2, -0.5, 1.5, "stage (0=all,1=accepted)"};
+    const AxisSpec candidateSpecies{4, -0.5, 3.5, "species (0=D0,1=D0bar,2=Dstar+,3=Dstar-)"};
+    registry.add("events", "Events", kTH1F, {eventStage});
+    registry.add("candidates", "Candidates", kTH1F, {candidateSpecies});
+  }
+
+  template <o2::hf_centrality::CentralityEstimator CentEstimator, typename CollType>
+  bool acceptCollision(CollType const& col)
+  {
+    registry.fill(HIST("events"), 0);
+    float cent = -1.f; // No centrality for pp MB.
+    const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask<true, CentEstimator, aod::BCsWithTimestamps>(col, cent, ccdb, registry);
+    hfEvSel.fillHistograms(col, rejectionMask, cent);
+    if (rejectionMask != 0) {
+      return false;
+    }
+    hasD0 = false;
+    hasDstar = false;
+    auto bc = col.template bc_as<aod::BCsWithTimestamps>();
+    initCCDB(bc, runNumber, ccdb, ccdbPathGrpMag.value, nullptr, false);
+    // Propagator field is in kG; FemtoDream uses Tesla.
+    const float fieldTesla = 0.1f * o2::base::Propagator::Instance()->getNominalBz();
+    // MultV0M is the common reduced percentile column: FT0C or -1 (no
+    // estimator).
+    collisions(col.posZ(), cent, col.multNTracksPV(), 2.f, fieldTesla);
+    registry.fill(HIST("events"), 1);
+    // Keep all accepted collisions, including zero/single-candidate events.
+    return true;
+  }
+
+  template <typename Scores>
+  std::array<float, 3> readScores(Scores const& scores)
+  {
+    if (scores.size() != NMlClasses) {
+      LOGP(fatal, "Expected three selector ML scores: background, prompt, nonprompt");
+      return {-1.f, -1.f, -1.f};
+    }
+    for (const auto& score : scores) {
+      if (!std::isfinite(score)) {
+        LOGP(fatal, "Non-finite selector ML score");
+      }
+    }
+    return {scores[0], scores[1], scores[2]};
+  }
+
+  template <bool WithMl, typename CollType, typename CandsType>
+  void fillD0(CollType const& col, CandsType const& candidates)
+  {
+    const auto timestamp = col.template bc_as<aod::BCsWithTimestamps>().timestamp();
+    for (const auto& cand : candidates) {
+      auto p0 = cand.template prong0_as<aod::Tracks>();
+      auto p1 = cand.template prong1_as<aod::Tracks>();
+      // The OR filter accepts the row if either hypothesis passes. Write only
+      // the passing hypotheses, each with its own flavour and ML scores.
+      for (int hypothesis = 0; hypothesis < ND0Hypotheses; ++hypothesis) {
+        if ((hypothesis == 0 ? cand.isSelD0() : cand.isSelD0bar()) < selectionFlagD0) {
+          continue;
+        }
+        std::array<float, 3> scores{-1.f, -1.f, -1.f};
+        if constexpr (WithMl) {
+          if (mlD0.applyMlMode != NoMl) {
+            scores = hypothesis == 0 ? readScores(cand.mlProbD0()) : readScores(cand.mlProbD0bar());
+          }
+          if (mlD0.applyMlMode == FillMlFromNewBDT) {
+            // Do not call the ML response with an out-of-range model index.
+            if (!std::isfinite(cand.pt()) || cand.pt() < mlD0.binsPtMl.value.front() || cand.pt() >= mlD0.binsPtMl.value.back()) {
+              continue;
+            }
+            const int pdgCode = hypothesis == 0 ? o2::constants::physics::kD0 : -o2::constants::physics::kD0;
+            auto features = hfMlResponseD0.getInputFeatures<true>(cand, pdgCode);
+            std::vector<float> output;
+            if (!hfMlResponseD0.isSelectedMl(features, cand.pt(), output)) {
+              continue; // Reject this hypothesis only.
+            }
+            scores = readScores(output);
+          }
+        }
+        d0Rows(collisions.lastIndex(), timestamp, hypothesis == 0 ? 1 : -1, p0.globalIndex(), p1.globalIndex(), p0.pt(), p1.pt(), p0.eta(), p1.eta(), p0.phi(), p1.phi(), 1 << hypothesis, scores[0], scores[1], scores[2]);
+        hasD0 = true;
+        registry.fill(HIST("candidates"), hypothesis);
+      }
+    }
+  }
+
+  template <bool WithMl, typename CollType, typename CandsType>
+  void fillDstar(CollType const& col, CandsType const& candidates)
+  {
+    const auto timestamp = col.template bc_as<aod::BCsWithTimestamps>().timestamp();
+    for (const auto& cand : candidates) {
+      auto p0 = cand.template prong0_as<aod::Tracks>();
+      auto p1 = cand.template prong1_as<aod::Tracks>();
+      auto soft = cand.template prongPi_as<aod::Tracks>();
+      std::array<float, 3> scores{-1.f, -1.f, -1.f};
+      if constexpr (WithMl) {
+        if (mlDstar.applyMlMode != NoMl) {
+          scores = readScores(cand.mlProbDstarToD0Pi());
+        }
+        if (mlDstar.applyMlMode == FillMlFromNewBDT) {
+          if (!std::isfinite(cand.pt()) || cand.pt() < mlDstar.binsPtMl.value.front() || cand.pt() >= mlDstar.binsPtMl.value.back()) {
+            continue;
+          }
+          // Match the unswapped daughter convention used by D+track.
+          auto features = hfMlResponseDstar.getInputFeatures(cand, false);
+          std::vector<float> output;
+          if (!hfMlResponseDstar.isSelectedMl(features, cand.pt(), output)) {
+            continue;
+          }
+          scores = readScores(output);
+        }
+      }
+      dstarRows(collisions.lastIndex(), timestamp, soft.sign(), p0.globalIndex(), p1.globalIndex(), soft.globalIndex(), p0.pt(), p1.pt(), soft.pt(), p0.eta(), p1.eta(), soft.eta(), p0.phi(), p1.phi(), soft.phi(), 1, scores[0], scores[1], scores[2]);
+      hasDstar = true;
+      registry.fill(HIST("candidates"), soft.sign() > 0 ? 2 : 3);
+    }
+  }
+
+  void processD0D0(Collisions::iterator const& col,
+                   aod::BCsWithTimestamps const&,
+                   aod::Tracks const&,
+                   soa::Filtered<D0s> const& d0s)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::None>(col)) {
+      fillD0<false>(col, d0s);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0D0, "D0 only, data", true);
+
+  void processD0Dstar(Collisions::iterator const& col,
+                      aod::BCsWithTimestamps const&,
+                      aod::Tracks const&,
+                      soa::Filtered<D0s> const& d0s,
+                      soa::Filtered<Dstars> const& dstars)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::None>(col)) {
+      fillD0<false>(col, d0s);
+      fillDstar<false>(col, dstars);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0Dstar, "D0 and Dstar, data", false);
+
+  void processD0D0Ml(Collisions::iterator const& col,
+                     aod::BCsWithTimestamps const&,
+                     aod::Tracks const&,
+                     soa::Filtered<D0sMl> const& d0s)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::None>(col)) {
+      fillD0<true>(col, d0s);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0D0Ml, "D0 with selector ML scores", false);
+
+  void processD0DstarMl(Collisions::iterator const& col,
+                        aod::BCsWithTimestamps const&,
+                        aod::Tracks const&,
+                        soa::Filtered<D0sMl> const& d0s,
+                        soa::Filtered<DstarsMl> const& dstars)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::None>(col)) {
+      fillD0<true>(col, d0s);
+      fillDstar<true>(col, dstars);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0DstarMl, "D0 and Dstar with selector ML scores", false);
+
+  void processD0D0WithFT0C(CollisionsWithFT0C::iterator const& col,
+                           aod::BCsWithTimestamps const&,
+                           aod::Tracks const&,
+                           soa::Filtered<D0s> const& d0s)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::FT0C>(col)) {
+      fillD0<false>(col, d0s);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0D0WithFT0C, "D0 only, data with FT0C centrality", false);
+
+  void processD0DstarWithFT0C(CollisionsWithFT0C::iterator const& col,
+                              aod::BCsWithTimestamps const&,
+                              aod::Tracks const&,
+                              soa::Filtered<D0s> const& d0s,
+                              soa::Filtered<Dstars> const& dstars)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::FT0C>(col)) {
+      fillD0<false>(col, d0s);
+      fillDstar<false>(col, dstars);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0DstarWithFT0C, "D0 and Dstar, data with FT0C centrality", false);
+
+  void processD0D0MlWithFT0C(CollisionsWithFT0C::iterator const& col,
+                             aod::BCsWithTimestamps const&,
+                             aod::Tracks const&,
+                             soa::Filtered<D0sMl> const& d0s)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::FT0C>(col)) {
+      fillD0<true>(col, d0s);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0D0MlWithFT0C, "D0 with selector ML scores and FT0C centrality", false);
+
+  void processD0DstarMlWithFT0C(CollisionsWithFT0C::iterator const& col,
+                                aod::BCsWithTimestamps const&,
+                                aod::Tracks const&,
+                                soa::Filtered<D0sMl> const& d0s,
+                                soa::Filtered<DstarsMl> const& dstars)
+  {
+    if (acceptCollision<o2::hf_centrality::CentralityEstimator::FT0C>(col)) {
+      fillD0<true>(col, d0s);
+      fillDstar<true>(col, dstars);
+      rowMasks(hasD0 ? 1 : 0, hasDstar ? 1 : 0, 0);
+    }
+  }
+  PROCESS_SWITCH(HfProducerCharmHadronsCharmFemtoDream, processD0DstarMlWithFT0C, "D0 and Dstar with selector ML scores and FT0C centrality", false);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<HfProducerCharmHadronsCharmFemtoDream>(cfgc)};
+}
diff --git a/PWGHF/HFC/Tasks/CMakeLists.txt b/PWGHF/HFC/Tasks/CMakeLists.txt
index 4fa2b6cd959..69e279424c9 100644
--- a/PWGHF/HFC/Tasks/CMakeLists.txt
+++ b/PWGHF/HFC/Tasks/CMakeLists.txt
@@ -70,3 +70,8 @@ o2physics_add_dpl_workflow(task-flow
                     COMPONENT_NAME Analysis)
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR}/PWGCF/Core)
+
+o2physics_add_dpl_workflow(task-charm-hadrons-charm-femto-dream
+                    SOURCES taskCharmHadronsCharmFemtoDream.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
diff --git a/PWGHF/HFC/Tasks/taskCharmHadronsCharmFemtoDream.cxx b/PWGHF/HFC/Tasks/taskCharmHadronsCharmFemtoDream.cxx
new file mode 100644
index 00000000000..f051d75cf27
--- /dev/null
+++ b/PWGHF/HFC/Tasks/taskCharmHadronsCharmFemtoDream.cxx
@@ -0,0 +1,457 @@
+// 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 taskCharmHadronsCharmFemtoDream.cxx
+/// \brief Build DD pair femtoscopy distributions
+/// \author Biao Zhang, Heidelberg University, biao.zhang@cern.ch
+
+#include "PWGCF/DataModel/FemtoDerived.h"
+#include "PWGCF/FemtoDream/Core/femtoDreamMath.h"
+
+#include <CommonConstants/PhysicsConstants.h>
+#include <Framework/ASoA.h>
+#include <Framework/ASoAHelpers.h>
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/BinningPolicy.h>
+#include <Framework/Configurable.h>
+#include <Framework/Expressions.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/HistogramSpec.h>
+#include <Framework/InitContext.h>
+#include <Framework/runDataProcessing.h>
+
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <string>
+#include <vector>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::constants::physics;
+using namespace o2::analysis::femtoDream;
+
+struct HfTaskCharmHadronsCharmFemtoDream {
+  // Like/unlike sign refers to charm flavour; D0 is electrically neutral.
+  // Keep charge-conjugate channels separate so they can also be combined
+  // offline.
+  enum PairChannel : int {
+    D0D0LikeSign = 0,
+    D0barD0barLikeSign = 1,
+    D0D0barUnlikeSign = 2,
+    D0DstarPlusLikeSign = 3,
+    D0barDstarMinusLikeSign = 4,
+    D0DstarMinusUnlikeSign = 5,
+    D0barDstarPlusUnlikeSign = 6,
+    NPairChannels = 7
+  };
+  static constexpr int MixingBinPolicyMax = 2;
+
+  Produces<aod::FDHfCharm2Prong> rowFemtoResultCharm2Prong;
+  Produces<aod::FDHfCharmDstar> rowFemtoResultCharmDstar;
+  Produces<aod::FDHfColl> rowFemtoResultColl;
+
+  Configurable<float> ptMinD0{"ptMinD0", 0.f, "Minimum D0 pT"};
+  Configurable<float> ptMaxD0{"ptMaxD0", 36.f, "Maximum D0 pT"};
+  Configurable<float> ptMinDstar{"ptMinDstar", 0.f, "Minimum Dstar pT"};
+  Configurable<float> ptMaxDstar{"ptMaxDstar", 36.f, "Maximum Dstar pT"};
+  Configurable<float> etaMax{"etaMax", 0.8f, "Maximum absolute candidate eta"};
+  Configurable<float> massMinD0{"massMinD0", 1.7f, "Minimum D0 mass (keep sidebands)"};
+  Configurable<float> massMaxD0{"massMaxD0", 2.0f, "Maximum D0 mass"};
+  Configurable<float> deltaMassMin{"deltaMassMin", MassPiPlus, "Minimum Dstar-D0 mass difference"};
+  Configurable<float> deltaMassMax{"deltaMassMax", 0.17f, "Maximum Dstar-D0 mass difference"};
+  Configurable<float> daughterMassMin{"daughterMassMin", 1.80f, "Minimum Dstar daughter D0 mass"};
+  Configurable<float> daughterMassMax{"daughterMassMax", 1.93f, "Maximum Dstar daughter D0 mass"};
+  Configurable<bool> useMl{"useMl", false, "Require valid ML scores and apply score cuts"};
+  Configurable<float> d0BkgBdtMax{"d0BkgBdtMax", 1.f, "Maximum D0 background score"};
+  Configurable<float> d0PromptBdtMin{"d0PromptBdtMin", 0.f, "Minimum D0 prompt score"};
+  Configurable<float> dstarBkgBdtMax{"dstarBkgBdtMax", 1.f, "Maximum Dstar background score"};
+  Configurable<float> dstarPromptBdtMin{"dstarPromptBdtMin", 0.f, "Minimum Dstar prompt score"};
+  Configurable<int> charmHadCandSel{"charmHadCandSel", 1, "Minimum reduced charm candidate selection flag"};
+  struct : ConfigurableGroup {
+    std::string prefix = "eventSel";
+    Configurable<bool> useCentrality{"useCentrality", false, "Apply percentile selection (requires a centrality-enabled producer)"};
+    Configurable<int> multMin{"multMin", 0, "Minimum multiplicity (MultNtr)"};
+    Configurable<int> multMax{"multMax", 99999, "Maximum multiplicity (MultNtr)"};
+    Configurable<float> multPercentileMin{"multPercentileMin", 0.f, "Minimum multiplicity percentile"};
+    Configurable<float> multPercentileMax{"multPercentileMax", 100.f, "Maximum multiplicity percentile"};
+  } eventSel;
+
+  struct : ConfigurableGroup {
+    Configurable<bool> doMixEvent{"doMixEvent", true, "Enable mixed events"};
+    Configurable<int> mixingBinPolicy{"mixingBinPolicy", 0, "0: multiplicity, 1: percentile, 2: both"};
+    Configurable<int> mixingDepth{"mixingDepth", 5, "Number of neighbours for O2 event mixing"};
+  } mixSetting;
+  SliceCache cache;
+  using FilteredCollisions = soa::Filtered<soa::Join<aod::FDCollisions, aod::FDColMasks>>;
+  using FilteredCharmCand2Prongs = soa::Filtered<aod::FDHfCand2Prong>;
+  using FilteredCharmCandDstars = soa::Filtered<aod::FDHfCandDstar>;
+
+  Filter eventMultiplicity = aod::femtodreamcollision::multNtr >= eventSel.multMin && aod::femtodreamcollision::multNtr <= eventSel.multMax;
+  Filter eventMultiplicityPercentile = ifnode(eventSel.useCentrality, aod::femtodreamcollision::multV0M >= eventSel.multPercentileMin && aod::femtodreamcollision::multV0M <= eventSel.multPercentileMax, Node{LiteralNode{true}});
+  Filter hfCandSelFilter = aod::fdhf::candidateSelFlag >= charmHadCandSel;
+
+  // Full-table process inputs do not automatically register grouping caches.
+  // Register both keys explicitly for the sliceByCached calls below.
+  Preslice<FilteredCharmCand2Prongs> perCollisionD0 = aod::femtodreamparticle::fdCollisionId;
+  Preslice<FilteredCharmCandDstars> perCollisionDstar = aod::femtodreamparticle::fdCollisionId;
+
+  Partition<FilteredCharmCand2Prongs> partitionCharmHadron2Prong = ifnode(useMl, aod::fdhf::bdtBkg <= d0BkgBdtMax && aod::fdhf::bdtPrompt >= d0PromptBdtMin, Node{LiteralNode{true}});
+  Partition<FilteredCharmCandDstars> partitionCharmHadronDstar = ifnode(useMl, aod::fdhf::bdtBkg <= dstarBkgBdtMax && aod::fdhf::bdtPrompt >= dstarPromptBdtMin, Node{LiteralNode{true}});
+
+  ConfigurableAxis mixingBinMult{"mixingBinMult", {VARIABLE_WIDTH, 0.f, 20.f, 60.f, 200.f}, "Mixing bins - multiplicity"};
+  ConfigurableAxis mixingBinMultPercentile{"mixingBinMultPercentile", {VARIABLE_WIDTH, 0.f, 100.f}, "Mixing bins - multiplicity percentile"};
+  ConfigurableAxis mixingBinVztx{"mixingBinVztx", {VARIABLE_WIDTH, -10.f, -4.f, 0.f, 4.f, 10.f}, "Mixing bins - z-vertex"};
+  ColumnBinningPolicy<aod::collision::PosZ, aod::femtodreamcollision::MultNtr> colBinningMult{{mixingBinVztx, mixingBinMult}, true};
+  ColumnBinningPolicy<aod::collision::PosZ, aod::femtodreamcollision::MultV0M> colBinningMultPercentile{{mixingBinVztx, mixingBinMultPercentile}, true};
+  ColumnBinningPolicy<aod::collision::PosZ, aod::femtodreamcollision::MultNtr, aod::femtodreamcollision::MultV0M> colBinningMultMultPercentile{{mixingBinVztx, mixingBinMult, mixingBinMultPercentile}, true};
+  aod::femtodreamcollision::BitMaskType bitMask = 1 << 0;
+  HistogramRegistry registry{"registry"};
+
+  void init(InitContext const&)
+  {
+    if (doprocessD0D0 == doprocessD0Dstar) {
+      LOGP(fatal, "Enable exactly one charm-charm analysis process");
+    }
+    if (mixSetting.mixingDepth < 0 || mixSetting.mixingBinPolicy < 0 || mixSetting.mixingBinPolicy > MixingBinPolicyMax || ptMinD0 < 0 || ptMinD0 >= ptMaxD0 || ptMinDstar < 0 || ptMinDstar >= ptMaxDstar || etaMax <= 0 || massMinD0 >= massMaxD0 || deltaMassMin >= deltaMassMax || daughterMassMin >= daughterMassMax || charmHadCandSel < 1 || eventSel.multMin > eventSel.multMax || eventSel.multPercentileMin > eventSel.multPercentileMax) {
+      LOGP(fatal, "Invalid charm-charm selection or mixing configuration");
+    }
+    colBinningMult = {{mixingBinVztx, mixingBinMult}, true};
+    colBinningMultPercentile = {{mixingBinVztx, mixingBinMultPercentile}, true};
+    colBinningMultMultPercentile = {{mixingBinVztx, mixingBinMult, mixingBinMultPercentile}, true};
+    const AxisSpec kstar{400, 0., 2., "k* (GeV/c)"};
+    const AxisSpec massD0{300, massMinD0.value, massMaxD0.value, "M(Kpi) (GeV/c2)"};
+    const AxisSpec deltaMass{310, deltaMassMin.value, deltaMassMax.value, "Delta M (GeV/c2)"};
+    const AxisSpec pt{72, ptMinD0.value, ptMaxD0.value, "D0 pT (GeV/c)"};
+    const AxisSpec ptStar{72, ptMinDstar.value, ptMaxDstar.value, "Dstar pT (GeV/c)"};
+    const AxisSpec mult{mixingBinMult, "NTracksPV"};
+    const AxisSpec channel{static_cast<int>(NPairChannels), -0.5, static_cast<double>(NPairChannels) - 0.5, "0=D0D0 LS,1=D0barD0bar LS,2=D0D0bar US,3=D0D*+ LS,4=D0barD*- LS,5=D0D*- US,6=D0barD*+ US"};
+    registry.add("SE/D0D0", "Same event", kTHnSparseF, {kstar, massD0, massD0, pt, pt, mult, channel});
+    registry.add("ME/D0D0", "Mixed event", kTHnSparseF, {kstar, massD0, massD0, pt, pt, mult, channel});
+    registry.add("SE/D0Dstar", "Same event", kTHnSparseF, {kstar, massD0, deltaMass, pt, ptStar, mult, channel});
+    registry.add("ME/D0Dstar", "Mixed event", kTHnSparseF, {kstar, massD0, deltaMass, pt, ptStar, mult, channel});
+    registry.add("QA/events", ";0=in mixing range,1=outside;events", kTH1F, {{2, -0.5, 1.5}});
+    registry.add("MixingQA/hSECollisionBins", ";mixing bin;events", kTH1F, {{1000, -0.5, 999.5}});
+    registry.add("MixingQA/hMECollisionBinsD0D0", ";mixing bin;event pairs", kTH1F, {{1000, -0.5, 999.5}});
+    registry.add("MixingQA/hMECollisionBinsD0Dstar", ";mixing bin;event pairs", kTH1F, {{1000, -0.5, 999.5}});
+    registry.add("QA/sharedDaughters", ";0=D0D0,1=D0Dstar;rejected pairs", kTH1F, {{2, -0.5, 1.5}});
+    registry.add("QA/d0MassPt", ";M(Kpi);pT", kTH2F, {massD0, pt});
+    registry.add("QA/dstarMassPt", ";Delta M;pT", kTH2F, {deltaMass, ptStar});
+  }
+
+  template <bool IsDstar, typename D0Row, typename OtherRow>
+  static bool sharesDaughter(D0Row const& a, OtherRow const& b)
+  {
+    for (const auto& id : std::array{a.prong0Id(), a.prong1Id()}) {
+      if (id < 0) {
+        continue;
+      }
+      if (id == b.prong0Id() || id == b.prong1Id()) {
+        return true;
+      }
+      if constexpr (IsDstar) {
+        if (id == b.prong2Id()) {
+          return true;
+        }
+      }
+    }
+    return false;
+  }
+
+  // D0 charge denotes flavour, not electric charge.
+  template <bool IsDstar, typename D0Row, typename OtherRow>
+  static PairChannel pairChannel(D0Row const& d0, OtherRow const& other)
+  {
+    if constexpr (!IsDstar) {
+      return d0.charge() != other.charge() ? D0D0barUnlikeSign : (d0.charge() > 0 ? D0D0LikeSign : D0barD0barLikeSign);
+    }
+    return d0.charge() > 0 ? (other.charge() > 0 ? D0DstarPlusLikeSign : D0DstarMinusUnlikeSign) : (other.charge() < 0 ? D0barDstarMinusLikeSign : D0barDstarPlusUnlikeSign);
+  }
+
+  // Canonical ordering keeps mass/pt axes independent of event ordering for
+  // D0D0.
+  template <typename FirstRow, typename SecondRow>
+  static bool reverseD0Order(FirstRow const& first, SecondRow const& second)
+  {
+    if (first.charge() != second.charge()) {
+      return first.charge() < second.charge(); // D0 first for D0-D0bar
+    }
+    return first.pt() < second.pt(); // leading pT first for equal flavours
+  }
+
+  template <typename CandsType>
+  void validateMlScores(CandsType const& candidates)
+  {
+    if (!useMl) {
+      return;
+    }
+    for (const auto& row : candidates) {
+      if (!std::isfinite(row.bdtBkg()) || !std::isfinite(row.bdtPrompt()) || row.bdtBkg() < 0.f || row.bdtPrompt() < 0.f) {
+        LOGP(fatal, "ML cuts requested on missing/invalid scores: use an ML producer process");
+      }
+    }
+  }
+
+  template <bool IsDstar, bool FillQa = false, typename Row>
+  bool select(Row const& row, float& mass)
+  {
+    if (std::abs(row.charge()) != 1 || !std::isfinite(row.pt()) || !std::isfinite(row.eta()) || !std::isfinite(row.phi()) || std::abs(row.eta()) >= etaMax || row.pt() < (IsDstar ? ptMinDstar.value : ptMinD0.value) || row.pt() >= (IsDstar ? ptMaxDstar.value : ptMaxD0.value)) {
+      return false;
+    }
+    const std::array<double, 2> masses = row.charge() > 0 ? std::array{MassPiPlus, MassKPlus} : std::array{MassKPlus, MassPiPlus};
+    if constexpr (IsDstar) {
+      const float daughterMass = row.mDaughD0(masses);
+      if (!std::isfinite(daughterMass) || daughterMass < daughterMassMin || daughterMass >= daughterMassMax) {
+        return false;
+      }
+      mass = row.m(std::array{masses[0], masses[1], MassPiPlus}) - daughterMass;
+    } else {
+      mass = row.m(masses);
+    }
+    if (!std::isfinite(mass) || mass < (IsDstar ? deltaMassMin.value : massMinD0.value) || mass >= (IsDstar ? deltaMassMax.value : massMaxD0.value)) {
+      return false;
+    }
+    if constexpr (IsDstar) {
+      if constexpr (FillQa) {
+        registry.fill(HIST("QA/dstarMassPt"), mass, row.pt());
+      }
+    } else {
+      if constexpr (FillQa) {
+        registry.fill(HIST("QA/d0MassPt"), mass, row.pt());
+      }
+    }
+    return true;
+  }
+
+  template <bool Mixed, bool IsDstar, typename FirstRow, typename SecondRow>
+  void fillPair(FirstRow const& first, SecondRow const& second, float firstMass, float secondMass, float mult)
+  {
+    if constexpr (!Mixed) {
+      if (sharesDaughter<IsDstar>(first, second)) {
+        registry.fill(HIST("QA/sharedDaughters"), IsDstar ? 1 : 0);
+        return;
+      }
+    }
+    const bool reverse = !IsDstar && reverseD0Order(first, second);
+    const float mass1 = reverse ? secondMass : firstMass;
+    const float mass2 = reverse ? firstMass : secondMass;
+    const float pt1 = reverse ? second.pt() : first.pt();
+    const float pt2 = reverse ? first.pt() : second.pt();
+    // Read kinematics directly from the candidate rows. Only D0D0 can reverse;
+    // its two nominal parent masses are equal. Never use Dstar delta mass here.
+    const float kstar = reverse ? FemtoDreamMath::getkstar(second, MassD0, first, MassD0) : FemtoDreamMath::getkstar(first, MassD0, second, IsDstar ? MassDStar : MassD0);
+    if (!std::isfinite(kstar)) {
+      return;
+    }
+    const int channel = pairChannel<IsDstar>(first, second);
+    if constexpr (IsDstar && Mixed) {
+      registry.fill(HIST("ME/D0Dstar"), kstar, mass1, mass2, pt1, pt2, mult, channel);
+    } else if constexpr (IsDstar) {
+      registry.fill(HIST("SE/D0Dstar"), kstar, mass1, mass2, pt1, pt2, mult, channel);
+    } else if constexpr (Mixed) {
+      registry.fill(HIST("ME/D0D0"), kstar, mass1, mass2, pt1, pt2, mult, channel);
+    } else {
+      registry.fill(HIST("SE/D0D0"), kstar, mass1, mass2, pt1, pt2, mult, channel);
+    }
+  }
+
+  template <typename CollType>
+  int mixingBin(CollType const& col)
+  {
+    switch (mixSetting.mixingBinPolicy) {
+      case aod::femtodreamcollision::kMult:
+        return colBinningMult.getBin({col.posZ(), col.multNtr()});
+      case aod::femtodreamcollision::kMultPercentile:
+        return colBinningMultPercentile.getBin({col.posZ(), col.multV0M()});
+      case aod::femtodreamcollision::kMultMultPercentile:
+        return colBinningMultMultPercentile.getBin({col.posZ(), col.multNtr(), col.multV0M()});
+      default:
+        return -1;
+    }
+  }
+
+  template <bool IsDstar, typename D0Table, typename OtherTable, typename Policy>
+  void doMixedEvent(FilteredCollisions const& cols, D0Table const& d0Rows, OtherTable const& otherRows, Policy const& policy)
+  {
+    Partition<FilteredCollisions> partitionMaskedCol1 = (aod::femtodreamcollision::bitmaskTrackOne & bitMask) == bitMask;
+    // TrackOne denotes D0 presence, TrackTwo denotes Dstar presence.
+    Partition<FilteredCollisions> partitionMaskedCol2 = (aod::femtodreamcollision::bitmaskTrackOne & bitMask) == bitMask;
+    Partition<FilteredCollisions> partitionMaskedColDstar = (aod::femtodreamcollision::bitmaskTrackTwo & bitMask) == bitMask;
+    partitionMaskedCol1.bindTable(cols);
+    if constexpr (IsDstar) {
+      partitionMaskedColDstar.bindTable(cols);
+    } else {
+      partitionMaskedCol2.bindTable(cols);
+    }
+    auto const& secondCollisions = IsDstar ? *partitionMaskedColDstar.mFiltered : *partitionMaskedCol2.mFiltered;
+    for (const auto& [collision1, collision2] : combinations(soa::CombinationsBlockFullIndexPolicy(policy, mixSetting.mixingDepth, -1, *partitionMaskedCol1.mFiltered, secondCollisions))) {
+      if (collision1.globalIndex() == collision2.globalIndex()) {
+        continue;
+      }
+      if constexpr (!IsDstar) {
+        // Full policy returns both orientations for identical inputs. Count
+        // once.
+        if (collision1.globalIndex() > collision2.globalIndex()) {
+          continue;
+        }
+      }
+      if (!std::isfinite(collision1.magField()) || collision1.magField() != collision2.magField()) {
+        continue;
+      }
+      if constexpr (IsDstar) {
+        registry.fill(HIST("MixingQA/hMECollisionBinsD0Dstar"), mixingBin(collision1));
+      } else {
+        registry.fill(HIST("MixingQA/hMECollisionBinsD0D0"), mixingBin(collision1));
+      }
+      auto first = d0Rows.sliceByCached(aod::femtodreamparticle::fdCollisionId, collision1.globalIndex(), cache);
+      auto second = otherRows.sliceByCached(aod::femtodreamparticle::fdCollisionId, collision2.globalIndex(), cache);
+      // Species roles remain fixed; the full policy supplies the reverse event
+      // orientation when eligible. Do not add a second manual reverse loop.
+      for (const auto& [row1, row2] : combinations(soa::CombinationsFullIndexPolicy(first, second))) {
+        float mass1{}, mass2{};
+        if (select<false>(row1, mass1) && select<IsDstar>(row2, mass2)) {
+          fillPair<true, IsDstar>(row1, row2, mass1, mass2, collision1.multNtr());
+        }
+      }
+    }
+  }
+
+  // Write candidates once per collision, independently of SE/ME pair counts.
+  // Reuse the same output tables and field conventions as the track/V0 tasks.
+  template <bool WithDstar, typename CollType, typename D0Slice, typename DstarSlice>
+  void fillTables(CollType const& col, D0Slice const& d0s, DstarSlice const& dstars)
+  {
+    int64_t timeStamp = -1;
+    bool hasCandidate = false;
+    auto recordTimeStamp = [&](int64_t value) {
+      if (hasCandidate && value != timeStamp) {
+        LOGP(fatal, "Inconsistent charm timestamps in one reduced collision");
+      }
+      timeStamp = value;
+      hasCandidate = true;
+    };
+    for (const auto& part : d0s) {
+      float mass{};
+      if (!select<false, true>(part, mass)) {
+        continue;
+      }
+      recordTimeStamp(part.timeStamp());
+      rowFemtoResultCharm2Prong(col.globalIndex(), timeStamp, mass, part.pt(), part.eta(), part.phi(), part.prong0Id(), part.prong1Id(), part.charge(), part.bdtBkg(), part.bdtPrompt(), part.bdtFD(), 0, 0);
+    }
+    if constexpr (WithDstar) {
+      for (const auto& part : dstars) {
+        float mass{};
+        if (!select<true, true>(part, mass)) {
+          continue;
+        }
+        recordTimeStamp(part.timeStamp());
+        const std::array<double, 2> daughterMasses = part.charge() > 0 ? std::array{MassPiPlus, MassKPlus} : std::array{MassKPlus, MassPiPlus};
+        // CharmM is delta mass, CharmDaughM is M(Kpi), as in D+Track.
+        rowFemtoResultCharmDstar(col.globalIndex(), timeStamp, mass, part.mDaughD0(daughterMasses), part.pt(), part.eta(), part.phi(), part.prong0Id(), part.prong1Id(), part.prong2Id(), part.charge(), part.bdtBkg(), part.bdtPrompt(), part.bdtFD(), 0, 0);
+      }
+    }
+    // Retain D0-only/Dstar-only events even when no clean SE pair exists.
+    // MC fields above are zero for these data processes, as in existing tasks.
+    if (hasCandidate) {
+      rowFemtoResultColl(col.globalIndex(), timeStamp, col.posZ(), col.multNtr());
+    }
+  }
+
+  template <bool WithDstar, typename D0Table, typename DstarTable, typename Policy>
+  void runWithPolicy(FilteredCollisions const& cols, D0Table const& d0Rows, DstarTable const& dstarRows, Policy const& policy)
+  {
+    for (const auto& col : cols) {
+      auto d0s = d0Rows.sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache);
+      // Export before online mixing-bin cuts, so offline mixing can change bins.
+      if constexpr (WithDstar) {
+        auto dstars = dstarRows.sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache);
+        fillTables<true>(col, d0s, dstars);
+      } else {
+        fillTables<false>(col, d0s, d0s);
+      }
+      const int bin = mixingBin(col);
+      if (bin < 0 || !std::isfinite(col.magField())) {
+        registry.fill(HIST("QA/events"), 1);
+        continue;
+      }
+      registry.fill(HIST("QA/events"), 0);
+      registry.fill(HIST("MixingQA/hSECollisionBins"), bin);
+      for (const auto& [row1, row2] : combinations(soa::CombinationsStrictlyUpperIndexPolicy(d0s, d0s))) {
+        float mass1{}, mass2{};
+        if (select<false>(row1, mass1) && select<false>(row2, mass2)) {
+          fillPair<false, false>(row1, row2, mass1, mass2, col.multNtr());
+        }
+      }
+      if constexpr (WithDstar) {
+        auto dstars = dstarRows.sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache);
+        for (const auto& [row1, row2] : combinations(soa::CombinationsFullIndexPolicy(d0s, dstars))) {
+          float mass1{}, mass2{};
+          if (select<false>(row1, mass1) && select<true>(row2, mass2)) {
+            fillPair<false, true>(row1, row2, mass1, mass2, col.multNtr());
+          }
+        }
+      }
+    }
+    if (mixSetting.doMixEvent && mixSetting.mixingDepth > 0) {
+      doMixedEvent<false>(cols, d0Rows, d0Rows, policy);
+      if constexpr (WithDstar) {
+        doMixedEvent<true>(cols, d0Rows, dstarRows, policy);
+      }
+    }
+  }
+
+  template <bool WithDstar, typename DstarTable>
+  void run(FilteredCollisions const& cols, FilteredCharmCand2Prongs const& d0Rows, DstarTable const& dstarRows)
+  {
+    validateMlScores(d0Rows);
+    partitionCharmHadron2Prong.bindTable(d0Rows);
+    if constexpr (WithDstar) {
+      validateMlScores(dstarRows);
+      partitionCharmHadronDstar.bindTable(dstarRows);
+    }
+    switch (mixSetting.mixingBinPolicy) {
+      case aod::femtodreamcollision::kMult:
+        runWithPolicy<WithDstar>(cols, partitionCharmHadron2Prong, partitionCharmHadronDstar, colBinningMult);
+        break;
+      case aod::femtodreamcollision::kMultPercentile:
+        runWithPolicy<WithDstar>(cols, partitionCharmHadron2Prong, partitionCharmHadronDstar, colBinningMultPercentile);
+        break;
+      case aod::femtodreamcollision::kMultMultPercentile:
+        runWithPolicy<WithDstar>(cols, partitionCharmHadron2Prong, partitionCharmHadronDstar, colBinningMultMultPercentile);
+        break;
+      default:
+        LOGP(fatal, "Invalid mixing binning policy");
+    }
+  }
+
+  void processD0D0(FilteredCollisions const& cols,
+                   FilteredCharmCand2Prongs const& d0s)
+  {
+    run<false>(cols, d0s, d0s);
+  }
+  PROCESS_SWITCH(HfTaskCharmHadronsCharmFemtoDream, processD0D0, "D0D0 data, SE and ME", true);
+
+  void processD0Dstar(FilteredCollisions const& cols,
+                      FilteredCharmCand2Prongs const& d0s,
+                      FilteredCharmCandDstars const& dstars)
+  {
+    run<true>(cols, d0s, dstars);
+  }
+  PROCESS_SWITCH(HfTaskCharmHadronsCharmFemtoDream, processD0Dstar, "D0D0 and D0Dstar data, SE and ME", false);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{adaptAnalysisTask<HfTaskCharmHadronsCharmFemtoDream>(cfgc)};
+}