[Trigger,PWGHF] Add new trigger for non-prompt antideuterons from beauty

EventFiltering/CMakeLists.txt

@@ -107,7 +107,12 @@ o2physics_add_dpl_workflow(lf-f1proton-filter
     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore KFParticle::KFParticle O2::DetectorsBase
     COMPONENT_NAME Analysis)
 
- o2physics_add_dpl_workflow(lf-doublephi-filter
+o2physics_add_dpl_workflow(lf-doublephi-filter
     SOURCES PWGLF/filterdoublephi.cxx
     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase
     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(h2-from-beauty-filter
+    SOURCES PWGHF/H2fromLbFilter.cxx
+    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase
+    COMPONENT_NAME Analysis)

EventFiltering/PWGHF/H2fromLbFilter.cxx

@@ -0,0 +1,246 @@
+// 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.
+//
+/// \brief A filter task for non prompt deuterons
+/// \author Marta Razza [email protected]
+/// \author Francesca Ercolessi [email protected]
+/// \since Dec 17, 2025
+
+#include "../filterTables.h"
+
+#include "PWGLF/DataModel/LFParticleIdentification.h"
+
+#include "Common/Core/RecoDecay.h"
+#include "Common/Core/TrackSelection.h"
+#include "Common/Core/trackUtilities.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/CollisionAssociationTables.h"
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/Multiplicity.h"
+#include "Common/DataModel/PIDResponseTOF.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+
+#include "CCDB/BasicCCDBManager.h"
+#include "CommonConstants/PhysicsConstants.h"
+#include "DCAFitter/DCAFitterN.h"
+#include "DataFormatsParameters/GRPMagField.h"
+#include "DataFormatsParameters/GRPObject.h"
+#include "DetectorsBase/Propagator.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+#include "ReconstructionDataFormats/Track.h"
+#include "ReconstructionDataFormats/TrackParametrization.h"
+
+#include "TVector3.h"
+
+#include <cmath>
+
+struct H2fromLbFilter {
+
+  // Recall the output table
+  o2::framework::Produces<o2::aod::H2fromLbFilters> table;
+
+  o2::framework::Configurable<int> nBinsPt{"nBinsPt", 100, "N bins in pT histo"};
+
+  o2::base::Propagator::MatCorrType noMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE;
+
+  // Define a histograms and registries
+  o2::framework::HistogramRegistry QAHistos{"QAHistos", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, true};
+  o2::framework::OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", "Event filtered;; Number of events", 6, 0., 6.)};
+
+  o2::framework::ConfigurableAxis pTAxis{"pTAxis", {200, -10.f, 10.f}, "p_{T} GeV/c"};
+  o2::framework::ConfigurableAxis nSigmaAxis{"nSigmaAxis", {200, -10.f, 10.f}, "p_{T} GeV/c"};
+  o2::framework::ConfigurableAxis DCAxyAxis{"DCAxyAxis", {1000, -0.2f, 0.2f}, "DCA xy (cm)"};
+  o2::framework::ConfigurableAxis DCAzAxis{"DCAzAxis", {1000, -0.2f, 0.2f}, "DCA z (cm)"};
+
+  o2::framework::Configurable<float> cutzvertex{"cutzvertex", 20.0f, "Accepted z-vertex range"}; // 20 cm
+  o2::framework::Configurable<bool> isTVX{"isTVX", true, "isTVX event selection"};
+  o2::framework::Configurable<bool> isNoTimeFrameBorder{"isNoTimeFrameBorder", true, "isNoTimeFrameBorder event selection"};
+  o2::framework::Configurable<bool> isNoITSROFrameBorder{"isNoITSROFrameBorder", true, "isNoITSROFrameBorder event selection"};
+  o2::framework::Configurable<float> cfgTPCNsigma{"cfgTPCNsigma", 4.0f, "TPC n sigma for deuteron PID"};
+  o2::framework::Configurable<float> cfgTOFNsigma{"cfgTOFNsigma", 4.0f, "TOF n sigma for deuteron PID"};
+  o2::framework::Configurable<float> cfgEta{"cfgEta", 0.8f, "Track eta selection"};
+  o2::framework::Configurable<int> cfgTPCNclsFound{"cfgTPCNclsFound", 100, "Minimum TPC clusters found"};
+  o2::framework::Configurable<float> cfgTPCChi2Ncl{"cfgTPCChi2Ncl", 4.0f, "Maximum TPC chi2 per N clusters"};
+  o2::framework::Configurable<float> cfgITSChi2Ncl{"cfgITSChi2Ncl", 36.0f, "Maximum ITS chi2 per N clusters"};
+  o2::framework::Configurable<int> cfgITScls{"cfgITScls", 2, "Minimum ITS clusters"};
+  o2::framework::Configurable<float> cfgMaxPt{"cfgMaxPt", 5.0f, "Maximum pT cut"};
+  o2::framework::Configurable<float> cfgMinPt{"cfgMinPt", 1.0f, "Maximum pT cut"};
+  o2::framework::Configurable<float> cfgDCAcut{"cfgDCAcut", 0.003f, "DCA cut for non prompt deuteron"};
+
+  o2::framework::Service<o2::ccdb::BasicCCDBManager> ccdb;
+
+  void init(o2::framework::InitContext&)
+  {
+    ccdb->setURL("http://alice-ccdb.cern.ch"); // Set CCDB URL to get magnetic field
+    ccdb->setCaching(true);
+    ccdb->setLocalObjectValidityChecking();
+    ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
+
+    o2::framework::AxisSpec ptAxis = {100, 0.0f, 10.0f, "#it{p}_{T} (GeV/#it{c})"};
+
+    // general QA histograms
+    QAHistos.add("hVtxZ", "Z-Vertex distribution after selection;Z (cm)", o2::framework::HistType::kTH1F, {{100, -50, 50}});
+    QAHistos.add("hDCAxyVsPt", "DCAxy #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAxyAxis}});
+    QAHistos.add("hDCAzVsPt", "DCAz #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAzAxis}});
+    QAHistos.add("hnSigmaTPCVsPt", "n#sigma TPC vs p_{T} for #bar{d} hypothesis; p_{T} (GeV/c); n#sigma TPC", {o2::framework::HistType::kTH2D, {pTAxis, nSigmaAxis}});
+    QAHistos.add("hnSigmaTOFVsPt", "n#sigma TOF vs p_{T} for #bar{d} hypothesis; p_{T} (GeV/c); n#sigma TOF", {o2::framework::HistType::kTH2D, {pTAxis, nSigmaAxis}});
+    QAHistos.add("ptAntiDeuteron", "ptAntiDeuteron", {o2::framework::HistType::kTH1F, {ptAxis}});
+    QAHistos.add("etaAntideuteron", "etaAntideuteron", {o2::framework::HistType::kTH1F, {{100, -1.0f, 1.0f, "eta #bar{d}"}}});
+    QAHistos.add("hDCAxyVsPt-pre_selection", "DCAxy #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAxyAxis}});
+    QAHistos.add("hDCAzVsPt-pre-selection", "DCAz #bar{d} vs p_{T}", {o2::framework::HistType::kTH2D, {pTAxis, DCAzAxis}});
+
+    // processed events
+    hProcessedEvents->GetXaxis()->SetBinLabel(1, "Events processed");
+    hProcessedEvents->GetXaxis()->SetBinLabel(2, "TVX");
+    hProcessedEvents->GetXaxis()->SetBinLabel(3, "TF border");
+    hProcessedEvents->GetXaxis()->SetBinLabel(4, "ITS-RO border");
+    hProcessedEvents->GetXaxis()->SetBinLabel(5, "z-vertex");
+    hProcessedEvents->GetXaxis()->SetBinLabel(6, o2::aod::filtering::H2fromLb::columnLabel());
+  }
+
+  // Tables
+  using CollisionCandidates = o2::soa::Join<o2::aod::Collisions, o2::aod::EvSels>;
+  using TrackCandidates = o2::soa::Join<o2::aod::Tracks, o2::aod::TracksCov, o2::aod::TracksExtra, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCFullDe, o2::aod::pidTOFFullDe>;
+
+  // Single-Track Selection
+  template <typename T1>
+  bool passedSingleTrackSelection(const T1& track)
+  {
+    // Single-Track Selections
+    if (std::abs(track.eta()) > cfgEta)
+      return false;
+    if (!track.hasITS())
+      return false;
+    if (!track.hasTPC())
+      return false;
+    if (!track.hasTOF())
+      return false;
+    if (track.tpcNClsFound() < cfgTPCNclsFound)
+      return false;
+    if (track.tpcChi2NCl() > cfgTPCChi2Ncl)
+      return false;
+    if (track.itsChi2NCl() > cfgITSChi2Ncl)
+      return false;
+    if (track.itsNCls() < cfgITScls)
+      return false;
+    if (track.pt() > cfgMaxPt)
+      return false;
+    if (track.pt() < cfgMinPt)
+      return false;
+    if (track.sign() > 0)
+      return false;
+
+    return true;
+  }
+
+  void fillTriggerTable(bool keepEvent[])
+  {
+    table(keepEvent[0]);
+  }
+
+  o2::framework::Preslice<o2::aod::TrackAssoc> trackIndicesPerCollision = o2::aod::track_association::collisionId;
+  int mCurrentRun = -1;
+
+  void process(CollisionCandidates const& collisions,
+               o2::aod::TrackAssoc const& trackIndices,
+               TrackCandidates const& tracks,
+               o2::aod::BCsWithTimestamps const&)
+  {
+    for (const auto& collision : collisions) // start loop over collisions
+    {
+      if (mCurrentRun != collision.bc_as<o2::aod::BCsWithTimestamps>().runNumber()) { // If the run is new then we need to initialize the propagator field
+        o2::parameters::GRPMagField* grpo = ccdb->getForTimeStamp<o2::parameters::GRPMagField>("GLO/Config/GRPMagField", collision.bc_as<o2::aod::BCsWithTimestamps>().timestamp());
+        o2::base::Propagator::initFieldFromGRP(grpo);
+        mCurrentRun = collision.bc_as<o2::aod::BCsWithTimestamps>().runNumber();
+      }
+
+      // Is event good? keepEvent[0] = non promp deuteron
+      bool keepEvent[1]{}; // explicitly zero-initialised
+      hProcessedEvents->Fill(0.5);
+      // Event selection cuts
+      if (isTVX && !collision.selection_bit(o2::aod::evsel::kIsTriggerTVX)) {
+        fillTriggerTable(keepEvent);
+        continue;
+      }
+      hProcessedEvents->Fill(1.5);
+      if (isNoTimeFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) {
+        fillTriggerTable(keepEvent);
+        continue;
+      }
+      hProcessedEvents->Fill(2.5);
+      if (isNoITSROFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
+        fillTriggerTable(keepEvent);
+        continue;
+      }
+      hProcessedEvents->Fill(3.5);
+      if (std::abs(collision.posZ()) > cutzvertex) {
+        fillTriggerTable(keepEvent);
+        continue;
+      }
+      hProcessedEvents->Fill(4.5);
+      QAHistos.fill(HIST("hVtxZ"), collision.posZ());
+
+      // Loop over tracks
+      const auto& trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex());
+
+      for (const auto& trackId : trackIdsThisCollision) { // start loop over tracks
+
+        const auto& track = tracks.rawIteratorAt(trackId.trackId());
+
+        std::array<float, 2> dca{track.dcaXY(), track.dcaZ()};
+        std::array<float, 3> pVec = track.pVector();
+
+        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);
+          getPxPyPz(trackPar, pVec);
+        }
+
+        if (!passedSingleTrackSelection(track)) {
+          continue;
+        }
+
+        const bool isTOFDe = std::abs(track.tofNSigmaDe()) < cfgTOFNsigma;
+        const bool isTPCDe = std::abs(track.tpcNSigmaDe()) < cfgTPCNsigma;
+
+        if (isTPCDe && isTOFDe) {
+
+          QAHistos.fill(HIST("hDCAxyVsPt-pre_selection"), track.pt(), dca[0]);
+          QAHistos.fill(HIST("hDCAzVsPt-pre-selection"), track.pt(), dca[1]);
+          if (std::abs(dca[0]) < cfgDCAcut) {
+            continue;
+          }
+          keepEvent[0] = true;
+          QAHistos.fill(HIST("ptAntiDeuteron"), track.pt());
+          QAHistos.fill(HIST("etaAntideuteron"), track.eta());
+          QAHistos.fill(HIST("hDCAxyVsPt"), track.pt(), dca[0]);
+          QAHistos.fill(HIST("hDCAzVsPt"), track.pt(), dca[1]);
+          QAHistos.fill(HIST("hnSigmaTPCVsPt"), track.pt(), track.tpcNSigmaDe());
+          QAHistos.fill(HIST("hnSigmaTOFVsPt"), track.pt(), track.tofNSigmaDe());
+        }
+
+      } // end track loop
+      if (keepEvent[0]) {
+        hProcessedEvents->Fill(5.5);
+      }
+
+      fillTriggerTable(keepEvent);
+    }
+  }
+};
+
+o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)
+{
+  return o2::framework::WorkflowSpec{o2::framework::adaptAnalysisTask<H2fromLbFilter>(cfgc)};
+}

EventFiltering/filterTables.h

@@ -210,6 +210,9 @@ DECLARE_SOA_COLUMN(EtaPrimeP, hasEtaPrimeP, bool);   //! eta'P meson candidate (
 DECLARE_SOA_COLUMN(EtaPrimePP, hasEtaPrimePP, bool); //! eta'PP meson candidate (3pi) in the collision
 DECLARE_SOA_COLUMN(EtaPrimed, hasEtaPrimed, bool);   //! eta'd meson candidate (3pi) in the collision
 
+// deuterons non prompt
+DECLARE_SOA_COLUMN(H2fromLb, hasH2fromLb, bool); //! deuteron from Lb decay
+
 } // namespace filtering
 
 namespace decision
