[PWGUD] MC analysis for flow

PWGUD/Tasks/CMakeLists.txt

@@ -14,27 +14,27 @@
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(upc
                            SOURCES upcAnalysis.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-spectra
                            SOURCES sgSpectraAnalyzer.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-pid-spectra-table
                            SOURCES sgPIDSpectraTable.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-pid-analyzer
                            SOURCES sgPIDAnalyzer.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-pid-spectra
                            SOURCES sgPIDSpectra.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
@@ -49,27 +49,27 @@
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-excl-universe
                            SOURCES sgExcUniverse.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-fourpi
                            SOURCES sgFourPiAnalyzer.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-sixpi
                            SOURCES sgSixPiAnalyzer.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-twopi
                            SOURCES sgTwoPiAnalyzer.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(sg-d0
                            SOURCES sgD0Analyzer.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase
                            COMPONENT_NAME Analysis)
@@ -254,6 +254,11 @@
 			               PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGCFCore
                            COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(flow-mc-upc
+                           SOURCES flowMcUpc.cxx
+			               PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore
+                           COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(analysis-mc-dpm-jet-sg-v3
                            SOURCES analysisMCDPMJetSGv3.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
@@ -272,4 +277,4 @@
 o2physics_add_dpl_workflow(fitbit-mapping
 													 SOURCES upcTestFITBitMapping.cxx
 													 PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
-													 COMPONENT_NAME Analysis)
+													 COMPONENT_NAME Analysis)

PWGUD/Tasks/flowCumulantsUpc.cxx

@@ -125,9 +125,18 @@ struct FlowCumulantsUpc {
   O2_DEFINE_CONFIGURABLE(cfgDcazCut, float, 10.0, "dcaz cut")
   O2_DEFINE_CONFIGURABLE(cfgItsClusterSize, unsigned int, 5, "ITS cluster size")
   O2_DEFINE_CONFIGURABLE(cfgMaxTPCChi2NCl, int, 4, "tpcchi2")
+  O2_DEFINE_CONFIGURABLE(cfgConsistentEventFlag, int, 0, "Flag to select consistent events - 0: off, 1: v2{2} gap calculable, 2: v2{4} full calculable, 4: v2{4} gap calculable, 8: v2{4} 3sub calculable")
+
   Configurable<std::vector<std::string>> cfgUserDefineGFWCorr{"cfgUserDefineGFWCorr", std::vector<std::string>{"refN02 {2} refP02 {-2}", "refN12 {2} refP12 {-2}"}, "User defined GFW CorrelatorConfig"};
   Configurable<std::vector<std::string>> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector<std::string>{"Ch02Gap22", "Ch12Gap22"}, "User defined GFW Name"};
   Configurable<std::vector<int>> cfgRunRemoveList{"cfgRunRemoveList", std::vector<int>{-1}, "excluded run numbers"};
+  Configurable<std::vector<float>> cfgConsistentEventVector{"cfgConsistentEventVector", std::vector<float>{-0.8, -0.5, -0.4, 0.4, 0.5, 0.8}, "eta regions: left(min,max), mid(min,max), right(min,max)"};
+  struct AcceptedTracks {
+    int nNeg;
+    int nMid;
+    int nPos;
+    int nFull;
+  };
 
   ConfigurableAxis axisPtHist{"axisPtHist", {100, 0., 10.}, "pt axis for histograms"};
   ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.5, 4, 5, 6, 8, 10}, "pt axis for histograms"};
@@ -217,13 +226,15 @@ struct FlowCumulantsUpc {
 
     // Add some output objects to the histogram registry
     // Event QA
-    registry.add("hEventCount", "Number of Event;; Count", {HistType::kTH1D, {{5, 0, 5}}});
+    registry.add("hEventCount", "Number of Event;; Count", {HistType::kTH1D, {{6, 0, 6}}});
     registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(1, "Filtered event");
     registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "after gapside selection");
     registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(3, "after its selection");
     registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(4, "after pt selection");
     registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(5, "after occupancy");
-    registry.add("hTrackCount", "Number of tracks;; Count", {HistType::kTH1D, {{5, 0, 5}}});
+    registry.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(6, "after consistency check");
+
+    registry.add("hTrackCount", "Number of tracks;; Count", {HistType::kTH1D, {{7, 0, 7}}});
     registry.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(1, "after event selection");
     registry.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(2, "PVContributor");
     registry.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(3, "dcaz");
