Project
Loading...
Searching...
No Matches
TrackExtension.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
12#include "CommonUtils/TreeStreamRedirector.h"
13#include "DataFormatsGlobalTracking/RecoContainer.h"
14#include "DataFormatsITS/TrackITS.h"
15#include "DataFormatsITSMFT/CompCluster.h"
16#include "DetectorsBase/GRPGeomHelper.h"
17#include "DetectorsBase/Propagator.h"
18#include "Framework/ConfigParamRegistry.h"
19#include "Framework/Task.h"
20#include "ITSBase/GeometryTGeo.h"
21#include "ITSStudies/Helpers.h"
22#include "ITSStudies/TrackExtension.h"
23#include "SimulationDataFormat/MCEventHeader.h"
24#include "SimulationDataFormat/MCTrack.h"
25#include "Steer/MCKinematicsReader.h"
26#include "ReconstructionDataFormats/Vertex.h"
27#include "ReconstructionDataFormats/DCA.h"
28
29#include <bitset>
30
31#include "TFile.h"
32#include "TH1D.h"
33#include "TH2D.h"
34#include "TEfficiency.h"
35
36namespace o2::its::study
37{
38using namespace o2::framework;
39using namespace o2::globaltracking;
40
41using GTrackID = o2::dataformats::GlobalTrackID;
42using o2::steer::MCKinematicsReader;
43class TrackExtensionStudy : public Task
44{
45 struct ParticleInfo {
46 float eventX;
47 float eventY;
48 float eventZ;
49 int pdg;
50 float pt;
51 float eta;
52 float phi;
53 int mother;
54 int first;
55 float vx;
56 float vy;
57 float vz;
58 uint8_t clusters = 0u;
59 uint8_t fakeClusters = 0u;
60 uint8_t isReco = 0u;
61 uint8_t isFake = 0u;
62 bool isPrimary = false;
63 unsigned char storedStatus = 2;
64 int prodProcess;
65 o2::its::TrackITS track;
66 MCTrack mcTrack;
67 };
68
69 public:
70 TrackExtensionStudy(std::shared_ptr<DataRequest> dr,
71 mask_t src,
72 std::shared_ptr<o2::steer::MCKinematicsReader> kineReader,
73 std::shared_ptr<o2::base::GRPGeomRequest> gr) : mDataRequest(dr), mTracksSrc(src), mKineReader(kineReader), mGGCCDBRequest(gr)
74 {
75 LOGP(info, "Read MCKine reader with {} sources", mKineReader->getNSources());
76 }
77
78 ~TrackExtensionStudy() final = default;
79 void init(InitContext& /*ic*/) final;
80 void run(ProcessingContext& /*pc*/) final;
81 void endOfStream(EndOfStreamContext& /*ec*/) final;
82 void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final;
83 void process();
84
85 private:
86 static constexpr std::array<uint8_t, 9> mBitPatternsBefore{15, 30, 31, 60, 62, 63, 120, 124, 126};
87 static constexpr std::array<uint8_t, 16> mBitPatternsAfter{31, 47, 61, 62, 63, 79, 94, 95, 111, 121, 122, 123, 124, 125, 126, 127};
88 const std::bitset<7> mTopMask{"1110000"};
89 const std::bitset<7> mBotMask{"0000111"};
90
91 void updateTimeDependentParams(ProcessingContext& pc);
92 std::string mOutFileName = "TrackExtensionStudy.root";
93 std::shared_ptr<MCKinematicsReader> mKineReader;
94 GeometryTGeo* mGeometry{};
95
96 gsl::span<const o2::itsmft::ROFRecord> mTracksROFRecords;
97 gsl::span<const o2::its::TrackITS> mTracks;
98 gsl::span<const o2::MCCompLabel> mTracksMCLabels;
99 gsl::span<const o2::itsmft::CompClusterExt> mClusters;
100 gsl::span<const int> mInputITSidxs;
101 const o2::dataformats::MCLabelContainer* mClustersMCLCont{};
102
103 GTrackID::mask_t mTracksSrc{};
104 std::shared_ptr<DataRequest> mDataRequest;
105 std::vector<std::vector<std::vector<ParticleInfo>>> mParticleInfo; // src/event/track
106 unsigned short mMask = 0x7f;
107
108 std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
109 std::unique_ptr<utils::TreeStreamRedirector> mStream;
110 bool mWithTree{false};
111
112 std::unique_ptr<TH1D> mHTrackCounts;
113 std::unique_ptr<TH1D> mHLengthAny, mHLengthGood, mHLengthFake;
114 std::unique_ptr<TH1D> mHChi2Any, mHChi2Good, mHChi2Fake;
115 std::unique_ptr<TH1D> mHPtAny, mHPtGood, mHPtFake;
116 std::unique_ptr<TH1D> mHExtensionAny, mHExtensionGood, mHExtensionFake;
117 std::unique_ptr<TH2D> mHExtensionPatternsAny, mHExtensionPatternsGood, mHExtensionPatternsFake, mHExtensionPatternsGoodMissed, mHExtensionPatternsGoodEmpty;
118 std::unique_ptr<TH1D> mEExtensionNum, mEExtensionDen, mEExtensionPurityNum, mEExtensionPurityDen, mEExtensionFakeNum, mEExtensionFakeDen;
119 std::unique_ptr<TH1D> mEExtensionFakeBeforeNum, mEExtensionFakeAfterNum, mEExtensionFakeMixNum;
120 std::unique_ptr<TH1D> mEExtensionTopNum, mEExtensionTopPurityNum, mEExtensionTopFakeNum;
121 std::unique_ptr<TH1D> mEExtensionBotNum, mEExtensionBotPurityNum, mEExtensionBotFakeNum;
122 std::unique_ptr<TH1D> mEExtensionMixNum, mEExtensionMixPurityNum, mEExtensionMixFakeNum;
123 std::array<std::unique_ptr<TH1D>, mBitPatternsBefore.size()> mEExtensionPatternGoodNum, mEExtensionPatternFakeNum;
124 std::array<std::array<std::unique_ptr<TH1D>, mBitPatternsAfter.size()>, mBitPatternsBefore.size()> mEExtensionPatternIndGoodNum, mEExtensionPatternIndFakeNum;
125 // DCA
126 std::unique_ptr<TH2D> mDCAxyVsPtPionsNormal, mDCAxyVsPtPionsExtended;
127 std::unique_ptr<TH2D> mDCAzVsPtPionsNormal, mDCAzVsPtPionsExtended;
128
129 template <class T, typename... C, typename... F>
130 std::unique_ptr<T> createHistogram(C... n, F... b)
131 {
132 auto t = std::make_unique<T>(n..., b...);
133 mHistograms.push_back(static_cast<TH1*>(t.get()));
134 return std::move(t);
135 }
136 std::vector<TH1*> mHistograms;
137};
138
139void TrackExtensionStudy::init(InitContext& ic)
140{
141 o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest);
142 mWithTree = ic.options().get<bool>("with-tree");
143
144 constexpr size_t effHistBins = 40;
145 constexpr float effPtCutLow = 0.01;
146 constexpr float effPtCutHigh = 10.;
147 auto xbins = helpers::makeLogBinning(effHistBins, effPtCutLow, effPtCutHigh);
148
149 // Track Counting
150 mHTrackCounts = createHistogram<TH1D>("hTrackCounts", "Track Stats", 10, 0, 10);
151 mHTrackCounts->GetXaxis()->SetBinLabel(1, "Total Tracks");
152 mHTrackCounts->GetXaxis()->SetBinLabel(2, "Normal ANY Tracks");
153 mHTrackCounts->GetXaxis()->SetBinLabel(3, "Normal GOOD Tracks");
154 mHTrackCounts->GetXaxis()->SetBinLabel(4, "Normal FAKE Tracks");
155 mHTrackCounts->GetXaxis()->SetBinLabel(5, "Extended ANY Tracks");
156 mHTrackCounts->GetXaxis()->SetBinLabel(6, "Extended GOOD Tracks");
157 mHTrackCounts->GetXaxis()->SetBinLabel(7, "Extended FAKE Tracks");
158 mHTrackCounts->GetXaxis()->SetBinLabel(8, "Extended FAKE BEFORE Tracks");
159 mHTrackCounts->GetXaxis()->SetBinLabel(9, "Extended FAKE AFTER Tracks");
160 mHTrackCounts->GetXaxis()->SetBinLabel(10, "Extended FAKE BEFORE&AFTER Tracks");
161
162 // Length
163 mHLengthAny = createHistogram<TH1D>("hLengthAny", "Extended Tracks Length (ANY);NCluster;Entries", 5, 3, 8);
164 mHLengthGood = createHistogram<TH1D>("hLengthGood", "Extended Tracks Length (GOOD);NCluster;Entries", 5, 3, 8);
165 mHLengthFake = createHistogram<TH1D>("hLengthFake", "Extended Tracks Length (FAKE);NCluster;Entries", 5, 3, 8);
166
167 // Chi2
168 mHChi2Any = createHistogram<TH1D>("hChi2Any", "Extended Tracks Length (ANY);#chi^{2};Entries", 50, 0, 100);
169 mHChi2Good = createHistogram<TH1D>("hChi2Good", "Extended Tracks Length (GOOD);#chi^{2};Entries", 50, 0, 100);
170 mHChi2Fake = createHistogram<TH1D>("hChi2Fake", "Extended Tracks Length (FAKE);#chi^{2};Entries", 50, 0, 100);
171
172 // Pt
173 mHPtAny = createHistogram<TH1D>("hPtAny", "Extended Tracks Length (ANY);#it{p}_{T};Entries", effHistBins, xbins.data());
174 mHPtGood = createHistogram<TH1D>("hPtGood", "Extended Tracks Length (GOOD);#it{p}_{T};Entries", effHistBins, xbins.data());
175 mHPtFake = createHistogram<TH1D>("hPtFake", "Extended Tracks Length (FAKE);#it{p}_{T};Entries", effHistBins, xbins.data());
176
177 // Length
178 mHExtensionAny = createHistogram<TH1D>("hExtensionAny", "Extended Tracks Length (ANY);Extended Layer;Entries", 7, 0, 7);
179 mHExtensionGood = createHistogram<TH1D>("hExtensionGood", "Extended Tracks Length (GOOD);Extended Layer;Entries", 7, 0, 7);
180 mHExtensionFake = createHistogram<TH1D>("hExtensionFake", "Extended Tracks Length (FAKE);Extended Layer;Entries", 7, 0, 7);
181
182 // Patterns
183 auto makePatternAxisLabels = [&](TH1* h, bool xBefore = true) {
184 for (int i{1}; i <= h->GetXaxis()->GetNbins(); ++i) {
185 if (xBefore) {
186 h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsBefore[i - 1]).c_str());
187 } else {
188 h->GetXaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str());
189 }
190 }
191 for (int i{1}; i <= h->GetYaxis()->GetNbins(); ++i) {
192 h->GetYaxis()->SetBinLabel(i, fmt::format("{:07b}", mBitPatternsAfter[i - 1]).c_str());
193 }
194 };
195 mHExtensionPatternsAny = createHistogram<TH2D>("hExtensionPatternsAny", "Extended Tracks Pattern (ANY);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
196 makePatternAxisLabels(mHExtensionPatternsAny.get());
197 mHExtensionPatternsGood = createHistogram<TH2D>("hExtensionPatternsGood", "Extended Tracks Pattern (GOOD);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
198 makePatternAxisLabels(mHExtensionPatternsGood.get());
199 mHExtensionPatternsFake = createHistogram<TH2D>("hExtensionPatternsFake", "Extended Tracks Pattern (FAKE);Before;After;Entries", mBitPatternsBefore.size(), 0, mBitPatternsBefore.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
200 makePatternAxisLabels(mHExtensionPatternsFake.get());
201 mHExtensionPatternsGoodMissed = createHistogram<TH2D>("hExtensionPatternsGoodMissed", "Extended Tracks Pattern (GOOD) Missed Clusters;After;Missed;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
202 makePatternAxisLabels(mHExtensionPatternsGoodMissed.get(), false);
203 mHExtensionPatternsGoodEmpty = createHistogram<TH2D>("hExtensionPatternsGoodEmpty", "Extended Tracks Pattern (GOOD) Empty Clusters;Before;After;Entries", mBitPatternsAfter.size(), 0, mBitPatternsAfter.size(), mBitPatternsAfter.size(), 0, mBitPatternsAfter.size());
204 makePatternAxisLabels(mHExtensionPatternsGoodEmpty.get(), false);
205
207 mEExtensionNum = createHistogram<TH1D>("hExtensionNum", "Extension Numerator", effHistBins, xbins.data());
208 mEExtensionDen = createHistogram<TH1D>("hExtensionDen", "Extension Dennominator", effHistBins, xbins.data());
209 // Purity
210 mEExtensionPurityNum = createHistogram<TH1D>("hExtensionPurityNum", "Extension Purity Numerator", effHistBins, xbins.data());
211 mEExtensionPurityDen = createHistogram<TH1D>("hExtensionPurityDen", "Extension Purity Denominator", effHistBins, xbins.data());
212 // Fake
213 mEExtensionFakeNum = createHistogram<TH1D>("hExtensionFakeNum", "Extension Fake Numerator", effHistBins, xbins.data());
214 mEExtensionFakeDen = createHistogram<TH1D>("hExtensionFakeDen", "Extension Fake Denominator", effHistBins, xbins.data());
215 mEExtensionFakeBeforeNum = createHistogram<TH1D>("hExtensionFakeBeforeNum", "Extension Fake Before Numerator", effHistBins, xbins.data());
216 mEExtensionFakeAfterNum = createHistogram<TH1D>("hExtensionFakeAfterNum", "Extension Fake After Numerator", effHistBins, xbins.data());
217 mEExtensionFakeMixNum = createHistogram<TH1D>("hExtensionFakeMixNum", "Extension Fake Mix Numerator", effHistBins, xbins.data());
218 // Top
219 mEExtensionTopNum = createHistogram<TH1D>("hExtensionTopNum", "Extension Top Numerator", effHistBins, xbins.data());
220 mEExtensionTopPurityNum = createHistogram<TH1D>("hExtensionTopPurityNum", "Extension Top Purity Numerator", effHistBins, xbins.data());
221 mEExtensionTopFakeNum = createHistogram<TH1D>("hExtensionTopFakeNum", "Extension Top Fake Numerator", effHistBins, xbins.data());
222 mEExtensionBotNum = createHistogram<TH1D>("hExtensionBotNum", "Extension Bot Numerator", effHistBins, xbins.data());
223 mEExtensionBotPurityNum = createHistogram<TH1D>("hExtensionBotPurityNum", "Extension Bot Purity Numerator", effHistBins, xbins.data());
224 mEExtensionBotFakeNum = createHistogram<TH1D>("hExtensionBotFakeNum", "Extension Bot Fake Numerator", effHistBins, xbins.data());
225 mEExtensionMixNum = createHistogram<TH1D>("hExtensionMixNum", "Extension Mix Numerator", effHistBins, xbins.data());
226 mEExtensionMixPurityNum = createHistogram<TH1D>("hExtensionMixPurityNum", "Extension Mix Purity Numerator", effHistBins, xbins.data());
227 mEExtensionMixFakeNum = createHistogram<TH1D>("hExtensionMixFakeNum", "Extension Mix Fake Numerator", effHistBins, xbins.data());
228 // Patterns
229 for (int i{0}; i < mBitPatternsBefore.size(); ++i) {
230 mEExtensionPatternGoodNum[i] = createHistogram<TH1D>(fmt::format("hExtensionPatternGood_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b}", mBitPatternsBefore[i]).c_str(), effHistBins, xbins.data());
231 mEExtensionPatternFakeNum[i] = createHistogram<TH1D>(fmt::format("hExtensionPatternFake_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b}", mBitPatternsBefore[i]).c_str(), effHistBins, xbins.data());
232 for (int j{0}; j < mBitPatternsAfter.size(); ++j) {
233 mEExtensionPatternIndGoodNum[i][j] = createHistogram<TH1D>(fmt::format("hExtensionPatternGood_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} -> {:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), effHistBins, xbins.data());
234 mEExtensionPatternIndFakeNum[i][j] = createHistogram<TH1D>(fmt::format("hExtensionPatternFake_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} -> {:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), effHistBins, xbins.data());
235 }
236 }
237
239 mDCAxyVsPtPionsNormal = createHistogram<TH2D>("hDCAxyVsPtResNormal", "DCA_{#it{xy}} NORMAL Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500);
240 mDCAxyVsPtPionsExtended = createHistogram<TH2D>("hDCAxyVsPtResExtended", "DCA_{#it{xy}} EXTENDED Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500);
241 mDCAzVsPtPionsNormal = createHistogram<TH2D>("hDCAzVsPtResNormal", "DCA_{#it{z}} NORMAL Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500);
242 mDCAzVsPtPionsExtended = createHistogram<TH2D>("hDCAzVsPtResExtended", "DCA_{#it{z}} EXTENDED Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", effHistBins, xbins.data(), 1000, -500, 500);
243
244 mStream = std::make_unique<utils::TreeStreamRedirector>(mOutFileName.c_str(), "RECREATE");
245}
246
247void TrackExtensionStudy::run(ProcessingContext& pc)
248{
249 o2::globaltracking::RecoContainer recoData;
250 recoData.collectData(pc, *mDataRequest);
251 updateTimeDependentParams(pc);
252
253 mTracksROFRecords = recoData.getITSTracksROFRecords();
254 mTracks = recoData.getITSTracks();
255 mTracksMCLabels = recoData.getITSTracksMCLabels();
256 mClusters = recoData.getITSClusters();
257 mClustersMCLCont = recoData.getITSClustersMCLabels();
258 mInputITSidxs = recoData.getITSTracksClusterRefs();
259
260 LOGP(info, "** Found in {} rofs:\n\t- {} clusters with {} labels\n\t- {} tracks with {} labels",
261 mTracksROFRecords.size(), mClusters.size(), mClustersMCLCont->getIndexedSize(), mTracks.size(), mTracksMCLabels.size());
262 LOGP(info, "** Found {} sources from kinematic files", mKineReader->getNSources());
263
264 process();
265}
266
267void TrackExtensionStudy::process()
268{
269 LOGP(info, "** Filling particle table ... ");
270 mParticleInfo.resize(mKineReader->getNSources()); // sources
271 for (int iSource{0}; iSource < mKineReader->getNSources(); ++iSource) {
272 mParticleInfo[iSource].resize(mKineReader->getNEvents(iSource)); // events
273 for (int iEvent{0}; iEvent < mKineReader->getNEvents(iSource); ++iEvent) {
274 const auto& mcEvent = mKineReader->getMCEventHeader(iSource, iEvent);
275 mParticleInfo[iSource][iEvent].resize(mKineReader->getTracks(iSource, iEvent).size()); // tracks
276 for (auto iPart{0}; iPart < mKineReader->getTracks(iEvent).size(); ++iPart) {
277 const auto& part = mKineReader->getTracks(iSource, iEvent)[iPart];
278 mParticleInfo[iSource][iEvent][iPart].eventX = mcEvent.GetX();
279 mParticleInfo[iSource][iEvent][iPart].eventY = mcEvent.GetY();
280 mParticleInfo[iSource][iEvent][iPart].eventZ = mcEvent.GetZ();
281 mParticleInfo[iSource][iEvent][iPart].pdg = part.GetPdgCode();
282 mParticleInfo[iSource][iEvent][iPart].pt = part.GetPt();
283 mParticleInfo[iSource][iEvent][iPart].phi = part.GetPhi();
284 mParticleInfo[iSource][iEvent][iPart].eta = part.GetEta();
285 mParticleInfo[iSource][iEvent][iPart].vx = part.Vx();
286 mParticleInfo[iSource][iEvent][iPart].vy = part.Vy();
287 mParticleInfo[iSource][iEvent][iPart].vz = part.Vz();
288 mParticleInfo[iSource][iEvent][iPart].isPrimary = part.isPrimary();
289 mParticleInfo[iSource][iEvent][iPart].mother = part.getMotherTrackId();
290 mParticleInfo[iSource][iEvent][iPart].prodProcess = part.getProcess();
291 }
292 }
293 }
294 LOGP(info, "** Creating particle/clusters correspondance ... ");
295 for (auto iSource{0}; iSource < mParticleInfo.size(); ++iSource) {
296 for (auto iCluster{0}; iCluster < mClusters.size(); ++iCluster) {
297 auto labs = mClustersMCLCont->getLabels(iCluster); // ideally I can have more than one label per cluster
298 for (auto& lab : labs) {
299 if (!lab.isValid()) {
300 continue; // We want to skip channels related to noise, e.g. sID = 99: QED
301 }
302 int trackID, evID, srcID;
303 bool fake;
304 lab.get(trackID, evID, srcID, fake);
305 auto& cluster = mClusters[iCluster];
306 auto layer = mGeometry->getLayer(cluster.getSensorID());
307 mParticleInfo[srcID][evID][trackID].clusters |= (1 << layer);
308 if (fake) {
309 mParticleInfo[srcID][evID][trackID].fakeClusters |= (1 << layer);
310 }
311 }
312 }
313 }
314
315 LOGP(info, "** Analysing tracks ... ");
316 int unaccounted{0}, good{0}, fakes{0}, extended{0};
317 for (auto iTrack{0}; iTrack < mTracks.size(); ++iTrack) {
318 const auto& lab = mTracksMCLabels[iTrack];
319 if (!lab.isValid()) {
320 unaccounted++;
321 continue;
322 }
323 int trackID, evID, srcID;
324 bool fake;
325 lab.get(trackID, evID, srcID, fake);
326
327 if (srcID == 99) { // skip QED
328 unaccounted++;
329 continue;
330 }
331
332 for (int iLayer{0}; iLayer < 7; ++iLayer) {
333 if (mTracks[iTrack].isExtendedOnLayer(iLayer)) {
334 ++extended;
335 break;
336 }
337 }
338
339 mParticleInfo[srcID][evID][trackID].isReco += !fake;
340 mParticleInfo[srcID][evID][trackID].isFake += fake;
341 if (mTracks[iTrack].isBetter(mParticleInfo[srcID][evID][trackID].track, 1.e9)) {
342 mParticleInfo[srcID][evID][trackID].storedStatus = fake;
343 mParticleInfo[srcID][evID][trackID].track = mTracks[iTrack];
344 mParticleInfo[srcID][evID][trackID].mcTrack = *mKineReader->getTrack(lab);
345 }
346 fakes += fake;
347 good += !fake;
348 }
349 LOGP(info, "** Some statistics:");
350 LOGP(info, "\t- Total number of tracks: {}", mTracks.size());
351 LOGP(info, "\t- Total number of tracks not corresponding to particles: {} ({:.2f} %)", unaccounted, unaccounted * 100. / mTracks.size());
352 LOGP(info, "\t- Total number of fakes: {} ({:.2f} %)", fakes, fakes * 100. / mTracks.size());
353 LOGP(info, "\t- Total number of good: {} ({:.2f} %)", good, good * 100. / mTracks.size());
354 LOGP(info, "\t- Total number of extensions: {} ({:.2f} %)", extended, extended * 100. / mTracks.size());
355
356 o2::dataformats::VertexBase collision;
357 o2::dataformats::DCA impactParameter;
358 LOGP(info, "** Filling histograms ... ");
359 for (auto iTrack{0}; iTrack < mTracks.size(); ++iTrack) {
360 auto& lab = mTracksMCLabels[iTrack];
361 if (!lab.isValid()) {
362 unaccounted++;
363 continue;
364 }
365 int trackID, evID, srcID;
366 bool fake;
367 lab.get(trackID, evID, srcID, fake);
368 const auto& part = mParticleInfo[srcID][evID][trackID];
369 if (!part.isPrimary) {
370 continue;
371 }
372 const auto& trk = part.track;
373 bool isGood = part.isReco && !part.isFake;
374 mHTrackCounts->Fill(0);
375
376 std::bitset<7> extPattern{0};
377 for (int iLayer{0}; iLayer < 7; ++iLayer) {
378 if (trk.isExtendedOnLayer(iLayer)) {
379 extPattern.set(iLayer);
380 }
381 }
382
383 // Tree
384 while (mWithTree) {
385 constexpr float refRadius{70.f};
386 constexpr float maxSnp{0.9f};
387 auto cTrk = trk;
388 if (!o2::base::Propagator::Instance()->PropagateToXBxByBz(cTrk, refRadius, maxSnp, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrTGeo)) {
389 break;
390 }
391 std::array<float, 3> xyz{(float)part.mcTrack.GetStartVertexCoordinatesX(), (float)part.mcTrack.GetStartVertexCoordinatesY(), (float)part.mcTrack.GetStartVertexCoordinatesZ()};
392 std::array<float, 3> pxyz{(float)part.mcTrack.GetStartVertexMomentumX(), (float)part.mcTrack.GetStartVertexMomentumY(), (float)part.mcTrack.GetStartVertexMomentumZ()};
393 auto pdg = O2DatabasePDG::Instance()->GetParticle(part.pdg);
394 if (pdg == nullptr) {
395 LOGP(error, "MC info not available");
396 break;
397 }
398 auto mcTrk = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pdg->Charge() / 3.), true);
399 if (!mcTrk.rotate(cTrk.getAlpha()) || !o2::base::Propagator::Instance()->PropagateToXBxByBz(mcTrk, refRadius, maxSnp, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrTGeo)) {
400 break;
401 }
402 (*mStream) << "tree"
403 << "trk=" << cTrk
404 << "mcTrk=" << mcTrk
405 << "isGood=" << isGood
406 << "isExtended=" << extPattern.any()
407 << "\n";
408 break;
409 }
410
411 // impact parameter
412 while (isGood && std::abs(part.pdg) == 211) {
413 auto trkC = part.track;
414 collision.setXYZ(part.eventX, part.eventY, part.eventZ);
415 if (!o2::base::Propagator::Instance()->propagateToDCA(collision, trkC, o2::base::Propagator::Instance()->getNominalBz(), 2.0, o2::base::Propagator::MatCorrType::USEMatCorrTGeo, &impactParameter)) {
416 break;
417 }
418
419 auto dcaXY = impactParameter.getY() * 1e4;
420 auto dcaZ = impactParameter.getZ() * 1e4;
421 if (!extPattern.any()) {
422 mDCAxyVsPtPionsNormal->Fill(part.pt, dcaXY);
423 mDCAzVsPtPionsNormal->Fill(part.pt, dcaZ);
424 } else {
425 mDCAxyVsPtPionsExtended->Fill(part.pt, dcaXY);
426 mDCAzVsPtPionsExtended->Fill(part.pt, dcaZ);
427 }
428 break;
429 }
430
431 mEExtensionDen->Fill(trk.getPt());
432
433 if (!extPattern.any()) {
434 mHTrackCounts->Fill(1);
435 if (part.isReco || !part.isFake) {
436 mHTrackCounts->Fill(2);
437 } else {
438 mHTrackCounts->Fill(3);
439 }
440 continue;
441 }
442
443 mHTrackCounts->Fill(4);
444 mHLengthAny->Fill(trk.getNClusters());
445 mHChi2Any->Fill(trk.getChi2());
446 mHPtAny->Fill(trk.getPt());
447 mEExtensionNum->Fill(trk.getPt());
448 mEExtensionPurityDen->Fill(trk.getPt());
449 mEExtensionFakeDen->Fill(trk.getPt());
450 if (isGood) {
451 mHTrackCounts->Fill(5);
452 mHLengthGood->Fill(trk.getNClusters());
453 mHChi2Good->Fill(trk.getChi2());
454 mHPtGood->Fill(trk.getPt());
455 mEExtensionPurityNum->Fill(trk.getPt());
456 } else {
457 mHTrackCounts->Fill(6);
458 mHLengthFake->Fill(trk.getNClusters());
459 mHChi2Fake->Fill(trk.getChi2());
460 mHPtFake->Fill(trk.getPt());
461 mEExtensionFakeNum->Fill(trk.getPt());
462 }
463
464 std::bitset<7> clusPattern{static_cast<uint8_t>(trk.getPattern())};
465 for (int iLayer{0}; iLayer < 7; ++iLayer) {
466 if (extPattern.test(iLayer)) {
467 extPattern.set(iLayer);
468 mHExtensionAny->Fill(iLayer);
469 if (isGood) {
470 mHExtensionGood->Fill(iLayer);
471 } else {
472 mHExtensionFake->Fill(iLayer);
473 }
474 }
475 }
476 std::bitset<7> oldPattern{clusPattern & ~extPattern}, holePattern{clusPattern};
477 holePattern.flip();
478 auto clusN = clusPattern.to_ulong();
479 auto clusIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), clusN));
480 auto oldN = oldPattern.to_ulong();
481 auto oldIdx = std::distance(std::begin(mBitPatternsBefore), std::find(std::begin(mBitPatternsBefore), std::end(mBitPatternsBefore), oldN));
482 mHExtensionPatternsAny->Fill(oldIdx, clusIdx);
483 if (isGood) {
484 mHExtensionPatternsGood->Fill(oldIdx, clusIdx);
485 mEExtensionPatternGoodNum[oldIdx]->Fill(trk.getPt());
486 mEExtensionPatternIndGoodNum[oldIdx][clusIdx]->Fill(trk.getPt());
487 } else {
488 mHExtensionPatternsFake->Fill(oldIdx, clusIdx);
489 mEExtensionPatternFakeNum[oldIdx]->Fill(trk.getPt());
490 mEExtensionPatternIndFakeNum[oldIdx][clusIdx]->Fill(trk.getPt());
491 }
492
493 // old pattern
494 bool oldFake{false}, newFake{false};
495 for (int iLayer{0}; iLayer < 7; ++iLayer) {
496 if (trk.isFakeOnLayer(iLayer)) {
497 if (oldPattern.test(iLayer)) {
498 oldFake = true;
499 } else if (extPattern.test(iLayer)) {
500 newFake = true;
501 }
502 }
503 }
504 if (oldFake && newFake) {
505 mHTrackCounts->Fill(9);
506 mEExtensionFakeMixNum->Fill(trk.getPt());
507 } else if (oldFake) {
508 mHTrackCounts->Fill(7);
509 mEExtensionFakeBeforeNum->Fill(trk.getPt());
510 } else if (newFake) {
511 mHTrackCounts->Fill(8);
512 mEExtensionFakeAfterNum->Fill(trk.getPt());
513 }
514
515 // Check if we missed some clusters
516 if (isGood && holePattern.any()) {
517 auto missPattern{clusPattern}, emptyPattern{clusPattern};
518 for (int iLayer{0}; iLayer < 7; ++iLayer) {
519 if (!holePattern.test(iLayer)) {
520 continue;
521 }
522
523 // Check if there was actually a cluster that we missed
524 if ((part.clusters & (1 << iLayer)) != 0) {
525 missPattern.set(iLayer);
526 } else {
527 emptyPattern.set(iLayer);
528 }
529 }
530
531 if (missPattern != clusPattern) {
532 auto missN = missPattern.to_ulong();
533 auto missIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), missN));
534 mHExtensionPatternsGoodMissed->Fill(clusIdx, missIdx);
535 }
536 if (emptyPattern != clusPattern) {
537 auto emptyN = emptyPattern.to_ulong();
538 auto emptyIdx = std::distance(std::begin(mBitPatternsAfter), std::find(std::begin(mBitPatternsAfter), std::end(mBitPatternsAfter), emptyN));
539 mHExtensionPatternsGoodEmpty->Fill(clusIdx, emptyIdx);
540 }
541 }
542
543 // Top/Bot/Mixed Extension
544 bool isTop = (extPattern & mTopMask).any();
545 bool isBot = (extPattern & mBotMask).any();
546 if (isTop && isBot) {
547 mEExtensionMixNum->Fill(trk.getPt());
548 if (isGood) {
549 mEExtensionMixPurityNum->Fill(trk.getPt());
550 } else {
551 mEExtensionMixFakeNum->Fill(trk.getPt());
552 }
553 } else if (isTop) {
554 mEExtensionTopNum->Fill(trk.getPt());
555 if (isGood) {
556 mEExtensionTopPurityNum->Fill(trk.getPt());
557 } else {
558 mEExtensionTopFakeNum->Fill(trk.getPt());
559 }
560 } else {
561 mEExtensionBotNum->Fill(trk.getPt());
562 if (isGood) {
563 mEExtensionBotPurityNum->Fill(trk.getPt());
564 } else {
565 mEExtensionBotFakeNum->Fill(trk.getPt());
566 }
567 }
568 }
569}
570
571void TrackExtensionStudy::updateTimeDependentParams(ProcessingContext& pc)
572{
573 static bool initOnceDone = false;
574 o2::base::GRPGeomHelper::instance().checkUpdates(pc);
575 if (!initOnceDone) { // this params need to be queried only once
576 initOnceDone = true;
577 mGeometry = GeometryTGeo::Instance();
578 mGeometry->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot, o2::math_utils::TransformType::T2G));
579 }
580}
581
582void TrackExtensionStudy::endOfStream(EndOfStreamContext& ec)
583{
584 LOGP(info, "Writing results to {}", mOutFileName);
585 mStream->GetFile()->cd();
586 for (const auto h : mHistograms) {
587 h->Write();
588 }
589
590 LOGP(info, "Calculating efficiencies");
591 auto makeEff = [](auto num, auto den, const char* name, const char* title) {
592 auto e = std::make_unique<TEfficiency>(*num, *den);
593 e->SetName(name);
594 e->SetTitle(title);
595 e->Write();
596 };
597 makeEff(mEExtensionNum.get(), mEExtensionDen.get(), "eExtension", "Track Extension EXT TRK/ALL");
598 makeEff(mEExtensionPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionPurity", "Track Extension Purity GOOD/EXT TRK");
599 makeEff(mEExtensionFakeNum.get(), mEExtensionFakeDen.get(), "eExtensionFake", "Track Extension Fake FAKE/EXT TRK");
600 makeEff(mEExtensionFakeBeforeNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeBefore", "Track Extension Fake FAKE BEF/FAKE EXT TRK");
601 makeEff(mEExtensionFakeAfterNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeAfter", "Track Extension Fake FAKE AFT/FAKE EXT TRK");
602 makeEff(mEExtensionFakeMixNum.get(), mEExtensionFakeNum.get(), "eExtensionFakeMix", "Track Extension Fake FAKE MIX/FAKE EXT TRK");
603 makeEff(mEExtensionTopNum.get(), mEExtensionDen.get(), "eExtensionTop", "Track Extension Top");
604 makeEff(mEExtensionTopPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionTopPurity", "Track Extension Purity GOOD TOP/EXT TRK");
605 makeEff(mEExtensionTopFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionTopFake", "Track Extension FAKE TOP/EXT FAKE TRK");
606 makeEff(mEExtensionBotNum.get(), mEExtensionDen.get(), "eExtensionBot", "Track Extension Bot");
607 makeEff(mEExtensionBotPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionBotPurity", "Track Extension Purity GOOD BOT/EXT TRK");
608 makeEff(mEExtensionBotFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionBotFake", "Track Extension FAKE BOT/EXT FAKE TRK");
609 makeEff(mEExtensionMixNum.get(), mEExtensionDen.get(), "eExtensionMix", "Track Extension Mix");
610 makeEff(mEExtensionMixPurityNum.get(), mEExtensionPurityDen.get(), "eExtensionMixPurity", "Track Extension Purity GOOD MIX/EXT TRK");
611 makeEff(mEExtensionMixFakeNum.get(), mEExtensionFakeNum.get(), "eExtensionMixFake", "Track Extension FAKE MIX/EXT FAKE TRK");
612 for (int i{0}; i < mBitPatternsBefore.size(); ++i) {
613 makeEff(mEExtensionPatternGoodNum[i].get(), mEExtensionPurityNum.get(), fmt::format("eExtensionPatternGood_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} GOOD EXT TRK/EXT TRK", mBitPatternsBefore[i]).c_str());
614 makeEff(mEExtensionPatternFakeNum[i].get(), mEExtensionFakeNum.get(), fmt::format("eExtensionPatternFake_{:07b}", mBitPatternsBefore[i]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} FAKE EXT TRK/EXT TRK", mBitPatternsBefore[i]).c_str());
615 for (int j{0}; j < mBitPatternsAfter.size(); ++j) {
616 makeEff(mEExtensionPatternIndGoodNum[i][j].get(), mEExtensionPatternGoodNum[i].get(), fmt::format("eExtensionPatternGood_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (GOOD) {:07b} -> {:07b} GOOD EXT TRK/EXT TRK", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str());
617 makeEff(mEExtensionPatternIndFakeNum[i][j].get(), mEExtensionPatternFakeNum[i].get(), fmt::format("eExtensionPatternFake_{:07b}_{:07b}", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str(), fmt::format("Extended Tracks Pattern (FAKE) {:07b} -> {:07b} FAKE EXT TRK/EXT TRK", mBitPatternsBefore[i], mBitPatternsAfter[j]).c_str());
618 }
619 }
620
621 mStream->Close();
622}
623
624void TrackExtensionStudy::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
625{
626 if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
627 return;
628 }
629}
630
631DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, std::shared_ptr<o2::steer::MCKinematicsReader> kineReader)
632{
633 std::vector<OutputSpec> outputs;
634 auto dataRequest = std::make_shared<DataRequest>();
635 dataRequest->requestTracks(srcTracksMask, true);
636 dataRequest->requestClusters(srcClustersMask, true);
637
638 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(false, // orbitResetTime
639 true, // GRPECS=true
640 false, // GRPLHCIF
641 true, // GRPMagField
642 true, // askMatLUT
643 o2::base::GRPGeomRequest::Aligned, // geometry
644 dataRequest->inputs,
645 true);
646
647 return DataProcessorSpec{
648 "its-study-track-extension",
649 dataRequest->inputs,
650 outputs,
651 AlgorithmSpec{adaptFromTask<TrackExtensionStudy>(dataRequest, srcTracksMask, kineReader, ggRequest)},
652 Options{{"with-tree", o2::framework::VariantType::Bool, false, {"Produce in addition a tree"}}}};
653}
654
655} // namespace o2::its::study
CompCluster.h
Definition of the ITSMFT compact cluster.
RecoContainer.h
Wrapper container for different reconstructed object types.
i
int32_t i
Definition GPUCommonAlgorithm.h:443
GRPGeomHelper.h
Helper for geometry and GRP related CCDB requests.
GeometryTGeo.h
Definition of the GeometryTGeo class.
MCTrack.h
Definition of the MCTrack class.
j
uint32_t j
Definition RawData.h:0
TrackITS.h
Definition of the ITS track.
num
double num
Definition bench_Clusterizer.cxx:122
h
Class for time synchronization of RawReader instances.
o2::O2DatabasePDG::Instance
static TDatabasePDG * Instance()
Definition O2DatabasePDG.h:35
o2::base::GRPGeomHelper::checkUpdates
void checkUpdates(o2::framework::ProcessingContext &pc)
Definition GRPGeomHelper.cxx:178
o2::base::GRPGeomHelper::instance
static GRPGeomHelper & instance()
Definition GRPGeomHelper.h:120
o2::base::GRPGeomHelper::setRequest
void setRequest(std::shared_ptr< GRPGeomRequest > req)
Definition GRPGeomHelper.cxx:101
o2::base::PropagatorImpl< float >::Instance
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
o2::dataformats::MCTruthContainer::getLabels
gsl::span< TruthElement > getLabels(uint32_t dataindex)
Definition MCTruthContainer.h:153
o2::framework::InitContext::options
ConfigParamRegistry const & options()
Definition InitContext.h:33
o2::its::GeometryTGeo::getLayer
int getLayer(int index) const
Get chip layer, from 0.
Definition GeometryTGeo.cxx:154
o2::its::GeometryTGeo::Instance
static GeometryTGeo * Instance()
Definition GeometryTGeo.h:55
o2::its::GeometryTGeo::fillMatrixCache
void fillMatrixCache(int mask) override
Definition GeometryTGeo.cxx:525
o2::its::study::TrackExtensionStudy::run
void run(ProcessingContext &) final
Definition TrackExtension.cxx:247
o2::its::study::TrackExtensionStudy::finaliseCCDB
void finaliseCCDB(ConcreteDataMatcher &matcher, void *obj) final
Definition TrackExtension.cxx:624
o2::its::study::TrackExtensionStudy::endOfStream
void endOfStream(EndOfStreamContext &) final
This is invoked whenever we have an EndOfStream event.
Definition TrackExtension.cxx:582
o2::its::study::TrackExtensionStudy::TrackExtensionStudy
TrackExtensionStudy(std::shared_ptr< DataRequest > dr, mask_t src, std::shared_ptr< o2::steer::MCKinematicsReader > kineReader, std::shared_ptr< o2::base::GRPGeomRequest > gr)
Definition TrackExtension.cxx:70
n
GLdouble n
Definition glcorearb.h:1982
src
GLenum src
Definition glcorearb.h:1767
size
GLsizeiptr size
Definition glcorearb.h:659
array
GLenum array
Definition glcorearb.h:4274
name
GLuint const GLchar * name
Definition glcorearb.h:781
b
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
layer
GLenum GLuint GLint GLint layer
Definition glcorearb.h:1310
o2::framework
Defining PrimaryVertex explicitly as messageable.
Definition TFIDInfo.h:20
o2::framework::Options
std::vector< ConfigParamSpec > Options
Definition DataProcessorSpec.h:31
o2::header::get
auto get(const std::byte *buffer, size_t=0)
Definition DataHeader.h:454
o2::its::study::helpers::makeLogBinning
std::vector< double > makeLogBinning(const int nbins, const double min, const double max)
get a vector containing binning info for constant sized bins on a log axis
Definition Helpers.cxx:21
o2::its::study::dcaZ
float dcaZ[3]
Definition Efficiency.h:49
o2::its::study::getTrackExtensionStudy
o2::framework::DataProcessorSpec getTrackExtensionStudy(mask_t srcTracksMask, mask_t srcClustersMask, std::shared_ptr< o2::steer::MCKinematicsReader > kineReader)
Definition TrackExtension.cxx:631
o2::its::study::dcaXY
float dcaXY[3]
Definition Efficiency.h:48
o2::math_utils::detail::T
T
Definition SMatrixGPU.h:747
o2::tpc::C
@ C
Definition Defs.h:36
o2::track::TrackPar
TrackParF TrackPar
Definition Track.h:29
std
Defining DataPointCompositeObject explicitly as copiable.
Definition CcdbObjectInfo.h:121
o2::globaltracking::RecoContainer::collectData
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
Definition RecoContainer.cxx:609
o2::math_utils::TransformType::T2L
static constexpr int T2L
Definition Cartesian.h:55
o2::math_utils::TransformType::T2GRot
static constexpr int T2GRot
Definition Cartesian.h:57
o2::math_utils::TransformType::T2G
static constexpr int T2G
Definition Cartesian.h:56
clusters
std::vector< Cluster > clusters
Definition test_ctf_io_cpv.cxx:41