[PWGCF] CF correlations and JCorran multiplicity correlations

PWGCF/DataModel/CorrelationsDerived.h

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -8,7 +8,7 @@
 // In applying this license CERN does not waive the privileges and immunities
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 #ifndef PWGCF_DATAMODEL_CORRELATIONSDERIVED_H_
 #define PWGCF_DATAMODEL_CORRELATIONSDERIVED_H_
 
 #include "Common/DataModel/Centrality.h"
@@ -31,7 +31,7 @@
 
 namespace cfmcparticle
 {
 DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision
 DECLARE_SOA_COLUMN(Pt, pt, float);                      //! pT (GeV/c)
 DECLARE_SOA_COLUMN(Eta, eta, float);                    //! Pseudorapidity
 DECLARE_SOA_COLUMN(Phi, phi, float);                    //! Phi angle
@@ -48,7 +48,7 @@
 
 namespace cfcollision
 {
 DECLARE_SOA_INDEX_COLUMN(CFMcCollision, cfMcCollision); //! Index to reduced MC collision
 DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float);  //! Centrality/multiplicity value
 } // namespace cfcollision
 DECLARE_SOA_TABLE(CFCollisions, "AOD", "CFCOLLISION", //! Reduced collision table
@@ -64,8 +64,8 @@
 
 namespace cftrack
 {
 DECLARE_SOA_INDEX_COLUMN(CFCollision, cfCollision);   //! Index to collision
 DECLARE_SOA_INDEX_COLUMN(CFMcParticle, cfMCParticle); //! Index to MC particle
 DECLARE_SOA_COLUMN(Pt, pt, float);                    //! pT (GeV/c)
 DECLARE_SOA_COLUMN(Eta, eta, float);                  //! Pseudorapidity
 DECLARE_SOA_COLUMN(Phi, phi, float);                  //! Phi angle
@@ -90,7 +90,14 @@
 namespace cfmultset
 {
 DECLARE_SOA_COLUMN(Multiplicities, multiplicities, std::vector<float>); //! List of auxiliary multiplicities
-}
+enum MultiplicityEstimators : uint8_t {
+  CentFT0C = 0x1,
+  MultFV0A = 0x2,
+  MultNTracksPV = 0x4,
+  MultNTracksGlobal = 0x8
+};
+
+} // namespace cfmultset
 DECLARE_SOA_TABLE(CFMultSets, "AOD", "CFMULTSET", cfmultset::Multiplicities); //! Auxilary multiplicity set table
 
 using CFMultSet = CFMultSets::iterator;
@@ -120,8 +127,8 @@
 
 namespace cf2prongtrack
 {
 DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng0, cfTrackProng0, int, CFTracks, "_0"); //! Index to prong 1 CFTrack
 DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng1, cfTrackProng1, int, CFTracks, "_1"); //! Index to prong 2 CFTrack
 DECLARE_SOA_COLUMN(Pt, pt, float);                                                //! pT (GeV/c)
 DECLARE_SOA_COLUMN(Eta, eta, float);                                              //! Pseudorapidity
 DECLARE_SOA_COLUMN(Phi, phi, float);                                              //! Phi angle

PWGCF/JCorran/Core/JFFlucAnalysis.cxx

@@ -12,10 +12,12 @@
 /// \author Jasper Parkkila (jparkkil@cern.ch)
 /// \since Sep 2022
 
-#include <TMath.h>
+#include "JFFlucAnalysis.h"
+
 #include <TComplex.h>
+#include <TMath.h>
+
 #include <algorithm>
-#include "JFFlucAnalysis.h"
 
 JFFlucAnalysis::JFFlucAnalysis() : TNamed(),
                                    fVertex(0),

PWGCF/JCorran/Core/JFFlucAnalysis.h

@@ -15,13 +15,16 @@
 #ifndef PWGCF_JCORRAN_CORE_JFFLUCANALYSIS_H_
 #define PWGCF_JCORRAN_CORE_JFFLUCANALYSIS_H_
 
-#include <experimental/type_traits>
 #include "JQVectors.h"
+
 #include <TComplex.h>
-#include <TNamed.h>
 #include <TH1.h>
 #include <THn.h>
 #include <THnSparse.h>
