[PWGCF] Nch mc correlation matrix

PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx

@@ -686,8 +686,6 @@ struct FlowGenericFramework {
         registryQA.add("trackQA/after/Nch_uncorrected", "; N_{ch}; Counts", {HistType::kTH1D, {nchAxis}});
         registryQA.add("trackQA/after/etaNch", "; #eta; Counts", {HistType::kTH1D, {etaAxis}});
         registryQA.add("trackQA/after/etaPtPt", "; #eta; Counts", {HistType::kTH1D, {etaAxis}});
-        registryQA.add("trackQA/after/etaV02", "; #eta; Counts", {HistType::kTH1D, {etaAxis}});
-        registryQA.add("trackQA/after/etaV0", "; #eta; Counts", {HistType::kTH1D, {etaAxis}});
         if (!cfgFill.cfgFillRunByRunQA) {
           if (cfgUsePID) {
             registryQA.add<TH3>("phi_eta_vtxz_ref", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
@@ -732,6 +730,10 @@ struct FlowGenericFramework {
       AxisSpec axisLambdaMass = {resoSwitchVals[MassBins][Lambda], resoCutVals[MassMin][Lambda], resoCutVals[MassMax][Lambda]};
       AxisSpec yAxis = {100, -1, 1};
       // QA histograms for V0s
+      if (cfgFill.cfgFillV0QA && (resoSwitchVals[UseParticle][K0] != 0 || resoSwitchVals[UseParticle][Lambda] != 0)) {
+        registryQA.add("trackQA/after/etaV02", "; #eta; Counts", {HistType::kTH1D, {etaAxis}});
+        registryQA.add("trackQA/after/etaV0", "; #eta; Counts", {HistType::kTH1D, {etaAxis}});
+      }
       if (resoSwitchVals[UseParticle][K0] != 0) {
         if (cfgFill.cfgFillV0QA) {
           registryQA.add("K0/PiPlusTPC_K0", "", {HistType::kTH2D, {{ptAxis, axisNsigmaTPC}}});
@@ -2373,7 +2375,9 @@ struct FlowGenericFramework {
     registryQA.fill(HIST("K0/hK0Count"), FillDaughterTrackSelected);
     selection.selected = true;
     selection.isK0 = true;
-    registryQA.fill(HIST("K0/hK0AP"), v0.alpha(), v0.qtarm());
+    if (cfgFill.cfgFillV0QA) {
+      registryQA.fill(HIST("K0/hK0AP"), v0.alpha(), v0.qtarm());
+    }
     return selection;
   }
 
@@ -2481,15 +2485,15 @@ struct FlowGenericFramework {
       if (!selectionV0Daughter(postrack, Protons) || !selectionV0Daughter(negtrack, Pions)) {
         return selection;
       }
-      if (fillSelectionQA) {
+      if (fillSelectionQA && cfgFill.cfgFillV0QA) {
         registryQA.fill(HIST("Lambda/hLambdaAP"), v0.alpha(), v0.qtarm());
       }
     }
     if (selection.isAL) {
       if (!selectionV0Daughter(postrack, Pions) || !selectionV0Daughter(negtrack, Protons)) {
         return selection;
       }
-      if (fillSelectionQA) {
+      if (fillSelectionQA && cfgFill.cfgFillV0QA) {
         registryQA.fill(HIST("Lambda/hAntiLambdaAP"), v0.alpha(), v0.qtarm());
       }
     }
@@ -2827,11 +2831,15 @@ struct FlowGenericFramework {
 
     if (cfgEventSelection.cfgDoOccupancySel) {
       int occupancy = collision.trackOccupancyInTimeRange();
-      registryQA.fill(HIST("eventQA/before/occ_mult_cent"), occupancy, tracks.size(), centrality);
+      if (cfgFill.cfgFillQA) {
+        registryQA.fill(HIST("eventQA/before/occ_mult_cent"), occupancy, tracks.size(), centrality);
+      }
       if (occupancy < 0 || occupancy > cfgEventSelection.cfgOccupancySelection) {
         return;
       }
-      registryQA.fill(HIST("eventQA/after/occ_mult_cent"), occupancy, tracks.size(), centrality);
+      if (cfgFill.cfgFillQA) {
+        registryQA.fill(HIST("eventQA/after/occ_mult_cent"), occupancy, tracks.size(), centrality);
+      }
     }
     registryQA.fill(HIST("eventQA/eventSel"), 2.5);
     if (cfgFill.cfgFillRunByRunQA) {
@@ -2874,7 +2882,9 @@ struct FlowGenericFramework {
   void processOnTheFly(soa::Filtered<aod::McCollisions>::iterator const& mcCollision, aod::McParticles const& mcParticles, aod::V0Datas const& v0s)
   {
     int run = 0;
-    registryQA.fill(HIST("MCGen/impactParameter"), mcCollision.impactParameter(), mcParticles.size());
+    if (cfgFill.cfgFillQA) {
+      registryQA.fill(HIST("MCGen/impactParameter"), mcCollision.impactParameter(), mcParticles.size());
+    }
     processCollision<Gen>(mcCollision, mcParticles, v0s, mcCollision.impactParameter(), -999, run);
   }
   PROCESS_SWITCH(FlowGenericFramework, processOnTheFly, "Process analysis for MC on-the-fly generated events", false);
@@ -3269,7 +3279,9 @@ struct FlowGenericFramework {
         continue;
       }
       fillGeneratedEfficiencyTrack(particle, selectedCentrality);
-      fillGeneratedLambdaFeeddownXi(particle, selectedCentrality);
+      if (cfgFill.cfgFillV0QA) {
+        fillGeneratedLambdaFeeddownXi(particle, selectedCentrality);
+      }
       if (isGeneratedEfficiencyV0(particle, PDG_t::kK0Short, K0) && resoSwitchVals[UseParticle][K0] != 0) {
         fillGeneratedEfficiencyV0(particle, EfficiencyK0, selectedCentrality);
       }
@@ -3294,7 +3306,9 @@ struct FlowGenericFramework {
         if (v0.collisionId() != bestCollisionIndex) {
           continue;
         }
-        fillLambdaFeeddownReco(v0, collision, tracks, selectedCentrality);
+        if (cfgFill.cfgFillV0QA) {
+          fillLambdaFeeddownReco(v0, collision, tracks, selectedCentrality);
+        }
         fillEfficiencyRecoV0(v0, collision, tracks, selectedCentrality);
       }
       break;

PWGCF/GenericFramework/Tasks/flowGfwNonflow.cxx

@@ -44,6 +44,7 @@
 
 #include <TF1.h>
 #include <TH1.h>
+#include <TH2.h>
 #include <TH3.h>
 #include <TNamed.h>
 #include <TObjArray.h>
@@ -78,14 +79,15 @@ struct FlowGfwNonflow {
   Configurable<int> cfgMpar{"cfgMpar", 4, "Highest order of pt-pt correlations"};
   Configurable<int> cfgCentEstimator{"cfgCentEstimator", 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A, 4:NTPV, 5:NGlobal, 6:MFT"};
   Configurable<bool> cfgUseNch{"cfgUseNch", false, "Do correlations as function of Nch"};
-  Configurable<int> cfgUseNchCorrection{"cfgUseNchCorrection", 1, "Use correction for Nch; 0: Use size of tracks table, 1: Use efficiency-corrected Nch values, 2: Use uncorrected Nch values"};
+  Configurable<int> cfgUseNchCorrection{"cfgUseNchCorrection", 1, "Nch used on the x-axis; 0: tracks table size, 1: efficiency-corrected, 2: accepted reconstructed, 3: Reco.-gen. response-matrix corrected"};
   Configurable<bool> cfgRunByRun{"cfgRunByRun", false, "Use run-by-run NUA"};
   Configurable<bool> cfgFillQA{"cfgFillQA", false, "Fill QA histograms"};
   Configurable<bool> cfgUseCentralMoments{"cfgUseCentralMoments", true, "Use central moments in vn-pt calculations"};
   Configurable<bool> cfgUseMultiplicityFlowWeights{"cfgUseMultiplicityFlowWeights", true, "Enable or disable the use of multiplicity-based event weighting"};
   struct : ConfigurableGroup {
     Configurable<std::string> cfgEfficiencyPath{"cfgEfficiencyPath", "", "CCDB path to efficiency object"};
     Configurable<bool> cfgUse2DEfficiency{"cfgUse2DEfficiency", false, "Toggle the use of 2D (pt, centrality) efficiency versus centrality integrated efficiency"};
+    Configurable<std::string> cfgNchResponsePath{"cfgNchResponsePath", "", "CCDB path to TH2 response matrix (reconstructed Nch on x, generated Nch on y)"};
     Configurable<std::string> cfgAcceptancePath{"cfgAcceptancePath", "", "CCDB path to acceptance object"};
   } cfgCorrections;
   struct : ConfigurableGroup {
@@ -181,6 +183,7 @@ struct FlowGfwNonflow {
 
   struct Config {
     TH1* mEfficiency = nullptr;
+    TH2* mNchResponse = nullptr;
     std::vector<GFWWeights*> mAcceptance;
     bool correctionsLoaded = false;
   } correctionsConfig;
@@ -240,6 +243,12 @@ struct FlowGfwNonflow {
     ProtonID,
     SpeciesCount
   };
+  enum NchSelector {
+    TableSize,
+    Corrected,
+    Uncorrected,
+    ResponseMatrixCorrected
+  };
 
   // Generic Framework
   GFW* fGFW = new GFW();
@@ -379,6 +388,9 @@ struct FlowGfwNonflow {
         registry.add("eventQA/before/occ_mult_cent", "; occupancy; N_{ch}; centrality (%)", {HistType::kTH3D, {occAxis, nchAxis, centAxis}});
       }
     }
+    if (doprocessMCReco) {
+      registry.add("MCReco/Nch_reco_gen", "; N_{ch}^{reco}; N_{ch}^{gen}", {HistType::kTH2D, {nchAxis, nchAxis}});
+    }
     registry.add("eventQA/before/centrality", "; centrality (%); Counts", {HistType::kTH1D, {centAxis}});
     registry.add("eventQA/before/multiplicity", "; N_{ch}; Counts", {HistType::kTH1D, {nchAxis}});
     registry.addClone("eventQA/before/", "eventQA/after/");
@@ -639,9 +651,43 @@ struct FlowGfwNonflow {
       }
       LOGF(info, "Loaded efficiency histogram from %s", cfgCorrections.cfgEfficiencyPath.value.c_str());
     }
+    if (!cfgCorrections.cfgNchResponsePath.value.empty()) {
+      correctionsConfig.mNchResponse = ccdb->getForTimeStamp<TH2D>(cfgCorrections.cfgNchResponsePath, timestamp);
+      if (correctionsConfig.mNchResponse == nullptr) {
+        LOGF(fatal, "Could not load Nch response matrix from %s", cfgCorrections.cfgNchResponsePath.value.c_str());
+      }
+      LOGF(info, "Loaded Nch response matrix from %s", cfgCorrections.cfgNchResponsePath.value.c_str());
+    } else if (cfgUseNchCorrection == NchSelector::ResponseMatrixCorrected) {
+      LOGF(fatal, "cfgUseNchCorrection=3 requires cfgNchResponsePath");
+    }
     correctionsConfig.correctionsLoaded = true;
   }
 
+  float getResponseCorrectedNch(const unsigned int multReconstructed) const
+  {
+    if (!correctionsConfig.mNchResponse) {
+      return multReconstructed;
+    }
+    const auto* response = correctionsConfig.mNchResponse;
+    const int recoBin = response->GetXaxis()->FindFixBin(multReconstructed);
+    if (recoBin < 1 || recoBin > response->GetNbinsX()) {
+      LOGF(warn, "Reconstructed Nch %u is outside the response matrix; using the uncorrected value", multReconstructed);
+      return multReconstructed;
+    }
+    double sumWeights = 0.;
+    double sumGeneratedNch = 0.;
+    for (int genBin = 1; genBin <= response->GetNbinsY(); ++genBin) {
+      const double weight = response->GetBinContent(recoBin, genBin);
+      sumWeights += weight;
+      sumGeneratedNch += weight * response->GetYaxis()->GetBinCenter(genBin);
+    }
+    if (sumWeights <= 0.) {
+      LOGF(warn, "Response matrix has no entries for reconstructed Nch %u; using the uncorrected value", multReconstructed);
+      return multReconstructed;
+    }
+    return sumGeneratedNch / sumWeights;
+  }
+
   template <typename TTrack>
   double getAcceptance(const TTrack& track, const double& vtxz)
   { // 0 ref, 1 ch, 2 pi, 3 ka, 4 pr
@@ -677,22 +723,22 @@ struct FlowGfwNonflow {
         return -1.;
       }
       return 1. / eff;
-    } else {
-      auto* effHist = dynamic_cast<TH1D*>(correctionsConfig.mEfficiency);
-      if (!effHist) {
-        LOGF(error, "Efficiency object at %s is not a TH1D", cfgCorrections.cfgEfficiencyPath.value.c_str());
-        return -1.;
-      }
-      bin = effHist->FindBin(track.pt());
-      if (!bin) {
-        return -1.;
-      }
-      const double eff = effHist->GetBinContent(bin);
-      if (!std::isfinite(eff) || eff <= 0.) {
-        return -1.;
-      }
-      return 1. / eff;
     }
+
+    auto* effHist = dynamic_cast<TH1D*>(correctionsConfig.mEfficiency);
+    if (!effHist) {
+      LOGF(error, "Efficiency object at %s is not a TH1D", cfgCorrections.cfgEfficiencyPath.value.c_str());
+      return -1.;
+    }
+    bin = effHist->FindBin(track.pt());
+    if (!bin) {
+      return -1.;
+    }
+    const double eff = effHist->GetBinContent(bin);
+    if (!std::isfinite(eff) || eff <= 0.) {
+      return -1.;
+    }
+    return 1. / eff;
   }
 
   template <typename TCollision>
@@ -1014,7 +1060,7 @@ struct FlowGfwNonflow {
   };
 
   template <DataType dt, typename TCollision, typename TTracks>
-  void processCollision(const TCollision& collision, const TTracks& tracks, const float& centrality, const float& field)
+  void processCollision(const TCollision& collision, const TTracks& tracks, const float& centrality, const float& field, const int generatedNch = -1)
   {
     if (tracks.size() < 1) {
       return;
@@ -1037,22 +1083,30 @@ struct FlowGfwNonflow {
     for (const auto& track : tracks) {
       processTrack(track, vtxz, field, centrality, acceptedTracks);
     }
+    if constexpr (dt == Reco) {
+      if (generatedNch >= 0) {
+        registry.fill(HIST("MCReco/Nch_reco_gen"), acceptedTracks.totaluncorr, generatedNch);
+      }
+    }
     if (dt != Gen && cfgFillQA) {
       registry.fill(HIST("trackQA/after/Nch_corrected"), acceptedTracks.total);
       registry.fill(HIST("trackQA/after/Nch_uncorrected"), acceptedTracks.totaluncorr);
     }
 
     float multiplicity = 0.f;
     switch (cfgUseNchCorrection) {
-      case 0:
+      case NchSelector::TableSize:
         multiplicity = tracks.size();
         break;
-      case 1:
+      case NchSelector::Corrected:
         multiplicity = acceptedTracks.total;
         break;
-      case 2:
+      case NchSelector::Uncorrected:
         multiplicity = acceptedTracks.totaluncorr;
         break;
+      case NchSelector::ResponseMatrixCorrected:
+        multiplicity = (dt == Gen) ? acceptedTracks.totaluncorr : getResponseCorrectedNch(acceptedTracks.totaluncorr);
+        break;
       default:
         multiplicity = tracks.size();
         break;
@@ -1155,8 +1209,10 @@ struct FlowGfwNonflow {
 
   using GFWCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs>>;
   using GFWMCCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs, aod::McCollisionLabels>;
+  using FilteredGFWMCCollisions = soa::Filtered<GFWMCCollisions>;
   using GFWTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA>>;
   using GFWMCTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::McTrackLabels>>;
+  Preslice<aod::McParticles> particlesPerMcCollision = aod::mcparticle::mcCollisionId;
 
   SliceCache cache;
   Partition<GFWTracks> posTracks = aod::track::signed1Pt > 0.0f;
@@ -1221,7 +1277,7 @@ struct FlowGfwNonflow {
   }
   PROCESS_SWITCH(FlowGfwNonflow, processData, "Process analysis for non-derived data", true);
 
-  void processMCReco(GFWCollisions::iterator const& collision, aod::BCsWithTimestamps const&, GFWMCTracks const& tracks, aod::McParticles const&)
+  void processMCReco(FilteredGFWMCCollisions::iterator const& collision, aod::BCsWithTimestamps const&, GFWMCTracks const& tracks, aod::McParticles const& particles)
   {
     auto bc = collision.bc_as<aod::BCsWithTimestamps>();
     int run = bc.runNumber();
@@ -1265,9 +1321,16 @@ struct FlowGfwNonflow {
     if (cfgFillQA) {
       fillEventQA<After>(collision, tracks);
     }
+    unsigned int generatedNch = 0;
+    const auto particlesThisCollision = particles.sliceBy(particlesPerMcCollision, collision.mcCollisionId());
+    for (const auto& particle : particlesThisCollision) {
+      if (particle.isPhysicalPrimary() && particle.eta() > cfgKinematics.cfgEtaNch->first && particle.eta() < cfgKinematics.cfgEtaNch->second && particle.pt() > gfwMemberCache.ptlow && particle.pt() < gfwMemberCache.ptup) {
+        ++generatedNch;
+      }
+    }
     loadCorrections(bc);
     auto field = (cfgEventSelection.cfgMagField == DefaultMagneticFieldCut) ? getMagneticField(bc.timestamp()) : static_cast<int>(cfgEventSelection.cfgMagField);
-    processCollision<Reco>(collision, tracks, centrality, field);
+    processCollision<Reco>(collision, tracks, centrality, field, generatedNch);
   }
   PROCESS_SWITCH(FlowGfwNonflow, processMCReco, "Process analysis for MC reconstructed events", false);