Merge branch 'master' into separate

Author: dsekihat

PWGCF/EbyEFluctuations/Tasks/partNumFluc.cxx

@@ -1634,16 +1634,33 @@ struct LambdaSpinPolarization {
     float dphi = RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PI);
     float dR = std::sqrt(drap * drap + dphi * dphi);
 
-    // Get Lambda-Proton four-momentum
     std::array<float, 4> l1 = {p1.px(), p1.py(), p1.pz(), MassLambda0};
     std::array<float, 4> l2 = {p2.px(), p2.py(), p2.pz(), MassLambda0};
+
+    std::array<float, 4> llpair = {
+      l1[0] + l2[0],
+      l1[1] + l2[1],
+      l1[2] + l2[2],
+      l1[3] + l2[3]};
+
     std::array<float, 4> pr1 = {p1.prPx(), p1.prPy(), p1.prPz(), MassProton};
     std::array<float, 4> pr2 = {p2.prPx(), p2.prPy(), p2.prPz(), MassProton};
-    std::array<float, 3> v1, v2;
-    getBoostVector(l1, v1, cInvBoostFlag);
-    getBoostVector(l2, v2, cInvBoostFlag);
-    boost(pr1, v1);
-    boost(pr2, v2);
+
+    std::array<float, 3> vPair;
+    getBoostVector(llpair, vPair, cInvBoostFlag);
+
+    boost(l1, vPair);
+    boost(l2, vPair);
+
+    boost(pr1, vPair);
+    boost(pr2, vPair);
+
+    std::array<float, 3> v1_pair, v2_pair;
+    getBoostVector(l1, v1_pair, cInvBoostFlag);
+    getBoostVector(l2, v2_pair, cInvBoostFlag);
+
+    boost(pr1, v1_pair);
+    boost(pr2, v2_pair);
 
     std::array<float, 3> pr1tv = {pr1[0], pr1[1], pr1[2]};
     std::array<float, 3> pr2tv = {pr2[0], pr2[1], pr2[2]};
@@ -1676,7 +1693,7 @@ struct LambdaSpinPolarization {
           continue;
         }
 
-        // Kinematic matching with named constants
+        // Kinematic matching
         float deltaPt = std::abs(trk_1.pt() - trk_2.pt());
         float deltaPhi = std::abs(RecoDecay::constrainAngle(trk_1.phi() - trk_2.phi(), -PI));
         float deltaRap = std::abs(trk_1.rap() - trk_2.rap());
@@ -1718,7 +1735,7 @@ struct LambdaSpinPolarization {
     analyzePairs<kAntiLambdaAntiLambda, true>(antiLambdaTracks, antiLambdaTracks);
   }
 
-  PROCESS_SWITCH(LambdaSpinPolarization, processDataReco, "Process for Data and MCReco", false);
+  PROCESS_SWITCH(LambdaSpinPolarization, processDataReco, "Process for Data and MCReco", true);
 
   struct GetMultiplicity {
     float operator()(auto const& col) const
@@ -1766,7 +1783,7 @@ struct LambdaSpinPolarization {
     }
   }
 
