Skip to content

Commit 50598af

Browse files
christianreckziegelchristianreckziegel
committed
analyseMC() function updated for more robust event selection handling and QA
1 parent 72501b7 commit 50598af

1 file changed

Lines changed: 132 additions & 37 deletions

File tree

‎PWGJE/Tasks/hfFragmentationFunction.cxx‎

Lines changed: 132 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,25 @@ double deltaPhi(double phi1, double phi2)
6565
return std::abs(dphi);
6666
}
6767

68+
//
69+
/// Collision counter selection indexes
70+
///
71+
/// The collision selection is done and stored in multiple steps, for later QA analysis.
72+
/// In order not to hard code which bins should be filled throughout different process
73+
/// function, this namespace with enums is create
74+
namespace collisionSelections
75+
{
76+
enum CollisionSelectionStep {
77+
kMCCollisions = 0, ///< raw mccollisions with no selection, starts with 0
78+
kMCCollisionsZCut, ///< mccollisions with z vtx selection
79+
kMCCollisionsZCutSel8, ///< mccollisions with z vtx and sel8 mc emulated selections
80+
kMCCollisionsZCutSel8HasCollisions, ///< mccollisions with z vtx and sel8 mc emulated selections, with at least one reconstructed collisions
81+
kMCCollisionsZCutSel8SplitCollisions, ///< mccollisions with z vtx and sel8 mc emulated selections, with no split reconstructed collisions
82+
kRecoCollisions, ///< raw reconstructed collisions after previous mccollisions selection
83+
kRecoCollisionsZcut, ///< reconstructed collisions with z vtx selection after previous mccollisions selection
84+
kRecoCollisionsZcutSel8 ///< reconstructed collisions with z vtx and sel8 selections after previous mccollisions selection
85+
};
86+
}
6887
// creating table for storing distance data
6988
namespace o2::aod
7089
{
@@ -199,6 +218,7 @@ struct HfFragmentationFunction {
199218
Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
200219
Configurable<bool> applyMcEventSelection{"applyMcEventSelection", false, "Choose a boolean value"};
201220
Configurable<bool> applyRecoEventSelection{"applyRecoEventSelection", true, "Choose a boolean value"};
221+
Configurable<bool> rejectSplitCollisions{"rejectSplitCollisions", true, "reject generated events associated to more than one reconstructed collision"};
202222

203223
std::vector<int> eventSelectionBits;
204224

@@ -208,12 +228,12 @@ struct HfFragmentationFunction {
208228
eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
209229

210230
// create histograms
211-
// collision system histograms
212-
std::vector<std::string> histLabels = {"mccollisions", "z_cut", "collisions", "sel8"};
231+
// collision counter histograms
232+
std::vector<std::string> histLabels = {"mccollisions", "mccollisions+z_cut", "mccollisions+z_{cut}+sel8", "mccollisions+z_{cut}+sel8+HasCollisions", "mccollisions+z_{cut}+sel8+NoSplitVtx", "collisions", "collisions+z_{cut}", "collisions+z_{cut}+sel8"};
213233
registry.add("h_collision_counter", ";# of collisions;", HistType::kTH1F, {{static_cast<int>(histLabels.size()), 0.0, static_cast<double>(histLabels.size())}});
214-
auto counter = registry.get<TH1>(HIST("h_collision_counter"));
234+
auto collCounter = registry.get<TH1>(HIST("h_collision_counter"));
215235
for (std::vector<std::string>::size_type iCounter = 0; iCounter < histLabels.size(); iCounter++) {
216-
counter->GetXaxis()->SetBinLabel(iCounter + 1, histLabels[iCounter].data());
236+
collCounter->GetXaxis()->SetBinLabel(iCounter + 1, histLabels[iCounter].data());
217237
}
218238
registry.add("h_jet_counter", ";# of jets;", {HistType::kTH1F, {{6, 0., 3.0}}});
219239
auto jetCounter = registry.get<TH1>(HIST("h_jet_counter"));
@@ -245,11 +265,11 @@ struct HfFragmentationFunction {
245265
aod::JetTracks const&)
246266
{
247267
// apply event selection and fill histograms for sanity check
248-
registry.fill(HIST("h_collision_counter"), 2.0);
268+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisions);
249269
if (applyRecoEventSelection && (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut))) {
250270
return;
251271
}
252-
registry.fill(HIST("h_collision_counter"), 3.0);
272+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisionsZcutSel8);
253273

254274
for (const auto& jet : jets) {
255275
// fill jet counter histogram
@@ -320,22 +340,22 @@ struct HfFragmentationFunction {
320340
{
321341
for (const auto& mccollision : mccollisions) {
322342

323-
registry.fill(HIST("h_collision_counter"), 0.0);
343+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisions);
324344
// skip collisions outside of |z| < vertexZCut
325345
if (applyMcEventSelection && (!jetderiveddatautilities::selectCollision(mccollision, eventSelectionBits) || !(std::abs(mccollision.posZ()) < vertexZCut))) {
326346
continue;
327347
}
328-
registry.fill(HIST("h_collision_counter"), 1.0);
348+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisionsZCutSel8);
329349

330350
// reconstructed collisions associated to same mccollision
331351
const auto collisionsPerMCCollision = collisions.sliceBy(collisionsPerMCCollisionPreslice, mccollision.globalIndex());
332352
for (const auto& collision : collisionsPerMCCollision) {
333353

334-
registry.fill(HIST("h_collision_counter"), 2.0);
354+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisions);
335355
if (applyRecoEventSelection && (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut))) {
336356
continue;
337357
}
338-
registry.fill(HIST("h_collision_counter"), 3.0);
358+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisionsZcutSel8);
339359

340360
// d0 detector level jets associated to the current same collision
341361
const auto d0mcdJetsPerCollision = mcdjets.sliceBy(d0MCDJetsPerCollisionPreslice, collision.globalIndex());
@@ -392,25 +412,65 @@ struct HfFragmentationFunction {
392412
}
393413
PROCESS_SWITCH(HfFragmentationFunction, processMcEfficiency, "non-matched and matched MC HF and jets", false);
394414

395-
template <typename TMCPJetsPerMCCollisionPreslice, typename TJetsMCD, typename TJetsMCP, typename TCandidatesMCD, typename TCandidatesMCP>
415+
template <typename TMCPJetsPerMCCollisionPreslice, typename TMCDJetsPerCollisionPreslice, typename TJetsMCP, typename TJetsMCD, typename TCandidatesMCP, typename TCandidatesMCD>
396416
void analyzeMC(TMCPJetsPerMCCollisionPreslice const& MCPJetsPerMCCollisionPreslice,
417+
TMCDJetsPerCollisionPreslice const& MCDJetsPerCollisionPreslice,
397418
aod::JetMcCollisions const& mccollisions,
398419
aod::JetCollisionsMCD const& collisions,
399-
TJetsMCD const&,
400420
TJetsMCP const& mcpjets,
401-
TCandidatesMCD const&,
421+
TJetsMCD const& mcdjets,
402422
TCandidatesMCP const&,
403-
aod::JetTracks const&,
404-
aod::JetParticles const&)
423+
TCandidatesMCD const&,
424+
aod::JetParticles const&,
425+
aod::JetTracks const&)
405426
{
406427
for (const auto& mccollision : mccollisions) {
407-
registry.fill(HIST("h_collision_counter"), 0.0);
428+
429+
// --- begin event selection
430+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisions);
408431
// skip collisions outside of |z| < vertexZCut
409-
if (applyMcEventSelection && (!jetderiveddatautilities::selectCollision(mccollision, eventSelectionBits) || !(std::abs(mccollision.posZ()) < vertexZCut))) {
432+
if (applyMcEventSelection && !(std::abs(mccollision.posZ()) < vertexZCut)) {
433+
continue;
434+
}
435+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisionsZCut);
436+
if (applyMcEventSelection && !jetderiveddatautilities::selectCollision(mccollision, eventSelectionBits)) {
437+
continue;
438+
}
439+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisionsZCutSel8);
440+
441+
// reconstructed collisions associated to same mccollision
442+
const auto collisionsPerMCCollision = collisions.sliceBy(collisionsPerMCCollisionPreslice, mccollision.globalIndex());
443+
// only consider events with at least one reconstructed collision
444+
if (collisionsPerMCCollision.size() == 0) {
445+
continue;
446+
}
447+
// only consider events with no split vertices (one mccollision-to-one collision)
448+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisionsZCutSel8HasCollisions);
449+
if (rejectSplitCollisions && collisionsPerMCCollision.size() > 1) {
410450
continue;
411451
}
412-
registry.fill(HIST("h_collision_counter"), 1.0);
452+
registry.fill(HIST("h_collision_counter"), collisionSelections::kMCCollisionsZCutSel8SplitCollisions);
413453

