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/flowCorrelationsUpc.cxx b/PWGUD/Tasks/flowCorrelationsUpc.cxx index 061eda25566..6bd5736ee48 100644 --- a/PWGUD/Tasks/flowCorrelationsUpc.cxx +++ b/PWGUD/Tasks/flowCorrelationsUpc.cxx @@ -377,12 +377,21 @@ struct FlowCorrelationsUpc { // event mixing SliceCache cache; - using MixedBinning = ColumnBinningPolicy; + // using MixedBinning = ColumnBinningPolicy; // the process for filling the mixed events void processMixed(UDCollisionsFull const& collisions, UdTracksFull const& tracks) { - MixedBinning binningOnVtxAndMult{{vtxMix, multMix}, true}; // true is for 'ignore overflows' (true by default) + 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/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 9624c1a5c86..7d831eb2bf8 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -404,33 +404,34 @@ struct FlowCumulantsUpc { delete oba; // eta region - fGFW->AddRegion("full", -0.8, 0.8, 1, 1); - fGFW->AddRegion("refN00", -0.8, 0., 1, 1); // gap0 negative region - fGFW->AddRegion("refP00", 0., 0.8, 1, 1); // gap0 positve region - fGFW->AddRegion("refN02", -0.8, -0.1, 1, 1); // gap2 negative region - fGFW->AddRegion("refP02", 0.1, 0.8, 1, 1); // gap2 positve region - fGFW->AddRegion("refN04", -0.8, -0.2, 1, 1); // gap4 negative region - fGFW->AddRegion("refP04", 0.2, 0.8, 1, 1); // gap4 positve region - fGFW->AddRegion("refN06", -0.8, -0.3, 1, 1); // gap6 negative region - fGFW->AddRegion("refP06", 0.3, 0.8, 1, 1); // gap6 positve region - fGFW->AddRegion("refN08", -0.8, -0.4, 1, 1); - fGFW->AddRegion("refP08", 0.4, 0.8, 1, 1); - fGFW->AddRegion("refN10", -0.8, -0.5, 1, 1); - fGFW->AddRegion("refP10", 0.5, 0.8, 1, 1); - fGFW->AddRegion("refN12", -0.8, -0.6, 1, 1); - fGFW->AddRegion("refP12", 0.6, 0.8, 1, 1); - fGFW->AddRegion("refN14", -0.8, -0.7, 1, 1); - fGFW->AddRegion("refP14", 0.7, 0.8, 1, 1); - fGFW->AddRegion("refN", -0.8, -0.4, 1, 1); - fGFW->AddRegion("refP", 0.4, 0.8, 1, 1); + fGFW->AddRegion("full", -0.9, 0.9, 1, 1); + fGFW->AddRegion("refN00", -0.9, 0., 1, 1); // gap0 negative region + fGFW->AddRegion("refP00", 0., 0.9, 1, 1); // gap0 positve region + fGFW->AddRegion("refN02", -0.9, -0.1, 1, 1); // gap2 negative region + fGFW->AddRegion("refP02", 0.1, 0.9, 1, 1); // gap2 positve region + fGFW->AddRegion("refN04", -0.9, -0.2, 1, 1); // gap4 negative region + fGFW->AddRegion("refP04", 0.2, 0.9, 1, 1); // gap4 positve region + fGFW->AddRegion("refN06", -0.9, -0.3, 1, 1); // gap6 negative region + fGFW->AddRegion("refP06", 0.3, 0.9, 1, 1); // gap6 positve region + fGFW->AddRegion("refN08", -0.9, -0.4, 1, 1); + fGFW->AddRegion("refP08", 0.4, 0.9, 1, 1); + fGFW->AddRegion("refN10", -0.9, -0.5, 1, 1); + fGFW->AddRegion("refP10", 0.5, 0.9, 1, 1); + fGFW->AddRegion("refN12", -0.9, -0.6, 1, 1); + fGFW->AddRegion("refP12", 0.6, 0.9, 1, 1); + fGFW->AddRegion("refN14", -0.9, -0.7, 1, 1); + fGFW->AddRegion("refP14", 0.7, 0.9, 1, 1); + fGFW->AddRegion("refN", -0.9, -0.4, 1, 1); + fGFW->AddRegion("refP", 0.4, 0.9, 1, 1); fGFW->AddRegion("refM", -0.4, 0.4, 1, 1); - fGFW->AddRegion("poiN", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 2); - fGFW->AddRegion("poiN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 2); - fGFW->AddRegion("poifull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 2); - fGFW->AddRegion("olN", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 4); - fGFW->AddRegion("olN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 4); - fGFW->AddRegion("olfull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 4); - + fGFW->AddRegion("poiN", -0.9, -0.4, 1 + fPtAxis->GetNbins(), 2); + fGFW->AddRegion("poiN10", -0.9, -0.5, 1 + fPtAxis->GetNbins(), 2); + fGFW->AddRegion("poifull", -0.9, 0.9, 1 + fPtAxis->GetNbins(), 2); + fGFW->AddRegion("olN", -0.9, -0.4, 1 + fPtAxis->GetNbins(), 4); + fGFW->AddRegion("olN10", -0.9, -0.5, 1 + fPtAxis->GetNbins(), 4); + fGFW->AddRegion("olfull", -0.9, 0.9, 1 + fPtAxis->GetNbins(), 4); + + // eta region for MC, can be different from data to study the effect of acceptance fGFWMC->AddRegion("full", -0.8, 0.8, 1, 1); fGFWMC->AddRegion("refN00", -0.8, 0., 1, 1); // gap0 negative region fGFWMC->AddRegion("refP00", 0., 0.8, 1, 1); // gap0 positve region diff --git a/PWGUD/Tasks/flowMcUpc.cxx b/PWGUD/Tasks/flowMcUpc.cxx new file mode 100644 index 00000000000..a927c075658 --- /dev/null +++ b/PWGUD/Tasks/flowMcUpc.cxx @@ -0,0 +1,196 @@ +// 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 +#include +#include +#include +#include +#include +#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 McParts = 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 const& track) + { + // auto momentum = std::array{track.px(), track.py(), track.pz()}; + auto pt = track.pt(); + 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, McParts 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(); + + if (imp >= minB && imp <= maxB) { + // event within range + histos.fill(HIST("hImpactParameter"), imp); + + for (auto const& mcParticle : mcParts) { + 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); + + 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; + + 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) { + // 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); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}