Update PiNucleiFemto.cxx

Author: blacwovie

PWGCF/Femto/FemtoNuclei/TableProducer/PiNucleiFemto.cxx

@@ -19,6 +19,7 @@
 #include "PWGLF/DataModel/EPCalibrationTables.h"
 #include "PWGLF/DataModel/LFHypernucleiTables.h"
 #include "PWGLF/Utils/svPoolCreator.h"
+#include "PWGCF/FemtoWorld/Core/FemtoWorldMath.h"
 
 #include "Common/Core/PID/PIDTOF.h"
 #include "Common/Core/PID/TPCPIDResponse.h"
@@ -76,7 +77,7 @@ using std::array;
 using CollBracket = o2::math_utils::Bracket<int>;
 using CollisionsFull = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::FT0Mults>;
 using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0Cs, aod::FT0Mults>;
-using TrackCandidates = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullDe, aod::pidTOFFullDe, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::TOFSignal, aod::TOFEvTime>;
+using TrackCandidates = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullDe, aod::pidTOFFullDe, aod::pidTPCFullPr, aod::pidTOFFullPr, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::TOFSignal, aod::TOFEvTime>;
 
 namespace
 {
@@ -125,6 +126,8 @@ struct PiNucandidate {
   float chi2TPCPi = -10.f;
   float nSigmaNu = -10.f;
   float nSigmaPi = -10.f;
+  float tpcPrnsigma = -10.f;
+  float tofPrnsigma = -10.f;
   uint32_t pidTrkNu = 0xFFFFF; // PID in tracking
   uint32_t pidTrkPi = 0xFFFFF;
   float massTOFNu = -10;
@@ -141,8 +144,12 @@ struct PiNucandidate {
   int trackIDNu = -1;
   int trackIDPi = -1;
 
+  float kstar = 1.f;
+  float mT = 1.f;
+
   // collision information
   int32_t collisionID = 0;
+  float cent =1.f;
 };
 
 struct PiNucleiFemto {
@@ -237,6 +244,7 @@ struct PiNucleiFemto {
      {"hNcontributor", "Number of primary vertex contributor", {HistType::kTH1F, {{2000, 0.0f, 2000.0f}}}},
      {"hCentrality", "Centrality", {HistType::kTH1F, {{100, 0.0f, 100.0f}}}},
      {"hTrackSel", "Accepted tracks", {HistType::kTH1F, {{Selections::kAll, -0.5, static_cast<double>(Selections::kAll) - 0.5}}}},
+     {"hSkipReasons", "Why storedEvent skipped;Reason;Counts", {HistType::kTH1F, {{5, -0.5, 4.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}}}},
      {"hdcaxyNu", ";DCA_{xy} (cm)", {HistType::kTH1F, {{200, -1.0f, 1.0f}}}},
@@ -247,6 +255,8 @@ struct PiNucleiFemto {
      {"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}}}},
+     {"hSingleNuPt", "#it{p}_{T} distribution; #it{p}_{T} (GeV/#it{c})", {HistType::kTH1F, {{240, -6.0f, 6.0f}}}},
+     {"hSinglePiPt", "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}}}},
@@ -268,16 +278,24 @@ struct PiNucleiFemto {
      {"h2NsigmaPiTOF", "NsigmaPi TOF distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(p)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}}},
      {"h2NsigmaNuTOF", "NsigmaNu TOF distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(Nu)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {200, -5.0f, 5.0f}}}},
      {"h2NsigmaPiTOF_preselection", "NsigmaPi TOF distribution; #iit{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(p)", {HistType::kTH2F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}}}},
+     {"hkStaVsmTVsCent_LS_M", ";kStar (GeV/c);mT (GeV/#it{c}^{2});Centrality", {HistType::kTH3F, {{300, 0.0f, 3.0f}, {100, 0.2f, 3.2f}, {100, 0.0f, 100.0f}}}},
+     {"hkStaVsmTVsCent_LS_A", ";kStar (GeV/c);mT (GeV/#it{c}^{2});Centrality", {HistType::kTH3F, {{300, 0.0f, 3.0f}, {100, 0.2f, 3.2f}, {100, 0.0f, 100.0f}}}},
+     {"hkStaVsmTVsCent_US_M", ";kStar (GeV/c);mT (GeV/#it{c}^{2});Centrality", {HistType::kTH3F, {{300, 0.0f, 3.0f}, {100, 0.2f, 3.2f}, {100, 0.0f, 100.0f}}}},
+     {"hkStaVsmTVsCent_US_A", ";kStar (GeV/c);mT (GeV/#it{c}^{2});Centrality", {HistType::kTH3F, {{300, 0.0f, 3.0f}, {100, 0.2f, 3.2f}, {100, 0.0f, 100.0f}}}},
+     {"hCollIDVsCentEachPion", ";CollisionID;Centrality", {HistType::kTH2F, {{4000, 0.0f, 4000.0f}, {100, 0.0f, 100.0f}}}},
+     {"hNHypsPerPrevColl", "Number of V0Hypers in previous collision used for mixing;N_{V0Hypers};Entries",{HistType::kTH2F, {{4000, 0.0f, 4000.0f}, {50, -0.5, 49.5}}}}, 
      {"hkStar_LS_M", ";kStar (GeV/c)", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
      {"hkStar_LS_A", ";kStar (GeV/c)", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
      {"hkStar_US_M", ";kStar (GeV/c)", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
      {"hkStar_US_A", ";kStar (GeV/c)", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
-     {"hkStar_All", ";kStar (GeV/c)", {HistType::kTH1F, {{300, 0.0f, 3.0f}}}},
+     {"h2NsigmaPiPrTPC", "NsigmaPi TPC distribution; #it{p}_{T}(GeV/#it{c}); n#sigma_{TPC}(p)", {HistType::kTH1F, {{200, -5.0f, 5.0f}}}},
+     {"h2NsigmaPiPrTOF", "NsigmaPi TOF distribution; #it{p}_{T}(GeV/#it{c}); n#sigma_{TPC}(p)", {HistType::kTH1F, {{200, -5.0f, 5.0f}}}},
      {"hisBkgEM", "; isBkgEM;", {HistType::kTH1F, {{3, -1, 2}}}}},
     OutputObjHandlingPolicy::AnalysisObject,
     false,
     true};
 
+ 
   int numOfCentBins = 40;
   int numOfVertexZBins = 30;
   float Vz_low = -10.0f;
@@ -460,10 +478,21 @@ struct PiNucleiFemto {
   template <typename Ttrack>
   bool selectionPIDPion(const Ttrack& candidate)
   {
+    if (abs(candidate.dcaXY()) > settingCutPiDCAxyMin || abs(candidate.dcaZ()) > settingCutPiDCAzMin)
+      return false;
+
     auto tpcNSigmaPi = candidate.tpcNSigmaPi();
     mQaRegistry.fill(HIST("h2NsigmaPiTPC_preselection"), candidate.tpcInnerParam(), tpcNSigmaPi);
     if (std::abs(candidate.pt()) < settingCutPiptMin || std::abs(candidate.pt()) > settingCutPiptMax)
       return false;
+      //reject protons
+    if (std::abs(candidate.tpcNSigmaPr()) < settingCutNsigTPCPrMin)
+      return false;
+    mQaRegistry.fill(HIST("h2NsigmaPiPrTPC"), candidate.tpcNSigmaPr());
+    if (candidate.hasTOF() && std::abs(candidate.tofNSigmaPr()) < settingCutNsigTOFPrMin)
+      return false;
+    mQaRegistry.fill(HIST("h2NsigmaPiPrTOF"), candidate.tofNSigmaPr());
+      
     if (candidate.hasTOF() && candidate.tpcInnerParam() >= settingCutPinMinTOFPi) {
       auto tofNSigmaPi = candidate.tofNSigmaPi();
       auto combNsigma = std::sqrt(tofNSigmaPi * tofNSigmaPi + tpcNSigmaPi * tpcNSigmaPi);
@@ -499,6 +528,10 @@ struct PiNucleiFemto {
   {
     float tpcInnerParam = candidate.tpcInnerParam();
     mQaRegistry.fill(HIST("h2dEdx"), candidate.sign() * tpcInnerParam, candidate.tpcSignal());
+    
+    float DeDCAxyMin = 0.015 + 0.0305 / TMath::Power(candidate.pt(), 1.1);
+    if (abs(candidate.dcaXY()) > DeDCAxyMin || abs(candidate.dcaXY()) > settingCutDeDCAzMin)
+      return false;
 
     if (std::abs(tpcInnerParam) < settingCutPinMinDe) {
       return false;
@@ -690,6 +723,9 @@ struct PiNucleiFemto {
       piNucand.massTOFPi = trackPi.tpcInnerParam() * std::sqrt(1.f / (beta * beta) - 1.f);
     }
 
+    piNucand.kstar = o2::analysis::femtoWorld::FemtoWorldMath::getkstar(trackPi, o2::constants::physics::MassPiPlus, trackDe, o2::constants::physics::MassDeuteron);
+    piNucand.mT = o2::analysis::femtoWorld::FemtoWorldMath::getmT(trackPi, o2::constants::physics::MassPiPlus, trackDe, o2::constants::physics::MassDeuteron);
+
     return true;
   }
 
@@ -763,11 +799,12 @@ struct PiNucleiFemto {
   }
 
   template <typename Ttrack>
-  void pairTracksSameEvent(const Ttrack& tracks)
+  void pairTracksSameEvent(const Ttrack& tracks, float cent)
   {
+    bool filledAllOnce = false;
     // LOG(info) << "Number of tracks: " << tracks.size();
     for (const auto& track0 : tracks) {
-
+    
       mQaRegistry.fill(HIST("hTrackSel"), Selections::kNoCuts);
 
       if (!selectTrack(track0)) {
@@ -779,6 +816,7 @@ struct PiNucleiFemto {
         continue;
       }
       mQaRegistry.fill(HIST("hTrackSel"), Selections::kPID);
+      mQaRegistry.fill(HIST("hSingleNuPt"), track0.pt() * track0.sign());
 
       for (const auto& track1 : tracks) {
         if (track0 == track1) {
@@ -797,6 +835,11 @@ struct PiNucleiFemto {
         if (!selectTrack(track1) || !selectionPIDPion(track1)) {
           continue;
         }
+   
+        if (!filledAllOnce) {
+          mQaRegistry.fill(HIST("hCollIDVsCentEachPion"), track1.collisionId(), cent);
+          mQaRegistry.fill(HIST("hSinglePiPt"), track1.pt() * track1.sign());
+        }
 
         SVCand trackPair;
         trackPair.tr0Idx = track0.globalIndex();
@@ -806,6 +849,7 @@ struct PiNucleiFemto {
         trackPair.collBracket = collBracket;
         mTrackPairs.push_back(trackPair);
       }
+      filledAllOnce = true; 
     }
   }
 
@@ -956,6 +1000,7 @@ struct PiNucleiFemto {
       piNucand.itsClSizePi,
       piNucand.itsClSizeNu,
       piNucand.sharedClustersPi,
+      piNucand.trackIDPi,
       piNucand.isBkgUS,
       piNucand.isBkgEM);
     if (settingFillMultiplicity) {
@@ -985,80 +1030,26 @@ struct PiNucleiFemto {
     mQaRegistry.fill(HIST("hisBkgEM"), piNucand.isBkgEM);
   }
 
-  double computePrTPCnsig(double InnerParamTPCHad, double SignalTPCHad)
-  {
-    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);
-  }
-
-  double tofNSigmaCalculation(double MassTOFHad, double ptHad)
-  {
-    double fExpTOFMassHad = 0.9487; // Proton mass in TOF
-    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;
-  }
-
-  double computeKstar(const PiNucandidate& piNucand)
-  {
-    constexpr double massDe = o2::constants::physics::MassDeuteron;
-    constexpr double massHad = o2::constants::physics::MassPiPlus;
-
-    const ROOT::Math::PtEtaPhiMVector De(std::abs(piNucand.recoPtNu()), piNucand.recoEtaNu(), piNucand.recoPhiNu(), massDe);
-    const ROOT::Math::PtEtaPhiMVector Had(std::abs(piNucand.recoPtPi()), piNucand.recoEtaPi(), piNucand.recoPhiPi(), massHad);
-    const ROOT::Math::PtEtaPhiMVector trackSum = De + Had;
-
-    const float beta = trackSum.Beta();
-    const float betax = beta * std::cos(trackSum.Phi()) * std::sin(trackSum.Theta());
-    const float betay = beta * std::sin(trackSum.Phi()) * std::sin(trackSum.Theta());
-    const float betaz = beta * std::cos(trackSum.Theta());
-
-    ROOT::Math::PxPyPzMVector DeCMS(De);
-    ROOT::Math::PxPyPzMVector HadCMS(Had);
-
-    const ROOT::Math::Boost boostPRF = ROOT::Math::Boost(-betax, -betay, -betaz);
-    DeCMS = boostPRF(DeCMS);
-    HadCMS = boostPRF(HadCMS);
-
-    const ROOT::Math::PxPyPzMVector RelKstar = DeCMS - HadCMS;
-    return 0.5 * RelKstar.P();
-  }
-
-  void fillKstar(const PiNucandidate& piNucand)
+  template <typename Tcoll>
+  void fillKstar(const PiNucandidate& piNucand, const Tcoll& collision)
   {
-    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);
-
-    double kstar = computeKstar(piNucand);
     if (piNucand.isBkgUS == 0) {
       if (piNucand.recoPtNu() > 0) {
-        mQaRegistry.fill(HIST("hkStar_LS_M"), kstar);
+        mQaRegistry.fill(HIST("hkStar_LS_M"), piNucand.kstar);
+        mQaRegistry.fill(HIST("hkStaVsmTVsCent_LS_M"), piNucand.kstar, piNucand.mT, collision.centFT0C());
       } else {
-        mQaRegistry.fill(HIST("hkStar_LS_A"), kstar);
+        mQaRegistry.fill(HIST("hkStar_LS_A"), piNucand.kstar);
+        mQaRegistry.fill(HIST("hkStaVsmTVsCent_LS_A"), piNucand.kstar, piNucand.mT, collision.centFT0C());
       }
     } else {
       if (piNucand.recoPtNu() > 0) {
-        mQaRegistry.fill(HIST("hkStar_US_M"), kstar);
+        mQaRegistry.fill(HIST("hkStar_US_M"), piNucand.kstar);
+        mQaRegistry.fill(HIST("hkStaVsmTVsCent_US_M"), piNucand.kstar, piNucand.mT, collision.centFT0C());
       } else {
-        mQaRegistry.fill(HIST("hkStar_US_A"), kstar);
+        mQaRegistry.fill(HIST("hkStar_US_A"), piNucand.kstar);
+        mQaRegistry.fill(HIST("hkStaVsmTVsCent_US_A"), piNucand.kstar, piNucand.mT, collision.centFT0C());
       }
     }
-    mQaRegistry.fill(HIST("hkStar_All"), kstar);
   }
 
   // ==================================================================================================================
@@ -1077,9 +1068,9 @@ struct PiNucleiFemto {
         continue;
       }
 
-      fillKstar(piNucand);
-
       auto collision = collisions.rawIteratorAt(piNucand.collisionID);
+      fillKstar(piNucand, collision);
+      fillHistograms(piNucand);
 
       if (settingFillTable) {
         fillTable(piNucand, collision);
@@ -1139,7 +1130,7 @@ struct PiNucleiFemto {
       auto trackTableThisCollision = tracks.sliceBy(mPerCol, collIdx);
       trackTableThisCollision.bindExternalIndices(&tracks);
 
-      pairTracksSameEvent(trackTableThisCollision);
+      pairTracksSameEvent(trackTableThisCollision, collision.centFT0C());
 
       if (mTrackPairs.size() == 0) {
         continue;
@@ -1234,42 +1225,48 @@ PROCESS_SWITCH(PiNucleiFemto, processMixedEventHyper, "Process Mixed event", fal
       mQaRegistry.fill(HIST("hCentrality"), collision.centFT0C());
       mQaRegistry.fill(HIST("hVtxZ"), collision.posZ());
       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;
       }
+      auto& pool = All_Event_pool[poolIndexPi];
 
       for (auto const& storedEvent : pool.events) {
-        auto c1 = collisions.iteratorAt(storedEvent.collisionId);
-        const auto& c2 = collision;
-        if (!c1.sel8() || !c2.sel8())
+        const auto& c2 = collisions.iteratorAt(storedEvent.collisionId);
+        if (!collision.sel8() || !c2.sel8()) {
+          mQaRegistry.fill(HIST("hSkipReasons"), 0);
           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);
+        
+        auto hypdTablepreviousCollision = V0Hypers.sliceBy(hypPerCol, c2.globalIndex());
+        hypdTablepreviousCollision.bindExternalIndices(&V0Hypers);
+        
+        if (hypdTablepreviousCollision.size() == 0) {
+          mQaRegistry.fill(HIST("hSkipReasons"), 1);
+          continue;
+        }
+        
+        auto firstHyp = hypdTablepreviousCollision.iteratorAt(0);
+        int poolIndexHyp = where_pool(firstHyp.zPrimVtx(), firstHyp.centralityFT0C());
+        
+        if (poolIndexHyp != poolIndexPi) {
+         mQaRegistry.fill(HIST("hSkipReasons"), 2);
+         continue;
         }
-        pairHyperEventMixing(tracks1, hypers2);
+        
+        auto trackTableThisCollision = pitracks.sliceBy(mPerCol, collision.globalIndex());
+        trackTableThisCollision.bindExternalIndices(&pitracks);
+
+        mQaRegistry.fill(HIST("hNHypsPerPrevColl"), c2.globalIndex(), hypdTablepreviousCollision.size());
+        pairHyperEventMixing(trackTableThisCollision, hypdTablepreviousCollision);
       }
-      fillPairsHyper(collisions, pitracks, V0Hypers, /*isMixedEvent*/ true);
 
       if (static_cast<int>(pool.events.size()) >= settingNoMixedEvents) {
         pool.events.pop_front();
       }
       pool.events.push_back({collision.globalIndex()});
-    }
+      }
+    fillPairsHyper(collisions, pitracks, V0Hypers, /*isMixedEvent*/ true);
   }
   PROCESS_SWITCH(PiNucleiFemto, processMixedEventHyperPool, "Process Mixed event", false);
 };