@@ -955,7 +966,7 @@ struct FlowCumulantsUpc {
     registry.fill(HIST("hMult"), tracks.size());
     registry.fill(HIST("hCent"), cent);
     fGFW->Clear();
-    if (cfgIfVertex && abs(vtxz) > cfgCutVertex) {
+    if (cfgIfVertex && std::abs(vtxz) > cfgCutVertex) {
       return;
     }
     registry.fill(HIST("hEventCount"), 3.5);
@@ -971,6 +982,8 @@ struct FlowCumulantsUpc {
     if (cfgUseNch) {
       independent = static_cast<float>(tracks.size());
     }
+    AcceptedTracks acceptedTracks{0, 0, 0, 0};
+    std::vector<float> consistentEventVector = cfgConsistentEventVector;
 
     for (const auto& track : tracks) {
       registry.fill(HIST("hChi2prTPCcls"), track.tpcChi2NCl());
@@ -996,6 +1009,16 @@ struct FlowCumulantsUpc {
         continue;
       }
       registry.fill(HIST("hPt"), track.pt());
+
+      if (cfgConsistentEventFlag && consistentEventVector.size() == 6) { // o2-linter: disable=magic-number (size match)
+        acceptedTracks.nFull += 1;
+        if (eta > consistentEventVector[0] && eta < consistentEventVector[1])
+          acceptedTracks.nNeg += 1;
+        if (eta > consistentEventVector[2] && eta < consistentEventVector[3])
+          acceptedTracks.nMid += 1;
+        if (eta > consistentEventVector[4] && eta < consistentEventVector[5])
+          acceptedTracks.nPos += 1;
+      }
       if (withinPtRef) {
         registry.fill(HIST("hPhi"), phi);
         registry.fill(HIST("hPhiWeighted"), phi, wacc);
@@ -1017,6 +1040,23 @@ struct FlowCumulantsUpc {
       registry.fill(HIST("hEtaNch2D"), eta, tracks.size());
     }
     registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected);
+    if (cfgConsistentEventFlag) {
+      if (cfgConsistentEventFlag & 1) {
+        if (!acceptedTracks.nPos || !acceptedTracks.nNeg)
+          return;
+      } else if (cfgConsistentEventFlag & 2) {
+        if (acceptedTracks.nFull < 4) // o2-linter: disable=magic-number (at least four tracks in full acceptance)
+          return;
+      } else if (cfgConsistentEventFlag & 4) {
+        if (acceptedTracks.nPos < 2 || acceptedTracks.nNeg < 2) // o2-linter: disable=magic-number (at least two tracks in each subevent)
+          return;
+      }
+      if (cfgConsistentEventFlag & 8) {
+        if (acceptedTracks.nPos < 2 || acceptedTracks.nMid < 2 || acceptedTracks.nNeg < 2) // o2-linter: disable=magic-number (at least two tracks in all three subevents)
+          return;
+      }
+    }
+    registry.fill(HIST("hEventCount"), 5.5);
 
     // Filling Flow Container
     for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) {

PWGUD/Tasks/flowMcUpc.cxx

@@ -0,0 +1,193 @@
+// 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.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// 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.
+
+/// \file   flowMcUpc.cxx
+/// \author Zhiyong Lu (zhiyong.lu@cern.ch), Yongxi Du (yongxi.du@cern.ch)
+/// \since  Apr/2/2026
+/// \brief  flow efficiency analysis on UPC MC
+
+#include "PWGUD/Core/SGSelector.h"
+#include "PWGUD/DataModel/UDTables.h"
+
+#include "Common/Core/RecoDecay.h"
+#include "Common/Core/TrackSelection.h"
+#include "Common/Core/TrackSelectionDefaults.h"
+#include "Common/Core/trackUtilities.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/RunningWorkflowInfo.h"
+#include "Framework/runDataProcessing.h"
+#include "ReconstructionDataFormats/Track.h"
+#include <CCDB/BasicCCDBManager.h>
+
+#include <TF1.h>
+#include <TPDGCode.h>
+#include <TProfile.h>
+#include <TRandom3.h>
+
+#include <string>
+#include <vector>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
+
+struct FlowMcUpc {
+  HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
+
+  Configurable<float> minB{"minB", 0.0f, "min impact parameter"};
+  Configurable<float> maxB{"maxB", 20.0f, "max impact parameter"};
+  O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
+  O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
+  O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.1f, "Minimal pT for tracks")
+  O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 1000.0f, "Maximal pT for tracks")
+  O2_DEFINE_CONFIGURABLE(cfgCutDCAxy, float, 0.2f, "DCAxy cut for tracks")
+  O2_DEFINE_CONFIGURABLE(cfgDcaxy, bool, true, "choose dcaxy")
+
+  ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""};
+  ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f}, "pt axis"};
+  // Connect to ccdb
+  Service<ccdb::BasicCCDBManager> ccdb;
+  Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+
+  double epsilon = 1e-6;
+
+  using McParticles = soa::Join<aod::UDMcParticles, aod::UDMcTrackLabels>;
+
+  void init(InitContext&)
+  {
+    ccdb->setURL(ccdbUrl.value);
+    ccdb->setCaching(true);
+    auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+    ccdb->setCreatedNotAfter(now);
+
+    const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
+    const AxisSpec axisEta{20, -1., 1., "#eta"};
+    const AxisSpec axisCounter{1, 0, +1, ""};
+    // QA histograms
+    histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {{5, 0, 5}});
+    histos.add<TH1>("RecoProcessEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}});
+    histos.add<TH1>("hImpactParameter", "hImpactParameter", HistType::kTH1D, {axisB});
+
+    histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
+    histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
+    histos.add<TH1>("hPtReco", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
+    histos.add<TH3>("hEtaPtVtxzMCReco", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
+  }
+
+  // template <typename TCollision>
+  // bool eventSelected(TCollision collision)
+  // {
+  //   return true;
+  // }
+
+  template <typename TTrack>
+  bool trackSelected(TTrack track)
+  {
+    auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
+    double pt = RecoDecay::pt(momentum);
+    if (pt < cfgPtCutMin || pt > cfgPtCutMax) {
+      return false;
+    }
+    double dcaLimit = 0.0105 + 0.035 / std::pow(pt, 1.1);
+    if (cfgDcaxy && !(std::fabs(track.dcaXY()) < dcaLimit)) {
+      return false;
+    }
+    return true;
+  }
+
+  void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, McParticles const& mcParticles, aod::BCs const& bcs)
+  {
+    if (bcs.size() == 0) {
+      return;
+    }
+    histos.fill(HIST("mcEventCounter"), 0.5);
+    float imp = mcCollision.impactParameter();
+    float vtxz = mcCollision.posZ();
+
+    if (imp >= minB && imp <= maxB) {
+      // event within range
+      histos.fill(HIST("hImpactParameter"), imp);
+
+      for (auto const& mcParticle : mcParticles) {
+        auto momentum = std::array<double, 3>{mcParticle.px(), mcParticle.py(), mcParticle.pz()};
+        double pt = RecoDecay::pt(momentum);
+        // double phi = RecoDecay::phi(momentum);
+        double eta = RecoDecay::eta(momentum);
+        // focus on bulk: e, mu, pi, k, p
+        int pdgCode = std::abs(mcParticle.pdgCode());
+        if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
+          continue;
+
+        if (!mcParticle.isPhysicalPrimary())
+          continue;
+        // if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
+        //   continue;
+
+        histos.fill(HIST("hPtMCGen"), pt);
+        histos.fill(HIST("hEtaPtVtxzMCGen"), eta, pt, vtxz);
+      }
+    }
+  }
+  PROCESS_SWITCH(FlowMcUpc, processMCTrue, "process pure simulation information", true);
+
+  using MCRecoTracks = soa::Join<aod::UDTracks, aod::UDTracksPID, aod::UDTracksExtra, aod::UDTracksFlags, aod::UDTracksDCA, aod::UDMcTrackLabels>;
+  using MCRecoCollisions = soa::Join<aod::UDCollisions, aod::SGCollisions, aod::UDCollisionSelExtras, aod::UDCollisionsSels, aod::UDZdcsReduced, aod::UDMcCollsLabels>;
+
+  void processReco(MCRecoCollisions::iterator const& collision, MCRecoTracks const& tracks)
+  {
+    histos.fill(HIST("RecoProcessEventCounter"), 0.5);
+    // if (!eventSelected(collision))
+    //   return;
+    histos.fill(HIST("RecoProcessEventCounter"), 1.5);
+    if (!collision.has_udMcCollision())
+      return;
+    histos.fill(HIST("RecoProcessEventCounter"), 2.5);
+    if (tracks.size() < 1)
+      return;
+    histos.fill(HIST("RecoProcessEventCounter"), 3.5);
+
+    float vtxz = collision.posZ();
+
+    for (const auto& track : tracks) {
+      auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
+      double pt = RecoDecay::pt(momentum);
+      // double phi = RecoDecay::phi(momentum);
+      double eta = RecoDecay::eta(momentum);
+      if (!trackSelected(track) || (!track.has_udMcParticle()))
+        continue;
+      auto mcParticle = track.udMcParticle();
+      int pdgCode = std::abs(mcParticle.pdgCode());
+      if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
+        continue;
+      // if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
+      //   continue;
+      if (!mcParticle.isPhysicalPrimary())
+        continue;
+
+      histos.fill(HIST("hPtReco"), pt);
+      histos.fill(HIST("hEtaPtVtxzMCReco"), eta, pt, vtxz);
+    }
+  }
+  PROCESS_SWITCH(FlowMcUpc, processReco, "process reconstructed information", true);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<FlowMcUpc>(cfgc)};
+}