diff --git a/PWGCF/DataModel/FemtoDerived.h b/PWGCF/DataModel/FemtoDerived.h index cac62250bbe..cf93d9f9b4f 100644 --- a/PWGCF/DataModel/FemtoDerived.h +++ b/PWGCF/DataModel/FemtoDerived.h @@ -411,6 +411,55 @@ DECLARE_SOA_DYNAMIC_COLUMN(Y, y, //! o2::constants::physics::MassXiCPlus); }); //! Rapidity distribution of Xic candidates } // namespace fdhf_xic +namespace fdhf_sigmac +{ +DECLARE_SOA_COLUMN(SoftPiId, softPiId, int); //! original soft-pion track ID +DECLARE_SOA_COLUMN(SoftPiPt, softPiPt, float); //! soft-pion pT +DECLARE_SOA_COLUMN(SoftPiEta, softPiEta, float); //! soft-pion eta +DECLARE_SOA_COLUMN(SoftPiPhi, softPiPhi, float); //! soft-pion phi +DECLARE_SOA_COLUMN(ChargeLc, chargeLc, int8_t); //! LambdaC charge, also particle/antiparticle tag for neutral SigmaC +DECLARE_SOA_COLUMN(MassLc, massLc, float); //! LambdaC mass for the accepted hypothesis +DECLARE_SOA_COLUMN(DeltaMass, deltaMass, float); //! M(Lc_reco pi) - M(Lc_reco) +DECLARE_SOA_COLUMN(SourceScId, sourceScId, int64_t); //! source HF SigmaC candidate, shared by its hypotheses +DECLARE_SOA_COLUMN(SoftPiDcaXY, softPiDcaXY, float); //! DCA to the SigmaC collision, including reassociation +DECLARE_SOA_COLUMN(SoftPiDcaZ, softPiDcaZ, float); //! DCA to the SigmaC collision, including reassociation +} // namespace fdhf_sigmac + +DECLARE_SOA_TABLE(FDHfCandSc, "AOD", "FDHFCANDSC", //! SigmaC(2455) -> Lc(p K pi) pi, one row per accepted Lc hypothesis + o2::soa::Index<>, + femtodreamparticle::FDCollisionId, + fdhf::TimeStamp, + fdhf::Charge, + fdhf::Prong0Id, + fdhf::Prong1Id, + fdhf::Prong2Id, + fdhf::Prong0Pt, + fdhf::Prong1Pt, + fdhf::Prong2Pt, + fdhf::Prong0Eta, + fdhf::Prong1Eta, + fdhf::Prong2Eta, + fdhf::Prong0Phi, + fdhf::Prong1Phi, + fdhf::Prong2Phi, + fdhf::CandidateSelFlag, + fdhf::BDTBkg, + fdhf::BDTPrompt, + fdhf::BDTFD, + femtodreamparticle::Pt, + femtodreamparticle::Eta, + femtodreamparticle::Phi, + fdhf_sigmac::SoftPiId, + fdhf_sigmac::SoftPiPt, + fdhf_sigmac::SoftPiEta, + fdhf_sigmac::SoftPiPhi, + fdhf_sigmac::ChargeLc, + fdhf_sigmac::MassLc, + fdhf_sigmac::DeltaMass, + fdhf_sigmac::SourceScId, + fdhf_sigmac::SoftPiDcaXY, + fdhf_sigmac::SoftPiDcaZ); + DECLARE_SOA_TABLE(FDHfCand3Prong, "AOD", "FDHFCAND3PRONG", //! Table to store the derived data for charm 3prong candidates o2::soa::Index<>, femtodreamparticle::FDCollisionId, @@ -605,6 +654,28 @@ DECLARE_SOA_TABLE(FDHfCharmDstar, "AOD", "FDHFCHARMDSTAR", //! table to store re fdhf::FlagMc, fdhf::OriginMcRec); +DECLARE_SOA_TABLE(FDHfCharmSc, "AOD", "FDHFCHARMSC", //! table to store results for HF femtoscopy + fdhf::GIndexCol, + fdhf::TimeStamp, + fdhf::CharmM, + fdhf::CharmDaughM, + fdhf::CharmPt, + fdhf::CharmEta, + fdhf::CharmPhi, + fdhf::Prong0Id, + fdhf::Prong1Id, + fdhf::Prong2Id, + fdhf_sigmac::SoftPiId, + fdhf_sigmac::ChargeLc, + fdhf_sigmac::SourceScId, + fdhf::CandidateSelFlag, + fdhf::Charge, + fdhf::BDTBkg, + fdhf::BDTPrompt, + fdhf::BDTFD, + fdhf::FlagMc, + fdhf::OriginMcRec); + DECLARE_SOA_TABLE(FDHfTrk, "AOD", "FDHFTRK", //! table to store results for HF femtoscopy fdhf::GIndexCol, fdhf::TimeStamp, diff --git a/PWGCF/FemtoDream/Core/femtoDreamUtils.h b/PWGCF/FemtoDream/Core/femtoDreamUtils.h index fde124f5db4..df48c5a57dc 100644 --- a/PWGCF/FemtoDream/Core/femtoDreamUtils.h +++ b/PWGCF/FemtoDream/Core/femtoDreamUtils.h @@ -71,6 +71,12 @@ inline float getMass(int pdgCode) case o2::constants::physics::Pdg::kDS: mass = o2::constants::physics::MassDS; break; + case o2::constants::physics::Pdg::kSigmaC0: + mass = o2::constants::physics::MassSigmaC0; + break; + case o2::constants::physics::Pdg::kSigmaCPlusPlus: + mass = o2::constants::physics::MassSigmaCPlusPlus; + break; case o2::constants::physics::Pdg::kLambdaCPlus: mass = o2::constants::physics::MassLambdaCPlus; break; diff --git a/PWGHF/HFC/TableProducer/producerCharmHadronsTrackFemtoDream.cxx b/PWGHF/HFC/TableProducer/producerCharmHadronsTrackFemtoDream.cxx index 0b1102163c7..7aa75dc27a5 100644 --- a/PWGHF/HFC/TableProducer/producerCharmHadronsTrackFemtoDream.cxx +++ b/PWGHF/HFC/TableProducer/producerCharmHadronsTrackFemtoDream.cxx @@ -22,12 +22,14 @@ #include "PWGHF/Core/CentralityEstimation.h" #include "PWGHF/Core/DecayChannels.h" #include "PWGHF/Core/DecayChannelsLegacy.h" +#include "PWGHF/Core/HfHelper.h" #include "PWGHF/Core/HfMlResponseD0ToKPi.h" #include "PWGHF/Core/HfMlResponseDplusToPiKPi.h" #include "PWGHF/Core/HfMlResponseDstarToD0Pi.h" #include "PWGHF/Core/HfMlResponseLcToPKPi.h" #include "PWGHF/Core/HfMlResponseXicToXiPiPi.h" #include "PWGHF/Core/SelectorCuts.h" +#include "PWGHF/D2H/Utils/utilsSigmac.h" #include "PWGHF/DataModel/AliasTables.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" @@ -67,6 +69,7 @@ #include #include +#include #include #include #include @@ -102,7 +105,8 @@ enum DecayChannel { DplusToPiKPi = 0, LcToPKPi, D0ToPiK, DstarToD0Pi, - XicToXiPiPi + XicToXiPiPi, + ScToLcPi }; enum class D0CandFlag : uint8_t { @@ -121,6 +125,7 @@ struct HfProducerCharmHadronsTrackFemtoDream { Produces rowCandCharm3ProngXicQa; Produces rowCandCharm2Prong; Produces rowCandCharmDstar; + Produces rowCandCharmSc; Produces rowCandMcCharmHad; Produces rowCandCharmHadGen; Produces outputPartsIndex; @@ -200,6 +205,15 @@ struct HfProducerCharmHadronsTrackFemtoDream { Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; } mlCfg; + struct : ConfigurableGroup { + std::string prefix = "sigmac"; + Configurable deltaMassMin{"deltaMassMin", static_cast(o2::constants::physics::MassPiPlus), "Minimum retained DeltaM, including sidebands (GeV/c2)"}; + Configurable deltaMassMax{"deltaMassMax", 0.190f, "Maximum retained DeltaM, including sidebands (GeV/c2)"}; + Configurable lcMassWindow{"lcMassWindow", 0.03f, "Maximum |M(Lc)-PDG|, also limited by the HF creator (GeV/c2)"}; + Configurable yMax{"yMax", 0.8f, "Maximum |y| using the SigmaC(2455) mass"}; + Configurable rejectAmbiguousHypothesis{"rejectAmbiguousHypothesis", false, "Reject candidates with both accepted Lc hypotheses"}; + } sigmacCfg; + FemtoDreamTrackSelection trackCuts; o2::analysis::HfMlResponseLcToPKPi hfMlResponseLc; @@ -261,7 +275,8 @@ struct HfProducerCharmHadronsTrackFemtoDream { void init(InitContext&) { - std::array processes = {doprocessDataDplusToPiKPi, doprocessMcDplusToPiKPi, doprocessDataDplusToPiKPiWithML, doprocessMcDplusToPiKPiWithML, doprocessMcDplusToPiKPiGen, + std::array processes = {doprocessDataScToLcPi, doprocessDataScToLcPiWithML, doprocessMcScToLcPi, doprocessMcScToLcPiWithML, doprocessMcScToLcPiGen, + doprocessDataDplusToPiKPi, doprocessMcDplusToPiKPi, doprocessDataDplusToPiKPiWithML, doprocessMcDplusToPiKPiWithML, doprocessMcDplusToPiKPiGen, doprocessDataLcToPKPi, doprocessMcLcToPKPi, doprocessDataLcToPKPiWithML, doprocessMcLcToPKPiWithML, doprocessMcLcToPKPiGen, doprocessDataD0ToPiK, doprocessMcD0ToPiK, doprocessDataD0ToPiKWithML, doprocessMcD0ToPiKWithML, doprocessMcD0ToPiKGen, doprocessDataDstarToD0Pi, doprocessMcDstarToD0Pi, doprocessDataDstarToD0PiWithML, doprocessMcDstarToD0PiWithML, doprocessMcDstarToD0PiGen, @@ -271,6 +286,16 @@ struct HfProducerCharmHadronsTrackFemtoDream { LOGP(fatal, "One and only one process function must be enabled at a time."); } + if ((doprocessDataScToLcPi || doprocessDataScToLcPiWithML || doprocessMcScToLcPi || doprocessMcScToLcPiWithML) && + (sigmacCfg.deltaMassMin.value >= sigmacCfg.deltaMassMax.value || sigmacCfg.lcMassWindow.value <= 0.f || sigmacCfg.yMax.value <= 0.f)) { + LOG(fatal) << "Invalid SigmaC mass or rapidity window."; + } + if (doprocessDataScToLcPi || doprocessDataScToLcPiWithML || doprocessMcScToLcPi || doprocessMcScToLcPiWithML) { + qaRegistry.add("Sigmac/hPtVsDeltaMass", ";pT(SigmaC) (GeV/c);DeltaM (GeV/c2)", kTH2F, {{180, 0., 36.}, {300, 0.138, 0.198}}); + qaRegistry.add("Sigmac/hLcMass", ";M(Lc) (GeV/c2);Entries", kTH1F, {{300, 2.15, 2.42}}); + qaRegistry.add("Sigmac/hSoftPiPt", ";pT(soft pi) (GeV/c);Entries", kTH1F, {{200, 0., 2.}}); + } + int const cutBits = 8 * sizeof(o2::aod::femtodreamparticle::cutContainerType); trackRegistry.add("AnalysisQA/CutCounter", "; Bit; Counter", kTH1F, {{cutBits + 1, -0.5, cutBits + 0.5}}); @@ -319,7 +344,7 @@ struct HfProducerCharmHadronsTrackFemtoDream { 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 useLcMl = doprocessDataLcToPKPiWithML || doprocessMcLcToPKPiWithML || doprocessDataScToLcPiWithML || doprocessMcScToLcPiWithML; bool useDplusMl = doprocessDataDplusToPiKPiWithML || doprocessMcDplusToPiKPiWithML; bool useD0Ml = doprocessDataD0ToPiKWithML || doprocessMcD0ToPiKWithML; bool useDstarMl = doprocessDataDstarToD0PiWithML || doprocessMcDstarToD0PiWithML; @@ -1084,10 +1109,10 @@ struct HfProducerCharmHadronsTrackFemtoDream { /// \param TC type of the femto track cuts /// \return whether or not the tracks fulfills the all selections template - bool isNoSelectedTracks(C const& /*col*/, T const& tracks, TC& trackCuts) + bool isNoSelectedTracks(C const& /*col*/, T const& tracks, TC& trackSelection) { for (auto const& track : tracks) { - if (trackCuts.isSelectedMinimal(track)) { + if (trackSelection.isSelectedMinimal(track)) { return false; } } @@ -1378,6 +1403,210 @@ struct HfProducerCharmHadronsTrackFemtoDream { } PROCESS_SWITCH(HfProducerCharmHadronsTrackFemtoDream, processMcLcToPKPiGen, "Provide Mc Generated lctopkpi", false); + template + void fillSigmacHadronTable(CollisionType const& col, TrackType const& tracks, ScType const& candidates, + LcType const& /*lcs*/, DaughterTracks const& /*daughterTracks*/, aod::McParticles const* mcParticles = nullptr) + { + 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 (generalCfg.isRun3) { + if (generalCfg.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 (isNoSelectedTracks(col, tracks, trackCuts) && sizeCand <= 0) { + qaRegistry.fill(HIST("hEventQA"), 1 + Event::RejNoTracksAndCharm); + return; + } + + tables.outputCollision(vtxZ, mult, multNtr, spher, magField); + if constexpr (IsMc) { + fillMcCollision(col); + } + + int nSelectedSc = 0; + constexpr std::size_t NumMlScores = 3; + constexpr int NumHypotheses = 2; + constexpr int ChargeDoublyCharged = 2; + auto bc = col.template bc_as(); + for (auto const& sc : candidates) { + if (sc.charge() != 0 && std::abs(sc.charge()) != ChargeDoublyCharged) { + continue; + } + auto const rapidity = sc.charge() == 0 ? HfHelper::ySc0(sc) : HfHelper::yScPlusPlus(sc); + if (!std::isfinite(rapidity) || std::abs(rapidity) > sigmacCfg.yMax.value) { + continue; + } + auto lc = sc.template prongLc_as(); + if (lc.collisionId() != sc.collisionId()) { + continue; + } + auto p0 = lc.template prong0_as(); + auto p1 = lc.template prong1_as(); + auto p2 = lc.template prong2_as(); + auto softPi = sc.template prong1_as(); + if (softPi.globalIndex() == p0.globalIndex() || softPi.globalIndex() == p1.globalIndex() || softPi.globalIndex() == p2.globalIndex()) { + continue; + } + std::array const lcMass{HfHelper::invMassLcToPKPi(lc), HfHelper::invMassLcToPiKP(lc)}; + std::array const scMass{HfHelper::invMassScRecoLcToPKPi(sc, lc), HfHelper::invMassScRecoLcToPiKP(sc, lc)}; + auto const decayMask = o2::hf_sigmac_utils::isDecayToPKPiToPiKP(lc, sc); + std::array accepted{}; + std::array, NumHypotheses> scores{std::vector(NumMlScores, -1.f), std::vector(NumMlScores, -1.f)}; + for (int h = 0; h < NumHypotheses; ++h) { + int const selected = h == 0 ? lc.isSelLcToPKPi() : lc.isSelLcToPiKP(); + accepted[h] = (decayMask & (1 << h)) && selected >= generalCfg.selectionFlagHadron.value && + std::isfinite(lcMass[h]) && std::abs(lcMass[h] - o2::constants::physics::MassLambdaCPlus) <= sigmacCfg.lcMassWindow.value; + if (!accepted[h]) { + continue; + } + if constexpr (UseCharmMl) { + if (mlCfg.applyMlMode == FillMlFromSelector) { + auto const& probabilities = h == 0 ? lc.mlProbLcToPKPi() : lc.mlProbLcToPiKP(); + accepted[h] = probabilities.size() >= NumMlScores; + if (accepted[h]) { + scores[h].assign(probabilities.begin(), probabilities.begin() + NumMlScores); + } + } else if (mlCfg.applyMlMode == FillMlFromNewBDT) { + auto features = hfMlResponseLc.getInputFeatures(lc, h == 0); + accepted[h] = hfMlResponseLc.isSelectedMl(features, lc.pt(), scores[h]); + accepted[h] = accepted[h] && scores[h].size() >= NumMlScores; + } else if (mlCfg.applyMlMode != NoMl) { + LOG(fatal) << "Invalid ML mode for SigmaC LambdaC daughter."; + } + } + } + if (sigmacCfg.rejectAmbiguousHypothesis.value && accepted[0] && accepted[1]) { + continue; + } + for (int h = 0; h < NumHypotheses; ++h) { + double const deltaMass = scMass[h] - lcMass[h]; + if (!accepted[h] || !std::isfinite(deltaMass) || deltaMass < sigmacCfg.deltaMassMin.value || deltaMass > sigmacCfg.deltaMassMax.value) { + continue; + } + tables.rowCandCharmSc(tables.outputCollision.lastIndex(), bc.timestamp(), sc.charge(), + p0.globalIndex(), p1.globalIndex(), p2.globalIndex(), p0.pt(), p1.pt(), p2.pt(), + p0.eta(), p1.eta(), p2.eta(), p0.phi(), p1.phi(), p2.phi(), + 1 << h, scores[h][0], scores[h][1], scores[h][2], + sc.pt(), sc.eta(), sc.phi(), softPi.globalIndex(), softPi.pt(), softPi.eta(), softPi.phi(), + p0.sign(), lcMass[h], deltaMass, sc.globalIndex(), sc.softPiDcaXY(), sc.softPiDcaZ()); + if constexpr (IsMc) { + // HF matches the complete SigmaC chain. Invalidate the wrong p/pi + // assignment separately so a reflection cannot masquerade as signal. + int flag = sc.flagMcMatchRec(); + int origin = sc.originMcRec(); + int const expectedPdg0 = h == 0 ? kProton : kPiPlus; + if (!p0.has_mcParticle() || std::abs(p0.mcParticle().pdgCode()) != expectedPdg0) { + flag = 0; + origin = 0; + } + tables.rowCandMcCharmHad(flag, origin); + } + ++nSelectedSc; + qaRegistry.fill(HIST("Sigmac/hPtVsDeltaMass"), sc.pt(), deltaMass); + qaRegistry.fill(HIST("Sigmac/hLcMass"), lcMass[h]); + qaRegistry.fill(HIST("Sigmac/hSoftPiPt"), softPi.pt()); + } + } + bool const isTrackFilled = fillTracksForCharmHadron(col, tracks, mcParticles); + aod::femtodreamcollision::BitMaskType const bitTrack = isTrackFilled ? 1 : 0; + aod::femtodreamcollision::BitMaskType const bitCand = nSelectedSc > 0 ? 1 : 0; + if (isTrackFilled) { + qaRegistry.fill(HIST("hEventQA"), 1 + Event::TrackSelected); + } + if (nSelectedSc > 0) { + qaRegistry.fill(HIST("hEventQA"), 1 + Event::CharmSelected); + } + if (isTrackFilled && nSelectedSc > 0) { + qaRegistry.fill(HIST("hEventQA"), 1 + Event::PairSelected); + } + tables.rowMasks(bitTrack, bitCand, 0); + } + + void processDataScToLcPi(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + aod::Tracks const& daughterTracks, + CandidateLc const& lcs, + aod::HfCandSc const& candidates) + { + getMagneticFieldTesla(col.bc_as()); + fillSigmacHadronTable(col, tracks, candidates, lcs, daughterTracks); + } + PROCESS_SWITCH(HfProducerCharmHadronsTrackFemtoDream, processDataScToLcPi, "SigmaC(2455) -> Lc pi femto data", false); + + void processDataScToLcPiWithML(FemtoFullCollision const& col, + aod::BCsWithTimestamps const&, + FemtoHFTracks const& tracks, + aod::Tracks const& daughterTracks, + soa::Join const& lcs, + aod::HfCandSc const& candidates) + { + getMagneticFieldTesla(col.bc_as()); + fillSigmacHadronTable(col, tracks, candidates, lcs, daughterTracks); + } + PROCESS_SWITCH(HfProducerCharmHadronsTrackFemtoDream, processDataScToLcPiWithML, "SigmaC(2455) -> Lc pi femto data with Lc ML", false); + + void processMcScToLcPi(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + aod::TracksWMc const& daughterTracks, + CandidateLcMc const& lcs, + soa::Join const& candidates, + aod::McParticles const& mcParticles) + { + getMagneticFieldTesla(col.bc_as()); + fillSigmacHadronTable(col, tracks, candidates, lcs, daughterTracks, &mcParticles); + } + PROCESS_SWITCH(HfProducerCharmHadronsTrackFemtoDream, processMcScToLcPi, "SigmaC(2455) -> Lc pi femto MC", false); + + void processMcScToLcPiWithML(FemtoFullCollisionMc const& col, + aod::BCsWithTimestamps const&, + FemtoHFMcTracks const& tracks, + aod::TracksWMc const& daughterTracks, + soa::Join const& lcs, + soa::Join const& candidates, + aod::McParticles const& mcParticles) + { + getMagneticFieldTesla(col.bc_as()); + fillSigmacHadronTable(col, tracks, candidates, lcs, daughterTracks, &mcParticles); + } + PROCESS_SWITCH(HfProducerCharmHadronsTrackFemtoDream, processMcScToLcPiWithML, "SigmaC(2455) -> Lc pi femto MC with Lc ML", false); + + void processMcScToLcPiGen(soa::Join const& particles) + { + for (auto const& particle : particles) { + int const flag = std::abs(particle.flagMcMatchGen()); + if (flag == o2::hf_decay::hf_cand_sigmac::DecayChannelMain::Sc0ToPKPiPi || + flag == o2::hf_decay::hf_cand_sigmac::DecayChannelMain::ScplusplusToPKPiPi) { + tables.rowCandCharmHadGen(particle.mcCollisionId(), particle.flagMcMatchGen(), particle.originMcGen()); + } + } + } + PROCESS_SWITCH(HfProducerCharmHadronsTrackFemtoDream, processMcScToLcPiGen, "Generated SigmaC(2455) -> Lc pi", false); + /// XicToXiPiPi void processDataXicToXiPiPi(FemtoFullCollision const& col, aod::BCsWithTimestamps const&, diff --git a/PWGHF/HFC/Tasks/taskCharmHadronsTrackFemtoDream.cxx b/PWGHF/HFC/Tasks/taskCharmHadronsTrackFemtoDream.cxx index cfbd216ec22..000e0479ffb 100644 --- a/PWGHF/HFC/Tasks/taskCharmHadronsTrackFemtoDream.cxx +++ b/PWGHF/HFC/Tasks/taskCharmHadronsTrackFemtoDream.cxx @@ -46,6 +46,7 @@ #include #include +#include #include #include #include @@ -87,20 +88,25 @@ struct HfTaskCharmHadronsTrackFemtoDream { LcToPKPi, D0ToPiK, DstarToD0Pi, - XicToXiPiPi + XicToXiPiPi, + ScToLcPi }; constexpr static int OriginRecPrompt = 1; constexpr static int OriginRecFD = 2; constexpr static int CutBitChargePositive = 2; constexpr static uint32_t PidTrackPass = 1u; - - Produces rowFemtoResultPairs; - Produces rowFemtoResultCharm3Prong; - Produces rowFemtoResultCharm2Prong; - Produces rowFemtoResultCharmDstar; - Produces rowFemtoResultTrk; - Produces rowFemtoResultColl; + static constexpr int ChargeDoublyCharged = 2; + + struct : ProducesGroup { + Produces rowFemtoResultPairs; + Produces rowFemtoResultCharm3Prong; + Produces rowFemtoResultCharm2Prong; + Produces rowFemtoResultCharmDstar; + Produces rowFemtoResultCharmSc; + Produces rowFemtoResultTrk; + Produces rowFemtoResultColl; + } tables; Configurable confTempFitVarMomentum{"confTempFitVarMomentum", 0, "Momentum used for binning: 0 -> pt; 1 -> preco; 2 -> ptpc"}; @@ -119,6 +125,14 @@ struct HfTaskCharmHadronsTrackFemtoDream { 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; + struct : ConfigurableGroup { + std::string prefix = "sigmac"; + Configurable particleSign{"particleSign", 0, "1: particle, -1: antiparticle, 0: both (also for neutral SigmaC)"}; + Configurable deltaMassMin{"deltaMassMin", 0.160f, "Lower DeltaM edge (GeV/c2), tune using mass fits"}; + Configurable deltaMassMax{"deltaMassMax", 0.175f, "Upper DeltaM edge (GeV/c2), configurable for sidebands"}; + Configurable lcMassWindow{"lcMassWindow", 0.03f, "Maximum |M(Lc)-PDG| (GeV/c2)"}; + ConfigurableAxis deltaMassBins{"deltaMassBins", {300, 0.138, 0.198}, "SigmaC DeltaM histogram binning"}; + } sigmacSel; /// General options struct : ConfigurableGroup { Configurable cprDeltaEtaMax{"cprDeltaEtaMax", 0.01, "Max. Delta Eta for Close Pair Rejection"}; @@ -173,6 +187,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { using FilteredCharmCandDstars = soa::Filtered; using FilteredCharmCandDstar = FilteredCharmCandDstars::iterator; + using FilteredCharmCandScs = soa::Filtered; + using FilteredCharmMcCandScs = soa::Filtered>; using FilteredCharmMcCand3Prongs = soa::Filtered>; using FilteredCharmMcCand3Prong = FilteredCharmMcCand3Prongs::iterator; @@ -201,7 +217,10 @@ struct HfTaskCharmHadronsTrackFemtoDream { 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); + Filter hfMcSelFilter = ifnode(o2::framework::expressions::PlaceholderNode{false, std::string{"processMcScTrk"}}, + (nabs(aod::fdhf::flagMc) == static_cast(o2::hf_decay::hf_cand_sigmac::DecayChannelMain::Sc0ToPKPiPi) || + nabs(aod::fdhf::flagMc) == static_cast(o2::hf_decay::hf_cand_sigmac::DecayChannelMain::ScplusplusToPKPiPi)), + nabs(aod::fdhf::flagMc) == charmSel.charmHadMcSel); Filter trackEtaFilterLow = ifnode(aod::femtodreamparticle::partType == static_cast(aod::femtodreamparticle::ParticleType::kTrack), aod::femtodreamparticle::eta < trackSel.etaTrack1Max, true); Filter trackEtaFilterUp = ifnode(aod::femtodreamparticle::partType == static_cast(aod::femtodreamparticle::ParticleType::kTrack), aod::femtodreamparticle::eta > trackSel.etaTrack1Min, true); Filter trackPtFilterLow = ifnode(aod::femtodreamparticle::partType == static_cast(aod::femtodreamparticle::ParticleType::kTrack), aod::femtodreamparticle::pt < trackSel.ptTrack1Max, true); @@ -213,6 +232,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { Preslice perHf3ProngXicByCol = aod::femtodreamparticle::fdCollisionId; Preslice perHf2ProngByCol = aod::femtodreamparticle::fdCollisionId; Preslice perHfDstarByCol = aod::femtodreamparticle::fdCollisionId; + Preslice perHfScByCol = aod::femtodreamparticle::fdCollisionId; } preslices; /// Partition for particle 1 @@ -230,11 +250,17 @@ struct HfTaskCharmHadronsTrackFemtoDream { Partition partitionCharmHadron3ProngXic = 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 partitionCharmHadronSc = ((aod::fdhf::bdtBkg < 0.f && aod::fdhf::bdtPrompt < 0.f && aod::fdhf::bdtFD < 0.f) || + (aod::fdhf::bdtBkg < charmSel.charmHadBkgBDTmax && aod::fdhf::bdtFD < charmSel.charmHadFdBDTmax && aod::fdhf::bdtFD > charmSel.charmHadFdBDTmin && + aod::fdhf::bdtPrompt < charmSel.charmHadPromptBDTmax && aod::fdhf::bdtPrompt > charmSel.charmHadPromptBDTmin)); Partition partitionMcCharmHadron3Prong = aod::fdhf::originMcRec == OriginRecPrompt || aod::fdhf::originMcRec == OriginRecFD; Partition partitionMcCharmHadron3ProngXic = 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; + Partition partitionMcCharmHadronSc = (aod::fdhf::originMcRec == OriginRecPrompt || aod::fdhf::originMcRec == OriginRecFD) && ((aod::fdhf::bdtBkg < 0.f && aod::fdhf::bdtPrompt < 0.f && aod::fdhf::bdtFD < 0.f) || + (aod::fdhf::bdtBkg < charmSel.charmHadBkgBDTmax && aod::fdhf::bdtFD < charmSel.charmHadFdBDTmax && aod::fdhf::bdtFD > charmSel.charmHadFdBDTmin && + aod::fdhf::bdtPrompt < charmSel.charmHadPromptBDTmax && aod::fdhf::bdtPrompt > charmSel.charmHadPromptBDTmin)); /// Axis configurables ConfigurableAxis dummy{"dummy", {1, 0, 1}, "dummy axis"}; @@ -300,7 +326,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { void init(InitContext& /*context*/) { - std::array processes = {doprocessDataLcTrk, doprocessDataDplusTrk, doprocessDataD0Trk, doprocessDataDstarTrk, doprocessDataXicTrk, doprocessMcLcTrk, doprocessMcDplusTrk, doprocessMcD0Trk, doprocessMcDstarTrk, doprocessMcXicTrk}; + std::array processes = {doprocessDataScTrk, doprocessMcScTrk, doprocessDataLcTrk, doprocessDataDplusTrk, doprocessDataD0Trk, doprocessDataDstarTrk, doprocessDataXicTrk, doprocessMcLcTrk, doprocessMcDplusTrk, doprocessMcD0Trk, doprocessMcDstarTrk, doprocessMcXicTrk}; if (std::accumulate(processes.begin(), processes.end(), 0) != 1) { LOGP(fatal, "One and only one process function must be enabled at a time."); } @@ -308,8 +334,17 @@ struct HfTaskCharmHadronsTrackFemtoDream { bool process2Prong = doprocessDataD0Trk || doprocessMcD0Trk; bool processDstar = doprocessDataDstarTrk || doprocessMcDstarTrk; + bool const processSc = doprocessDataScTrk || doprocessMcScTrk; + if (processSc) { + isMc.value = doprocessMcScTrk.value; + if (std::abs(sigmacSel.particleSign.value) > 1 || sigmacSel.deltaMassMin.value >= sigmacSel.deltaMassMax.value || sigmacSel.lcMassWindow.value <= 0.f) { + LOG(fatal) << "Invalid SigmaC sign or mass window."; + } + binInvMass.value = sigmacSel.deltaMassBins.value; + } + massOne = o2::analysis::femtoDream::getMass(trackSel.pdgCodeTrack1.value); - massTwo = o2::analysis::femtoDream::getMass(charmSel.charmHadPDGCode.value); + massTwo = processSc ? 0.f : o2::analysis::femtoDream::getMass(charmSel.charmHadPDGCode.value); // setup columnpolicy for binning colBinningMult = {{mixingBinVztx, mixingBinMult}}; @@ -319,6 +354,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { allTrackHisto.init(®istry, binmultTempFit, binMulPercentile, binpTTrack, binEta, binPhi, binTempFitVarTrack, binNSigmaTPC, binNSigmaTOF, binNSigmaTPCTOF, binTPCClusters, dummy, dummy, isMc, trackSel.pdgCodeTrack1, true); selectedTrackHisto.init(®istry, binmultTempFit, binMulPercentile, binpTTrack, binEta, binPhi, binTempFitVarTrack, binNSigmaTPC, binNSigmaTOF, binNSigmaTPCTOF, binTPCClusters, dummy, dummy, isMc, trackSel.pdgCodeTrack1, true); + // SigmaC masses are updated from each candidate's charge when filling pairs. + int const charmPdg = processSc ? Pdg::kSigmaCPlusPlus : charmSel.charmHadPDGCode.value; sameEventCont.init(®istry, binkstar, binpTTrack, binkT, binmT, mixingBinMult, mixingBinMultPercentile, bin4Dkstar, bin4DmT, bin4DMult, bin4DmultPercentile, @@ -326,7 +363,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { highkstarCut, smearingByOrigin, binInvMass); - sameEventCont.setPDGCodes(trackSel.pdgCodeTrack1.value, charmSel.charmHadPDGCode.value); + sameEventCont.setPDGCodes(trackSel.pdgCodeTrack1.value, charmPdg); mixedEventCont.init(®istry, binkstar, binpTTrack, binkT, binmT, mixingBinMult, mixingBinMultPercentile, bin4Dkstar, bin4DmT, bin4DMult, bin4DmultPercentile, @@ -334,7 +371,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { highkstarCut, smearingByOrigin, binInvMass); - mixedEventCont.setPDGCodes(trackSel.pdgCodeTrack1.value, charmSel.charmHadPDGCode.value); + mixedEventCont.setPDGCodes(trackSel.pdgCodeTrack1.value, charmPdg); 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}}); @@ -368,7 +405,9 @@ struct HfTaskCharmHadronsTrackFemtoDream { float getCharmHadronMass(const Candidate& cand, bool ReturnDaughMass = false) { float invMass = 0.0f; - if constexpr (Channel == DecayChannel::LcToPKPi) { + if constexpr (Channel == DecayChannel::ScToLcPi) { + return ReturnDaughMass ? cand.massLc() : cand.deltaMass(); + } else if constexpr (Channel == DecayChannel::LcToPKPi) { if (cand.candidateSelFlag() == 1) { invMass = cand.m(std::array{MassProton, MassKPlus, MassPiPlus}); return invMass; @@ -434,6 +473,11 @@ struct HfTaskCharmHadronsTrackFemtoDream { LOG(fatal) << "Invalid PDG code for track mass hypothesis: " << trkPDGCode; } + if constexpr (Channel == DecayChannel::ScToLcPi) { + auto pVecSc = RecoDecayPtEtaPhi::pVector(cand.pt(), cand.eta(), cand.phi()); + return RecoDecay::m(std::array{pVecSc, pVecTrack}, std::array{static_cast(cand.massLc() + cand.deltaMass()), trackMassHyp}); + } + // D0 → K π + track (2-prong) if constexpr (Channel == DecayChannel::D0ToPiK) { const auto pVecCharmTrk = std::array{pVecProng0, pVecProng1, pVecTrack}; @@ -486,6 +530,39 @@ struct HfTaskCharmHadronsTrackFemtoDream { return -1.f; } + template + bool isCharmCandidateSelected(Candidate const& candidate, float mass) const + { + if constexpr (Channel == DecayChannel::ScToLcPi) { + int const chargeAbs = std::abs(candidate.charge()); + return (chargeAbs == 0 || chargeAbs == ChargeDoublyCharged) && + (sigmacSel.particleSign.value == 0 || candidate.chargeLc() == sigmacSel.particleSign.value) && + mass >= sigmacSel.deltaMassMin.value && mass <= sigmacSel.deltaMassMax.value && + std::abs(candidate.massLc() - MassLambdaCPlus) <= sigmacSel.lcMassWindow.value && + candidate.pt() >= charmSel.charmHadMinPt.value && candidate.pt() <= charmSel.charmHadMaxPt.value; + } + return mass >= charmSel.charmHadMinInvMass.value && mass <= charmSel.charmHadMaxInvMass.value; + } + + template + float getCharmHadronPdgMass(Candidate const& candidate) const + { + if constexpr (Channel == DecayChannel::ScToLcPi) { + return candidate.charge() == 0 ? MassSigmaC0 : MassSigmaCPlusPlus; + } + return massTwo; + } + + template + void fillPairHistogram(Container& container, Track const& track, Candidate const& candidate, Collision const& col) + { + if constexpr (Channel == DecayChannel::ScToLcPi) { + int const pdg = candidate.charge() == 0 ? Pdg::kSigmaC0 : Pdg::kSigmaCPlusPlus; + container.setPDGCodes(trackSel.pdgCodeTrack1.value, pdg); + } + container.template setPair(track, candidate, col.multNtr(), col.multV0M(), pairQASetting.use4D, pairQASetting.extendedPlots, smearingByOrigin); + } + /// This function processes the same event and takes care of all the histogramming template void doSameEvent(CandType& sliceCharmHad, PartitionType& sliceTrk1, TableTracks const& parts, Collision const& col) @@ -538,6 +615,13 @@ struct HfTaskCharmHadronsTrackFemtoDream { } } + if constexpr (Channel == DecayChannel::ScToLcPi) { + if (p1.trackId() == p2.prong0Id() || p1.trackId() == p2.prong1Id() || + p1.trackId() == p2.prong2Id() || p1.trackId() == p2.softPiId()) { + continue; + } + } + if constexpr (Channel == DecayChannel::DstarToD0Pi) { if (p1.trackId() == p2.prong0Id() || p1.trackId() == p2.prong1Id() || p1.trackId() == p2.prong2Id()) { continue; @@ -553,13 +637,14 @@ struct HfTaskCharmHadronsTrackFemtoDream { } } - float kstar = FemtoDreamMath::getkstar(p1, massOne, p2, massTwo); + float const charmMass = getCharmHadronPdgMass(p2); + float kstar = FemtoDreamMath::getkstar(p1, massOne, p2, charmMass); if (kstar > highkstarCut) { continue; } float invMass = getCharmHadronMass(p2); - if (invMass < charmSel.charmHadMinInvMass || invMass > charmSel.charmHadMaxInvMass) { + if (!isCharmCandidateSelected(p2, invMass)) { continue; } @@ -567,7 +652,11 @@ struct HfTaskCharmHadronsTrackFemtoDream { continue; } - float deltaInvMassPair = getCharmHadronTrackMass(p2, p1, trackSel.pdgCodeTrack1.value) - invMass; + float massForPairDifference = invMass; + if constexpr (Channel == DecayChannel::ScToLcPi) { + massForPairDifference += p2.massLc(); + } + float deltaInvMassPair = getCharmHadronTrackMass(p2, p1, trackSel.pdgCodeTrack1.value) - massForPairDifference; // associated track charge float chargeTrack = 0.; @@ -577,7 +666,9 @@ struct HfTaskCharmHadronsTrackFemtoDream { chargeTrack = NegativeCharge; } int pairSign = 0; - if (chargeTrack == p2.charge()) { + if constexpr (Channel == DecayChannel::ScToLcPi) { + pairSign = p2.charge() == 0 ? PairNotDefined : (chargeTrack * p2.charge() > 0 ? LikeSignPair : UnLikeSignPair); + } else if (chargeTrack == p2.charge()) { pairSign = LikeSignPair; } else if (chargeTrack == -p2.charge()) { pairSign = UnLikeSignPair; @@ -595,7 +686,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { originType = p2.originMcRec(); } - rowFemtoResultPairs( + tables.rowFemtoResultPairs( invMass, p2.pt(), p1.pt(), @@ -603,8 +694,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { p2.bdtPrompt(), p2.bdtFD(), kstar, - FemtoDreamMath::getkT(p1, massOne, p2, massTwo), - FemtoDreamMath::getmT(p1, massOne, p2, massTwo), + FemtoDreamMath::getkT(p1, massOne, p2, charmMass), + FemtoDreamMath::getmT(p1, massOne, p2, charmMass), col.multNtr(), col.multV0M(), p2.charge(), @@ -614,7 +705,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { charmHadMc, originType); - sameEventCont.setPair(p1, p2, col.multNtr(), col.multV0M(), pairQASetting.use4D, pairQASetting.extendedPlots, smearingByOrigin); + fillPairHistogram(sameEventCont, p1, p2, col); } } @@ -676,6 +767,13 @@ struct HfTaskCharmHadronsTrackFemtoDream { } } + if constexpr (Channel == DecayChannel::ScToLcPi) { + if (p1.trackId() == p2.prong0Id() || p1.trackId() == p2.prong1Id() || + p1.trackId() == p2.prong2Id() || p1.trackId() == p2.softPiId()) { + continue; + } + } + if constexpr (Channel == DecayChannel::DstarToD0Pi) { if (pairQASetting.useCPR.value) { @@ -689,14 +787,15 @@ struct HfTaskCharmHadronsTrackFemtoDream { } } - float kstar = FemtoDreamMath::getkstar(p1, massOne, p2, massTwo); + float const charmMass = getCharmHadronPdgMass(p2); + float kstar = FemtoDreamMath::getkstar(p1, massOne, p2, charmMass); if (kstar > highkstarCut) { continue; } float invMass = getCharmHadronMass(p2); - if (invMass < charmSel.charmHadMinInvMass || invMass > charmSel.charmHadMaxInvMass) { + if (!isCharmCandidateSelected(p2, invMass)) { continue; } @@ -704,7 +803,11 @@ struct HfTaskCharmHadronsTrackFemtoDream { continue; } - float deltaInvMassPair = getCharmHadronTrackMass(p2, p1, trackSel.pdgCodeTrack1.value) - invMass; + float massForPairDifference = invMass; + if constexpr (Channel == DecayChannel::ScToLcPi) { + massForPairDifference += p2.massLc(); + } + float deltaInvMassPair = getCharmHadronTrackMass(p2, p1, trackSel.pdgCodeTrack1.value) - massForPairDifference; // associated track charge float chargeTrack = 0.; @@ -714,7 +817,9 @@ struct HfTaskCharmHadronsTrackFemtoDream { chargeTrack = NegativeCharge; } int pairSign = 0; - if (chargeTrack == p2.charge()) { + if constexpr (Channel == DecayChannel::ScToLcPi) { + pairSign = p2.charge() == 0 ? PairNotDefined : (chargeTrack * p2.charge() > 0 ? LikeSignPair : UnLikeSignPair); + } else if (chargeTrack == p2.charge()) { pairSign = LikeSignPair; } else if (chargeTrack == -p2.charge()) { pairSign = UnLikeSignPair; @@ -729,7 +834,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { originType = p2.originMcRec(); } - rowFemtoResultPairs( + tables.rowFemtoResultPairs( invMass, p2.pt(), p1.pt(), @@ -737,8 +842,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { p2.bdtPrompt(), p2.bdtFD(), kstar, - FemtoDreamMath::getkT(p1, massOne, p2, massTwo), - FemtoDreamMath::getmT(p1, massOne, p2, massTwo), + FemtoDreamMath::getkT(p1, massOne, p2, charmMass), + FemtoDreamMath::getmT(p1, massOne, p2, charmMass), collision1.multNtr(), collision1.multV0M(), p2.charge(), @@ -748,7 +853,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { charmHadMc, originType); - mixedEventCont.setPair(p1, p2, collision1.multNtr(), collision1.multV0M(), pairQASetting.use4D, pairQASetting.extendedPlots, smearingByOrigin); + fillPairHistogram(mixedEventCont, p1, p2, collision1); } } } @@ -766,6 +871,11 @@ struct HfTaskCharmHadronsTrackFemtoDream { HIST("CharmHadronQA/hPtVsMass"), part.pt(), invMass); + if constexpr (Channel == DecayChannel::ScToLcPi) { + if (!isCharmCandidateSelected(part, invMass)) { + continue; + } + } timeStamp = part.timeStamp(); int charmHadMc = 0; @@ -777,7 +887,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { if constexpr (Channel == DecayChannel::DplusToPiKPi || Channel == DecayChannel::LcToPKPi) { - rowFemtoResultCharm3Prong( + tables.rowFemtoResultCharm3Prong( col.globalIndex(), timeStamp, invMass, @@ -793,8 +903,29 @@ struct HfTaskCharmHadronsTrackFemtoDream { part.bdtFD(), charmHadMc, originType); + } else if constexpr (Channel == DecayChannel::ScToLcPi) { + tables.rowFemtoResultCharmSc(col.globalIndex(), + timeStamp, + invMass, + part.massLc(), + part.pt(), + part.eta(), + part.phi(), + part.prong0Id(), + part.prong1Id(), + part.prong2Id(), + part.softPiId(), + part.chargeLc(), + part.sourceScId(), + part.candidateSelFlag(), + part.charge(), + part.bdtBkg(), + part.bdtPrompt(), + part.bdtFD(), + charmHadMc, + originType); } else if constexpr (Channel == DecayChannel::XicToXiPiPi) { - rowFemtoResultCharm3Prong( + tables.rowFemtoResultCharm3Prong( col.globalIndex(), timeStamp, invMass, @@ -811,7 +942,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { charmHadMc, originType); } else if constexpr (Channel == DecayChannel::D0ToPiK) { - rowFemtoResultCharm2Prong( + tables.rowFemtoResultCharm2Prong( col.globalIndex(), timeStamp, invMass, @@ -828,7 +959,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { originType); } else if constexpr (Channel == DecayChannel::DstarToD0Pi) { float invMassD0 = getCharmHadronMass(part, true); - rowFemtoResultCharmDstar( + tables.rowFemtoResultCharmDstar( col.globalIndex(), timeStamp, invMass, @@ -884,7 +1015,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { << trackSel.pdgCodeTrack1.value; break; } - rowFemtoResultTrk( + tables.rowFemtoResultTrk( col.globalIndex(), timeStamp, part.pt(), @@ -901,7 +1032,7 @@ struct HfTaskCharmHadronsTrackFemtoDream { // ---- Fill Collision Table ---- if (sliceCharmHad.size() > 0 || sliceTrk1.size() > 0) { - rowFemtoResultColl( + tables.rowFemtoResultColl( col.globalIndex(), timeStamp, col.posZ(), @@ -913,8 +1044,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { FilteredFDParticles const& parts, FilteredCharmCand3Prongs const& candidates) { - rowFemtoResultCharm3Prong.reserve(2 * candidates.size() + 1); - rowFemtoResultTrk.reserve(parts.size() + 1); + tables.rowFemtoResultCharm3Prong.reserve(2 * candidates.size() + 1); + tables.rowFemtoResultTrk.reserve(parts.size() + 1); for (const auto& col : cols) { eventHisto.fillQA(col); auto* partitionTrk1Selected = &partitionTrk1; @@ -957,8 +1088,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { FilteredFDParticles const& parts, FilteredCharmCand3Prongs const& candidates) { - rowFemtoResultCharm3Prong.reserve(candidates.size() + 1); - rowFemtoResultTrk.reserve(parts.size() + 1); + tables.rowFemtoResultCharm3Prong.reserve(candidates.size() + 1); + tables.rowFemtoResultTrk.reserve(parts.size() + 1); for (const auto& col : cols) { eventHisto.fillQA(col); auto* partitionTrk1Selected = &partitionTrk1; @@ -1002,8 +1133,8 @@ struct HfTaskCharmHadronsTrackFemtoDream { FilteredFDParticles const& parts, FilteredCharmCand2Prongs const& candidates) { - rowFemtoResultCharm2Prong.reserve(candidates.size() * 2 + 1); - rowFemtoResultTrk.reserve(parts.size() + 1); + tables.rowFemtoResultCharm2Prong.reserve(candidates.size() * 2 + 1); + tables.rowFemtoResultTrk.reserve(parts.size() + 1); for (const auto& col : cols) { eventHisto.fillQA(col); auto* partitionTrk1Selected = &partitionTrk1; @@ -1042,12 +1173,56 @@ struct HfTaskCharmHadronsTrackFemtoDream { } PROCESS_SWITCH(HfTaskCharmHadronsTrackFemtoDream, processDataD0Trk, "Enable processing D0ToPiK and Tracks correlation", false); + void processDataScTrk(FilteredCollisions const& cols, + FilteredFDParticles const& parts, + FilteredCharmCandScs const& candidates) + { + tables.rowFemtoResultCharmSc.reserve(candidates.size() + 1); + tables.rowFemtoResultTrk.reserve(parts.size() + 1); + for (const auto& col : cols) { + eventHisto.fillQA(col); + auto* partitionTrk1Selected = &partitionTrk1; + if (trackSel.pdgCodeTrack1.value == kKPlus) { + partitionTrk1Selected = &partitionTrk1Ka; + } + auto sliceTrk1 = partitionTrk1Selected->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + auto sliceCharmHad = partitionCharmHadronSc->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + if (fillTableWithCharm.value && sliceCharmHad.size() == 0) { + continue; + } + fillTables(col, sliceTrk1, sliceCharmHad); + if (sliceCharmHad.size() > 0 && sliceTrk1.size() > 0) { + doSameEvent(sliceCharmHad, sliceTrk1, parts, col); + } + } + if (mixSetting.doMixEvent) { + auto* partitionTrk1Selected = &partitionTrk1; + if (trackSel.pdgCodeTrack1.value == kKPlus) { + partitionTrk1Selected = &partitionTrk1Ka; + } + switch (mixSetting.mixingBinPolicy) { + case femtodreamcollision::kMult: + doMixedEvent(cols, partitionCharmHadronSc, *partitionTrk1Selected, parts, colBinningMult); + break; + case femtodreamcollision::kMultPercentile: + doMixedEvent(cols, partitionCharmHadronSc, *partitionTrk1Selected, parts, colBinningMultPercentile); + break; + case femtodreamcollision::kMultMultPercentile: + doMixedEvent(cols, partitionCharmHadronSc, *partitionTrk1Selected, parts, colBinningMultMultPercentile); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + PROCESS_SWITCH(HfTaskCharmHadronsTrackFemtoDream, processDataScTrk, "Enable processing ScToLcPi and Tracks correlation", false); + void processDataDstarTrk(FilteredCollisions const& cols, FilteredFDParticles const& parts, FilteredCharmCandDstars const& candidates) { - rowFemtoResultCharmDstar.reserve(candidates.size() + 1); - rowFemtoResultTrk.reserve(parts.size() + 1); + tables.rowFemtoResultCharmDstar.reserve(candidates.size() + 1); + tables.rowFemtoResultTrk.reserve(parts.size() + 1); for (const auto& col : cols) { eventHisto.fillQA(col); auto* partitionTrk1Selected = &partitionTrk1; @@ -1263,6 +1438,51 @@ struct HfTaskCharmHadronsTrackFemtoDream { } PROCESS_SWITCH(HfTaskCharmHadronsTrackFemtoDream, processMcD0Trk, "Enable processing D0ToPiK and Tracks correlation for Monte Carlo", false); + void processMcScTrk(FilteredMcColisions const& cols, + FilteredFDMcParts const& parts, + o2::aod::FDMCParticles const&, + o2::aod::FDExtMCParticles const&, + FilteredCharmMcCandScs const&) + { + for (const auto& col : cols) { + eventHisto.fillQA(col); + auto* partitionTrk1Selected = &partitionMcTrk1; + if (trackSel.pdgCodeTrack1.value == kKPlus) { + partitionTrk1Selected = &partitionMcTrk1Ka; + } + auto sliceMcTrk1 = partitionTrk1Selected->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + auto sliceMcCharmHad = partitionMcCharmHadronSc->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); + if ((col.bitmaskTrackOne() & bitMask) != bitMask || (col.bitmaskTrackTwo() & bitMask) != bitMask) { + continue; + } + if (fillTableWithCharm.value && sliceMcCharmHad.size() == 0) { + continue; + } + fillTables(col, sliceMcTrk1, sliceMcCharmHad); + doSameEvent(sliceMcCharmHad, sliceMcTrk1, parts, col); + } + if (mixSetting.doMixEvent) { + auto* partitionTrk1Selected = &partitionMcTrk1; + if (trackSel.pdgCodeTrack1.value == kKPlus) { + partitionTrk1Selected = &partitionMcTrk1Ka; + } + switch (mixSetting.mixingBinPolicy) { + case femtodreamcollision::kMult: + doMixedEvent(cols, partitionMcCharmHadronSc, *partitionTrk1Selected, parts, colBinningMult); + break; + case femtodreamcollision::kMultPercentile: + doMixedEvent(cols, partitionMcCharmHadronSc, *partitionTrk1Selected, parts, colBinningMultPercentile); + break; + case femtodreamcollision::kMultMultPercentile: + doMixedEvent(cols, partitionMcCharmHadronSc, *partitionTrk1Selected, parts, colBinningMultMultPercentile); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + PROCESS_SWITCH(HfTaskCharmHadronsTrackFemtoDream, processMcScTrk, "Enable processing ScToLcPi and Tracks correlation for Monte Carlo", false); + void processMcDstarTrk(FilteredMcColisions const& cols, FilteredFDMcParts const& parts, o2::aod::FDMCParticles const&,