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