-  PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixed, "Process for Data and MCReco for Mixed events", true);
+  PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixed, "Process for Data and MCReco for Mixed events", false);
 
   void processDataRecoMixEvent(LambdaCollisions::iterator const& collision, LambdaTracks const& tracks)
   {

PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx

@@ -423,8 +423,66 @@ std::tuple<float, float, float, float, float, int> VarManager::BimodalityCoeffic
   // trim the distribution if requested, by requiring a minimum of "trim" counts in each bin
   if (trim > 0) {
     for (int i = 0; i < nBins; ++i) {
+      // if the count in the bin is less than the trim value, set it to zero
       if (counts[i] < trim) {
-        counts[i] = 0;
+        // set the count to zero only if this is an isolated bin,
+        // i.e. if this count belongs to a peak, we want to keep it even if it has counts below the trim limit, as long as the whole peak is not below the trim limit
+        // check the previous bins until we find an empty bin or we reach the beginning of the histogram
+        int localPeakCount = counts[i];
+        for (int j = i - 1; j >= 0; --j) {
+          if (counts[j] == 0) {
+            break;
+          }
+          localPeakCount += counts[j];
+        }
+        // check the next bins until we find an empty bin or we reach the end of the histogram
+        for (int j = i + 1; j < nBins; ++j) {
+          if (counts[j] == 0) {
+            break;
+          }
+          localPeakCount += counts[j];
+        }
+        if (localPeakCount < trim) {
+          counts[i] = 0;
+        }
+      }
+    }
+  }
+  if (trim < 0) {
+    // if trim is negative, then we remove all counts belonging to local peaks with an integrated count below 1/abs(trim)
+    for (int i = 0; i < nBins; ++i) {
+      if (counts[i] == 0) {
+        continue; // skip empty bins
+      }
+      // check the previous bins until we find an empty bin or we reach the beginning of the histogram
+      int localPeakCount = counts[i];
+      for (int j = i - 1; j >= 0; --j) {
+        if (counts[j] == 0) {
+          break;
+        }
+        localPeakCount += counts[j];
+      }
+      // check the next bins until we find an empty bin or we reach the end of the histogram
+      for (int j = i + 1; j < nBins; ++j) {
+        if (counts[j] == 0) {
+          break;
+        }
+        localPeakCount += counts[j];
+      }
+      if (localPeakCount < (1.0 / std::abs(trim)) * data.size()) {
+        // set all bins belonging to this local peak to zero
+        for (int j = i; j >= 0; --j) {
+          if (counts[j] == 0) {
+            break;
+          }
+          counts[j] = 0;
+        }
+        for (int j = i + 1; j < nBins; ++j) {
+          if (counts[j] == 0) {
+            break;
+          }
+          counts[j] = 0;
+        }
       }
     }
   }
@@ -453,7 +511,7 @@ std::tuple<float, float, float, float, float, int> VarManager::BimodalityCoeffic
   }
 
   if (totalCounts == 0) {
-    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0, -1);
+    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0, nPeaks);
   }
   mean /= totalCounts;
 
@@ -477,7 +535,7 @@ std::tuple<float, float, float, float, float, int> VarManager::BimodalityCoeffic
   m4 /= totalCounts;
 
   if (m2 == 0.0) {
-    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0, -1);
+    return std::make_tuple(-1.0, -1.0, -1.0, -1.0, -1.0, nPeaks);
   }
 
   float stddev = std::sqrt(m2);

PWGDQ/Core/HistogramsLibrary.cxx

@@ -2364,6 +2364,36 @@ void VarManager::FillEvent(T const& event, float* values)
     values[kBGFDDCpf] = event.bgFDDCpf();
   }
 
+  // check if the template T type includes the O2remerge table and if so, fill the corresponding variables
+  if constexpr (requires { event.dcazBimodalityCoefficient(); }) {
+    values[kDCAzBimodalityCoefficient] = event.dcazBimodalityCoefficient();
+    values[kDCAzBimodalityCoefficientBinned] = event.dcazBimodalityCoefficientBinned();
+    values[kDCAzBimodalityCoefficientBinnedTrimmed1] = event.dcazBimodalityCoefficientBinnedTrimmed1();
+    values[kDCAzBimodalityCoefficientBinnedTrimmed2] = event.dcazBimodalityCoefficientBinnedTrimmed2();
+    values[kDCAzBimodalityCoefficientBinnedTrimmed3] = event.dcazBimodalityCoefficientBinnedTrimmed3();
+    values[kDCAzMean] = event.dcazMean();
+    values[kDCAzMeanBinnedTrimmed1] = event.dcazMeanBinnedTrimmed1();
+    values[kDCAzMeanBinnedTrimmed2] = event.dcazMeanBinnedTrimmed2();
+    values[kDCAzMeanBinnedTrimmed3] = event.dcazMeanBinnedTrimmed3();
+    values[kDCAzRMS] = event.dcazRMS();
+    values[kDCAzRMSBinnedTrimmed1] = event.dcazRMSBinnedTrimmed1();
+    values[kDCAzRMSBinnedTrimmed2] = event.dcazRMSBinnedTrimmed2();
+    values[kDCAzRMSBinnedTrimmed3] = event.dcazRMSBinnedTrimmed3();
+    values[kDCAzSkewness] = event.dcazSkewness();
+    values[kDCAzKurtosis] = event.dcazKurtosis();
+    values[kDCAzNPeaks] = event.dcazNPeaks();
+    values[kDCAzNPeaksTrimmed1] = event.dcazNPeaksTrimmed1();
+    values[kDCAzNPeaksTrimmed2] = event.dcazNPeaksTrimmed2();
+    values[kDCAzNPeaksTrimmed3] = event.dcazNPeaksTrimmed3();
+    values[kDCAzFracAbove100um] = event.dcazFracAbove100um();
+    values[kDCAzFracAbove200um] = event.dcazFracAbove200um();
+    values[kDCAzFracAbove500um] = event.dcazFracAbove500um();
+    values[kDCAzFracAbove1mm] = event.dcazFracAbove1mm();
+    values[kDCAzFracAbove2mm] = event.dcazFracAbove2mm();
+    values[kDCAzFracAbove5mm] = event.dcazFracAbove5mm();
+    values[kDCAzFracAbove10mm] = event.dcazFracAbove10mm();
+  }
+
   // FillEventDerived(values);
 }
 
