Skip to content

Commit 3786573

Browse files
LuzhiyonggLuzhiyongg
authored
[PWGCF] flowMc: add Ev/Tr sel; flowTask/Dihadron: add multCorr config, roll back to wo PID (#12245)
1 parent 02cb80b commit 3786573

File tree

3 files changed

+241
-141
lines changed

3 files changed

+241
-141
lines changed

‎PWGCF/Flow/Tasks/flowMc.cxx‎

Lines changed: 89 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -54,13 +54,19 @@ struct FlowMc {
5454

5555
Configurable<float> minB{"minB", 0.0f, "min impact parameter"};
5656
Configurable<float> maxB{"maxB", 20.0f, "max impact parameter"};
57+
O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
58+
O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for tracks")
59+
O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 1000.0f, "Maximal pT for tracks")
60+
O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
5761
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
5862
O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks")
5963
O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks")
6064
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks")
6165
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks")
62-
O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 70.0f, "minimum TPC clusters")
63-
O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 6.0f, "minimum ITS clusters")
66+
O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters")
67+
O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters")
68+
O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters")
69+
O2_DEFINE_CONFIGURABLE(cfgCutTPCcrossedrows, float, 70.0f, "minimum TPC crossed rows")
6470
O2_DEFINE_CONFIGURABLE(cfgFlowAcceptance, std::string, "", "CCDB path to acceptance object")
6571
O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object")
6672
O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth")
@@ -70,6 +76,10 @@ struct FlowMc {
7076
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples")
7177
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")
7278
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrSlopeFactor, float, 1.0f, "A factor to scale the track density efficiency slope")
79+
O2_DEFINE_CONFIGURABLE(cfgRecoEvRejectMC, bool, false, "reject both MC and Reco events when reco do not pass")
80+
O2_DEFINE_CONFIGURABLE(cfgRecoEvSel8, bool, false, "require sel8 for reconstruction events")
81+
O2_DEFINE_CONFIGURABLE(cfgRecoEvkIsGoodITSLayersAll, bool, false, "require kIsGoodITSLayersAll for reconstruction events")
82+
O2_DEFINE_CONFIGURABLE(cfgRecoEvkNoSameBunchPileup, bool, false, "require kNoSameBunchPileup for reconstruction events")
7383
Configurable<std::vector<double>> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector<double>{0.6003720411, 0.6152630970, 0.6288860646, 0.6360694031, 0.6409494798, 0.6450540203, 0.6482117301, 0.6512592056, 0.6640008690, 0.6862631416, 0.7005738691, 0.7106567432, 0.7170728333}, "parameter 0 for track density efficiency correction"};
7484
Configurable<std::vector<double>> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector<double>{-1.007592e-05, -8.932635e-06, -9.114538e-06, -1.054818e-05, -1.220212e-05, -1.312304e-05, -1.376433e-05, -1.412813e-05, -1.289562e-05, -1.050065e-05, -8.635725e-06, -7.380821e-06, -6.201250e-06}, "parameter 1 for track density efficiency correction"};
7585
float maxEta = 0.8;
@@ -81,6 +91,18 @@ struct FlowMc {
8191

8292
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"};
8393

94+
// Filter for MCcollisions
95+
Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex;
96+
using FilteredMcCollisions = soa::Filtered<aod::McCollisions>;
97+
// Filter for MCParticle
98+
Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax);
99+
using FilteredMcParticles = soa::Filtered<soa::Join<aod::McParticles, aod::ParticlesToTracks>>;
100+
// Filter for reco tracks
101+
Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls);
102+
using FilteredTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::McTrackLabels>>;
103+
104+
// using FilteredTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TrackSelection>;
105+
84106
// Cent vs IP
85107
TH1D* mCentVsIPTruth = nullptr;
86108
bool centVsIPTruthLoaded = false;
@@ -98,7 +120,6 @@ struct FlowMc {
98120

99121
// Connect to ccdb
100122
Service<ccdb::BasicCCDBManager> ccdb;
101-
Configurable<int64_t> ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
102123
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
103124

104125
OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
@@ -110,10 +131,22 @@ struct FlowMc {
110131
std::vector<GFW::CorrConfig> corrconfigsTruth;
111132
std::vector<GFW::CorrConfig> corrconfigsReco;
112133
TRandom3* fRndm = new TRandom3(0);
134+
double epsilon = 1e-6;
113135

114136
void init(InitContext&)
115137
{
138+
ccdb->setURL(ccdbUrl.value);
139+
ccdb->setCaching(true);
140+
auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
141+
ccdb->setCreatedNotAfter(now);
142+
143+
const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
144+
const AxisSpec axisEta{20, -1., 1., "#eta"};
145+
const AxisSpec axisCounter{1, 0, +1, ""};
116146
// QA histograms
147+
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter});
148+
histos.add<TH1>("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}});
149+
histos.add<TH1>("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter});
117150
histos.add<TH1>("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}});
118151
histos.add<TH1>("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}});
119152
// pT histograms
@@ -156,7 +189,9 @@ struct FlowMc {
156189
histos.add<TH2>("hPtNchGeneratedLambda", "Reco production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
157190
histos.add<TH2>("hPtNchGlobalLambda", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
158191
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
192+
histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
159193
histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
194+
histos.add<TH3>("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
160195
histos.add<TH1>("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}});
161196

162197
o2::framework::AxisSpec axis = axisPt;
@@ -318,9 +353,35 @@ struct FlowMc {
318353
}
319354
}
320355

321-
using RecoTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TrackSelection>;
356+
template <typename TCollision>
357+
bool eventSelected(TCollision collision)
358+
{
359+
if (std::fabs(collision.posZ()) > cfgCutVertex) {
360+
return 0;
361+
}
362+
if (cfgRecoEvSel8 && !collision.sel8()) {
363+
return 0;
364+
}
365+
if (cfgRecoEvkNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
366+
// rejects collisions which are associated with the same "found-by-T0" bunch crossing
367+
// https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof
368+
return 0;
369+
}
370+
if (cfgRecoEvkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
371+
// from Jan 9 2025 AOT meeting
372+
// cut time intervals with dead ITS staves
373+
return 0;
374+
}
375+
return 1;
376+
}
377+
378+
template <typename TTrack>
379+
bool trackSelected(TTrack track)
380+
{
381+
return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgCutTPCcrossedrows) && (track.itsNCls() >= cfgCutITSclu));
382+
}
322383

323-
void process(aod::McCollision const& mcCollision, aod::BCsWithTimestamps const&, soa::Join<aod::McParticles, aod::ParticlesToTracks> const& mcParticles, RecoTracks const&)
384+
void process(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels>> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&)
324385
{
325386

326387
float imp = mcCollision.impactParameter();
@@ -337,6 +398,21 @@ struct FlowMc {
337398
auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
338399
loadCorrections(bc.timestamp());
339400

401+
if (collisions.size() > -1) {
402+
histos.fill(HIST("mcEventCounter"), 0.5);
403+
histos.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed
404+
if (cfgRecoEvRejectMC) {
405+
if (collisions.size() != 1) { // only pass those have one reconstruction event
406+
return;
407+
}
408+
for (auto const& collision : collisions) {
409+
if (!eventSelected(collision))
410+
return;
411+
}
412+
}
413+
histos.fill(HIST("RecoEventCounter"), 0.5);
414+
}
415+
340416
if (imp > minB && imp < maxB) {
341417
// event within range
342418
histos.fill(HIST("hImpactParameter"), imp);
@@ -363,9 +439,9 @@ struct FlowMc {
363439
if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance
364440
continue;
365441
if (mcParticle.has_tracks()) {
366-
auto const& tracks = mcParticle.tracks_as<RecoTracks>();
442+
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
367443
for (auto const& track : tracks) {
368-
if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) {
444+
if (!trackSelected(track)) {
369445
continue;
370446
}
371447
if (cfgIsGlobalTrack && track.isGlobalTrack()) {
@@ -418,6 +494,7 @@ struct FlowMc {
418494
histos.fill(HIST("hBVsPtVsPhiGenerated"), imp, deltaPhi, mcParticle.pt());
419495
histos.fill(HIST("hPtNchGenerated"), mcParticle.pt(), nChGlobal);
420496
histos.fill(HIST("hPtMCGen"), mcParticle.pt());
497+
histos.fill(HIST("hEtaPtVtxzMCGen"), mcParticle.eta(), mcParticle.pt(), vtxz);
421498
if (pdgCode == PDG_t::kPiPlus)
422499
histos.fill(HIST("hPtNchGeneratedPion"), mcParticle.pt(), nChGlobal);
423500
if (pdgCode == PDG_t::kKPlus)
@@ -437,9 +514,9 @@ struct FlowMc {
437514
bool validITSTrack = false;
438515
bool validITSABTrack = false;
439516
if (mcParticle.has_tracks()) {
440-
auto const& tracks = mcParticle.tracks_as<RecoTracks>();
517+
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
441518
for (auto const& track : tracks) {
442-
if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) {
519+
if (!trackSelected(track)) {
443520
continue;
444521
}
445522
histos.fill(HIST("hnTPCClu"), track.tpcNClsFound());
@@ -456,10 +533,10 @@ struct FlowMc {
456533
if (track.hasTPC()) {
457534
validTPCTrack = true;
458535
}
459-
if (track.hasITS() && track.itsChi2NCl() > -1e-6) {
536+
if (track.hasITS() && track.itsChi2NCl() > -1. * epsilon) {
460537
validITSTrack = true;
461538
}
462-
if (track.hasITS() && track.itsChi2NCl() < -1e-6) {
539+
if (track.hasITS() && track.itsChi2NCl() < -1. * epsilon) {
463540
validITSABTrack = true;
464541
}
465542
}
@@ -519,6 +596,7 @@ struct FlowMc {
519596
histos.fill(HIST("hBVsPtVsPhiGlobal"), imp, deltaPhi, mcParticle.pt(), wacc * weff);
520597
histos.fill(HIST("hPtNchGlobal"), mcParticle.pt(), nChGlobal);
521598
histos.fill(HIST("hPtMCGlobal"), mcParticle.pt());
599+
histos.fill(HIST("hEtaPtVtxzMCGlobal"), mcParticle.eta(), mcParticle.pt(), vtxz);
522600
if (pdgCode == PDG_t::kPiPlus)
523601
histos.fill(HIST("hPtNchGlobalPion"), mcParticle.pt(), nChGlobal);
524602
if (pdgCode == PDG_t::kKPlus)

0 commit comments

Comments
 (0)