Skip to content

Commit 73f63a7

Browse files
y-songy-song
authored
[PWGJE] Add new function and histograms to check number of tracks according to their MC associations (#17778)
1 parent 94c99d7 commit 73f63a7

1 file changed

Lines changed: 209 additions & 3 deletions

File tree

‎PWGJE/Tasks/trackEfficiency.cxx‎

Lines changed: 209 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -258,7 +258,7 @@ struct TrackEfficiency {
258258
AxisSpec dcaxyAxis = {1000, -1.0, 1.0, "dca_{xy}"};
259259
AxisSpec dcazAxis = {4000, -4.0, 4.0, "dca_{z}"};
260260

261-
if (doprocessEFficiencyPurity || doprocessEFficiencyPurityWeighted) {
261+
if (doprocessEFficiencyPurity || doprocessEFficiencyPurityWeighted || doprocessQcCheck) {
262262

263263
registry.add("hMcCollCutsCounts", "McColl cuts count checks", {HistType::kTH1F, {{10, 0., 10.}}});
264264
registry.get<TH1>(HIST("hMcCollCutsCounts"))->GetXaxis()->SetBinLabel(1, "allMcColl");
@@ -283,7 +283,7 @@ struct TrackEfficiency {
283283
registry.get<TH1>(HIST("hTrackCutsCounts"))->GetXaxis()->SetBinLabel(2, "trackSel");
284284
registry.get<TH1>(HIST("hTrackCutsCounts"))->GetXaxis()->SetBinLabel(3, "hasMcParticle");
285285

286-
if (doprocessEFficiencyPurity) {
286+
if (doprocessEFficiencyPurity || doprocessQcCheck) {
287287
registry.get<TH1>(HIST("hTrackCutsCounts"))->GetXaxis()->SetBinLabel(4, "mcPartIsPrimary");
288288
registry.get<TH1>(HIST("hTrackCutsCounts"))->GetXaxis()->SetBinLabel(5, "etaAcc"); // not actually applied here but it will give an idea of what will be done in the post processing
289289
}
@@ -292,7 +292,13 @@ struct TrackEfficiency {
292292
registry.get<TH1>(HIST("hTrackCutsCounts"))->GetXaxis()->SetBinLabel(5, "mcPartIsPrimary");
293293
registry.get<TH1>(HIST("hTrackCutsCounts"))->GetXaxis()->SetBinLabel(6, "etaAcc"); // not actually applied here but it will give an idea of what will be done in the post processing
294294
}
295-
295+
if (doprocessQcCheck) {
296+
registry.add("h_ntrack_nonassociatedtrack", "Non-associated tracks;N_{tracks};counts", {HistType::kTH1I, {nTracksAxis}});
297+
registry.add("h_ntrack_associatedtrack_primary", "Associated tracks, primary;N_{tracks};counts", {HistType::kTH1I, {nTracksAxis}});
298+
registry.add("h_ntrack_associatedtrack_nonprimary", "Associated tracks, non-primary;N_{tracks};counts", {HistType::kTH1I, {nTracksAxis}});
299+
registry.add("h_ntrack_associatedtrack_split_primary", "Associated split tracks, primary;N_{tracks};counts", {HistType::kTH1I, {nTracksAxis}});
300+
registry.add("h_ntrack_associatedtrack_split_nonprimary", "Associated split tracks, non-primary;N_{tracks};counts", {HistType::kTH1I, {nTracksAxis}});
301+
}
296302
// ptAxisLow
297303
registry.add("h3_particle_pt_particle_eta_particle_phi_mcpartofinterest", "#it{p}_{T, mcpart} (GeV/#it{c}); #eta_{mcpart}; #phi_{mcpart}", {HistType::kTH3F, {ptAxisEff, etaAxisEff, phiAxisEff}});
298304
registry.add("h3_particle_pt_particle_eta_particle_phi_mcpart_nonprimary", "#it{p}_{T, mcpart} (GeV/#it{c}); #eta_{mcpart}; #phi_{mcpart}", {HistType::kTH3F, {ptAxisEff, etaAxisEff, phiAxisEff}});
@@ -1500,6 +1506,206 @@ struct TrackEfficiency {
15001506
}
15011507
}
15021508
PROCESS_SWITCH(TrackEfficiency, processItsTpcMatchingMC, "fills histograms for ITS-TPC matching analysis - MC study, true primary and true secondary separated", false);
1509+
1510+
void processQcCheck(aod::JetMcCollisions::iterator const& mcCollision,
1511+
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions, // smallgroups gives only the collisions associated to the current mccollision, thanks to the mccollisionlabel pre-integrated in jetcollisionsmcd
1512+
soa::Join<aod::JetTracksMCD, aod::JTrackExtras, aod::JTrackPIs> const& jetTracks,
1513+
soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA> const&,
1514+
JetParticlesWithOriginal const& jMcParticles)
1515+
{
1516+
registry.fill(HIST("hMcCollCutsCounts"), 0.5); // all mcCollisions
1517+
1518+
if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
1519+
return;
1520+
}
1521+
registry.fill(HIST("hMcCollCutsCounts"), 1.5); // mcCollision.posZ() condition
1522+
1523+
if (collisions.size() < 1) {
1524+
return;
1525+
}
1526+
registry.fill(HIST("hMcCollCutsCounts"), 2.5); // mcCollisions with at least one reconstructed collision
1527+
1528+
if (acceptSplitCollisions == NonSplitOnly && collisions.size() > 1) {
1529+
return;
1530+
}
1531+
registry.fill(HIST("hMcCollCutsCounts"), 3.5); // split mcCollisions condition
1532+
1533+
float centrality = -1;
1534+
bool hasSel8Coll = false;
1535+
bool centralityCheck = false;
1536+
bool occupancyCheck = false;
1537+
if (acceptSplitCollisions == SplitOkCheckFirstAssocCollOnly || acceptSplitCollisions == NonSplitOnly) { // check only that the first reconstructed collision passes the check (for the NonSplitOnly case, there's only one associated collision)
1538+
if (jetderiveddatautilities::selectCollision(collisions.begin(), eventSelectionBits, skipMBGapEvents, applyRCTSelections)) { // Skipping MC events that have their first associated collision not reconstructed
1539+
hasSel8Coll = true;
1540+
}
1541+
if (!checkOccupancy || ((trackOccupancyInTimeRangeMin < collisions.begin().trackOccupancyInTimeRange()) && (collisions.begin().trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMax))) { // check occupancy only in GP Pb-Pb MC
1542+
occupancyCheck = true;
1543+
}
1544+
centrality = checkCentFT0M ? collisions.begin().centFT0M() : collisions.begin().centFT0C();
1545+
if (!cutCentrality || ((centralityMin < centrality) && (centrality < centralityMax))) { // mcCollision.centFT0C() isn't filled at the moment; can use it instead when it is added to O2Physics
1546+
centralityCheck = true;
1547+
}
1548+
} else if (acceptSplitCollisions == SplitOkCheckAnyAssocColl) { // check that at least one of the reconstructed collisions passes the checks
1549+
for (auto const& collision : collisions) {
1550+
if (jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRCTSelections)) { // Skipping MC events that have not a single selected reconstructed collision ; effect unclear if mcColl is split
1551+
hasSel8Coll = true;
1552+
}
1553+
if (!checkOccupancy || ((trackOccupancyInTimeRangeMin < collision.trackOccupancyInTimeRange()) && (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMax))) { // check occupancy only in GP Pb-Pb MC
1554+
occupancyCheck = true;
1555+
}
1556+
centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C();
1557+
if (!cutCentrality || ((centralityMin < centrality) && (centrality < centralityMax))) { // effect unclear if mcColl is split
1558+
centralityCheck = true;
1559+
}
1560+
}
1561+
}
1562+
if (!hasSel8Coll) {
1563+
return;
1564+
}
1565+
registry.fill(HIST("hMcCollCutsCounts"), 4.5); // at least one of the reconstructed collisions associated with this mcCollision is selected
1566+
1567+
// float centrality = checkCentFT0M ? mcCollision.centFT0M() : mcCollision.centFT0C(); mcCollision.centFT0C() isn't filled at the moment; can be added back when it is
1568+
// if (cutCentrality && (centrality < centralityMin || centralityMax < centrality)) {
1569+
// return;
1570+
// }
1571+
if (!centralityCheck) {
1572+
return;
1573+
}
1574+
registry.fill(HIST("hMcCollCutsCounts"), 5.5); // at least one of the reconstructed collisions associated with this mcCollision is selected with regard to centrality
1575+
1576+
float pTHat = mcCollision.ptHard() < pTHatSettingSentinelValue ? mcCollision.ptHard() : simPtRef / (std::pow(mcCollision.weight(), 1.0 / pTHatExponent));
1577+
if (pTHat < ptHatMin || pTHat > ptHatMax) { // only allows mcCollisions with weight in between min and max
1578+
return;
1579+
}
1580+
registry.fill(HIST("hMcCollCutsCounts"), 6.5); // ptHat condition
1581+
1582+
if (checkOccupancy) {
1583+
if (!occupancyCheck) {
1584+
return;
1585+
}
1586+
registry.fill(HIST("hMcCollCutsCounts"), 7.5);
1587+
}
1588+
1589+
for (auto const& jMcParticle : jMcParticles) {
1590+
registry.fill(HIST("hMcPartCutsCounts"), 0.5); // allPartsInSelMcColl
1591+
1592+
if (!isChargedParticle(jMcParticle.pdgCode())) {
1593+
continue;
1594+
}
1595+
registry.fill(HIST("hMcPartCutsCounts"), 1.5); // isCharged
1596+
1597+
registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_mcpart_nonprimary"), jMcParticle.pt(), jMcParticle.eta(), jMcParticle.phi());
1598+
1599+
if (checkPrimaryPart && !jMcParticle.isPhysicalPrimary()) { // global tracks should be mostly primaries
1600+
continue;
1601+
}
1602+
registry.fill(HIST("hMcPartCutsCounts"), 2.5); // isPrimary
1603+
1604+
registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_mcpartofinterest"), jMcParticle.pt(), jMcParticle.eta(), jMcParticle.phi());
1605+
1606+
registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_mcpartofinterest"), jMcParticle.pt(), jMcParticle.eta(), jMcParticle.phi());
1607+
1608+
if ((std::abs(jMcParticle.eta()) < trackEtaAcceptanceCountQA)) { // removed from actual cuts for now because all the histograms have an eta axis
1609+
registry.fill(HIST("hMcPartCutsCounts"), 3.5); // etaAccept // not actually applied here but it will give an idea of what will be done in the post processing
1610+
}
1611+
}
1612+
1613+
std::vector<int> seenMcParticlesVector; // is reset every mc collision
1614+
1615+
int splitCollCounter = 0;
1616+
for (auto const& collision : collisions) {
1617+
splitCollCounter++;
1618+
if (acceptSplitCollisions == SplitOkCheckFirstAssocCollOnly && splitCollCounter > 1) {
1619+
return;
1620+
}
1621+
1622+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRCTSelections) || !(std::abs(collision.posZ()) < vertexZCut)) {
1623+
continue;
1624+
}
1625+
1626+
auto collTracks = jetTracks.sliceBy(tracksPerJCollision, collision.globalIndex());
1627+
int ntrack_nonassociatedtrack = 0;
1628+
int ntrack_associatedtrack_nonprimary = 0;
1629+
int ntrack_associatedtrack_primary = 0;
1630+
int ntrack_associatedtrack_split_nonprimary = 0;
1631+
int ntrack_associatedtrack_split_primary = 0;
1632+
for (auto const& track : collTracks) {
1633+
registry.fill(HIST("hTrackCutsCounts"), 0.5);
1634+
1635+
if (!isAcceptedTrack(track)) {
1636+
continue;
1637+
}
1638+
registry.fill(HIST("hTrackCutsCounts"), 1.5);
1639+
1640+
if (!track.has_mcParticle()) {
1641+
ntrack_nonassociatedtrack += 1;
1642+
1643+
registry.fill(HIST("h3_track_pt_track_eta_track_phi_nonassociatedtrack"), track.pt(), track.eta(), track.phi());
1644+
1645+
registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_nonassociatedtrack"), track.pt(), track.eta(), track.phi());
1646+
continue;
1647+
}
1648+
registry.fill(HIST("hTrackCutsCounts"), 2.5);
1649+
1650+
auto jMcParticleFromTrack = track.mcParticle_as<JetParticlesWithOriginal>();
1651+
if (!jMcParticleFromTrack.isPhysicalPrimary()) {
1652+
ntrack_associatedtrack_nonprimary += 1;
1653+
1654+
registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_nonprimary"), track.pt(), track.eta(), track.phi());
1655+
registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1656+
1657+
registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_nonprimary"), track.pt(), track.eta(), track.phi());
1658+
registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1659+
1660+
if (std::find(seenMcParticlesVector.begin(), seenMcParticlesVector.end(), jMcParticleFromTrack.globalIndex()) != seenMcParticlesVector.end()) {
1661+
ntrack_associatedtrack_split_nonprimary += 1;
1662+
1663+
registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_split_nonprimary"), track.pt(), track.eta(), track.phi());
1664+
registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_split_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1665+
1666+
registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_split_nonprimary"), track.pt(), track.eta(), track.phi());
1667+
registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_split_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1668+
} else {
1669+
seenMcParticlesVector.push_back(jMcParticleFromTrack.globalIndex());
1670+
}
1671+
1672+
continue;
1673+
}
1674+
1675+
registry.fill(HIST("hTrackCutsCounts"), 3.5);
1676+
1677+
ntrack_associatedtrack_primary += 1;
1678+
registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_primary"), track.pt(), track.eta(), track.phi());
1679+
registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1680+
registry.fill(HIST("h2_particle_pt_track_pt_residual_associatedtrack_primary"), jMcParticleFromTrack.pt(), (jMcParticleFromTrack.pt() - track.pt()) / jMcParticleFromTrack.pt());
1681+
1682+
registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_primary"), track.pt(), track.eta(), track.phi());
1683+
registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1684+
registry.fill(HIST("h2_particle_pt_high_track_pt_high_residual_associatedtrack_primary"), jMcParticleFromTrack.pt(), (jMcParticleFromTrack.pt() - track.pt()) / jMcParticleFromTrack.pt());
1685+
1686+
if (std::find(seenMcParticlesVector.begin(), seenMcParticlesVector.end(), jMcParticleFromTrack.globalIndex()) != seenMcParticlesVector.end()) {
1687+
ntrack_associatedtrack_split_primary += 1;
1688+
registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_split_primary"), track.pt(), track.eta(), track.phi());
1689+
registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_split_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1690+
1691+
registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_split_primary"), track.pt(), track.eta(), track.phi());
1692+
registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_split_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi());
1693+
} else {
1694+
seenMcParticlesVector.push_back(jMcParticleFromTrack.globalIndex());
1695+
}
1696+
1697+
if (std::abs(jMcParticleFromTrack.eta()) < trackEtaAcceptanceCountQA) { // not actually applied here but it will give an idea of what will be done in the post processing
1698+
registry.fill(HIST("hTrackCutsCounts"), 4.5);
1699+
}
1700+
}
1701+
registry.fill(HIST("h_ntrack_nonassociatedtrack"), ntrack_nonassociatedtrack);
1702+
registry.fill(HIST("h_ntrack_associatedtrack_nonprimary"), ntrack_associatedtrack_nonprimary);
1703+
registry.fill(HIST("h_ntrack_associatedtrack_split_nonprimary"), ntrack_associatedtrack_split_nonprimary);
1704+
registry.fill(HIST("h_ntrack_associatedtrack_primary"), ntrack_associatedtrack_primary);
1705+
registry.fill(HIST("h_ntrack_associatedtrack_split_primary"), ntrack_associatedtrack_split_primary);
1706+
}
1707+
}
1708+
PROCESS_SWITCH(TrackEfficiency, processQcCheck, "Histograms for QC checks", false);
15031709
};
15041710

15051711
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)