[PWGCF] pi-hypertriton ME process

PWGCF/Femto/FemtoNuclei/TableProducer/PiNucleiFemto.cxx

@@ -63,6 +63,7 @@
 #include <array>
 #include <cmath>
 #include <cstdlib>
+#include <deque>
 #include <iterator> // std::prev
 #include <string>
 #include <vector>
@@ -204,12 +205,13 @@
   PresliceUnsorted<o2::aod::DataHypCandsWColl> hypPerCol = o2::aod::hyperrec::collisionId;
 
   // binning for EM background
-  ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "Binning for multiplicity"};
+  ConfigurableAxis axisVertex{"axisVertex", {30, -10, 10}, "Binning for vtxz"};
   ConfigurableAxis axisCentrality{"axisCentrality", {40, 0, 100}, "Binning for centrality"};
   using BinningType = ColumnBinningPolicy<aod::collision::PosZ, aod::cent::CentFT0C>;
   BinningType binningPolicy{{axisVertex, axisCentrality}, true};
   SliceCache cache;
   SameKindPair<CollisionsFull, TrackCandidates, BinningType> mPair{binningPolicy, settingNoMixedEvents, -1, &cache};
+  // Pair<CollisionsFull, TrackCandidates, o2::aod::DataHypCandsWColl, BinningType> hyperPair{binningPolicy, settingNoMixedEvents, -1, &cache};
 
   std::array<float, 6> mBBparamsDe;
 