@@ -2414,51 +2444,50 @@ void VarManager::FillEventTracks(T const& tracks, float* values)
   auto [bimodalityCoefficientBin, meanBin, stddevBin, skewnessBin, kurtosisBin, nPeaksBin] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005);
   if (stddevBin > -1.0) {
     values[kDCAzBimodalityCoefficientBinned] = bimodalityCoefficientBin;
-    values[kDCAzNPeaks] = nPeaksBin;
   } else {
     values[kDCAzBimodalityCoefficientBinned] = -9999.0;
-    values[kDCAzNPeaks] = -9999.0;
   }
+  values[kDCAzNPeaks] = nPeaksBin;
   cout << "Bimodality coefficient binned: " << bimodalityCoefficientBin << ", mean: " << mean << ", stddev: " << stddev << ", skewness: " << skewness << ", kurtosis: " << kurtosis << ", nPeaks: " << nPeaksBin << endl;
-  // compute the binned bimodality coefficient and related statistics with a bin width of 50um and trimming of 1, 2, 3 entries per bin
-  auto [bimodalityCoefficientBinTrimmed1, meanBinTrimmed1, stddevBinTrimmed1, skewnessBinTrimmed1, kurtosisBinTrimmed1, nPeaksBinTrimmed1] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 2);
+  // compute the binned bimodality coefficient and related statistics with a bin width of 50um and different trimming versions
+  // more then 3 counts per bin
+  auto [bimodalityCoefficientBinTrimmed1, meanBinTrimmed1, stddevBinTrimmed1, skewnessBinTrimmed1, kurtosisBinTrimmed1, nPeaksBinTrimmed1] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 4);
   if (stddevBinTrimmed1 > -1.0) {
     values[kDCAzBimodalityCoefficientBinnedTrimmed1] = bimodalityCoefficientBinTrimmed1;
     values[kDCAzMeanBinnedTrimmed1] = meanBinTrimmed1;
     values[kDCAzRMSBinnedTrimmed1] = stddevBinTrimmed1;
-    values[kDCAzNPeaksTrimmed1] = nPeaksBinTrimmed1;
   } else {
     values[kDCAzBimodalityCoefficientBinnedTrimmed1] = -9999.0;
     values[kDCAzMeanBinnedTrimmed1] = -9999.0;
     values[kDCAzRMSBinnedTrimmed1] = -9999.0;
-    values[kDCAzNPeaksTrimmed1] = -9999.0;
   }
+  values[kDCAzNPeaksTrimmed1] = nPeaksBinTrimmed1;
   cout << "Bimodality coefficient (trimmed 1): " << bimodalityCoefficientBinTrimmed1 << ", mean: " << meanBinTrimmed1 << ", stddev: " << stddevBinTrimmed1 << ", skewness: " << skewnessBinTrimmed1 << ", kurtosis: " << kurtosisBinTrimmed1 << ", nPeaks: " << nPeaksBinTrimmed1 << endl;
-  auto [bimodalityCoefficientBinTrimmed2, meanBinTrimmed2, stddevBinTrimmed2, skewnessBinTrimmed2, kurtosisBinTrimmed2, nPeaksBinTrimmed2] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 3);
+  // more than 3% of the entries per peak
+  auto [bimodalityCoefficientBinTrimmed2, meanBinTrimmed2, stddevBinTrimmed2, skewnessBinTrimmed2, kurtosisBinTrimmed2, nPeaksBinTrimmed2] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, -100);
   if (stddevBinTrimmed2 > -1.0) {
     values[kDCAzBimodalityCoefficientBinnedTrimmed2] = bimodalityCoefficientBinTrimmed2;
     values[kDCAzMeanBinnedTrimmed2] = meanBinTrimmed2;
     values[kDCAzRMSBinnedTrimmed2] = stddevBinTrimmed2;
-    values[kDCAzNPeaksTrimmed2] = nPeaksBinTrimmed2;
   } else {
     values[kDCAzBimodalityCoefficientBinnedTrimmed2] = -9999.0;
     values[kDCAzMeanBinnedTrimmed2] = -9999.0;
     values[kDCAzRMSBinnedTrimmed2] = -9999.0;
-    values[kDCAzNPeaksTrimmed2] = -9999.0;
   }
