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

EventFiltering/PWGLF/nucleiFilter.cxx

@@ -1,3 +1,3 @@
 // 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.
@@ -8,11 +8,13 @@
 // 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.
 // O2 includes
 
 #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 @@
 {
 
 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 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,17 +101,18 @@
   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;
   float d_bz;
   Service<o2::ccdb::BasicCCDBManager> ccdb;
   using TrackCandidates = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TrackSelection, aod::TracksDCA, aod::EvTimeTOFFT0ForTrack, aod::pidTPCFullPi, aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTPCFullTr, aod::pidTPCFullHe, aod::pidTPCFullAl, aod::pidTOFFullDe, aod::pidTOFFullTr, aod::pidTOFFullHe, aod::pidTOFFullAl>; // FIXME: positio has been changed
   o2::aod::pidtofgeneric::TofPidNewCollision<TrackCandidates::iterator> bachelorTOFPID;
@@ -160,9 +162,9 @@
   } 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.};
 
@@ -211,7 +213,7 @@
     }
 
     // In case override, don't proceed, please - no CCDB access required
     if (trgH3L3Body.d_bz_input > -990) {
       d_bz = trgH3L3Body.d_bz_input;
       fitter3body.setBz(d_bz);
       o2::parameters::GRPMagField grpmag;
@@ -246,7 +248,7 @@
     // Set magnetic field value once known
     fitter3body.setBz(d_bz);
 
     if (trgH3L3Body.useMatCorrType == 2) {
       // setMatLUT only after magfield has been initalized
       // (setMatLUT has implicit and problematic init field call if not)
       o2::base::Propagator::Instance()->setMatLUT(lut);
@@ -313,59 +315,80 @@
     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()))};
         for (int iC{0}; iC < 7; ++iC) {
           avgClsSize += track.itsClsSizeInLayer(iC);
         }
         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) {
         continue;
       }
 
       if (std::abs(track.tpcNSigmaDe()) < 5) {
         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 @@
         }
 
         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 @@
         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
       //
@@ -474,7 +544,7 @@
 
       std::array<float, 3> pos = {0.};
       const auto& vtxXYZ = fitter3body.getPCACandidate();
       for (int i = 0; i < 3; i++) {
         pos[i] = vtxXYZ[i];
       }
 
@@ -531,7 +601,7 @@
           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 @@
           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]);
   }
 };
 

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;