diff --git a/PWGHF/DataModel/DerivedTables.h b/PWGHF/DataModel/DerivedTables.h index e3602a0ba63..3917627c547 100644 --- a/PWGHF/DataModel/DerivedTables.h +++ b/PWGHF/DataModel/DerivedTables.h @@ -356,6 +356,9 @@ DECLARE_SOA_COLUMN(PtProngPi1, ptProngPi1, float); DECLARE_SOA_COLUMN(RSecondaryVertex, rSecondaryVertex, float); //! distance of the secondary vertex from the z axis // D*± → D0(bar) π± DECLARE_SOA_COLUMN(SignProng1, signProng1, int8_t); +DECLARE_SOA_COLUMN(AbsEtaTrackMin, absEtaTrackMin, float); +DECLARE_SOA_COLUMN(NumItsClsMin, numItsClsMin, int); +DECLARE_SOA_COLUMN(NumTpcClsMin, numTpcClsMin, int); // TOF DECLARE_SOA_COLUMN(NSigTofKa0, nSigTofKa0, float); DECLARE_SOA_COLUMN(NSigTofKa1, nSigTofKa1, float); @@ -963,6 +966,9 @@ DECLARE_SOA_TABLE_STAGED(HfDstarPars, "HFDSTPAR", //! Table with candidate prope hf_cand_par::NSigTpcPi1, hf_cand_par::NSigTofPi1, hf_cand_par::NSigTpcTofPi1, + hf_cand_par::AbsEtaTrackMin, + hf_cand_par::NumItsClsMin, + hf_cand_par::NumTpcClsMin, o2::soa::Marker); DECLARE_SOA_TABLE_STAGED(HfDstarParD0s, "HFDSTPARD0", //! Table with candidate properties used for selection diff --git a/PWGHF/HFJ/CMakeLists.txt b/PWGHF/HFJ/CMakeLists.txt index 62f88c324cf..d1ccd071b54 100644 --- a/PWGHF/HFJ/CMakeLists.txt +++ b/PWGHF/HFJ/CMakeLists.txt @@ -7,4 +7,6 @@ # # 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. \ No newline at end of file +# or submit itself to any jurisdiction. + +add_subdirectory(Tasks) \ No newline at end of file diff --git a/PWGHF/HFJ/Tasks/CMakeLists.txt b/PWGHF/HFJ/Tasks/CMakeLists.txt new file mode 100644 index 00000000000..acdf2c92eb1 --- /dev/null +++ b/PWGHF/HFJ/Tasks/CMakeLists.txt @@ -0,0 +1,17 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2physics_add_dpl_workflow(task-dstar-polarisation-in-jet + SOURCES taskDstarPolarisationInJet.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/PWGJE/Core) \ No newline at end of file diff --git a/PWGHF/HFJ/Tasks/taskDstarPolarisationInJet.cxx b/PWGHF/HFJ/Tasks/taskDstarPolarisationInJet.cxx new file mode 100644 index 00000000000..d8119be808d --- /dev/null +++ b/PWGHF/HFJ/Tasks/taskDstarPolarisationInJet.cxx @@ -0,0 +1,980 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file taskDstarPolarisationInJet.cxx +/// \brief Analysis task for Dstar spin alignment in jet (rho00 vs. z||) +/// +/// \author F. Grosa (CERN) fabrizio.grosa@cern.ch +/// \author Mingze Li (CCNU) Mingze.li@cern.ch + +#include "PWGHF/Core/DecayChannels.h" +#include "PWGJE/Core/FastJetUtilities.h" +#include "PWGJE/Core/JetFinder.h" +#include "PWGJE/Core/JetHFUtilities.h" +#include "PWGJE/Core/JetUtilities.h" +#include "PWGJE/DataModel/Jet.h" + +#include +#include +#include + +#include +#include +#include + +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace o2 +{ +namespace charm_polarisation +{ +enum CosThetaStarType : uint8_t { + Helicity = 0, + JetAxis, + Production, + Beam, + Random, + NTypes +}; +enum DecayChannel : uint8_t { + DstarToDzeroPi = 0, + NChannels +}; +} // namespace charm_polarisation +} // namespace o2 + +struct HfTaskDstarPolarisationInJet { + + Configurable selectionFlagDstarToD0Pi{"selectionFlagDstarToD0Pi", true, "Selection Flag for D* decay to D0 Pi"}; + Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; + + /// activate rotational background + Configurable nBkgRotations{"nBkgRotations", 0, "Number of rotated copies (background) per each original candidate"}; + Configurable minRotAngleMultByPi{"minRotAngleMultByPi", 5. / 6, "Minimum angle rotation for track rotation, to be multiplied by pi"}; + Configurable maxRotAngleMultByPi{"maxRotAngleMultByPi", 7. / 6, "Maximum angle rotation for track rotation, to be multiplied by pi"}; + + // activate study of systematic uncertainties of tracking + Configurable activateTrackingSys{"activateTrackingSys", false, "Activate the study of systematic uncertainties of tracking"}; + + /// output THnSparses + Configurable activateTHnSparseCosThStarHelicity{"activateTHnSparseCosThStarHelicity", true, "Activate the THnSparse with cosThStar w.r.t. helicity axis"}; + Configurable activateTHnSparseCosThStarJetAxis{"activateTHnSparseCosThStarJetAxis", true, "Activate the THnSparse with cosThStar w.r.t. production axis"}; + Configurable activateTHnSparseCosThStarProduction{"activateTHnSparseCosThStarProduction", true, "Activate the THnSparse with cosThStar w.r.t. production axis"}; + Configurable activatePartRecoDstar{"activatePartRecoDstar", false, "Activate the study of partly reconstructed D*+ -> D0 (-> KPiPi0) Pi decays"}; + + /// Application of rapidity cut for reconstructed candidates + Configurable absRapidityCandMax{"absRapidityCandMax", 999.f, "Max. value of reconstructed candidate rapidity (abs. value)"}; + + float invMassMin{0.f}; + float invMassMax{1000.f}; + float bkgRotationAngleStep{0.f}; + uint8_t nMassHypos{0u}; + struct HfHistoInput { + float invMassCharmHad; + float ptCharmHad; + float rapCharmHad; + float invMassD0; + float cosThetaStar; + const std::array& outputMl; + int isRotatedCandidate; + int8_t origin; + float ptBhadMother; + float absEtaMin; + int numItsClsMin; + int numTpcClsMin; + int8_t nMuons; + float zParallel; + float jetPt; + }; + + // Tables for MC jet matching + using DstarJets = soa::Join; + using JetMCDTable = soa::Join; + using TracksWithExtra = soa::Join; + + // slices for accessing proper HF mcdjets collision associated to mccollisions + Preslice dstarMCDJetsPerCollisionPreslice = aod::jet::collisionId; + Preslice dstarJetsPerCollision = aod::jet::collisionId; + PresliceUnsorted collisionsPerMCCollisionPreslice = aod::jmccollisionlb::mcCollisionId; + + ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {200, 0.139f, 0.179f}, "#it{M} (GeV/#it{c}^{2})"}; // o2-linter: disable=pdg/explicit-mass (false positive) + ConfigurableAxis configThnAxisPt{"configThnAxisPt", {100, 0.f, 100.f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis configThnAxisY{"configThnAxisY", {20, -1.f, 1.f}, "#it{y}"}; + ConfigurableAxis configThnAxisCosThetaStar{"configThnAxisCosThetaStar", {20, -1.f, 1.f}, "cos(#vartheta_{helicity})"}; + ConfigurableAxis configThnAxisMlBkg{"configThnAxisMlBkg", {100, 0.f, 1.f}, "ML bkg"}; + ConfigurableAxis configThnAxisInvMassD0{"configThnAxisInvMassD0", {250, 1.65f, 2.15f}, "#it{M}(D^{0}) (GeV/#it{c}^{2})"}; // only for D*+ + // ConfigurableAxis configThnAxisMlPrompt{"configThnAxisMlPrompt", {100, 0.f, 1.f}, "ML prompt"}; + ConfigurableAxis configThnAxisMlNonPrompt{"configThnAxisMlNonPrompt", {100, 0.f, 1.f}, "ML non-prompt"}; + ConfigurableAxis configThnAxisNumCandidates{"configThnAxisNumCandidates", {1, -0.5f, 0.5f}, "num candidates"}; + ConfigurableAxis configThnAxisPtB{"configThnAxisPtB", {3000, 0.f, 300.f}, "#it{p}_{T}(B mother) (GeV/#it{c})"}; + ConfigurableAxis configThnAxisAbsEtaTrackMin{"configThnAxisAbsEtaTrackMin", {3, 0.f, 0.3f}, "min |#it{#eta_{track}}|"}; + ConfigurableAxis configThnAxisNumItsClsMin{"configThnAxisNumItsClsMin", {4, 3.5f, 7.5f}, "min #it{N}_{cls ITS}"}; + ConfigurableAxis configThnAxisNumTpcClsMin{"configThnAxisNumTpcClsMin", {3, 79.5f, 140.5f}, "min #it{N}_{cls TPC}"}; + ConfigurableAxis configThnAxisCharge{"configThnAxisCharge", {2, -2.f, 2.f}, "electric charge"}; + ConfigurableAxis configThnAxisProjection{"configThnAxisProjection", {300, 0.f, 10.f}, "z^{D^{*+},jet}_{||}"}; + ConfigurableAxis configThnAxisJetPt{"configThnAxisJetPt", {100, 0.f, 100.f}, "#it{p}_{T} (GeV/#it{c})"}; + + HistogramRegistry registry{"registry", {}}; + + std::vector eventSelectionBits; + + void init(InitContext&) + { + /// check process functions + if (doprocessDstar + doprocessDstarWithMl + doprocessDstarMc + doprocessDstarMcWithMl > 1) { + LOGP(fatal, "Only one process function should be enabled at a time, please check your configuration"); + } + if (doprocessDstar + doprocessDstarWithMl + doprocessDstarMc + doprocessDstarMcWithMl == 0) { + LOGP(fatal, "No process function enabled"); + } + + // initialise event selection: + eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(eventSelections.value); + + /// check output THnSparses + std::array sparses = {activateTHnSparseCosThStarHelicity, activateTHnSparseCosThStarJetAxis, activateTHnSparseCosThStarProduction}; + if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) { + LOGP(fatal, "No output THnSparses enabled"); + } + if (activateTHnSparseCosThStarHelicity) { + LOGP(info, "THnSparse with cosThStar w.r.t. helicity axis active."); + } + if (activateTHnSparseCosThStarJetAxis) { + LOGP(info, "THnSparse with cosThStar w.r.t. jet axis active."); + } + if (activateTHnSparseCosThStarProduction) { + LOGP(info, "THnSparse with cosThStar w.r.t. production axis active."); + } + + if (activatePartRecoDstar && !(doprocessDstarMc || doprocessDstarMcWithMl)) { + LOGP(fatal, "Check on partly reconstructed D* mesons only possible for processDstarMc and processDstarMcWithMl"); + } + + // check bkg rotation for MC (not supported currently) + if (nBkgRotations > 0 && (doprocessDstarMc || doprocessDstarMcWithMl)) { + LOGP(fatal, "No background rotation supported for MC."); + } + + bkgRotationAngleStep = (nBkgRotations > 1) ? (maxRotAngleMultByPi - minRotAngleMultByPi) * constants::math::PI / (nBkgRotations - 1) : 0.; + + const AxisSpec thnAxisInvMass{configThnAxisInvMass, "#it{M} (GeV/#it{c}^{2})"}; + const AxisSpec thnAxisInvMassD0{configThnAxisInvMassD0, "#it{M}(D^{0}) (GeV/#it{c}^{2})"}; + const AxisSpec thnAxisPt{configThnAxisPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec thnAxisY{configThnAxisY, "#it{y}"}; + const AxisSpec thnAxisCosThetaStar{configThnAxisCosThetaStar, "cos(#vartheta)"}; + const AxisSpec thnAxisMlBkg{configThnAxisMlBkg, "ML bkg"}; + const AxisSpec thnAxisMlNonPrompt{configThnAxisMlNonPrompt, "ML non-prompt"}; + const AxisSpec thnAxisIsRotatedCandidate{2, -0.5f, 1.5f, "rotated bkg"}; + const AxisSpec thnAxisNumcandidates{configThnAxisNumCandidates, "num candidates"}; + const AxisSpec thnAxisPtB{configThnAxisPtB, "#it{p}_{T}(B mother) (GeV/#it{c})"}; + const AxisSpec thnAxisDausAcc{2, -0.5f, 1.5f, "daughters in acceptance"}; + const AxisSpec thnAxisDauToMuons{4, -0.5f, 3.5f, "daughters decayed to muons"}; + const AxisSpec thnAxisAbsEtaTrackMin{configThnAxisAbsEtaTrackMin, "min |#it{#eta_{track}}|"}; + const AxisSpec thnAxisNumItsClsMin{configThnAxisNumItsClsMin, "min #it{N}_{cls ITS}"}; + const AxisSpec thnAxisNumTpcClsMin{configThnAxisNumTpcClsMin, "min #it{N}_{cls TPC}"}; + const AxisSpec thnAxisCharge{configThnAxisCharge, "charge"}; + const AxisSpec thnAxisProjection{configThnAxisProjection, "z_{||}"}; + const AxisSpec thnAxisJetPt{configThnAxisJetPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec thnAxisSelFlag{2, -0.5f, 1.5f, "Sel flag"}; + + auto invMassBins = thnAxisInvMass.binEdges; + invMassMin = invMassBins.front(); + invMassMax = invMassBins.back(); + + std::vector thnRecDataAxes = {thnAxisInvMass, thnAxisPt, thnAxisY, thnAxisInvMassD0, thnAxisCosThetaStar}; + if (activateTrackingSys) { + thnRecDataAxes.insert(thnRecDataAxes.end(), {thnAxisAbsEtaTrackMin, thnAxisNumItsClsMin, thnAxisNumTpcClsMin}); + } + if (nBkgRotations > 0) { + thnRecDataAxes.push_back(thnAxisIsRotatedCandidate); + } + if (doprocessDstarWithMl || doprocessDstarMcWithMl) { + thnRecDataAxes.insert(thnRecDataAxes.end(), {thnAxisMlBkg, thnAxisMlNonPrompt}); + } + std::vector thnRecPromptAxes = thnRecDataAxes; + std::vector thnRecNonPromptAxes = thnRecDataAxes; + + if (doprocessDstar || doprocessDstarWithMl) { + thnRecDataAxes.insert(thnRecDataAxes.end(), {thnAxisProjection, thnAxisJetPt}); + // Data histos + if (activateTHnSparseCosThStarHelicity) { + registry.add("hHelicity", "THn for polarisation studies with cosThStar ", HistType::kTHnSparseF, thnRecDataAxes); + } + if (activateTHnSparseCosThStarProduction) { + registry.add("hProduction", "THn for polarisation studies with cosThStar w.r.t. production axis", HistType::kTHnSparseF, thnRecDataAxes); + } + if (activateTHnSparseCosThStarJetAxis) { + registry.add("hJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis", HistType::kTHnSparseF, thnRecDataAxes); + } + } else if (doprocessDstarMc || doprocessDstarMcWithMl) { + // MC Reco histos + thnRecPromptAxes.insert(thnRecPromptAxes.end(), {thnAxisDauToMuons, thnAxisProjection, thnAxisJetPt}); + thnRecNonPromptAxes.insert(thnRecNonPromptAxes.end(), {thnAxisPtB, thnAxisDauToMuons, thnAxisProjection, thnAxisJetPt}); + if (activateTHnSparseCosThStarHelicity) { + registry.add("hRecoPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores -- reco prompt signal", HistType::kTHnSparseF, thnRecPromptAxes); + registry.add("hRecoNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores -- reco non-prompt signal", HistType::kTHnSparseF, thnRecNonPromptAxes); + if (activatePartRecoDstar) { + registry.add("hPartRecoPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores -- partially reco prompt signal", HistType::kTHnSparseF, thnRecPromptAxes); + registry.add("hPartRecoNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis and BDT scores -- partially reco non-prompt signal", HistType::kTHnSparseF, thnRecNonPromptAxes); + } + } + if (activateTHnSparseCosThStarProduction) { + registry.add("hRecoPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores -- reco prompt signal", HistType::kTHnSparseF, thnRecPromptAxes); + registry.add("hRecoNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores -- reco non-prompt signal", HistType::kTHnSparseF, thnRecNonPromptAxes); + if (activatePartRecoDstar) { + registry.add("hPartRecoPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores -- partially reco prompt signal", HistType::kTHnSparseF, thnRecPromptAxes); + registry.add("hPartRecoNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis and BDT scores -- partially reco non-prompt signal", HistType::kTHnSparseF, thnRecNonPromptAxes); + } + } + if (activateTHnSparseCosThStarJetAxis) { + registry.add("hRecoPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis and BDT scores -- reco prompt signal", HistType::kTHnSparseF, thnRecPromptAxes); + registry.add("hRecoNonPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis and BDT scores -- reco non-prompt signal", HistType::kTHnSparseF, thnRecNonPromptAxes); + if (activatePartRecoDstar) { + registry.add("hPartRecoPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis and BDT scores -- partially reco prompt signal", HistType::kTHnSparseF, thnRecPromptAxes); + registry.add("hPartRecoNonPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis and BDT scores -- partially reco non-prompt signal", HistType::kTHnSparseF, thnRecNonPromptAxes); + } + } + + // MC Gen histos + std::vector thnGenPromptAxes = {thnAxisPt, thnAxisY, thnAxisCosThetaStar, thnAxisDausAcc, thnAxisCharge, thnAxisProjection, thnAxisJetPt}; + std::vector thnGenNonPromptAxes = {thnAxisPt, thnAxisY, thnAxisCosThetaStar, thnAxisPtB, thnAxisDausAcc, thnAxisCharge, thnAxisProjection, thnAxisJetPt}; + + if (activateTHnSparseCosThStarHelicity) { + registry.add("hGenPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis -- gen prompt signal", HistType::kTHnSparseF, thnGenPromptAxes); + registry.add("hGenNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis -- gen non-prompt signal", HistType::kTHnSparseF, thnGenNonPromptAxes); + if (activatePartRecoDstar && (doprocessDstarMc || doprocessDstarMcWithMl)) { + registry.add("hGenPartRecoPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis -- gen prompt partly reco signal", HistType::kTHnSparseF, thnGenPromptAxes); + registry.add("hGenPartRecoNonPromptHelicity", "THn for polarisation studies with cosThStar w.r.t. helicity axis -- gen non-prompt partly reco signal", HistType::kTHnSparseF, thnGenNonPromptAxes); + } + } + if (activateTHnSparseCosThStarProduction) { + registry.add("hGenPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis -- gen prompt signal", HistType::kTHnSparseF, thnGenPromptAxes); + registry.add("hGenNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis -- gen non-prompt signal", HistType::kTHnSparseF, thnGenNonPromptAxes); + if (activatePartRecoDstar && (doprocessDstarMc || doprocessDstarMcWithMl)) { + registry.add("hGenPartRecoPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis -- gen prompt partly reco signal", HistType::kTHnSparseF, thnGenPromptAxes); + registry.add("hGenPartRecoNonPromptProduction", "THn for polarisation studies with cosThStar w.r.t. production axis -- gen non-prompt partly reco signal", HistType::kTHnSparseF, thnGenNonPromptAxes); + } + } + if (activateTHnSparseCosThStarJetAxis) { + registry.add("hGenPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis -- gen prompt signal", HistType::kTHnSparseF, thnGenPromptAxes); + registry.add("hGenNonPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis -- gen non-prompt signal", HistType::kTHnSparseF, thnGenNonPromptAxes); + if (activatePartRecoDstar && (doprocessDstarMc || doprocessDstarMcWithMl)) { + registry.add("hGenPartRecoPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis -- gen prompt partly reco signal", HistType::kTHnSparseF, thnGenPromptAxes); + registry.add("hGenPartRecoNonPromptJetAxis", "THn for polarisation studies with cosThStar w.r.t. jet axis -- gen non-prompt partly reco signal", HistType::kTHnSparseF, thnGenNonPromptAxes); + } + } + } + // inv. mass hypothesis to loop over + nMassHypos = 1; + }; // end init + + /// \param invMassCharmHad is the invariant-mass of the candidate + /// \param ptCharmHad is the pt of the candidate + /// \param rapCharmHad is the rapidity of the candidate + /// \param invMassD0 is the invariant-mass of the D0 daugher (only for D*+) + /// \param cosThetaStar is the cosThetaStar of the candidate + /// \param outputMl is the array with ML output scores + /// \param isRotatedCandidate is a flag that keeps the info of the rotation of the candidate for bkg studies + /// \param origin is the MC origin + /// \param ptBhadMother is the pt of the b-hadron mother (only in case of non-prompt) + /// \param absEtaMin is the minimum absolute eta of the daughter tracks + /// \param numItsClsMin is the minimum number of ITS clusters of the daughter tracks + /// \param numTpcClsMin is the minimum number of TPC clusters of the daughter tracks + /// \param nMuons is the number of muons from daughter decays + /// \param isPartRecoDstar is a flag indicating if it is a partly reconstructed Dstar meson (MC only) + void fillRecoHistos(charm_polarisation::CosThetaStarType cosThetaStarType, bool withMl, bool doMc, bool isPartRecoDstar, HfHistoInput& recoHistoInput) + { + + if (cosThetaStarType == charm_polarisation::CosThetaStarType::Helicity) { // Helicity + if (!doMc) { // data + if (withMl) { // with ML + if (activateTrackingSys) { + if (nBkgRotations > 0) { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.isRotatedCandidate, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (nBkgRotations > 0) { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.isRotatedCandidate, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } else { // without ML + if (activateTrackingSys) { + if (nBkgRotations > 0) { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.isRotatedCandidate, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (nBkgRotations > 0) { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.isRotatedCandidate, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } else { // MC --> no distinction among channels, since rotational bkg not supported + if (withMl) { // with ML + if (recoHistoInput.origin == RecoDecay::OriginType::Prompt) { // prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } else { // non-prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } else { // without ML + if (recoHistoInput.origin == RecoDecay::OriginType::Prompt) { // prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + + } else { // non-prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptHelicity"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } + } + } else if (cosThetaStarType == charm_polarisation::CosThetaStarType::Production) { // Production + if (!doMc) { // data + if (withMl) { // with ML + if (activateTrackingSys) { + if (nBkgRotations > 0) { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.isRotatedCandidate, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (nBkgRotations > 0) { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.isRotatedCandidate, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } else { // without ML + if (activateTrackingSys) { + if (nBkgRotations > 0) { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.isRotatedCandidate, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (nBkgRotations > 0) { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.isRotatedCandidate, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } else { // MC --> no distinction among channels, since rotational bkg not supported + if (withMl) { // with ML + if (recoHistoInput.origin == RecoDecay::OriginType::Prompt) { // prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } else { // non-prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } else { // without ML + if (recoHistoInput.origin == RecoDecay::OriginType::Prompt) { // prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + + } else { // non-prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptProduction"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } + } + } else if (cosThetaStarType == charm_polarisation::CosThetaStarType::JetAxis) { // JetAxis + if (!doMc) { // data + if (withMl) { // with ML + if (activateTrackingSys) { + if (nBkgRotations > 0) { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.isRotatedCandidate, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (nBkgRotations > 0) { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.isRotatedCandidate, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } else { // without ML + if (activateTrackingSys) { + if (nBkgRotations > 0) { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.isRotatedCandidate, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (nBkgRotations > 0) { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.isRotatedCandidate, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } else { // MC --> no distinction among channels, since rotational bkg not supported + if (withMl) { // with ML + if (recoHistoInput.origin == RecoDecay::OriginType::Prompt) { // prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } else { // non-prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.outputMl[0], /*recoHistoInput.outputMl[1],*/ recoHistoInput.outputMl[2], recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } else { // without ML + if (recoHistoInput.origin == RecoDecay::OriginType::Prompt) { // prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + + } else { // non-prompt + if (activateTrackingSys) { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.absEtaMin, recoHistoInput.numItsClsMin, recoHistoInput.numTpcClsMin, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } else { + if (!isPartRecoDstar) { + registry.fill(HIST("hRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } else { + registry.fill(HIST("hPartRecoNonPromptJetAxis"), recoHistoInput.invMassCharmHad, recoHistoInput.ptCharmHad, recoHistoInput.rapCharmHad, recoHistoInput.invMassD0, recoHistoInput.cosThetaStar, recoHistoInput.ptBhadMother, recoHistoInput.nMuons, recoHistoInput.zParallel, recoHistoInput.jetPt); + } + } + } + } + } + } + } + + /// TODO: To be implemented + /// \param ptCharmHad is the pt of the particle + /// \param rapCharmHad is the rapidity of the particle + /// \param cosThetaStar is the cosThetaStar of the particle + /// \param origin is the MC origin + /// \param ptBhadMother is the pt of the b-hadron mother (only in case of non-prompt) + /// \param areDausInAcc is a flag indicating whether the daughters are in acceptance or not + /// \param isPartRecoDstar is a flag indicating if it is a partly reconstructed Dstar->D0pi->Kpipipi0 meson (MC only) + void fillGenHistos(charm_polarisation::CosThetaStarType cosThetaStarType, float ptCharmHad, float rapCharmHad, float cosThetaStar, int8_t origin, float ptBhadMother, bool areDausInAcc, int8_t charge, bool isPartRecoDstar, float zParallel, float jetPt) + { + if (cosThetaStarType == charm_polarisation::CosThetaStarType::Helicity) { // Helicity + if (origin == RecoDecay::OriginType::Prompt) { // prompt + if (!isPartRecoDstar) { + registry.fill(HIST("hGenPromptHelicity"), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt); + } else { + registry.fill(HIST("hGenPartRecoPromptHelicity"), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt); + } + } else { // non-prompt + if (!isPartRecoDstar) { + registry.fill(HIST("hGenNonPromptHelicity"), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt); + } else { + registry.fill(HIST("hGenPartRecoNonPromptHelicity"), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt); + } + } + } else if (cosThetaStarType == charm_polarisation::CosThetaStarType::Production) { // Production + if (origin == RecoDecay::OriginType::Prompt) { // prompt + if (!isPartRecoDstar) { + registry.fill(HIST("hGenPromptProduction"), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt); + } else { + registry.fill(HIST("hGenPartRecoPromptProduction"), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt); + } + } else { // non-prompt + if (!isPartRecoDstar) { + registry.fill(HIST("hGenNonPromptProduction"), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt); + } else { + registry.fill(HIST("hGenPartRecoNonPromptProduction"), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt); + } + } + } else if (cosThetaStarType == charm_polarisation::CosThetaStarType::JetAxis) { // JetAxis + if (origin == RecoDecay::OriginType::Prompt) { // prompt + if (!isPartRecoDstar) { + registry.fill(HIST("hGenPromptJetAxis"), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt); + } else { + registry.fill(HIST("hGenPartRecoPromptJetAxis"), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt); + } + } else { // non-prompt + if (!isPartRecoDstar) { + registry.fill(HIST("hGenNonPromptJetAxis"), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt); + } else { + registry.fill(HIST("hGenPartRecoNonPromptJetAxis"), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt); + } + } + } + } + + /// \param invMass is the invariant mass + /// \return true if candidate in signal region + bool isInSignalRegion(charm_polarisation::DecayChannel channel, float invMass) + { + if (channel == charm_polarisation::DecayChannel::DstarToDzeroPi) { // D*+ + float invMassSigMin = 0.142f; + float invMassSigMax = 0.15f; + if (invMassSigMin < invMass && invMass < invMassSigMax) { + return true; + } + } + + return false; + } + + /// TODO: To be implemented + /// \param daughter is the daughter particle + /// \param ptMin is the minimum pt + /// \param etaMax is the maximum eta + /// \return true if daughter is in acceptance + template + bool isDaughterInAcceptance(Part const& daughter, float ptMin, float etaMax) + { + if (daughter.pt() < ptMin) { + return false; + } + if (std::abs(daughter.eta()) > etaMax) { + return false; + } + return true; + } + + /// \param jet is the jet containing the candidates + /// \param candidates are the selected candidates + /// \param bkgRotationId is the id for the background rotation + /// \return true if candidate in signal region + template + bool runPolarisationAnalysis(charm_polarisation::DecayChannel channel, Jet const& jet, Cand const& candidate, int bkgRotationId) + { + bool isCandidateInSignalRegion{false}; + int8_t origin{RecoDecay::OriginType::None}; + float ptBhadMother{-1.f}; + bool partRecoDstar{false}; + if constexpr (doMc) { + if (channel == charm_polarisation::DecayChannel::DstarToDzeroPi) { + partRecoDstar = std::abs(candidate.flagMcMatchRec()) == hf_decay::hf_cand_dstar::DecayChannelMain::DstarToPiKPiPi0 && std::abs(candidate.flagMcMatchRecCharm()) == hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiKPi0; + bool signalDstar = std::abs(candidate.flagMcMatchRec()) == hf_decay::hf_cand_dstar::DecayChannelMain::DstarToPiKPi && std::abs(candidate.flagMcMatchRecCharm()) == hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK; + if (!signalDstar && (!partRecoDstar || !activatePartRecoDstar)) { // this candidate is not signal and not partially reconstructed signal, skip + return isCandidateInSignalRegion; + } + origin = candidate.originMcRec(); + ptBhadMother = candidate.ptBhadMotherPart(); + int pdgBhadMother = candidate.pdgBhadMotherPart(); + // For unknown reasons there are charm hadrons coming directly from beauty diquarks without an intermediate B-hadron which have an unreasonable correlation between the pT of the charm hadron and the beauty mother. We also remove charm hadrons from quarkonia. + if (origin == RecoDecay::OriginType::NonPrompt && (pdgBhadMother == 5101 || pdgBhadMother == 5103 || pdgBhadMother == 5201 || pdgBhadMother == 5203 || pdgBhadMother == 5301 || pdgBhadMother == 5303 || pdgBhadMother == 5401 || pdgBhadMother == 5403 || pdgBhadMother == 5503 || pdgBhadMother == 553 || pdgBhadMother == 555 || pdgBhadMother == 557)) { // o2-linter: disable=pdg/explicit-code, magic-number (constants not in the PDG header) + return isCandidateInSignalRegion; + } + } + } + + // loop over mass hypotheses + for (uint8_t iMass = 0u; iMass < nMassHypos; iMass++) { + + // variable definition + float pxDau{-1000.f}, pyDau{-1000.f}, pzDau{-1000.f}; + float pxCharmHad{-1000.f}, pyCharmHad{-1000.f}, pzCharmHad{-1000.f}; + float massDau{0.f}, invMassCharmHad{0.f}, invMassCharmHadForSparse{0.f}, invMassD0{0.f}; + float rapidity{-999.f}; + std::array outputMl{-1.f, -1.f, -1.f}; + std::array threeVecCand{}; + int isRotatedCandidate = 0; // currently meaningful only for Lc->pKpi + + if (channel == charm_polarisation::DecayChannel::DstarToDzeroPi) { + // Dstar analysis + // polarization measured from the soft-pion daughter (*) + + massDau = o2::constants::physics::MassPiPlus; // (*) + const float bkgRotAngle = (bkgRotationId > 0) ? minRotAngleMultByPi * constants::math::PI + bkgRotationAngleStep * (bkgRotationId - 1) : 0; + + std::array threeVecSoftPi{candidate.pxProng1() * std::cos(bkgRotAngle) - candidate.pyProng1() * std::sin(bkgRotAngle), candidate.pxProng1() * std::sin(bkgRotAngle) + candidate.pyProng1() * std::cos(bkgRotAngle), candidate.pzProng1()}; // we rotate the soft pion + std::array threeVecD0Prong0{candidate.pxProng0Charm(), candidate.pyProng0Charm(), candidate.pzProng0Charm()}; + std::array threeVecD0Prong1{candidate.pxProng1Charm(), candidate.pyProng1Charm(), candidate.pzProng1Charm()}; + if (bkgRotationId > 0) { + isRotatedCandidate = 1; + pxDau = threeVecSoftPi[0]; + pyDau = threeVecSoftPi[1]; + pzDau = threeVecSoftPi[2]; + threeVecCand = RecoDecay::pVec(threeVecSoftPi, threeVecD0Prong0, threeVecD0Prong1); + pxCharmHad = threeVecCand[0]; + pyCharmHad = threeVecCand[1]; + pzCharmHad = threeVecCand[2]; + if (candidate.signProng1() > 0) { + invMassCharmHad = RecoDecay::m(std::array{threeVecD0Prong0, threeVecD0Prong1, threeVecSoftPi}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKaonCharged, o2::constants::physics::MassPiPlus}); + invMassD0 = RecoDecay::m(std::array{threeVecD0Prong0, threeVecD0Prong1}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKaonCharged}); + } else { + invMassCharmHad = RecoDecay::m(std::array{threeVecD0Prong0, threeVecD0Prong1, threeVecSoftPi}, std::array{o2::constants::physics::MassKaonCharged, o2::constants::physics::MassPiPlus, o2::constants::physics::MassPiPlus}); + invMassD0 = RecoDecay::m(std::array{threeVecD0Prong0, threeVecD0Prong1}, std::array{o2::constants::physics::MassKaonCharged, o2::constants::physics::MassPiPlus}); + } + rapidity = RecoDecay::y(threeVecCand, o2::constants::physics::MassDStar); + } else { + isRotatedCandidate = 0; + pxDau = candidate.pxProng1(); + pyDau = candidate.pyProng1(); + pzDau = candidate.pzProng1(); + // threeVecCand = RecoDecay::pVec(std::array{candidate.pxProng1(), candidate.pyProng1(), candidate.pzProng1()}, + // std::array{candidate.pyProng0Charm(), candidate.pxProng0Charm(), candidate.pzProng0Charm()}, + // std::array{candidate.pxProng1Charm(), candidate.pyProng1Charm(), candidate.pzProng1Charm()}); + threeVecCand = {candidate.px(), candidate.py(), candidate.pz()}; + pxCharmHad = threeVecCand[0]; + pyCharmHad = threeVecCand[1]; + pzCharmHad = threeVecCand[2]; + invMassCharmHad = candidate.m(); + invMassD0 = candidate.invMassCharm(); + rapidity = candidate.y(); + } + invMassCharmHadForSparse = invMassCharmHad - invMassD0; + + if constexpr (withMl) { + std::copy_n(candidate.mlScores().begin(), outputMl.size(), outputMl.begin()); + } + } + if (invMassCharmHadForSparse < invMassMin || invMassCharmHadForSparse > invMassMax) { + continue; + } + + /// apply rapidity selection on the reconstructed candidate + if (std::abs(rapidity) > absRapidityCandMax) { + continue; + } + + ROOT::Math::PxPyPzMVector fourVecDau = ROOT::Math::PxPyPzMVector(pxDau, pyDau, pzDau, massDau); + ROOT::Math::PxPyPzMVector fourVecMother = ROOT::Math::PxPyPzMVector(pxCharmHad, pyCharmHad, pzCharmHad, invMassCharmHad); + ROOT::Math::Boost boost{fourVecMother.BoostToCM()}; + ROOT::Math::PxPyPzMVector fourVecDauCM = boost(fourVecDau); + ROOT::Math::XYZVector threeVecDauCM = fourVecDauCM.Vect(); + ROOT::Math::PxPyPzMVector fourVecJet = ROOT::Math::PxPyPzMVector(jet.px(), jet.py(), jet.pz(), jet.mass()); + ROOT::Math::PxPyPzMVector fourVecJetCM = boost(fourVecJet); + + float ptCharmHad = RecoDecay::pt(threeVecCand); // this definition is valid for both rotated and original candidates + + if (!isCandidateInSignalRegion) { // it could be that only one mass hypothesis is in signal region + isCandidateInSignalRegion = isInSignalRegion(channel, invMassCharmHadForSparse); + } + float absEtaTrackMin = candidate.absEtaTrackMin(); + int numItsClsMin = candidate.numItsClsMin(); + int numTpcClsMin = candidate.numTpcClsMin(); + + // helicity + ROOT::Math::XYZVector helicityVec = fourVecMother.Vect(); + float cosThetaStarHelicity = -10.f; + // production + ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector(pyCharmHad, -pxCharmHad, 0.f); + float cosThetaStarProduction = -10.f; + // jet axis + ROOT::Math::XYZVector jetaxisVec = fourVecJetCM.Vect(); + float cosThetaStarJet = -10.f; + + float zParallel = helicityVec.Dot(jetaxisVec) / std::sqrt(jetaxisVec.Mag2()); + float jetPt = jet.pt(); + + int8_t nMuons{0u}; + if constexpr (doMc) { + nMuons = candidate.nTracksDecayed(); + } + if (activateTHnSparseCosThStarHelicity) { + // helicity + cosThetaStarHelicity = helicityVec.Dot(threeVecDauCM) / std::sqrt(threeVecDauCM.Mag2()) / std::sqrt(helicityVec.Mag2()); + HfHistoInput helicityInput{ + .invMassCharmHad = invMassCharmHadForSparse, + .ptCharmHad = ptCharmHad, + .rapCharmHad = rapidity, + .invMassD0 = invMassD0, + .cosThetaStar = cosThetaStarHelicity, + .outputMl = outputMl, + .isRotatedCandidate = isRotatedCandidate, + .origin = origin, + .ptBhadMother = ptBhadMother, + .absEtaMin = absEtaTrackMin, + .numItsClsMin = numItsClsMin, + .numTpcClsMin = numTpcClsMin, + .nMuons = nMuons, + .zParallel = zParallel, + .jetPt = jetPt}; + fillRecoHistos(charm_polarisation::CosThetaStarType::Helicity, withMl, doMc, partRecoDstar, helicityInput); + } + if (activateTHnSparseCosThStarProduction) { + // production + cosThetaStarProduction = normalVec.Dot(threeVecDauCM) / std::sqrt(threeVecDauCM.Mag2()) / std::sqrt(normalVec.Mag2()); + HfHistoInput productionInput{ + .invMassCharmHad = invMassCharmHadForSparse, + .ptCharmHad = ptCharmHad, + .rapCharmHad = rapidity, + .invMassD0 = invMassD0, + .cosThetaStar = cosThetaStarProduction, + .outputMl = outputMl, + .isRotatedCandidate = isRotatedCandidate, + .origin = origin, + .ptBhadMother = ptBhadMother, + .absEtaMin = absEtaTrackMin, + .numItsClsMin = numItsClsMin, + .numTpcClsMin = numTpcClsMin, + .nMuons = nMuons, + .zParallel = zParallel, + .jetPt = jetPt}; + fillRecoHistos(charm_polarisation::CosThetaStarType::Production, withMl, doMc, partRecoDstar, productionInput); + } + if (activateTHnSparseCosThStarJetAxis) { + // jet axis + cosThetaStarJet = jetaxisVec.Dot(threeVecDauCM) / std::sqrt(threeVecDauCM.Mag2()) / std::sqrt(jetaxisVec.Mag2()); + HfHistoInput jetAxisInput{ + .invMassCharmHad = invMassCharmHadForSparse, + .ptCharmHad = ptCharmHad, + .rapCharmHad = rapidity, + .invMassD0 = invMassD0, + .cosThetaStar = cosThetaStarJet, + .outputMl = outputMl, + .isRotatedCandidate = isRotatedCandidate, + .origin = origin, + .ptBhadMother = ptBhadMother, + .absEtaMin = absEtaTrackMin, + .numItsClsMin = numItsClsMin, + .numTpcClsMin = numTpcClsMin, + .nMuons = nMuons, + .zParallel = zParallel, + .jetPt = jetPt}; + fillRecoHistos(charm_polarisation::CosThetaStarType::JetAxis, withMl, doMc, partRecoDstar, jetAxisInput); + } + } /// end loop over mass hypotheses + + return isCandidateInSignalRegion; + } + + ///////////////////////// + // Dstar analysis /// + ///////////////////////// + + // Dstar with rectangular cuts + void processDstar(aod::JetCollisions const& collisions, + DstarJets const& jets, + aod::CandidatesDstarData const&) + { + for (const auto& collision : collisions) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + continue; + } + auto thisCollId = collision.globalIndex(); + auto groupedDstarjets = jets.sliceBy(dstarJetsPerCollision, thisCollId); + for (const auto& jet : groupedDstarjets) { + for (const auto& dstarCandidate : jet.candidates_as()) { + runPolarisationAnalysis(charm_polarisation::DecayChannel::DstarToDzeroPi, jet, dstarCandidate, 0); + for (int iRotation{1}; iRotation <= nBkgRotations; ++iRotation) { + runPolarisationAnalysis(charm_polarisation::DecayChannel::DstarToDzeroPi, jet, dstarCandidate, iRotation); + } + break; // hf jet should have only one Dstar candidate but for safety + } + } + } + } + PROCESS_SWITCH(HfTaskDstarPolarisationInJet, processDstar, "Process Dstar candidates without ML", true); + + // Dstar with ML cuts + void processDstarWithMl(aod::JetCollisions const& collisions, + DstarJets const& jets, + aod::CandidatesDstarData const&) + { + for (const auto& collision : collisions) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + continue; + } + auto thisCollId = collision.globalIndex(); + auto groupedDstarJets = jets.sliceBy(dstarJetsPerCollision, thisCollId); + for (const auto& jet : groupedDstarJets) { + for (const auto& dstarCandidate : jet.candidates_as()) { + runPolarisationAnalysis(charm_polarisation::DecayChannel::DstarToDzeroPi, jet, dstarCandidate, 0); + for (int iRotation{1}; iRotation <= nBkgRotations; ++iRotation) { + runPolarisationAnalysis(charm_polarisation::DecayChannel::DstarToDzeroPi, jet, dstarCandidate, iRotation); + } + break; // hf jet should have only one Dstar candidate but for safety + } + } + } + } + PROCESS_SWITCH(HfTaskDstarPolarisationInJet, processDstarWithMl, "Process Dstar candidates with ML", false); + + // Dstar in MC with rectangular cuts + void processDstarMc(aod::JetMcCollisions const& mcCollisions, + aod::JetCollisions const& collisions, + JetMCDTable const& mcdJets, + aod::CandidatesDstarMCD const&) + { + for (const auto& mcCollision : mcCollisions) { + const auto collisionsPerMCCollision = collisions.sliceBy(collisionsPerMCCollisionPreslice, mcCollision.globalIndex()); + for (const auto& collision : collisionsPerMCCollision) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + continue; + } + const auto dstarMcDJetsPerCollision = mcdJets.sliceBy(dstarMCDJetsPerCollisionPreslice, collision.globalIndex()); + for (const auto& mcdJet : dstarMcDJetsPerCollision) { + for (const auto& mcdDstarCand : mcdJet.candidates_as()) { + runPolarisationAnalysis(charm_polarisation::DecayChannel::DstarToDzeroPi, mcdJet, mcdDstarCand, 0); + break; + } + } + } + } + } + PROCESS_SWITCH(HfTaskDstarPolarisationInJet, processDstarMc, "Process Dstar candidates in MC without ML", false); + + // Dstar in MC with ML cuts + void processDstarMcWithMl(aod::JetMcCollisions const& mcCollisions, + aod::JetCollisions const& collisions, + JetMCDTable const& mcdJets, + aod::CandidatesDstarMCD const&) + { + for (const auto& mcCollision : mcCollisions) { + const auto collisionsPerMCCollision = collisions.sliceBy(collisionsPerMCCollisionPreslice, mcCollision.globalIndex()); + for (const auto& collision : collisionsPerMCCollision) { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + continue; + } + const auto dstarMcdJetsPerCollision = mcdJets.sliceBy(dstarMCDJetsPerCollisionPreslice, collision.globalIndex()); + for (const auto& mcdJet : dstarMcdJetsPerCollision) { + for (const auto& mcdDstarCand : mcdJet.candidates_as()) { + runPolarisationAnalysis(charm_polarisation::DecayChannel::DstarToDzeroPi, mcdJet, mcdDstarCand, 0); + } + } + } + } + } + PROCESS_SWITCH(HfTaskDstarPolarisationInJet, processDstarMcWithMl, "Process Dstar candidates in MC with ML", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/TableProducer/derivedDataCreatorDstarToD0Pi.cxx b/PWGHF/TableProducer/derivedDataCreatorDstarToD0Pi.cxx index dfbfce30e68..1587c99fb71 100644 --- a/PWGHF/TableProducer/derivedDataCreatorDstarToD0Pi.cxx +++ b/PWGHF/TableProducer/derivedDataCreatorDstarToD0Pi.cxx @@ -88,7 +88,7 @@ struct HfDerivedDataCreatorDstarToD0Pi { using CollisionsWCentMult = soa::Join; using CollisionsWMcCentMult = soa::Join; - using TracksWPid = soa::Join; + using TracksWPid = soa::Join; using SelectedCandidates = soa::Filtered>; using SelectedCandidatesMc = soa::Filtered>; using SelectedCandidatesMl = soa::Filtered>; @@ -125,12 +125,34 @@ struct HfDerivedDataCreatorDstarToD0Pi { rowsCommon.init(confDerData); } + /// prongTracks is the vector of daughter tracks + /// etaMin is the minimum eta + /// nItsClsMin is the minumum number of clusters in ITS + /// nTpcClsMin is the minumum number of clusters in TPC + template + void getTrackingInfos(std::array const& prongTracks, float& etaMin, uint8_t& nItsClsMin, int16_t& nTpcClsMin) + { + etaMin = 10.f; + nItsClsMin = 100; + nTpcClsMin = 1000; + + for (const auto& track : prongTracks) { + etaMin = std::min(etaMin, std::abs(track.eta())); + nItsClsMin = std::min(nItsClsMin, track.itsNCls()); + nTpcClsMin = std::min(nTpcClsMin, track.tpcNClsCrossedRows()); + } + } + template void fillTablesCandidate(const T& candidate, const U& prong0, const U& prong1, const U& prongSoftPi, int candFlag, double invMass, double invMassD0, double y, int8_t flagMc, int8_t flagMcD0, int8_t origin, int8_t nTracksDecayed, double ptBhad, int pdgBhad, const std::vector& mlScores) { rowsCommon.fillTablesCandidate(candidate, invMass, y); if (fillCandidatePar) { + float absEtaTrackMin{-1.f}; + uint8_t numItsClsMin{200u}; + int16_t numTpcClsMin{1000}; + getTrackingInfos(std::array{prong0, prong1, prongSoftPi}, absEtaTrackMin, numItsClsMin, numTpcClsMin); rowCandidatePar( candidate.pxD0(), candidate.pyD0(), @@ -143,7 +165,10 @@ struct HfDerivedDataCreatorDstarToD0Pi { candidate.normalisedImpParamSoftPi(), prongSoftPi.tpcNSigmaPi(), prongSoftPi.tofNSigmaPi(), - prongSoftPi.tpcTofNSigmaPi()); + prongSoftPi.tpcTofNSigmaPi(), + absEtaTrackMin, + numItsClsMin, + numTpcClsMin); } if (fillCandidateParD0) { rowCandidateParD0( @@ -250,11 +275,6 @@ struct HfDerivedDataCreatorDstarToD0Pi { double ptBhadMotherPart = 0; int pdgBhadMotherPart = 0; for (const auto& candidate : candidatesThisColl) { - if constexpr (IsMl) { - if (!TESTBIT(candidate.isSelDstarToD0Pi(), aod::SelectionStep::RecoMl)) { - continue; - } - } if constexpr (IsMc) { flagMcRec = candidate.flagMcMatchRec(); flagMcRecD0 = candidate.flagMcMatchRecD0(); diff --git a/PWGJE/Core/JetHFUtilities.h b/PWGJE/Core/JetHFUtilities.h index 6e25d808d37..afefa3ea0d9 100644 --- a/PWGJE/Core/JetHFUtilities.h +++ b/PWGJE/Core/JetHFUtilities.h @@ -1065,7 +1065,10 @@ void fillDstarCandidateTable(T const& candidate, U& DstarParTable, V& DstarParDa candidate.impactParameterNormalised1(), candidate.nSigTpcPi1(), candidate.nSigTofPi1(), - candidate.nSigTpcTofPi1()); + candidate.nSigTpcTofPi1(), + candidate.absEtaTrackMin(), + candidate.numItsClsMin(), + candidate.numTpcClsMin()); DstarParDaughterTable( candidate.chi2PCACharm(),