From b6008adb3d637b33948699e8181a53bdbcefdc4e Mon Sep 17 00:00:00 2001 From: mkruger Date: Fri, 24 Oct 2025 15:22:36 +0200 Subject: [PATCH 1/2] use PCC in charged particle analysis --- PWGLF/Tasks/Nuspex/chargedParticles.cxx | 104 ++++++++++++++++++------ 1 file changed, 79 insertions(+), 25 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/chargedParticles.cxx b/PWGLF/Tasks/Nuspex/chargedParticles.cxx index 42bbb437a89..9408c94c7ce 100644 --- a/PWGLF/Tasks/Nuspex/chargedParticles.cxx +++ b/PWGLF/Tasks/Nuspex/chargedParticles.cxx @@ -18,6 +18,7 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "PWGLF/DataModel/particleCompositionCorrectionTable.h" #include #include @@ -27,6 +28,7 @@ #include #include +#include using namespace o2; using namespace o2::framework; @@ -38,6 +40,19 @@ using aod::track::TrackSelectionFlags; //-------------------------------------------------------------------------------------------------- //-------------------------------------------------------------------------------------------------- struct ChargedParticles { + std::random_device rnd; + std::mt19937 gen{rnd()}; + std::uniform_real_distribution dist{0.f, 1.f}; + uint32_t getNRepetitons(float scalingFactor) + { + uint32_t nRepetitions = static_cast(scalingFactor); + float rest = scalingFactor - nRepetitions; + if (rest) { + nRepetitions += (dist(gen) <= rest) ? 1u : 0u; + // LOGP(info, "scalingFactor: {} -> {}", scalingFactor, nRepetitions); + } + return nRepetitions; + }; HistogramRegistry histos; Service pdg; @@ -62,9 +77,11 @@ struct ChargedParticles { kSystUpDCAxy, kSystDownDCAz, kSystUpDCAz, - kSystITSHits, + kSystITSHits, // only relevant for converted data kSystDownTPCCrossedRows, - kSystUpTPCCrossedRows + kSystUpTPCCrossedRows, + kSystDownPCC = 120, + kSystUpPCC }; Configurable systMode{"systMode", kSystNominal, "variation for systematic uncertainties"}; uint16_t trackSelMask{TrackSelectionFlags::kGlobalTrackWoPtEta}; // track selection bitmask (without cut that is being varied) @@ -80,6 +97,7 @@ struct ChargedParticles { bool isAcceptedEventMC{false}; bool isINELGT0EventMC{false}; bool isChargedPrimary{false}; + uint32_t nRepetitions{1u}; }; VarContainer vars; static constexpr float kMaxVtxZ = 10.f; @@ -92,7 +110,7 @@ struct ChargedParticles { template bool initTrack(const T& track); - template + template void initEvent(const C& collision, const T& tracks); template @@ -112,9 +130,9 @@ struct ChargedParticles { using CollisionTableMCTrue = aod::McCollisions; using CollisionTableMC = soa::SmallGroups>; using TrackTableMC = soa::Join; - using ParticleTableMC = aod::McParticles; + using ParticleTableMC = soa::Join; Preslice perCollision = aod::track::collisionId; - void processMC(CollisionTableMCTrue::iterator const& mcCollision, CollisionTableMC const& collisions, TrackTableMC const& tracks, ParticleTableMC const& particles); + void processMC(CollisionTableMCTrue::iterator const& mcCollision, TrackTableMC const& tracks, CollisionTableMC const& collisions, ParticleTableMC const& particles); PROCESS_SWITCH(ChargedParticles, processMC, "process mc", true); }; @@ -241,7 +259,7 @@ void ChargedParticles::init(InitContext const&) //************************************************************************************************** void ChargedParticles::processData(CollisionTableData::iterator const& collision, TrackTableData const& tracks) { - initEvent(collision, tracks); + initEvent(collision, tracks); processMeas(tracks); } @@ -250,7 +268,7 @@ void ChargedParticles::processData(CollisionTableData::iterator const& collision * Entrypoint to processes mc. */ //************************************************************************************************** -void ChargedParticles::processMC(CollisionTableMCTrue::iterator const& mcCollision, CollisionTableMC const& collisions, TrackTableMC const& tracks, ParticleTableMC const& particles) +void ChargedParticles::processMC(CollisionTableMCTrue::iterator const& mcCollision, TrackTableMC const& tracks, CollisionTableMC const& collisions, ParticleTableMC const& particles) { histos.fill(HIST("collision_ambiguity"), collisions.size()); @@ -270,7 +288,7 @@ void ChargedParticles::processMC(CollisionTableMCTrue::iterator const& mcCollisi } else { for (const auto& collision : collisions) { auto curTracks = tracks.sliceBy(perCollision, collision.globalIndex()); - initEvent(collision, curTracks); + initEvent(collision, curTracks); processMeas(curTracks); break; // for now look only at first collision... } @@ -302,6 +320,15 @@ bool ChargedParticles::initParticle(const P& particle) if (particle.pt() <= ptMinCut || particle.pt() >= ptMaxCut) { return false; } + + float pccWeight = particle.pccWeight(); + if (systMode == kSystDownPCC) { + pccWeight = particle.pccWeightSysDown(); + } else if (systMode == kSystUpPCC) { + pccWeight = particle.pccWeightSysUp(); + } + vars.nRepetitions = getNRepetitons(pccWeight); + // FIXME: in case of nRepetitions = 0 INELGT0 can be wrong return true; } @@ -334,20 +361,43 @@ bool ChargedParticles::initTrack(const T& track) * Check if event is good. */ //************************************************************************************************** -template +template void ChargedParticles::initEvent(const C& collision, const T& tracks) { vars.multMeas = 0; for (const auto& track : tracks) { if (initTrack(track)) { - ++vars.multMeas; + if constexpr (IS_MC) { + if (!track.has_mcParticle()) { + continue; + } + const auto& particle = track.template mcParticle_as(); + if (!initParticle(particle)) { + continue; + } + vars.multMeas += vars.nRepetitions; + } else { + ++vars.multMeas; + } } } vars.isAcceptedEvent = false; if (std::abs(collision.posZ()) < kMaxVtxZ) { - if (isRun3 ? collision.sel8() : collision.sel7()) { - if ((isRun3 || doprocessMC) ? true : collision.alias_bit(kINT7)) { + if (isRun3) { + if (collision.sel8() && + collision.selection_bit(aod::evsel::kNoSameBunchPileup) && + collision.selection_bit(aod::evsel::kIsVertexITSTPC) && + collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + // !collision.selection_bit(o2::aod::evsel::kIsVertexTRDmatched) && + // collision.selection_bit(aod::evsel::kNoCollInTimeRangeStandard) && + // collision.selection_bit(aod::evsel::kNoCollInRofStandard) && + // collision.selection_bit(aod::evsel::kIsVertexTOFmatched) + { + vars.isAcceptedEvent = true; + } + } else { + if (collision.sel7() && ((doprocessMC) ? true : collision.alias_bit(kINT7))) { vars.isAcceptedEvent = true; } } @@ -368,7 +418,7 @@ void ChargedParticles::initEventMC(const C& collision, const P& particles) if (!initParticle(particle) || !vars.isChargedPrimary) { continue; } - ++vars.multTrue; + vars.multTrue += vars.nRepetitions; } bool isGoodEventClass = (normINELGT0) ? vars.isINELGT0EventMC : (vars.multTrue > 0); vars.isAcceptedEventMC = isGoodEventClass && (std::abs(collision.posZ()) < kMaxVtxZ); @@ -390,9 +440,11 @@ void ChargedParticles::processTrue(const P& particles) for (const auto& particle : particles) { if (initParticle(particle) && vars.isChargedPrimary) { - histos.fill(HIST("multPtSpec_prim_gen"), vars.multTrue, particle.pt()); - if (!vars.isAcceptedEvent) { - histos.fill(HIST("multPtSpec_prim_gen_evtloss"), vars.multTrue, particle.pt()); + for (auto i = 0u; i < vars.nRepetitions; ++i) { + histos.fill(HIST("multPtSpec_prim_gen"), vars.multTrue, particle.pt()); + if (!vars.isAcceptedEvent) { + histos.fill(HIST("multPtSpec_prim_gen_evtloss"), vars.multTrue, particle.pt()); + } } } } @@ -441,7 +493,7 @@ void ChargedParticles::processMeas(const T& tracks) foundParticles.push_back(track.mcParticleId()); - const auto& particle = track.template mcParticle_as(); + const auto& particle = track.template mcParticle_as(); if (!vars.isAcceptedEventMC) { histos.fill(HIST("multPtSpec_trk_meas_evtcont"), vars.multMeas, track.pt()); @@ -450,19 +502,21 @@ void ChargedParticles::processMeas(const T& tracks) histos.fill(HIST("multPtSpec_trk_inter"), vars.multTrue, track.pt()); if (initParticle(particle)) { - if (!vars.isChargedPrimary) { - histos.fill(HIST("multPtSpec_trk_sec_meas"), vars.multMeas, track.pt()); - } else { - histos.fill(HIST("multCorrel_prim"), vars.multMeas, vars.multTrue); - histos.fill(HIST("ptCorrel_prim"), track.pt(), particle.pt()); - histos.fill(HIST("multPtSpec_prim_meas"), vars.multTrue, particle.pt()); - histos.fill(HIST("multPtSpec_trk_prim_meas"), vars.multMeas, track.pt()); + for (auto i = 0u; i < vars.nRepetitions; ++i) { + if (!vars.isChargedPrimary) { + histos.fill(HIST("multPtSpec_trk_sec_meas"), vars.multMeas, track.pt()); + } else { + histos.fill(HIST("multCorrel_prim"), vars.multMeas, vars.multTrue); + histos.fill(HIST("ptCorrel_prim"), track.pt(), particle.pt()); + histos.fill(HIST("multPtSpec_prim_meas"), vars.multTrue, particle.pt()); + histos.fill(HIST("multPtSpec_trk_prim_meas"), vars.multMeas, track.pt()); + } } } } } - std::unordered_set uniqueIndices(foundParticles.begin(), foundParticles.end()); + std::unordered_set uniqueIndices(foundParticles.begin(), foundParticles.end()); for (const auto& mcParticleID : uniqueIndices) { histos.fill(HIST("track_ambiguity"), std::count(foundParticles.begin(), foundParticles.end(), mcParticleID)); } From f99a2263d17d6ff9a715d57333ef9740a871d0c5 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 24 Oct 2025 13:24:33 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGLF/Tasks/Nuspex/chargedParticles.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/chargedParticles.cxx b/PWGLF/Tasks/Nuspex/chargedParticles.cxx index 9408c94c7ce..9da546afd57 100644 --- a/PWGLF/Tasks/Nuspex/chargedParticles.cxx +++ b/PWGLF/Tasks/Nuspex/chargedParticles.cxx @@ -13,12 +13,13 @@ /// \brief Task for analysis of charged particle pt spectra vs multiplicity with 2d unfolding. /// \author Mario Krüger +#include "PWGLF/DataModel/particleCompositionCorrectionTable.h" + #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "PWGLF/DataModel/particleCompositionCorrectionTable.h" #include #include @@ -26,9 +27,9 @@ #include #include +#include #include #include -#include using namespace o2; using namespace o2::framework;