added scalar product method and efficiency calculation for polarization

Author: Prottay Das

PWGLF/Tasks/Strangeness/lambdapolsp.cxx

@@ -109,6 +109,7 @@ struct lambdapolsp {
     Configurable<double> etaMix{"etaMix", 0.1, "eta difference in mixing"};
     Configurable<double> ptMix{"ptMix", 0.1, "pt difference in mixing"};
     Configurable<double> phiMix{"phiMix", 0.1, "phi difference in mixing"};
+    Configurable<bool> useSP{"useSP", false, "use scalar product"};
   } randGrp;
   // events
   Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"};
@@ -232,6 +233,7 @@ struct lambdapolsp {
     std::vector<AxisSpec> runaxes = {thnAxisInvMass, axisGrp.configthnAxispT, axisGrp.configthnAxisPol, axisGrp.configcentAxis};
     if (needetaaxis)
       runaxes.insert(runaxes.end(), {axisGrp.configbinAxis});
+    std::vector<AxisSpec> runaxes2 = {thnAxisInvMass, axisGrp.configthnAxispT, axisGrp.configcentAxis};
 
     if (checkwithpub) {
       if (useprofile == 2) {
@@ -428,6 +430,11 @@ struct lambdapolsp {
       histos.add("hptnegantilambda", "hptnegantilambda", HistType::kTH1D, {distGrp.axispt}, true);
     }
 
+    histos.add("hSparseGenLambda", "hSparseGenLambda", HistType::kTHnSparseF, runaxes2, true);
+    histos.add("hSparseGenAntiLambda", "hSparseGenAntiLambda", HistType::kTHnSparseF, runaxes2, true);
+    histos.add("hSparseRecLambda", "hSparseRecLambda", HistType::kTHnSparseF, runaxes2, true);
+    histos.add("hSparseRecAntiLambda", "hSparseRecAntiLambda", HistType::kTHnSparseF, runaxes2, true);
+
     ccdb->setURL(cfgCcdbParam.cfgURL);
     ccdbApi.init("http://alice-ccdb.cern.ch");
     ccdb->setCaching(true);
@@ -672,6 +679,10 @@ struct lambdapolsp {
     if (randGrp.doRandomPhi) {
       phiangle = randPhi.Uniform(0, 2 * TMath::Pi());
     }
+
+    auto ux = TMath::Cos(phiangle);
+    auto uy = TMath::Sin(phiangle);
+
     auto phiminuspsiC = GetPhiInRange(phiangle - psiZDCC);
     auto phiminuspsiA = GetPhiInRange(phiangle - psiZDCA);
     auto phiminuspsi = GetPhiInRange(phiangle - psiZDC);
@@ -680,6 +691,7 @@ struct lambdapolsp {
     auto PolC = TMath::Sin(phiminuspsiC);
     auto PolA = TMath::Sin(phiminuspsiA);
     auto Pol = TMath::Sin(phiminuspsi);
+    auto PolSP = uy * TMath::Cos(psiZDC) - ux * TMath::Sin(psiZDC);
 
     auto sinPhiStar = TMath::Sin(GetPhiInRange(phiangle));
     auto cosPhiStar = TMath::Cos(GetPhiInRange(phiangle));
@@ -694,6 +706,9 @@ struct lambdapolsp {
     auto PolAwgt = PolA / acvalue;
     auto PolCwgt = PolC / acvalue;
 
+    if (randGrp.useSP)
+      Polwgt = PolSP / acvalue;
+
     // Fill histograms using constructed names
     if (tag2) {
       if (needetaaxis) {
@@ -807,6 +822,7 @@ struct lambdapolsp {
   Filter dcaCutFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
 
   using EventCandidates = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::FV0Mults, aod::TPCMults, aod::CentFV0As, aod::CentFT0Ms, aod::CentFT0Cs, aod::CentFT0As, aod::SPCalibrationTables, aod::Mults>>;
+  using EventCandidatesMC = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs>>;
   using AllTrackCandidates = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullPi, aod::pidTPCFullPr, aod::pidTPCFullKa>>;
   using ResoV0s = aod::V0Datas;
 
@@ -895,6 +911,8 @@ struct lambdapolsp {
     if (!checkwithpub) {
       // histos.fill(HIST("hVtxZ"), collision.posZ());
       histos.fill(HIST("hpRes"), centrality, (TMath::Cos(GetPhiInRange(psiZDCA - psiZDCC))));
+      if (randGrp.useSP)
+        histos.fill(HIST("hpRes"), centrality, ((modqxZDCA * modqxZDCC) + (modqyZDCA * modqyZDCC)));
       // histos.fill(HIST("hpResSin"), centrality, (TMath::Sin(GetPhiInRange(psiZDCA - psiZDCC))));
       /*histos.fill(HIST("hpCosPsiA"), centrality, (TMath::Cos(GetPhiInRange(psiZDCA))));
       histos.fill(HIST("hpCosPsiC"), centrality, (TMath::Cos(GetPhiInRange(psiZDCC))));
@@ -1505,7 +1523,142 @@ struct lambdapolsp {
   }
   PROCESS_SWITCH(lambdapolsp, processDerivedData, "Process derived data", false);
 
+  using TrackMCTrueTable = aod::McParticles;
+  ROOT::Math::PxPyPzMVector lambdamc, antiLambdamc, lambdadummymc, antiLambdadummymc, protonmc, pionmc, antiProtonmc, antiPionmc;
+
+  void processMC(EventCandidatesMC::iterator const& collision, AllTrackCandidates const& tracks, TrackMCTrueTable const& GenParticles, ResoV0s const& V0s)
+  {
+    if (!collision.sel8()) {
+      return;
+    }
+    double centrality = -999.;
+    centrality = collision.centFT0C();
+
+    if (additionalEvSel && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) {
+      return;
+    }
+
+    if (additionalEvSel2 && (collision.trackOccupancyInTimeRange() > cfgMaxOccupancy || collision.trackOccupancyInTimeRange() < cfgMinOccupancy)) {
+      return;
+    }
+
+    if (additionalEvSel3 && (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder))) {
+      return;
+    }
+    if (additionalEvSel4 && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
+      return;
+    }
+
+    if (rctCut.requireRCTFlagChecker && !rctChecker(collision)) {
+      return;
+    }
+
+    histos.fill(HIST("hCentrality"), centrality);
+
+    for (const auto& v0 : V0s) {
+
+      auto postrack = v0.template posTrack_as<AllTrackCandidates>();
+      auto negtrack = v0.template negTrack_as<AllTrackCandidates>();
+
+      int LambdaTag = 0;
+      int aLambdaTag = 0;
+
+      const auto signpos = postrack.sign();
+      const auto signneg = negtrack.sign();
+
+      if (signpos < 0 || signneg > 0) {
+        continue;
+      }
+
+      if (isSelectedV0Daughter(v0, postrack, 0, 0) && isSelectedV0Daughter(v0, negtrack, 1, 0)) {
+        LambdaTag = 1;
+      }
+      if (isSelectedV0Daughter(v0, negtrack, 0, 1) && isSelectedV0Daughter(v0, postrack, 1, 1)) {
+        aLambdaTag = 1;
+      }
+
+      if (!LambdaTag && !aLambdaTag)
+        continue;
+
+      if (!SelectionV0(collision, v0)) {
+        continue;
+      }
+
+      if (LambdaTag) {
+        Proton = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPr);
+        AntiPion = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPi);
+        Lambdadummy = Proton + AntiPion;
+      }
+      if (aLambdaTag) {
+        AntiProton = ROOT::Math::PxPyPzMVector(v0.pxneg(), v0.pyneg(), v0.pzneg(), massPr);
+        Pion = ROOT::Math::PxPyPzMVector(v0.pxpos(), v0.pypos(), v0.pzpos(), massPi);
+        AntiLambdadummy = AntiProton + Pion;
+      }
+
+      if (shouldReject(LambdaTag, aLambdaTag, Lambdadummy, AntiLambdadummy)) {
+        continue;
+      }
+
+      if (TMath::Abs(v0.eta()) > 0.8)
+        continue;
+
+      if (LambdaTag) {
+        Lambda = Proton + AntiPion;
+        histos.fill(HIST("hSparseRecLambda"), v0.mLambda(), v0.pt(), centrality);
+      }
+      if (aLambdaTag) {
+        AntiLambda = AntiProton + Pion;
+        histos.fill(HIST("hSparseRecAntiLambda"), v0.mAntiLambda(), v0.pt(), centrality);
+      }
+    }
+
+    for (const auto& mcParticle : GenParticles) {
+      if (std::abs(mcParticle.pdgCode()) != PDG_t::kLambda0) {
+        continue;
+      }
+      if (std::abs(mcParticle.y()) > ConfV0Rap) {
+        continue;
+      }
+      auto pdg1 = mcParticle.pdgCode();
+      auto kDaughters = mcParticle.daughters_as<aod::McParticles>();
+      int daughsize = 2;
+      if (kDaughters.size() != daughsize) {
+        continue;
+      }
+      for (const auto& kCurrentDaughter : kDaughters) {
+
+        if (std::abs(kCurrentDaughter.pdgCode()) != PDG_t::kProton && std::abs(kCurrentDaughter.pdgCode()) != PDG_t::kPiPlus) {
+          continue;
+        }
+        if (kCurrentDaughter.pdgCode() == PDG_t::kProton) {
+          protonmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassProton);
+        }
+        if (kCurrentDaughter.pdgCode() == PDG_t::kPiMinus) {
+          antiPionmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassPionCharged);
+        }
+
+        if (kCurrentDaughter.pdgCode() == PDG_t::kProtonBar) {
+          antiProtonmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassProton);
+        }
+        if (kCurrentDaughter.pdgCode() == PDG_t::kPiPlus) {
+          pionmc = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassPionCharged);
+        }
+      }
+      if (pdg1 == PDG_t::kLambda0) {
+        lambdadummymc = protonmc + antiPionmc;
+        histos.fill(HIST("hSparseGenLambda"), lambdadummymc.M(), lambdadummymc.Pt(), centrality);
+      }
+
+      if (pdg1 == PDG_t::kLambda0Bar) {
+        antiLambdadummymc = antiProtonmc + pionmc;
+        histos.fill(HIST("hSparseGenAntiLambda"), antiLambdadummymc.M(), antiLambdadummymc.Pt(), centrality);
+      }
+    }
+  }
+  PROCESS_SWITCH(lambdapolsp, processMC, "Process MC", false);
+
   // Processing Event Mixing
+  /*
   using BinningType = ColumnBinningPolicy<aod::collision::PosZ, aod::cent::CentFT0C>;
   BinningType colBinning{{meGrp.axisVertex, meGrp.axisMultiplicityClass}, true};
   Preslice<v0Candidates> tracksPerCollisionV0Mixed = o2::aod::v0data::straCollisionId; // for derived data only
@@ -1644,6 +1797,7 @@ struct lambdapolsp {
   }
   PROCESS_SWITCH(lambdapolsp, processDerivedDataMixed, "Process mixed event using derived data", false);
 
+
   void processDerivedDataMixed2(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraStamps, aod::StraZDCSP> const& collisions, v0Candidates const& V0s, dauTracks const&)
   {
     TRandom3 randGen(0);
@@ -1795,150 +1949,9 @@ struct lambdapolsp {
     }
   }
   PROCESS_SWITCH(lambdapolsp, processDerivedDataMixed2, "Process mixed event2 using derived data", false);
-
-  void processDerivedDataMixedFIFO(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraStamps, aod::StraZDCSP> const& collisions, v0Candidates const& V0s, dauTracks const&)
-  {
-
-    auto nBins = colBinning.getAllBinsCount();
-    std::vector<std::deque<int>> eventPools(nBins); // Pool per bin holding just event indices
-
-    for (auto& collision1 : collisions) {
-
-      if (!collision1.sel8()) {
-        continue;
-      }
-      if (!collision1.triggereventsp()) { // provided by StraZDCSP
-        continue;
-      }
-      if (rctCut.requireRCTFlagChecker && !rctChecker(collision1)) {
-        continue;
-      }
-
-      if (additionalEvSel && (!collision1.selection_bit(aod::evsel::kNoSameBunchPileup) || !collision1.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) {
-        continue;
-      }
-      if (additionalEvSel2 && (collision1.trackOccupancyInTimeRange() > cfgMaxOccupancy || collision1.trackOccupancyInTimeRange() < cfgMinOccupancy)) {
-        continue;
-      }
-      if (additionalEvSel3 && (!collision1.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision1.selection_bit(aod::evsel::kNoITSROFrameBorder))) {
-        continue;
-      }
-      if (additionalEvSel4 && !collision1.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
-        continue;
-      }
-
-      int bin = colBinning.getBin(std::make_tuple(collision1.posZ(), collision1.centFT0C()));
-      auto groupV0_evt1 = V0s.sliceBy(tracksPerCollisionV0Mixed, collision1.index());
-      float centrality = collision1.centFT0C();
-      auto qxZDCA = collision1.qxZDCA();
-      auto qxZDCC = collision1.qxZDCC();
-      auto qyZDCA = collision1.qyZDCA();
-      auto qyZDCC = collision1.qyZDCC();
-      auto psiZDCC = collision1.psiZDCC();
-      auto psiZDCA = collision1.psiZDCA();
-      double modqxZDCA;
-      double modqyZDCA;
-      double modqxZDCC;
-      double modqyZDCC;
-
-      if (bin < 0)
-        continue;
-      modqxZDCA = TMath::Sqrt((qxZDCA * qxZDCA) + (qyZDCA * qyZDCA)) * TMath::Cos(psiZDCA);
-      modqyZDCA = TMath::Sqrt((qxZDCA * qxZDCA) + (qyZDCA * qyZDCA)) * TMath::Sin(psiZDCA);
-      modqxZDCC = TMath::Sqrt((qxZDCC * qxZDCC) + (qyZDCC * qyZDCC)) * TMath::Cos(psiZDCC);
-      modqyZDCC = TMath::Sqrt((qxZDCC * qxZDCC) + (qyZDCC * qyZDCC)) * TMath::Sin(psiZDCC);
-
-      auto psiZDC = TMath::ATan2((modqyZDCC - modqyZDCA), (modqxZDCC - modqxZDCA)); // full event plane from collision
-
-      histos.fill(HIST("hCentrality"), centrality);
-      histos.fill(HIST("hpRes"), centrality, (TMath::Cos(GetPhiInRange(psiZDCA - psiZDCC))));
-
-      // For deduplication of (v0_evt1, v0_evt2) pairs per mixed event
-      std::unordered_map<int, std::set<std::pair<int, int>>> seenMap;
-
-      // Loop over Λ candidates in collision1 (keep psi from here)
-
-      for (auto& v0_evt1 : groupV0_evt1) {
-        if (!SelectionV0(collision1, v0_evt1))
-          continue;
-        bool LambdaTag1 = isCompatible(v0_evt1, 0);
-        bool aLambdaTag1 = isCompatible(v0_evt1, 1);
-        ROOT::Math::PxPyPzMVector proton1, pion1, antiproton1, antipion1, LambdaTag1dummy, AntiLambdaTag1dummy;
-        if (LambdaTag1) {
-          proton1 = {v0_evt1.pxpos(), v0_evt1.pypos(), v0_evt1.pzpos(), massPr};
-          antipion1 = {v0_evt1.pxneg(), v0_evt1.pyneg(), v0_evt1.pzneg(), massPi};
-          LambdaTag1dummy = proton1 + antipion1;
-        }
-        if (aLambdaTag1) {
-          antiproton1 = {v0_evt1.pxneg(), v0_evt1.pyneg(), v0_evt1.pzneg(), massPr};
-          pion1 = {v0_evt1.pxpos(), v0_evt1.pypos(), v0_evt1.pzpos(), massPi};
-          AntiLambdaTag1dummy = antiproton1 + pion1;
-        }
-        if (shouldReject(LambdaTag1, aLambdaTag1, LambdaTag1dummy, AntiLambdaTag1dummy)) {
-          continue;
-        }
-        if (TMath::Abs(v0_evt1.eta()) > 0.8)
-          continue;
-
-        // Loop over all FIFO pool events (mixed events) for this centrality bin
-        int nMixedEvents = 0;
-        for (auto it = eventPools[bin].rbegin(); it != eventPools[bin].rend() && nMixedEvents < meGrp.nMix; ++it, ++nMixedEvents) {
-          int collision2idx = *it;
-          if (collision1.index() == collision2idx)
-            continue;
-          auto groupV0_evt2 = V0s.sliceBy(tracksPerCollisionV0Mixed, collision2idx);
-
-          // Now loop over Λ candidates in collision2 to randomize proton phi* (randomize decay angle)
-          for (auto& v0_evt2 : groupV0_evt2) {
-            if (!SelectionV0(collision1, v0_evt2))
-              continue;
-            bool LambdaTag2 = isCompatible(v0_evt2, 0);
-            bool aLambdaTag2 = isCompatible(v0_evt2, 1);
-            if (!LambdaTag2 && !aLambdaTag2)
-              continue;
-
-            // Deduplicate (v0_evt1, v0_evt2) pairs per collision2idx
-            auto key = std::make_pair(v0_evt1.index(), v0_evt2.index());
-            if (!seenMap[collision2idx].insert(key).second)
-              continue;
-
-            ROOT::Math::PxPyPzMVector proton_mix, antiproton_mix, pion_mix, antipion_mix, LambdaTag2dummy, AntiLambdaTag2dummy;
-            if (LambdaTag2) {
-              proton_mix = {v0_evt2.pxpos(), v0_evt2.pypos(), v0_evt2.pzpos(), massPr};
-              antipion_mix = {v0_evt2.pxneg(), v0_evt2.pyneg(), v0_evt2.pzneg(), massPi};
-              LambdaTag2dummy = proton_mix + antipion_mix;
-            }
-            if (aLambdaTag2) {
-              antiproton_mix = {v0_evt2.pxneg(), v0_evt2.pyneg(), v0_evt2.pzneg(), massPr};
-              pion_mix = {v0_evt2.pxpos(), v0_evt2.pypos(), v0_evt2.pzpos(), massPi};
-              AntiLambdaTag2dummy = antiproton_mix + pion_mix;
-            }
-            if (shouldReject(LambdaTag2, aLambdaTag2, LambdaTag2dummy, AntiLambdaTag2dummy)) {
-              continue;
-            }
-            if (TMath::Abs(v0_evt2.eta()) > 0.8)
-              continue;
-            if (LambdaTag1) {
-              double acvalue = 1.0;
-              fillHistograms(1, 0, LambdaTag1dummy, proton_mix, psiZDCC, psiZDCA, psiZDC, centrality, v0_evt1.mLambda(), v0_evt1.pt(), v0_evt1.eta(), acvalue, 1.0);
-            }
-            if (aLambdaTag1) {
-              double acvalue = 1.0;
-              fillHistograms(0, 1, AntiLambdaTag1dummy, antiproton_mix, psiZDCC, psiZDCA, psiZDC, centrality, v0_evt1.mAntiLambda(), v0_evt1.pt(), v0_evt1.eta(), acvalue, 1.0);
-            }
-          }
-        }
-      }
-      // After processing all mixes, add current event V0s to pool for future mixing
-      eventPools[bin].push_back(collision1.index());
-      // Keep only N last events in FIFO queue
-      if (static_cast<int>(eventPools[bin].size()) > meGrp.nMix) {
-        eventPools[bin].pop_front();
-      }
-    }
-  }
-  PROCESS_SWITCH(lambdapolsp, processDerivedDataMixedFIFO, "Process mixed event using derived data with FIFO method", false);
+ */
 };
+
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   return WorkflowSpec{