454+
bool hasSelectedCollision = false;
455+
for (const auto& collision : collisionsPerMCCollision) {
456+
457+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisions);
458+
if (applyRecoEventSelection && !(std::abs(collision.posZ()) < vertexZCut)) {
459+
continue;
460+
}
461+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisionsZcut);
462+
if (applyRecoEventSelection && !jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
463+
continue;
464+
}
465+
registry.fill(HIST("h_collision_counter"), collisionSelections::kRecoCollisionsZcutSel8);
466+
hasSelectedCollision = true;
467+
} // end of collisions loop
468+
469+
if (!hasSelectedCollision) {
470+
continue;
471+
}
472+
473+
// --- begin particle level jets storage
414474
// hf particle level jets associated to same mccollision
415475
const auto mcpJetsPerMCCollision = mcpjets.sliceBy(MCPJetsPerMCCollisionPreslice, mccollision.globalIndex());
416476
for (const auto& mcpjet : mcpJetsPerMCCollision) {
@@ -427,14 +487,6 @@ struct HfFragmentationFunction {
427487
for (const auto& mcdjet : mcpjet.template matchedJetCand_as<TJetsMCD>()) {
428488
registry.fill(HIST("h_jet_counter"), 2.0);
429489

430-
// apply collision sel8 selection on detector level jet's collision
431-
const auto& collision = collisions.iteratorAt(mcdjet.collisionId());
432-
registry.fill(HIST("h_collision_counter"), 2.0);
433-
if (applyRecoEventSelection && (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut))) {
434-
continue;
435-
}
436-
registry.fill(HIST("h_collision_counter"), 3.0);
437-
438490
// obtain leading HF candidate in jet
439491
auto mcdcand = mcdjet.template candidates_first_as<TCandidatesMCD>();
440492

@@ -460,37 +512,80 @@ struct HfFragmentationFunction {
460512
matchJetTable(jetutilities::deltaR(mcpjet, mcpcand), mcpjet.pt(), mcpjet.eta(), mcpjet.phi(), mcpjet.template tracks_as<aod::JetParticles>().size() + mcpjet.template candidates_as<TCandidatesMCP>().size(), // particle level jet
461513
mcpcand.pt(), mcpcand.eta(), mcpcand.phi(), mcpcand.y(), (mcpcand.originMcGen() == RecoDecay::OriginType::Prompt), // particle level HF
462514
-2, -2, -2, -2, -2, // no detector-level jet found
463-
-2, -2, -2, -2, -2, -2, // no detector-level jet found
515+
-2, -2, -2, -2, -2, false, // no detector-level jet found
464516
-2, -2, -2, // no detector-level jet found
465517
-2, -2); // no detector-level jet found
466518
}
467519
} // end of mcpjets loop
520+
521+
// --- begin non-matched detector level jets storage (fake candidates and correlated background if present)
522+
// reconstructed collisions associated to same mccollision
523+
for (const auto& collision : collisionsPerMCCollision) {
524+
525+
// Also apply the reconstructed level collisions selections
526+
if (applyRecoEventSelection && (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut))) {
527+
continue;
528+
}
529+
530+
// d0 detector level jets associated to the current same collision
531+
const auto mcdJetsPerCollision = mcdjets.sliceBy(MCDJetsPerCollisionPreslice, collision.globalIndex());
532+
for (const auto& mcdjet : mcdJetsPerCollision) {
533+
534+
registry.fill(HIST("h_jet_counter"), 0.5);
535+
536+
// obtain leading HF candidate in jet
537+
auto mcdcand = mcdjet.template candidates_first_as<TCandidatesMCD>();
538+
539+
if (mcdjet.has_matchedJetCand()) {
540+
registry.fill(HIST("h_jet_counter"), 1.5);
541+
} else { // store the detector level non-matched candidates
542+
543+
// reflection information for storage: D0 = +1, D0bar = -1, neither = 0
544+
int selectedAs = 0;
545+
546+
// bitwise AND operation: Checks whether BIT(i) is set, regardless of other bits
547+
if (mcdcand.candidateSelFlag() & BIT(0)) { // CandidateSelFlag == BIT(0) -> selected as D0
548+
selectedAs = 1;
549+
} else if (mcdcand.candidateSelFlag() & BIT(1)) { // CandidateSelFlag == BIT(1) -> selected as D0bar
550+
selectedAs = -1;
551+
}
552+
553+
// store matched particle and detector level data in one single table (calculate angular distance in eta-phi plane on the fly)
554+
matchJetTable(-2, -2, -2, -2, -2, // particle level jet
555+
-2, -2, -2, -2, false, // particle level HF
556+
jetutilities::deltaR(mcdjet, mcdcand), mcdjet.pt(), mcdjet.eta(), mcdjet.phi(), mcdjet.template tracks_as<aod::JetTracks>().size() + mcdjet.template candidates_as<TCandidatesMCD>().size(), // detector level jet
557+
mcdcand.pt(), mcdcand.eta(), mcdcand.phi(), mcdcand.m(), mcdcand.y(), (mcdcand.originMcRec() == RecoDecay::OriginType::Prompt), // detector level HF
558+
mcdcand.mlScores()[0], mcdcand.mlScores()[1], mcdcand.mlScores()[2], // Machine Learning PID scores: background, prompt, non-prompt
559+
static_cast<int>(mcdcand.flagMcMatchRec()), selectedAs); // HF = +1, HFbar = -1, neither = 0
560+
}
561+
} // end of non-matched detector level jets loop
562+
} // end of collisions loop
468563
} // end of mccollisions loop
469564
} // end of analyzeMC function
470565

