Skip to content

Commit 556996a

Browse files
mario-kruegermario-krueger
committed
use PCC in charged particle analysis
1 parent 7c3d72f commit 556996a

File tree

1 file changed

+80
-25
lines changed

1 file changed

+80
-25
lines changed

‎PWGLF/Tasks/Nuspex/chargedParticles.cxx‎

Lines changed: 80 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@
1313
/// \brief Task for analysis of charged particle pt spectra vs multiplicity with 2d unfolding.
1414
/// \author Mario Krüger <mario.kruger@cern.ch>
1515

16+
#include "PWGLF/DataModel/particleCompositionCorrectionTable.h"
17+
1618
#include "Common/Core/TrackSelection.h"
1719
#include "Common/Core/TrackSelectionDefaults.h"
1820
#include "Common/DataModel/Centrality.h"
@@ -25,6 +27,7 @@
2527
#include <Framework/runDataProcessing.h>
2628
#include <ReconstructionDataFormats/Track.h>
2729

30+
#include <random>
2831
#include <unordered_set>
2932
#include <vector>
3033

@@ -38,6 +41,19 @@ using aod::track::TrackSelectionFlags;
3841
//--------------------------------------------------------------------------------------------------
3942
//--------------------------------------------------------------------------------------------------
4043
struct ChargedParticles {
44+
std::random_device rnd;
45+
std::mt19937 gen{rnd()};
46+
std::uniform_real_distribution<float> dist{0.f, 1.f};
47+
uint32_t getNRepetitons(float scalingFactor)
48+
{
49+
uint32_t nRepetitions = static_cast<uint32_t>(scalingFactor);
50+
float rest = scalingFactor - nRepetitions;
51+
if (rest) {
52+
nRepetitions += (dist(gen) <= rest) ? 1u : 0u;
53+
// LOGP(info, "scalingFactor: {} -> {}", scalingFactor, nRepetitions);
54+
}
55+
return nRepetitions;
56+
};
4157

4258
HistogramRegistry histos;
4359
Service<o2::framework::O2DatabasePDG> pdg;
@@ -62,9 +78,11 @@ struct ChargedParticles {
6278
kSystUpDCAxy,
6379
kSystDownDCAz,
6480
kSystUpDCAz,
65-
kSystITSHits,
81+
kSystITSHits, // only relevant for converted data
6682
kSystDownTPCCrossedRows,
67-
kSystUpTPCCrossedRows
83+
kSystUpTPCCrossedRows,
84+
kSystDownPCC = 120,
85+
kSystUpPCC
6886
};
6987
Configurable<uint32_t> systMode{"systMode", kSystNominal, "variation for systematic uncertainties"};
7088
uint16_t trackSelMask{TrackSelectionFlags::kGlobalTrackWoPtEta}; // track selection bitmask (without cut that is being varied)
@@ -80,6 +98,7 @@ struct ChargedParticles {
8098
bool isAcceptedEventMC{false};
8199
bool isINELGT0EventMC{false};
82100
bool isChargedPrimary{false};
101+
uint32_t nRepetitions{1u};
83102
};
84103
VarContainer vars;
85104
static constexpr float kMaxVtxZ = 10.f;
@@ -92,7 +111,7 @@ struct ChargedParticles {
92111
template <typename T>
93112
bool initTrack(const T& track);
94113

95-
template <typename C, typename T>
114+
template <bool IS_MC, typename C, typename T>
96115
void initEvent(const C& collision, const T& tracks);
97116

98117
template <typename C, typename P>
@@ -112,9 +131,9 @@ struct ChargedParticles {
112131
using CollisionTableMCTrue = aod::McCollisions;
113132
using CollisionTableMC = soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels>>;
114133
using TrackTableMC = soa::Join<aod::FullTracks, aod::TracksDCA, aod::TrackSelection, aod::McTrackLabels>;
115-
using ParticleTableMC = aod::McParticles;
134+
using ParticleTableMC = soa::Join<aod::McParticles, aod::ParticleCompositionCorrection>;
116135
Preslice<TrackTableMC> perCollision = aod::track::collisionId;
117-
void processMC(CollisionTableMCTrue::iterator const& mcCollision, CollisionTableMC const& collisions, TrackTableMC const& tracks, ParticleTableMC const& particles);
136+
void processMC(CollisionTableMCTrue::iterator const& mcCollision, TrackTableMC const& tracks, CollisionTableMC const& collisions, ParticleTableMC const& particles);
118137
PROCESS_SWITCH(ChargedParticles, processMC, "process mc", true);
119138
};
120139

@@ -241,7 +260,7 @@ void ChargedParticles::init(InitContext const&)
241260
//**************************************************************************************************
242261
void ChargedParticles::processData(CollisionTableData::iterator const& collision, TrackTableData const& tracks)
243262
{
244-
initEvent(collision, tracks);
263+
initEvent<false>(collision, tracks);
245264
processMeas<false>(tracks);
246265
}
247266

@@ -250,7 +269,7 @@ void ChargedParticles::processData(CollisionTableData::iterator const& collision
250269
* Entrypoint to processes mc.
251270
*/
252271
//**************************************************************************************************
253-
void ChargedParticles::processMC(CollisionTableMCTrue::iterator const& mcCollision, CollisionTableMC const& collisions, TrackTableMC const& tracks, ParticleTableMC const& particles)
272+
void ChargedParticles::processMC(CollisionTableMCTrue::iterator const& mcCollision, TrackTableMC const& tracks, CollisionTableMC const& collisions, ParticleTableMC const& particles)
254273
{
255274
histos.fill(HIST("collision_ambiguity"), collisions.size());
256275

@@ -270,7 +289,7 @@ void ChargedParticles::processMC(CollisionTableMCTrue::iterator const& mcCollisi
270289
} else {
271290
for (const auto& collision : collisions) {
272291
auto curTracks = tracks.sliceBy(perCollision, collision.globalIndex());
273-
initEvent(collision, curTracks);
292+
initEvent<true>(collision, curTracks);
274293
processMeas<true>(curTracks);
275294
break; // for now look only at first collision...
276295
}
@@ -302,6 +321,15 @@ bool ChargedParticles::initParticle(const P& particle)
302321
if (particle.pt() <= ptMinCut || particle.pt() >= ptMaxCut) {
303322
return false;
304323
}
324+
325+
float pccWeight = particle.pccWeight();
326+
if (systMode == kSystDownPCC) {
327+
pccWeight = particle.pccWeightSysDown();
328+
} else if (systMode == kSystUpPCC) {
329+
pccWeight = particle.pccWeightSysUp();
330+
}
331+
vars.nRepetitions = getNRepetitons(pccWeight);
332+
// FIXME: in case of nRepetitions = 0 INELGT0 can be wrong
305333
return true;
306334
}
307335

@@ -334,20 +362,43 @@ bool ChargedParticles::initTrack(const T& track)
334362
* Check if event is good.
335363
*/
336364
//**************************************************************************************************
337-
template <typename C, typename T>
365+
template <bool IS_MC, typename C, typename T>
338366
void ChargedParticles::initEvent(const C& collision, const T& tracks)
339367
{
340368
vars.multMeas = 0;
341369
for (const auto& track : tracks) {
342370
if (initTrack(track)) {
343-
++vars.multMeas;
371+
if constexpr (IS_MC) {
372+
if (!track.has_mcParticle()) {
373+
continue;
374+
}
375+
const auto& particle = track.template mcParticle_as<ParticleTableMC>();
376+
if (!initParticle(particle)) {
377+
continue;
378+
}
379+
vars.multMeas += vars.nRepetitions;
380+
} else {
381+
++vars.multMeas;
382+
}
344383
}
345384
}
346385

347386
vars.isAcceptedEvent = false;
348387
if (std::abs(collision.posZ()) < kMaxVtxZ) {
349-
if (isRun3 ? collision.sel8() : collision.sel7()) {
350-
if ((isRun3 || doprocessMC) ? true : collision.alias_bit(kINT7)) {
388+
if (isRun3) {
389+
if (collision.sel8() &&
390+
collision.selection_bit(aod::evsel::kNoSameBunchPileup) &&
391+
collision.selection_bit(aod::evsel::kIsVertexITSTPC) &&
392+
collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))
393+
// !collision.selection_bit(o2::aod::evsel::kIsVertexTRDmatched) &&
394+
// collision.selection_bit(aod::evsel::kNoCollInTimeRangeStandard) &&
395+
// collision.selection_bit(aod::evsel::kNoCollInRofStandard) &&
396+
// collision.selection_bit(aod::evsel::kIsVertexTOFmatched)
397+
{
398+
vars.isAcceptedEvent = true;
399+
}
400+
} else {
401+
if (collision.sel7() && ((doprocessMC) ? true : collision.alias_bit(kINT7))) {
351402
vars.isAcceptedEvent = true;
352403
}
353404
}
@@ -368,7 +419,7 @@ void ChargedParticles::initEventMC(const C& collision, const P& particles)
368419
if (!initParticle(particle) || !vars.isChargedPrimary) {
369420
continue;
370421
}
371-
++vars.multTrue;
422+
vars.multTrue += vars.nRepetitions;
372423
}
373424
bool isGoodEventClass = (normINELGT0) ? vars.isINELGT0EventMC : (vars.multTrue > 0);
374425
vars.isAcceptedEventMC = isGoodEventClass && (std::abs(collision.posZ()) < kMaxVtxZ);
@@ -390,9 +441,11 @@ void ChargedParticles::processTrue(const P& particles)
390441

391442
for (const auto& particle : particles) {
392443
if (initParticle(particle) && vars.isChargedPrimary) {
393-
histos.fill(HIST("multPtSpec_prim_gen"), vars.multTrue, particle.pt());
394-
if (!vars.isAcceptedEvent) {
395-
histos.fill(HIST("multPtSpec_prim_gen_evtloss"), vars.multTrue, particle.pt());
444+
for (auto i = 0u; i < vars.nRepetitions; ++i) {
445+
histos.fill(HIST("multPtSpec_prim_gen"), vars.multTrue, particle.pt());
446+
if (!vars.isAcceptedEvent) {
447+
histos.fill(HIST("multPtSpec_prim_gen_evtloss"), vars.multTrue, particle.pt());
448+
}
396449
}
397450
}
398451
}
@@ -441,7 +494,7 @@ void ChargedParticles::processMeas(const T& tracks)
441494

442495
foundParticles.push_back(track.mcParticleId());
443496

444-
const auto& particle = track.template mcParticle_as<aod::McParticles>();
497+
const auto& particle = track.template mcParticle_as<ParticleTableMC>();
445498

446499
if (!vars.isAcceptedEventMC) {
447500
histos.fill(HIST("multPtSpec_trk_meas_evtcont"), vars.multMeas, track.pt());
@@ -450,19 +503,21 @@ void ChargedParticles::processMeas(const T& tracks)
450503

451504
histos.fill(HIST("multPtSpec_trk_inter"), vars.multTrue, track.pt());
452505
if (initParticle(particle)) {
453-
if (!vars.isChargedPrimary) {
454-
histos.fill(HIST("multPtSpec_trk_sec_meas"), vars.multMeas, track.pt());
455-
} else {
456-
histos.fill(HIST("multCorrel_prim"), vars.multMeas, vars.multTrue);
457-
histos.fill(HIST("ptCorrel_prim"), track.pt(), particle.pt());
458-
histos.fill(HIST("multPtSpec_prim_meas"), vars.multTrue, particle.pt());
459-
histos.fill(HIST("multPtSpec_trk_prim_meas"), vars.multMeas, track.pt());
506+
for (auto i = 0u; i < vars.nRepetitions; ++i) {
507+
if (!vars.isChargedPrimary) {
508+
histos.fill(HIST("multPtSpec_trk_sec_meas"), vars.multMeas, track.pt());
509+
} else {
510+
histos.fill(HIST("multCorrel_prim"), vars.multMeas, vars.multTrue);
511+
histos.fill(HIST("ptCorrel_prim"), track.pt(), particle.pt());
512+
histos.fill(HIST("multPtSpec_prim_meas"), vars.multTrue, particle.pt());
513+
histos.fill(HIST("multPtSpec_trk_prim_meas"), vars.multMeas, track.pt());
514+
}
460515
}
461516
}
462517
}
463518
}
464519

465-
std::unordered_set<int> uniqueIndices(foundParticles.begin(), foundParticles.end());
520+
std::unordered_set<int32_t> uniqueIndices(foundParticles.begin(), foundParticles.end());
466521
for (const auto& mcParticleID : uniqueIndices) {
467522
histos.fill(HIST("track_ambiguity"), std::count(foundParticles.begin(), foundParticles.end(), mcParticleID));
468523
}

0 commit comments

Comments
 (0)