From d616170a841e8b1ee1a1e5a3a4a8aa00dbae1c3d Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Thu, 2 Apr 2026 12:52:23 +0800 Subject: [PATCH 1/8] add event consistency flag --- PWGUD/Tasks/flowCumulantsUpc.cxx | 96 +++++++++++++++++++++----------- 1 file changed, 63 insertions(+), 33 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index ad75ae1a3c3..49ef99e88e4 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -14,51 +14,41 @@ /// \since Mar/2025 /// \brief jira: , task to measure flow observables with cumulant method +#include "FlowContainer.h" +#include "GFW.h" +#include "GFWCumulant.h" +#include "GFWPowerArray.h" +#include "GFWWeights.h" + #include "PWGUD/Core/SGSelector.h" +#include "PWGUD/DataModel/SGTables.h" #include "PWGUD/DataModel/UDTables.h" -// -#include "PWGCF/GenericFramework/Core/FlowContainer.h" -#include "PWGCF/GenericFramework/Core/GFW.h" -#include "PWGCF/GenericFramework/Core/GFWWeights.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" - +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include + +#include "TList.h" +#include "TVector3.h" #include -#include -#include -#include #include #include #include -#include - -#include - -#include #include -#include #include -#include #include #include #include @@ -125,9 +115,17 @@ 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> cfgUserDefineGFWCorr{"cfgUserDefineGFWCorr", std::vector{"refN02 {2} refP02 {-2}", "refN12 {2} refP12 {-2}"}, "User defined GFW CorrelatorConfig"}; Configurable> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector{"Ch02Gap22", "Ch12Gap22"}, "User defined GFW Name"}; Configurable> cfgRunRemoveList{"cfgRunRemoveList", std::vector{-1}, "excluded run numbers"}; + Configurable> cfgConsistentEventVector{"cfgConsistentEventVector", std::vector{-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 +215,14 @@ 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(HIST("hEventCount"))->GetXaxis()->SetBinLabel(1, "Filtered event"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "after gapside selection"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(3, "after its selection"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(4, "after pt selection"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(5, "after occupancy"); - registry.add("hTrackCount", "Number of tracks;; Count", {HistType::kTH1D, {{5, 0, 5}}}); + registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(6, "after consistency check"); + registry.add("hTrackCount", "Number of tracks;; Count", {HistType::kTH1D, {{7, 0, 7}}}); registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(1, "after event selection"); registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(2, "PVContributor"); registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(3, "dcaz"); @@ -955,7 +954,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 +970,8 @@ struct FlowCumulantsUpc { if (cfgUseNch) { independent = static_cast(tracks.size()); } + AcceptedTracks acceptedTracks{0, 0, 0, 0}; + std::vector consistentEventVector = cfgConsistentEventVector; for (const auto& track : tracks) { registry.fill(HIST("hChi2prTPCcls"), track.tpcChi2NCl()); @@ -996,6 +997,17 @@ 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); @@ -1018,6 +1030,24 @@ struct FlowCumulantsUpc { } 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++) { fillFC(corrconfigs.at(l_ind), independent, lRandom); From 913605ee5cae54f8762f56482ae9ab1849e5128b Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Thu, 2 Apr 2026 17:33:37 +0800 Subject: [PATCH 2/8] add event consistency flag --- PWGUD/Tasks/flowCumulantsUpc.cxx | 54 +++++++++++++++++++------------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 49ef99e88e4..9624c1a5c86 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -14,41 +14,51 @@ /// \since Mar/2025 /// \brief jira: , task to measure flow observables with cumulant method -#include "FlowContainer.h" -#include "GFW.h" -#include "GFWCumulant.h" -#include "GFWPowerArray.h" -#include "GFWWeights.h" - #include "PWGUD/Core/SGSelector.h" -#include "PWGUD/DataModel/SGTables.h" #include "PWGUD/DataModel/UDTables.h" +// +#include "PWGCF/GenericFramework/Core/FlowContainer.h" +#include "PWGCF/GenericFramework/Core/GFW.h" +#include "PWGCF/GenericFramework/Core/GFWWeights.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisTask.h" -#include "Framework/HistogramRegistry.h" -#include "Framework/RunningWorkflowInfo.h" -#include "Framework/runDataProcessing.h" -#include -#include "TList.h" -#include "TVector3.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include #include +#include +#include +#include #include #include #include +#include + +#include + +#include #include +#include #include +#include #include #include #include @@ -116,6 +126,7 @@ struct FlowCumulantsUpc { 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> cfgUserDefineGFWCorr{"cfgUserDefineGFWCorr", std::vector{"refN02 {2} refP02 {-2}", "refN12 {2} refP12 {-2}"}, "User defined GFW CorrelatorConfig"}; Configurable> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector{"Ch02Gap22", "Ch12Gap22"}, "User defined GFW Name"}; Configurable> cfgRunRemoveList{"cfgRunRemoveList", std::vector{-1}, "excluded run numbers"}; @@ -222,6 +233,7 @@ struct FlowCumulantsUpc { registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(4, "after pt selection"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(5, "after occupancy"); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(6, "after consistency check"); + registry.add("hTrackCount", "Number of tracks;; Count", {HistType::kTH1D, {{7, 0, 7}}}); registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(1, "after event selection"); registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(2, "PVContributor"); @@ -1007,7 +1019,6 @@ struct FlowCumulantsUpc { if (eta > consistentEventVector[4] && eta < consistentEventVector[5]) acceptedTracks.nPos += 1; } - if (withinPtRef) { registry.fill(HIST("hPhi"), phi); registry.fill(HIST("hPhiWeighted"), phi, wacc); @@ -1029,7 +1040,6 @@ 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) From 7290de79c16209828cb328e5dae0b34481debb6c Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Thu, 2 Apr 2026 21:19:02 +0800 Subject: [PATCH 3/8] MC analysis for flow --- PWGUD/Tasks/CMakeLists.txt | 7 +- PWGUD/Tasks/flowMcUpc.cxx | 193 +++++++++++++++++++++++++++++++++++++ 2 files changed, 199 insertions(+), 1 deletion(-) create mode 100644 PWGUD/Tasks/flowMcUpc.cxx diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index d8aa575b906..ab54549251c 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -254,6 +254,11 @@ o2physics_add_dpl_workflow(flow-correlations-upc 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(sg-exclusive-jpsi-midrapidity o2physics_add_dpl_workflow(fitbit-mapping SOURCES upcTestFITBitMapping.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore - COMPONENT_NAME Analysis) \ No newline at end of file + COMPONENT_NAME Analysis) diff --git a/PWGUD/Tasks/flowMcUpc.cxx b/PWGUD/Tasks/flowMcUpc.cxx new file mode 100644 index 00000000000..457367cf97d --- /dev/null +++ b/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 + +#include +#include +#include +#include + +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; + +struct FlowMcUpc { + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Configurable minB{"minB", 0.0f, "min impact parameter"}; + Configurable 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; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + double epsilon = 1e-6; + + using McParticles = soa::Join; + + void init(InitContext&) + { + ccdb->setURL(ccdbUrl.value); + ccdb->setCaching(true); + auto now = std::chrono::duration_cast(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("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {{5, 0, 5}}); + histos.add("RecoProcessEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}}); + histos.add("hImpactParameter", "hImpactParameter", HistType::kTH1D, {axisB}); + + histos.add("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); + histos.add("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}}); + histos.add("hPtReco", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); + histos.add("hEtaPtVtxzMCReco", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}}); + } + + // template + // bool eventSelected(TCollision collision) + // { + // return true; + // } + + template + bool trackSelected(TTrack track) + { + auto momentum = std::array{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{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; + using MCRecoCollisions = soa::Join; + + 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{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(cfgc)}; +} From 9c3c8206f727db2f45132abd9f66029d9fbbcd18 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Wed, 8 Apr 2026 14:32:14 +0800 Subject: [PATCH 4/8] use get track size --- PWGUD/Tasks/flowCorrelationsUpc.cxx | 9 ++++ PWGUD/Tasks/flowMcUpc.cxx | 77 ++++++++++++++++++++--------- 2 files changed, 64 insertions(+), 22 deletions(-) diff --git a/PWGUD/Tasks/flowCorrelationsUpc.cxx b/PWGUD/Tasks/flowCorrelationsUpc.cxx index 061eda25566..bebb14f1158 100644 --- a/PWGUD/Tasks/flowCorrelationsUpc.cxx +++ b/PWGUD/Tasks/flowCorrelationsUpc.cxx @@ -382,6 +382,15 @@ struct FlowCorrelationsUpc { // the process for filling the mixed events void processMixed(UDCollisionsFull const& collisions, UdTracksFull const& tracks) { + auto getTracksSize = [&tracks, this](UDCollisionsFull::iterator const& collision) { + auto associatedTracks = tracks.sliceByCached(o2::aod::udtrack::udCollisionId, collision.udCollisionId(), this->cache); + + auto mult = associatedTracks.size(); + return mult; + }; + + // using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; + // MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default) auto tracksTuple = std::make_tuple(tracks); SameKindPair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip diff --git a/PWGUD/Tasks/flowMcUpc.cxx b/PWGUD/Tasks/flowMcUpc.cxx index 457367cf97d..ffa9455d0dd 100644 --- a/PWGUD/Tasks/flowMcUpc.cxx +++ b/PWGUD/Tasks/flowMcUpc.cxx @@ -23,16 +23,17 @@ #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 +#include +#include +#include +#include +#include +#include +#include +#include +#include -#include #include #include #include @@ -43,6 +44,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using LorentzVectorM = ROOT::Math::LorentzVector>; #define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; @@ -96,10 +98,10 @@ struct FlowMcUpc { // } template - bool trackSelected(TTrack track) + bool trackSelected(TTrack const& track) { - auto momentum = std::array{track.px(), track.py(), track.pz()}; - double pt = RecoDecay::pt(momentum); + // auto momentum = std::array{track.px(), track.py(), track.pz()}; + double pt = track.pt(); if (pt < cfgPtCutMin || pt > cfgPtCutMax) { return false; } @@ -124,15 +126,28 @@ struct FlowMcUpc { histos.fill(HIST("hImpactParameter"), imp); for (auto const& mcParticle : mcParticles) { - auto momentum = std::array{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 + // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; + LorentzVectorM pMC(mcParticle.px(), mcParticle.py(), mcParticle.pz(), 0); // double phi = RecoDecay::phi(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) + if (pdgCode == PDG_t::kElectron) + pMC.SetM(o2::constants::physics::MassElectron); + else if (pdgCode == PDG_t::kMuonMinus) + pMC.SetM(o2::constants::physics::MassMuon); + else if (pdgCode == PDG_t::kPiPlus) + pMC.SetM(o2::constants::physics::MassPionCharged); + else if (pdgCode == PDG_t::kKPlus) + pMC.SetM(o2::constants::physics::MassKaonCharged); + else if (pdgCode == PDG_t::kProton) + pMC.SetM(o2::constants::physics::MassProton); + else continue; + double pt = pMC.Pt(); + double eta = pMC.Eta(); + + // if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton) + // continue; + if (!mcParticle.isPhysicalPrimary()) continue; // if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance @@ -164,16 +179,34 @@ struct FlowMcUpc { float vtxz = collision.posZ(); for (const auto& track : tracks) { - auto momentum = std::array{track.px(), track.py(), track.pz()}; - double pt = RecoDecay::pt(momentum); - double phi = RecoDecay::phi(momentum); - double eta = RecoDecay::eta(momentum); + LorentzVectorM recoMC(track.px(), track.py(), track.pz(), 0); + // double phi = RecoDecay::phi(momentum); + // focus on bulk: e, mu, pi, k, p + // auto momentum = std::array{track.px(), track.py(), track.pz()}; + // double pt = track.pt(); + // double phi = RecoDecay::phi(momentum); + // double eta = track.eta(); 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) + if (pdgCode == PDG_t::kElectron) + recoMC.SetM(o2::constants::physics::MassElectron); + else if (pdgCode == PDG_t::kMuonMinus) + recoMC.SetM(o2::constants::physics::MassMuon); + else if (pdgCode == PDG_t::kPiPlus) + recoMC.SetM(o2::constants::physics::MassPionCharged); + else if (pdgCode == PDG_t::kKPlus) + recoMC.SetM(o2::constants::physics::MassKaonCharged); + else if (pdgCode == PDG_t::kProton) + recoMC.SetM(o2::constants::physics::MassProton); + else continue; + + double pt = recoMC.Pt(); + double eta = recoMC.Eta(); + // if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton) + // continue; // if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance // continue; if (!mcParticle.isPhysicalPrimary()) From 1a853697307fa74ae90ffae0e26ed9ec1937827f Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Wed, 15 Apr 2026 15:33:43 +0800 Subject: [PATCH 5/8] preslice particles --- PWGUD/Tasks/flowMcUpc.cxx | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/PWGUD/Tasks/flowMcUpc.cxx b/PWGUD/Tasks/flowMcUpc.cxx index a927c075658..a068b297d9f 100644 --- a/PWGUD/Tasks/flowMcUpc.cxx +++ b/PWGUD/Tasks/flowMcUpc.cxx @@ -109,6 +109,8 @@ struct FlowMcUpc { return true; } + PresliceUnsorted partPerMcCollision = aod::udmcparticle::udMcCollisionId; + void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, McParts const& mcParts, aod::BCs const& bcs) { if (bcs.size() == 0) { @@ -122,7 +124,9 @@ struct FlowMcUpc { // event within range histos.fill(HIST("hImpactParameter"), imp); - for (auto const& mcParticle : mcParts) { + auto const& tempParts = mcParts.sliceBy(partPerMcCollision, static_cast(mcCollision.globalIndex())); + + for (auto const& mcParticle : tempParts) { auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; int pdgCode = std::abs(mcParticle.pdgCode()); @@ -147,6 +151,9 @@ struct FlowMcUpc { using MCRecoTracks = soa::Join; using MCRecoCollisions = soa::Join; + // PresliceUnsorted trackPerMcParticle = aod::udmctracklabel::udMcParticleId; + Preslice trackPerCollision = aod::udtrack::udCollisionId; // sorted preslice used because the pair track-collision is already sorted in processDataSG function + void processReco(MCRecoCollisions::iterator const& collision, MCRecoTracks const& tracks) { histos.fill(HIST("RecoProcessEventCounter"), 0.5); @@ -162,7 +169,9 @@ struct FlowMcUpc { float vtxz = collision.posZ(); - for (const auto& track : tracks) { + auto const& tempTracks = tracks.sliceBy(trackPerCollision, static_cast(collision.globalIndex())); + + for (const auto& track : tempTracks) { // focus on bulk: e, mu, pi, k, p auto momentum = std::array{track.px(), track.py(), track.pz()}; double pt = RecoDecay::pt(momentum); From 0dc9aab497545466996d5a8d0725ecf677d47150 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Thu, 16 Apr 2026 17:12:42 +0800 Subject: [PATCH 6/8] not joining particllabel --- PWGUD/Tasks/flowMcUpc.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGUD/Tasks/flowMcUpc.cxx b/PWGUD/Tasks/flowMcUpc.cxx index a068b297d9f..b35511573df 100644 --- a/PWGUD/Tasks/flowMcUpc.cxx +++ b/PWGUD/Tasks/flowMcUpc.cxx @@ -65,7 +65,7 @@ struct FlowMcUpc { double epsilon = 1e-6; - using McParts = soa::Join; + // using McParts = soa::Join; void init(InitContext&) { @@ -111,7 +111,7 @@ struct FlowMcUpc { PresliceUnsorted partPerMcCollision = aod::udmcparticle::udMcCollisionId; - void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, McParts const& mcParts, aod::BCs const& bcs) + void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, aod::UDMcParticles const& mcParts, aod::BCs const& bcs) { if (bcs.size() == 0) { return; From 751a111a2199fe0226d3c6fdc84cbcb4f85bd830 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Fri, 17 Apr 2026 13:32:29 +0800 Subject: [PATCH 7/8] iterate collsion manually --- PWGUD/Tasks/flowMcUpc.cxx | 124 ++++++++++++++++++++------------------ 1 file changed, 64 insertions(+), 60 deletions(-) diff --git a/PWGUD/Tasks/flowMcUpc.cxx b/PWGUD/Tasks/flowMcUpc.cxx index b35511573df..136626fb7a6 100644 --- a/PWGUD/Tasks/flowMcUpc.cxx +++ b/PWGUD/Tasks/flowMcUpc.cxx @@ -111,88 +111,92 @@ struct FlowMcUpc { PresliceUnsorted partPerMcCollision = aod::udmcparticle::udMcCollisionId; - void processMCTrue(aod::UDMcCollisions::iterator const& mcCollision, aod::UDMcParticles const& mcParts, aod::BCs const& bcs) + void processMCTrue(aod::UDMcCollisions const& mcCollisions, aod::UDMcParticles const& mcParts, aod::BCs const& bcs) { if (bcs.size() == 0) { return; } - histos.fill(HIST("mcEventCounter"), 0.5); - float imp = mcCollision.impactParameter(); - float vtxz = mcCollision.posZ(); + for (const auto& mcCollision : mcCollisions) { + 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); + if (imp >= minB && imp <= maxB) { + // event within range + histos.fill(HIST("hImpactParameter"), imp); - auto const& tempParts = mcParts.sliceBy(partPerMcCollision, static_cast(mcCollision.globalIndex())); + auto const& tempParts = mcParts.sliceBy(partPerMcCollision, static_cast(mcCollision.globalIndex())); - for (auto const& mcParticle : tempParts) { - auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; - int pdgCode = std::abs(mcParticle.pdgCode()); + for (auto const& mcParticle : tempParts) { + auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; + int pdgCode = std::abs(mcParticle.pdgCode()); - double pt = RecoDecay::pt(momentum); - double eta = RecoDecay::eta(momentum); + double pt = RecoDecay::pt(momentum); + double eta = RecoDecay::eta(momentum); - if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton) - continue; + if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton) + continue; - if (!mcParticle.isPhysicalPrimary()) - continue; - if (std::fabs(eta) > cfgCutEta) // main acceptance - continue; + if (!mcParticle.isPhysicalPrimary()) + continue; + if (std::fabs(eta) > cfgCutEta) // main acceptance + continue; - histos.fill(HIST("hPtMCGen"), pt); - histos.fill(HIST("hEtaPtVtxzMCGen"), eta, pt, vtxz); + 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; - using MCRecoCollisions = soa::Join; + using MCRecoCollisions = soa::Join; // PresliceUnsorted trackPerMcParticle = aod::udmctracklabel::udMcParticleId; Preslice trackPerCollision = aod::udtrack::udCollisionId; // sorted preslice used because the pair track-collision is already sorted in processDataSG function - void processReco(MCRecoCollisions::iterator const& collision, MCRecoTracks const& tracks) + void processReco(MCRecoCollisions const& collisions, 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(); - - auto const& tempTracks = tracks.sliceBy(trackPerCollision, static_cast(collision.globalIndex())); - - for (const auto& track : tempTracks) { - // focus on bulk: e, mu, pi, k, p - auto momentum = std::array{track.px(), track.py(), track.pz()}; - double pt = RecoDecay::pt(momentum); - double eta = RecoDecay::eta(momentum); - // double phi = RecoDecay::phi(momentum); - if (!trackSelected(track) || (!track.has_udMcParticle())) - continue; - auto mcParticle = track.udMcParticle(); - int pdgCode = std::abs(mcParticle.pdgCode()); - - // double pt = recoMC.Pt(); - // double eta = recoMC.Eta(); - if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton) - continue; - if (std::fabs(eta) > cfgCutEta) // main acceptance - continue; - if (!mcParticle.isPhysicalPrimary()) - continue; - - histos.fill(HIST("hPtReco"), pt); - histos.fill(HIST("hEtaPtVtxzMCReco"), eta, pt, vtxz); + for (const auto& collision : collisions) { + 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(); + + auto const& tempTracks = tracks.sliceBy(trackPerCollision, static_cast(collision.globalIndex())); + + for (const auto& track : tempTracks) { + // focus on bulk: e, mu, pi, k, p + auto momentum = std::array{track.px(), track.py(), track.pz()}; + double pt = RecoDecay::pt(momentum); + double eta = RecoDecay::eta(momentum); + // double phi = RecoDecay::phi(momentum); + if (!trackSelected(track) || (!track.has_udMcParticle())) + continue; + auto mcParticle = track.udMcParticle(); + int pdgCode = std::abs(mcParticle.pdgCode()); + + // double pt = recoMC.Pt(); + // double eta = recoMC.Eta(); + if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton) + continue; + if (std::fabs(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); From 6b7bd53c65fe1575448f8920b6f4e9fbd45e13c1 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Mon, 20 Apr 2026 16:09:30 +0800 Subject: [PATCH 8/8] add eff & independent axis --- PWGUD/Tasks/CMakeLists.txt | 2 +- PWGUD/Tasks/flowCorrelationsUpc.cxx | 407 ++++++++++++++++++++++++---- 2 files changed, 354 insertions(+), 55 deletions(-) diff --git a/PWGUD/Tasks/CMakeLists.txt b/PWGUD/Tasks/CMakeLists.txt index ab54549251c..bf1c0e7be72 100644 --- a/PWGUD/Tasks/CMakeLists.txt +++ b/PWGUD/Tasks/CMakeLists.txt @@ -251,7 +251,7 @@ o2physics_add_dpl_workflow(flow-cumulants-upc o2physics_add_dpl_workflow(flow-correlations-upc SOURCES flowCorrelationsUpc.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGCFCore + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::GFWCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(flow-mc-upc diff --git a/PWGUD/Tasks/flowCorrelationsUpc.cxx b/PWGUD/Tasks/flowCorrelationsUpc.cxx index 6bd5736ee48..8ed0140fa95 100644 --- a/PWGUD/Tasks/flowCorrelationsUpc.cxx +++ b/PWGUD/Tasks/flowCorrelationsUpc.cxx @@ -15,11 +15,17 @@ /// copied from Thor Jensen (thor.kjaersgaard.jensen@cern.ch) and Debojit Sarkar (debojit.sarkar@cern.ch) #include "PWGCF/Core/CorrelationContainer.h" +#include "PWGCF/GenericFramework/Core/GFW.h" +#include "PWGCF/GenericFramework/Core/GFWWeights.h" #include "PWGUD/Core/SGSelector.h" #include "PWGUD/DataModel/UDTables.h" +#include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include #include #include #include @@ -31,14 +37,33 @@ #include #include #include +#include #include +#include +#include +#include +#include +#include +#include +#include #include #include +#include + +#include + +#include +#include #include #include #include +#include +#include +#include +#include +#include #include namespace o2::aod @@ -129,8 +154,6 @@ struct FlowCorrelationsUpc { O2_DEFINE_CONFIGURABLE(cfgRadiusLow, float, 0.8, "Low radius for merging cut") O2_DEFINE_CONFIGURABLE(cfgRadiusHigh, float, 2.5, "High radius for merging cut") O2_DEFINE_CONFIGURABLE(cfgIsGoodItsLayers, bool, false, "whether choose itslayers") - O2_DEFINE_CONFIGURABLE(cfgGapSide, int, 1, "0: gapside A;1:C") - O2_DEFINE_CONFIGURABLE(cfgGapSideMerge, bool, false, "whether merge A and C side together") O2_DEFINE_CONFIGURABLE(cfgDcaxy, bool, true, "choose dcaxy") O2_DEFINE_CONFIGURABLE(cfgDcaz, bool, false, "choose dcaz") O2_DEFINE_CONFIGURABLE(cfgDcazCut, float, 10.0, "dcaz cut") @@ -139,6 +162,19 @@ struct FlowCorrelationsUpc { O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, true, "Occupancy cut") O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 1000, "High cut on TPC occupancy") O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") + O2_DEFINE_CONFIGURABLE(cfgCutTPCCrossedRows, float, 70.0f, "minimum number of crossed TPC Rows") + O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum number of found TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum number of ITS clusters") + O2_DEFINE_CONFIGURABLE(cfgCutTPCChi2NCl, int, 4, "max chi2 per TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgGlobalTrack, bool, true, "require TPC+ITS track") + O2_DEFINE_CONFIGURABLE(cfgUseNchCorrected, bool, true, "use corrected Nch for X axis") + O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights") + O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeightsRefPt, bool, false, "NUA weights are filled in ref pt bins") + O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeightsRunbyRun, bool, false, "NUA weights are filled run-by-run") + O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") + O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object") + O2_DEFINE_CONFIGURABLE(cfgUseEventWeights, bool, false, "Use event weights for mixed event") + O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object") ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"}; ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"}; @@ -163,11 +199,31 @@ struct FlowCorrelationsUpc { Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; + ConfigurableAxis axisIndependent{"axisIndependent", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60}, "X axis for histograms"}; + ConfigurableAxis axisNch{"axisNch", {300, 0, 300}, "N_{ch}"}; + + // Corrections + TH1D* mEfficiency = nullptr; + GFWWeights* mAcceptance = nullptr; + bool correctionsLoaded = false; // make the filters and cuts. Filter trackFilter = (aod::udtrack::isPVContributor == true); Filter collisionFilter = ((aod::udcollision::gapSide == (uint8_t)1 || aod::udcollision::gapSide == (uint8_t)0) && (cfgIfVertex == false || aod::collision::posZ < cfgZVtxCut) && (aod::udcollision::occupancyInTime > 0 && aod::udcollision::occupancyInTime < cfgCutOccupancyHigh) && (aod::flowcorrupc::truegapside == 1 || aod::flowcorrupc::truegapside == 0)); + // Connect to ccdb + Service ccdb; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + OutputObj fWeights{GFWWeights("weights")}; + OutputObj fWeightsMc{GFWWeights("weightsMC")}; + + TAxis* fPtAxis = nullptr; + int lastRunNumber = -1; + std::vector runNumbers; + std::map> th3sPerRun; // map of TH3 histograms for all runs + std::vector efficiencyCache; + using UdTracks = soa::Filtered>; using UdTracksFull = soa::Filtered>; @@ -181,12 +237,12 @@ struct FlowCorrelationsUpc { void init(InitContext&) { - LOG(info) << "cfgGapSide = " << cfgGapSide; - LOG(info) << "cfgGapSide value type: " << typeid(cfgGapSide).name(); LOGF(info, "Starting init"); + ccdb->setURL(ccdbUrl.value); + ccdb->setCaching(true); + auto now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); // Make histograms to check the distributions after cuts - registry.add("deltaEta_deltaPhi_same", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}}); registry.add("Phi", "Phi", {HistType::kTH1D, {axisPhi}}); registry.add("Eta", "Eta", {HistType::kTH1D, {axisEta}}); registry.add("pT", "pT", {HistType::kTH1D, {axisPtTrigger}}); @@ -195,17 +251,38 @@ struct FlowCorrelationsUpc { registry.add("Nch_vs_zVtx", "Nch vs zVtx", {HistType::kTH2D, {axisVertex, axisMultiplicity}}); registry.add("Nch_same", "Nch same event", {HistType::kTH1D, {axisMultiplicity}}); registry.add("Nch_mixed", "Nch mixed event", {HistType::kTH1D, {axisMultiplicity}}); - - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - - registry.add("eventcount_same", "bin", {HistType::kTH1F, {{10, 0, 10, "bin"}}}); // histogram to see how many events are in the same and mixed event - registry.add("eventcount_mixed", "bin", {HistType::kTH1F, {{10, 0, 10, "bin"}}}); // histogram to see how many events are in the same and mixed event - - registry.add("trackcount_same", "bin", {HistType::kTH1F, {{10, 0, 10, "bin"}}}); // histogram to see how many tracks are in the same and mixed event - registry.add("trackcount_mixed", "bin", {HistType::kTH1F, {{10, 0, 10, "bin"}}}); // histogram to see how many tracks are in the same and mixed event + registry.add("EtaCorrected", "Eta corrected", {HistType::kTH1D, {axisEta}}); + registry.add("pTCorrected", "pT corrected", {HistType::kTH1D, {axisPtTrigger}}); + + registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisIndependent, axisPtTrigger}}}); + + registry.add("eventcont", "bin", {HistType::kTH1F, {{10, 0, 10, "bin"}}}); // histogram to see how many events are in the same and mixed event + registry.add("deltaEta_deltaPhi_same", "deltaeta-deltaphi", {HistType::kTH2D, {axisDeltaEta, axisDeltaPhi}}); // histogram to check the delta eta and delta phi distribution + registry.add("deltaEta_deltaPhi_mixed", "deltaeta-deltaphi", {HistType::kTH2D, {axisDeltaEta, axisDeltaPhi}}); // histogram to check the delta eta and delta phi distribution + registry.add("Nch_raw_vs_independent", "Raw vs Independent", {HistType::kTH2D, {axisMultiplicity, axisIndependent}}); + + // if (doprocessSim) { + // registry.add("eventCounterMC", "Number of MC Events;; Count", {HistType::kTH1D, {{5, 0, 5}}}); + // registry.add("hVtxZMC", "Vexter Z distribution (MC)", {HistType::kTH1D, {axisVertex}}); + // registry.add("hMultMC", "Multiplicity distribution (MC)", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); + // registry.add("numberOfTracksMC", "Number of MC tracks;; Count", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); + // } + + o2::framework::AxisSpec axis = axisPtTrigger; + int nPtBins = axis.binEdges.size() - 1; + double* ptBins = &(axis.binEdges)[0]; + fPtAxis = new TAxis(nPtBins, ptBins); + + if (cfgOutputNUAWeights) { + fWeights->setPtBins(nPtBins, ptBins); + fWeights->init(true, false); + fWeightsMc->setPtBins(nPtBins, ptBins); + fWeightsMc->init(true, false); + } std::vector corrAxis = {{axisSample, "Sample"}, {axisVertex, "z-vtx (cm)"}, + {axisIndependent, "Independent (N_{ch} corrected)"}, {axisPtTrigger, "p_{T} (GeV/c)"}, {axisPtAssoc, "p_{T} (GeV/c)"}, {axisDeltaPhi, "#Delta#varphi (rad)"}, @@ -266,46 +343,165 @@ struct FlowCorrelationsUpc { return false; } // registry.fill(HIST("hTrackCount"), 3.5); - if (track.itsClusterSizes() <= cfgItsClusterSize) { + if (track.itsNCls() <= cfgCutITSclu) { return false; } // registry.fill(HIST("hTrackCount"), 4.5); if (track.tpcChi2NCl() >= cfgMaxTPCChi2NCl) { return false; } + if (track.tpcNClsCrossedRows() < cfgCutTPCCrossedRows) { + return false; + } + auto tpcClu = track.tpcNClsFindable() - track.tpcNClsFindableMinusFound(); + if (tpcClu < cfgCutTPCclu) { + return false; + } // registry.fill(HIST("hTrackCount"), 5.5); return true; } + void createOutputObjectsForRun(int runNumber) + { + const AxisSpec axisPhi{60, 0.0, constants::math::TwoPI, "#varphi"}; + std::shared_ptr histPhiEtaVtxz = registry.add(Form("%d/hPhiEtaVtxz", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.8, 1.8}, {40, -10, 10}}}); + th3sPerRun.insert(std::make_pair(runNumber, histPhiEtaVtxz)); + } + + void loadCorrections(uint64_t timestamp) + { + if (correctionsLoaded) { + return; + } + if (cfgAcceptance.value.empty() == false) { + mAcceptance = ccdb->getForTimeStamp(cfgAcceptance, timestamp); + if (mAcceptance) { + LOGF(info, "Loaded acceptance weights from %s (%p)", cfgAcceptance.value.c_str(), (void*)mAcceptance); + } else { + LOGF(warning, "Could not load acceptance weights from %s (%p)", cfgAcceptance.value.c_str(), (void*)mAcceptance); + } + } + if (cfgEfficiency.value.empty() == false) { + mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); + if (mEfficiency == nullptr) { + LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiency.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)mEfficiency); + } + correctionsLoaded = true; + } + + bool getEfficiencyCorrection(float& weight_nue, float eta, float pt, float posZ) + { + float eff = 1.; + if (mEfficiency) { + int etaBin = mEfficiency->GetXaxis()->FindBin(eta); + int ptBin = mEfficiency->GetYaxis()->FindBin(pt); + int zBin = mEfficiency->GetZaxis()->FindBin(posZ); + eff = mEfficiency->GetBinContent(etaBin, ptBin, zBin); + } else { + eff = 1.0; + } + if (eff <= 0) + return false; + weight_nue = 1. / eff; + return true; + } + + bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float phi, float eta, float pt, float vtxz) + { + float eff = 1.; + if (mEfficiency) + eff = mEfficiency->GetBinContent(mEfficiency->FindBin(pt)); + else + eff = 1.0; + if (eff == 0) + return false; + weight_nue = 1. / eff; + + if (mAcceptance) + weight_nua = mAcceptance->getNUA(phi, eta, vtxz); + else + weight_nua = 1; + return true; + } + // fill multiple histograms template - void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms. + void fillYield(TCollision collision, TTracks tracks, int runNumber, float vtxz) // function to fill the yield and etaphi histograms. { registry.fill(HIST("Nch"), tracks.size()); registry.fill(HIST("zVtx"), collision.posZ()); for (auto const& track1 : tracks) { - auto momentum1 = std::array{track1.px(), track1.py(), track1.pz()}; - registry.fill(HIST("Phi"), RecoDecay::phi(momentum1)); - registry.fill(HIST("Eta"), RecoDecay::eta(momentum1)); - registry.fill(HIST("pT"), track1.pt()); + if (!trackSelected(track1)) + continue; + auto momentum = std::array{track1.px(), track1.py(), track1.pz()}; + double pt = RecoDecay::pt(momentum); + double phi = RecoDecay::phi(momentum); + double eta = RecoDecay::eta(momentum); + float weff = 1.; + if (!getEfficiencyCorrection(weff, eta, pt, vtxz)) + continue; + + registry.fill(HIST("Phi"), phi); + registry.fill(HIST("Eta"), eta); + registry.fill(HIST("pT"), pt); + registry.fill(HIST("EtaCorrected"), eta, weff); + registry.fill(HIST("pTCorrected"), pt, weff); + + if (cfgOutputNUAWeights) { + if (cfgOutputNUAWeightsRefPt) { + if (pt > cfgPtCutMin && pt < cfgPtCutMax) { + fWeights->fill(phi, eta, vtxz, pt, 0., 0); + if (cfgOutputNUAWeightsRunbyRun) + th3sPerRun[runNumber]->Fill(phi, eta, vtxz); + } + } else { + fWeights->fill(phi, eta, vtxz, pt, 0., 0); + if (cfgOutputNUAWeightsRunbyRun) + th3sPerRun[runNumber]->Fill(phi, eta, vtxz); + } + } } } template - void fillCorrelations(TTracks tracks1, TTracks tracks2, float posZ, int system, int runnum) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms + void fillCorrelations(TTracks tracks1, TTracks tracks2, float posZ, int system, int runnum, float vtxz, float eventWeight, double independent) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms { + if (mEfficiency) { + efficiencyCache.clear(); + efficiencyCache.reserve(static_cast(tracks2.size())); + for (const auto& track2 : tracks2) { + auto momentum = std::array{track2.px(), track2.py(), track2.pz()}; + double pt = RecoDecay::pt(momentum); + double eta = RecoDecay::eta(momentum); + float weff = 1.; + getEfficiencyCorrection(weff, eta, pt, vtxz); + efficiencyCache.push_back(weff); + } + } + int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); // loop over all tracks for (auto const& track1 : tracks1) { - if (!trackSelected(track1)) continue; + auto momentum = std::array{track1.px(), track1.py(), track1.pz()}; + double pt1 = RecoDecay::pt(momentum); + double phi1 = RecoDecay::phi(momentum); + double eta1 = RecoDecay::eta(momentum); + + // 计算track1的权重 + float weff1 = 1., wacc1 = 1.; + if (!setCurrentParticleWeights(weff1, wacc1, phi1, eta1, pt1, vtxz)) + continue; + if (system == SameEvent) { - registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt()); + registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, independent, pt1, eventWeight * weff1 * wacc1); } for (auto const& track2 : tracks2) { @@ -313,24 +509,42 @@ struct FlowCorrelationsUpc { continue; if (track1.globalIndex() == track2.globalIndex()) - continue; // For pt-differential correlations, skip if the trigger and associate are the same track - if (system == SameEvent && track1.pt() <= track2.pt()) - continue; // Without pt-differential correlations, skip if the trigger pt is less than the associate pt - - auto momentum1 = std::array{track1.px(), track1.py(), track1.pz()}; - auto momentum2 = std::array{track2.px(), track2.py(), track2.pz()}; - double phi1 = RecoDecay::phi(momentum1); - double phi2 = RecoDecay::phi(momentum2); + continue; + if (system == SameEvent && cfgUsePtOrder && pt1 <= track2.pt()) + continue; + if (system == MixedEvent && cfgUsePtOrderInMixEvent && pt1 <= track2.pt()) + continue; + + auto momentum = std::array{track2.px(), track2.py(), track2.pz()}; + double pt2 = RecoDecay::pt(momentum); + double phi2 = RecoDecay::phi(momentum); + double eta2 = RecoDecay::eta(momentum); + + // 计算track2的权重 + float weff2 = 1., wacc2 = 1.; + if (mEfficiency) { + weff2 = efficiencyCache[track2.filteredIndex()]; + } else { + getEfficiencyCorrection(weff2, eta2, pt2, vtxz); + } + + if (mAcceptance) { + wacc2 = mAcceptance->getNUA(phi2, eta2, vtxz); + } else { + wacc2 = 1; + } + float deltaPhi = RecoDecay::constrainAngle(phi1 - phi2, -PIHalf); - float deltaEta = RecoDecay::eta(momentum1) - RecoDecay::eta(momentum2); + float deltaEta = eta1 - eta2; - if (std::abs(deltaEta) < cfgCutMerging) { + // 计算组合权重 + float weight = eventWeight * weff1 * weff2 * wacc1 * wacc2; + // Merging cut + if (std::abs(deltaEta) < cfgCutMerging) { double dPhiStarHigh = getDPhiStar(track1, track2, cfgRadiusHigh, runnum, phi1, phi2); double dPhiStarLow = getDPhiStar(track1, track2, cfgRadiusLow, runnum, phi1, phi2); - const double kLimit = 3.0 * cfgCutMerging; - bool bIsBelow = false; if (std::abs(dPhiStarLow) < kLimit || std::abs(dPhiStarHigh) < kLimit || dPhiStarLow * dPhiStarHigh < 0) { @@ -346,13 +560,13 @@ struct FlowCorrelationsUpc { } } - // fill the right sparse and histograms + // fill the right sparse and histograms with weights if (system == SameEvent) { - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta); + same->getPairHist()->Fill(step, fSampleIndex, posZ, independent, pt1, pt2, deltaPhi, deltaEta, weight); + registry.fill(HIST("deltaEta_deltaPhi_same"), deltaPhi, deltaEta, weight); } else if (system == MixedEvent) { - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta); + mixed->getPairHist()->Fill(step, fSampleIndex, posZ, independent, pt1, pt2, deltaPhi, deltaEta, weight); + registry.fill(HIST("deltaEta_deltaPhi_mixed"), deltaPhi, deltaEta, weight); } } } @@ -364,13 +578,61 @@ struct FlowCorrelationsUpc { if (tracks.size() < cfgMinMult || tracks.size() > cfgMaxMult) { return; } - registry.fill(HIST("eventcount_same"), 3.5); - int runIndex = collision.runNumber(); - // registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - registry.fill(HIST("Nch_vs_zVtx"), collision.posZ(), tracks.size()); - fillYield(collision, tracks); - fillCorrelations(tracks, tracks, collision.posZ(), SameEvent, runIndex); // fill the SE histogram and Sparse + float vtxz = collision.posZ(); + auto currentRunNumber = collision.runNumber(); + auto runDuration = ccdb->getRunDuration(currentRunNumber); + + loadCorrections(runDuration.first); + + if (cfgOutputNUAWeightsRunbyRun && currentRunNumber != lastRunNumber) { + lastRunNumber = currentRunNumber; + if (std::find(runNumbers.begin(), runNumbers.end(), currentRunNumber) == runNumbers.end()) { + createOutputObjectsForRun(currentRunNumber); + runNumbers.push_back(currentRunNumber); + } + + if (th3sPerRun.find(currentRunNumber) == th3sPerRun.end()) { + LOGF(fatal, "RunNumber %d not found in th3sPerRun", currentRunNumber); + return; + } + } + + registry.fill(HIST("eventcont"), 3.5); + + //-----------independent--------------- + double nTracksRaw = 0.; + double nTracksCorrected = 0.; + + for (const auto& track : tracks) { + if (!trackSelected(track)) + continue; + + auto momentum = std::array{track.px(), track.py(), track.pz()}; + double pt = RecoDecay::pt(momentum); + double eta = RecoDecay::eta(momentum); + + nTracksRaw += 1.; + + if (cfgUseNchCorrected) { + float weff = 1.; + if (getEfficiencyCorrection(weff, eta, pt, vtxz)) { + nTracksCorrected += weff; + } + } + } + registry.fill(HIST("Nch_raw_vs_independent"), nTracksRaw, nTracksCorrected); + + double independent = nTracksRaw; + if (cfgUseNchCorrected) { + independent = nTracksCorrected; + } + + fillYield(collision, tracks, currentRunNumber, vtxz); + + fillCorrelations( + tracks, tracks, collision.posZ(), SameEvent, + currentRunNumber, vtxz, 1.0f, independent); } PROCESS_SWITCH(FlowCorrelationsUpc, processSame, "Process same event", true); @@ -384,24 +646,61 @@ struct FlowCorrelationsUpc { { auto getTracksSize = [&tracks, this](UDCollisionsFull::iterator const& collision) { auto associatedTracks = tracks.sliceByCached(o2::aod::udtrack::udCollisionId, collision.udCollisionId(), this->cache); - auto mult = associatedTracks.size(); return mult; }; using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; - // MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default) auto tracksTuple = std::make_tuple(tracks); - SameKindPair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + SameKindPair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; - for (auto const& [collision1, tracks1, collision2, tracks2] : pairs) { - // registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin - if (tracks1.size() < cfgMinMult || tracks1.size() > cfgMaxMult || tracks2.size() < cfgMinMult || tracks2.size() > cfgMaxMult) { + for (auto it = pairs.begin(); it != pairs.end(); it++) { + auto& [collision1, tracks1, collision2, tracks2] = *it; + if (tracks1.size() < cfgMinMult || tracks1.size() > cfgMaxMult || + tracks2.size() < cfgMinMult || tracks2.size() > cfgMaxMult) { continue; } - registry.fill(HIST("eventcount_same"), 4.5); - fillCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, collision1.runNumber()); // fill the ME histogram and Sparse + + auto runDuration1 = ccdb->getRunDuration(collision1.runNumber()); + loadCorrections(runDuration1.first); + + registry.fill(HIST("eventcont"), 4.5); + + double nTracksRaw = 0.; + double nTracksCorrected = 0.; + + for (const auto& track : tracks1) { + if (!trackSelected(track)) + continue; + + auto momentum = std::array{track.px(), track.py(), track.pz()}; + double pt = RecoDecay::pt(momentum); + double eta = RecoDecay::eta(momentum); + + nTracksRaw += 1.; + + if (cfgUseNchCorrected) { + float weff = 1.; + if (getEfficiencyCorrection(weff, eta, pt, collision1.posZ())) { + nTracksCorrected += weff; + } + } + } + + double independent = nTracksRaw; + if (cfgUseNchCorrected) { + independent = nTracksCorrected; + } + + float eventWeight = 1.0f; + if (cfgUseEventWeights) { + eventWeight = 1.0f / it.currentWindowNeighbours(); + } + + fillCorrelations( + tracks1, tracks2, collision1.posZ(), MixedEvent, + collision1.runNumber(), collision1.posZ(), eventWeight, independent); } } PROCESS_SWITCH(FlowCorrelationsUpc, processMixed, "Process mixed events", true);