[PWGCF] JCorran updates (#9754)

Author: jaelpark

PWGCF/JCorran/Core/JFFlucAnalysis.cxx

@@ -116,6 +116,13 @@ TComplex JFFlucAnalysis::Q(int n, int p)
   return n >= 0 ? pqvecs->QvectorQC[n][p] : C(pqvecs->QvectorQC[-n][p]);
 }
 
+TComplex JFFlucAnalysis::Q(const JQVectorsT& qvecs, int n, int p)
+{
+  // Return QvectorQC
+  // Q{-n, p} = Q{n, p}*
+  return n >= 0 ? qvecs.QvectorQC[n][p] : C(qvecs.QvectorQC[-n][p]);
+}
+
 TComplex JFFlucAnalysis::Two(int n1, int n2)
 {
   // two-particle correlation <exp[i(n1*phi1 + n2*phi2)]>
@@ -124,9 +131,26 @@ TComplex JFFlucAnalysis::Two(int n1, int n2)
 
 TComplex JFFlucAnalysis::Four(int n1, int n2, int n3, int n4)
 {
-
   return Q(n1, 1) * Q(n2, 1) * Q(n3, 1) * Q(n4, 1) - Q(n1 + n2, 2) * Q(n3, 1) * Q(n4, 1) - Q(n2, 1) * Q(n1 + n3, 2) * Q(n4, 1) - Q(n1, 1) * Q(n2 + n3, 2) * Q(n4, 1) + 2. * Q(n1 + n2 + n3, 3) * Q(n4, 1) - Q(n2, 1) * Q(n3, 1) * Q(n1 + n4, 2) + Q(n2 + n3, 2) * Q(n1 + n4, 2) - Q(n1, 1) * Q(n3, 1) * Q(n2 + n4, 2) + Q(n1 + n3, 2) * Q(n2 + n4, 2) + 2. * Q(n3, 1) * Q(n1 + n2 + n4, 3) - Q(n1, 1) * Q(n2, 1) * Q(n3 + n4, 2) + Q(n1 + n2, 2) * Q(n3 + n4, 2) + 2. * Q(n2, 1) * Q(n1 + n3 + n4, 3) + 2. * Q(n1, 1) * Q(n2 + n3 + n4, 3) - 6. * Q(n1 + n2 + n3 + n4, 4);
 }
+
+TComplex JFFlucAnalysis::TwoDiff(int n1, int n2)
+{
+#define dp(n, p) Q(*pqvecs, n, p)    // POI
+#define dQ(n, p) Q(*pqvecsRef, n, p) // REF
+#define dq(n, p) dp(n, p)            //(dp(n,p)+dQ(n,p)) //POI+REF in narrow bin. Since there is no mass for ref, q = POI
+                                     // #define dq(n,p) (dp(n,p)+dQ(n,p)) //POI+REF in narrow bin. Since there is no mass for ref, q = POI
+  return dp(n1, 1) * dQ(n2, 1) - dq(n1 + n2, 2);
+}
+
+TComplex JFFlucAnalysis::FourDiff(int n1, int n2, int n3, int n4)
+{
+  return dp(n1, 1) * dQ(n2, 1) * dQ(n3, 1) * dQ(n4, 1) - dq(n1 + n2, 2) * dQ(n3, 1) * dQ(n4, 1) - dq(n1 + n3, 2) * dQ(n2, 1) * dQ(n4, 1) - dp(n1, 1) * dQ(n2 + n3, 2) * dQ(n4, 1) + 2. * dq(n1 + n2 + n3, 3) * dQ(n4, 1) - dQ(n2, 1) * dQ(n3, 1) * dq(n1 + n4, 2) + dQ(n2 + n3, 2) * dq(n1 + n4, 2) - dp(n1, 1) * dQ(n3, 1) * dQ(n2 + n4, 2) + dq(n1 + n3, 2) * dQ(n2 + n4, 2) + 2. * dQ(n3, 1) * dq(n1 + n2 + n4, 3) - dp(n1, 1) * dQ(n2, 1) * dQ(n3 + n4, 2) + dq(n1 + n2, 2) * dQ(n3 + n4, 2) + 2. * dQ(n2, 1) * dq(n1 + n3 + n4, 3) + 2. * dp(n1, 1) * dQ(n2 + n3 + n4, 3) - 6. * dq(n1 + n2 + n3 + n4, 4);
+}
+
+#undef dp
+#undef dQ
+#undef dq
 #undef C
 
 //________________________________________________________________________
@@ -296,19 +320,21 @@ void JFFlucAnalysis::UserExec(Option_t* /*popt*/) // NOLINT(readability/casting)
     }
   }
 
