Skip to content

Commit caa29a2

Browse files
f3schf3sch
committed
ITS: TrackExtensionStudy add missed/empty cluster patterns
Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>
1 parent 7a1e4e4 commit caa29a2

File tree

1 file changed

+88
-81
lines changed

1 file changed

+88
-81
lines changed

‎Detectors/ITSMFT/ITS/postprocessing/studies/src/TrackExtension.cxx‎

Lines changed: 88 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -101,12 +101,22 @@ class TrackExtensionStudy : public Task
101101
bool mWithTree{false};
102102

103103
std::unique_ptr<TH1D> mHTrackCounts;
104-
std::unique_ptr<TH1D> mHClustersCounts;
105104
std::unique_ptr<TH1D> mHLengthAny, mHLengthGood, mHLengthFake;
106105
std::unique_ptr<TH1D> mHChi2Any, mHChi2Good, mHChi2Fake;
107106
std::unique_ptr<TH1D> mHPtAny, mHPtGood, mHPtFake;
108107
std::unique_ptr<TH1D> mHExtensionAny, mHExtensionGood, mHExtensionFake;
109-
std::unique_ptr<TH2D> mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake;
108+
std::unique_ptr<TH2D> mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake, mHExtensionPatternsGoodMissed, mHExtensionPatternsGoodEmpty;
109+
;
110+
111+
template <class T, typename... C, typename... F>
112+
std::unique_ptr<T> createHistogram(C... n, F... b)
113+
{
114+
auto t = std::make_unique<T>(n..., b...);
115+
mHistograms.push_back(static_cast<TH1*>(t.get()));
116+
return std::move(t);
117+
}
118+
119+
std::vector<TH1*> mHistograms;
110120
};
111121

112122
void TrackExtensionStudy::init(InitContext& ic)
@@ -120,7 +130,7 @@ void TrackExtensionStudy::init(InitContext& ic)
120130
auto xbins = helpers::makeLogBinning(effHistBins, effPtCutLow, effPtCutHigh);
121131

122132
// Track Counting
123-
mHTrackCounts = std::make_unique<TH1D>("hTrackCounts", "Track Stats", 10, 0, 10);
133+
mHTrackCounts = createHistogram<TH1D>("hTrackCounts", "Track Stats", 10, 0, 10);
124134
mHTrackCounts->GetXaxis()->SetBinLabel(1, "Total Tracks");
125135
mHTrackCounts->GetXaxis()->SetBinLabel(2, "Normal ANY Tracks");
126136
mHTrackCounts->GetXaxis()->SetBinLabel(3, "Normal GOOD Tracks");
@@ -132,49 +142,49 @@ void TrackExtensionStudy::init(InitContext& ic)
132142
mHTrackCounts->GetXaxis()->SetBinLabel(9, "Extended FAKE AFTER Tracks");
133143
mHTrackCounts->GetXaxis()->SetBinLabel(10, "Extended FAKE BEFORE&AFTER Tracks");
134144

