[PWGCF] JCorran: Add nonflow SPC correlators and a smaller Q-vector fill grid

Author: Neelkamal Mallick

PWGCF/JCorran/Core/FlowJSPCAnalysis.cxx

@@ -125,6 +125,19 @@ void FlowJSPCAnalysis::calculateCorrelators(const int fCentBin)
     correlationDenom = 0.;
     weightCorrelationDenom = 0.;
   }
+
+  // N_m = Re(Q_{0,1}); weight M_m. 3SPC → fN3, 4SPC → fN4.
+  if (mHistRegistry && qvecs) {
+    const double nSel = qvecs->QvectorQC[0][1].Re();
+    const float centX = static_cast<float>(fCentBin) + 0.5f;
+    if (nSel > 0.0 && std::isfinite(nSel)) {
+      if (mWhichSPC == 0 && fCorrelDenoms[2] > 0.0) {
+        mHistRegistry->fill(HIST("fN3"), centX, nSel, fCorrelDenoms[2]);
+      } else if (mWhichSPC == 1 && fCorrelDenoms[3] > 0.0) {
+        mHistRegistry->fill(HIST("fN4"), centX, nSel, fCorrelDenoms[3]);
+      }
+    }
+  }
 }
 
 void FlowJSPCAnalysis::fillHistograms(const int fCentBin, int ind, double cNum, double cDenom, double wNum, double wDenom)

PWGCF/JCorran/Core/FlowJSPCAnalysis.h

@@ -29,6 +29,7 @@
 
 #include <Rtypes.h>
 
+#include <cstdint>
 #include <cstring>
 #include <string_view>
 #include <vector>
@@ -42,6 +43,27 @@ class FlowJSPCAnalysis
   int getCentBin(float cValue);
 
   using JQVectorsT = JQVectors<TComplex, 113, 15, false>;
+  /// Fill-grid size for Q_{n,p}: (v8 * nPartDen)+1 harmonics, nPartDen+1 powers.
+  static constexpr uint32_t Nh3p = 49;  // 3-particle SPC, 6-particle denominator
+  static constexpr uint32_t Nk3p = 7;
+  static constexpr uint32_t Nh4p = 65;  // 4-particle SPC, 8-particle denominator
+  static constexpr uint32_t Nk4p = 9;
+  static constexpr uint32_t NhFull = 113;
+  static constexpr uint32_t NkFull = 15;
+  void qVectorGrid(int whichSPC, uint32_t& nhUse, uint32_t& nkUse)
+  {
+    mWhichSPC = whichSPC;
+    if (whichSPC == 0) {
+      nhUse = Nh3p;
+      nkUse = Nk3p;
+    } else if (whichSPC == 1) {
+      nhUse = Nh4p;
+      nkUse = Nk4p;
+    } else {
+      nhUse = NhFull;
+      nkUse = NkFull;
+    }
+  }
   inline void setQvectors(const JQVectorsT* _qvecs) { qvecs = _qvecs; }
   void correlation(int c_nPart, int c_nHarmo, int* harmo, double* correlData);
   void calculateCorrelators(const int fCentBin);
@@ -57,6 +79,9 @@ class FlowJSPCAnalysis
       return;
     }
     mHistRegistry->add("FullCentrality", "FullCentrality", o2::framework::HistType::kTH1D, {{100, 0., 100.}}, true);
+    // Effective N_m per centrality class, weighted by M_m (arXiv:2606.10258 c0).
+    mHistRegistry->add("fN3", "Effective N_{3};centrality class;N_{3}", {o2::framework::HistType::kTProfile, {{9, 0., 9.}}}, true);
+    mHistRegistry->add("fN4", "Effective N_{4};centrality class;N_{4}", {o2::framework::HistType::kTProfile, {{9, 0., 9.}}}, true);
     mHistRegistry->add("Centrality_0/fResults", "Numerators and denominators", {o2::framework::HistType::kTProfile, {{24, 0., 24.}}}, true);
     mHistRegistry->add("Centrality_0/fCovResults", "Covariance N*D", {o2::framework::HistType::kTProfile, {{48, 0., 48.}}}, true);
     mHistRegistry->add("Centrality_0/phiBefore", "Phi before", {o2::framework::HistType::kTH1D, {{100, 0., o2::constants::math::TwoPI}}}, true);
