[PWGLF,Trigger] Add V0 with He3, triton femto, and tracked 3 body triggers (#10834)

mpuccio · web-flow · commit 2989f9046f8f · 2025-04-14T21:26:19.000+02:00
diff --git a/EventFiltering/PWGLF/nucleiFilter.cxx b/EventFiltering/PWGLF/nucleiFilter.cxx
@@ -13,6 +13,8 @@
 #include <cmath>
 #include <string>
 
+#include "Math/Vector4D.h"
+#include "Math/GenVector/Boost.h"
 #include "Common/DataModel/EventSelection.h"
 #include "Common/DataModel/PIDResponse.h"
 #include "Common/DataModel/TrackSelectionTables.h"
@@ -37,6 +39,7 @@
 #include "CCDB/BasicCCDBManager.h"
 #include "DCAFitter/DCAFitterN.h"
 #include "PWGLF/DataModel/pidTOFGeneric.h"
+#include "PWGLF/DataModel/LFStrangenessTables.h"
 #include "Common/Core/PID/PIDTOF.h"
 
 using namespace o2;
@@ -47,8 +50,6 @@ namespace
 {
 
 static constexpr int nNuclei{3};
-static constexpr int nHyperNuclei{1};
-static constexpr int nITStriggers{2};
 static constexpr int nCutsPID{5};
 static constexpr std::array<float, nNuclei> masses{
   constants::physics::MassDeuteron, constants::physics::MassTriton,
@@ -57,7 +58,7 @@ static constexpr std::array<int, nNuclei> charges{1, 1, 2};
 static const std::vector<std::string> matterOrNot{"Matter", "Antimatter"};
 static const std::vector<std::string> nucleiNames{"H2", "H3", "Helium"};
 static const std::vector<std::string> hypernucleiNames{"H3L"}; // 3-body decay case
-static const std::vector<std::string> columnsNames{o2::aod::filtering::H2::columnLabel(), "fH3", o2::aod::filtering::He::columnLabel(), o2::aod::filtering::H3L3Body::columnLabel(), o2::aod::filtering::ITSmildIonisation::columnLabel(), o2::aod::filtering::ITSextremeIonisation::columnLabel()};
+static const std::vector<std::string> columnsNames{o2::aod::filtering::H2::columnLabel(), o2::aod::filtering::He::columnLabel(), o2::aod::filtering::HeV0::columnLabel(), o2::aod::filtering::TritonFemto::columnLabel(), o2::aod::filtering::H3L3Body::columnLabel(), o2::aod::filtering::Tracked3Body::columnLabel(), o2::aod::filtering::ITSmildIonisation::columnLabel(), o2::aod::filtering::ITSextremeIonisation::columnLabel()};
 static const std::vector<std::string> cutsNames{
   "TPCnSigmaMin", "TPCnSigmaMax", "TOFnSigmaMin", "TOFnSigmaMax", "TOFpidStartPt"};
 constexpr double betheBlochDefault[nNuclei][6]{
@@ -100,13 +101,14 @@ struct nucleiFilter {
   Configurable<float> cfgCutNclusTPC{"cfgCutNclusTPC", 80, "Minimum number of TPC clusters"};
   Configurable<float> cfgCutDCAxy{"cfgCutDCAxy", 3, "Max DCAxy"};
   Configurable<float> cfgCutDCAz{"cfgCutDCAz", 10, "Max DCAz"};
+  Configurable<float> cfgCutKstar{"cfgCutKstar", 1.f, "Kstar cut for triton femto trigger"};
 
   Configurable<LabeledArray<double>> cfgBetheBlochParams{"cfgBetheBlochParams", {betheBlochDefault[0], nNuclei, 6, nucleiNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"};
   Configurable<LabeledArray<double>> cfgMomentumScalingBetheBloch{"cfgMomentumScalingBetheBloch", {bbMomScalingDefault[0], nNuclei, 2, nucleiNames, matterOrNot}, "TPC Bethe-Bloch momentum scaling for light nuclei"};
   Configurable<LabeledArray<double>> cfgMinTPCmom{"cfgMinTPCmom", {minTPCmom[0], nNuclei, 2, nucleiNames, matterOrNot}, "Minimum TPC p/Z for nuclei PID"};
 
   Configurable<LabeledArray<float>> cfgCutsPID{"nucleiCutsPID", {cutsPID[0], nNuclei, nCutsPID, nucleiNames, cutsNames}, "Nuclei PID selections"};
-  Configurable<bool> fixTPCinnerParam{"fixTPCinnerParam", false, "Fix TPC inner param"};
+  Configurable<bool> cfgFixTPCinnerParam{"cfgFixTPCinnerParam", false, "Fix TPC inner param"};
 
   // variable/tool for hypertriton 3body decay
   int mRunNumber;
@@ -160,9 +162,9 @@ struct nucleiFilter {
   } trgH3L3Body;
 
   HistogramRegistry qaHists{"qaHists", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
-  OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", ";;Number of filtered events", nNuclei + nHyperNuclei + nITStriggers + 1, -0.5, nNuclei + nHyperNuclei + nITStriggers + 0.5)};
+  OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", ";;Number of filtered events", kNtriggers + 1, -0.5, static_cast<double>(kNtriggers) + 0.5)};
 
-  void init(o2::framework::InitContext&)
+  void init(InitContext&)
   {
     std::vector<double> ptBinning = {0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4., 5.};
 
@@ -313,27 +315,50 @@ struct nucleiFilter {
     bachelorTOFPID.SetParams(mRespParamsV2);
   }
 
+  enum {
+    kH2 = 0,
+    kHe,
+    kHeV0,
+    kTritonFemto,
+    kH3L3Body,
+    kTracked3Body,
+    kITSmildIonisation,
+    kITSextremeIonisation,
+    kNtriggers
+  } TriggerType;
   // void process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision, aod::Vtx3BodyDatas const& vtx3bodydatas, TrackCandidates const& tracks)
   using ColWithEvTime = soa::Join<aod::Collisions, aod::EvSels, aod::EvTimeTOFFT0>;
-  void process(ColWithEvTime::iterator const& collision, aod::Decay3Bodys const& decay3bodys, TrackCandidates const& tracks, aod::BCsWithTimestamps const&)
+  void process(ColWithEvTime::iterator const& collision, aod::Decay3Bodys const& decay3bodys, TrackCandidates const& tracks, aod::AssignedTracked3Bodys const& tracked3Bodys, aod::V0s const& v0s, aod::BCsWithTimestamps const&)
   {
     // collision process loop
-    bool keepEvent[nNuclei + nHyperNuclei + nITStriggers]{false};
+    std::array<bool, kNtriggers> keepEvent{false};
     //
     qaHists.fill(HIST("fCollZpos"), collision.posZ());
     hProcessedEvents->Fill(0);
     //
     if (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
-      tags(keepEvent[0], keepEvent[2], keepEvent[3], keepEvent[4], keepEvent[5]);
+      tags(keepEvent[kH2], keepEvent[kHe], keepEvent[kHeV0], keepEvent[kTritonFemto], keepEvent[kH3L3Body], keepEvent[kTracked3Body], keepEvent[kITSmildIonisation], keepEvent[kITSextremeIonisation]);
       return;
     }
+
     //
     const double bgScalings[nNuclei][2]{
       {charges[0] * cfgMomentumScalingBetheBloch->get(0u, 0u) / masses[0], charges[0] * cfgMomentumScalingBetheBloch->get(0u, 1u) / masses[0]},
       {charges[1] * cfgMomentumScalingBetheBloch->get(1u, 0u) / masses[1], charges[1] * cfgMomentumScalingBetheBloch->get(1u, 1u) / masses[1]},
       {charges[2] * cfgMomentumScalingBetheBloch->get(2u, 0u) / masses[2], charges[2] * cfgMomentumScalingBetheBloch->get(2u, 1u) / masses[2]}};
 
-    for (auto& track : tracks) { // start loop over tracks
+    constexpr int nucleusIndex[nNuclei]{kH2, -1, kHe}; /// remap for nuclei triggers
+    std::vector<int> h3indices;
+    std::vector<ROOT::Math::PtEtaPhiMVector> h3vectors;
+
+    auto getNsigma = [&](const auto& track, int iN, int iC) {
+      float fixTPCrigidity{(cfgFixTPCinnerParam && (track.pidForTracking() == track::PID::Helium3 || track.pidForTracking() == track::PID::Alpha)) ? 0.5f : 1.f};
+      double expBethe{tpc::BetheBlochAleph(static_cast<double>(track.tpcInnerParam() * fixTPCrigidity * bgScalings[iN][iC]), cfgBetheBlochParams->get(iN, 0u), cfgBetheBlochParams->get(iN, 1u), cfgBetheBlochParams->get(iN, 2u), cfgBetheBlochParams->get(iN, 3u), cfgBetheBlochParams->get(iN, 4u))};
+      double expSigma{expBethe * cfgBetheBlochParams->get(iN, 5u)};
+      return static_cast<float>((track.tpcSignal() - expBethe) / expSigma);
+    };
+
+    for (const auto& track : tracks) { // start loop over tracks
       if (track.itsNCls() >= cfgCutNclusExtremeIonisationITS) {
         double avgClsSize{0.};
         double cosL{std::sqrt(1. / (1. + track.tgl() * track.tgl()))};
@@ -342,8 +367,8 @@ struct nucleiFilter {
         }
         avgClsSize = avgClsSize * cosL / track.itsNCls();
         qaHists.fill(HIST("fExtremeIonisationITS"), track.itsNCls(), avgClsSize, track.p());
-        keepEvent[4] = track.p() > cfgMomentumCutExtremeIonisation && avgClsSize > cfgCutClsSizeMildIonisation;
-        keepEvent[5] = track.p() > cfgMomentumCutExtremeIonisation && avgClsSize > cfgCutClsSizeExtremeIonisation;
+        keepEvent[kITSmildIonisation] = track.p() > cfgMomentumCutExtremeIonisation && avgClsSize > cfgCutClsSizeMildIonisation;
+        keepEvent[kITSextremeIonisation] = track.p() > cfgMomentumCutExtremeIonisation && avgClsSize > cfgCutClsSizeExtremeIonisation;
       }
       if (track.itsNCls() < cfgCutNclusITS ||
           track.tpcNClsFound() < cfgCutNclusTPC) {
@@ -354,18 +379,16 @@ struct nucleiFilter {
         qaHists.fill(HIST("fDeuTOFNsigma"), track.p() * track.sign(), track.tofNSigmaDe());
       }
 
-      if (track.sign() > 0 && (std::abs(track.dcaXY()) > cfgCutDCAxy ||
-                               std::abs(track.dcaZ()) > cfgCutDCAz)) {
-        continue;
-      }
+      bool passesDCAselection{(track.sign() < 0 || (std::abs(track.dcaXY()) < cfgCutDCAxy &&
+                                                    std::abs(track.dcaZ()) < cfgCutDCAz))};
 
       float nSigmaTPC[nNuclei]{
         track.tpcNSigmaDe(), track.tpcNSigmaTr(), track.tpcNSigmaHe()};
       const float nSigmaTOF[nNuclei]{
         track.tofNSigmaDe(), track.tofNSigmaTr(), track.tofNSigmaHe()};
       const int iC{track.sign() < 0};
 
-      float fixTPCrigidity{(fixTPCinnerParam && (track.pidForTracking() == track::PID::Helium3 || track.pidForTracking() == track::PID::Alpha)) ? 0.5f : 1.f};
+      float fixTPCrigidity{(cfgFixTPCinnerParam && (track.pidForTracking() == track::PID::Helium3 || track.pidForTracking() == track::PID::Alpha)) ? 0.5f : 1.f};
 
       // fill QA hist: dEdx for all charged tracks
       qaHists.fill(HIST("fTPCsignalAll"), track.sign() * track.tpcInnerParam() * fixTPCrigidity, track.tpcSignal());
@@ -377,9 +400,7 @@ struct nucleiFilter {
         }
 
         if (cfgBetheBlochParams->get(iN, 5u) > 0.f) {
-          double expBethe{tpc::BetheBlochAleph(static_cast<double>(track.tpcInnerParam() * fixTPCrigidity * bgScalings[iN][iC]), cfgBetheBlochParams->get(iN, 0u), cfgBetheBlochParams->get(iN, 1u), cfgBetheBlochParams->get(iN, 2u), cfgBetheBlochParams->get(iN, 3u), cfgBetheBlochParams->get(iN, 4u))};
-          double expSigma{expBethe * cfgBetheBlochParams->get(iN, 5u)};
-          nSigmaTPC[iN] = static_cast<float>((track.tpcSignal() - expBethe) / expSigma);
+          nSigmaTPC[iN] = getNsigma(track, iN, iC);
         }
         h2TPCnSigma[iN]->Fill(track.sign() * track.tpcInnerParam() * fixTPCrigidity, nSigmaTPC[iN]);
         if (nSigmaTPC[iN] < cfgCutsPID->get(iN, 0u) || nSigmaTPC[iN] > cfgCutsPID->get(iN, 1u)) {
@@ -388,12 +409,61 @@ struct nucleiFilter {
         if (track.p() > cfgCutsPID->get(iN, 4u) && (nSigmaTOF[iN] < cfgCutsPID->get(iN, 2u) || nSigmaTOF[iN] > cfgCutsPID->get(iN, 3u))) {
           continue;
         }
-        keepEvent[iN] = true;
-        if (keepEvent[iN]) {
+        if (iN == 1 && passesDCAselection) {
+          h3indices.push_back(track.globalIndex());
+          h3vectors.emplace_back(track.pt(), track.eta(), track.phi(), masses[iN]);
+        }
+        if (nucleusIndex[iN] < 0) {
+          continue;
+        }
+        keepEvent[nucleusIndex[iN]] = passesDCAselection;
+        if (keepEvent[nucleusIndex[iN]]) {
           h2TPCsignal[iN]->Fill(track.sign() * track.tpcInnerParam() * fixTPCrigidity, track.tpcSignal());
         }
       }
 
+      for (const auto& track : tracks) {
+        if (track.itsNCls() < cfgCutNclusITS ||
+            track.tpcNClsFound() < cfgCutNclusTPC ||
+            std::abs(track.dcaXY()) > cfgCutDCAxy ||
+            std::abs(track.dcaZ()) > cfgCutDCAz ||
+            std::abs(track.eta()) > cfgCutEta) {
+          continue;
+        }
+        const ROOT::Math::PtEtaPhiMVector trackVector(track.pt(), track.eta(), track.phi(), constants::physics::MassPiMinus);
+        for (size_t iH3{0}; iH3 < h3vectors.size(); ++iH3) {
+          if (h3indices[iH3] == track.globalIndex()) {
+            continue;
+          }
+          const auto& h3vector = h3vectors[iH3];
+          auto pivector = trackVector;
+          auto cm = h3vector + trackVector;
+          const ROOT::Math::Boost boost(cm.BoostToCM());
+          boost(pivector);
+          if (pivector.P() < cfgCutKstar) {
+            keepEvent[kTritonFemto] = true;
+            break;
+          }
+        }
+      }
+
+      for (const auto& v0 : v0s) {
+        const auto& posTrack = tracks.rawIteratorAt(v0.posTrackId());
+        const auto& negTrack = tracks.rawIteratorAt(v0.negTrackId());
+        if ((posTrack.itsNCls() < cfgCutNclusITS || posTrack.tpcNClsFound() < cfgCutNclusTPC) &&
+            (negTrack.itsNCls() < cfgCutNclusITS || negTrack.tpcNClsFound() < cfgCutNclusTPC)) {
+          continue;
+        }
+        float nSigmas[2]{
+          cfgBetheBlochParams->get(2, 5u) > 0.f ? getNsigma(posTrack, 2, 0) : posTrack.tpcNSigmaHe(),
+          cfgBetheBlochParams->get(2, 5u) > 0.f ? getNsigma(negTrack, 2, 1) : negTrack.tpcNSigmaHe()};
+        if ((nSigmas[0] > cfgCutsPID->get(2, 0u) && nSigmas[0] < cfgCutsPID->get(2, 1u)) ||
+            (nSigmas[1] > cfgCutsPID->get(2, 0u) && nSigmas[1] < cfgCutsPID->get(2, 1u))) {
+          keepEvent[kHeV0] = true;
+          break;
+        }
+      }
+
       //
       // fill QA histograms
       //
@@ -531,7 +601,7 @@ struct nucleiFilter {
           qaHists.fill(HIST("fBachDeuTOFNsigma"), track2.p() * track2.sign(), tofNSigmaDeuteron);
           qaHists.fill(HIST("fH3LDcaVsPt"), pt3B, dcaDaughters);
           qaHists.fill(HIST("fH3LCosPAVsPt"), pt3B, vtxCosPA);
-          keepEvent[3] = true;
+          keepEvent[kH3L3Body] = true;
         }
       }
       if (invmassAntiH3L >= trgH3L3Body.h3LMassLowerlimit && invmassAntiH3L <= trgH3L3Body.h3LMassUpperlimit) {
@@ -541,17 +611,20 @@ struct nucleiFilter {
           qaHists.fill(HIST("fBachDeuTOFNsigma"), track2.p() * track2.sign(), tofNSigmaDeuteron);
           qaHists.fill(HIST("fH3LDcaVsPt"), pt3B, dcaDaughters);
           qaHists.fill(HIST("fH3LCosPAVsPt"), pt3B, vtxCosPA);
-          keepEvent[3] = true;
+          keepEvent[kH3L3Body] = true;
         }
       }
     }
 
-    for (int iDecision{0}; iDecision < nNuclei + nHyperNuclei + nITStriggers; ++iDecision) {
+    keepEvent[kTracked3Body] = tracked3Bodys.size() > 0;
+
+    for (int iDecision{0}; iDecision < kNtriggers; ++iDecision) {
       if (keepEvent[iDecision]) {
         hProcessedEvents->Fill(iDecision + 1);
       }
     }
-    tags(keepEvent[0], keepEvent[2], keepEvent[3], keepEvent[4], keepEvent[5]);
+
+    tags(keepEvent[kH2], keepEvent[kHe], keepEvent[kHeV0], keepEvent[kTritonFemto], keepEvent[kH3L3Body], keepEvent[kTracked3Body], keepEvent[kITSmildIonisation], keepEvent[kITSextremeIonisation]);
   }
 };
 
diff --git a/EventFiltering/filterTables.h b/EventFiltering/filterTables.h
@@ -48,6 +48,8 @@ namespace filtering
 {
 DECLARE_SOA_COLUMN(H2, hasH2, bool);                                     //! deuteron trigger for the helium normalisation (to be downscaled)
 DECLARE_SOA_COLUMN(He, hasHe, bool);                                     //! helium
+DECLARE_SOA_COLUMN(HeV0, hasHeV0, bool);                                 //! V0 containing a V0
+DECLARE_SOA_COLUMN(TritonFemto, hasTritonFemto, bool);                   //! Triton hadron femtoscopy
 DECLARE_SOA_COLUMN(H3L3Body, hasH3L3Body, bool);                         //! hypertriton 3body
 DECLARE_SOA_COLUMN(ITSextremeIonisation, hasITSextremeIonisation, bool); //! ITS extreme ionisation
 DECLARE_SOA_COLUMN(ITSmildIonisation, hasITSmildIonisation, bool);       //! ITS mild ionisation (normalisation of the extreme ionisation), to be downscaled
@@ -215,7 +217,7 @@ DECLARE_SOA_COLUMN(BCend, hasBCend, uint64_t);     //! CEFP bcrange
 
 // nuclei
 DECLARE_SOA_TABLE(NucleiFilters, "AOD", "NucleiFilters", //!
-                  filtering::H2, filtering::He, filtering::H3L3Body, filtering::ITSmildIonisation,
+                  filtering::H2, filtering::He, filtering::HeV0, filtering::TritonFemto, filtering::H3L3Body, filtering::Tracked3Body, filtering::ITSmildIonisation,
                   filtering::ITSextremeIonisation);
 using NucleiFilter = NucleiFilters::iterator;
 