471566
void processD0MC(aod::JetMcCollisions const& mccollisions,
472567
aod::JetCollisionsMCD const& collisions,
473-
JetD0MCDTable const& mcdjets,
474568
JetD0MCPTable const& mcpjets,
475-
aod::CandidatesD0MCD const& mcdcands,
569+
JetD0MCDTable const& mcdjets,
476570
aod::CandidatesD0MCP const& mcpcands,
477-
aod::JetTracks const& jettracks,
478-
aod::JetParticles const& jetparticles)
571+
aod::CandidatesD0MCD const& mcdcands,
572+
aod::JetParticles const& jetparticles,
573+
aod::JetTracks const& jettracks)
479574
{
480-
analyzeMC<Preslice<JetD0MCPTable>, JetD0MCDTable, JetD0MCPTable, aod::CandidatesD0MCD, aod::CandidatesD0MCP>(d0MCPJetsPerMCCollisionPreslice, mccollisions, collisions, mcdjets, mcpjets, mcdcands, mcpcands, jettracks, jetparticles);
575+
analyzeMC<Preslice<JetD0MCPTable>, Preslice<JetD0MCDTable>, JetD0MCPTable, JetD0MCDTable, aod::CandidatesD0MCP, aod::CandidatesD0MCD>(d0MCPJetsPerMCCollisionPreslice, d0MCDJetsPerCollisionPreslice, mccollisions, collisions, mcpjets, mcdjets, mcpcands, mcdcands, jetparticles, jettracks);
481576
}
482577
PROCESS_SWITCH(HfFragmentationFunction, processD0MC, "Store all simulated D0 jets information with matched candidate (if any found)", false);
483578