@@ -67,13 +92,13 @@ class FlowJSPCAnalysis
     }
   }
 
-  void setCorrSet(int obsInd, int harmo[8])
+  void setCorrSet(int obsInd, int const harmo[8])
   {
     for (int i = 0; i < 8; i++) {
       fHarmosArray[obsInd][i] = harmo[i];
     }
   }
-  void setFullCorrSet(int harmo[12][8])
+  void setFullCorrSet(int const harmo[12][8])
   {
     memcpy(fHarmosArray, harmo, sizeof(int) * 12 * 8);
   }
@@ -92,13 +117,14 @@ class FlowJSPCAnalysis
  private:
   const int mNqHarmos = 113; ///< Highest harmo for Q(n,p): (v8*14part)+1.
   const int mNqPowers = 15;  ///< Max power for Q(n,p): 14part+1.
-  const JQVectorsT* qvecs;
+  const JQVectorsT* qvecs = nullptr;
 
   o2::framework::HistogramRegistry* mHistRegistry = nullptr;
 
-  int fHarmosArray[12][8];
+  int fHarmosArray[12][8] = {{0}};
 
-  double fCorrelDenoms[14];
+  double fCorrelDenoms[14] = {0};
+  int mWhichSPC = 0;
 
   ClassDefNV(FlowJSPCAnalysis, 1);
 };

PWGCF/JCorran/Core/FlowJSPCObservables.h