@@ -229,6 +231,7 @@
     "QA",
     {{"hVtxZ", "Vertex distribution in Z;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}},
      {"hNcontributor", "Number of primary vertex contributor", {HistType::kTH1F, {{2000, 0.0f, 2000.0f}}}},
+     {"hCentrality", "Centrality", {HistType::kTH1F, {{200, 0.0f, 100.0f}}}},
      {"hTrackSel", "Accepted tracks", {HistType::kTH1F, {{Selections::kAll, -0.5, static_cast<double>(Selections::kAll) - 0.5}}}},
      {"hEvents", "; Events;", {HistType::kTH1F, {{3, -0.5, 2.5}}}},
      {"hEmptyPool", "svPoolCreator did not find track pairs false/true", {HistType::kTH1F, {{2, -0.5, 1.5}}}},
@@ -240,6 +243,9 @@
      {"hNuPitInvMass", "; M(Nu + p) (GeV/#it{c}^{2})", {HistType::kTH1F, {{300, 3.74f, 4.34f}}}},
      {"hNuPt", "#it{p}_{T} distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}},
      {"hPiPt", "Pt distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}},
+     {"hHe3TPCnsigma", "NsigmaHe3 TPC distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TPC}(He3)", {HistType::kTH2F, {{100, -2.0f, 2.0f}, {200, -5.0f, 5.0f}}}},
+     {"hHe3P", "Pin distribution; p (GeV/#it{c})", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}},
+     {"hHe3P_preselected", "Pin distribution_preselected; p (GeV/#it{c})", {HistType::kTH1F, {{120, -3.0f, 3.0f}}}},
      {"hNuEta", "eta distribution; #eta(Nu)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}},
      {"hPiEta", "eta distribution; #eta(#pi)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}},
      {"hNuPhi", "phi distribution; phi(Nu)", {HistType::kTH1F, {{600, -4.0f, 4.0f}}}},
@@ -268,6 +274,52 @@
     false,
     true};
 
+  int numOfCentBins = 40;
+  int numOfVertexZBins = 30;
+  float Vz_low = -10.0f;
+  float Vz_high = 10.0f;
+  float Vz_step = (Vz_high - Vz_low) / numOfVertexZBins;
+
+  struct EventRef {
+    int64_t collisionId;
+  };
+
+  struct PoolBin {
+    std::deque<EventRef> events;
+  };
+
+  std::vector<PoolBin> All_Event_pool;
+  bool isInitialized = false;
+
+  int nPoolBins() const { return numOfVertexZBins * numOfCentBins; }
+
+  void initializePools()
+  {
+    All_Event_pool.clear();
+    All_Event_pool.resize(nPoolBins());
+    isInitialized = true;
+  }
+
+  int where_pool(float vz, float v0Centr) const
+  {
+    float CentBinWidth = 100.0 / numOfCentBins; // = 2.5
+
+    int iy = static_cast<int>(std::floor(v0Centr / CentBinWidth));
+    if (iy < 0)
+      iy = 0;
+    if (iy >= numOfCentBins)
+      iy = numOfCentBins - 1;
+
+    int ix = static_cast<int>(std::floor((vz - Vz_low) / Vz_step));
+    if (ix < 0)
+      ix = 0;
+    if (ix >= numOfVertexZBins)
+      ix = numOfVertexZBins - 1;
+
+    int bin = ix + numOfVertexZBins * iy;
+    return bin;
+  }
+
   void init(o2::framework::InitContext&)
   {
     mZorroSummary.setObject(mZorro.getZorroSummary());
@@ -488,6 +540,41 @@
     return false;
   }
 
+  float averageClusterSizeCosl(uint32_t itsClusterSizes, float eta)
+  {
+    float average = 0;
+    int nclusters = 0;
+    const float cosl = 1. / std::cosh(eta);
+    const int nlayerITS = 7;
+
+    for (int layer = 0; layer < nlayerITS; layer++) {
+      if ((itsClusterSizes >> (layer * 4)) & 0xf) {
+        nclusters++;
+        average += (itsClusterSizes >> (layer * 4)) & 0xf;
+      }
+    }
+    if (nclusters == 0) {
+      return 0;
+    }
+    return average * cosl / nclusters;
+  };
+
+  bool selectionPIDHyper(const aod::DataHypCandsWColl::iterator& V0Hyper)
+  {
+    mQaRegistry.fill(HIST("hHe3P_preselected"), V0Hyper.tpcMomHe());
+    float averClusSizeHe = averageClusterSizeCosl(V0Hyper.itsClusterSizesHe(), V0Hyper.etaHe3());
+    if (averClusSizeHe <= 4) {
+      return false;
+    }
+    if (V0Hyper.tpcChi2He() <= 0.5) {
+      return false;
+    }
+    mQaRegistry.fill(HIST("hHe3P"), V0Hyper.tpcMomHe());
+    mQaRegistry.fill(HIST("hHe3TPCnsigma"), V0Hyper.nSigmaHe());
+
+    return true;
+  }
+
   // ==================================================================================================================
 
   template <typename Ttrack, typename Tcollisions, typename Ttracks>
@@ -720,21 +807,23 @@
   template <typename Ttrack, typename Thypers>
   void pairTracksSameEventHyper(const Ttrack& piTracks, const Thypers& V0Hypers)
   {
-    for (const auto& piTrack : piTracks) {
-
-      mQaRegistry.fill(HIST("hTrackSel"), Selections::kNoCuts);
-
-      if (!selectTrack(piTrack)) {
+    for (const auto& V0Hyper : V0Hypers) {
+      if (!selectionPIDHyper(V0Hyper)) {
         continue;
       }
-      mQaRegistry.fill(HIST("hTrackSel"), Selections::kTrackCuts);
+      for (const auto& piTrack : piTracks) {
 
-      if (!selectionPIDPion(piTrack)) {
-        continue;
-      }
-      mQaRegistry.fill(HIST("hTrackSel"), Selections::kPID);
+        mQaRegistry.fill(HIST("hTrackSel"), Selections::kNoCuts);
 
-      for (const auto& V0Hyper : V0Hypers) {
+        if (!selectTrack(piTrack)) {
+          continue;
+        }
+        mQaRegistry.fill(HIST("hTrackSel"), Selections::kTrackCuts);
+
+        if (!selectionPIDPion(piTrack)) {
+          continue;
+        }
+        mQaRegistry.fill(HIST("hTrackSel"), Selections::kPID);
 
         SVCand pair;
         pair.tr0Idx = piTrack.globalIndex();
@@ -770,6 +859,29 @@
     }
   }
 
+  template <typename T1, typename T2>
+  void pairHyperEventMixing(T1& pionCands, T2& hypCands)
+  {
+    for (const auto& hypCand : hypCands) {
+      if (!selectionPIDHyper(hypCand)) {
+        continue;
+      }
+      for (const auto& pionCand : pionCands) {
+        if (!selectTrack(pionCand) || !selectionPIDPion(pionCand)) {
+          continue;
+        }
+
+        SVCand pair;
+        pair.tr0Idx = hypCand.globalIndex();
+        pair.tr1Idx = pionCand.globalIndex();
+        const int collIdx = hypCand.collisionId();
+        CollBracket collBracket{collIdx, collIdx};
+        pair.collBracket = collBracket;
+        mTrackHypPairs.push_back(pair);
+      }
+    }
+  }
+
   template <typename Tcoll>
   void fillTable(const PiNucandidate& piNucand, const Tcoll& collision)
   {
@@ -856,7 +968,7 @@
     mQaRegistry.fill(HIST("hNuPitInvMass"), piNucand.invMass);
     mQaRegistry.fill(HIST("hdcaxyNu"), piNucand.dcaxyNu);
     mQaRegistry.fill(HIST("hdcazNu"), piNucand.dcazNu);
     mQaRegistry.fill(HIST("hdcazNu_min"), (abs(piNucand.dcazNu) - settingCutDeDCAzMin));
     mQaRegistry.fill(HIST("hNClsNuITS"), piNucand.nClsItsNu);
     mQaRegistry.fill(HIST("hNClsPiITS"), piNucand.nClsItsPi);
     mQaRegistry.fill(HIST("hisBkgEM"), piNucand.isBkgEM);
@@ -867,7 +979,7 @@
     double m_BBparamsProton[6] = {-54.42066571222577, 0.2857381250239097, 1.247140602468868, 0.6297483918147729, 2.985438833884555, 0.09};
 
     float TPCinnerParam = InnerParamTPCHad;
     float expTPCSignal = o2::tpc::BetheBlochAleph((TPCinnerParam / 0.9382721), m_BBparamsProton[0], m_BBparamsProton[1], m_BBparamsProton[2], m_BBparamsProton[3], m_BBparamsProton[4]);
     double resoTPC{expTPCSignal * m_BBparamsProton[5]};
     return ((SignalTPCHad - expTPCSignal) / resoTPC);
   }
@@ -878,7 +990,7 @@
     const float kp0 = 1.22204e-02;
     const float kp1 = 7.48467e-01;
 
     double fSigmaTOFMassHad = (kp0 * TMath::Exp(kp1 * TMath::Abs(ptHad))) * fExpTOFMassHad;
     double fNSigmaTOFHad = (MassTOFHad - fExpTOFMassHad) / fSigmaTOFMassHad;
     return fNSigmaTOFHad;
   }
@@ -912,12 +1024,12 @@
   {
     double PrTPCnsigma = computePrTPCnsig(piNucand.momPiTPC, piNucand.tpcSignalPi);
     double PrTOFnsigma = tofNSigmaCalculation(piNucand.massTOFPi, piNucand.recoPtPi());
     if (abs(PrTPCnsigma) < settingCutNsigTPCPrMin)
       return;
     if (abs(PrTOFnsigma) < settingCutNsigTOFPrMin)
       return;
     float DeDCAxyMin = 0.015 + 0.0305 / TMath::Power(piNucand.recoPtNu(), 1.1);
     if (abs(piNucand.dcaxyNu) > DeDCAxyMin || abs(piNucand.dcazNu) > settingCutDeDCAzMin || abs(piNucand.dcaxyPi) > settingCutPiDCAxyMin || abs(piNucand.dcazPi) > settingCutPiDCAzMin)
       return;
     fillHistograms(piNucand);
 
@@ -977,6 +1089,15 @@
       if (!fillCandidateInfoHyper(v0hyper, piTrack, piNucand, isMixedEvent)) {
         continue;
       }
+      mQaRegistry.fill(HIST("hNuPt"), piNucand.recoPtNu());
+      mQaRegistry.fill(HIST("hPiPt"), piNucand.recoPtPi());
+      mQaRegistry.fill(HIST("hNuEta"), piNucand.recoEtaNu());
+      mQaRegistry.fill(HIST("hPiEta"), piNucand.recoEtaPi());
+      mQaRegistry.fill(HIST("hNuPhi"), piNucand.recoPhiNu());
+      mQaRegistry.fill(HIST("hPiPhi"), piNucand.recoPhiPi());
+      mQaRegistry.fill(HIST("hNuPitInvMass"), piNucand.invMass);
+      mQaRegistry.fill(HIST("hNClsPiITS"), piNucand.nClsItsPi);
+      mQaRegistry.fill(HIST("hisBkgEM"), piNucand.isBkgEM);
 
       auto collision = collisions.rawIteratorAt(piNucand.collisionID);
 
@@ -1021,6 +1142,7 @@
   {
     mGoodCollisions.clear();
     mGoodCollisions.resize(collisions.size(), false);
+    // LOG(info) << "Number of hyperCandidates read = " << V0Hypers.size();
 
     for (const auto& collision : collisions) {
 
@@ -1068,6 +1190,76 @@
     fillPairs(collisions, tracks, /*isMixedEvent*/ true);
   }
   PROCESS_SWITCH(PiNucleiFemto, processMixedEvent, "Process Mixed event", false);
+
+  /*void processMixedEventHyper(const CollisionsFull& collisions, o2::aod::DataHypCandsWColl const& V0Hypers, const TrackCandidates& pitracks)
+  {
+    LOG(debug) << "Processing mixed event for hypertriton";
+    mTrackHypPairs.clear();
+
+    for (const auto& [c1, tracks1, c2, V0Hypers2] : hyperPair) {
+      if (!c1.sel8() || !c2.sel8()) {
+        continue;
+      }
+
+      mQaRegistry.fill(HIST("hNcontributor"), c2.numContrib());
+      //mQaRegistry.fill(HIST("hCentrality"), c2.centFT0C());
+      mQaRegistry.fill(HIST("hVtxZ"), c2.posZ());
+
+      pairHyperEventMixing(tracks1, V0Hypers2);
+    }
+}
+PROCESS_SWITCH(PiNucleiFemto, processMixedEventHyper, "Process Mixed event", false);*/
+
+  void processMixedEventHyperPool(const CollisionsFull& collisions, o2::aod::DataHypCandsWColl const& V0Hypers, const TrackCandidates& pitracks)
+  {
+    LOG(info) << "Processing mixed event for hypertriton";
+
+    mTrackHypPairs.clear();
+    if (!isInitialized) {
+      initializePools();
+      LOG(info) << "Initialized event pool with size = " << All_Event_pool.size();
+    }
+    for (auto const& collision : collisions) {
+      int poolIndexPi = where_pool(collision.posZ(), collision.centFT0C());
+      auto& pool = All_Event_pool[poolIndexPi];
+
+      if (poolIndexPi < 0 || static_cast<size_t>(poolIndexPi) >= All_Event_pool.size()) {
+        continue;
+      }
+
+      for (auto const& storedEvent : pool.events) {
+        auto c1 = collisions.iteratorAt(storedEvent.collisionId);
+        const auto& c2 = collision;
+        if (!c1.sel8() || !c2.sel8())
+          continue;
+
+        std::vector<TrackCandidates::iterator> tracks1;
+        for (auto const& t : pitracks) {
+          if (t.collisionId() == c1.globalIndex()) {
+            tracks1.push_back(t);
+          }
+        }
+
+        std::vector<o2::aod::DataHypCandsWColl::iterator> hypers2;
+        for (auto const& h : V0Hypers) {
+          if (h.collisionId() != c2.globalIndex())
+            continue;
+          int poolIndexHyp = where_pool(h.zPrimVtx(), h.centralityFT0C());
+          if (poolIndexHyp != poolIndexPi)
+            continue;
+          hypers2.push_back(h);
+        }
+        pairHyperEventMixing(tracks1, hypers2);
+      }
+      fillPairsHyper(collisions, pitracks, V0Hypers, /*isMixedEvent*/ true);
+
+      if (static_cast<int>(pool.events.size()) >= settingNoMixedEvents) {
+        pool.events.pop_front();
+      }
+      pool.events.push_back({collision.globalIndex()});
+    }
+  }
+  PROCESS_SWITCH(PiNucleiFemto, processMixedEventHyperPool, "Process Mixed event", false);
 };
 
 WorkflowSpec defineDataProcessing(const ConfigContext& cfgc)