+#include <TNamed.h>
+
+#include <experimental/type_traits>
+#include <vector>
 
 class JFFlucAnalysis : public TNamed
 {
@@ -61,7 +64,6 @@ class JFFlucAnalysis : public TNamed
   enum HIST_THN {
     HIST_THN_PHIETAZ,
     HIST_THN_PTETA,
-    HIST_THN_PHIETA,
     HIST_THN_SC_with_QC_4corr,
     HIST_THN_SC_with_QC_2corr,
     HIST_THN_SC_with_QC_2corr_gap,
@@ -98,6 +100,7 @@ class JFFlucAnalysis : public TNamed
   enum HIST_THN_SPARSE {
     HIST_THN_SPARSE_VN,
     HIST_THN_SPARSE_VN_VN,
+    HIST_THN_SPARSE_MULTCORR,
     HIST_THN_SPARSE_COUNT
   };
   enum {
@@ -144,12 +147,17 @@ class JFFlucAnalysis : public TNamed
   using hasWeightEff = decltype(std::declval<T&>().weightEff());
   template <class T>
   using hasSign = decltype(std::declval<T&>().sign());
+  template <class T>
+  using hasMultSet = decltype(std::declval<T&>().multiplicities());
 
   template <class JInputClass>
   inline void FillQA(JInputClass& inputInst, UInt_t type = 0u)
   {
-    ph1[HIST_TH1_CENTRALITY]->Fill(fCent);
-    ph1[HIST_TH1_IMPACTPARAM]->Fill(fImpactParameter);
+    if (type == 0u) {
+      ph1[HIST_TH1_CENTRALITY]->Fill(fCent);
+      ph1[HIST_TH1_ZVERTEX]->Fill(fVertex);
+      ph1[HIST_TH1_IMPACTPARAM]->Fill(fImpactParameter);
+    }
 
     for (auto& track : inputInst) {
       Double_t corrInv = 1.0;
@@ -162,11 +170,19 @@ class JFFlucAnalysis : public TNamed
         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);
       pht[HIST_THN_PHIETAZ]->Fill(fCent, static_cast<Double_t>(type), track.phi(), track.eta(), fVertex, corrInv);
     }
+  }
 
-    ph1[HIST_TH1_ZVERTEX]->Fill(fVertex);
+  template <class JEventClass>
+  inline void FillMultSet(JEventClass& event)
+  {
+    //
+    if constexpr (std::experimental::is_detected<hasMultSet, const JEventClass>::value) {
+      // need to convert to vec of doubles since THnSparse has no way to fill vec of floats directly
+      std::vector<double> v(event.multiplicities().begin(), event.multiplicities().end());
+      phs[HIST_THN_SPARSE_MULTCORR]->Fill(v.data());
+    }
   }
 
 #define kcNH kH4 // max second dimension + 1

PWGCF/JCorran/Core/JFFlucAnalysisO2Hist.cxx

@@ -12,24 +12,41 @@
 /// \since Sep 2024
 
 #include "JFFlucAnalysisO2Hist.h"
-#include <vector>
-#include <memory>
+
+#include "PWGCF/DataModel/CorrelationsDerived.h"
+
 #include "CommonConstants/MathConstants.h"
 
+#include <memory>
+#include <vector>
+
 using namespace o2;
 
-JFFlucAnalysisO2Hist::JFFlucAnalysisO2Hist(HistogramRegistry& registry, AxisSpec& axisMultiplicity, AxisSpec& phiAxis, AxisSpec& etaAxis, AxisSpec& zvtAxis, AxisSpec& ptAxis, AxisSpec& massAxis, const TString& folder) : JFFlucAnalysis()
+JFFlucAnalysisO2Hist::JFFlucAnalysisO2Hist(HistogramRegistry& registry, AxisSpec& axisMultiplicity, AxisSpec& phiAxis, AxisSpec& etaAxis, AxisSpec& zvtAxis, AxisSpec& ptAxis, AxisSpec& massAxis, uint16_t multCorrMask, 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();
 
+  if (multCorrMask != 0) {
+    std::vector<AxisSpec> multAxes;
+    if (multCorrMask & aod::cfmultset::CentFT0C)
+      multAxes.emplace_back(100, 0, 100, "FT0C centrality");
+    if (multCorrMask & aod::cfmultset::MultFV0A)
+      multAxes.emplace_back(100, 0, 100000, "V0A multiplicity");
+    if (multCorrMask & aod::cfmultset::MultNTracksPV)
+      multAxes.emplace_back(100, 0, 1000, "Nch PV");
+    if (multCorrMask & aod::cfmultset::MultNTracksGlobal)
+      multAxes.emplace_back(100, 0, 1000, "Nch Global");
+    registry.add("multCorrelations", "Multiplicity correlations", {HistType::kTHnSparseF, multAxes});
+    phs[HIST_THN_SPARSE_MULTCORR] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_multcorr", folder.Data()), "multiplicity/centrality correlations", {HistType::kTHnSparseF, multAxes})).get();
+  }
+
   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<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_PHIETAZ] = std::get<std::shared_ptr<THnSparse>>(registry.add(Form("%s/h_phietaz", folder.Data()), "(corrected) 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}"};

PWGCF/JCorran/Core/JFFlucAnalysisO2Hist.h

@@ -15,6 +15,7 @@
 #define PWGCF_JCORRAN_CORE_JFFLUCANALYSISO2HIST_H_
 
 #include "JFFlucAnalysis.h"
+
 #include "Framework/HistogramRegistry.h"
 
 using namespace o2;
@@ -23,7 +24,7 @@ using namespace o2::framework;
 class JFFlucAnalysisO2Hist : public JFFlucAnalysis
 {
  public:
-  JFFlucAnalysisO2Hist(HistogramRegistry&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, const TString&);
+  JFFlucAnalysisO2Hist(HistogramRegistry&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, AxisSpec&, uint16_t, const TString&);
   ~JFFlucAnalysisO2Hist();
 };
 

PWGCF/JCorran/Tasks/jflucAnalysisTask.cxx

@@ -55,6 +55,7 @@ struct jflucAnalysisTask {
   O2_DEFINE_CONFIGURABLE(ptmin, float, 0.2, "Minimum pt for tracks");
   O2_DEFINE_CONFIGURABLE(ptmax, float, 5.0, "Maximum pt for tracks");
   O2_DEFINE_CONFIGURABLE(cfgCentBinsForMC, int, 0, "0 = OFF and 1 = ON for data like multiplicity/centrality bins for MC process");
+  O2_DEFINE_CONFIGURABLE(cfgMultCorrelationsMask, uint16_t, 0, "Selection bitmask for the multiplicity correlations. This should match the filter selection cfgEstimatorBitMask.")
 
   ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 100.1}, "multiplicity / centrality axis for histograms"};
   ConfigurableAxis phiAxis{"axisPhi", {50, 0.0, o2::constants::math::TwoPI}, "phi axis for histograms"};
@@ -78,15 +79,15 @@ struct jflucAnalysisTask {
     auto axisSpecZvt = AxisSpec(zvtAxis);
     auto axisSpecPt = AxisSpec(ptAxis);
     auto axisSpecMass = AxisSpec(massAxis);
-    if (doprocessJDerived || doprocessJDerivedCorrected || doprocessCFDerived || doprocessCFDerivedCorrected || doprocessMCCFDerived) {
-      pcf = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecPt, axisSpecMass, "jfluc");
+    if (doprocessJDerived || doprocessJDerivedCorrected || doprocessCFDerived || doprocessCFDerivedCorrected || doprocessCFDerivedMultSet || doprocessCFDerivedMultSetCorrected || doprocessMCCFDerived) {
+      pcf = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecPt, axisSpecMass, cfgMultCorrelationsMask, "jfluc");
       pcf->AddFlags(JFFlucAnalysis::kFlucEbEWeighting);
       pcf->UserCreateOutputObjects();
     } else {
       pcf = 0;
     }
     if (doprocessCF2ProngDerived || doprocessCF2ProngDerivedCorrected) {
-      pcf2Prong = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecPt, axisSpecMass, "jfluc2prong");
+      pcf2Prong = new JFFlucAnalysisO2Hist(registry, axisSpecMult, axisSpecPhi, axisSpecEta, axisSpecZvt, axisSpecPt, axisSpecMass, cfgMultCorrelationsMask, "jfluc2prong");
       pcf2Prong->AddFlags(JFFlucAnalysis::kFlucEbEWeighting);
       pcf2Prong->UserCreateOutputObjects();
 
@@ -95,6 +96,8 @@ struct jflucAnalysisTask {
     } else {
       pcf2Prong = 0;
     }
+    if ((doprocessCFDerivedMultSet || doprocessCFDerivedMultSetCorrected) && cfgMultCorrelationsMask == 0)
+      LOGF(fatal, "cfgMultCorrelationsMask can not be 0 when MultSet process functions are in use.");
   }
 
   template <class T>
@@ -107,6 +110,7 @@ struct jflucAnalysisTask {
     pcf->SetEventCentrality(collision.multiplicity());
     pcf->SetEventVertex(collision.posZ());
     pcf->FillQA(tracks);
+    pcf->FillMultSet(collision);
     qvecs.Calculate(tracks, etamin, etamax);
     pcf->SetJQVectors(&qvecs);
     pcf->UserExec("");
@@ -158,6 +162,22 @@ struct jflucAnalysisTask {
   }
   PROCESS_SWITCH(jflucAnalysisTask, processCFDerivedCorrected, "Process CF derived data with corrections", true);
 
+  void processCFDerivedMultSet(soa::Join<aod::CFCollisions, aod::CFMultSets>::iterator const& collision, soa::Filtered<aod::CFTracks> const& tracks)
+  {
+    if (std::popcount(cfgMultCorrelationsMask.value) != static_cast<int>(collision.multiplicities().size()))
+      LOGF(fatal, "Multiplicity selections (cfgMultCorrelationsMask = 0x%x) do not match the size of the table column (%ld). The histogram filling relies on the preservation of order.", cfgMultCorrelationsMask.value, collision.multiplicities().size());
+    analyze(collision, tracks);
+  }
+  PROCESS_SWITCH(jflucAnalysisTask, processCFDerivedMultSet, "Process CF derived data with multiplicity sets", false);
+
+  void processCFDerivedMultSetCorrected(soa::Join<aod::CFCollisions, aod::CFMultSets>::iterator const& collision, soa::Filtered<soa::Join<aod::CFTracks, aod::JWeights>> const& tracks)
+  {
+    if (std::popcount(cfgMultCorrelationsMask.value) != static_cast<int>(collision.multiplicities().size()))
+      LOGF(fatal, "Multiplicity selections (cfgMultCorrelationsMask = 0x%x) do not match the size of the table column (%ld). The histogram filling relies on the preservation of order.", cfgMultCorrelationsMask.value, collision.multiplicities().size());
+    analyze(collision, tracks);
+  }
+  PROCESS_SWITCH(jflucAnalysisTask, processCFDerivedMultSetCorrected, "Process CF derived data with corrections and multiplicity sets", false);
+
   void processCF2ProngDerived(aod::CFCollision const& collision, soa::Filtered<aod::CFTracks> const& tracks, soa::Filtered<aod::CF2ProngTracks> const& p2tracks)
   {
     analyze(collision, p2tracks, tracks);

PWGCF/TableProducer/filterCorrelations.cxx

@@ -64,13 +64,6 @@ struct FilterCF {
     kPIDProton = BIT(1)
   };
 
-  enum MultiplicityEstimators : uint8_t {
-    kCentFT0C = BIT(0),
-    kMultFV0A = BIT(1),
-    kMultNTracksPV = BIT(2),
-    kMultNTracksGlobal = BIT(3),
-  };
-
   // Configuration
   O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 7.0f, "Accepted z-vertex range")
   O2_DEFINE_CONFIGURABLE(cfgCutPt, float, 0.5f, "Minimal pT for tracks")
@@ -271,13 +264,13 @@ struct FilterCF {
 
     if constexpr (std::experimental::is_detected<HasMultTables, C1>::value) {
       multiplicities.clear();
-      if (cfgEstimatorBitMask & kCentFT0C)
+      if (cfgEstimatorBitMask & aod::cfmultset::CentFT0C)
         multiplicities.push_back(collision.centFT0C());
-      if (cfgEstimatorBitMask & kMultFV0A)
+      if (cfgEstimatorBitMask & aod::cfmultset::MultFV0A)
         multiplicities.push_back(collision.multFV0A());
-      if (cfgEstimatorBitMask & kMultNTracksPV)
+      if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksPV)
         multiplicities.push_back(collision.multNTracksPV());
-      if (cfgEstimatorBitMask & kMultNTracksGlobal)
+      if (cfgEstimatorBitMask & aod::cfmultset::MultNTracksGlobal)
         multiplicities.push_back(collision.multNTracksGlobal());
       outputMultSets(multiplicities);
     }