[PWGLF] added option to equalise width of Q distributions vs centrality and a possible downscale factor for pions (#17482)

Co-authored-by: Prottay Das <prottay@alipap1.cern.ch>

Author: prottayCMT

PWGLF/TableProducer/Common/spvector.cxx

@@ -247,6 +247,22 @@ struct spvector {
     Configurable<std::string> confRecentereVzTimeSp3{"confRecentereVzTimeSp3", "Users/p/prottay/My/Object/Testingwithsparse/NewPbPbpass4_17092024/recenter", "TH3F path for (Vz,time) recentering, cycle 3"};
     Configurable<std::string> confRecentereVzTimeSp4{"confRecentereVzTimeSp4", "Users/p/prottay/My/Object/Testingwithsparse/NewPbPbpass4_17092024/recenter", "TH3F path for (Vz,time) recentering, cycle 4"};
     Configurable<std::string> confRecentereVzTimeSp5{"confRecentereVzTimeSp5", "Users/p/prottay/My/Object/Testingwithsparse/NewPbPbpass4_17092024/recenter", "TH3F path for (Vz,time) recentering, cycle 5"};
+
+    Configurable<bool> useWidthEqualization{
+      "useWidthEqualization",
+      false,
+      "Apply Q-vector width equalization"};
+
+    Configurable<std::string> confWidthSigma{
+      "confWidthSigma",
+      "Users/p/prottay/My/Object/ZDC/widthEq",
+      "TH2F path for Q-vector width equalization"};
+
+    Configurable<std::string> confWidthMean{
+      "confWidthMean",
+      "Users/p/prottay/My/Object/ZDC/widthEq",
+      "TH2F path for Q-vector mean equalization"};
+
   } confignewpro;
 
   struct : ConfigurableGroup {
@@ -503,6 +519,59 @@ struct spvector {
   std::array<TH3F*, 6> hrecenterecentvySpA{};
   std::array<TH3F*, 6> hrecenterecentvzSpA{};
   std::array<TH3F*, 5> hrecenterevztimeSpA{};
+  TH2F* hWidthMean = nullptr;
+  TH2F* hWidthSigma = nullptr;
+
+  bool CorrectWidthEqualization(TH2F* hMean,
+                                TH2F* hWidth,
+                                auto centrality,
+                                auto& qxZDCA,
+                                auto& qyZDCA,
+                                auto& qxZDCC,
+                                auto& qyZDCC)
+  {
+
+    if (!hMean || !hWidth) {
+      return false;
+    }
+
+    const double meanQxA = hMean->GetBinContent(
+      hMean->FindBin(centrality + 0.00001, 0.5));
+
+    const double meanQyA = hMean->GetBinContent(
+      hMean->FindBin(centrality + 0.00001, 1.5));
+
+    const double meanQxC = hMean->GetBinContent(
+      hMean->FindBin(centrality + 0.00001, 2.5));
+
+    const double meanQyC = hMean->GetBinContent(
+      hMean->FindBin(centrality + 0.00001, 3.5));
+
+    const double sigmaQxA = hWidth->GetBinContent(
+      hWidth->FindBin(centrality + 0.00001, 0.5));
+
+    const double sigmaQyA = hWidth->GetBinContent(
+      hWidth->FindBin(centrality + 0.00001, 1.5));
+
+    const double sigmaQxC = hWidth->GetBinContent(
+      hWidth->FindBin(centrality + 0.00001, 2.5));
+
+    const double sigmaQyC = hWidth->GetBinContent(
+      hWidth->FindBin(centrality + 0.00001, 3.5));
+
+    if (sigmaQxA <= 0.0 || sigmaQyA <= 0.0 ||
+        sigmaQxC <= 0.0 || sigmaQyC <= 0.0) {
+      return false;
+    }
+
+    qxZDCA = (qxZDCA - meanQxA) / sigmaQxA;
+    qyZDCA = (qyZDCA - meanQyA) / sigmaQyA;
+
+    qxZDCC = (qxZDCC - meanQxC) / sigmaQxC;
+    qyZDCC = (qyZDCC - meanQyC) / sigmaQyC;
+
+    return true;
+  }
 
   bool CorrectfineCent(TH2F* hrecenterecentSp,
                        auto centrality,
@@ -1404,6 +1473,19 @@ struct spvector {
             LOGF(fatal, "Cannot apply fine (cent,Vz) recentering, stage 6");
           }
         }
+
+        // -------------------- Width equalization --------------------
+        if (confignewpro.useWidthEqualization) {
+          if (currentRunNumber != lastRunNumber) {
+
+            hWidthMean = ccdb->getForTimeStamp<TH2F>(confignewpro.confWidthMean.value, bc.timestamp());
+            hWidthSigma = ccdb->getForTimeStamp<TH2F>(confignewpro.confWidthSigma.value, bc.timestamp());
+          }
+
+          if (!CorrectWidthEqualization(hWidthMean, hWidthSigma, centrality, qxZDCA, qyZDCA, qxZDCC, qyZDCC)) {
+            LOGF(fatal, "Cannot apply Q-vector width equalization");
+          }
+        }
       }
 
       histos.fill(HIST("hpQxZDCAvstime"), timeMin, qxZDCA);

PWGLF/TableProducer/Resonances/cksspinalignment.cxx

@@ -120,6 +120,8 @@ struct cksspinalignment {
   Configurable<float> cfgK0sMassMin{"cfgK0sMassMin", 0.45f, "Minimum K0s invariant mass"};
   Configurable<float> cfgK0sMassMax{"cfgK0sMassMax", 0.55f, "Maximum K0s invariant mass"};
 
+  Configurable<float> cfgKeepFraction{"cfgKeepFraction", 0.4f, "Fraction of events to keep"};
+
   HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
   RCTFlagsChecker rctChecker;
 
@@ -136,6 +138,26 @@ struct cksspinalignment {
     histos.add("hNStoredPions", "hNStoredPions;N_{#pi};Events", kTH1F, {{500, 0.0f, 500.0f}});
   }
 
+  bool keepEvent(uint64_t eventIndex) const
+  { /*
+     if (cfgKeepFraction >= 1.0f) {
+       return true;
+     }
+    */
+    // SplitMix64 pseudo-random hash
+    uint64_t x = eventIndex + 0x9e3779b97f4a7c15ULL;
+    x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
+    x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
+    x = x ^ (x >> 31);
+
+    // Convert hash to uniform random number in [0,1)
+    const double u =
+      static_cast<double>(x >> 11) *
+      (1.0 / 9007199254740992.0);
+
+    return u < cfgKeepFraction.value;
+  }
+
   template <typename T>
   bool selectionTrack(const T& candidate)
   {
@@ -386,6 +408,10 @@ struct cksspinalignment {
       return;
     }
 
+    if (!keepEvent(static_cast<uint64_t>(collision.globalIndex()))) {
+      return;
+    }
+
     histos.fill(HIST("hEvtSelInfo"), 2.5);
     histos.fill(HIST("hCent"), centrality);