+  auto four = [&](int a, int b, int c, int d) -> TComplex { return pqvecsRef ? FourDiff(a, b, c, d) : Four(a, b, c, d); };
+  auto two = [&](int a, int b) -> TComplex { return pqvecsRef ? TwoDiff(a, b) : Two(a, b); };
   Double_t event_weight_four = 1.0;
   Double_t event_weight_two = 1.0;
   if (flags & kFlucEbEWeighting) {
-    event_weight_four = Four(0, 0, 0, 0).Re();
-    event_weight_two = Two(0, 0).Re();
+    event_weight_four = four(0, 0, 0, 0).Re();
+    event_weight_two = two(0, 0).Re();
   }
 
   for (UInt_t ih = 2; ih < kNH; ih++) {
     for (UInt_t ihh = 2, mm = (ih < kcNH ? ih : static_cast<UInt_t>(kcNH)); ihh < mm; ihh++) {
-      TComplex scfour = Four(ih, ihh, -ih, -ihh) / Four(0, 0, 0, 0).Re();
+      TComplex scfour = four(ih, ihh, -ih, -ihh) / four(0, 0, 0, 0).Re();
       pht[HIST_THN_SC_with_QC_4corr]->Fill(fCent, fAvgInvariantMass, ih, ihh, scfour.Re(), event_weight_four);
     }
-    TComplex sctwo = Two(ih, -ih) / Two(0, 0).Re();
+    TComplex sctwo = two(ih, -ih) / two(0, 0).Re();
     pht[HIST_THN_SC_with_QC_2corr]->Fill(fCent, fAvgInvariantMass, ih, sctwo.Re(), event_weight_two);
   }
 }

PWGCF/JCorran/Core/JFFlucAnalysis.h

@@ -37,6 +37,8 @@ class JFFlucAnalysis : public TNamed
   TComplex Q(int n, int p);
   TComplex Two(int n1, int n2);
   TComplex Four(int n1, int n2, int n3, int n4);
+  TComplex TwoDiff(int n1, int n2);
+  TComplex FourDiff(int n1, int n2, int n3, int n4);
   void UserExec(Option_t* option);
   void Terminate(Option_t*);
 
@@ -124,6 +126,7 @@ class JFFlucAnalysis : public TNamed
          kK4,
          nKL }; // order
   using JQVectorsT = JQVectors<TComplex, kNH, nKL, true>;
+  TComplex Q(const JQVectorsT& qvecs, int n, int p);
   inline void SetJQVectors(const JQVectorsT* _pqvecs)
   {
     pqvecs = _pqvecs;
@@ -140,10 +143,10 @@ class JFFlucAnalysis : public TNamed
   template <class T>
   using hasWeightEff = decltype(std::declval<T&>().weightEff());
   template <class T>
-  using hasType = decltype(std::declval<T&>().particleType());
+  using hasSign = decltype(std::declval<T&>().sign());
 
   template <class JInputClass>
-  inline void FillQA(JInputClass& inputInst, UInt_t type = 0)
+  inline void FillQA(JInputClass& inputInst, UInt_t type = 0u)
   {
     ph1[HIST_TH1_CENTRALITY]->Fill(fCent);
     ph1[HIST_TH1_IMPACTPARAM]->Fill(fImpactParameter);
@@ -153,19 +156,20 @@ class JFFlucAnalysis : public TNamed
       using JInputClassIter = typename JInputClass::iterator;
       if constexpr (std::experimental::is_detected<hasWeightEff, const JInputClassIter>::value)
         corrInv /= track.weightEff();
-      pht[HIST_THN_PTETA]->Fill(fCent, track.pt(), track.eta(), corrInv);
+      if constexpr (std::experimental::is_detected<hasSign, const JInputClassIter>::value)
+        pht[HIST_THN_PTETA]->Fill(fCent, track.pt(), track.eta(), track.sign(), corrInv);
+      else
+        pht[HIST_THN_PTETA]->Fill(fCent, track.pt(), track.eta(), 0.0, corrInv);
       if constexpr (std::experimental::is_detected<hasWeightNUA, const JInputClassIter>::value)
         corrInv /= track.weightNUA();
       pht[HIST_THN_PHIETA]->Fill(fCent, track.phi(), track.eta(), corrInv);
-      if constexpr (std::experimental::is_detected<hasType, const JInputClassIter>::value)
-        type = track.particleType();
       pht[HIST_THN_PHIETAZ]->Fill(fCent, static_cast<Double_t>(type), track.phi(), track.eta(), fVertex, corrInv);
     }
 
     ph1[HIST_TH1_ZVERTEX]->Fill(fVertex);
   }
 
