diff --git a/PWGHF/HFC/TableProducer/CMakeLists.txt b/PWGHF/HFC/TableProducer/CMakeLists.txt index fbf7a06c6dc..be466d7e86e 100644 --- a/PWGHF/HFC/TableProducer/CMakeLists.txt +++ b/PWGHF/HFC/TableProducer/CMakeLists.txt @@ -83,3 +83,8 @@ o2physics_add_dpl_workflow(producer-charm-hadrons-track-femto-dream SOURCES producerCharmHadronsTrackFemtoDream.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) + +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) diff --git a/PWGHF/HFC/TableProducer/producerCharmHadronsV0FemtoDream.cxx b/PWGHF/HFC/TableProducer/producerCharmHadronsV0FemtoDream.cxx new file mode 100644 index 00000000000..61eb7181061 --- /dev/null +++ b/PWGHF/HFC/TableProducer/producerCharmHadronsV0FemtoDream.cxx @@ -0,0 +1,1290 @@ +// Copyright 2019-2025 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 producerCharmHadronsV0FemtoDream.cxx +/// \brief Tasks that produces the V0 and CharmHadrons tables used for the pairing +/// \author Biao Zhang, Heidelberg University, biao.zhang@cern.ch + +#include "PWGCF/DataModel/FemtoDerived.h" +#include "PWGCF/FemtoDream/Core/femtoDreamSelection.h" +#include "PWGCF/FemtoDream/Core/femtoDreamTrackSelection.h" +#include "PWGCF/FemtoDream/Core/femtoDreamUtils.h" +#include "PWGCF/FemtoDream/Core/femtoDreamV0SelectionK0Short.h" +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/Core/DecayChannels.h" +#include "PWGHF/Core/HfMlResponseD0ToKPi.h" +#include "PWGHF/Core/HfMlResponseDplusToPiKPi.h" +#include "PWGHF/Core/HfMlResponseDstarToD0Pi.h" +#include "PWGHF/Core/HfMlResponseLcToPKPi.h" +#include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" +#include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::analysis; +using namespace o2::framework::expressions; +using namespace o2::analysis::femtoDream; +using namespace o2::hf_evsel; +using namespace o2::hf_centrality; + +// event types +enum Event : uint8_t { + All = 0, + RejEveSel, + RejNoV0sAndCharm, + V0Selected, + CharmSelected, + PairSelected +}; + +// ml modes +enum MlMode : uint8_t { + NoMl = 0, + FillMlFromSelector, + FillMlFromNewBDT +}; + +// decay channels +enum DecayChannel { DplusToPiKPi = 0, + LcToPKPi, + D0ToPiK, + DstarToD0Pi +}; + +enum V0Channel { K0S = 0, + Lambda +}; + +struct HfProducerCharmHadronsV0FemtoDream { + + Produces outputCollision; + Produces rowMasks; + Produces rowCandCharm3Prong; + Produces rowCandCharm2Prong; + Produces rowCandCharmDstar; + Produces rowCandMcCharmHad; + Produces rowCandCharmHadGen; + Produces outputPartsIndex; + Produces outputPartsTime; + Produces outputMcCollision; + Produces outputCollsMcLabels; + Produces outputParts; + Produces outputPartsMc; + Produces outputDebugParts; + Produces outputPartsMcLabels; + Produces outputDebugPartsMc; + Produces outputPartsExtMcLabels; + + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbPathLut{"ccdbPathLut", "GLO/Param/MatLUT", "Path for LUT parametrization"}; + Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "Path of the grp file (Run 2)"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"EventFiltering/PWGHF/BDTLc"}, "Paths of models on CCDB"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"ModelHandler_onnx_LcToPKPi.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + + // Configurable isForceGRP{"isForceGRP", false, "Set true if the magnetic field configuration is not available in the usual CCDB directory (e.g. for Run 2 converted data or unanchorad Monte Carlo)"}; + + Configurable isDebug{"isDebug", true, "Enable Debug tables"}; + Configurable isRun3{"isRun3", true, "Running on Run3 or pilot"}; + Configurable selectionFlagV0{"selectionFlagV0", 0, "Activate filling of V0: 0 for K0S, 1 for Lambda"}; + /// Charm hadron table + Configurable selectionFlagCharmHadron{"selectionFlagCharmHadron", 1, "Selection Flag for Charm Hadron: 1 for Lc, 7 for Dplus (Topologic and PID cuts)"}; + Configurable useCent{"useCent", false, "Enable centrality for Charm Hadron"}; + + Configurable v0PDGCode{"v0PDGCode", 2212, "PDG code of the selected track for Monte Carlo truth"}; + Configurable trkPIDnSigmaOffsetTPC{"trkPIDnSigmaOffsetTPC", 0., "Offset for TPC nSigma because of bad calibration"}; // set to zero for run3 or so + Configurable trkPIDnSigmaOffsetTOF{"trkPIDnSigmaOffsetTOF", 0., "Offset for TOF nSigma because of bad calibration"}; + Configurable trkRejectNotPropagated{"trkRejectNotPropagated", false, "True: reject not propagated tracks"}; + + struct : o2::framework::ConfigurableGroup { + Configurable> confLambdaSign{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0Sign, "confLambda"), std::vector{-1, 1}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0Sign, "V0 selection: ")}; + Configurable> confLambdaPtMin{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0pTMin, "confLambda"), std::vector{0.3f, 0.4f, 0.5f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0pTMin, "V0 selection: ")}; + Configurable> confLambdaPtMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0pTMax, "confLambda"), std::vector{3.3f, 3.4f, 3.5f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0pTMax, "V0 selection: ")}; + Configurable> confLambdaEtaMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0etaMax, "confLambda"), std::vector{0.8f, 0.7f, 0.9f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0etaMax, "V0 selection: ")}; + Configurable> confLambdaDCADaughMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0DCADaughMax, "confLambda"), std::vector{1.2f, 1.5f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0DCADaughMax, "V0 selection: ")}; + Configurable> confLambdaCPAMin{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0CPAMin, "confLambda"), std::vector{0.99f, 0.995f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0CPAMin, "V0 selection: ")}; + Configurable> confLambdaTranRadMin{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0TranRadMin, "confLambda"), std::vector{0.2f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0TranRadMin, "V0 selection: ")}; + Configurable> confLambdaTranRadMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0TranRadMax, "confLambda"), std::vector{100.f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0TranRadMax, "V0 selection: ")}; + Configurable> confLambdaDecVtxMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0DecVtxMax, "confLambda"), std::vector{100.f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0DecVtxMax, "V0 selection: ")}; + + Configurable confLambdaInvMassLowLimit{"confLambdaInvMassLowLimit", 1.05, "Lower limit of the V0 invariant mass"}; + Configurable confLambdaInvMassUpLimit{"confLambdaInvMassUpLimit", 1.30, "Upper limit of the V0 invariant mass"}; + Configurable confLambdaRejectKaons{"confLambdaRejectKaons", false, "Switch to reject kaons"}; + Configurable confLambdaRejectLambdas{"confLambdaRejectLambdas", false, "Switch to reject lambdas (if mother is kaon)"}; + Configurable confLambdaInvKaonMassLowLimit{"confLambdaInvKaonMassLowLimit", 0.48, "Lower limit of the V0 invariant mass for Kaon rejection"}; + Configurable confLambdaInvKaonMassUpLimit{"confLambdaInvKaonMassUpLimit", 0.515, "Upper limit of the V0 invariant mass for Kaon rejection"}; + + Configurable> confLambdaChildSign{"confLambdaChildSign", std::vector{-1, 1}, "V0 Child sel: Charge"}; + Configurable> confLambdaChildEtaMax{"confLambdaChildEtaMax", std::vector{0.8f}, "V0 Child sel: max eta"}; + Configurable> confLambdaChildTPCnClsMin{"confLambdaChildTPCnClsMin", std::vector{80.f, 70.f, 60.f}, "V0 Child sel: Min. nCls TPC"}; + Configurable> confLambdaChildDCAMin{"confLambdaChildDCAMin", std::vector{0.05f, 0.06f}, "V0 Child sel: Max. DCA Daugh to PV (cm)"}; + Configurable> confLambdaChildPIDnSigmaMax{"confLambdaChildPIDnSigmaMax", std::vector{5.f, 4.f}, "V0 Child sel: Max. PID nSigma TPC"}; + Configurable> confLambdaChildPIDspecies{"confLambdaChildPIDspecies", std::vector{o2::track::PID::Pion, o2::track::PID::Proton}, "V0 Child sel: Particles species for PID"}; + + // cuts and object for v0 2 + Configurable> confK0shortSign{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0Sign, "confK0short"), std::vector{-1, 1}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0Sign, "V0 selection: ")}; + Configurable> confK0shortPtMin{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0pTMin, "confK0short"), std::vector{0.3f, 0.4f, 0.5f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0pTMin, "V0 selection: ")}; + Configurable> confK0shortPtMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0pTMax, "confK0short"), std::vector{3.3f, 3.4f, 3.5f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0pTMax, "V0 selection: ")}; + Configurable> confK0shortEtaMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0etaMax, "confK0short"), std::vector{0.8f, 0.7f, 0.9f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0etaMax, "V0 selection: ")}; + Configurable> confK0shortDCADaughMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0DCADaughMax, "confK0short"), std::vector{1.2f, 1.5f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0DCADaughMax, "V0 selection: ")}; + Configurable> confK0shortCPAMin{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0CPAMin, "confK0short"), std::vector{0.99f, 0.995f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0CPAMin, "V0 selection: ")}; + Configurable> confK0shortTranRadMin{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0TranRadMin, "confK0short"), std::vector{0.2f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0TranRadMin, "V0 selection: ")}; + Configurable> confK0shortTranRadMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0TranRadMax, "confK0short"), std::vector{100.f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0TranRadMax, "V0 selection: ")}; + Configurable> confK0shortDecVtxMax{FemtoDreamV0Selection::getSelectionName(femto_dream_v0_selection::kV0DecVtxMax, "confK0short"), std::vector{100.f}, FemtoDreamV0Selection::getSelectionHelper(femto_dream_v0_selection::kV0DecVtxMax, "V0 selection: ")}; + + Configurable confK0shortInvMassLowLimit{"confK0shortInvMassLowLimit", 0.45, "Lower limit of the V0 invariant mass"}; + Configurable confK0shortInvMassUpLimit{"confK0shortInvMassUpLimit", 0.55, "Upper limit of the V0 invariant mass"}; + Configurable confK0shortRejectKaons{"confK0shortRejectKaons", false, "Switch to reject kaons"}; + Configurable confK0shortRejectLambdas{"confK0shortRejectLambdas", false, "Switch to reject lambdas (if mother is kaon)"}; + Configurable confK0shortInvKaonMassLowLimit{"confK0shortInvKaonMassLowLimit", 0.48, "Lower limit of the V0 invariant mass for Kaon rejection"}; + Configurable confK0shortInvKaonMassUpLimit{"confK0shortInvKaonMassUpLimit", 0.515, "Upper limit of the V0 invariant mass for Kaon rejection"}; + + Configurable> confK0shortChildSign{"confK0shortChildSign", std::vector{-1, 1}, "V0 Child sel: Charge"}; + Configurable> confK0shortChildEtaMax{"confK0shortChildEtaMax", std::vector{0.8f}, "V0 Child sel: max eta"}; + Configurable> confK0shortChildTPCnClsMin{"confK0shortChildTPCnClsMin", std::vector{80.f, 70.f, 60.f}, "V0 Child sel: Min. nCls TPC"}; + Configurable> confK0shortChildDCAMin{"confK0shortChildDCAMin", std::vector{0.05f, 0.06f}, "V0 Child sel: Max. DCA Daugh to PV (cm)"}; + Configurable> confK0shortChildPIDnSigmaMax{"confK0shortChildPIDnSigmaMax", std::vector{5.f, 4.f}, "V0 Child sel: Max. PID nSigma TPC"}; + Configurable> confK0shortChildPIDspecies{"confK0shortChildPIDspecies", std::vector{o2::track::PID::Pion, o2::track::PID::Proton}, "V0 Child sel: Particles species for PID"}; + } V0Sel; + + // ML inference + Configurable applyMlMode{"applyMlMode", 1, "None: 0, BDT model from candidate selector: 1, New BDT model on Top of candidate selector: 2"}; + Configurable> binsPtMl{"binsPtMl", std::vector{hf_cuts_ml::vecBinsPt}, "pT bin limits for ML application"}; + Configurable> cutDirMl{"cutDirMl", std::vector{hf_cuts_ml::vecCutDir}, "Whether to reject score values greater or smaller than the threshold"}; + Configurable> cutsMl{"cutsMl", {hf_cuts_ml::Cuts[0], hf_cuts_ml::NBinsPt, hf_cuts_ml::NCutScores, hf_cuts_ml::labelsPt, hf_cuts_ml::labelsCutScore}, "ML selections per pT bin"}; + Configurable nClassesMl{"nClassesMl", static_cast(hf_cuts_ml::NCutScores), "Number of classes in ML model"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; + + o2::analysis::HfMlResponseLcToPKPi hfMlResponseLc; + o2::analysis::HfMlResponseDplusToPiKPi hfMlResponseDplus; + o2::analysis::HfMlResponseD0ToKPi hfMlResponseD0; + o2::analysis::HfMlResponseDstarToD0Pi hfMlResponseDstar; + + std::vector outputMlD0; + std::vector outputMlD0bar; + std::vector outputMlDstar; + std::vector outputMlDplus; + std::vector outputMlPKPi; + std::vector outputMlPiKP; + o2::ccdb::CcdbApi ccdbApi; + o2::hf_evsel::HfEventSelection hfEvSel; + Service ccdb; /// Accessing the CCDB + o2::base::MatLayerCylSet* lut{}; + // if (doPvRefit){ lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(ccdbPathLut));} //! may be it useful, will check later + + float magField{}; + int runNumber{}; + using CandidateD0 = soa::Join; + using CandidateD0Mc = soa::Join; + using CandidateDstar = soa::Join; + using CandidateDstarMc = soa::Join; + using CandidateDplus = soa::Join; + using CandidateDplusMc = soa::Join; + using CandidateLc = soa::Join; + using CandidateLcMc = soa::Join; + using FemtoHFTracks = soa::Join; + using FemtoHFMcTracks = soa::Join; + using FemtoFullCollision = soa::Join::iterator; + using FemtoFullCollisionMc = soa::Join::iterator; + using FemtoFullMcgenCollisions = soa::Join; + using FemtoFullMcgenCollision = FemtoFullMcgenCollisions::iterator; + + using GeneratedMc = soa::Filtered>; + + Filter filterSelectCandidateD0 = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagCharmHadron || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagCharmHadron); + Filter filterSelectCandidateDstar = aod::hf_sel_candidate_dstar::isSelDstarToD0Pi == true; + Filter filterSelectCandidateDplus = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagCharmHadron; + Filter filterSelectCandidateLc = (aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlagCharmHadron || aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlagCharmHadron); + + HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry qaRegistryV0{"V0", {}, OutputObjHandlingPolicy::AnalysisObject}; + FemtoDreamV0Selection K0SCuts; + FemtoDreamV0Selection LambdaCuts; + void init(InitContext&) + { + std::array processes = {doprocessDataDplusToPiKPi, doprocessMcDplusToPiKPi, doprocessDataDplusToPiKPiWithML, doprocessMcDplusToPiKPiWithML, doprocessMcDplusToPiKPiGen, + doprocessDataLcToPKPi, doprocessMcLcToPKPi, doprocessDataLcToPKPiWithML, doprocessMcLcToPKPiWithML, doprocessMcLcToPKPiGen, doprocessDataD0ToPiK, doprocessMcD0ToPiK, doprocessDataD0ToPiKWithML, doprocessMcD0ToPiKWithML, doprocessMcD0ToPiKGen, doprocessDataDstarToD0Pi, doprocessMcDstarToD0Pi, doprocessDataDstarToD0PiWithML, doprocessMcDstarToD0PiWithML, doprocessMcDstarToD0PiGen}; + if (std::accumulate(processes.begin(), processes.end(), 0) != 1) { + LOGP(fatal, "One and only one process function must be enabled at a time."); + } + + int const cutBits = 8 * sizeof(o2::aod::femtodreamparticle::cutContainerType); + qaRegistryV0.add("AnalysisQA/CutCounter", "; Bit; Counter", kTH1F, {{cutBits + 1, -0.5, cutBits + 0.5}}); + + // event QA histograms + constexpr int kEventTypes = PairSelected + 1; + std::string labels[kEventTypes]; + labels[Event::All] = "All events"; + labels[Event::RejEveSel] = "rejected by event selection"; + labels[Event::RejNoV0sAndCharm] = "rejected by no V0s and charm"; + labels[Event::V0Selected] = "with V0 "; + labels[Event::CharmSelected] = "with charm hadrons "; + labels[Event::PairSelected] = "with pairs"; + + static const AxisSpec axisEvents = {kEventTypes, 0.5, kEventTypes + 0.5, ""}; + qaRegistry.add("hEventQA", "Events;;entries", HistType::kTH1F, {axisEvents}); + for (int iBin = 0; iBin < kEventTypes; iBin++) { + qaRegistry.get(HIST("hEventQA"))->GetXaxis()->SetBinLabel(iBin + 1, labels[iBin].data()); + } + + if (selectionFlagV0 == V0Channel::K0S) { + K0SCuts.setSelection(V0Sel.confK0shortSign, femto_dream_v0_selection::kV0Sign, femtoDreamSelection::kEqual); + K0SCuts.setSelection(V0Sel.confK0shortPtMin, femto_dream_v0_selection::kV0pTMin, femtoDreamSelection::kLowerLimit); + K0SCuts.setSelection(V0Sel.confK0shortPtMax, femto_dream_v0_selection::kV0pTMax, femtoDreamSelection::kUpperLimit); + K0SCuts.setSelection(V0Sel.confK0shortEtaMax, femto_dream_v0_selection::kV0etaMax, femtoDreamSelection::kAbsUpperLimit); + K0SCuts.setSelection(V0Sel.confK0shortDCADaughMax, femto_dream_v0_selection::kV0DCADaughMax, femtoDreamSelection::kUpperLimit); + K0SCuts.setSelection(V0Sel.confK0shortCPAMin, femto_dream_v0_selection::kV0CPAMin, femtoDreamSelection::kLowerLimit); + K0SCuts.setSelection(V0Sel.confK0shortTranRadMin, femto_dream_v0_selection::kV0TranRadMin, femtoDreamSelection::kLowerLimit); + K0SCuts.setSelection(V0Sel.confK0shortTranRadMax, femto_dream_v0_selection::kV0TranRadMax, femtoDreamSelection::kUpperLimit); + K0SCuts.setSelection(V0Sel.confK0shortDecVtxMax, femto_dream_v0_selection::kV0DecVtxMax, femtoDreamSelection::kUpperLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confK0shortChildSign, femtoDreamTrackSelection::kSign, femtoDreamSelection::kEqual); + K0SCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confK0shortChildEtaMax, femtoDreamTrackSelection::kEtaMax, femtoDreamSelection::kAbsUpperLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confK0shortChildTPCnClsMin, femtoDreamTrackSelection::kTPCnClsMin, femtoDreamSelection::kLowerLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confK0shortChildDCAMin, femtoDreamTrackSelection::kDCAMin, femtoDreamSelection::kAbsLowerLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confK0shortChildPIDnSigmaMax, femtoDreamTrackSelection::kPIDnSigmaMax, femtoDreamSelection::kAbsUpperLimit); + + K0SCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confK0shortChildSign, femtoDreamTrackSelection::kSign, femtoDreamSelection::kEqual); + K0SCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confK0shortChildEtaMax, femtoDreamTrackSelection::kEtaMax, femtoDreamSelection::kAbsUpperLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confK0shortChildTPCnClsMin, femtoDreamTrackSelection::kTPCnClsMin, femtoDreamSelection::kLowerLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confK0shortChildDCAMin, femtoDreamTrackSelection::kDCAMin, femtoDreamSelection::kAbsLowerLimit); + K0SCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confK0shortChildPIDnSigmaMax, femtoDreamTrackSelection::kPIDnSigmaMax, femtoDreamSelection::kAbsUpperLimit); + K0SCuts.setChildPIDSpecies(femto_dream_v0_selection::kPosTrack, V0Sel.confK0shortChildPIDspecies); + K0SCuts.setChildPIDSpecies(femto_dream_v0_selection::kNegTrack, V0Sel.confK0shortChildPIDspecies); + K0SCuts.init(&qaRegistryV0, &qaRegistryV0); + K0SCuts.setInvMassLimits(V0Sel.confK0shortInvMassLowLimit, V0Sel.confK0shortInvMassUpLimit); + K0SCuts.setIsMother(false); + + K0SCuts.setChildRejectNotPropagatedTracks(femto_dream_v0_selection::kPosTrack, trkRejectNotPropagated); + K0SCuts.setChildRejectNotPropagatedTracks(femto_dream_v0_selection::kNegTrack, trkRejectNotPropagated); + + K0SCuts.setnSigmaPIDOffsetTPC(trkPIDnSigmaOffsetTPC); + K0SCuts.setChildnSigmaPIDOffset(femto_dream_v0_selection::kPosTrack, trkPIDnSigmaOffsetTPC, trkPIDnSigmaOffsetTOF); + K0SCuts.setChildnSigmaPIDOffset(femto_dream_v0_selection::kNegTrack, trkPIDnSigmaOffsetTPC, trkPIDnSigmaOffsetTOF); + + K0SCuts.setRejectLambda(V0Sel.confK0shortRejectLambdas); + K0SCuts.setKaonInvMassLimits(V0Sel.confK0shortInvKaonMassLowLimit, V0Sel.confK0shortInvKaonMassUpLimit); + } + + if (selectionFlagV0 == V0Channel::Lambda) { + LambdaCuts.setSelection(V0Sel.confLambdaSign, femto_dream_v0_selection::kV0Sign, femtoDreamSelection::kEqual); + LambdaCuts.setSelection(V0Sel.confLambdaPtMin, femto_dream_v0_selection::kV0pTMin, femtoDreamSelection::kLowerLimit); + LambdaCuts.setSelection(V0Sel.confLambdaPtMax, femto_dream_v0_selection::kV0pTMax, femtoDreamSelection::kUpperLimit); + LambdaCuts.setSelection(V0Sel.confLambdaEtaMax, femto_dream_v0_selection::kV0etaMax, femtoDreamSelection::kAbsUpperLimit); + LambdaCuts.setSelection(V0Sel.confLambdaDCADaughMax, femto_dream_v0_selection::kV0DCADaughMax, femtoDreamSelection::kUpperLimit); + LambdaCuts.setSelection(V0Sel.confLambdaCPAMin, femto_dream_v0_selection::kV0CPAMin, femtoDreamSelection::kLowerLimit); + LambdaCuts.setSelection(V0Sel.confLambdaTranRadMin, femto_dream_v0_selection::kV0TranRadMin, femtoDreamSelection::kLowerLimit); + LambdaCuts.setSelection(V0Sel.confLambdaTranRadMax, femto_dream_v0_selection::kV0TranRadMax, femtoDreamSelection::kUpperLimit); + LambdaCuts.setSelection(V0Sel.confLambdaDecVtxMax, femto_dream_v0_selection::kV0DecVtxMax, femtoDreamSelection::kUpperLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confLambdaChildSign, femtoDreamTrackSelection::kSign, femtoDreamSelection::kEqual); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confLambdaChildEtaMax, femtoDreamTrackSelection::kEtaMax, femtoDreamSelection::kAbsUpperLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confLambdaChildTPCnClsMin, femtoDreamTrackSelection::kTPCnClsMin, femtoDreamSelection::kLowerLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confLambdaChildDCAMin, femtoDreamTrackSelection::kDCAMin, femtoDreamSelection::kAbsLowerLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kPosTrack, V0Sel.confLambdaChildPIDnSigmaMax, femtoDreamTrackSelection::kPIDnSigmaMax, femtoDreamSelection::kAbsUpperLimit); + + LambdaCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confLambdaChildSign, femtoDreamTrackSelection::kSign, femtoDreamSelection::kEqual); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confLambdaChildEtaMax, femtoDreamTrackSelection::kEtaMax, femtoDreamSelection::kAbsUpperLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confLambdaChildTPCnClsMin, femtoDreamTrackSelection::kTPCnClsMin, femtoDreamSelection::kLowerLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confLambdaChildDCAMin, femtoDreamTrackSelection::kDCAMin, femtoDreamSelection::kAbsLowerLimit); + LambdaCuts.setChildCuts(femto_dream_v0_selection::kNegTrack, V0Sel.confLambdaChildPIDnSigmaMax, femtoDreamTrackSelection::kPIDnSigmaMax, femtoDreamSelection::kAbsUpperLimit); + LambdaCuts.setChildPIDSpecies(femto_dream_v0_selection::kPosTrack, V0Sel.confLambdaChildPIDspecies); + LambdaCuts.setChildPIDSpecies(femto_dream_v0_selection::kNegTrack, V0Sel.confLambdaChildPIDspecies); + LambdaCuts.init(&qaRegistryV0, &qaRegistryV0); + LambdaCuts.setInvMassLimits(V0Sel.confLambdaInvMassLowLimit, V0Sel.confLambdaInvMassUpLimit); + LambdaCuts.setIsMother(true); + + LambdaCuts.setChildRejectNotPropagatedTracks(femto_dream_v0_selection::kPosTrack, trkRejectNotPropagated); + LambdaCuts.setChildRejectNotPropagatedTracks(femto_dream_v0_selection::kNegTrack, trkRejectNotPropagated); + + LambdaCuts.setnSigmaPIDOffsetTPC(trkPIDnSigmaOffsetTPC); + LambdaCuts.setChildnSigmaPIDOffset(femto_dream_v0_selection::kPosTrack, trkPIDnSigmaOffsetTPC, trkPIDnSigmaOffsetTOF); + LambdaCuts.setChildnSigmaPIDOffset(femto_dream_v0_selection::kNegTrack, trkPIDnSigmaOffsetTPC, trkPIDnSigmaOffsetTOF); + + if (V0Sel.confLambdaRejectKaons) { + LambdaCuts.setKaonInvMassLimits(V0Sel.confLambdaInvKaonMassLowLimit, V0Sel.confLambdaInvKaonMassUpLimit); + } + } + + runNumber = 0; + magField = 0.0; + /// Initializing CCDB + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + + hfEvSel.addHistograms(qaRegistry); // collision monitoring + + int64_t const now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); + + bool useLcMl = doprocessDataLcToPKPiWithML || doprocessMcLcToPKPiWithML; + bool useDplusMl = doprocessDataDplusToPiKPiWithML || doprocessMcDplusToPiKPiWithML; + bool useD0Ml = doprocessDataD0ToPiKWithML || doprocessMcD0ToPiKWithML; + bool useDstarMl = doprocessDataDstarToD0PiWithML || doprocessMcDstarToD0PiWithML; + + if (applyMlMode == FillMlFromNewBDT) { + if (useLcMl) { + hfMlResponseLc.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + hfMlResponseLc.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponseLc.init(); + } + if (useDplusMl) { + hfMlResponseDplus.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + hfMlResponseDplus.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponseDplus.init(); + } + if (useD0Ml) { + hfMlResponseD0.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + hfMlResponseD0.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponseD0.init(); + } + if (useDstarMl) { + hfMlResponseDstar.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + hfMlResponseDstar.cacheInputFeaturesIndices(namesInputFeatures); + hfMlResponseDstar.init(); + } + + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + if (useLcMl) { + hfMlResponseLc.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } + if (useDplusMl) { + hfMlResponseDplus.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } + if (useD0Ml) { + hfMlResponseD0.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } + if (useDstarMl) { + hfMlResponseDstar.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + } + + } else { + if (useLcMl) { + hfMlResponseLc.setModelPathsLocal(onnxFileNames); + } + if (useDplusMl) { + hfMlResponseDplus.setModelPathsLocal(onnxFileNames); + } + if (useD0Ml) { + hfMlResponseD0.setModelPathsLocal(onnxFileNames); + } + if (useDstarMl) { + hfMlResponseDstar.setModelPathsLocal(onnxFileNames); + } + } + } + } + + /// Function to retrieve the nominal magnetic field in kG (0.1T) and convert it directly to T + void getMagneticFieldTesla(const aod::BCsWithTimestamps::iterator& bc) + { + initCCDB(bc, runNumber, ccdb, !isRun3 ? ccdbPathGrp : ccdbPathGrpMag, lut, !isRun3); + } + + template + void fillDebugParticle(ParticleType const& particle) + { + if constexpr (isTrackOrV0) { + outputDebugParts(particle.sign(), + (uint8_t)particle.tpcNClsFound(), + particle.tpcNClsFindable(), + (uint8_t)particle.tpcNClsCrossedRows(), + particle.tpcNClsShared(), + particle.tpcInnerParam(), + particle.itsNCls(), + particle.itsNClsInnerBarrel(), + particle.dcaXY(), + particle.dcaZ(), + particle.tpcSignal(), + -999., + particle.tpcNSigmaPi(), + particle.tpcNSigmaKa(), + particle.tpcNSigmaPr(), + particle.tpcNSigmaDe(), + -999., + -999., + -999., + particle.tofNSigmaPi(), + particle.tofNSigmaKa(), + particle.tofNSigmaPr(), + particle.tofNSigmaDe(), + -999., -999., -999., -999., + -999., -999., -999., -999., + -999., -999., -999., -999., + -999., -999., -999., -999., + -999., -999., -999., -999., + -999., -999., -999.); + } else { + outputDebugParts(-999., // sign + -999., -999., -999., -999., -999., -999., -999., -999., -999., // track properties (DCA, NCls, crossed rows, etc.) + -999., -999., -999., -999., -999., -999., -999., -999., // TPC PID (TPC signal + particle hypothesis) + -999., -999., -999., -999., -999., -999., -999., // TOF PID + -999., -999., -999., -999., -999., -999., -999., -999., // ITS PID + particle.dcaV0daughters(), + particle.v0radius(), + particle.x(), + particle.y(), + particle.z(), + particle.mK0Short(), + -999., -999., -999., -999., -999., -999., -999.); + } + } + + template + void fillMcParticle(CollisionType const& col, ParticleType const& particle, o2::aod::femtodreamparticle::ParticleType fdparttype) + { + if (particle.has_mcParticle()) { + // get corresponding MC particle and its info + auto particleMc = particle.mcParticle(); + auto pdgCode = particleMc.pdgCode(); + int particleOrigin = 99; + int pdgCodeMother = -1; + constexpr int GenFromTransport = -1; // -1 if a particle produced during transport + // get list of mothers, but it could be empty (for example in case of injected light nuclei) + auto motherparticlesMc = particleMc.template mothers_as(); + // check pdg code + // if this fails, the particle is a fake + if (std::abs(pdgCode) == std::abs(v0PDGCode.value)) { + // check first if particle is from pile up + // check if the collision associated with the particle is the same as the analyzed collision by checking their Ids + if ((col.has_mcCollision() && (particleMc.mcCollisionId() != col.mcCollisionId())) || !col.has_mcCollision()) { + particleOrigin = aod::femtodreamMCparticle::ParticleOriginMCTruth::kWrongCollision; + // check if particle is primary + } else if (particleMc.isPhysicalPrimary()) { + particleOrigin = aod::femtodreamMCparticle::ParticleOriginMCTruth::kPrimary; + // check if particle is secondary + // particle is from a decay -> getProcess() == 4 + // particle is generated during transport -> getGenStatusCode() == -1 + // list of mothers is not empty + } else if (particleMc.getProcess() == TMCProcess::kPDecay && particleMc.getGenStatusCode() == GenFromTransport && !motherparticlesMc.empty()) { + // get direct mother + auto motherparticleMc = motherparticlesMc.front(); + pdgCodeMother = motherparticleMc.pdgCode(); + particleOrigin = checkDaughterType(fdparttype, motherparticleMc.pdgCode()); + // check if particle is material + // particle is from inelastic hadronic interaction -> getProcess() == 23 + // particle is generated during transport -> getGenStatusCode() == -1 + } else if (particleMc.getProcess() == TMCProcess::kPHInhelastic && particleMc.getGenStatusCode() == GenFromTransport) { + particleOrigin = aod::femtodreamMCparticle::ParticleOriginMCTruth::kMaterial; + // cross check to see if we missed a case + } else { + particleOrigin = aod::femtodreamMCparticle::ParticleOriginMCTruth::kElse; + } + // if pdg code is wrong, particle is fake + } else { + particleOrigin = aod::femtodreamMCparticle::ParticleOriginMCTruth::kFake; + } + + outputPartsMc(particleOrigin, pdgCode, particleMc.pt(), particleMc.eta(), particleMc.phi()); + outputPartsMcLabels(outputPartsMc.lastIndex()); + if (isDebug) { + outputPartsExtMcLabels(outputPartsMc.lastIndex()); + outputDebugPartsMc(pdgCodeMother); + } + } else { + outputPartsMcLabels(-1); + if (isDebug) { + outputPartsExtMcLabels(-1); + } + } + } + + template + void fillMcCollision(CollisionType const& col) + { + if (col.has_mcCollision()) { + // auto genMCcol = col.template mcCollision_as(); + // outputMcCollision(genMCcol.multMCNParticlesEta08()); + outputCollsMcLabels(outputMcCollision.lastIndex()); + } else { + outputCollsMcLabels(-1); + } + } + + template + bool fillV0sForCharmHadron(CollisionType const& col, V0Type const& fullV0s, TrackType const&) + { + const bool isK0S = (selectionFlagV0 == V0Channel::K0S); + const bool isLam = (selectionFlagV0 == V0Channel::Lambda); + if (!isK0S && !isLam) { + LOG(fatal) << "Invalid V0 particle !! Please check the configuration"; + } + + std::vector childIDs = {0, 0}; + std::vector tmpIDtrack; + bool fIsV0Filled = false; + + auto bc = col.template bc_as(); + int64_t timeStamp = bc.timestamp(); + + for (const auto& v0 : fullV0s) { + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + + float massV0 = 0.f; + float antiMassV0 = 0.f; + + std::array cutContainerV0{}; + + if (isK0S) { + K0SCuts.fillLambdaQA(col, v0, postrack, negtrack); + if (!K0SCuts.isSelectedMinimal(col, v0, postrack, negtrack)) { + continue; + } + K0SCuts.fillQA(col, v0, postrack, negtrack); + + cutContainerV0 = K0SCuts.getCutContainer(col, v0, postrack, negtrack); + massV0 = v0.mK0Short(); + antiMassV0 = v0.mK0Short(); + + } else { // Lambda + LambdaCuts.fillLambdaQA(col, v0, postrack, negtrack); + if (!LambdaCuts.isSelectedMinimal(col, v0, postrack, negtrack)) { + continue; + } + LambdaCuts.fillQA(col, v0, postrack, negtrack); + cutContainerV0 = LambdaCuts.getCutContainer(col, v0, postrack, negtrack); + massV0 = v0.mLambda(); + antiMassV0 = v0.mAntiLambda(); + } + + int rowPos = getRowDaughters(v0.posTrackId(), tmpIDtrack); + int rowNeg = getRowDaughters(v0.negTrackId(), tmpIDtrack); + if (rowPos < 0 || rowNeg < 0) + continue; + + // --- pos child + childIDs[0] = rowPos; + childIDs[1] = 0; + + auto daughType = isK0S ? aod::femtodreamparticle::ParticleType::kV0K0ShortChild + : aod::femtodreamparticle::ParticleType::kV0Child; + + outputParts(outputCollision.lastIndex(), + v0.positivept(), v0.positiveeta(), v0.positivephi(), + daughType, + cutContainerV0.at(femto_dream_v0_selection::V0ContainerPosition::kPosCuts), + cutContainerV0.at(femto_dream_v0_selection::V0ContainerPosition::kPosPID), + postrack.dcaXY(), + childIDs, + 0, + 0); + const int rowOfPosTrack = outputParts.lastIndex(); + + // --- neg child + childIDs[0] = 0; + childIDs[1] = rowNeg; + + outputParts(outputCollision.lastIndex(), + v0.negativept(), v0.negativeeta(), v0.negativephi(), + daughType, + cutContainerV0.at(femto_dream_v0_selection::V0ContainerPosition::kNegCuts), + cutContainerV0.at(femto_dream_v0_selection::V0ContainerPosition::kNegPID), + negtrack.dcaXY(), + childIDs, + 0, + 0); + const int rowOfNegTrack = outputParts.lastIndex(); + + // --- mother + std::vector indexChildID = {rowOfPosTrack, rowOfNegTrack}; + auto motherType = isK0S ? aod::femtodreamparticle::ParticleType::kV0K0Short + : aod::femtodreamparticle::ParticleType::kV0; + + outputParts(outputCollision.lastIndex(), + v0.pt(), v0.eta(), v0.phi(), + motherType, + cutContainerV0.at(femto_dream_v0_selection::V0ContainerPosition::kV0), + 0, + v0.v0cosPA(), + indexChildID, + massV0, + antiMassV0); + + outputPartsTime(timeStamp); + outputPartsIndex(v0.globalIndex()); + + if (isDebug.value) { + fillDebugParticle(postrack); + fillDebugParticle(negtrack); + fillDebugParticle(v0); + } + + if constexpr (IsMc) { + fillMcParticle(col, v0, o2::aod::femtodreamparticle::ParticleType::kV0); + } + + fIsV0Filled = true; + } + return fIsV0Filled; + } + + template + void fillCharmHadronTable(CollisionType const& col, TrackType const& tracks, V0Type const& fullV0s, CandType const& candidates) + { + const auto vtxZ = col.posZ(); + const auto sizeCand = candidates.size(); + const auto spher = 2.; // dummy value for the moment + float mult = 0; + int multNtr = 0; + if (isRun3) { + if (useCent) { + mult = col.centFT0M(); + } else { + mult = 0; + } + multNtr = col.multNTracksPV(); + } else { + mult = 1; // multiplicity percentile is know in Run 2 + multNtr = col.multTracklets(); + } + + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(col, mult, ccdb, qaRegistry); + + qaRegistry.fill(HIST("hEventQA"), 1 + Event::All); + + /// monitor the satisfied event selections + hfEvSel.fillHistograms(col, rejectionMask, mult); + if (rejectionMask != 0) { + /// at least one event selection not satisfied --> reject the candidate + qaRegistry.fill(HIST("hEventQA"), 1 + Event::RejEveSel); + return; + } + + if (isNoSelectedV0s(col, tracks, fullV0s, K0SCuts) && sizeCand <= 0) { + qaRegistry.fill(HIST("hEventQA"), 1 + Event::RejNoV0sAndCharm); + return; + } + + outputCollision(vtxZ, mult, multNtr, spher, magField); + if constexpr (IsMc) { + fillMcCollision(col); + } + + // Filling candidate properties + if constexpr (Channel == DecayChannel::DplusToPiKPi || Channel == DecayChannel::LcToPKPi) { + rowCandCharm3Prong.reserve(sizeCand); + } else if constexpr (Channel == DecayChannel::D0ToPiK) { + rowCandCharm2Prong.reserve(sizeCand); + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { + rowCandCharmDstar.reserve(sizeCand); + } + bool isV0Filled = false; + bool isSelectedMlLcToPKPi = true; + bool isSelectedMlLcToPiKP = true; + bool isSelectedMlDplusToPiKPi = true; + bool isSelectedMlD0ToPiK = true; + bool isSelectedMlD0barToKPi = true; + bool isSelectedMlDstarToD0Pi = true; + + for (const auto& candidate : candidates) { + outputMlD0 = {-1.0f, -1.0f, -1.0f}; + outputMlD0bar = {-1.0f, -1.0f, -1.0f}; + outputMlDplus = {-1.0f, -1.0f, -1.0f}; + outputMlDstar = {-1.0f, -1.0f, -1.0f}; + outputMlPKPi = {-1.0f, -1.0f, -1.0f}; + outputMlPiKP = {-1.0f, -1.0f, -1.0f}; + auto trackPos1 = candidate.template prong0_as(); // positive daughter (negative for the antiparticles) + auto trackNeg = candidate.template prong1_as(); // negative daughter (positive for the antiparticles) + + auto bc = col.template bc_as(); + int64_t timeStamp = bc.timestamp(); + + auto fillTable = [&](int candFlag, + int functionSelection, + float bdtScoreBkg, + float bdtScorePrompt, + float bdtScoreFd) { + if (functionSelection >= 1) { + if constexpr (Channel == DecayChannel::DplusToPiKPi || Channel == DecayChannel::LcToPKPi) { + auto trackPos2 = candidate.template prong2_as(); + rowCandCharm3Prong( + outputCollision.lastIndex(), + timeStamp, + trackPos1.sign() + trackNeg.sign() + trackPos2.sign(), + trackPos1.globalIndex(), + trackNeg.globalIndex(), + trackPos2.globalIndex(), + trackPos1.pt(), + trackNeg.pt(), + trackPos2.pt(), + trackPos1.eta(), + trackNeg.eta(), + trackPos2.eta(), + trackPos1.phi(), + trackNeg.phi(), + trackPos2.phi(), + 1 << candFlag, + bdtScoreBkg, + bdtScorePrompt, + bdtScoreFd); + + } else if constexpr (Channel == DecayChannel::D0ToPiK) { + rowCandCharm2Prong( + outputCollision.lastIndex(), + timeStamp, + trackPos1.sign() + trackNeg.sign(), + trackPos1.globalIndex(), + trackNeg.globalIndex(), + trackPos1.pt(), + trackNeg.pt(), + trackPos1.eta(), + trackNeg.eta(), + trackPos1.phi(), + trackNeg.phi(), + 1 << candFlag, + bdtScoreBkg, + bdtScorePrompt, + bdtScoreFd); + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { + auto trackPos2 = candidate.template prongPi_as(); + rowCandCharmDstar( + outputCollision.lastIndex(), + timeStamp, + candidate.template prongPi_as().sign(), + trackPos1.globalIndex(), + trackNeg.globalIndex(), + trackPos2.globalIndex(), + trackPos1.pt(), + trackNeg.pt(), + trackPos2.pt(), + trackPos1.eta(), + trackNeg.eta(), + trackPos2.eta(), + trackPos1.phi(), + trackNeg.phi(), + trackPos2.phi(), + 1 << candFlag, + bdtScoreBkg, + bdtScorePrompt, + bdtScoreFd); + } + + if constexpr (IsMc) { + rowCandMcCharmHad( + candidate.flagMcMatchRec(), + candidate.originMcRec()); + } + } + }; + + if constexpr (Channel == DecayChannel::DplusToPiKPi) { + if constexpr (UseCharmMl) { + /// fill with ML information + /// BDT index 0: bkg score; BDT index 1: prompt score; BDT index 2: non-prompt score + if (applyMlMode == FillMlFromSelector) { + if (candidate.mlProbDplusToPiKPi().size() > 0) { + outputMlDplus.at(0) = candidate.mlProbDplusToPiKPi()[0]; /// bkg score + outputMlDplus.at(1) = candidate.mlProbDplusToPiKPi()[1]; /// prompt score + outputMlDplus.at(2) = candidate.mlProbDplusToPiKPi()[2]; /// non-prompt score + } + } else if (applyMlMode == FillMlFromNewBDT) { + isSelectedMlDplusToPiKPi = false; + if (candidate.mlProbDplusToPiKPi().size() > 0) { + std::vector inputFeaturesDplusToPiKPi = hfMlResponseDplus.getInputFeatures(candidate); + isSelectedMlDplusToPiKPi = hfMlResponseDplus.isSelectedMl(inputFeaturesDplusToPiKPi, candidate.pt(), outputMlDplus); + } + if (!isSelectedMlDplusToPiKPi) { + continue; + } + } else { + LOGF(fatal, "Please check your Ml configuration!!"); + } + } + fillTable(2, candidate.isSelDplusToPiKPi(), outputMlDplus.at(0), outputMlDplus.at(1), outputMlDplus.at(2)); + + } else if constexpr (Channel == DecayChannel::LcToPKPi) { + if constexpr (UseCharmMl) { + /// fill with ML information + /// BDT index 0: bkg score; BDT index 1: prompt score; BDT index 2: non-prompt score + if (applyMlMode == FillMlFromSelector) { + if (candidate.mlProbLcToPKPi().size() > 0) { + outputMlPKPi.at(0) = candidate.mlProbLcToPKPi()[0]; /// bkg score + outputMlPKPi.at(1) = candidate.mlProbLcToPKPi()[1]; /// prompt score + outputMlPKPi.at(2) = candidate.mlProbLcToPKPi()[2]; /// non-prompt score + } + if (candidate.mlProbLcToPiKP().size() > 0) { + outputMlPiKP.at(0) = candidate.mlProbLcToPiKP()[0]; /// bkg score + outputMlPiKP.at(1) = candidate.mlProbLcToPiKP()[1]; /// prompt score + outputMlPiKP.at(2) = candidate.mlProbLcToPiKP()[2]; /// non-prompt score + } + } else if (applyMlMode == FillMlFromNewBDT) { + isSelectedMlLcToPKPi = false; + isSelectedMlLcToPiKP = false; + if (candidate.mlProbLcToPKPi().size() > 0) { + std::vector inputFeaturesLcToPKPi = hfMlResponseLc.getInputFeatures(candidate, true); + isSelectedMlLcToPKPi = hfMlResponseLc.isSelectedMl(inputFeaturesLcToPKPi, candidate.pt(), outputMlPKPi); + } + if (candidate.mlProbLcToPiKP().size() > 0) { + std::vector inputFeaturesLcToPiKP = hfMlResponseLc.getInputFeatures(candidate, false); + isSelectedMlLcToPiKP = hfMlResponseLc.isSelectedMl(inputFeaturesLcToPiKP, candidate.pt(), outputMlPKPi); + } + if (!isSelectedMlLcToPKPi && !isSelectedMlLcToPiKP) { + continue; + } + } else { + LOGF(fatal, "Please check your Ml configuration!!"); + } + } + fillTable(0, candidate.isSelLcToPKPi(), outputMlPKPi.at(0), outputMlPKPi.at(1), outputMlPKPi.at(2)); + fillTable(1, candidate.isSelLcToPiKP(), outputMlPiKP.at(0), outputMlPiKP.at(1), outputMlPiKP.at(2)); + } else if constexpr (Channel == DecayChannel::D0ToPiK) { + if constexpr (UseCharmMl) { + + /// fill with ML information + /// BDT index 0: bkg score; BDT index 1: prompt score; BDT index 2: non-prompt score + if (applyMlMode == FillMlFromSelector) { + if (candidate.mlProbD0().size() > 0) { + outputMlD0.at(0) = candidate.mlProbD0()[0]; /// bkg score + outputMlD0.at(1) = candidate.mlProbD0()[1]; /// prompt score + outputMlD0.at(2) = candidate.mlProbD0()[2]; /// non-prompt score + } + if (candidate.mlProbD0bar().size() > 0) { + outputMlD0bar.at(0) = candidate.mlProbD0bar()[0]; /// bkg score + outputMlD0bar.at(1) = candidate.mlProbD0bar()[1]; /// prompt score + outputMlD0bar.at(2) = candidate.mlProbD0bar()[2]; /// non-prompt score + } + + } else if (applyMlMode == FillMlFromNewBDT) { + isSelectedMlD0ToPiK = false; + isSelectedMlD0barToKPi = false; + + if (candidate.mlProbD0().size() > 0) { + std::vector inputFeaturesD0ToPiK = hfMlResponseD0.getInputFeatures(candidate, o2::constants::physics::kD0); + isSelectedMlD0ToPiK = hfMlResponseD0.isSelectedMl(inputFeaturesD0ToPiK, candidate.pt(), outputMlD0); + } + if (candidate.mlProbD0bar().size() > 0) { + std::vector inputFeaturesD0barToKPi = hfMlResponseD0.getInputFeatures(candidate, o2::constants::physics::kD0); + isSelectedMlD0barToKPi = hfMlResponseD0.isSelectedMl(inputFeaturesD0barToKPi, candidate.pt(), outputMlD0bar); + } + if (!isSelectedMlD0ToPiK && !isSelectedMlD0barToKPi) { + continue; + } + } else { + LOGF(fatal, "Please check your Ml configuration!!"); + } + } + fillTable(0, candidate.isSelD0(), outputMlD0.at(0), outputMlD0.at(1), outputMlD0.at(2)); + fillTable(1, candidate.isSelD0bar(), outputMlD0bar.at(0), outputMlD0bar.at(1), outputMlD0bar.at(2)); + + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { + if constexpr (UseCharmMl) { + /// fill with ML information + /// BDT index 0: bkg score; BDT index 1: prompt score; BDT index 2: non-prompt score + if (applyMlMode == FillMlFromSelector) { + if (candidate.mlProbDstarToD0Pi().size() > 0) { + outputMlDstar.at(0) = candidate.mlProbDstarToD0Pi()[0]; /// bkg score + outputMlDstar.at(1) = candidate.mlProbDstarToD0Pi()[1]; /// prompt score + outputMlDstar.at(2) = candidate.mlProbDstarToD0Pi()[2]; /// non-prompt score + } + } else if (applyMlMode == FillMlFromNewBDT) { + isSelectedMlDstarToD0Pi = false; + if (candidate.mlProbDstarToD0Pi().size() > 0) { + std::vector inputFeaturesDstarToD0Pi = hfMlResponseDstar.getInputFeatures(candidate, false); + isSelectedMlDstarToD0Pi = hfMlResponseDstar.isSelectedMl(inputFeaturesDstarToD0Pi, candidate.pt(), outputMlDstar); + } + if (!isSelectedMlDstarToD0Pi) { + continue; + } + } else { + LOGF(fatal, "Please check your Ml configuration!!"); + } + } + fillTable(2, candidate.isSelDstarToD0Pi(), outputMlDstar.at(0), outputMlDstar.at(1), outputMlDstar.at(2)); + } + } + isV0Filled = fillV0sForCharmHadron(col, fullV0s, tracks); + + aod::femtodreamcollision::BitMaskType bitV0 = 0; + if (isV0Filled) { + bitV0 |= 1 << 0; + qaRegistry.fill(HIST("hEventQA"), 1 + Event::V0Selected); + } + + aod::femtodreamcollision::BitMaskType bitCand = 0; + if (sizeCand > 0) { + bitCand |= 1 << 0; + qaRegistry.fill(HIST("hEventQA"), 1 + Event::CharmSelected); + } + + if (isV0Filled && (sizeCand > 0)) { + qaRegistry.fill(HIST("hEventQA"), 1 + Event::PairSelected); + } + + rowMasks(bitV0, + bitCand, + 0); + } + + // check if there is no selected v0 + /// \param C type of the collision + /// \param T type of the V0s + /// \param TC type of the femto V0s cuts + /// \return whether or not the V0s fulfills the all selections + template + bool isNoSelectedV0s(CollType const& col, TrackType const& /*col*/, TV0 const& V0s, CutType& v0Cuts) + { + for (auto const& v0 : V0s) { + auto postrack = v0.template posTrack_as(); + auto negtrack = v0.template negTrack_as(); + if (v0Cuts.isSelectedMinimal(col, v0, postrack, negtrack)) { + return false; + } + } + return true; + } + + template + int getRowDaughters(int daughID, T const& vecID) + { + int rowInPrimaryTrackTableDaugh = -1; + for (size_t i = 0; i < vecID.size(); i++) { + if (vecID.at(i) == daughID) { + rowInPrimaryTrackTableDaugh = i; + break; + } + } + return rowInPrimaryTrackTableDaugh; + } + + template + void fillCharmHadMcGen(ParticleType particles) + { + // Filling particle properties + rowCandCharmHadGen.reserve(particles.size()); + if constexpr (Channel == DecayChannel::DplusToPiKPi) { + for (const auto& particle : particles) { + if (std::abs(particle.flagMcMatchGen()) == hf_decay::hf_cand_3prong::DecayChannelMain::DplusToPiKPi) { + rowCandCharmHadGen( + particle.mcCollisionId(), + particle.flagMcMatchGen(), + particle.originMcGen()); + } + } + } else if constexpr (Channel == DecayChannel::LcToPKPi) { + for (const auto& particle : particles) { + if (std::abs(particle.flagMcMatchGen()) == hf_decay::hf_cand_3prong::DecayChannelMain::LcToPKPi) { + rowCandCharmHadGen( + particle.mcCollisionId(), + particle.flagMcMatchGen(), + particle.originMcGen()); + } + } + } else if constexpr (Channel == DecayChannel::D0ToPiK) { + for (const auto& particle : particles) { + if (std::abs(particle.flagMcMatchGen()) == hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK) { + rowCandCharmHadGen( + particle.mcCollisionId(), + particle.flagMcMatchGen(), + particle.originMcGen()); + } + } + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { + for (const auto& particle : particles) { + if (std::abs(particle.flagMcMatchGen()) == hf_decay::hf_cand_dstar::DecayChannelMain::DstarToPiKPi) { + rowCandCharmHadGen( + particle.mcCollisionId(), + particle.flagMcMatchGen(), + particle.originMcGen()); + } + } + } + } + + /// D0ToPiK + void processDataD0ToPiK(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataD0ToPiK, "Provide experimental data for D0ToPiK femto", false); + + void processDataD0ToPiKWithML(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered> const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataD0ToPiKWithML, "Provide experimental data for D0ToPiK with ml", false); + + void processMcD0ToPiK(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + CandidateD0Mc const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcD0ToPiK, "Provide Mc for D0ToPiK", false); + + void processMcD0ToPiKWithML(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + soa::Join const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcD0ToPiKWithML, "Provide Mc for D0ToPiK with ml", false); + + void processMcD0ToPiKGen(GeneratedMc const& particles) + { + fillCharmHadMcGen(particles); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcD0ToPiKGen, "Provide Mc Generated D0ToPiK", false); + + /// DstarToD0Pi + void processDataDstarToD0Pi(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataDstarToD0Pi, "Provide experimental data for DstarToD0Pi femto", false); + + void processDataDstarToD0PiWithML(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered> const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataDstarToD0PiWithML, "Provide experimental data for DstarToD0Pi with ml", false); + + void processMcDstarToD0Pi(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + CandidateDstarMc const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcDstarToD0Pi, "Provide Mc for DstarToD0Pi", false); + + void processMcDstarToD0PiWithML(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + soa::Join const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcDstarToD0PiWithML, "Provide Mc for DstarToD0Pi with ml", false); + + void processMcDstarToD0PiGen(GeneratedMc const& particles) + { + + fillCharmHadMcGen(particles); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcDstarToD0PiGen, "Provide Mc Generated DstarToD0Pi", false); + + /// DplusToPiKPi + void processDataDplusToPiKPi(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataDplusToPiKPi, "Provide experimental data for DplusToPiKPi femto", false); + + void processDataDplusToPiKPiWithML(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered> const& candidates) + { + + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataDplusToPiKPiWithML, "Provide experimental data for DplusToPiKPi with ml", false); + + void processMcDplusToPiKPi(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + CandidateDplusMc const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcDplusToPiKPi, "Provide Mc for DplusToPiKPi", false); + + void processMcDplusToPiKPiWithML(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + soa::Join const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcDplusToPiKPiWithML, "Provide Mc for DplusToPiKPi with ml", false); + + void processMcDplusToPiKPiGen(GeneratedMc const& particles) + { + + fillCharmHadMcGen(particles); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcDplusToPiKPiGen, "Provide Mc Generated DplusToPiKPi", false); + + /// LcToPKPi + void processDataLcToPKPi(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataLcToPKPi, "Provide experimental data for Lc(PKPi)-proton femto", false); + + void processDataLcToPKPiWithML(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + soa::Filtered> const& candidates) + { + + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processDataLcToPKPiWithML, "Provide experimental data for Lc(PKPi)-proton femto with ml", false); + + void processMcLcToPKPi(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + CandidateLcMc const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcLcToPKPi, "Provide Mc for lctopkpi", false); + + void processMcLcToPKPiWithML(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + soa::Join const& fullV0s, + aod::McParticles const&, + soa::Join const& candidates) + { + // get magnetic field for run + getMagneticFieldTesla(col.bc_as()); + fillCharmHadronTable(col, tracks, fullV0s, candidates); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcLcToPKPiWithML, "Provide Mc for lctopkpi with ml", false); + + void processMcLcToPKPiGen(GeneratedMc const& particles) + { + + fillCharmHadMcGen(particles); + } + PROCESS_SWITCH(HfProducerCharmHadronsV0FemtoDream, processMcLcToPKPiGen, "Provide Mc Generated lctopkpi", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/HFC/Tasks/CMakeLists.txt b/PWGHF/HFC/Tasks/CMakeLists.txt index 9db29e2030f..5132b7174b7 100644 --- a/PWGHF/HFC/Tasks/CMakeLists.txt +++ b/PWGHF/HFC/Tasks/CMakeLists.txt @@ -14,6 +14,11 @@ o2physics_add_dpl_workflow(task-charm-hadrons-track-femto-dream PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(task-charm-hadrons-v0-femto-dream + SOURCES taskCharmHadronsV0FemtoDream.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(task-correlation-d0-hadrons SOURCES taskCorrelationD0Hadrons.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGHF/HFC/Tasks/taskCharmHadronsV0FemtoDream.cxx b/PWGHF/HFC/Tasks/taskCharmHadronsV0FemtoDream.cxx new file mode 100644 index 00000000000..25b9941e274 --- /dev/null +++ b/PWGHF/HFC/Tasks/taskCharmHadronsV0FemtoDream.cxx @@ -0,0 +1,951 @@ +// Copyright 2019-2025 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 taskCharmHadronsTrackFemtoDream.cxx +/// \brief Tasks that reads the V0 and CharmHadrons tables used for the pairing and builds pairs +/// \author Biao Zhang, Heidelberg University, biao.zhang@cern.ch + +#include "PWGCF/DataModel/FemtoDerived.h" +#include "PWGCF/FemtoDream/Core/femtoDreamContainer.h" +#include "PWGCF/FemtoDream/Core/femtoDreamDetaDphiStar.h" +#include "PWGCF/FemtoDream/Core/femtoDreamEventHisto.h" +#include "PWGCF/FemtoDream/Core/femtoDreamMath.h" +#include "PWGCF/FemtoDream/Core/femtoDreamPairCleaner.h" +#include "PWGCF/FemtoDream/Core/femtoDreamParticleHisto.h" +#include "PWGCF/FemtoDream/Core/femtoDreamUtils.h" +#include "PWGHF/Core/DecayChannels.h" + +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::soa; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::analysis::femtoDream; +using namespace o2::hf_decay; +using namespace o2::hf_decay::hf_cand_3prong; +using namespace o2::constants::physics; + +struct HfTaskCharmHadronsV0FemtoDream { + + enum TrackCharge { + PositiveCharge = 1, + NegativeCharge = -1 + }; + + enum PairSign { + PairNotDefined = 0, + LikeSignPair = 1, + UnLikeSignPair = 2 + }; + // decay channels + enum DecayChannel { DplusToPiKPi = 0, + LcToPKPi, + D0ToPiK, + DstarToD0Pi + }; + + constexpr static int OriginRecPrompt = 1; + constexpr static int OriginRecFD = 2; + constexpr static int CutBitChargePositive = 2; + + Produces rowFemtoResultPairs; + Produces rowFemtoResultCharm3Prong; + Produces rowFemtoResultCharm2Prong; + Produces rowFemtoResultCharmDstar; + Produces rowFemtoResultV0; + Produces rowFemtoResultColl; + + Configurable confTempFitVarMomentum{"confTempFitVarMomentum", 0, "Momentum used for binning: 0 -> pt; 1 -> preco; 2 -> ptpc"}; + + /// Particle 2 (Charm Hadrons) + struct : ConfigurableGroup { + Configurable charmHadBkgBDTmax{"charmHadBkgBDTmax", 1., "Maximum background bdt score for Charm Hadron (particle 2)"}; + Configurable charmHadCandSel{"charmHadCandSel", 1, "candidate selection for charm hadron"}; + Configurable charmHadMcSel{"charmHadMcSel", DecayChannelMain::LcToPKPi, "charm hadron selection for mc, DplusToPiKPi = 1, LcToPKPi = 17"}; + Configurable charmHadFdBDTmin{"charmHadFdBDTmin", 0., "Minimum feed-down bdt score Charm Hadron (particle 2)"}; + Configurable charmHadFdBDTmax{"charmHadFdBDTmax", 1., "Maximum feed-down bdt score Charm Hadron (particle 2)"}; + Configurable charmHadMaxInvMass{"charmHadMaxInvMass", 2.45, "Maximum invariant mass of Charm Hadron (particle 2)"}; + Configurable charmHadMinInvMass{"charmHadMinInvMass", 2.15, "Minimum invariant mass of Charm Hadron (particle 2)"}; + Configurable charmHadMinPt{"charmHadMinPt", 0., "Minimum pT of Charm Hadron (particle 2)"}; + Configurable charmHadMaxPt{"charmHadMaxPt", 999., "Maximum pT of Charm Hadron (particle 2)"}; + Configurable charmHadPDGCode{"charmHadPDGCode", 4122, "PDG code of particle 2 Charm Hadron"}; + Configurable charmHadPromptBDTmin{"charmHadPromptBDTmin", 0., "Minimum prompt bdt score Charm Hadron (particle 2)"}; + Configurable charmHadPromptBDTmax{"charmHadPromptBDTmax", 1., "Maximum prompt bdt score Charm Hadron (particle 2)"}; + } charmSel; + /// General options + Configurable cprDeltaEtaMax{"cprDeltaEtaMax", 0.01, "Max. Delta Eta for Close Pair Rejection"}; + Configurable cprDeltaPhiMax{"cprDeltaPhiMax", 0.01, "Max. Delta Phi for Close Pair Rejection"}; + Configurable cprPlotPerRadii{"cprPlotPerRadii", false, "Plot CPR per radii"}; + Configurable extendedPlots{"extendedPlots", false, "Enable additional three dimensional histogramms. High memory consumption. Use for debugging"}; + Configurable highkstarCut{"highkstarCut", 100000., "Set a cut for high k*, above which the pairs are rejected"}; + Configurable isMc{"isMc", false, "Set true in the case of a MonteCarlo Run"}; + Configurable smearingByOrigin{"smearingByOrigin", false, "Obtain the smearing matrix differential in the MC origin of particle 1 and particle 2. High memory consumption. Use with care!"}; + Configurable use4D{"use4D", false, "Enable four dimensional histogramms (to be used only for analysis with high statistics): k* vs multiplicity vs multiplicity percentil vs mT"}; + Configurable useCPR{"useCPR", false, "Close Pair Rejection"}; + Configurable fillTableWithCharm{"fillTableWithCharm", true, "Write charm/tracks/collision table only if >=1 charm hadron in this collision"}; + + // Mixing configurables + struct : ConfigurableGroup { + Configurable doMixEvent{"doMixEvent", false, "choose do mix-event online"}; + Configurable mixingBinPolicy{"mixingBinPolicy", 0, "Binning policy for mixing - 0: multiplicity, 1: multipliciy percentile, 2: both"}; + Configurable mixingDepth{"mixingDepth", 5, "Number of events for mixing"}; + } mixSetting; + /// Event selection + struct : ConfigurableGroup { + std::string prefix = "eventSel"; + Configurable multMin{"multMin", 0, "Minimum Multiplicity (MultNtr)"}; + Configurable multMax{"multMax", 99999, "Maximum Multiplicity (MultNtr)"}; + Configurable multPercentileMin{"multPercentileMin", 0, "Maximum Multiplicity Percentile"}; + Configurable multPercentileMax{"multPercentileMax", 100, "Minimum Multiplicity Percentile"}; + } eventSel; + + /// particle 1 (V0), Λ or Ks0 + struct : ConfigurableGroup { + std::string prefix = std::string("v0Sel"); + Configurable pdgCodeV0{"pdgCodeV0", 310, "PDG code of V0 (310: K0S, 3122: Lambda)"}; + Configurable cutBit{"cutBit", 7518, "Selection bit for particle 1 (V0)"}; + Configurable childPosCutBit{"childPosCutBit", 210, "Selection bit for positive child of V01"}; + Configurable childPosTPCBit{"childPosTPCBit", 64, "PID TPC bit for positive child of V01"}; + Configurable childNegCutBit{"childNegCutBit", 209, "Selection bit for negative child of V01"}; + Configurable childNegTPCBit{"childNegTPCBit", 256, "PID TPC bit for negative child of V01"}; + Configurable invMassV0Min{"invMassV0Min", 0.45, "Minimum invariant mass of Partricle 1 (particle) (V0)"}; + Configurable invMassV0Max{"invMassV0Max", 0.55, "Maximum invariant mass of Partricle 1 (particle) (V0)"}; + Configurable invMassAntiV0Min{"invMassAntiV0Min", 0.45, "Minimum invariant mass of Partricle 1 (antiparticle) (V0)"}; + Configurable invMassAntiV0Max{"invMassAntiV0Max", 0.55, "Maximum invariant mass of Partricle 1 (antiparticle) (V0)"}; + Configurable ptV0Min{"ptV0Min", 0., "Minimum pT of Partricle 1 (V0)"}; + Configurable ptV0Max{"ptV0Max", 999., "Maximum pT of Partricle 1 (V0)"}; + Configurable etaV0Min{"etaV0Min", -10., "Minimum eta of Partricle 1 (V0)"}; + Configurable etaV0Max{"etaV0Max", 10., "Maximum eta of Partricle 1 (V0)"}; + } v0Sel; + + SliceCache cache; + + using FilteredCharmCand3Prongs = soa::Filtered; + using FilteredCharmCand3Prong = FilteredCharmCand3Prongs::iterator; + + using FilteredCharmCand2Prongs = soa::Filtered; + using FilteredCharmCand2Prong = FilteredCharmCand2Prongs::iterator; + + using FilteredCharmCandDstars = soa::Filtered; + using FilteredCharmCandDstar = FilteredCharmCandDstars::iterator; + + using FilteredCharmMcCand3Prongs = soa::Filtered>; + using FilteredCharmMcCand3Prong = FilteredCharmMcCand3Prongs::iterator; + + using FilteredCharmMcCand2Prongs = soa::Filtered>; + using FilteredCharmMcCand2Prong = FilteredCharmMcCand2Prongs::iterator; + + using FilteredCharmMcCandDstars = soa::Filtered>; + using FilteredCharmMcCandDstar = FilteredCharmMcCandDstars::iterator; + + using FilteredCollisions = soa::Filtered>; + using FilteredCollision = FilteredCollisions::iterator; + + using FilteredMcColisions = soa::Filtered>; + using FilteredMcColision = FilteredMcColisions::iterator; + + using FilteredFDMcParts = soa::Filtered>; + using FilteredFDMcPart = FilteredFDMcParts::iterator; + + using FilteredFDParticles = soa::Filtered>; + using FilteredFDParticle = FilteredFDParticles::iterator; + + Filter eventMultiplicity = aod::femtodreamcollision::multNtr >= eventSel.multMin && aod::femtodreamcollision::multNtr <= eventSel.multMax; + Filter eventMultiplicityPercentile = aod::femtodreamcollision::multV0M >= eventSel.multPercentileMin && aod::femtodreamcollision::multV0M <= eventSel.multPercentileMax; + Filter hfCandSelFilter = aod::fdhf::candidateSelFlag >= charmSel.charmHadCandSel; + Filter hfMcSelFilter = (nabs(aod::fdhf::flagMc) == charmSel.charmHadMcSel); + + Preslice perCol = aod::femtodreamparticle::fdCollisionId; + Preslice perHf3ProngByCol = aod::femtodreamparticle::fdCollisionId; + Preslice perHf2ProngByCol = aod::femtodreamparticle::fdCollisionId; + Preslice perHfDstarByCol = aod::femtodreamparticle::fdCollisionId; + + /// Partition for particle Lambda + Partition partitionLambda = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kV0)) && + ((aod::femtodreamparticle::cut & v0Sel.cutBit) == v0Sel.cutBit) && + (aod::femtodreamparticle::pt > v0Sel.ptV0Min) && + (aod::femtodreamparticle::pt < v0Sel.ptV0Max) && + (aod::femtodreamparticle::eta > v0Sel.etaV0Min) && + (aod::femtodreamparticle::eta < v0Sel.etaV0Max) && + (aod::femtodreamparticle::mLambda > v0Sel.invMassV0Min) && + (aod::femtodreamparticle::mLambda < v0Sel.invMassV0Max) && + (aod::femtodreamparticle::mAntiLambda > v0Sel.invMassAntiV0Min) && + (aod::femtodreamparticle::mAntiLambda < v0Sel.invMassAntiV0Max); + + /// Partition for particle K0Short + Partition partitionK0Short = (aod::femtodreamparticle::partType == uint8_t(aod::femtodreamparticle::ParticleType::kV0K0Short)) && + ((aod::femtodreamparticle::cut & v0Sel.cutBit) == v0Sel.cutBit) && + (aod::femtodreamparticle::pt > v0Sel.ptV0Min) && + (aod::femtodreamparticle::pt < v0Sel.ptV0Max) && + (aod::femtodreamparticle::eta > v0Sel.etaV0Min) && + (aod::femtodreamparticle::eta < v0Sel.etaV0Max) && + (aod::femtodreamparticle::mLambda > v0Sel.invMassV0Min) && + (aod::femtodreamparticle::mLambda < v0Sel.invMassV0Max) && + (aod::femtodreamparticle::mAntiLambda > v0Sel.invMassAntiV0Min) && + (aod::femtodreamparticle::mAntiLambda < v0Sel.invMassAntiV0Max); + + /// Partition for particle 2 + Partition partitionCharmHadron3Prong = aod::fdhf::bdtBkg < charmSel.charmHadBkgBDTmax && aod::fdhf::bdtFD < charmSel.charmHadFdBDTmax && aod::fdhf::bdtFD > charmSel.charmHadFdBDTmin&& aod::fdhf::bdtPrompt charmSel.charmHadPromptBDTmin; + Partition partitionCharmHadron2Prong = aod::fdhf::bdtBkg < charmSel.charmHadBkgBDTmax && aod::fdhf::bdtFD < charmSel.charmHadFdBDTmax && aod::fdhf::bdtFD > charmSel.charmHadFdBDTmin&& aod::fdhf::bdtPrompt charmSel.charmHadPromptBDTmin; + Partition partitionCharmHadronDstar = aod::fdhf::bdtBkg < charmSel.charmHadBkgBDTmax && aod::fdhf::bdtFD < charmSel.charmHadFdBDTmax && aod::fdhf::bdtFD > charmSel.charmHadFdBDTmin&& aod::fdhf::bdtPrompt charmSel.charmHadPromptBDTmin; + + Partition partitionMcCharmHadron3Prong = aod::fdhf::originMcRec == OriginRecPrompt || aod::fdhf::originMcRec == OriginRecFD; + Partition partitionMcCharmHadron2Prong = aod::fdhf::originMcRec == OriginRecPrompt || aod::fdhf::originMcRec == OriginRecFD; + Partition partitionMcCharmHadronDstar = aod::fdhf::originMcRec == OriginRecPrompt || aod::fdhf::originMcRec == OriginRecFD; + struct : ConfigurableGroup { + /// Axis configurables + ConfigurableAxis dummy{"dummy", {1, 0, 1}, "dummy axis"}; + /// Binning configurables + ConfigurableAxis bin4Dkstar{"bin4Dkstar", {1500, 0., 6.}, "binning kstar for the 4Dimensional plot: k* vs multiplicity vs multiplicity percentile vs mT (set <> to true in order to use)"}; + ConfigurableAxis bin4DMult{"bin4DMult", {VARIABLE_WIDTH, 0.0f, 4.0f, 8.0f, 12.0f, 16.0f, 20.0f, 24.0f, 28.0f, 32.0f, 36.0f, 40.0f, 44.0f, 48.0f, 52.0f, 56.0f, 60.0f, 64.0f, 68.0f, 72.0f, 76.0f, 80.0f, 84.0f, 88.0f, 92.0f, 96.0f, 100.0f, 200.0f}, "multiplicity Binning for the 4Dimensional plot: k* vs multiplicity vs multiplicity percentile vs mT (set <> to true in order to use)"}; + ConfigurableAxis bin4DmT{"bin4DmT", {VARIABLE_WIDTH, 1.02f, 1.14f, 1.20f, 1.26f, 1.38f, 1.56f, 1.86f, 4.50f}, "mT Binning for the 4Dimensional plot: k* vs multiplicity vs multiplicity percentile vs mT (set <> to true in order to use)"}; + ConfigurableAxis bin4DmultPercentile{"bin4DmultPercentile", {10, 0.0f, 100.0f}, "multiplicity percentile Binning for the 4Dimensional plot: k* vs multiplicity vs multiplicity percentile vs mT (set <> to true in order to use)"}; + ConfigurableAxis binInvMass{"binInvMass", {400, 2.10, 2.50}, "InvMass binning"}; + ConfigurableAxis binpTCharm{"binpTCharm", {360, 0, 36}, "pT binning of charm hadron"}; + ConfigurableAxis binTempFitVarV0{"binTempFitVarV0", {300, 0.9, 1}, "binning of the TempFitVar in the pT vs. TempFitVar plot (charm))"}; + ConfigurableAxis binTempFitVarV0Child{"binTempFitVarV0Child", {300, -0.15, 0.15}, "binning of the TempFitVar in the pT vs. TempFitVar plot (V0 child)"}; + ConfigurableAxis binmT{"binmT", {225, 0., 7.5}, "binning mT"}; + ConfigurableAxis binmultTempFit{"binmultTempFit", {1, 0, 1}, "multiplicity Binning for the TempFitVar plot"}; + ConfigurableAxis binMulPercentile{"binMulPercentile", {10, 0.0f, 100.0f}, "multiplicity percentile Binning"}; + ConfigurableAxis binpTV0{"binpTV0", {50, 0.5, 10.05}, "pT binning of the pT vs. TempFitVar plot (Track)"}; + ConfigurableAxis binpTV0Child{"binpTV0Child", {20, 0.5, 4.05}, "pT binning of the pT vs. TempFitVar plot (V0 Child)"}; + ConfigurableAxis binEta{"binEta", {{200, -1.5, 1.5}}, "eta binning"}; + ConfigurableAxis binPhi{"binPhi", {{200, 0, o2::constants::math::TwoPI}}, "phi binning"}; + ConfigurableAxis binkT{"binkT", {150, 0., 9.}, "binning kT"}; + ConfigurableAxis binkstar{"binkstar", {1500, 0., 6.}, "binning kstar"}; + ConfigurableAxis binNSigmaTPC{"binNSigmaTPC", {1600, -8, 8}, "Binning of Nsigma TPC plot"}; + ConfigurableAxis binNSigmaTOF{"binNSigmaTOF", {3000, -15, 15}, "Binning of the Nsigma TOF plot"}; + ConfigurableAxis binNSigmaTPCTOF{"binNSigmaTPCTOF", {3000, -15, 15}, "Binning of the Nsigma TPC+TOF plot"}; + ConfigurableAxis binTPCClusters{"binTPCClusters", {163, -0.5, 162.5}, "Binning of TPC found clusters plot"}; + ConfigurableAxis mixingBinMult{"mixingBinMult", {VARIABLE_WIDTH, 0.0f, 20.0f, 60.0f, 200.0f}, "Mixing bins - multiplicity"}; + ConfigurableAxis mixingBinMultPercentile{"mixingBinMultPercentile", {VARIABLE_WIDTH, 0.0f, 100.f}, "Mixing bins - multiplicity percentile"}; + ConfigurableAxis mixingBinVztx{"mixingBinVztx", {VARIABLE_WIDTH, -10.0f, -4.f, 0.f, 4.f, 10.f}, "Mixing bins - z-vertex"}; + } AxisBinning; + + ColumnBinningPolicy colBinningMult{{AxisBinning.mixingBinVztx, AxisBinning.mixingBinMult}, true}; + ColumnBinningPolicy colBinningMultPercentile{{AxisBinning.mixingBinVztx, AxisBinning.mixingBinMultPercentile}, true}; + ColumnBinningPolicy colBinningMultMultPercentile{{AxisBinning.mixingBinVztx, AxisBinning.mixingBinMult, AxisBinning.mixingBinMultPercentile}, true}; + + FemtoDreamContainer sameEventCont; + FemtoDreamContainer mixedEventCont; + FemtoDreamPairCleaner pairCleaner3Prong; + FemtoDreamDetaDphiStar pairCloseRejectionSE3Prong; + FemtoDreamDetaDphiStar pairCloseRejectionME3Prong; + + FemtoDreamPairCleaner pairCleaner2Prong; + FemtoDreamDetaDphiStar pairCloseRejectionSE2Prong; + FemtoDreamDetaDphiStar pairCloseRejectionME2Prong; + + FemtoDreamPairCleaner pairCleanerDstar; + FemtoDreamDetaDphiStar pairCloseRejectionSEDstar; + FemtoDreamDetaDphiStar pairCloseRejectionMEDstar; + + femtodreamcollision::BitMaskType bitMask = 1 << 0; + + /// Histogramming for particle 1 + FemtoDreamParticleHisto v0HistoPartOne; + FemtoDreamParticleHisto v0HistoPartOneSelected; + FemtoDreamParticleHisto posChildHistos; + FemtoDreamParticleHisto negChildHistos; + /// Histogramming for Event + FemtoDreamEventHisto eventHisto; + /// Histogram output + HistogramRegistry registry{"CorrelationsAndQA", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry registryMixQa{"registryMixQa"}; + HistogramRegistry registryCharmHadronQa{"registryCharmHadronQa"}; + + float massOne = o2::analysis::femtoDream::getMass(v0Sel.pdgCodeV0); + float massTwo = o2::analysis::femtoDream::getMass(charmSel.charmHadPDGCode); + int8_t partSign = 0; + int64_t processType = 0; + + void init(InitContext& /*context*/) + { + std::array processes = {doprocessDataLcK0s, doprocessDataDplusK0s, doprocessDataD0K0s, doprocessDataDstarK0s}; + if (std::accumulate(processes.begin(), processes.end(), 0) != 1) { + LOGP(fatal, "One and only one process function must be enabled at a time."); + } + bool process3Prong = doprocessDataLcK0s || doprocessDataDplusK0s; + bool process2Prong = doprocessDataD0K0s; + bool processDstar = doprocessDataDstarK0s; + + // setup columnpolicy for binning + colBinningMult = {{AxisBinning.mixingBinVztx, AxisBinning.mixingBinMult}, true}; + colBinningMultPercentile = {{AxisBinning.mixingBinVztx, AxisBinning.mixingBinMultPercentile}, true}; + colBinningMultMultPercentile = {{AxisBinning.mixingBinVztx, AxisBinning.mixingBinMult, AxisBinning.mixingBinMultPercentile}, true}; + eventHisto.init(®istry); + v0HistoPartOne.init(®istry, AxisBinning.binmultTempFit, AxisBinning.dummy, AxisBinning.binpTV0, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.binTempFitVarV0, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, isMc, v0Sel.pdgCodeV0); + v0HistoPartOneSelected.init(®istry, AxisBinning.binmultTempFit, AxisBinning.dummy, AxisBinning.binpTV0, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.binTempFitVarV0, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, isMc, v0Sel.pdgCodeV0); + posChildHistos.init(®istry, AxisBinning.binmultTempFit, AxisBinning.dummy, AxisBinning.binpTV0Child, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.binTempFitVarV0Child, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, false, 0); + negChildHistos.init(®istry, AxisBinning.binmultTempFit, AxisBinning.dummy, AxisBinning.binpTV0Child, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.binTempFitVarV0Child, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, AxisBinning.dummy, false, 0); + sameEventCont.init(®istry, + AxisBinning.binkstar, AxisBinning.binpTV0, AxisBinning.binkT, AxisBinning.binmT, AxisBinning.mixingBinMult, AxisBinning.mixingBinMultPercentile, + AxisBinning.bin4Dkstar, AxisBinning.bin4DmT, AxisBinning.bin4DMult, AxisBinning.bin4DmultPercentile, + isMc, use4D, extendedPlots, + highkstarCut, + smearingByOrigin, AxisBinning.binInvMass); + + sameEventCont.setPDGCodes(v0Sel.pdgCodeV0, charmSel.charmHadPDGCode); + mixedEventCont.init(®istry, + AxisBinning.binkstar, AxisBinning.binpTV0, AxisBinning.binkT, AxisBinning.binmT, AxisBinning.mixingBinMult, AxisBinning.mixingBinMultPercentile, + AxisBinning.bin4Dkstar, AxisBinning.bin4DmT, AxisBinning.bin4DMult, AxisBinning.bin4DmultPercentile, + isMc, use4D, extendedPlots, + highkstarCut, + smearingByOrigin, AxisBinning.binInvMass); + + mixedEventCont.setPDGCodes(v0Sel.pdgCodeV0, charmSel.charmHadPDGCode); + registryMixQa.add("MixingQA/hSECollisionBins", "; bin; Entries", kTH1F, {{120, -0.5, 119.5}}); + registryMixQa.add("MixingQA/hSECollisionPool", "; Vz (cm); Mul", kTH2F, {{100, -10, 10}, {200, 0, 200}}); + registryMixQa.add("MixingQA/hMECollisionBins", "; bin; Entries", kTH1F, {{120, -0.5, 119.5}}); + registryCharmHadronQa.add("CharmHadronQA/hPtVsMass", "; #it{p}_{T} (GeV/#it{c}); inv. mass (GeV/#it{c}^{2})", kTH2F, {AxisBinning.binpTCharm, AxisBinning.binInvMass}); + + if (useCPR.value && process3Prong) { + pairCleaner3Prong.init(®istry); + pairCloseRejectionSE3Prong.init(®istry, ®istry, cprDeltaPhiMax.value, cprDeltaEtaMax.value, cprPlotPerRadii.value, 1); + pairCloseRejectionME3Prong.init(®istry, ®istry, cprDeltaPhiMax.value, cprDeltaEtaMax.value, cprPlotPerRadii.value, 2); + } else if (useCPR.value && process2Prong) { + pairCleaner2Prong.init(®istry); + pairCloseRejectionSE2Prong.init(®istry, ®istry, cprDeltaPhiMax.value, cprDeltaEtaMax.value, cprPlotPerRadii.value, 1); + pairCloseRejectionME2Prong.init(®istry, ®istry, cprDeltaPhiMax.value, cprDeltaEtaMax.value, cprPlotPerRadii.value, 2); + } else if (useCPR.value && processDstar) { + pairCleanerDstar.init(®istry); + pairCloseRejectionSEDstar.init(®istry, ®istry, cprDeltaPhiMax.value, cprDeltaEtaMax.value, cprPlotPerRadii.value, 1); + pairCloseRejectionMEDstar.init(®istry, ®istry, cprDeltaPhiMax.value, cprDeltaEtaMax.value, cprPlotPerRadii.value, 2); + } + } + + template + void fillCollision(CollisionType const& col) + { + registryMixQa.fill(HIST("MixingQA/hSECollisionBins"), colBinningMult.getBin({col.posZ(), col.multNtr()})); + registryMixQa.fill(HIST("MixingQA/hSECollisionPool"), col.posZ(), col.multNtr()); + } + + /// Compute the charm hadron candidates mass with the daughter masses + /// assumes the candidate is either a D+ or Λc+ or D0 or Dstar + template + float getCharmHadronMass(const Candidate& cand, bool ReturnDaughMass = false) + { + float invMass = 0.0f; + if constexpr (Channel == DecayChannel::LcToPKPi) { + if (cand.candidateSelFlag() == 1) { + invMass = cand.m(std::array{MassProton, MassKPlus, MassPiPlus}); + return invMass; + } + invMass = cand.m(std::array{MassPiPlus, MassKPlus, MassProton}); + return invMass; + } else if constexpr (Channel == DecayChannel::DplusToPiKPi) { // D+ → π K π (PDG: 411) + invMass = cand.m(std::array{MassPiPlus, MassKPlus, MassPiPlus}); + return invMass; + } else if constexpr (Channel == DecayChannel::D0ToPiK) { // D0 → π K (PDG: 421) + if (cand.candidateSelFlag() == 1) { + invMass = cand.m(std::array{MassPiPlus, MassKPlus}); + return invMass; + } else { + invMass = cand.m(std::array{MassKPlus, MassPiPlus}); + return invMass; + } + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { // D* → D0π (PDG: 413) + float mDstar = 0.f; + float mD0 = 0.f; + if (cand.charge() > 0.f) { + mDstar = cand.m(std::array{MassPiPlus, MassKPlus, MassPiPlus}); + mD0 = cand.mDaughD0(std::array{MassPiPlus, MassKPlus}); + } else { + mDstar = cand.m(std::array{MassKPlus, MassPiPlus, MassPiPlus}); + mD0 = cand.mDaughD0(std::array{MassKPlus, MassPiPlus}); + } + if (ReturnDaughMass) { + return mD0; + } else { + return mDstar - mD0; + } + } + // Add more channels as needed + return 0.f; + } + + /// This function processes the same event and takes care of all the histogramming + template + void doSameEvent(CandType& sliceCharmHad, PartitionType& sliceV01, FDParticles const& femtoParts, Collision const& col) + { + fillCollision(col); + processType = 1; // for same event + for (auto const& [p1, p2] : combinations(CombinationsFullIndexPolicy(sliceV01, sliceCharmHad))) { + + if constexpr (Channel == DecayChannel::D0ToPiK) { + if (p1.childrenIds()[0] == p2.prong0Id() || p1.childrenIds()[0] == p2.prong1Id() || p1.childrenIds()[1] == p2.prong0Id() || p1.childrenIds()[1] == p2.prong1Id()) + continue; + + if (useCPR.value) { + if (pairCloseRejectionSE2Prong.isClosePair(p1, p2, femtoParts, col.magField())) { + continue; + } + } + if (!pairCleaner2Prong.isCleanPair(p1, p2, femtoParts)) { + continue; + } + } else if constexpr (Channel == DecayChannel::LcToPKPi || Channel == DecayChannel::DplusToPiKPi) { + if (p1.childrenIds()[0] == p2.prong0Id() || p1.childrenIds()[0] == p2.prong1Id() || p1.childrenIds()[0] == p2.prong2Id() || p1.childrenIds()[1] == p2.prong0Id() || p1.childrenIds()[1] == p2.prong1Id() || p1.childrenIds()[1] == p2.prong2Id()) + continue; + if (useCPR.value) { + if (pairCloseRejectionSE3Prong.isClosePair(p1, p2, femtoParts, col.magField())) { + continue; + } + } + if (!pairCleaner3Prong.isCleanPair(p1, p2, femtoParts)) { + continue; + } + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { + if (p1.childrenIds()[0] == p2.prong0Id() || p1.childrenIds()[0] == p2.prong1Id() || p1.childrenIds()[0] == p2.prong2Id() || p1.childrenIds()[1] == p2.prong0Id() || p1.childrenIds()[1] == p2.prong1Id() || p1.childrenIds()[1] == p2.prong2Id()) + continue; + if (useCPR.value) { + if (pairCloseRejectionSEDstar.isClosePair(p1, p2, femtoParts, col.magField())) { + continue; + } + } + + if (!pairCleanerDstar.isCleanPair(p1, p2, femtoParts)) { + continue; + } + } + // v0 daughters selection + const auto& posChild = femtoParts.iteratorAt(p1.index() - 2); + const auto& negChild = femtoParts.iteratorAt(p1.index() - 1); + if (((posChild.cut() & v0Sel.childPosCutBit) == v0Sel.childPosCutBit) && + ((posChild.pidcut() & v0Sel.childPosTPCBit) == v0Sel.childPosTPCBit) && + ((negChild.cut() & v0Sel.childNegCutBit) == v0Sel.childNegCutBit) && + ((negChild.pidcut() & v0Sel.childNegTPCBit) == v0Sel.childNegTPCBit)) { + + v0HistoPartOneSelected.fillQA(p1, static_cast(confTempFitVarMomentum.value), col.multNtr(), col.multV0M()); + posChildHistos.fillQA(posChild, static_cast(confTempFitVarMomentum.value), col.multNtr(), col.multV0M()); + negChildHistos.fillQA(negChild, static_cast(confTempFitVarMomentum.value), col.multNtr(), col.multV0M()); + } else { + continue; + } + + float kstar = FemtoDreamMath::getkstar(p1, massOne, p2, massTwo); + if (kstar > highkstarCut) { + continue; + } + float invMass = getCharmHadronMass(p2); + + if (invMass < charmSel.charmHadMinInvMass || invMass > charmSel.charmHadMaxInvMass) { + continue; + } + + if (p2.pt() < charmSel.charmHadMinPt || p2.pt() > charmSel.charmHadMaxPt) { + continue; + } + + float invMassV0 = 0.f; + if (p1.sign() > 0) + invMassV0 = p1.mLambda(); + else + invMassV0 = p1.mAntiLambda(); + + float chargeV0 = 0.; + if ((p1.cut() & p1.sign()) == CutBitChargePositive) { + chargeV0 = PositiveCharge; + } else { + chargeV0 = NegativeCharge; + } + int pairSign = 0; + if (chargeV0 == p2.charge()) { + pairSign = LikeSignPair; + } else { + pairSign = UnLikeSignPair; + } + + int charmHadMc = 0; + int originType = 0; + if constexpr (IsMc) { + charmHadMc = p2.flagMc(); + originType = p2.originMcRec(); + } + + rowFemtoResultPairs( + invMass, + p2.pt(), + p1.pt(), + p2.bdtBkg(), + p2.bdtPrompt(), + p2.bdtFD(), + kstar, + FemtoDreamMath::getkT(p1, massOne, p2, massTwo), + FemtoDreamMath::getmT(p1, massOne, p2, massTwo), + col.multNtr(), + col.multV0M(), + p2.charge(), + pairSign, + invMassV0, // For the moment, use the delta-mass column to represent the V0 invariant mass + processType, + charmHadMc, + originType); + + sameEventCont.setPair(p1, p2, col.multNtr(), col.multV0M(), use4D, extendedPlots, smearingByOrigin); + } + } + + template + void doMixedEvent(CollisionType const& cols, PartitionType1& charms, PartitionType2& v0s, FDParticles const& femtoParts, BinningType policy) + { + processType = 2; // for mixed event + // Mixed events that contain the pair of interest + Partition partitionMaskedCol1 = (aod::femtodreamcollision::bitmaskTrackOne & bitMask) == bitMask; + partitionMaskedCol1.bindTable(cols); + + Partition partitionMaskedCol2 = (aod::femtodreamcollision::bitmaskTrackTwo & bitMask) == bitMask; + partitionMaskedCol2.bindTable(cols); + + for (auto const& [collision1, collision2] : combinations(soa::CombinationsBlockFullIndexPolicy(policy, mixSetting.mixingDepth, -1, *partitionMaskedCol1.mFiltered, *partitionMaskedCol2.mFiltered))) { + // make sure that tracks in the same events are not mixed + if (collision1.globalIndex() == collision2.globalIndex()) { + continue; + } + + const int multiplicityCol = collision1.multNtr(); + registryMixQa.fill(HIST("MixingQA/hMECollisionBins"), colBinningMult.getBin({collision1.posZ(), multiplicityCol})); + + auto sliceV01 = v0s->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision1.globalIndex(), cache); + auto sliceCharmHad = charms->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision2.globalIndex(), cache); + for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(sliceV01, sliceCharmHad))) { + + if constexpr (Channel == DecayChannel::D0ToPiK) { + + if (useCPR.value) { + if (pairCloseRejectionME2Prong.isClosePair(p1, p2, femtoParts, collision1.magField())) { + continue; + } + } + + if (!pairCleaner2Prong.isCleanPair(p1, p2, femtoParts)) { + continue; + } + } + + if constexpr (Channel == DecayChannel::DplusToPiKPi || Channel == DecayChannel::LcToPKPi) { + + if (useCPR.value) { + if (pairCloseRejectionME3Prong.isClosePair(p1, p2, femtoParts, collision1.magField())) { + continue; + } + } + + if (!pairCleaner3Prong.isCleanPair(p1, p2, femtoParts)) { + continue; + } + } + + if constexpr (Channel == DecayChannel::DstarToD0Pi) { + + if (useCPR.value) { + if (pairCloseRejectionME3Prong.isClosePair(p1, p2, femtoParts, collision1.magField())) { + continue; + } + } + + if (!pairCleanerDstar.isCleanPair(p1, p2, femtoParts)) { + continue; + } + } + + float kstar = FemtoDreamMath::getkstar(p1, massOne, p2, massTwo); + if (kstar > highkstarCut) { + continue; + } + + float invMass = getCharmHadronMass(p2); + + if (invMass < charmSel.charmHadMinInvMass || invMass > charmSel.charmHadMaxInvMass) { + continue; + } + + if (p2.pt() < charmSel.charmHadMinPt || p2.pt() > charmSel.charmHadMaxPt) { + continue; + } + + float invMassV0 = 0.f; + if (p1.sign() > 0) + invMassV0 = p1.mLambda(); + else + invMassV0 = p1.mAntiLambda(); + + float chargeV0 = 0.; + if ((p1.cut() & p1.sign()) == CutBitChargePositive) { + chargeV0 = PositiveCharge; + } else { + chargeV0 = NegativeCharge; + } + int pairSign = 0; + if (chargeV0 == p2.charge()) { + pairSign = LikeSignPair; + } else { + pairSign = UnLikeSignPair; + } + + int charmHadMc = 0; + int originType = 0; + if constexpr (IsMc) { + charmHadMc = p2.flagMc(); + originType = p2.originMcRec(); + } + + rowFemtoResultPairs( + invMass, + p2.pt(), + p1.pt(), + p2.bdtBkg(), + p2.bdtPrompt(), + p2.bdtFD(), + kstar, + FemtoDreamMath::getkT(p1, massOne, p2, massTwo), + FemtoDreamMath::getmT(p1, massOne, p2, massTwo), + collision1.multNtr(), + collision1.multV0M(), + p2.charge(), + pairSign, + invMassV0, + processType, + charmHadMc, + originType); + + mixedEventCont.setPair(p1, p2, collision1.multNtr(), collision1.multV0M(), use4D, extendedPlots, smearingByOrigin); + } + } + } + template + void fillTables(const CollType& col, + const V0Type& sliceV01, + const CandType& sliceCharmHad, + const FDParticles& femtoParts) + { + int64_t timeStamp = -999; + + // ---- Fill Charm-Hadron Table ---- + for (auto const& part : sliceCharmHad) { + float invMass = getCharmHadronMass(part); + registryCharmHadronQa.fill(HIST("CharmHadronQA/hPtVsMass"), part.pt(), invMass); + timeStamp = part.timeStamp(); + + if constexpr (Channel == DecayChannel::DplusToPiKPi || Channel == DecayChannel::LcToPKPi) { + + rowFemtoResultCharm3Prong( + col.globalIndex(), + timeStamp, + invMass, + part.pt(), + part.eta(), + part.phi(), + part.prong0Id(), + part.prong1Id(), + part.prong2Id(), + part.charge(), + part.bdtBkg(), + part.bdtPrompt(), + part.bdtFD()); + } else if constexpr (Channel == DecayChannel::D0ToPiK) { + rowFemtoResultCharm2Prong( + col.globalIndex(), + timeStamp, + invMass, + part.pt(), + part.eta(), + part.phi(), + part.prong0Id(), + part.prong1Id(), + part.charge(), + part.bdtBkg(), + part.bdtPrompt(), + part.bdtFD()); + } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { + float invMassD0 = getCharmHadronMass(part, true); + rowFemtoResultCharmDstar( + col.globalIndex(), + timeStamp, + invMass, + invMassD0, + part.pt(), + part.eta(), + part.phi(), + part.prong0Id(), + part.prong1Id(), + part.prong2Id(), + part.charge(), + part.bdtBkg(), + part.bdtPrompt(), + part.bdtFD()); + } + } + + // ---- Fill V0 Table ---- + for (auto const& part : sliceV01) { + + // v0 daughters selection + const auto& posChild = femtoParts.iteratorAt(part.index() - 2); + const auto& negChild = femtoParts.iteratorAt(part.index() - 1); + if (((posChild.cut() & v0Sel.childPosCutBit) == v0Sel.childPosCutBit) && + ((posChild.pidcut() & v0Sel.childPosTPCBit) == v0Sel.childPosTPCBit) && + ((negChild.cut() & v0Sel.childNegCutBit) == v0Sel.childNegCutBit) && + ((negChild.pidcut() & v0Sel.childNegTPCBit) == v0Sel.childNegTPCBit)) { + v0HistoPartOne.fillQA(part, static_cast(confTempFitVarMomentum.value), col.multNtr(), col.multV0M()); + } else { + continue; + } + + std::vector childIds = {part.childrenIds()[0], part.childrenIds()[1]}; + timeStamp = part.timeStamp(); + rowFemtoResultV0( + col.globalIndex(), + timeStamp, + part.pt(), + part.eta(), + part.phi(), + childIds, + part.sign(), + part.mLambda(), + part.mAntiLambda(), + part.tpcNClsFound(), + part.tpcNClsFindable(), + part.tpcNClsCrossedRows(), + part.tpcNSigmaPr(), + part.tofNSigmaPr()); + } + + // ---- Fill Collision Table ---- + if (sliceCharmHad.size() > 0 || sliceV01.size() > 0) { + rowFemtoResultColl( + col.globalIndex(), + timeStamp, + col.posZ(), + col.multNtr()); + } + } + + void processDataLcK0s(FilteredCollisions const& cols, + FilteredFDParticles const& parts, + FilteredCharmCand3Prongs const&) + { + for (const auto& col : cols) { + eventHisto.fillQA(col); + auto sliceCharmHad = partitionCharmHadron3Prong->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + if (fillTableWithCharm.value && sliceCharmHad.size() == 0) { + continue; + } + if (v0Sel.pdgCodeV0 == kLambda0) { + auto sliceV0 = partitionLambda->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + + } else if (v0Sel.pdgCodeV0 == kK0Short) { + auto sliceV0 = partitionK0Short->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + } else { + LOG(fatal) << "Unsupported V0 PDG: " << v0Sel.pdgCodeV0 << " (allowed: 3122, 310)"; + } + } + + if (mixSetting.doMixEvent) { + switch (mixSetting.mixingBinPolicy) { + case femtodreamcollision::kMult: + doMixedEvent(cols, partitionCharmHadron3Prong, partitionK0Short, parts, colBinningMult); + break; + case femtodreamcollision::kMultPercentile: + doMixedEvent(cols, partitionCharmHadron3Prong, partitionK0Short, parts, colBinningMultPercentile); + break; + case femtodreamcollision::kMultMultPercentile: + doMixedEvent(cols, partitionCharmHadron3Prong, partitionK0Short, parts, colBinningMultMultPercentile); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + PROCESS_SWITCH(HfTaskCharmHadronsV0FemtoDream, processDataLcK0s, "Enable processing LcToPKPi and V0 correlation", false); + + void processDataDplusK0s(FilteredCollisions const& cols, + FilteredFDParticles const& parts, + FilteredCharmCand3Prongs const&) + { + for (const auto& col : cols) { + eventHisto.fillQA(col); + auto sliceCharmHad = partitionCharmHadron3Prong->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + if (fillTableWithCharm.value && sliceCharmHad.size() == 0) { + continue; + } + if (v0Sel.pdgCodeV0 == kLambda0) { + auto sliceV0 = partitionLambda->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + + } else if (v0Sel.pdgCodeV0 == kK0Short) { + auto sliceV0 = partitionK0Short->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + } else { + LOG(fatal) << "Unsupported V0 PDG: " << v0Sel.pdgCodeV0 << " (allowed: 3122, 310)"; + } + } + + if (mixSetting.doMixEvent) { + switch (mixSetting.mixingBinPolicy) { + case femtodreamcollision::kMult: + doMixedEvent(cols, partitionCharmHadron3Prong, partitionK0Short, parts, colBinningMult); + break; + case femtodreamcollision::kMultPercentile: + doMixedEvent(cols, partitionCharmHadron3Prong, partitionK0Short, parts, colBinningMultPercentile); + break; + case femtodreamcollision::kMultMultPercentile: + doMixedEvent(cols, partitionCharmHadron3Prong, partitionK0Short, parts, colBinningMultMultPercentile); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + PROCESS_SWITCH(HfTaskCharmHadronsV0FemtoDream, processDataDplusK0s, "Enable processing DplusToPiKPi and V0 correlation", false); + + void processDataD0K0s(FilteredCollisions const& cols, + FilteredFDParticles const& parts, + FilteredCharmCand2Prongs const&) + { + for (const auto& col : cols) { + eventHisto.fillQA(col); + auto sliceCharmHad = partitionCharmHadron2Prong->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + if (fillTableWithCharm.value && sliceCharmHad.size() == 0) { + continue; + } + if (v0Sel.pdgCodeV0 == kLambda0) { + auto sliceV0 = partitionLambda->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + + } else if (v0Sel.pdgCodeV0 == kK0Short) { + auto sliceV0 = partitionK0Short->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + } else { + LOG(fatal) << "Unsupported V0 PDG: " << v0Sel.pdgCodeV0 << " (allowed: 3122, 310)"; + } + } + + if (mixSetting.doMixEvent) { + switch (mixSetting.mixingBinPolicy) { + case femtodreamcollision::kMult: + doMixedEvent(cols, partitionCharmHadron2Prong, partitionK0Short, parts, colBinningMult); + break; + case femtodreamcollision::kMultPercentile: + doMixedEvent(cols, partitionCharmHadron2Prong, partitionK0Short, parts, colBinningMultPercentile); + break; + case femtodreamcollision::kMultMultPercentile: + doMixedEvent(cols, partitionCharmHadron2Prong, partitionK0Short, parts, colBinningMultMultPercentile); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + PROCESS_SWITCH(HfTaskCharmHadronsV0FemtoDream, processDataD0K0s, "Enable processing D0ToPiK and V0 correlation", false); + + void processDataDstarK0s(FilteredCollisions const& cols, + FilteredFDParticles const& parts, + FilteredCharmCandDstars const&) + { + for (const auto& col : cols) { + eventHisto.fillQA(col); + auto sliceCharmHad = partitionCharmHadronDstar->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + if (fillTableWithCharm.value && sliceCharmHad.size() == 0) { + continue; + } + if (v0Sel.pdgCodeV0 == kLambda0) { + auto sliceV0 = partitionLambda->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + + } else if (v0Sel.pdgCodeV0 == kK0Short) { + auto sliceV0 = partitionK0Short->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + fillTables(col, sliceV0, sliceCharmHad, parts); + + if (sliceCharmHad.size() > 0 && sliceV0.size() > 0) { + doSameEvent(sliceCharmHad, sliceV0, parts, col); + } + } else { + LOG(fatal) << "Unsupported V0 PDG: " << v0Sel.pdgCodeV0 << " (allowed: 3122, 310)"; + } + } + + if (mixSetting.doMixEvent) { + switch (mixSetting.mixingBinPolicy) { + case femtodreamcollision::kMult: + doMixedEvent(cols, partitionCharmHadronDstar, partitionK0Short, parts, colBinningMult); + break; + case femtodreamcollision::kMultPercentile: + doMixedEvent(cols, partitionCharmHadronDstar, partitionK0Short, parts, colBinningMultPercentile); + break; + case femtodreamcollision::kMultMultPercentile: + doMixedEvent(cols, partitionCharmHadronDstar, partitionK0Short, parts, colBinningMultMultPercentile); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + PROCESS_SWITCH(HfTaskCharmHadronsV0FemtoDream, processDataDstarK0s, "Enable processing DstarToD0Pi and V0 correlation", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}