Skip to content

Commit 1b15091

Browse files
LuzhiyonggLuzhiyongg
committed
add event/track selections
1 parent 52ced55 commit 1b15091

File tree

1 file changed

+86
-9
lines changed

1 file changed

+86
-9
lines changed

‎PWGCF/Flow/Tasks/flowMc.cxx‎

Lines changed: 86 additions & 9 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")};
@@ -113,7 +134,18 @@ struct FlowMc {
113134

114135
void init(InitContext&)
115136
{
137+
ccdb->setURL(ccdbUrl.value);
138+
ccdb->setCaching(true);
139+
auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
140+
ccdb->setCreatedNotAfter(now);
141+
142+
const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
143+
const AxisSpec axisEta{20, -1., 1., "#eta"};
144+
const AxisSpec axisCounter{1, 0, +1, ""};
116145
// QA histograms
146+
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter});
147+
histos.add<TH1>("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}});
148+
histos.add<TH1>("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter});
117149
histos.add<TH1>("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}});
118150
histos.add<TH1>("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}});
119151
// pT histograms
@@ -156,7 +188,9 @@ struct FlowMc {
156188
histos.add<TH2>("hPtNchGeneratedLambda", "Reco production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
157189
histos.add<TH2>("hPtNchGlobalLambda", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
158190
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
191+
histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
159192
histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
193+
histos.add<TH3>("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
160194
histos.add<TH1>("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}});
161195

162196
o2::framework::AxisSpec axis = axisPt;
@@ -318,9 +352,35 @@ struct FlowMc {
318352
}
319353
}
320354

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

323-
void process(aod::McCollision const& mcCollision, aod::BCsWithTimestamps const&, soa::Join<aod::McParticles, aod::ParticlesToTracks> const& mcParticles, RecoTracks const&)
383+
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&)
324384
{
325385

326386
float imp = mcCollision.impactParameter();
@@ -337,6 +397,21 @@ struct FlowMc {
337397
auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
338398
loadCorrections(bc.timestamp());
339399

400+
if (collisions.size() > -1) {
401+
histos.fill(HIST("mcEventCounter"), 0.5);
402+
histos.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed
403+
if (cfgRecoEvRejectMC) {
404+
if (collisions.size() != 1) { // only pass those have one reconstruction event
405+
return;
406+
}
407+
for (auto const& collision : collisions) {
408+
if (!eventSelected(collision))
409+
return;
410+
}
411+
}
412+
histos.fill(HIST("RecoEventCounter"), 0.5);
413+
}
414+
340415
if (imp > minB && imp < maxB) {
341416
// event within range
342417
histos.fill(HIST("hImpactParameter"), imp);
@@ -363,9 +438,9 @@ struct FlowMc {
363438
if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance
364439
continue;
365440
if (mcParticle.has_tracks()) {
366-
auto const& tracks = mcParticle.tracks_as<RecoTracks>();
441+
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
367442
for (auto const& track : tracks) {
368-
if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) {
443+
if (!trackSelected(track)) {
369444
continue;
370445
}
371446
if (cfgIsGlobalTrack && track.isGlobalTrack()) {
@@ -418,6 +493,7 @@ struct FlowMc {
418493
histos.fill(HIST("hBVsPtVsPhiGenerated"), imp, deltaPhi, mcParticle.pt());
419494
histos.fill(HIST("hPtNchGenerated"), mcParticle.pt(), nChGlobal);
420495
histos.fill(HIST("hPtMCGen"), mcParticle.pt());
496+
histos.fill(HIST("hEtaPtVtxzMCGen"), mcParticle.eta(), mcParticle.pt(), vtxz);
421497
if (pdgCode == PDG_t::kPiPlus)
422498
histos.fill(HIST("hPtNchGeneratedPion"), mcParticle.pt(), nChGlobal);
423499
if (pdgCode == PDG_t::kKPlus)
@@ -437,9 +513,9 @@ struct FlowMc {
437513
bool validITSTrack = false;
438514
bool validITSABTrack = false;
439515
if (mcParticle.has_tracks()) {
440-
auto const& tracks = mcParticle.tracks_as<RecoTracks>();
516+
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
441517
for (auto const& track : tracks) {
442-
if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) {
518+
if (!trackSelected(track)) {
443519
continue;
444520
}
445521
histos.fill(HIST("hnTPCClu"), track.tpcNClsFound());
@@ -519,6 +595,7 @@ struct FlowMc {
519595
histos.fill(HIST("hBVsPtVsPhiGlobal"), imp, deltaPhi, mcParticle.pt(), wacc * weff);
520596
histos.fill(HIST("hPtNchGlobal"), mcParticle.pt(), nChGlobal);
521597
histos.fill(HIST("hPtMCGlobal"), mcParticle.pt());
598+
histos.fill(HIST("hEtaPtVtxzMCGlobal"), mcParticle.eta(), mcParticle.pt(), vtxz);
522599
if (pdgCode == PDG_t::kPiPlus)
523600
histos.fill(HIST("hPtNchGlobalPion"), mcParticle.pt(), nChGlobal);
524601
if (pdgCode == PDG_t::kKPlus)

0 commit comments

Comments
 (0)