-#define kcNH kH6 // max second dimension + 1
+#define kcNH kH4 // max second dimension + 1
  protected:
   Float_t fVertex;           //!
   Float_t fAvgInvariantMass; //!
@@ -178,7 +182,7 @@ class JFFlucAnalysis : public TNamed
   const JQVectorsT* pqvecsRef; //!
 
   TH1* ph1[HIST_TH1_COUNT];              //!
-  THn* pht[HIST_THN_COUNT];              //!
+  THnSparse* pht[HIST_THN_COUNT];        //!
   THnSparse* phs[HIST_THN_SPARSE_COUNT]; //!
 
   ClassDef(JFFlucAnalysis, 1)

PWGCF/JCorran/Core/JFFlucAnalysisO2Hist.cxx

@@ -18,28 +18,26 @@
 
 using namespace o2;
 
-JFFlucAnalysisO2Hist::JFFlucAnalysisO2Hist(HistogramRegistry& registry, AxisSpec& axisMultiplicity, AxisSpec& phiAxis, AxisSpec& etaAxis, AxisSpec& zvtAxis, AxisSpec& massAxis, const TString& folder) : JFFlucAnalysis()
+JFFlucAnalysisO2Hist::JFFlucAnalysisO2Hist(HistogramRegistry& registry, AxisSpec& axisMultiplicity, AxisSpec& phiAxis, AxisSpec& etaAxis, AxisSpec& zvtAxis, AxisSpec& ptAxis, AxisSpec& massAxis, const TString& folder) : JFFlucAnalysis()
 {
 
   ph1[HIST_TH1_CENTRALITY] = std::get<std::shared_ptr<TH1>>(registry.add(Form("%s/h_cent", folder.Data()), "multiplicity/centrality", {HistType::kTH1F, {axisMultiplicity}})).get();
   ph1[HIST_TH1_IMPACTPARAM] = std::get<std::shared_ptr<TH1>>(registry.add(Form("%s/h_IP", folder.Data()), "impact parameter", {HistType::kTH1F, {{400, -2.0, 20.0}}})).get();
   ph1[HIST_TH1_ZVERTEX] = std::get<std::shared_ptr<TH1>>(registry.add(Form("%s/h_vertex", folder.Data()), "z vertex", {HistType::kTH1F, {{100, -20.0, 20.0}}})).get();
-  //
-  // TODO: these shall be configurable
-  std::vector<double> ptBinning = {0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0};
-  AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/c)"};
+
+  AxisSpec chgAxis = {3, -1.5, 1.5, "charge"};
   AxisSpec typeAxis = {2, -0.5, 1.5, "type"};
-  pht[HIST_THN_PHIETAZ] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_phietaz", folder.Data()), "multiplicity/centrality, type, phi, eta, z", {HistType::kTHnF, {axisMultiplicity, typeAxis, phiAxis, etaAxis, zvtAxis}})).get();
-  pht[HIST_THN_PTETA] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_pteta", folder.Data()), "(corrected) multiplicity/centrality, pT, eta", {HistType::kTHnF, {axisMultiplicity, ptAxis, etaAxis}})).get();
-  pht[HIST_THN_PHIETA] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_phieta", folder.Data()), "(corrected) multiplicity/centrality, phi, eta", {HistType::kTHnF, {axisMultiplicity, phiAxis, etaAxis}})).get();
+  pht[HIST_THN_PHIETAZ] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_phietaz", folder.Data()), "multiplicity/centrality, type, phi, eta, z", {HistType::kTHnSparseF, {axisMultiplicity, typeAxis, phiAxis, etaAxis, zvtAxis}})).get();
+  pht[HIST_THN_PTETA] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_pteta", folder.Data()), "(corrected) multiplicity/centrality, pT, eta, charge", {HistType::kTHnSparseF, {axisMultiplicity, ptAxis, etaAxis, chgAxis}})).get();
+  pht[HIST_THN_PHIETA] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_phieta", folder.Data()), "(corrected) multiplicity/centrality, phi, eta", {HistType::kTHnSparseF, {axisMultiplicity, phiAxis, etaAxis}})).get();
   AxisSpec hAxis = {kNH, -0.5, static_cast<double>(kNH - 1) + 0.5, "#it{n}"};
   AxisSpec kAxis = {nKL, -0.5, static_cast<double>(nKL - 1) + 0.5, "#it{k}"};
   AxisSpec vnAxis = {2048, -0.1, 0.1, "#it{V}_#it{n}"};
