Skip to content

Commit be2e37d

Browse files
y-songy-song
committed
Add new function and histograms to check number of tracks according to their MC associations
1 parent eb8b678 commit be2e37d

1 file changed

Lines changed: 210 additions & 4 deletions

File tree

‎PWGJE/Tasks/trackEfficiency.cxx‎

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

260-
if (doprocessEFficiencyPurity || doprocessEFficiencyPurityWeighted) {
260+
if (doprocessEFficiencyPurity || doprocessEFficiencyPurityWeighted || doprocessQcCheck) {
261261

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

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

15041710
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
15051711
{
15061712
return WorkflowSpec{adaptAnalysisTask<TrackEfficiency>(cfgc)};
1507-
}
1713+
}

0 commit comments

Comments
 (0)