@@ -362,6 +365,11 @@ DECLARE_SOA_TABLE(HeavyNeutralMesonFilters, "AOD", "HNMesonFilters", //!
 
 using HeavyNeutralMesonFilter = HeavyNeutralMesonFilters::iterator;
 
+// beauty non prompt deuterons
+DECLARE_SOA_TABLE(H2fromLbFilters, "AOD", "H2LbFilters", //!
+                  filtering::H2fromLb);
+using H2fromLbFilter = H2fromLbFilters::iterator;
+
 // cefp decision
 DECLARE_SOA_TABLE(CefpDecisions, "AOD", "CefpDecision", //!
                   decision::BCId, decision::GlobalBCId, decision::EvSelBC, decision::CollisionTime, decision::CollisionTimeRes, decision::CefpTriggered0, decision::CefpTriggered1, decision::CefpSelected0, decision::CefpSelected1);
@@ -373,11 +381,11 @@ DECLARE_SOA_TABLE(BCRanges, "AOD", "BCRanges", //!
 using BCRange = BCRanges::iterator;
 
 /// List of the available filters, the description of their tables and the name of the tasks
-constexpr int NumberOfFilters{15};
-constexpr std::array<char[32], NumberOfFilters> AvailableFilters{"NucleiFilters", "DiffractionFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "StrangenessFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "DoublePhiFilters", "HeavyNeutralMesonFilters", "GlobalDimuonFilters"};
-constexpr std::array<char[16], NumberOfFilters> FilterDescriptions{"NucleiFilters", "DiffFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "LFStrgFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "LF2PhiFilters", "HNMesonFilters", "GlDimuonFilters"};
-constexpr std::array<char[128], NumberOfFilters> FilteringTaskNames{"o2-analysis-nuclei-filter", "o2-analysis-diffraction-filter", "o2-analysis-dq-filter-pp-with-association", "o2-analysis-hf-filter", "o2-analysis-cf-filter", "o2-analysis-je-filter", "o2-analysis-je-hf-filter", "o2-analysis-fje-filter", "o2-analysis-lf-strangeness-filter", "o2-analysis-mult-filter", "o2-analysis-em-photon-filter", "o2-analysis-lf-f1proton-filter", "o2-analysis-lf-doublephi-filter", "o2-analysis-heavy-neutral-meson-filter", "o2-analysis-em-global-dimuon-filter"};
-constexpr o2::framework::pack<NucleiFilters, DiffractionFilters, DqFilters, HfFilters, CFFilters, JetFilters, JetHFFilters, FullJetFilters, StrangenessFilters, MultFilters, PhotonFilters, F1ProtonFilters, DoublePhiFilters, HeavyNeutralMesonFilters, GlobalDimuonFilters> FiltersPack;
+constexpr int NumberOfFilters{16};
+constexpr std::array<char[32], NumberOfFilters> AvailableFilters{"NucleiFilters", "DiffractionFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "StrangenessFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "DoublePhiFilters", "HeavyNeutralMesonFilters", "GlobalDimuonFilters", "H2fromLbFilters"};
+constexpr std::array<char[16], NumberOfFilters> FilterDescriptions{"NucleiFilters", "DiffFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "LFStrgFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "LF2PhiFilters", "HNMesonFilters", "GlDimuonFilters", "H2LbFilters"};
+constexpr std::array<char[128], NumberOfFilters> FilteringTaskNames{"o2-analysis-nuclei-filter", "o2-analysis-diffraction-filter", "o2-analysis-dq-filter-pp-with-association", "o2-analysis-hf-filter", "o2-analysis-cf-filter", "o2-analysis-je-filter", "o2-analysis-je-hf-filter", "o2-analysis-fje-filter", "o2-analysis-lf-strangeness-filter", "o2-analysis-mult-filter", "o2-analysis-em-photon-filter", "o2-analysis-lf-f1proton-filter", "o2-analysis-lf-doublephi-filter", "o2-analysis-heavy-neutral-meson-filter", "o2-analysis-em-global-dimuon-filter", "o2-analysis-h2-from-beauty-filter"};
+constexpr o2::framework::pack<NucleiFilters, DiffractionFilters, DqFilters, HfFilters, CFFilters, JetFilters, JetHFFilters, FullJetFilters, StrangenessFilters, MultFilters, PhotonFilters, F1ProtonFilters, DoublePhiFilters, HeavyNeutralMesonFilters, GlobalDimuonFilters, H2fromLbFilters> FiltersPack;
 static_assert(o2::framework::pack_size(FiltersPack) == NumberOfFilters);
 
 template <typename T, typename C>