-  pht[HIST_THN_SC_with_QC_4corr] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_SC_with_QC_4corr", folder.Data()), "SC_with_QC_4corr", {HistType::kTHnF, {axisMultiplicity, massAxis, hAxis, hAxis, {2048, -0.001, 0.001, "correlation"}}})).get();
-  pht[HIST_THN_SC_with_QC_2corr] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_SC_with_QC_2corr", folder.Data()), "SC_with_QC_2corr", {HistType::kTHnF, {axisMultiplicity, massAxis, hAxis, {2048, -0.1, 0.1, "correlation"}}})).get();
-  pht[HIST_THN_SC_with_QC_2corr_gap] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_SC_with_QC_2corr_gap", folder.Data()), "SC_with_QC_2corr_gap", {HistType::kTHnF, {axisMultiplicity, massAxis, hAxis, {2048, -0.1, 0.1, "correlation"}}})).get();
+  pht[HIST_THN_SC_with_QC_4corr] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_SC_with_QC_4corr", folder.Data()), "SC_with_QC_4corr", {HistType::kTHnSparseF, {axisMultiplicity, massAxis, hAxis, hAxis, {2048, -0.001, 0.001, "correlation"}}})).get();
+  pht[HIST_THN_SC_with_QC_2corr] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_SC_with_QC_2corr", folder.Data()), "SC_with_QC_2corr", {HistType::kTHnSparseF, {axisMultiplicity, massAxis, hAxis, {2048, -0.1, 0.1, "correlation"}}})).get();
+  pht[HIST_THN_SC_with_QC_2corr_gap] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_SC_with_QC_2corr_gap", folder.Data()), "SC_with_QC_2corr_gap", {HistType::kTHnSparseF, {axisMultiplicity, massAxis, hAxis, {2048, -0.1, 0.1, "correlation"}}})).get();
   for (UInt_t i = HIST_THN_V4V2star_2; i < HIST_THN_COUNT; ++i)
-    pht[i] = std::get<std::shared_ptr<THn>>(registry.add(Form("%s/h_corrC%02u", folder.Data(), i - HIST_THN_V4V2star_2), "correlator", {HistType::kTHnF, {axisMultiplicity, massAxis, {2048, -3.0, 3.0, "correlation"}}})).get();
+    pht[i] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_corrC%02u", folder.Data(), i - HIST_THN_V4V2star_2), "correlator", {HistType::kTHnSparseF, {axisMultiplicity, massAxis, {2048, -3.0, 3.0, "correlation"}}})).get();
   for (UInt_t i = 0; i < HIST_THN_COUNT; ++i)
     pht[i]->Sumw2();
 

PWGCF/JCorran/Core/JFFlucAnalysisO2Hist.h

@@ -23,7 +23,7 @@ using namespace o2::framework;
 class JFFlucAnalysisO2Hist : public JFFlucAnalysis
 {
  public:
-  JFFlucAnalysisO2Hist(HistogramRegistry&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, const TString&);
+  JFFlucAnalysisO2Hist(HistogramRegistry&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, const TString&);
   ~JFFlucAnalysisO2Hist();
 };
 

PWGCF/JCorran/Core/JQVectors.h

@@ -65,7 +65,7 @@ class JQVectors : public std::conditional_t<gap, JQVectorsGapBase<Q, nh, nk>, JQ
         continue;
       using JInputClassIter = typename JInputClass::iterator;
       if constexpr (std::experimental::is_detected<hasInvMass, const JInputClassIter>::value) {
-        if (track.invMass() < massMin || track.invMass() > massMax)
+        if (track.invMass() < massMin || track.invMass() >= massMax)
           continue;
       }
 

PWGCF/JCorran/Tasks/jflucAnalysisTask.cxx

