diff --git a/PWGHF/D2H/Tasks/taskDeuteronFromLb.cxx b/PWGHF/D2H/Tasks/taskDeuteronFromLb.cxx index d50ee5c5c1b..496ee9c4176 100644 --- a/PWGHF/D2H/Tasks/taskDeuteronFromLb.cxx +++ b/PWGHF/D2H/Tasks/taskDeuteronFromLb.cxx @@ -59,7 +59,6 @@ struct HfTaskDeuteronFromLb { Configurable applySkimming{"applySkimming", false, "Skimmed dataset processing"}; Configurable cfgSkimming{"cfgSkimming", "fH2fromLb", "Configurable for skimming"}; Configurable sel8{"sel8", true, "sel8 event selection"}; - Configurable separateAntideuterons{"separateAntideuterons", true, "Fill FromLb histos for primary deuterons whose mother is Lb, Primary histos for all others"}; Configurable cfgEta{"cfgEta", 0.8f, "Track eta selection"}; Configurable cfgTPCNclsFound{"cfgTPCNclsFound", 100, "Minimum TPC clusters found"}; Configurable cfgTPCChi2Ncl{"cfgTPCChi2Ncl", 4.0f, "Maximum TPC chi2 per N clusters"}; @@ -73,29 +72,31 @@ struct HfTaskDeuteronFromLb { Configurable cfgDCAmin{"cfgDCAmin", 0.05f, "Minimum DCA for deuteron PID"}; Configurable cfgDCAmax{"cfgDCAmax", 1000.0f, "Maximum DCA for deuteron PID"}; Configurable rapidityCut{"rapidityCut", 0.5f, "Rapidity cut"}; - Configurable ptThresholdforPID{"ptThresholdforPID", 1.0f, "pT threshold for PID, above this value TPC+TOF PID is used, below only TPC PID is used"}; + Configurable ptThresholdforPID{"ptThresholdforPID", 1.0f, "pT threshold for PID, above this value TPC+TOF PID is used, below TPC+ITS PID is used"}; // PDG codes Configurable pdgCodeBeautyMeson{"pdgCodeBeautyMeson", -521, "PDG code of the beauty meson mother particle (default: B-)"}; Configurable pdgCodeBeautyBaryon{"pdgCodeBeautyBaryon", -5122, "PDG code of the beauty baryon mother particle (default: anti-Lambda_b)"}; Configurable pdgCodeDaughter{"pdgCodeDaughter", -1000010020, "PDG code of the daughter particle (default: anti-deuteron)"}; + Configurable pdgCodeHypertriton{"pdgCodeHypertriton", -1010010030, "PDG code of the hypertriton mother particle (default: anti-hypertriton)"}; int mRunNumber = 0; - float dBz = 0.f; int mCurrentRun = -1; framework::Service ccdb{}; using CollisionCandidates = o2::soa::Join; - using MCTrackCandidates = o2::soa::Join; + using TrackCandidates = o2::soa::Join; + using MCTrackCandidates = o2::soa::Join; using MCCollisionCandidates = o2::soa::Join; - using TrackCandidates = o2::soa::Join; Preslice trackIndicesPerCollision = o2::aod::track_association::collisionId; + Preslice mcParticlesPerCollision = o2::aod::mcparticle::mcCollisionId; ConfigurableAxis ptAxis{"ptAxis", {100, 0.f, 10.f}, "p_{T} GeV/c"}; ConfigurableAxis nSigmaAxis{"nSigmaAxis", {200, -10.f, 10.f}, "nSigma"}; - ConfigurableAxis dcaXyAxis{"dcaXyAxis", {1000, -0.2f, 0.2f}, "DCA xy (cm)"}; - ConfigurableAxis dcaZAxis{"dcaZAxis", {1000, -0.2f, 0.2f}, "DCA z (cm)"}; + ConfigurableAxis dcaXyAxis{"dcaXyAxis", {4000, -2.0f, 2.0f}, "DCA xy (cm)"}; + ConfigurableAxis dcaZAxis{"dcaZAxis", {4000, -2.0f, 2.0f}, "DCA z (cm)"}; + ConfigurableAxis ctauHypertritonAxis{"ctauHypertritonAxis", {100, 0.f, 30.f}, "c#tau (cm)"}; HistogramRegistry qaHistos{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -103,12 +104,12 @@ struct HfTaskDeuteronFromLb { OutputObj zorroSummary{"zorroSummary"}; OutputObj hProcessedEvents{TH1D("hProcessedEvents", "Event filtered;; Number of events", 4, 0., 4.)}; - bool isSelectedBeautyHadron(int pdgCode) + bool isSelectedMother(int pdgCode) { - return pdgCode == pdgCodeBeautyMeson || pdgCode == pdgCodeBeautyBaryon; + return pdgCode == pdgCodeBeautyMeson || pdgCode == pdgCodeBeautyBaryon || pdgCode == pdgCodeHypertriton; } - double getBeautyHadronMass(int pdgCode) + double getMotherMass(int pdgCode) { if (pdgCode == pdgCodeBeautyMeson) { return o2::constants::physics::MassBPlus; @@ -117,6 +118,9 @@ struct HfTaskDeuteronFromLb { if (pdgCode == pdgCodeBeautyBaryon) { return o2::constants::physics::MassLambdaB0; } + if (pdgCode == pdgCodeHypertriton) { + return o2::constants::physics::MassHyperTriton; + } return -1.; } @@ -139,54 +143,76 @@ struct HfTaskDeuteronFromLb { qaHistos.add("Data/etaAntideuteron", "etaAntideuteron", {HistType::kTH1F, {{100, -1.0f, 1.0f, "eta #bar{d}"}}}); qaHistos.add("Data/hVtxZ", "Z-Vertex distribution after selection;Z (cm)", HistType::kTH1F, {{100, -50, 50}}); qaHistos.add("Data/hnSigmaTOFVsPtPurity", "n#sigma TOF vs p_{T} for #bar{d} hypothesis for Data-driven purity check; p_{T} (GeV/c); n#sigma TOF", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); - qaHistos.add("Data/hnSigmaTOFVsPtPurityVeto", "n#sigma TOF vs p_{T} for #bar{d} hypothesis for Data-driven purity check with TPC veto for pi,K,p; p_{T} (GeV/c); n#sigma TOF", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); + qaHistos.add("Data/hnSigmaTOFVsPtPurityVeto", "n#sigma TOF vs p_{T} for #bar{d} hypothesis for Data-driven purity check with TPC veto for e,pi,K,p; p_{T} (GeV/c); n#sigma TOF", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); + qaHistos.add("Data/hnSigmaTOFVsPtPurityVetoTPCTOFClean", "n#sigma TOF vs p_{T} for #bar{d} hypothesis for Data-driven purity check with TPC and TOF veto for e,pi,K,p; p_{T} (GeV/c); n#sigma TOF", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); qaHistos.add("Data/hnSigmaTPCVsPtPurity", "n#sigma TPC vs p_{T} for #bar{d} hypothesis for Data-driven purity check; p_{T} (GeV/c); n#sigma TPC", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); qaHistos.add("Data/hnSigmaITSVsPt", "n#sigma ITS vs p_{T} for #bar{d} hypothesis, used to clean up the TPC/TOF purity sample; p_{T} (GeV/c); n#sigma ITS", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); // MC generated-level histograms qaHistos.add("MCGen/ptGeneratedBminus", "p_{T} generated B^{-};p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCGen/ptGeneratedAntiLambdaB", "p_{T} generated #bar{#Lambda}_{b};p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); - qaHistos.add("MCGen/ptAntiDeuteronPrimary", "p_{T} #bar{d} primary gen;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); - qaHistos.add("MCGen/ptAntiDeuteronFromBeautyHadron", "p_{T} #bar{d} from beauty gen hadrons;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); - qaHistos.add("MCGen/hMotherPdgCode", "PDG code of mother, gen level;PDG code;Counts", HistType::kTH1I, {{12000, -6000, 6000}}); qaHistos.add("MCGen/ctauBminus", "ctau of B^{-}, gen level;ctau (#mu m);Counts", HistType::kTH1F, {{25, 0., 2000.f}}); qaHistos.add("MCGen/ctauAntiLambdaB", "ctau of #bar{#Lambda}_{b}, gen level;ctau (#mu m);Counts", HistType::kTH1F, {{25, 0., 2000.f}}); - - // MC reco/MC-anchored histograms - qaHistos.add("MCReco/ptAntiDeuteronFromBminus", - "p_{T} #bar{d} from B^{-} reco/MC anchored;p_{T} (GeV/c);Counts", - HistType::kTH1F, {ptAxis}); - - qaHistos.add("MCReco/ptAntiDeuteronFromAntiLambdaB", - "p_{T} #bar{d} from #bar{#Lambda}_{b} reco/MC anchored;p_{T} (GeV/c);Counts", - HistType::kTH1F, {ptAxis}); - qaHistos.add( - "MCReco/ptAntiDeuteronPrimary", - "p_{T} #bar{d} primary reco/MC anchored;p_{T} (GeV/c);Counts", - HistType::kTH1F, - {ptAxis}); qaHistos.add("MCGen/ptAntiDeuteronFromBminus", "p_{T} #bar{d} from B^{-} gen;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCGen/ptAntiDeuteronFromAntiLambdaB", "p_{T} #bar{d} from #bar{#Lambda}_{b} gen;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCGen/ptAntiDeuteronAll", "p_{T} #bar{d} all gen;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); - qaHistos.add("MCReco/hDCAxy-Primary", "DCAxy primary reco/MC anchored;DCA xy (cm);Counts", {HistType::kTH1D, {{400, -0.2f, 0.2f, "DCA xy (cm)"}}}); - qaHistos.add("MCReco/hDCAxy-FromBeautyHadron", "DCAxy from beauty reco/MC anchored;DCA xy (cm);Counts", {HistType::kTH1D, {{400, -0.2f, 0.2f, "DCA xy (cm)"}}}); + qaHistos.add("MCGen/hEvents", "Events MC gen;;Counts", HistType::kTH1F, {{2, 0., 2.}}); + qaHistos.add("MCGen/ptGeneratedAntiHypertriton", "p_{T} generated ^{3}_{#bar{#Lambda}}#bar{H};p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCGen/ctauAntiHypertriton", "ctau of ^{3}_{#bar{#Lambda}}#bar{H}, gen level;c#tau (cm);Counts", HistType::kTH1F, {ctauHypertritonAxis}); + qaHistos.add("MCGen/ptAntiDeuteronFromAntiHypertriton", "p_{T} #bar{d} from ^{3}_{#bar{#Lambda}}#bar{H} gen;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + // MC reco/MC-anchored histograms + qaHistos.add("MCReco/hDCAxyVsPtDeFromBminus", "DCAxy #bar{d} from B^{-} vs p_{T};p_{T} (GeV/c);DCA_{xy} (cm)", {HistType::kTH2D, {ptAxis, dcaXyAxis}}); + qaHistos.add("MCReco/hDCAxyVsPtDeFromAntiLambdaB", "DCAxy #bar{d} from #bar{#Lambda}_{b} vs p_{T};p_{T} (GeV/c);DCA_{xy} (cm)", {HistType::kTH2D, {ptAxis, dcaXyAxis}}); + qaHistos.add("MCReco/hDCAxyVsPtDePrimary", "DCAxy #bar{d} primary vs p_{T};p_{T} (GeV/c);DCA_{xy} (cm)", {HistType::kTH2D, {ptAxis, dcaXyAxis}}); + qaHistos.add("MCReco/hDCAzVsPtDeFromBminus", "DCAz #bar{d} from B^{-} vs p_{T};p_{T} (GeV/c);DCA_{z} (cm)", {HistType::kTH2D, {ptAxis, dcaZAxis}}); + qaHistos.add("MCReco/hDCAzVsPtDeFromAntiLambdaB", "DCAz #bar{d} from #bar{#Lambda}_{b} vs p_{T};p_{T} (GeV/c);DCA_{z} (cm)", {HistType::kTH2D, {ptAxis, dcaZAxis}}); + qaHistos.add("MCReco/hDCAzVsPtDePrimary", "DCAz #bar{d} primary vs p_{T};p_{T} (GeV/c);DCA_{z} (cm)", {HistType::kTH2D, {ptAxis, dcaZAxis}}); + qaHistos.add("MCReco/ptAntiDeuteronFromBminus", "p_{T} #bar{d} from B^{-} reco/MC anchored;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/ptAntiDeuteronFromAntiLambdaB", "p_{T} #bar{d} from #bar{#Lambda}_{b} reco/MC anchored;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/ptAntiDeuteronPrimary", "p_{T} #bar{d} primary reco/MC anchored;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/hEvents", "Events MC reco;;Counts", HistType::kTH1F, {{3, 0., 3.}}); + qaHistos.add("MCReco/hDCAxy-Primary", "DCAxy primary reco/MC anchored;DCA xy (cm);Counts", {HistType::kTH1D, {dcaXyAxis}}); + qaHistos.add("MCReco/hDCAxy-FromBMeson", "DCAxy from B^{-} reco/MC anchored;DCA xy (cm);Counts", {HistType::kTH1D, {dcaXyAxis}}); + qaHistos.add("MCReco/hDCAxy-FromAntiLambdaB", "DCAxy from #bar{#Lambda}_{b} reco/MC anchored;DCA xy (cm);Counts", {HistType::kTH1D, {dcaXyAxis}}); qaHistos.add("MCReco/hMotherPdgCode", "PDG code of mother, reco/MC anchored;PDG code;Counts", HistType::kTH1I, {{12000, -6000, 6000}}); qaHistos.add("MCReco/ctauBminus", "ctau of B^{-}, reco/MC anchored;ctau (#mu m);Counts", HistType::kTH1F, {{25, 0., 2000.f}}); qaHistos.add("MCReco/ctauAntiLambdaB", "ctau of #bar{#Lambda}_{b}, reco/MC anchored;ctau (#mu m);Counts", HistType::kTH1F, {{25, 0., 2000.f}}); qaHistos.add("MCReco/ptAntiDeuteronPIDSelected", "p_{T} #bar{d} reco/MC anchored PID selection;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCReco/ptAntiDeuteronPIDSelectedTrue", "p_{T} #bar{d} reco/MC anchored PID and PDG;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); - qaHistos.add("MCReco/hnSigmaTPCVsPt", "n#sigma TPC vs p_{T} for #bar{d} hypothesis for MC; p_{T} (GeV/c); n#sigma TPC", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); - qaHistos.add("MCReco/hnSigmaTOFVsPt", "n#sigma TOF vs p_{T} for #bar{d} hypothesis for MC; p_{T} (GeV/c); n#sigma TOF", {HistType::kTH2D, {ptAxis, nSigmaAxis}}); qaHistos.add("MCReco/ptAntiDeuteronTPCOnlyPIDSelected", "p_{T} #bar{d} reco/MC anchored TPC only PID selection;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCReco/ptAntiDeuteronTPCOnlyPIDSelectedTrue", "p_{T} #bar{d} reco/MC anchored TPC only PID and PDG;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCReco/ptAntiDeuteronTPCorTPCTOFPIDSelected", "p_{T} #bar{d} reco/MC anchored TPC or TPC+TOF PID selection;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); qaHistos.add("MCReco/ptAntiDeuteronTPCorTPCTOFPIDSelectedTrue", "p_{T} #bar{d} reco/MC anchored TPC or TPC+TOF PID and PDG;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/ptAntiDeuteronGenFromBminus", "p_{T} #bar{d} from B^{-} gen in reco events;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/ptAntiDeuteronGenFromAntiLambdaB", "p_{T} #bar{d} from #bar{#Lambda}_{b} gen in reco events;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/hDCAxyVsPtDeFromAntiHypertriton", "DCAxy #bar{d} from ^{3}_{#bar{#Lambda}}#bar{H} vs p_{T};p_{T} (GeV/c);DCA_{xy} (cm)", {HistType::kTH2D, {ptAxis, dcaXyAxis}}); + qaHistos.add("MCReco/hDCAzVsPtDeFromAntiHypertriton", "DCAz #bar{d} from ^{3}_{#bar{#Lambda}}#bar{H} vs p_{T};p_{T} (GeV/c);DCA_{z} (cm)", {HistType::kTH2D, {ptAxis, dcaZAxis}}); + qaHistos.add("MCReco/ptAntiDeuteronFromAntiHypertriton", "p_{T} #bar{d} from ^{3}_{#bar{#Lambda}}#bar{H} reco/MC anchored;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/hDCAxy-FromAntiHypertriton", "DCAxy from ^{3}_{#bar{#Lambda}}#bar{H} reco/MC anchored;DCA xy (cm);Counts", {HistType::kTH1D, {dcaXyAxis}}); + qaHistos.add("MCReco/ctauAntiHypertriton", "ctau of ^{3}_{#bar{#Lambda}}#bar{H}, reco/MC anchored;c#tau (cm);Counts", HistType::kTH1F, {ctauHypertritonAxis}); + qaHistos.add("MCReco/ptAntiDeuteronGenFromAntiHypertriton", "p_{T} #bar{d} from ^{3}_{#bar{#Lambda}}#bar{H} gen in reco events;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); + qaHistos.add("MCReco/ptAntiDeuteronGenPrimary", "p_{T} #bar{d} primary gen in reco events;p_{T} (GeV/c);Counts", HistType::kTH1F, {ptAxis}); hProcessedEvents->GetXaxis()->SetBinLabel(1, "Events processed"); hProcessedEvents->GetXaxis()->SetBinLabel(2, "ZORRO"); hProcessedEvents->GetXaxis()->SetBinLabel(3, "sel8"); hProcessedEvents->GetXaxis()->SetBinLabel(4, "z-vertex"); + auto hEvents = qaHistos.get(HIST("MCReco/hEvents")); + hEvents->GetXaxis()->SetBinLabel(1, "All"); + hEvents->GetXaxis()->SetBinLabel(2, "sel8"); + hEvents->GetXaxis()->SetBinLabel(3, "z-vertex"); + auto hEventsGen = qaHistos.get(HIST("MCGen/hEvents")); + hEventsGen->GetXaxis()->SetBinLabel(1, "All"); + hEventsGen->GetXaxis()->SetBinLabel(2, "z-vertex"); + + qaHistos.add("MCReco/hMotherCategory", "Mother of #bar{d};;Counts", HistType::kTH1F, {{5, 0., 5.}}); + auto hMotherCategory = qaHistos.get(HIST("MCReco/hMotherCategory")); + hMotherCategory->GetXaxis()->SetBinLabel(1, "B^{-}"); + hMotherCategory->GetXaxis()->SetBinLabel(2, "#bar{#Lambda}_{b}"); + hMotherCategory->GetXaxis()->SetBinLabel(3, "^{3}_{#bar{#Lambda}}#bar{H}"); + hMotherCategory->GetXaxis()->SetBinLabel(4, "primary"); + hMotherCategory->GetXaxis()->SetBinLabel(5, "other"); } void initCCDB(o2::aod::BCsWithTimestamps::iterator const& bc) @@ -194,12 +220,12 @@ struct HfTaskDeuteronFromLb { if (mRunNumber == bc.runNumber()) { return; } - dBz = 0.f; if (applySkimming) { zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), cfgSkimming.value); zorro.populateHistRegistry(registry, bc.runNumber()); } + mRunNumber = bc.runNumber(); } template @@ -208,6 +234,9 @@ struct HfTaskDeuteronFromLb { if (std::abs(track.eta()) > cfgEta) { return false; } + if (std::abs(track.rapidity(o2::constants::physics::MassDeuteron)) > rapidityCut) { + return false; + } if (!track.hasITS()) { return false; } @@ -244,6 +273,7 @@ struct HfTaskDeuteronFromLb { TrackCandidates const& tracks, o2::aod::BCsWithTimestamps const&) { + auto tracksWithItsPid = o2::soa::Attach(tracks); for (const auto& collision : collisions) { if (mCurrentRun != collision.bc_as().runNumber()) { auto* grpo = ccdb->getForTimeStamp("GLO/Config/GRPMagField", collision.bc_as().timestamp()); @@ -271,7 +301,6 @@ struct HfTaskDeuteronFromLb { qaHistos.fill(HIST("Data/hVtxZ"), collision.posZ()); const auto& trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); - auto tracksWithItsPid = o2::soa::Attach(tracks); for (const auto& trackId : trackIdsThisCollision) { const auto& track = tracksWithItsPid.rawIteratorAt(trackId.trackId()); @@ -296,10 +325,18 @@ struct HfTaskDeuteronFromLb { qaHistos.fill(HIST("Data/hnSigmaITSVsPt"), track.pt(), track.itsNSigmaDe()); - const bool isTPCDeClean = isTPCDe && - std::abs(track.tpcNSigmaPi()) > cfgTPCNsigma && - std::abs(track.tpcNSigmaKa()) > cfgTPCNsigma && - std::abs(track.tpcNSigmaPr()) > cfgTPCNsigma; + // TPC veto against pi/K/p/e + const bool passTPCVeto = std::abs(track.tpcNSigmaPi()) > cfgTPCNsigma && + std::abs(track.tpcNSigmaKa()) > cfgTPCNsigma && + std::abs(track.tpcNSigmaPr()) > cfgTPCNsigma && + std::abs(track.tpcNSigmaEl()) > cfgTPCNsigma; + // TOF veto against pi/K/p/e + const bool passTOFVeto = std::abs(track.tofNSigmaPi()) > cfgTofNsigma && + std::abs(track.tofNSigmaKa()) > cfgTofNsigma && + std::abs(track.tofNSigmaPr()) > cfgTofNsigma && + std::abs(track.tofNSigmaEl()) > cfgTofNsigma; + + const bool isTPCDeClean = isTPCDe && passTPCVeto; if (track.pt() < ptThresholdforPID) { if (isITSDe) { @@ -312,17 +349,15 @@ struct HfTaskDeuteronFromLb { if (isTPCDeClean && track.hasTOF()) { qaHistos.fill(HIST("Data/hnSigmaTOFVsPtPurityVeto"), track.pt(), track.tofNSigmaDe()); } + if (isTPCDeClean && passTOFVeto && track.hasTOF()) { + qaHistos.fill(HIST("Data/hnSigmaTOFVsPtPurityVetoTPCTOFClean"), track.pt(), track.tofNSigmaDe()); + } } - if (track.pt() < ptThresholdforPID && isTPCDe && isITSDe) { - qaHistos.fill(HIST("Data/ptAntiDeuteron"), track.pt()); - qaHistos.fill(HIST("Data/etaAntideuteron"), track.eta()); - qaHistos.fill(HIST("Data/hDCAxyVsPt"), track.pt(), dca[0]); - qaHistos.fill(HIST("Data/hDCAzVsPt"), track.pt(), dca[1]); - qaHistos.fill(HIST("Data/hnSigmaTPCVsPt"), track.pt(), track.tpcNSigmaDe()); - qaHistos.fill(HIST("Data/hnSigmaTOFVsPt"), track.pt(), track.tofNSigmaDe()); - } - if (track.pt() >= ptThresholdforPID && isTPCDe && isTOFDe && track.hasTOF()) { + const bool passFinalSel = (track.pt() < ptThresholdforPID && isTPCDe && isITSDe) || + (track.pt() >= ptThresholdforPID && isTPCDeClean && passTOFVeto && isTOFDe && track.hasTOF()); + + if (passFinalSel) { qaHistos.fill(HIST("Data/ptAntiDeuteron"), track.pt()); qaHistos.fill(HIST("Data/etaAntideuteron"), track.eta()); qaHistos.fill(HIST("Data/hDCAxyVsPt"), track.pt(), dca[0]); @@ -335,95 +370,200 @@ struct HfTaskDeuteronFromLb { } PROCESS_SWITCH(HfTaskDeuteronFromLb, processData, "processData", true); - void processMC(MCCollisionCandidates::iterator const& collision, MCTrackCandidates const& tracks, o2::aod::McParticles const&) + void processMC(MCCollisionCandidates const& collisions, + o2::aod::TrackAssoc const& trackIndices, + MCTrackCandidates const& tracks, + o2::aod::McParticles const& mcParticles, + o2::aod::BCsWithTimestamps const&) { - for (const auto& track : tracks) { - if (std::abs(track.dcaXY()) < cfgDCAmin || std::abs(track.dcaXY()) > cfgDCAmax) { - continue; + auto tracksWithItsPid = o2::soa::Attach(tracks); + tracksWithItsPid.bindExternalIndices(&mcParticles); + + for (const auto& collision : collisions) { + const auto& bc = collision.bc_as(); + if (mCurrentRun != bc.runNumber()) { + auto* grpo = ccdb->getForTimeStamp("GLO/Config/GRPMagField", bc.timestamp()); + o2::base::Propagator::initFieldFromGRP(grpo); + mCurrentRun = bc.runNumber(); } - if (!passedSingleTrackSelection(track)) { + + qaHistos.fill(HIST("MCReco/hEvents"), 0.5); + if (sel8 && !collision.sel8()) { continue; } - - if (!track.has_mcParticle()) { + qaHistos.fill(HIST("MCReco/hEvents"), 1.5); + if (std::abs(collision.posZ()) > cutzvertex) { continue; } + qaHistos.fill(HIST("MCReco/hEvents"), 2.5); - auto mcParticle = track.mcParticle(); - const bool isTrueDeuteron = (mcParticle.pdgCode() == pdgCodeDaughter); - const bool isTPCDe = std::abs(track.tpcNSigmaDe()) < cfgTPCNsigma; - const bool isTOFDe = std::abs(track.tofNSigmaDe()) < cfgTofNsigma; - const bool selTpcOrTpcTof = isTPCDe && (!track.hasTOF() || isTOFDe); + if (collision.has_mcCollision()) { + const auto mcParticlesThisEvent = mcParticles.sliceBy(mcParticlesPerCollision, collision.mcCollisionId()); - if (isTPCDe) { - qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCOnlyPIDSelected"), track.pt()); - if (isTrueDeuteron) { - qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCOnlyPIDSelectedTrue"), track.pt()); + for (const auto& mcParticle : mcParticlesThisEvent) { + + if (mcParticle.pdgCode() != pdgCodeDaughter) { + continue; + } + if (std::abs(mcParticle.eta()) > cfgEta || std::abs(mcParticle.y()) > rapidityCut) { + continue; + } + if (mcParticle.pt() < cfgMinPt || mcParticle.pt() > cfgMaxPt) { + continue; + } + bool fromMother = false; + for (const auto& mom : mcParticle.mothers_as()) { + if (mom.pdgCode() == pdgCodeBeautyMeson) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronGenFromBminus"), mcParticle.pt()); + fromMother = true; + break; + } + if (mom.pdgCode() == pdgCodeBeautyBaryon) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronGenFromAntiLambdaB"), mcParticle.pt()); + fromMother = true; + break; + } + if (mom.pdgCode() == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronGenFromAntiHypertriton"), mcParticle.pt()); + fromMother = true; + break; + } + } + if (!fromMother && mcParticle.isPhysicalPrimary()) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronGenPrimary"), mcParticle.pt()); + } } } - if (track.hasTOF()) { - qaHistos.fill(HIST("MCReco/hnSigmaTPCVsPt"), track.pt(), track.tpcNSigmaDe()); - qaHistos.fill(HIST("MCReco/hnSigmaTOFVsPt"), track.pt(), track.tofNSigmaDe()); - } + const auto& trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); - if (selTpcOrTpcTof) { - qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCorTPCTOFPIDSelected"), track.pt()); - if (isTrueDeuteron) { - qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCorTPCTOFPIDSelectedTrue"), track.pt()); + for (const auto& trackId : trackIdsThisCollision) { + const auto& track = tracksWithItsPid.rawIteratorAt(trackId.trackId()); + std::array dca{track.dcaXY(), track.dcaZ()}; + + if (track.collisionId() != collision.globalIndex()) { + auto trackPar = getTrackParCov(track); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackPar, 2.f, noMatCorr, &dca); } - } - if (isTPCDe && isTOFDe && track.hasTOF()) { - qaHistos.fill(HIST("MCReco/ptAntiDeuteronPIDSelected"), track.pt()); - if (isTrueDeuteron) { - qaHistos.fill(HIST("MCReco/ptAntiDeuteronPIDSelectedTrue"), track.pt()); + if (std::abs(dca[0]) < cfgDCAmin || std::abs(dca[0]) > cfgDCAmax) { + continue; + } + if (!passedSingleTrackSelection(track)) { + continue; } - } else { - continue; - } - if (isTrueDeuteron) { - if (std::abs(mcParticle.y()) > rapidityCut) { + if (!track.has_mcParticle()) { continue; } - if (mcParticle.isPhysicalPrimary()) { + + auto mcParticle = track.mcParticle(); + const bool isTrueDeuteron = (mcParticle.pdgCode() == pdgCodeDaughter); + const bool isTPCDe = std::abs(track.tpcNSigmaDe()) < cfgTPCNsigma; + const bool isTOFDe = std::abs(track.tofNSigmaDe()) < cfgTofNsigma; + const bool isITSDe = track.itsNSigmaDe() > cfgITSNsigma; + + const bool selTpcOrTpcTof = isTPCDe && (!track.hasTOF() || isTOFDe); + + if (isTPCDe) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCOnlyPIDSelected"), track.pt()); + if (isTrueDeuteron) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCOnlyPIDSelectedTrue"), track.pt()); + } + } + + if (selTpcOrTpcTof) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCorTPCTOFPIDSelected"), track.pt()); + if (isTrueDeuteron) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronTPCorTPCTOFPIDSelectedTrue"), track.pt()); + } + } + // Purity study: TPC && TOF + if (isTPCDe && isTOFDe && track.hasTOF()) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronPIDSelected"), track.pt()); + if (isTrueDeuteron) { + qaHistos.fill(HIST("MCReco/ptAntiDeuteronPIDSelectedTrue"), track.pt()); + } + } + + // TPC veto against pi/K/p/e + const bool passTPCVeto = std::abs(track.tpcNSigmaPi()) > cfgTPCNsigma && + std::abs(track.tpcNSigmaKa()) > cfgTPCNsigma && + std::abs(track.tpcNSigmaPr()) > cfgTPCNsigma && + std::abs(track.tpcNSigmaEl()) > cfgTPCNsigma; + // TOF veto against pi/K/p/e + const bool passTOFVeto = std::abs(track.tofNSigmaPi()) > cfgTofNsigma && + std::abs(track.tofNSigmaKa()) > cfgTofNsigma && + std::abs(track.tofNSigmaPr()) > cfgTofNsigma && + std::abs(track.tofNSigmaEl()) > cfgTofNsigma; + + const bool isTPCDeClean = isTPCDe && passTPCVeto; + + if (isTrueDeuteron) { int motherPdg = 0; - if (separateAntideuterons) { - for (const auto& mom : mcParticle.mothers_as()) { - qaHistos.fill(HIST("MCReco/hMotherPdgCode"), mom.pdgCode()); - if (mom.pdgCode() == pdgCodeBeautyMeson || mom.pdgCode() == pdgCodeBeautyBaryon) { - motherPdg = mom.pdgCode(); - double dx = mcParticle.vx() - collision.posX(); - double dy = mcParticle.vy() - collision.posY(); - double dz = mcParticle.vz() - collision.posZ(); - - double flightDistance = std::sqrt(dx * dx + dy * dy + dz * dz); - double massBeauty = getBeautyHadronMass(mom.pdgCode()); - - if (massBeauty > 0.) { - double ctauBeauty = flightDistance / mom.p() * massBeauty; - ctauBeauty *= 1.e4; - - if (motherPdg == pdgCodeBeautyMeson) { - qaHistos.fill(HIST("MCReco/ctauBminus"), ctauBeauty); - } else if (motherPdg == pdgCodeBeautyBaryon) { - qaHistos.fill(HIST("MCReco/ctauAntiLambdaB"), ctauBeauty); - } + for (const auto& mom : mcParticle.mothers_as()) { + qaHistos.fill(HIST("MCReco/hMotherPdgCode"), mom.pdgCode()); + if (mom.pdgCode() == pdgCodeBeautyMeson || mom.pdgCode() == pdgCodeBeautyBaryon || mom.pdgCode() == pdgCodeHypertriton) { + motherPdg = mom.pdgCode(); + double dx = mcParticle.vx() - mom.vx(); + double dy = mcParticle.vy() - mom.vy(); + double dz = mcParticle.vz() - mom.vz(); + + double flightDistance = std::sqrt(dx * dx + dy * dy + dz * dz); + double massBeauty = getMotherMass(mom.pdgCode()); + + if (massBeauty > 0.) { + double ctauBeauty = flightDistance / mom.p() * massBeauty; + + if (motherPdg == pdgCodeBeautyMeson) { + qaHistos.fill(HIST("MCReco/ctauBminus"), ctauBeauty * 1.e4); + } else if (motherPdg == pdgCodeBeautyBaryon) { + qaHistos.fill(HIST("MCReco/ctauAntiLambdaB"), ctauBeauty * 1.e4); + } else if (motherPdg == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCReco/ctauAntiHypertriton"), ctauBeauty); // cm } - break; } + break; } } + if (motherPdg == pdgCodeBeautyMeson) { - qaHistos.fill(HIST("MCReco/hDCAxy-FromBeautyHadron"), track.dcaXY()); - qaHistos.fill(HIST("MCReco/ptAntiDeuteronFromBminus"), track.pt()); + qaHistos.fill(HIST("MCReco/hMotherCategory"), 0.5); } else if (motherPdg == pdgCodeBeautyBaryon) { - qaHistos.fill(HIST("MCReco/hDCAxy-FromBeautyHadron"), track.dcaXY()); - qaHistos.fill(HIST("MCReco/ptAntiDeuteronFromAntiLambdaB"), track.pt()); + qaHistos.fill(HIST("MCReco/hMotherCategory"), 1.5); + } else if (motherPdg == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCReco/hMotherCategory"), 2.5); + } else if (mcParticle.isPhysicalPrimary()) { + qaHistos.fill(HIST("MCReco/hMotherCategory"), 3.5); } else { - qaHistos.fill(HIST("MCReco/hDCAxy-Primary"), track.dcaXY()); - qaHistos.fill(HIST("MCReco/ptAntiDeuteronPrimary"), track.pt()); + qaHistos.fill(HIST("MCReco/hMotherCategory"), 4.5); + } + + // same final selection as in data + const bool passFinalSel = (track.pt() < ptThresholdforPID && isTPCDe && isITSDe) || + (track.pt() >= ptThresholdforPID && isTPCDeClean && passTOFVeto && isTOFDe && track.hasTOF()); + if (passFinalSel) { + if (motherPdg == pdgCodeBeautyMeson) { + qaHistos.fill(HIST("MCReco/hDCAxyVsPtDeFromBminus"), track.pt(), dca[0]); + qaHistos.fill(HIST("MCReco/hDCAzVsPtDeFromBminus"), track.pt(), dca[1]); + qaHistos.fill(HIST("MCReco/hDCAxy-FromBMeson"), dca[0]); + qaHistos.fill(HIST("MCReco/ptAntiDeuteronFromBminus"), track.pt()); + } else if (motherPdg == pdgCodeBeautyBaryon) { + qaHistos.fill(HIST("MCReco/hDCAxyVsPtDeFromAntiLambdaB"), track.pt(), dca[0]); + qaHistos.fill(HIST("MCReco/hDCAzVsPtDeFromAntiLambdaB"), track.pt(), dca[1]); + qaHistos.fill(HIST("MCReco/hDCAxy-FromAntiLambdaB"), dca[0]); + qaHistos.fill(HIST("MCReco/ptAntiDeuteronFromAntiLambdaB"), track.pt()); + } else if (motherPdg == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCReco/hDCAxyVsPtDeFromAntiHypertriton"), track.pt(), dca[0]); + qaHistos.fill(HIST("MCReco/hDCAzVsPtDeFromAntiHypertriton"), track.pt(), dca[1]); + qaHistos.fill(HIST("MCReco/hDCAxy-FromAntiHypertriton"), dca[0]); + qaHistos.fill(HIST("MCReco/ptAntiDeuteronFromAntiHypertriton"), track.pt()); + } else if (mcParticle.isPhysicalPrimary()) { + qaHistos.fill(HIST("MCReco/hDCAxyVsPtDePrimary"), track.pt(), dca[0]); + qaHistos.fill(HIST("MCReco/hDCAzVsPtDePrimary"), track.pt(), dca[1]); + qaHistos.fill(HIST("MCReco/hDCAxy-Primary"), dca[0]); + qaHistos.fill(HIST("MCReco/ptAntiDeuteronPrimary"), track.pt()); + } } } } @@ -433,38 +573,39 @@ struct HfTaskDeuteronFromLb { void processGen(o2::aod::McCollision const& collision, o2::aod::McParticles const& mcParticles) { - hProcessedEvents->Fill(0.5); + qaHistos.fill(HIST("MCGen/hEvents"), 0.5); + if (std::abs(collision.posZ()) > cutzvertex) { + return; + } + qaHistos.fill(HIST("MCGen/hEvents"), 1.5); for (const auto& mcParticle : mcParticles) { - - // Beauty hadron mother - if (isSelectedBeautyHadron(mcParticle.pdgCode())) { - if (std::abs(mcParticle.y()) <= rapidityCut) { - if (mcParticle.pdgCode() == pdgCodeBeautyMeson) { - qaHistos.fill(HIST("MCGen/ptGeneratedBminus"), mcParticle.pt()); - } else if (mcParticle.pdgCode() == pdgCodeBeautyBaryon) { - qaHistos.fill(HIST("MCGen/ptGeneratedAntiLambdaB"), mcParticle.pt()); - } + if (isSelectedMother(mcParticle.pdgCode())) { + if (mcParticle.pdgCode() == pdgCodeBeautyMeson) { + qaHistos.fill(HIST("MCGen/ptGeneratedBminus"), mcParticle.pt()); + } else if (mcParticle.pdgCode() == pdgCodeBeautyBaryon) { + qaHistos.fill(HIST("MCGen/ptGeneratedAntiLambdaB"), mcParticle.pt()); + } else if (mcParticle.pdgCode() == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCGen/ptGeneratedAntiHypertriton"), mcParticle.pt()); } - if (mcParticle.has_daughters()) { for (const auto& daughter : mcParticle.daughters_as()) { - double dx = daughter.vx() - collision.posX(); - double dy = daughter.vy() - collision.posY(); - double dz = daughter.vz() - collision.posZ(); - + double dx = daughter.vx() - mcParticle.vx(); + double dy = daughter.vy() - mcParticle.vy(); + double dz = daughter.vz() - mcParticle.vz(); double flightDistance = std::sqrt(dx * dx + dy * dy + dz * dz); - double massBeauty = getBeautyHadronMass(mcParticle.pdgCode()); + double massBeauty = getMotherMass(mcParticle.pdgCode()); if (massBeauty > 0.) { double ctauBeauty = flightDistance / mcParticle.p() * massBeauty; - ctauBeauty *= 1.e4; // cm -> micrometers if (mcParticle.pdgCode() == pdgCodeBeautyMeson) { - qaHistos.fill(HIST("MCGen/ctauBminus"), ctauBeauty); + qaHistos.fill(HIST("MCGen/ctauBminus"), ctauBeauty * 1.e4); } else if (mcParticle.pdgCode() == pdgCodeBeautyBaryon) { - qaHistos.fill(HIST("MCGen/ctauAntiLambdaB"), ctauBeauty); + qaHistos.fill(HIST("MCGen/ctauAntiLambdaB"), ctauBeauty * 1.e4); + } else if (mcParticle.pdgCode() == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCGen/ctauAntiHypertriton"), ctauBeauty); // cm } } break; @@ -473,14 +614,17 @@ struct HfTaskDeuteronFromLb { } if (mcParticle.pdgCode() == pdgCodeDaughter) { - if (std::abs(mcParticle.eta()) < cfgEta) { - qaHistos.fill(HIST("MCGen/ptAntiDeuteronAll"), mcParticle.pt()); + + if (std::abs(mcParticle.eta()) > cfgEta) { + continue; } if (std::abs(mcParticle.y()) > rapidityCut) { continue; } + qaHistos.fill(HIST("MCGen/ptAntiDeuteronAll"), mcParticle.pt()); + int motherPdg = 0; if (mcParticle.has_mothers()) { for (const auto& mom : mcParticle.mothers_as()) { @@ -496,6 +640,11 @@ struct HfTaskDeuteronFromLb { motherPdg = mom.pdgCode(); break; } + + if (mom.pdgCode() == pdgCodeHypertriton) { + motherPdg = mom.pdgCode(); + break; + } } } @@ -503,6 +652,8 @@ struct HfTaskDeuteronFromLb { qaHistos.fill(HIST("MCGen/ptAntiDeuteronFromBminus"), mcParticle.pt()); } else if (motherPdg == pdgCodeBeautyBaryon) { qaHistos.fill(HIST("MCGen/ptAntiDeuteronFromAntiLambdaB"), mcParticle.pt()); + } else if (motherPdg == pdgCodeHypertriton) { + qaHistos.fill(HIST("MCGen/ptAntiDeuteronFromAntiHypertriton"), mcParticle.pt()); } else if (mcParticle.isPhysicalPrimary()) { qaHistos.fill(HIST("MCGen/ptAntiDeuteronPrimary"), mcParticle.pt()); }