+  values[kDCAzNPeaksTrimmed2] = nPeaksBinTrimmed2;
   cout << "Bimodality coefficient (trimmed 2): " << bimodalityCoefficientBinTrimmed2 << ", mean: " << meanBinTrimmed2 << ", stddev: " << stddevBinTrimmed2 << ", skewness: " << skewnessBinTrimmed2 << ", kurtosis: " << kurtosisBinTrimmed2 << ", nPeaks: " << nPeaksBinTrimmed2 << endl;
-  auto [bimodalityCoefficientBinTrimmed3, meanBinTrimmed3, stddevBinTrimmed3, skewnessBinTrimmed3, kurtosisBinTrimmed3, nPeaksBinTrimmed3] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, 4);
+  // more than 5% of the entries per peak
+  auto [bimodalityCoefficientBinTrimmed3, meanBinTrimmed3, stddevBinTrimmed3, skewnessBinTrimmed3, kurtosisBinTrimmed3, nPeaksBinTrimmed3] = BimodalityCoefficientAndNPeaks(dcazValues, 0.005, -20);
   if (stddevBinTrimmed3 > -1.0) {
     values[kDCAzBimodalityCoefficientBinnedTrimmed3] = bimodalityCoefficientBinTrimmed3;
     values[kDCAzMeanBinnedTrimmed3] = meanBinTrimmed3;
     values[kDCAzRMSBinnedTrimmed3] = stddevBinTrimmed3;
-    values[kDCAzNPeaksTrimmed3] = nPeaksBinTrimmed3;
   } else {
     values[kDCAzBimodalityCoefficientBinnedTrimmed3] = -9999.0;
     values[kDCAzMeanBinnedTrimmed3] = -9999.0;
     values[kDCAzRMSBinnedTrimmed3] = -9999.0;
-    values[kDCAzNPeaksTrimmed3] = -9999.0;
   }
+  values[kDCAzNPeaksTrimmed3] = nPeaksBinTrimmed3;
   cout << "Bimodality coefficient (trimmed 3): " << bimodalityCoefficientBinTrimmed3 << ", mean: " << meanBinTrimmed3 << ", stddev: " << stddevBinTrimmed3 << ", skewness: " << skewnessBinTrimmed3 << ", kurtosis: " << kurtosisBinTrimmed3 << ", nPeaks: " << nPeaksBinTrimmed3 << endl;
 
   // compute fraction of tracks with |DCAz| > 100um, 200um, 500um, 1mm, 2mm, 5mm, 10mm

PWGDQ/Core/VarManager.cxx

@@ -194,8 +194,8 @@ DECLARE_SOA_TABLE(JPsieeCandidates, "AOD", "DQPSEUDOPROPER", dqanalysisflags::Ma
 } // namespace o2::aod
 
 // Declarations of various short names
-using MyEvents = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended>;
-using MyEventsMultExtra = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedEventsMultPV, aod::ReducedEventsMultAll>;
+using MyEvents = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedEventsMultPV, aod::ReducedEventsMultAll>;
+using MyEventsMultExtra = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedEventsMultPV, aod::ReducedEventsMultAll, aod::ReducedEventsQvectorCentr, aod::ReducedEventsMergingTable>;
 using MyEventsMultExtraQVector = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedEventsMultPV, aod::ReducedEventsMultAll, aod::ReducedEventsQvectorCentr, aod::ReducedEventsQvectorCentrExtra>;
 using MyEventsZdc = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedZdcs>;
 using MyEventsMultExtraZdc = soa::Join<aod::ReducedEvents, aod::ReducedEventsExtended, aod::ReducedEventsMultPV, aod::ReducedEventsMultAll, aod::ReducedZdcs>;

PWGDQ/Core/VarManager.h

@@ -41,6 +41,20 @@ static constexpr bool HasPrimaries = !std::is_same_v<TPrimaries, std::nullptr_t>
 template <typename TSecondaries>
 static constexpr bool HasSecondaries = !std::is_same_v<TSecondaries, std::nullptr_t>;
 
