alice3Multicharm.cxx

// 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 alice3Multicharm.cxx
/// \brief consumer task for alice 3 multicharm studies
/// \author Jesper Karlsson Gumprecht <jesper.gumprecht@cern.ch>

//  *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
//   Decay finder task for ALICE 3
//  *+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
//
//    Uses specific ALICE 3 PID and performance for studying
//    HF decays. Work in progress: use at your own risk!

#include "ALICE3/DataModel/OTFMulticharm.h"
#include "Tools/ML/model.h"

#include <CCDB/CcdbApi.h>
#include <Framework/AnalysisDataModel.h>
#include <Framework/AnalysisTask.h>
#include <Framework/Configurable.h>
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>
#include <Framework/InitContext.h>
#include <Framework/OutputObjHeader.h>
#include <Framework/runDataProcessing.h>

#include <TH1.h>

#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <map>
#include <string>
#include <vector>

using namespace o2;
using namespace o2::ml;
using namespace o2::framework;
using namespace o2::framework::expressions;

struct Alice3Multicharm {
  HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};

  o2::ml::OnnxModel bdtMCharm{};
  std::map<std::string, std::string> metadata{};
  o2::ccdb::CcdbApi ccdbApi{};
  static constexpr float ToMicrons = 1e+4;

  struct : ConfigurableGroup {
    std::string prefix = "bdt"; // JSON group name
    Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
    Configurable<std::string> localPath{"localPath", "MCharm_BDTModel.onnx", "(std::string) Path to the local .onnx file."};
    Configurable<std::string> pathCCDB{"pathCCDB", "Users/j/jekarlss/MLModels", "Path on CCDB"};
    Configurable<int64_t> timestampCCDB{"timestampCCDB", 1695750420200, "timestamp of the ONNX file for ML model used to query in CCDB. Please use 1695750420200"};
    Configurable<bool> loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"};
    Configurable<bool> enableOptimizations{"enableOptimizations", false, "Enables the ONNX extended model-optimization: sessionOptions.SetGraphOptimizationLevel(GraphOptimizationLevel::ORT_ENABLE_EXTENDED)"};
    Configurable<bool> enableML{"enableML", true, "Enables bdt model"};
  } bdt;

  ConfigurableAxis axisEta{"axisEta", {80, -4.0f, +4.0f}, "#eta"};
  ConfigurableAxis axisXicMass{"axisXicMass", {200, 2.368f, 2.568f}, "Xic Inv Mass (GeV/c^{2})"};
  ConfigurableAxis axisXiccMass{"axisXiccMass", {200, 3.521f, 3.721f}, "Xicc Inv Mass (GeV/c^{2})"};
  ConfigurableAxis axisDCA{"axisDCA", {400, 0, 400}, "DCA (#mum)"};
  ConfigurableAxis axisRadiusLarge{"axisRadiusLarge", {1000, 0, 20}, "Decay radius (cm)"};
  ConfigurableAxis axisRadius{"axisRadius", {10000, 0, 10000}, "Decay radius (#mum)"};
  ConfigurableAxis axisTofTrackDelta{"axisTofTrackDelta", {200, 0, 1000}, "TOF track time"};
  ConfigurableAxis axisNSigma{"axisNSigma", {21, -10, 10}, "nsigma"};
  ConfigurableAxis axisDecayLength{"axisDecayLength", {2000, 0, 2000}, "Decay lenght (#mum)"};
  ConfigurableAxis axisDcaDaughters{"axisDcaDaughters", {200, 0, 100}, "DCA (mum)"};
  ConfigurableAxis axisBDTScore{"axisBDTScore", {100, 0, 1}, "BDT Score"};
  ConfigurableAxis axisBDTScoreFine{"axisBDTScoreFine", {1000, 0, 1}, "BDT Score for 1D histogram"};
  ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"};

  Configurable<float> xiMinDCAxy{"xiMinDCAxy", -1, "[0] in |DCAxy| > [0]+[1]/pT"};
  Configurable<float> xiMinDCAz{"xiMinDCAz", -1, "[0] in |DCAz| > [0]+[1]/pT"};
  Configurable<float> xiMinRadius{"xiMinRadius", -1, "Minimum R2D for Xic decay (cm)"};

  Configurable<float> picMinDCAxy{"picMinDCAxy", -1, "[0] in |DCAz| > [0]+[1]/pT"};
  Configurable<float> picMinDCAz{"picMinDCAz", -1, "[0] in |DCAxy| > [0]+[1]/pT"};
  Configurable<float> picMinPt{"picMinPt", -1, "Minimum pT for Xic pions"};

  Configurable<float> piccMinDCAxy{"piccMinDCAxy", -1, "[0] in |DCAxy| > [0]+[1]/pT"};
  Configurable<float> piccMinDCAz{"piccMinDCAz", -1, "[0] in |DCAz| > [0]+[1]/pT"};
  Configurable<float> piccMinPt{"piccMinPt", -1, "Minimum pT for Xicc pions"};

  Configurable<float> xicMaxDauDCA{"xicMaxDauDCA", 1e+4, "DCA between Xic daughters (cm)"};
  Configurable<float> xicMinDCAxy{"xicMinDCAxy", -1, "[0] in |DCAz| > [0]+[1]/pT"};
  Configurable<float> xicMinDCAz{"xicMinDCAz", -1, "[0] in |DCAxy| > [0]+[1]/pT"};
  Configurable<float> xiccMaxDCAxy{"xiccMaxDCAxy", 1e+4, "Maximum DCAxy"};
  Configurable<float> xiccMaxDCAz{"xiccMaxDCAz", 1e+4, "Maximum DCAz"};
  Configurable<float> xicMinRadius{"xicMinRadius", -1, "Minimum R2D for Xic decay (cm)"};
  Configurable<float> xicMinDecayDistanceFromPV{"xicMinDecayDistanceFromPV", -1, "Minimum distance for Xic decay from PV (cm)"};
  Configurable<float> xicMinProperLength{"xicMinProperLength", -1, "Minimum proper length for Xic decay (cm)"};
  Configurable<float> xicMaxProperLength{"xicMaxProperLength", 1e+4, "Minimum proper length for Xic decay (cm)"};

  Configurable<float> xiccMaxDauDCA{"xiccMaxDauDCA", 1e+4, "DCA between Xicc daughters (cm)"};
  Configurable<float> xiccMinRadius{"xiccMinRadius", -1, "Minimum R2D for Xicc decay (cm)"};
  Configurable<float> xiccMinProperLength{"xiccMinProperLength", -1, "Minimum proper length for Xicc decay (cm)"};
  Configurable<float> xiccMaxProperLength{"xiccMaxProperLength", 1e+4, "Minimum proper length for Xicc decay (cm)"};
  Configurable<int> otfConfig{"otfConfig", 0, "OTF configuration flag"};

  Filter configFilter = (aod::otfmulticharm::lutConfigId == otfConfig);

  void init(InitContext&)
  {
    histos.add("CandidateQA/hDCAXicDaughters", "hDCAXicDaughters; DCA between Xic daughters (#mum)", kTH1D, {axisDcaDaughters});
    histos.add("CandidateQA/hDCAXiccDaughters", "hDCAXiccDaughters; DCA between Xicc daughters (#mum)", kTH1D, {axisDcaDaughters});
    histos.add("CandidateQA/hDCAxyXi", "hDCAxyXi; Xi DCAxy to PV (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hDCAzXi", "hDCAzXi; Xi DCAz to PV (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hDCAxyXic", "hDCAxyXic; Xic DCAxy to PV (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hDCAzXic", "hDCAzXic; Xic DCAz to PV (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hDCAxyXicc", "hDCAxyXicc; Xicc DCAxy to PV (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hDCAzXicc", "hDCAzXicc; Xicc DCAz to PV (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hDecayRadiusXic", "hDecayRadiusXic; Distance (#mum)", kTH1D, {axisRadius});
    histos.add("CandidateQA/hDecayRadiusXicc", "hDecayRadiusXicc; Distance (#mum)", kTH1D, {axisRadius});
    histos.add("CandidateQA/hDecayDistanceFromPVXic", "hDecayDistanceFromPVXic; Distance (#mum)", kTH1D, {axisDecayLength});
    histos.add("CandidateQA/hProperLengthXic", "hProperLengthXic; Distance (#mum)", kTH1D, {axisDecayLength});
    histos.add("CandidateQA/hProperLengthXicc", "hProperLengthXicc; Distance (#mum)", kTH1D, {axisDecayLength});
    histos.add("CandidateQA/hPi1cDCAxy", "hPi1cDCAxy; Pi1c DCAxy (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hPi1cDCAz", "hPi1cDCAz; Pi1c DCAz (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hPi2cDCAxy", "hPi2cDCAxy; Pi2c DCAxy (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hPi2cDCAz", "hPi2cDCAz; Pi2c DCAz (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hPiccDCAxy", "hPiccDCAxy; Picc DCAxy (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hPiccDCAz", "hPiccDCAz; Picc DCAz (#mum)", kTH1D, {axisDCA});
    histos.add("CandidateQA/hPi1cPt", "hPi1cPt; Pi1c pT (Gev/#it(c))", kTH1D, {axisPt});
    histos.add("CandidateQA/hPi2cPt", "hPi2cPt; Pi2c pT (Gev/#it(c))", kTH1D, {axisPt});
    histos.add("CandidateQA/hPiccPt", "hPiccPt; Picc pT (Gev/#it(c))", kTH1D, {axisPt});

    histos.add("SelectionQA/hDCAXicDaughters", "hDCAXicDaughters; DCA between Xic daughters (#mum)", kTH1D, {axisDcaDaughters});
    histos.add("SelectionQA/hDCAXiccDaughters", "hDCAXiccDaughters; DCA between Xicc daughters (#mum)", kTH1D, {axisDcaDaughters});
    histos.add("SelectionQA/hDCAxyXi", "hDCAxyXi; Xi DCAxy to PV (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hDCAzXi", "hDCAzXi; Xi DCAz to PV (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hDCAxyXic", "hDCAxyXic; Xic DCAxy to PV (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hDCAzXic", "hDCAzXic; Xic DCAz to PV (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hDCAxyXicc", "hDCAxyXicc; Xicc DCAxy to PV (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hDCAzXicc", "hDCAzXicc; Xicc DCAz to PV (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hDecayRadiusXic", "hDecayRadiusXic; Distance (#mum)", kTH1D, {axisRadius});
    histos.add("SelectionQA/hDecayRadiusXicc", "hDecayRadiusXicc; Distance (#mum)", kTH1D, {axisRadius});
    histos.add("SelectionQA/hDecayDistanceFromPVXic", "hDecayDistanceFromPVXic; Distance (#mum)", kTH1D, {axisDecayLength});
    histos.add("SelectionQA/hProperLengthXic", "hProperLengthXic; Distance (#mum)", kTH1D, {axisDecayLength});
    histos.add("SelectionQA/hProperLengthXicc", "hProperLengthXicc; Distance (#mum)", kTH1D, {axisDecayLength});
    histos.add("SelectionQA/hPi1cDCAxy", "hPi1cDCAxy; Pi1c DCAxy (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hPi1cDCAz", "hPi1cDCAz; Pi1c DCAz (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hPi2cDCAxy", "hPi2cDCAxy; Pi2c DCAxy (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hPi2cDCAz", "hPi2cDCAz; Pi2c DCAz (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hPiccDCAxy", "hPiccDCAxy; Picc DCAxy (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hPiccDCAz", "hPiccDCAz; Picc DCAz (#mum)", kTH1D, {axisDCA});
    histos.add("SelectionQA/hPi1cPt", "hPi1cPt; Pi1c pT (Gev/#it(c))", kTH1D, {axisPt});
    histos.add("SelectionQA/hPi2cPt", "hPi2cPt; Pi2c pT (Gev/#it(c))", kTH1D, {axisPt});
    histos.add("SelectionQA/hPiccPt", "hPiccPt; Picc pT (Gev/#it(c))", kTH1D, {axisPt});

    auto hMCharmBuilding = histos.add<TH1>("hMCharmBuilding", "hMCharmBuilding", kTH1D, {{22, -0.5, 21.5}});
    hMCharmBuilding->GetXaxis()->SetBinLabel(1, "nTotalCandidates");
    hMCharmBuilding->GetXaxis()->SetBinLabel(2, "xicMaxDauDCA");
    hMCharmBuilding->GetXaxis()->SetBinLabel(3, "xiccMaxDauDCA");
    hMCharmBuilding->GetXaxis()->SetBinLabel(4, "xiMinDCAxy");
    hMCharmBuilding->GetXaxis()->SetBinLabel(5, "xiMinDCAz");
    hMCharmBuilding->GetXaxis()->SetBinLabel(6, "pi1cMinDCAxy");
    hMCharmBuilding->GetXaxis()->SetBinLabel(7, "pi1cMinDCAz");
    hMCharmBuilding->GetXaxis()->SetBinLabel(8, "pi2cMinDCAxy");
    hMCharmBuilding->GetXaxis()->SetBinLabel(9, "pi2cMinDCAz");
    hMCharmBuilding->GetXaxis()->SetBinLabel(10, "piccMinDCAxy");
    hMCharmBuilding->GetXaxis()->SetBinLabel(11, "piccMinDCAz");
    hMCharmBuilding->GetXaxis()->SetBinLabel(12, "xicMinDCAxy");
    hMCharmBuilding->GetXaxis()->SetBinLabel(13, "xicMinDCAz");
    hMCharmBuilding->GetXaxis()->SetBinLabel(14, "xiccMaxDCAxy");
    hMCharmBuilding->GetXaxis()->SetBinLabel(15, "xiccMaxDCAz");
    hMCharmBuilding->GetXaxis()->SetBinLabel(16, "xicMinRadius");
    hMCharmBuilding->GetXaxis()->SetBinLabel(17, "xiccMinRadius");
    hMCharmBuilding->GetXaxis()->SetBinLabel(18, "xicMinProperLength");
    hMCharmBuilding->GetXaxis()->SetBinLabel(19, "xicMaxProperLength");
    hMCharmBuilding->GetXaxis()->SetBinLabel(20, "xiccMinProperLength");
    hMCharmBuilding->GetXaxis()->SetBinLabel(21, "xiccMaxProperLength");
    hMCharmBuilding->GetXaxis()->SetBinLabel(22, "xicMinDecayDistanceFromPV");

    histos.add("hXiccMass", "hXiccMass", kTH1D, {axisXiccMass});
    histos.add("hXicMass", "hXicMass", kTH1D, {axisXicMass});
    histos.add("hXiccPt", "hXiccPt", kTH1D, {axisPt});
    histos.add("hXicPt", "hXicPt", kTH1D, {axisPt});
    histos.add("h3dXicc", "h3dXicc; Xicc pT (GeV/#it(c)); Xicc #eta; Xicc mass (GeV/#it(c)^{2})", kTH3D, {axisPt, axisEta, axisXiccMass});
    histos.add("hConfigId", "hConfigId", kTH1D, {{11, -0.5, 10.5}});

    if (bdt.enableML) {
      if (bdt.loadModelsFromCCDB) {
        ccdbApi.init(bdt.ccdbUrl);
        LOG(info) << "Fetching model for timestamp: " << bdt.timestampCCDB.value;
        bool retrieveSuccessMCharm = ccdbApi.retrieveBlob(bdt.pathCCDB.value, ".", metadata, bdt.timestampCCDB.value, false, bdt.localPath.value);

        if (retrieveSuccessMCharm) {
          bdtMCharm.initModel(bdt.localPath.value, bdt.enableOptimizations.value);
        } else {
          LOG(fatal) << "Error encountered while fetching/loading the MCharm model from CCDB! Maybe the model doesn't exist yet for this runnumber/timestamp?";
        }
      } else {
        bdtMCharm.initModel(bdt.localPath.value, bdt.enableOptimizations.value);
      }

      histos.add("BDT/hBDTScoreSignalFine", "hBDTScoreSignalFine", kTH1D, {axisBDTScoreFine});
      histos.add("BDT/hBDTScoreSignal", "hBDTScoreSignal", kTH1D, {axisBDTScore});
      histos.add("BDT/hBDTScoreVsXiccMassSignal", "hBDTScoreVsXiccMassSignal", kTH2D, {axisXiccMass, axisBDTScore});
      histos.add("BDT/hBDTScoreVsXiccPtSignal", "hBDTScoreVsXiccPtSignal", kTH2D, {axisPt, axisBDTScore});
      histos.add("BDT/h3dBDTScoreSignal", "h3dBDTScoreSignal", kTH3D, {axisPt, axisXiccMass, axisBDTScore});
      histos.add("BDT/hBDTScoreBackgroundFine", "hBDTScoreBackgroundFine", kTH1D, {axisBDTScoreFine});
      histos.add("BDT/hBDTScoreBackground", "hBDTScoreBackground", kTH1D, {axisBDTScore});
      histos.add("BDT/hBDTScoreVsXiccMassBackground", "hBDTScoreVsXiccMassBackground", kTH2D, {axisXiccMass, axisBDTScore});
      histos.add("BDT/hBDTScoreVsXiccPtBackground", "hBDTScoreVsXiccPtBackground", kTH2D, {axisPt, axisBDTScore});
      histos.add("BDT/h3dBDTScoreBackground", "h3dBDTScoreBackground", kTH3D, {axisPt, axisXiccMass, axisBDTScore});

      histos.add("BDT/hDCAXicDaughters", "hDCAXicDaughters", kTH2D, {{axisBDTScore, axisDcaDaughters}});
      histos.add("BDT/hDCAXiccDaughters", "hDCAXiccDaughters", kTH2D, {{axisBDTScore, axisDcaDaughters}});
      histos.add("BDT/hDCAxyXi", "hDCAxyXi", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hDCAzXi", "hDCAzXi", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hDCAxyXic", "hDCAxyXic", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hDCAzXic", "hDCAzXic", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hDCAxyXicc", "hDCAxyXicc", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hDCAzXicc", "hDCAzXicc", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hDecayRadiusXic", "hDecayRadiusXic", kTH2D, {{axisBDTScore, axisRadius}});
      histos.add("BDT/hDecayRadiusXicc", "hDecayRadiusXicc", kTH2D, {{axisBDTScore, axisRadius}});
      histos.add("BDT/hDecayDistanceFromPVXic", "hDecayDistanceFromPVXic", kTH2D, {{axisBDTScore, axisDecayLength}});
      histos.add("BDT/hProperLengthXic", "hProperLengthXic", kTH2D, {{axisBDTScore, axisDecayLength}});
      histos.add("BDT/hProperLengthXicc", "hProperLengthXicc", kTH2D, {{axisBDTScore, axisDecayLength}});
      histos.add("BDT/hPi1cDCAxy", "hPi1cDCAxy", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hPi1cDCAz", "hPi1cDCAz", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hPi2cDCAxy", "hPi2cDCAxy", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hPi2cDCAz", "hPi2cDCAz", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hPiccDCAxy", "hPiccDCAxy", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hPiccDCAz", "hPiccDCAz", kTH2D, {{axisBDTScore, axisDCA}});
      histos.add("BDT/hPi1cPt", "hPi1cPt", kTH2D, {{axisBDTScore, axisPt}});
      histos.add("BDT/hPi2cPt", "hPi2cPt", kTH2D, {{axisBDTScore, axisPt}});
      histos.add("BDT/hPiccPt", "hPiccPt", kTH2D, {{axisBDTScore, axisPt}});
    }

    histos.print();
  }

  template <typename TMCharmCands>
  void genericProcessXicc(TMCharmCands const& xiccCands)
  {
    for (const auto& xiccCand : xiccCands) {
      const int icfg = xiccCand.lutConfigId();
      histos.fill(HIST("hConfigId"), icfg);

      histos.fill(HIST("CandidateQA/hDCAXicDaughters"), xiccCand.xicDauDCA() * ToMicrons);
      histos.fill(HIST("CandidateQA/hDCAXiccDaughters"), xiccCand.xiccDauDCA() * ToMicrons);
      histos.fill(HIST("CandidateQA/hDCAxyXi"), std::fabs(xiccCand.xiDCAxy() * ToMicrons));
      histos.fill(HIST("CandidateQA/hDCAzXi"), std::fabs(xiccCand.xiDCAz() * ToMicrons));
      histos.fill(HIST("CandidateQA/hDCAxyXic"), std::fabs(xiccCand.xicDCAxy() * ToMicrons));
      histos.fill(HIST("CandidateQA/hDCAzXic"), std::fabs(xiccCand.xicDCAz() * ToMicrons));
      histos.fill(HIST("CandidateQA/hDCAxyXicc"), std::fabs(xiccCand.xiccDCAxy() * ToMicrons));
      histos.fill(HIST("CandidateQA/hDCAzXicc"), std::fabs(xiccCand.xiccDCAz() * ToMicrons));
      histos.fill(HIST("CandidateQA/hDecayRadiusXic"), xiccCand.xicDecayRadius2D() * ToMicrons);
      histos.fill(HIST("CandidateQA/hDecayRadiusXicc"), xiccCand.xiccDecayRadius2D() * ToMicrons);
      histos.fill(HIST("CandidateQA/hDecayDistanceFromPVXic"), xiccCand.xicDistanceFromPV() * ToMicrons);
      histos.fill(HIST("CandidateQA/hProperLengthXic"), xiccCand.xicProperLength() * ToMicrons);
      histos.fill(HIST("CandidateQA/hProperLengthXicc"), xiccCand.xiccProperLength() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPi1cDCAxy"), xiccCand.pi1cDCAxy() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPi1cDCAz"), xiccCand.pi1cDCAz() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPi2cDCAxy"), xiccCand.pi2cDCAxy() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPi2cDCAz"), xiccCand.pi2cDCAz() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPiccDCAxy"), xiccCand.piccDCAxy() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPiccDCAz"), xiccCand.piccDCAz() * ToMicrons);
      histos.fill(HIST("CandidateQA/hPi1cPt"), xiccCand.pi1cPt());
      histos.fill(HIST("CandidateQA/hPi2cPt"), xiccCand.pi2cPt());
      histos.fill(HIST("CandidateQA/hPiccPt"), xiccCand.piccPt());

      if (bdt.enableML) {
        std::vector<float> inputFeatures{
          xiccCand.xicDauDCA(),
          xiccCand.xiccDauDCA(),
          xiccCand.xiDCAxy(),
          xiccCand.xicDCAxy(),
          xiccCand.xiccDCAxy(),
          xiccCand.xiDCAz(),
          xiccCand.xicDCAz(),
          xiccCand.xiccDCAz(),
          xiccCand.pi1cDCAxy(),
          xiccCand.pi2cDCAxy(),
          xiccCand.piccDCAxy(),
          xiccCand.pi1cDCAz(),
          xiccCand.pi2cDCAz(),
          xiccCand.piccDCAz(),
          xiccCand.xicDecayRadius2D(),
          xiccCand.xiccDecayRadius2D(),
          xiccCand.xicProperLength(),
          xiccCand.xicDistanceFromPV(),
          xiccCand.xiccProperLength()};

        float* probabilityMCharm = bdtMCharm.evalModel(inputFeatures);
        const float bdtPredictedBackground = probabilityMCharm[0];
        const float bdtPredictedSignal = probabilityMCharm[1];

        histos.fill(HIST("BDT/hBDTScoreSignal"), bdtPredictedSignal);
        histos.fill(HIST("BDT/hBDTScoreSignalFine"), bdtPredictedSignal);
        histos.fill(HIST("BDT/hBDTScoreVsXiccMassSignal"), xiccCand.xiccMass(), bdtPredictedSignal);
        histos.fill(HIST("BDT/hBDTScoreVsXiccPtSignal"), xiccCand.xiccPt(), bdtPredictedSignal);
        histos.fill(HIST("BDT/h3dBDTScoreSignal"), xiccCand.xiccPt(), xiccCand.xiccMass(), bdtPredictedSignal);
        histos.fill(HIST("BDT/hBDTScoreBackground"), bdtPredictedBackground);
        histos.fill(HIST("BDT/hBDTScoreBackgroundFine"), bdtPredictedBackground);
        histos.fill(HIST("BDT/hBDTScoreVsXiccMassBackground"), xiccCand.xiccMass(), bdtPredictedBackground);
        histos.fill(HIST("BDT/hBDTScoreVsXiccPtBackground"), xiccCand.xiccPt(), bdtPredictedBackground);
        histos.fill(HIST("BDT/h3dBDTScoreBackground"), xiccCand.xiccPt(), xiccCand.xiccMass(), bdtPredictedBackground);

        histos.fill(HIST("BDT/hDCAXicDaughters"), bdtPredictedSignal, xiccCand.xicDauDCA() * ToMicrons);
        histos.fill(HIST("BDT/hDCAXiccDaughters"), bdtPredictedSignal, xiccCand.xiccDauDCA() * ToMicrons);
        histos.fill(HIST("BDT/hDCAxyXi"), bdtPredictedSignal, std::fabs(xiccCand.xiDCAxy() * ToMicrons));
        histos.fill(HIST("BDT/hDCAzXi"), bdtPredictedSignal, std::fabs(xiccCand.xiDCAz() * ToMicrons));
        histos.fill(HIST("BDT/hDCAxyXic"), bdtPredictedSignal, std::fabs(xiccCand.xicDCAxy() * ToMicrons));
        histos.fill(HIST("BDT/hDCAzXic"), bdtPredictedSignal, std::fabs(xiccCand.xicDCAz() * ToMicrons));
        histos.fill(HIST("BDT/hDCAxyXicc"), bdtPredictedSignal, std::fabs(xiccCand.xiccDCAxy() * ToMicrons));
        histos.fill(HIST("BDT/hDCAzXicc"), bdtPredictedSignal, std::fabs(xiccCand.xiccDCAz() * ToMicrons));
        histos.fill(HIST("BDT/hDecayRadiusXic"), bdtPredictedSignal, xiccCand.xicDecayRadius2D() * ToMicrons);
        histos.fill(HIST("BDT/hDecayRadiusXicc"), bdtPredictedSignal, xiccCand.xiccDecayRadius2D() * ToMicrons);
        histos.fill(HIST("BDT/hDecayDistanceFromPVXic"), bdtPredictedSignal, xiccCand.xicDistanceFromPV() * ToMicrons);
        histos.fill(HIST("BDT/hProperLengthXic"), bdtPredictedSignal, xiccCand.xicProperLength() * ToMicrons);
        histos.fill(HIST("BDT/hProperLengthXicc"), bdtPredictedSignal, xiccCand.xiccProperLength() * ToMicrons);
        histos.fill(HIST("BDT/hPi1cDCAxy"), bdtPredictedSignal, xiccCand.pi1cDCAxy() * ToMicrons);
        histos.fill(HIST("BDT/hPi1cDCAz"), bdtPredictedSignal, xiccCand.pi1cDCAz() * ToMicrons);
        histos.fill(HIST("BDT/hPi2cDCAxy"), bdtPredictedSignal, xiccCand.pi2cDCAxy() * ToMicrons);
        histos.fill(HIST("BDT/hPi2cDCAz"), bdtPredictedSignal, xiccCand.pi2cDCAz() * ToMicrons);
        histos.fill(HIST("BDT/hPiccDCAxy"), bdtPredictedSignal, xiccCand.piccDCAxy() * ToMicrons);
        histos.fill(HIST("BDT/hPiccDCAz"), bdtPredictedSignal, xiccCand.piccDCAz() * ToMicrons);
        histos.fill(HIST("BDT/hPi1cPt"), bdtPredictedSignal, xiccCand.pi1cPt());
        histos.fill(HIST("BDT/hPi2cPt"), bdtPredictedSignal, xiccCand.pi2cPt());
        histos.fill(HIST("BDT/hPiccPt"), bdtPredictedSignal, xiccCand.piccPt());
      }

      histos.fill(HIST("hMCharmBuilding"), 0);
      if (xiccCand.xicDauDCA() > xicMaxDauDCA) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 1);
      if (xiccCand.xiccDauDCA() > xiccMaxDauDCA) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 2);
      if (std::fabs(xiccCand.xiDCAxy()) < xiMinDCAxy) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 3);
      if (std::fabs(xiccCand.xiDCAz()) < xiMinDCAz) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 4);
      if (std::fabs(xiccCand.pi1cDCAxy()) < picMinDCAxy) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 5);
      if (std::fabs(xiccCand.pi1cDCAz()) < picMinDCAz) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 6);
      if (std::fabs(xiccCand.pi2cDCAxy()) < picMinDCAxy) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 7);
      if (std::fabs(xiccCand.pi2cDCAz()) < picMinDCAz) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 8);
      if (std::fabs(xiccCand.piccDCAxy()) < piccMinDCAxy) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 9);
      if (std::fabs(xiccCand.piccDCAz()) < piccMinDCAz) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 10);
      if (std::fabs(xiccCand.xicDCAxy()) < xicMinDCAxy) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 11);
      if (std::fabs(xiccCand.xicDCAz()) < xicMinDCAz) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 12);
      if (std::fabs(xiccCand.xiccDCAxy()) > xiccMaxDCAxy) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 13);
      if (std::fabs(xiccCand.xiccDCAz()) > xiccMaxDCAz) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 14);
      if (xiccCand.xicDecayRadius2D() < xicMinRadius) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 15);
      if (xiccCand.xiccDecayRadius2D() < xiccMinRadius) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 16);
      if (xiccCand.xicProperLength() < xicMinProperLength) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 17);
      if (xiccCand.xicProperLength() > xicMaxProperLength) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 18);
      if (xiccCand.xiccProperLength() < xiccMinProperLength) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 19);
      if (xiccCand.xiccProperLength() > xiccMaxProperLength) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 20);
      if (xiccCand.xicDistanceFromPV() < xicMinDecayDistanceFromPV) {
        continue;
      }

      histos.fill(HIST("hMCharmBuilding"), 21);
      histos.fill(HIST("SelectionQA/hDCAXicDaughters"), xiccCand.xicDauDCA() * ToMicrons);
      histos.fill(HIST("SelectionQA/hDCAXiccDaughters"), xiccCand.xiccDauDCA() * ToMicrons);
      histos.fill(HIST("SelectionQA/hDCAxyXi"), std::fabs(xiccCand.xiDCAxy() * ToMicrons));
      histos.fill(HIST("SelectionQA/hDCAzXi"), std::fabs(xiccCand.xiDCAz() * ToMicrons));
      histos.fill(HIST("SelectionQA/hDCAxyXic"), std::fabs(xiccCand.xicDCAxy() * ToMicrons));
      histos.fill(HIST("SelectionQA/hDCAzXic"), std::fabs(xiccCand.xicDCAz() * ToMicrons));
      histos.fill(HIST("SelectionQA/hDCAxyXicc"), std::fabs(xiccCand.xiccDCAxy() * ToMicrons));
      histos.fill(HIST("SelectionQA/hDCAzXicc"), std::fabs(xiccCand.xiccDCAz() * ToMicrons));
      histos.fill(HIST("SelectionQA/hDecayRadiusXic"), xiccCand.xicDecayRadius2D() * ToMicrons);
      histos.fill(HIST("SelectionQA/hDecayRadiusXicc"), xiccCand.xiccDecayRadius2D() * ToMicrons);
      histos.fill(HIST("SelectionQA/hDecayDistanceFromPVXic"), xiccCand.xicDistanceFromPV() * ToMicrons);
      histos.fill(HIST("SelectionQA/hProperLengthXic"), xiccCand.xicProperLength() * ToMicrons);
      histos.fill(HIST("SelectionQA/hProperLengthXicc"), xiccCand.xiccProperLength() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPi1cDCAxy"), xiccCand.pi1cDCAxy() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPi1cDCAz"), xiccCand.pi1cDCAz() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPi2cDCAxy"), xiccCand.pi2cDCAxy() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPi2cDCAz"), xiccCand.pi2cDCAz() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPiccDCAxy"), xiccCand.piccDCAxy() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPiccDCAz"), xiccCand.piccDCAz() * ToMicrons);
      histos.fill(HIST("SelectionQA/hPi1cPt"), xiccCand.pi1cPt());
      histos.fill(HIST("SelectionQA/hPi2cPt"), xiccCand.pi2cPt());
      histos.fill(HIST("SelectionQA/hPiccPt"), xiccCand.piccPt());
      histos.fill(HIST("hXiccMass"), xiccCand.xiccMass());
      histos.fill(HIST("hXicMass"), xiccCand.xicMass());
      histos.fill(HIST("hXiccPt"), xiccCand.xiccPt());
      histos.fill(HIST("hXicPt"), xiccCand.xicPt());
      histos.fill(HIST("h3dXicc"), xiccCand.xiccPt(), xiccCand.xiccEta(), xiccCand.xiccMass());
    }
  }

  void processXicc(soa::Filtered<aod::MCharmCores> const& multiCharmTracks)
  {
    genericProcessXicc(multiCharmTracks);
  }

  PROCESS_SWITCH(Alice3Multicharm, processXicc, "find Xicc baryons", true);
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
  return WorkflowSpec{
    adaptAnalysisTask<Alice3Multicharm>(cfgc)};
}