135-
// Cluster Counting
136-
mHClustersCounts = std::make_unique<TH1D>("hClusterCounts", "Cluster Stats", 5, 0, 5);
137-
mHClustersCounts->GetXaxis()->SetBinLabel(1, "Total Clusters");
138-
mHClustersCounts->GetXaxis()->SetBinLabel(2, "Tracking");
139-
mHClustersCounts->GetXaxis()->SetBinLabel(3, "Extension");
140-
mHClustersCounts->GetXaxis()->SetBinLabel(4, "Good Extension");
141-
mHClustersCounts->GetXaxis()->SetBinLabel(5, "Fake Extension");
142-
143145
// Length
144-
mHLengthAny = std::make_unique<TH1D>("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8);
145-
mHLengthGood = std::make_unique<TH1D>("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8);
146-
mHLengthFake = std::make_unique<TH1D>("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8);
146+
mHLengthAny = createHistogram<TH1D>("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8);
147+
mHLengthGood = createHistogram<TH1D>("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8);
148+
mHLengthFake = createHistogram<TH1D>("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8);
147149

148150
// Chi2
149-
mHChi2Any = std::make_unique<TH1D>("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100);
150-
mHChi2Good = std::make_unique<TH1D>("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100);
151-
mHChi2Fake = std::make_unique<TH1D>("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100);
151+
mHChi2Any = createHistogram<TH1D>("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100);
152+
mHChi2Good = createHistogram<TH1D>("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100);
153+
mHChi2Fake = createHistogram<TH1D>("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100);
152154

153155
// Pt
154-
mHPtAny = std::make_unique<TH1D>("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh);
155-
mHPtGood = std::make_unique<TH1D>("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh);
156-
mHPtFake = std::make_unique<TH1D>("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh);
156+
mHPtAny = createHistogram<TH1D>("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh);
157+
mHPtGood = createHistogram<TH1D>("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh);
158+
mHPtFake = createHistogram<TH1D>("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", xbins.size(), effPtCutLow, effPtCutHigh);
157159

158160
// Length
159-
mHExtensionAny = std::make_unique<TH1D>("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7);
160-
mHExtensionGood = std::make_unique<TH1D>("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7);
161-
mHExtensionFake = std::make_unique<TH1D>("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7);
161+
mHExtensionAny = createHistogram<TH1D>("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7);
162+
mHExtensionGood = createHistogram<TH1D>("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7);
163+
mHExtensionFake = createHistogram<TH1D>("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7);
162164

163165
// Patterns
164-
auto makePatternAxisLabels = [&](TH1* h) {
166+
auto makePatternAxisLabels = [&](TH1* h, bool xBefore = true) {
165167
for (int i{1}; i <= h->GetXaxis()->GetNbins(); ++i) {
166-
h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str());
168+
if (xBefore) {
169+
h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str());
170+
} else {
171+
h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str());
172+
}
167173
}
168174
for (int i{1}; i <= h->GetYaxis()->GetNbins(); ++i) {
169175
h->GetYaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str());
170176
}
171177
};
172-
mHExtensionPatternsAny = std::make_unique<TH2D>("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
178+
mHExtensionPatternsAny = createHistogram<TH2D>("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
173179
makePatternAxisLabels(mHExtensionPatternsAny.get());
174-
mHExtensionPatternsGood = std::make_unique<TH2D>("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
180+
mHExtensionPatternsGood = createHistogram<TH2D>("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
175181
makePatternAxisLabels(mHExtensionPatternsGood.get());
176-
mHExtensionPatternsFake = std::make_unique<TH2D>("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
182+
mHExtensionPatternsFake = createHistogram<TH2D>("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
177183
makePatternAxisLabels(mHExtensionPatternsFake.get());
184+
mHExtensionPatternsGoodMissed = createHistogram<TH2D>("hExtensionPatternsGoodMissed", "Extended Tracks Pattern (GOOD) Missed Clusters;After;Missed;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
185+
makePatternAxisLabels(mHExtensionPatternsGoodMissed.get(), false);
186+
mHExtensionPatternsGoodEmpty = createHistogram<TH2D>("hExtensionPatternsGoodEmpty", "Extended Tracks Pattern (GOOD) Empty Clusters;Before;After;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
187+
makePatternAxisLabels(mHExtensionPatternsGoodEmpty.get(), false);
178188

179189
mStream = std::make_unique<utils::TreeStreamRedirector>(mOutFileName.c_str(), "RECREATE");
180190
}
@@ -192,8 +202,6 @@ void TrackExtensionStudy::run(ProcessingContext& pc)
192202
mClustersMCLCont = recoData.getITSClustersMCLabels();
193203
mInputITSidxs = recoData.getITSTracksClusterRefs();
194204

195-
mHClustersCounts->SetBinContent(1, mClusters.size());
196-
197205
LOGP(info, "** Found in {} rofs:\n\t- {} clusters with {} labels\n\t- {} tracks with {} labels",
198206
mTracksROFRecords.size(), mClusters.size(), mClustersMCLCont->getIndexedSize(), mTracks.size(), mTracksMCLabels.size());
199207
LOGP(info, "** Found {} sources from kinematic files", mKineReader->getNSources());
@@ -302,6 +310,7 @@ void TrackExtensionStudy::process()
302310
continue;
303311
}
304312
const auto& trk = part.track;
313+
bool isGood = part.isReco && !part.isFake;
305314
mHTrackCounts->Fill(0);
306315

307316
std::bitset<7> extPattern{0};
@@ -321,25 +330,11 @@ void TrackExtensionStudy::process()
321330
continue;
322331
}
323332

324-
if (mWithTree) {
325-
std::array<float, 3> xyz{(float)part.mcTrack.GetStartVertexCoordinatesX(), (float)part.mcTrack.GetStartVertexCoordinatesY(), (float)part.mcTrack.GetStartVertexCoordinatesZ()};
326-
std::array<float, 3> pxyz{(float)part.mcTrack.GetStartVertexMomentumX(), (float)part.mcTrack.GetStartVertexMomentumY(), (float)part.mcTrack.GetStartVertexMomentumZ()};
327-
auto pdg = O2DatabasePDG::Instance()->GetParticle(part.pdg);
328-
if (pdg == nullptr) {
329-
LOGP(fatal, "MC info not available");
330-
}
331-
auto mcTrk = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pdg->Charge() / 3.), true);
332-
(*mStream) << "tree"
333-
<< "trk=" << trk
334-
<< "mcTrk=" << mcTrk
335-
<< "\n";
336-
}
337-
338333
mHTrackCounts->Fill(4);
339334
mHLengthAny->Fill(trk.getNClusters());
340335
mHChi2Any->Fill(trk.getChi2());
341336
mHPtAny->Fill(trk.getPt());
342-
if (part.isReco || !part.isFake) {
337+
if (isGood) {
343338
mHTrackCounts->Fill(5);
344339
mHLengthGood->Fill(trk.getNClusters());
345340
mHChi2Good->Fill(trk.getChi2());
@@ -356,26 +351,21 @@ void TrackExtensionStudy::process()
356351
if (extPattern.test(iLayer)) {
357352
extPattern.set(iLayer);
358353
mHExtensionAny->Fill(iLayer);
359-
if (part.isReco || !part.isFake) {
354+
if (isGood) {
360355
mHExtensionGood->Fill(iLayer);
361356
} else {
362357
mHExtensionFake->Fill(iLayer);
363358
}
364-
365-
if ((part.fakeClusters & (0x1 << iLayer)) == 0) {
366-
mHClustersCounts->Fill(4);
367-
} else {
368-
mHClustersCounts->Fill(5);
369-
}
370359
}
371360
}
372-
std::bitset<7> oldPattern{clusPattern & ~extPattern};
361+
std::bitset<7> oldPattern{clusPattern & ~extPattern}, holePattern{clusPattern};
362+
holePattern.flip();
373363
auto clusN = clusPattern.to_ulong();
374364
auto clusIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), clusN));
375365
auto oldN = oldPattern.to_ulong();
376366
auto oldIdx = std::distance(std::begin(mBitPatternsBefore), std::find(std::begin(mBitPatternsBefore), std::end(mBitPatternsBefore), oldN));
377367
mHExtensionPatternsAny->Fill(oldIdx, clusIdx);
378-
if (part.isReco || !part.isFake) {
368+
if (isGood) {
379369
mHExtensionPatternsGood->Fill(oldIdx, clusIdx);
380370
} else {
381371
mHExtensionPatternsFake->Fill(oldIdx, clusIdx);
@@ -392,7 +382,6 @@ void TrackExtensionStudy::process()
392382
}
393383
}
394384
}
395-
396385
if (oldFake && newFake) {
397386
mHTrackCounts->Fill(9);
398387
} else if (oldFake) {
@@ -401,8 +390,47 @@ void TrackExtensionStudy::process()
401390
mHTrackCounts->Fill(8);
402391
}
403392

404-
mHClustersCounts->SetBinContent(2, mHClustersCounts->GetBinContent(2) + oldPattern.count());
405-
mHClustersCounts->SetBinContent(3, mHClustersCounts->GetBinContent(3) + extPattern.count());
393+
// Check if we missed some clusters
394+
if (isGood && holePattern.any()) {
395+
auto missPattern{clusPattern}, emptyPattern{clusPattern};
396+
for (int iLayer{0}; iLayer < 7; ++iLayer) {
397+
if (!holePattern.test(iLayer)) {
398+
continue;
399+
}
400+
401+
// Check if there was actually a cluster that we missed
402+
if ((part.clusters & (1 << iLayer)) != 0) {
403+
missPattern.set(iLayer);
404+
} else {
405+
emptyPattern.set(iLayer);
406+
}
407+
}
408+
409+
if (missPattern != clusPattern) {
410+
auto missN = missPattern.to_ulong();
411+
auto missIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), missN));
412+
mHExtensionPatternsGoodMissed->Fill(clusIdx, missIdx);
413+
}
414+
if (emptyPattern != clusPattern) {
415+
auto emptyN = emptyPattern.to_ulong();
416+
auto emptyIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), emptyN));
417+
mHExtensionPatternsGoodEmpty->Fill(clusIdx, emptyIdx);
418+
}
419+
}
420+
421+
if (mWithTree) {
422+
std::array<float, 3> xyz{(float)part.mcTrack.GetStartVertexCoordinatesX(), (float)part.mcTrack.GetStartVertexCoordinatesY(), (float)part.mcTrack.GetStartVertexCoordinatesZ()};
423+
std::array<float, 3> pxyz{(float)part.mcTrack.GetStartVertexMomentumX(), (float)part.mcTrack.GetStartVertexMomentumY(), (float)part.mcTrack.GetStartVertexMomentumZ()};
424+
auto pdg = O2DatabasePDG::Instance()->GetParticle(part.pdg);
425+
if (pdg == nullptr) {
426+
LOGP(fatal, "MC info not available");
427+
}
428+
auto mcTrk = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pdg->Charge() / 3.), true);
429+
(*mStream) << "tree"
430+
<< "trk=" << trk
431+
<< "mcTrk=" << mcTrk
432+
<< "\n";
433+
}
406434
}
407435
}
408436

@@ -421,30 +449,9 @@ void TrackExtensionStudy::endOfStream(EndOfStreamContext& ec)
421449
{
422450
LOGP(info, "Writing results to {}", mOutFileName);
423451
mStream->GetFile()->cd();
424-
425-
mHTrackCounts->Write();
426-
mHClustersCounts->Write();
427-
428-
mHLengthAny->Write();
429-
mHLengthGood->Write();
430-
mHLengthFake->Write();
431-
432-
mHChi2Any->Write();
433-
mHChi2Good->Write();
434-
mHChi2Fake->Write();
435-
436-
mHPtAny->Write();
437-
mHPtGood->Write();
438-
mHPtFake->Write();
439-
440-
mHExtensionAny->Write();
441-
mHExtensionGood->Write();
442-
mHExtensionFake->Write();
443-
444-
mHExtensionPatternsAny->Write();
445-
mHExtensionPatternsGood->Write();
446-
mHExtensionPatternsFake->Write();
447-
452+
for (const auto h : mHistograms) {
453+
h->Write();
454+
}
448455
mStream->Close();
449456
}
450457

0 commit comments

Comments
 (0)