@@ -57,6 +57,7 @@ struct jflucAnalysisTask {
   ConfigurableAxis phiAxis{"axisPhi", {50, 0.0, o2::constants::math::TwoPI}, "phi axis for histograms"};
   ConfigurableAxis etaAxis{"axisEta", {40, -2.0, 2.0}, "eta axis for histograms"};
   ConfigurableAxis zvtAxis{"axisZvt", {20, -10.0, 10.0}, "zvertex axis for histograms"};
+  ConfigurableAxis ptAxis{"axisPt", {60, 0.0, 300.0}, "pt axis for histograms"};
   ConfigurableAxis massAxis{"axisMass", {1, 0.0, 10.0}, "mass axis for histograms"};
 
   Filter jtrackFilter = (aod::jtrack::pt > ptmin) && (aod::jtrack::pt < ptmax);    // eta cuts done by jfluc
@@ -71,16 +72,17 @@ struct jflucAnalysisTask {
     auto axisSpecPhi = AxisSpec(phiAxis);
     auto axisSpecEta = AxisSpec(etaAxis);
     auto axisSpecZvt = AxisSpec(zvtAxis);
+    auto axisSpecPt = AxisSpec(ptAxis);
     auto axisSpecMass = AxisSpec(massAxis);
     if (doprocessJDerived || doprocessJDerivedCorrected || doprocessCFDerived || doprocessCFDerivedCorrected) {
-      pcf = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecMass, "jfluc");
+      pcf = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecPt, axisSpecMass, "jfluc");
       pcf->AddFlags(JFFlucAnalysis::kFlucEbEWeighting);
       pcf->UserCreateOutputObjects();
     } else {
       pcf = 0;
     }
     if (doprocessCF2ProngDerived || doprocessCF2ProngDerivedCorrected) {
-      pcf2Prong = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecMass, "jfluc2prong");
+      pcf2Prong = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecPt, axisSpecMass, "jfluc2prong");
       pcf2Prong->AddFlags(JFFlucAnalysis::kFlucEbEWeighting);
       pcf2Prong->UserCreateOutputObjects();
 
@@ -116,8 +118,8 @@ struct jflucAnalysisTask {
     pcf2Prong->Init();
     pcf2Prong->SetEventCentrality(collision.multiplicity());
     pcf2Prong->SetEventVertex(collision.posZ());
-    pcf2Prong->FillQA(poiTracks, 1u); // type = 1, all POI tracks in this list are of the same type
     pcf2Prong->FillQA(refTracks, 0u);
+    pcf2Prong->FillQA(poiTracks, 1u); // type = 1, all POI tracks in this list are of the same type
     qvecsRef.Calculate(refTracks, etamin, etamax);
     pcf2Prong->SetJQVectors(&qvecs, &qvecsRef);
     const AxisSpec& a = AxisSpec(massAxis);

PWGCF/JCorran/Tasks/jflucWeightsLoader.cxx

@@ -12,6 +12,7 @@
 /// \since May 2024
 
 #include <experimental/type_traits>
+#include <string>
 #include <TFile.h>
 #include <THn.h>
 
@@ -94,7 +95,8 @@ struct jflucWeightsLoader {
       float phiWeight, effWeight;
       if (ph) {
         UInt_t partType = 0; // partType 0 = all charged hadrons
-        if constexpr (std::experimental::is_detected<hasDecay, typename TrackT::iterator>::value) {
+        // TODO: code below to be enabled
+        /*if constexpr (std::experimental::is_detected<hasDecay, typename TrackT::iterator>::value) {
           switch (track.decay()) {
             case aod::cf2prongtrack::D0ToPiK:
             case aod::cf2prongtrack::D0barToKPi:
@@ -103,7 +105,7 @@ struct jflucWeightsLoader {
             default:
               break;
           }
-        }
+        }*/
         const Double_t coords[] = {collision.multiplicity(), static_cast<Double_t>(partType), track.phi(), track.eta(), collision.posZ()};
         phiWeight = ph->GetBinContent(ph->GetBin(coords));
       } else {

PWGCF/TableProducer/filterCorrelations.cxx

@@ -22,7 +22,6 @@
 #include "Common/DataModel/Centrality.h"
 
 #include <TH3F.h>
-#include <TDatabasePDG.h>
 
 using namespace o2;
 using namespace o2::framework;