+template <typename T>
+concept IsNonLinIterator = o2::soa::is_iterator<T> && requires(T t) {
+  // Check that the *elements* of the container have the required methods:
+  { t.corrE() } -> std::same_as<float>;
+  { t.corrPt() } -> std::same_as<float>;
+};
+
+template <typename T>
+concept IsNonLinContainer = o2::soa::is_table<T> && requires(T t) {
+  // Check that the *elements* of the container have the required methods:
+  { t.begin().corrE() } -> std::same_as<float>;
+  { t.begin().corrPt() } -> std::same_as<float>;
+};
+
 template <typename T>
 concept IsTrackIterator = o2::soa::is_iterator<T> && requires(T t) {
   // Check that the *elements* of the container have the required methods:
@@ -114,7 +128,7 @@ class EMCPhotonCut : public TNamed
 
   static const char* mCutNames[static_cast<int>(EMCPhotonCuts::kNCuts)];
 
-  constexpr auto getClusterId(o2::soa::is_iterator auto const& t) const
+  static constexpr auto getClusterId(o2::soa::is_iterator auto const& t)
   {
     if constexpr (requires { t.emEmcClusterId(); }) {
       return t.emEmcClusterId();
@@ -384,7 +398,11 @@ class EMCPhotonCut : public TNamed
         return cluster.definition() == mDefinition;
 
       case EMCPhotonCuts::kEnergy:
-        return cluster.e() > mMinE;
+        if constexpr (IsNonLinIterator<std::decay_t<decltype(cluster)>>) {
+          return cluster.corrE() > mMinE;
+        } else {
+          return cluster.e() > mMinE;
+        }
 
       case EMCPhotonCuts::kNCell:
         return cluster.nCells() >= mMinNCell;
@@ -478,7 +496,11 @@ class EMCPhotonCut : public TNamed
         return cluster.definition() == mDefinition;
 
       case EMCPhotonCuts::kEnergy:
-        return cluster.e() > mMinE;
+        if constexpr (IsNonLinIterator<Cluster>) {
+          return cluster.corrE() > mMinE;
+        } else {
+          return cluster.e() > mMinE;
+        }
 
       case EMCPhotonCuts::kNCell:
         return cluster.nCells() >= mMinNCell;
@@ -496,7 +518,9 @@ class EMCPhotonCut : public TNamed
             auto dPhi = std::fabs(emcmatchedtrack.deltaPhi());
             auto trackpt = emcmatchedtrack.trackPt();
             auto trackp = emcmatchedtrack.trackP();
-            bool result = (dEta > GetTrackMatchingEta(trackpt)) || (dPhi > GetTrackMatchingPhi(trackpt)) || (cluster.e() / trackp >= mMinEoverP);
+            bool result = (dEta > GetTrackMatchingEta(trackpt)) ||
+                          (dPhi > GetTrackMatchingPhi(trackpt)) ||
+                          (cluster.e() / trackp >= mMinEoverP);
             if (!result) {
               return false;
             }
@@ -506,7 +530,7 @@ class EMCPhotonCut : public TNamed
           auto dPhis = cluster.deltaPhi();   // std:vector<float>
           auto trackspt = cluster.trackpt(); // std:vector<float>
           auto tracksp = cluster.trackp();   // std:vector<float>
-          int ntrack = tracksp.size();
+          int ntrack = trackspt.size();
           for (int itr = 0; itr < ntrack; itr++) {
             float dEta = std::fabs(dEtas[itr]);
             float dPhi = std::fabs(dPhis[itr]);
@@ -526,8 +550,8 @@ class EMCPhotonCut : public TNamed
             auto dEta = std::fabs(emcmatchedtrack.deltaEta());
             auto dPhi = std::fabs(emcmatchedtrack.deltaPhi());
             auto trackpt = emcmatchedtrack.trackPt();
-            auto trackp = emcmatchedtrack.trackP();
-            bool result = (dEta > GetSecTrackMatchingEta(trackpt)) || (dPhi > GetSecTrackMatchingPhi(trackpt)) || (cluster.e() / trackp >= mMinEoverP);
+            bool result = (dEta > GetSecTrackMatchingEta(trackpt)) ||
+                          (dPhi > GetSecTrackMatchingPhi(trackpt));
             if (!result) {
               return false;
             }
@@ -536,8 +560,7 @@ class EMCPhotonCut : public TNamed
           auto dEtas = cluster.deltaEtaSec();   // std:vector<float>
           auto dPhis = cluster.deltaPhiSec();   // std:vector<float>
           auto trackspt = cluster.trackptSec(); // std:vector<float>
-          auto tracksp = cluster.trackpSec();   // std:vector<float>
-          int ntrack = tracksp.size();
+          int ntrack = trackspt.size();
           for (int itr = 0; itr < ntrack; itr++) {
             float dEta = std::fabs(dEtas[itr]);
             float dPhi = std::fabs(dPhis[itr]);