[PWGLF,Trigger] [Event Filtering] adding trigger condition for lambdalambda

EventFiltering/PWGLF/strangenessFilter.cxx

@@ -15,6 +15,7 @@
 /// \since June 1, 2021
 
 #include <cmath>
+#include "TVector3.h"
 #include "CCDB/BasicCCDBManager.h"
 #include "DataFormatsParameters/GRPMagField.h"
 #include "DataFormatsParameters/GRPObject.h"
@@ -25,6 +26,7 @@
 #include "Framework/AnalysisDataModel.h"
 #include "Framework/ASoAHelpers.h"
 #include "ReconstructionDataFormats/Track.h"
+#include "ReconstructionDataFormats/TrackParametrization.h"
 #include "Common/Core/RecoDecay.h"
 #include "Common/Core/trackUtilities.h"
 #include "PWGLF/DataModel/LFStrangenessTables.h"
@@ -73,10 +75,38 @@
   HistogramRegistry QAHistosTriggerParticles{"QAHistosTriggerParticles", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
   HistogramRegistry QAHistosStrangenessTracking{"QAHistosStrangenessTracking", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
   HistogramRegistry EventsvsMultiplicity{"EventsvsMultiplicity", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
-  OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", "Strangeness - event filtered;; Number of events", 16, -1., 15.)};
+  OutputObj<TH1D> hProcessedEvents{TH1D("hProcessedEvents", "Strangeness - event filtered;; Number of events", 17, -1., 16.)};
   OutputObj<TH1F> hCandidate{TH1F("hCandidate", "; Candidate pass selection; Number of events", 30, 0., 30.)};
   OutputObj<TH1F> hEvtvshMinPt{TH1F("hEvtvshMinPt", " Number of h-Omega events with pT_h higher than thrd; min p_{T, trigg} (GeV/c); Number of events", 11, 0., 11.)};
 
+  // Dedicated selection criteria for lambda-lambda
+  struct : ConfigurableGroup {
+    Configurable<float> cfgv0radiusMin{"cfgv0radiusMin", 1.2, "minimum decay radius"};
+    Configurable<float> cfgDCAPosToPVMin{"cfgDCAPosToPVMin", 0.05, "minimum DCA to PV for positive track"};
+    Configurable<float> cfgDCANegToPVMin{"cfgDCANegToPVMin", 0.2, "minimum DCA to PV for negative track"};
+    Configurable<float> cfgv0CosPA{"cfgv0CosPA", 0.995, "minimum v0 cosine"};
+    Configurable<float> cfgDCAV0Dau{"cfgDCAV0Dau", 1.0, "maximum DCA between daughters"};
+    Configurable<float> cfgV0PtMin{"cfgV0PtMin", 0, "minimum pT for lambda"};
+    Configurable<float> cfgV0RapMin{"cfgV0RapMin", -0.5, "maximum rapidity"};
+    Configurable<float> cfgV0RapMax{"cfgV0RapMax", 0.5, "maximum rapidity"};
+    Configurable<float> cfgV0LifeTime{"cfgV0LifeTime", 30., "maximum lambda lifetime"};
+    Configurable<int16_t> cfgDaughTPCnclsMin{"cfgDaughTPCnclsMin", 70, "minimum fired crossed rows"};
+    Configurable<uint8_t> cfgITSNclus{"cfgITSNclus", 1, "minimum its cluster"};
+    Configurable<float> cfgRCrossedFindable{"cfgRCrossedFindable", 0.0, "minimum ratio of crossed rows over findable clusters"};
+    Configurable<float> cfgDaughPIDCutsTPCPr{"cfgDaughPIDCutsTPCPr", 5, "proton nsigma for TPC"};
+    Configurable<float> cfgDaughPIDCutsTPCPi{"cfgDaughPIDCutsTPCPi", 5, "pion nsigma for TPC"};
+    Configurable<float> cfgDaughEtaMin{"cfgDaughEtaMin", -0.8, "minimum daughter eta"};
+    Configurable<float> cfgDaughEtaMax{"cfgDaughEtaMax", 0.8, "maximum daughter eta"};
+    Configurable<float> cfgDaughPrPt{"cfgDaughPrPt", 0.5, "minimum daughter proton pt"};
+    Configurable<float> cfgDaughPiPt{"cfgDaughPiPt", 0.5, "minimum daughter pion pt"};
+    Configurable<float> cfgLambdaMassWindow{"cfgLambdaMassWindow", 0.01, "window for lambda mass selection"};
+    Configurable<float> cfgCompV0Rej{"cfgCompV0Rej", 0.01, "competing V0 rejection"};
+    Configurable<float> cfgMinCPAV0V0{"cfgMinCPAV0V0", 0.8, "minimum CPA of v0v0"};
+    Configurable<float> cfgMaxRadiusV0V0{"cfgMaxRadiusV0V0", 10.0, "maximum radius of v0v0"};
+    Configurable<float> cfgMaxDistanceV0V0{"cfgMaxDistanceV0V0", 5.0, "maximum distance of v0v0"};
+    Configurable<float> cfgMaxDCAV0V0{"cfgMaxDCAV0V0", 5.0, "maximum DCA of v0v0"};
+  } cfgLLCuts;
+
   // Selection criteria for cascades
   Configurable<bool> useCascadeMomentumAtPrimVtx{"useCascadeMomentumAtPrimVtx", false, "use cascade momentum at PV"};
   Configurable<bool> doextraQA{"doextraQA", 1, "do extra QA"};
@@ -178,6 +208,81 @@
     return track.pt() > hMinPt && std::abs(track.eta()) < hEta && track.tpcNClsCrossedRows() >= tpcmincrossedrows && track.tpcCrossedRowsOverFindableCls() >= 0.8f && track.tpcChi2NCl() <= 4.f && track.itsChi2NCl() <= 36.f && (track.itsClusterMap() & 0x7) != 0;
   }
 
+  float getV0V0DCA(TVector3 v01pos, TVector3 v01mom, TVector3 v02pos, TVector3 v02mom)
+  {
+    TVector3 posdiff = v02pos - v01pos;
+    TVector3 cross = v01mom.Cross(v02mom);
+    TVector3 dcaVec = (posdiff.Dot(cross) / cross.Mag2()) * cross;
+    return dcaVec.Mag();
+  }
+  float getV0V0CPA(TVector3 v01mom, TVector3 v02mom)
+  {
+    return v01mom.Dot(v02mom) / (v01mom.Mag() * v02mom.Mag());
+  }
+  float getV0V0Distance(TVector3 v01pos, TVector3 v02pos)
+  {
+    TVector3 posdiff = v02pos - v01pos;
+    return posdiff.Mag();
+  }
+  float getV0V0Radius(TVector3 v01pos, TVector3 v01mom, TVector3 v02pos, TVector3 v02mom)
+  {
+    TVector3 posdiff = v02pos - v01pos;
+    v01mom *= 1. / v01mom.Mag();
+    v02mom *= 1. / v02mom.Mag();
+    float dd = 1. - TMath::Power(v01mom.Dot(v02mom), 2);
+    if (dd < 1e-5)
+      return 999;
+    float tt = posdiff.Dot(v01mom - v01mom.Dot(v02mom) * v02mom) / dd;
+    float ss = -posdiff.Dot(v02mom - v01mom.Dot(v02mom) * v01mom) / dd;
+    TVector3 radVec = v01pos + v02pos + tt * v01mom + ss * v02mom;
+    radVec *= 0.5;
+    return radVec.Mag();
+  }
+  bool isSelectedV0V0(TVector3 v01pos, TVector3 v01mom, TVector3 v02pos, TVector3 v02mom)
+  {
+    if (getV0V0DCA(v01pos, v01mom, v02pos, v02mom) > cfgLLCuts.cfgMaxDCAV0V0)
+      return false;
+    if (getV0V0CPA(v01mom, v02mom) < cfgLLCuts.cfgMinCPAV0V0)
+      return false;
+    if (getV0V0Distance(v01pos, v02pos) > cfgLLCuts.cfgMaxDistanceV0V0)
+      return false;
+    if (getV0V0Radius(v01pos, v01mom, v02pos, v02mom) > cfgLLCuts.cfgMaxRadiusV0V0)
+      return false;
+
+    return true;
+  }
+
+  template <typename T>
+  bool isSelectedV0Daughter(T const& track)
+  {
+    if (track.tpcNClsCrossedRows() < cfgLLCuts.cfgDaughTPCnclsMin)
+      return false;
+    if (track.tpcCrossedRowsOverFindableCls() < cfgLLCuts.cfgRCrossedFindable)
+      return false;
+    if (track.itsNCls() < cfgLLCuts.cfgITSNclus)
+      return false;
+    if (track.eta() > cfgLLCuts.cfgDaughEtaMax)
+      return false;
+    if (track.eta() < cfgLLCuts.cfgDaughEtaMin)
+      return false;
+
+    return true;
+  }
+  template <typename T>
+  bool isSelectedV0DaughterPID(T const& track, int pid) // pid 0: proton, pid 1: pion
+  {
+    if (pid == 0 && std::abs(track.tpcNSigmaPr()) > cfgLLCuts.cfgDaughPIDCutsTPCPr)
+      return false;
+    if (pid == 1 && std::abs(track.tpcNSigmaPi()) > cfgLLCuts.cfgDaughPIDCutsTPCPi)
+      return false;
+    if (pid == 0 && track.pt() < cfgLLCuts.cfgDaughPrPt)
+      return false;
+    if (pid == 1 && track.pt() < cfgLLCuts.cfgDaughPiPt)
+      return false;
+
+    return true;
+  }
+
   void init(o2::framework::InitContext&)
   {
     // set V0 parameters in the helper
@@ -223,6 +328,7 @@
     hProcessedEvents->GetXaxis()->SetBinLabel(14, aod::filtering::OmegaHighMult::columnLabel());
     hProcessedEvents->GetXaxis()->SetBinLabel(15, aod::filtering::DoubleOmega::columnLabel());
     hProcessedEvents->GetXaxis()->SetBinLabel(16, aod::filtering::OmegaXi::columnLabel());
+    hProcessedEvents->GetXaxis()->SetBinLabel(17, "LL");
 
     hCandidate->GetXaxis()->SetBinLabel(1, "All");
     hCandidate->GetXaxis()->SetBinLabel(2, "Has_V0");
@@ -439,7 +545,7 @@
   {
     const auto sigma = parSigmaMass->get(0u, s) * exp(parSigmaMass->get(1, s) * pt) + parSigmaMass->get(2, s) * exp(parSigmaMass->get(3, s) * pt);
     return nsigma * sigma;
   }
 
   ////////////////////////////////////////////////////////
   ////////// Strangeness Filter //////////////////////////
@@ -447,14 +553,14 @@
 
   void fillTriggerTable(bool keepEvent[])
   {
-    strgtable(keepEvent[0], keepEvent[1], keepEvent[2], keepEvent[3], keepEvent[4], keepEvent[5], keepEvent[6], keepEvent[7], keepEvent[8], keepEvent[9], keepEvent[10], keepEvent[11]);
+    strgtable(keepEvent[0], keepEvent[1], keepEvent[2], keepEvent[3], keepEvent[4], keepEvent[5], keepEvent[6], keepEvent[7], keepEvent[8], keepEvent[9], keepEvent[10], keepEvent[11], keepEvent[12]);
   }
 
-  void process(CollisionCandidates const& collision, TrackCandidates const& tracks, aod::Cascades const& cascadesBase, aod::AssignedTrackedCascades const& trackedCascades, aod::AssignedTrackedV0s const& /*trackedV0s*/, aod::AssignedTracked3Bodys const& /*tracked3Bodys*/, aod::V0s const&, aod::BCs const&, aod::FT0s const& /*ft0s*/)
+  void process(CollisionCandidates const& collision, TrackCandidates const& tracks, aod::Cascades const& cascadesBase, aod::AssignedTrackedCascades const& trackedCascades, aod::AssignedTrackedV0s const& /*trackedV0s*/, aod::AssignedTracked3Bodys const& /*tracked3Bodys*/, aod::V0s const& v0Base, aod::BCs const&, aod::FT0s const& /*ft0s*/)
   {
     // Is event good? [0] = Omega, [1] = high-pT hadron + Omega, [2] = 2Xi, [3] = 3Xi, [4] = 4Xi, [5] single-Xi, [6] Omega with high radius
     // [7] tracked Xi, [8] tracked Omega, [9] Omega + high mult event
-    bool keepEvent[12]{}; // explicitly zero-initialised
+    bool keepEvent[13]{}; // explicitly zero-initialised
     std::vector<std::array<int64_t, 2>> v0sFromOmegaID;
     std::vector<std::array<int64_t, 2>> v0sFromXiID;
 
@@ -486,7 +592,7 @@
 
     Bool_t isHighMultEvent = 0;
     float multFT0MNorm = 0.f;
     EventsvsMultiplicity.fill(HIST("AllEventsvsMultiplicityFT0M"), collision.multFT0M());
     if (!useNormalisedMult) {
       if (collision.multFT0M() > LowLimitFT0MMult) {
         isHighMultEvent = 1;
@@ -571,6 +677,102 @@
     const auto primaryVertex = getPrimaryVertex(collision);
     o2::dataformats::DCA impactParameterTrk;
 
+    std::vector<std::tuple<int64_t, int64_t, TVector3, TVector3>> v0sSelTuple;
+    for (auto& v00 : v0Base) { // loop over v0 for pre selection
+      hCandidate->Fill(0.5);   // All candidates
+
+      if (v00.v0Type() != 1) {
+        continue;
+      }
+
+      const auto posTrack0 = v00.posTrack_as<TrackCandidates>();
+      const auto negTrack0 = v00.negTrack_as<TrackCandidates>();
+
+      if (!isSelectedV0Daughter(posTrack0) || !isSelectedV0Daughter(negTrack0)) {
+        continue;
+      }
+
+      auto trackParPos0 = getTrackParCov(posTrack0);
+      auto trackParNeg0 = getTrackParCov(negTrack0);
+
+      if (!mStraHelper.buildV0Candidate(v00.collisionId(), pvPos[0], pvPos[1], pvPos[2], posTrack0, negTrack0, trackParPos0, trackParNeg0)) {
+        continue;
+      }
+
+      if (std::hypot(mStraHelper.v0.position[0], mStraHelper.v0.position[1]) < cfgLLCuts.cfgv0radiusMin) {
+        continue;
+      }
+      if (std::fabs(mStraHelper.v0.positiveDCAxy) < cfgLLCuts.cfgDCAPosToPVMin) {
+        continue;
+      }
+      if (std::fabs(mStraHelper.v0.negativeDCAxy) < cfgLLCuts.cfgDCANegToPVMin) {
+        continue;
+      }
+      if (TMath::Cos(mStraHelper.v0.pointingAngle) < cfgLLCuts.cfgv0CosPA) {
+        continue;
+      }
+      if (std::fabs(mStraHelper.v0.daughterDCA) > cfgLLCuts.cfgDCAV0Dau) {
+        continue;
+      }
+      if (std::hypot(mStraHelper.v0.momentum[0], mStraHelper.v0.momentum[1]) < cfgLLCuts.cfgV0PtMin) {
+        continue;
+      }
+      double yLambda = RecoDecay::y(array{mStraHelper.v0.momentum[0], mStraHelper.v0.momentum[1], mStraHelper.v0.momentum[2]}, o2::constants::physics::MassLambda0);
+      if (yLambda < cfgLLCuts.cfgV0RapMin) {
+        continue;
+      }
+      if (yLambda > cfgLLCuts.cfgV0RapMax) {
+        continue;
+      }
+      double distovertotmom = std::hypot(mStraHelper.v0.position[0] - collision.posX(), mStraHelper.v0.position[1] - collision.posY(), mStraHelper.v0.position[2] - collision.posZ()) / (std::hypot(mStraHelper.v0.momentum[0], mStraHelper.v0.momentum[1], mStraHelper.v0.momentum[2]) + 1e-13);
+      if (distovertotmom * o2::constants::physics::MassLambda0 > cfgLLCuts.cfgV0LifeTime) {
+        continue;
+      }
+
+      int Tag = 0;
+      if (isSelectedV0DaughterPID(posTrack0, 0) && isSelectedV0DaughterPID(negTrack0, 1)) {
+        if (cfgLLCuts.cfgLambdaMassWindow > std::fabs(mStraHelper.v0.massLambda - o2::constants::physics::MassLambda0)) {
+          if (cfgLLCuts.cfgCompV0Rej < std::fabs(mStraHelper.v0.massK0Short - o2::constants::physics::MassLambda0)) {
+            Tag++;
+          }
+        }
+      } // lambda
+      if (isSelectedV0DaughterPID(posTrack0, 1) && isSelectedV0DaughterPID(negTrack0, 0)) {
+        if (cfgLLCuts.cfgLambdaMassWindow > std::fabs(mStraHelper.v0.massAntiLambda - o2::constants::physics::MassLambda0)) {
+          if (cfgLLCuts.cfgCompV0Rej < std::fabs(mStraHelper.v0.massK0Short - o2::constants::physics::MassLambda0)) {
+            Tag++;
+          }
+        }
+      } // anti lambda
+      if (Tag != 1) { // Select when only one hypothesis is satisfied
+        continue;
+      }
+
+      TVector3 v0pos(mStraHelper.v0.position[0], mStraHelper.v0.position[1], mStraHelper.v0.position[2]);
+      TVector3 v0mom(mStraHelper.v0.momentum[0], mStraHelper.v0.momentum[1], mStraHelper.v0.momentum[2]);
+
+      v0sSelTuple.emplace_back(posTrack0.globalIndex(), negTrack0.globalIndex(), v0pos, v0mom);
+    }
+
+    for (size_t i = 0; i < v0sSelTuple.size(); ++i) {
+      for (size_t j = i + 1; j < v0sSelTuple.size(); ++j) {
+        auto d00 = std::get<0>(v0sSelTuple[i]);
+        auto d01 = std::get<1>(v0sSelTuple[i]);
+        auto d10 = std::get<0>(v0sSelTuple[j]);
+        auto d11 = std::get<1>(v0sSelTuple[j]);
+        if (d00 == d10 || d00 == d11 || d01 == d10 || d01 == d11) {
+          continue;
+        }
+        auto v00pos = std::get<2>(v0sSelTuple[i]);
+        auto v00mom = std::get<3>(v0sSelTuple[i]);
+        auto v01pos = std::get<2>(v0sSelTuple[j]);
+        auto v01mom = std::get<3>(v0sSelTuple[j]);
+        if (isSelectedV0V0(v00pos, v00mom, v01pos, v01mom)) {
+          keepEvent[12] = true;
+        }
+      }
+    }
+
     for (auto& casc : cascadesBase) { // loop over cascades
       hCandidate->Fill(0.5); // All candidates
 
@@ -610,7 +812,7 @@
       // pointing angle
       double v0DauCPA = RecoDecay::cpa(pvPos, array{mStraHelper.cascade.v0Position[0], mStraHelper.cascade.v0Position[1], mStraHelper.cascade.v0Position[2]}, array{mStraHelper.cascade.positiveMomentum[0] + mStraHelper.cascade.negativeMomentum[0], mStraHelper.cascade.positiveMomentum[1] + mStraHelper.cascade.negativeMomentum[1], mStraHelper.cascade.positiveMomentum[2] + mStraHelper.cascade.negativeMomentum[2]});
       double cascCPA = RecoDecay::cpa(
         pvPos,
         array{mStraHelper.cascade.cascadePosition[0], mStraHelper.cascade.cascadePosition[1], mStraHelper.cascade.cascadePosition[2]},
         array{mStraHelper.cascade.positiveMomentum[0] + mStraHelper.cascade.negativeMomentum[0] + mStraHelper.cascade.bachelorMomentum[0], mStraHelper.cascade.positiveMomentum[1] + mStraHelper.cascade.negativeMomentum[1] + mStraHelper.cascade.bachelorMomentum[1], mStraHelper.cascade.positiveMomentum[2] + mStraHelper.cascade.negativeMomentum[2] + mStraHelper.cascade.bachelorMomentum[2]});
       // dca V0 to PV
@@ -625,14 +827,14 @@
       if (mStraHelper.cascade.charge < 0) {
         LambdaMass = RecoDecay::m(array{array{mStraHelper.cascade.positiveMomentum[0], mStraHelper.cascade.positiveMomentum[1], mStraHelper.cascade.positiveMomentum[2]}, array{mStraHelper.cascade.negativeMomentum[0], mStraHelper.cascade.negativeMomentum[1], mStraHelper.cascade.negativeMomentum[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged});
       } else {
         LambdaMass = RecoDecay::m(array{array{mStraHelper.cascade.positiveMomentum[0], mStraHelper.cascade.positiveMomentum[1], mStraHelper.cascade.positiveMomentum[2]}, array{mStraHelper.cascade.negativeMomentum[0], mStraHelper.cascade.negativeMomentum[1], mStraHelper.cascade.negativeMomentum[2]}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton});
       }
 
       // rapidity
       double yXi = RecoDecay::y(array{mStraHelper.cascade.bachelorMomentum[0] + mStraHelper.cascade.positiveMomentum[0] + mStraHelper.cascade.negativeMomentum[0], mStraHelper.cascade.bachelorMomentum[1] + mStraHelper.cascade.positiveMomentum[1] + mStraHelper.cascade.negativeMomentum[1], mStraHelper.cascade.bachelorMomentum[2] + mStraHelper.cascade.positiveMomentum[2] + mStraHelper.cascade.negativeMomentum[2]}, o2::constants::physics::MassXiMinus);
       double yOmega = RecoDecay::y(array{mStraHelper.cascade.bachelorMomentum[0] + mStraHelper.cascade.positiveMomentum[0] + mStraHelper.cascade.negativeMomentum[0], mStraHelper.cascade.bachelorMomentum[1] + mStraHelper.cascade.positiveMomentum[1] + mStraHelper.cascade.negativeMomentum[1], mStraHelper.cascade.bachelorMomentum[2] + mStraHelper.cascade.positiveMomentum[2] + mStraHelper.cascade.negativeMomentum[2]}, o2::constants::physics::MassOmegaMinus);
 
       if (mStraHelper.cascade.charge > 0) {
         if (std::fabs(mStraHelper.cascade.positiveDCAxy) < dcamesontopv) {
           continue;
         }
@@ -1121,7 +1323,9 @@
     if (keepEvent[11]) {
       hProcessedEvents->Fill(14.5);
     }
-
+    if (keepEvent[12]) {
+      hProcessedEvents->Fill(15.5);
+    }
     // Filling the table
     fillTriggerTable(keepEvent);
   }

EventFiltering/filterTables.h

@@ -154,6 +154,7 @@ DECLARE_SOA_COLUMN(TrackedXi, hasTrackedXi, bool);               //! at least 1
 DECLARE_SOA_COLUMN(TrackedOmega, hasTrackedOmega, bool);         //! at least 1 tracked Omega
 DECLARE_SOA_COLUMN(Tracked3Body, hasTracked3Body, bool);         //! at least 1 tracked 3Body
 DECLARE_SOA_COLUMN(OmegaHighMult, hasOmegaHighMult, bool);       //! at least 1 Omega + high-mult event
+DECLARE_SOA_COLUMN(LambdaLambda, lambdaLambda, bool);            //! at least 2 lambda satisfying selection
 
 // F1-proton
 DECLARE_SOA_COLUMN(TriggerEventF1Proton, triggereventf1proton, bool); //! F1 - proton femto trigger event
@@ -295,7 +296,7 @@ using FullJetFilter = FullJetFilters::iterator;
 
 // strangeness (lf)
 DECLARE_SOA_TABLE(StrangenessFilters, "AOD", "LFStrgFilters", //!
-                  filtering::Omega, filtering::hadronOmega, filtering::DoubleXi, filtering::TripleXi, filtering::QuadrupleXi, filtering::SingleXiYN, filtering::OmegaLargeRadius, filtering::TrackedXi, filtering::TrackedOmega, filtering::OmegaHighMult, filtering::DoubleOmega, filtering::OmegaXi);
+                  filtering::Omega, filtering::hadronOmega, filtering::DoubleXi, filtering::TripleXi, filtering::QuadrupleXi, filtering::SingleXiYN, filtering::OmegaLargeRadius, filtering::TrackedXi, filtering::TrackedOmega, filtering::OmegaHighMult, filtering::DoubleOmega, filtering::OmegaXi, filtering::LambdaLambda);
 
 using StrangenessFilter = StrangenessFilters::iterator;