@@ -38,42 +38,53 @@ class FlowJSPCObservables
     switch (index) {
       case 0: {
         LOGF(info, "Computing three harmonic SPC");
-        int harmonicArray01[maxNrComb][8] = {
-          {3, 6, -3, -3, 0, 0, 0, 0},
-          {3, 4, -2, -2, 0, 0, 0, 0},
+        // fResults slot j: num at 2j+0.5, denom at 2j+1.5. arXiv:2606.10258 nonflow refs in unused slots.
+        int const harmonicArray01[maxNrComb][8] = {
+          {3, 6, -3, -3, 0, 0, 0, 0}, // 0: C633 = <V3 V3 V6*>
+          {3, 4, -2, -2, 0, 0, 0, 0}, // 1: C422 = <V2 V2 V4*>
           {3, 8, -4, -4, 0, 0, 0, 0},
-          {3, 2, 4, -6, 0, 0, 0, 0},
-          {3, 2, 3, -5, 0, 0, 0, 0},
+          {3, 2, 4, -6, 0, 0, 0, 0}, // 3: C246 = <V2 V4 V6*>
+          {3, 2, 3, -5, 0, 0, 0, 0}, // 4: C235 = <V2 V3 V5*>
           {3, 3, 4, -7, 0, 0, 0, 0}, // These are three harmonic SPC!!
           {3, 2, 5, -7, 0, 0, 0, 0}, // These are three harmonic SPC!!
           {3, 3, 5, -8, 0, 0, 0, 0}, // These are three harmonic SPC!!
-          {0, 6, -2, -2, -2, 0, 0, 0},
-          {0, 2, -3, -4, 5, 0, 0, 0},
-          {0, 2, -3, -3, 4, 0, 0, 0},
-          {0, 3, 3, -2, -2, -2, 0, 0}};
+          // {0, 6, -2, -2, -2, 0, 0, 0},
+          // {0, 2, -3, -4, 5, 0, 0, 0},
+          // {0, 2, -3, -3, 4, 0, 0, 0},
+          // {0, 3, 3, -2, -2, -2, 0, 0},
+          {3, 1, 1, -2, 0, 0, 0, 0}, // 8: C112 = <V1 V1 V2*>, Eqs. (IV.7), (IV.18)
+          {3, 1, 2, -3, 0, 0, 0, 0}, // 9: C123 = <V1 V2 V3*>, Eqs. (IV.8), (IV.18)
+          {0, 0, 0, 0, 0, 0, 0, 0},
+          {0, 0, 0, 0, 0, 0, 0, 0}};
 
         memcpy(harmonicArray, harmonicArray01, sizeof(int) * maxNrComb * 8);
       } break;
       case 1: {
         LOGF(info, "Computing four harmonic SPC");
-        int harmonicArray02[maxNrComb][8] = {
-          {4, 6, -2, -2, -2, 0, 0, 0},
+        // fResults slot j: num at 2j+0.5, denom at 2j+1.5. arXiv:2606.10258: c1{4}=<<4>>-2<<2>>^2 after averaging, Eq. (IV.6).
+        int const harmonicArray02[maxNrComb][8] = {
+          {4, 6, -2, -2, -2, 0, 0, 0}, // 0: C6222 = <V2 V2 V2 V6*>
           {4, 2, -3, -4, 5, 0, 0, 0},
           {4, 2, -3, -3, 4, 0, 0, 0},
           {4, 2, 2, 3, -7, 0, 0, 0}, // These are three harmonic SPC!!
           {4, 2, 2, 4, -8, 0, 0, 0}, // These are three harmonic SPC!!
           {4, 2, 7, -4, -5, 0, 0, 0},
           {4, 3, -4, -4, 5, 0, 0, 0},
-          {0, 0, 0, 0, 0, 0, 0, 0},
-          {0, 0, 0, 0, 0, 0, 0, 0},
-          {0, 0, 0, 0, 0, 0, 0, 0},
+          // {0, 0, 0, 0, 0, 0, 0, 0},
+          // {0, 0, 0, 0, 0, 0, 0, 0},
+          // {0, 0, 0, 0, 0, 0, 0, 0},
+          // {0, 0, 0, 0, 0, 0, 0, 0},
+          // {0, 0, 0, 0, 0, 0, 0, 0},
+          {4, 1, 1, -1, -1, 0, 0, 0}, // 7: <<4>>_{1,1,-1,-1} for c1{4}
+          {2, 1, -1, 0, 0, 0, 0, 0}, // 8: <<2>>_{1,-1} = <V1 V1*>
+          {3, 1, 1, -2, 0, 0, 0, 0}, // 9: C112 on the 4-particle sample (mixed-order Eq. (IV.19))
           {0, 0, 0, 0, 0, 0, 0, 0},
           {0, 0, 0, 0, 0, 0, 0, 0}};
         memcpy(harmonicArray, harmonicArray02, sizeof(int) * maxNrComb * 8);
       } break;
       case 2: {
         LOGF(info, "Computing five and six harmonic SPC");
-        int harmonicArray03[maxNrComb][8] = {
+        int const harmonicArray03[maxNrComb][8] = {
           {5, 3, 3, -2, -2, -2, 0, 0},
           {5, 2, 2, -3, 4, -5, 0, 0},
           {5, 2, 3, 3, -4, -4, 0, 0},
@@ -90,7 +101,7 @@ class FlowJSPCObservables
       } break;
       case 3: {
         LOGF(info, "Computing slected five harmonic SPC");
-        int harmonicArray04[maxNrComb][8] = {
+        int const harmonicArray04[maxNrComb][8] = {
           {5, 3, 3, -2, -2, -2, 0, 0},
           {0, 2, 2, -3, 4, -5, 0, 0},
           {5, 2, 3, 3, -4, -4, 0, 0},

PWGCF/JCorran/Core/JQVectors.h

@@ -19,6 +19,8 @@
 
 #include <RtypesCore.h>
 
+#include <algorithm>
+#include <cstdint>
 #include <experimental/type_traits>
 #include <type_traits>
 
@@ -52,19 +54,24 @@ class JQVectors : public std::conditional_t<gap, JQVectorsGapBase<Q, nh, nk>, JQ
   using hasInvMass = decltype(std::declval<T&>().invMass());
 
   template <class JInputClass>
-  inline void Calculate(JInputClass& inputInst, float etamin, float etamax, float massMin = 0.0f, float massMax = 999.9f)
+  inline void Calculate(JInputClass& inputInst, float etamin, float etamax, float massMin = 0.0f, float massMax = 999.9f,
+                        uint32_t nhUse = nh, uint32_t nkUse = nk)
   {
+    // nhUse/nkUse limit the filled (harmonic, power) grid. Defaults keep the full template size.
+    const uint32_t nH = std::min(nhUse, nh);
+    const uint32_t nK = std::min(nkUse, nk);
+
     // calculate Q-vector for QC method ( no subgroup )
-    for (UInt_t ih = 0; ih < nh; ++ih) {
-      for (UInt_t ik = 0; ik < nk; ++ik) {
+    for (UInt_t ih = 0; ih < nH; ++ih) {
+      for (UInt_t ik = 0; ik < nK; ++ik) {
         QvectorQC[ih][ik] = Q(0, 0);
         if constexpr (gap) {
           for (UInt_t isub = 0; isub < 2; ++isub)
             this->QvectorQCgap[isub][ih][ik] = Q(0, 0);
         }
       }
     }
-    for (auto& track : inputInst) {
+    for (auto const& track : inputInst) {
       if (track.eta() < -etamax || track.eta() > etamax)
         continue;
       using JInputClassIter = typename JInputClass::iterator;
@@ -74,10 +81,15 @@ class JQVectors : public std::conditional_t<gap, JQVectorsGapBase<Q, nh, nk>, JQ
       }
 
       UInt_t isub = (UInt_t)(track.eta() > 0.0);
-      for (UInt_t ih = 0; ih < nh; ++ih) {
+      const Double_t phi = track.phi();
+      const Double_t c1 = TMath::Cos(phi);
+      const Double_t s1 = TMath::Sin(phi);
+      Double_t cn = 1.0; // cos(ih * phi), ih = 0
+      Double_t sn = 0.0; // sin(ih * phi)
+      for (UInt_t ih = 0; ih < nH; ++ih) {
         Double_t tf = 1.0;
-        for (UInt_t ik = 0; ik < nk; ++ik) {
-          Q q(tf * TMath::Cos(ih * track.phi()), tf * TMath::Sin(ih * track.phi()));
+        for (UInt_t ik = 0; ik < nK; ++ik) {
+          Q q(tf * cn, tf * sn);
           QvectorQC[ih][ik] += q;
 
           if constexpr (gap) {
@@ -90,6 +102,10 @@ class JQVectors : public std::conditional_t<gap, JQVectorsGapBase<Q, nh, nk>, JQ
           if constexpr (std::experimental::is_detected<hasWeightEff, const JInputClassIter>::value)
             tf *= track.weightEff();
         }
+        const Double_t cnNext = cn * c1 - sn * s1;
+        const Double_t snNext = cn * s1 + sn * c1;
+        cn = cnNext;
+        sn = snNext;
       }
     }
   }

PWGCF/JCorran/Tasks/flowJSPCAnalysis.cxx

@@ -109,6 +109,8 @@ struct flowJSPCAnalysis {
 
   std::unique_ptr<TFormula> multCutFormula;
   std::array<uint, aod::cfmultset::NMultiplicityEstimators> multCutFormulaParamIndex;
+  uint32_t mNhUse = FlowJSPCAnalysis::NhFull;
+  uint32_t mNkUse = FlowJSPCAnalysis::NkFull;
 
   void init(InitContext const&)
   {
@@ -118,6 +120,8 @@ struct flowJSPCAnalysis {
 
     spcObservables.setSPCObservables(cfgWhichSPC);
     spcAnalysis.setFullCorrSet(spcObservables.harmonicArray);
+    spcAnalysis.qVectorGrid(cfgWhichSPC, mNhUse, mNkUse);
+    LOGF(info, "Q-vector fill grid: nh=%u nk=%u (cfgWhichSPC=%d)", mNhUse, mNkUse, cfgWhichSPC.value);
 
     histManager.setHistRegistryQA(&qaHistRegistry);
     histManager.setDebugLog(false);
@@ -186,7 +190,7 @@ struct flowJSPCAnalysis {
     if (cfgFillQA)
       histManager.fillEventQA<1>(collision, cBin, cent, nTracks);
 
-    jqvecs.Calculate(tracks, 0.0, cfgTrackCuts.cfgEtaMax);
+    jqvecs.Calculate(tracks, 0.0, cfgTrackCuts.cfgEtaMax, 0.0f, 999.9f, mNhUse, mNkUse);
     spcAnalysis.setQvectors(&jqvecs);
     spcAnalysis.calculateCorrelators(cBin);
   }