484579
void processLcMC(aod::JetMcCollisions const& mccollisions,
485580
aod::JetCollisionsMCD const& collisions,
486-
JetLcMCDTable const& mcdjets,
487581
JetLcMCPTable const& mcpjets,
488-
aod::CandidatesLcMCD const& mcdcands,
582+
JetLcMCDTable const& mcdjets,
489583
aod::CandidatesLcMCP const& mcpcands,
490-
aod::JetTracks const& jettracks,
491-
aod::JetParticles const& jetparticles)
584+
aod::CandidatesLcMCD const& mcdcands,
585+
aod::JetParticles const& jetparticles,
586+
aod::JetTracks const& jettracks)
492587
{
493-
analyzeMC<Preslice<JetLcMCPTable>, JetLcMCDTable, JetLcMCPTable, aod::CandidatesLcMCD, aod::CandidatesLcMCP>(lcMCPJetsPerMCCollisionPreslice, mccollisions, collisions, mcdjets, mcpjets, mcdcands, mcpcands, jettracks, jetparticles);
588+
analyzeMC<Preslice<JetLcMCPTable>, Preslice<JetLcMCDTable>, JetLcMCPTable, JetLcMCDTable, aod::CandidatesLcMCP, aod::CandidatesLcMCD>(lcMCPJetsPerMCCollisionPreslice, lcMCDJetsPerCollisionPreslice, mccollisions, collisions, mcpjets, mcdjets, mcpcands, mcdcands, jetparticles, jettracks);
494589
}
495590
PROCESS_SWITCH(HfFragmentationFunction, processLcMC, "Store all simulated Lc jets information with matched candidate (if any found)", false);
496591
};

0 commit comments

Comments
 (0)