Project
Loading...
Searching...
No Matches
GPUQA.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
14
15#define QA_DEBUG 0
16#define QA_TIMING 0
17
18#include "Rtypes.h" // Include ROOT header first, to use ROOT and disable replacements
19
20#include "TH1F.h"
21#include "TH2F.h"
22#include "TH1D.h"
23#include "TGraphAsymmErrors.h"
24#include "TCanvas.h"
25#include "TPad.h"
26#include "TLegend.h"
27#include "TColor.h"
28#include "TPaveText.h"
29#include "TF1.h"
30#include "TFile.h"
31#include "TTree.h"
32#include "TStyle.h"
33#include "TLatex.h"
34#include "TObjArray.h"
35#include <sys/stat.h>
36
37#include "GPUQA.h"
38#include "GPUTPCDef.h"
39#include "GPUTPCTrackingData.h"
40#include "GPUChainTracking.h"
41#include "GPUChainTrackingGetters.inc"
42#include "GPUTPCTrack.h"
43#include "GPUTPCTracker.h"
44#include "GPUTPCGMMergedTrack.h"
45#include "GPUTPCGMPropagator.h"
47#include "GPUTPCMCInfo.h"
48#include "GPUO2DataTypes.h"
49#include "GPUParam.inc"
51#include "GPUTPCConvertImpl.h"
52#include "TPCFastTransformPOD.h"
53#include "GPUROOTDump.h"
56#include "GPUSettings.h"
57#include "GPUDefMacros.h"
58#ifndef GPUCA_STANDALONE
68#include "TPDGCode.h"
69#include "TParticlePDG.h"
70#include "TDatabasePDG.h"
71#endif
72#include "GPUQAHelper.h"
73#include <algorithm>
74#include <cstdio>
75#include <cinttypes>
76#include <fstream>
77
78#include "utils/timer.h"
79
80#include <oneapi/tbb.h>
81
82using namespace o2::gpu;
83
84namespace o2::gpu
85{
87 bool unattached = false;
88 float qpt = 0.f;
89 bool lowPt = false;
90 bool mev200 = false;
93 int32_t id = 0;
94 bool physics = false, protect = false;
95};
96} // namespace o2::gpu
97
98template <bool COUNT, class T>
99inline checkClusterStateResult GPUQA::checkClusterState(uint32_t attach, T* counts) const
100{
102 r.unattached = attach == 0;
104 if (!r.unattached && !(attach & gputpcgmmergertypes::attachProtect)) {
105 r.qpt = fabsf(mTracking->mIOPtrs.mergedTracks[r.id].GetParam().GetQPt());
106 r.lowPt = r.qpt * mTracking->GetParam().qptB5Scaler > mTracking->GetParam().rec.tpc.rejectQPtB5;
107 r.mev200 = r.qpt > 5;
108 r.mergedLooperUnconnected = mTracking->mIOPtrs.mergedTracks[r.id].MergedLooperUnconnected();
109 r.mergedLooperConnected = mTracking->mIOPtrs.mergedTracks[r.id].MergedLooperConnected();
110 }
111 if (r.mev200) {
112 if constexpr (COUNT) {
113 counts->n200MeV++;
114 }
115 }
116 if (r.lowPt) {
117 if constexpr (COUNT) {
118 counts->nLowPt++;
119 }
120 } else if (r.mergedLooperUnconnected) {
121 if constexpr (COUNT) {
122 counts->nMergedLooperUnconnected++;
123 }
124 } else if (r.mergedLooperConnected) {
125 if constexpr (COUNT) {
126 counts->nMergedLooperConnected++;
127 }
128 } else if (attach) {
129 r.protect = !GPUTPCClusterRejection::GetRejectionStatus<COUNT>(attach, r.physics, counts, &r.mev200) && ((attach & gputpcgmmergertypes::attachProtect) || !GPUTPCClusterRejection::IsTrackRejected(mTracking->mIOPtrs.mergedTracks[r.id], mTracking->GetParam()));
130 }
131 return r;
132}
133
134static const GPUSettingsQA& GPUQA_GetConfig(GPUChainTracking* chain)
135{
136 static GPUSettingsQA defaultConfig;
137 if (chain && chain->mConfigQA) {
138 return *chain->mConfigQA;
139 } else {
140 return defaultConfig;
141 }
142}
143
144static const constexpr float LOG_PT_MIN = -1.;
145
146static constexpr float Y_MAX = 40;
147static constexpr float Z_MAX = 100;
148static constexpr float PT_MIN = 0.01; // TODO: Take from Param
149static constexpr float PT_MIN_PRIM = 0.1;
150static constexpr float PT_MIN_CLUST = 0.01;
151static constexpr float PT_MAX = 20;
152static constexpr float ETA_MAX = 1.5;
153static constexpr float ETA_MAX2 = 0.9;
154static constexpr int32_t PADROW_CHECK_MINCLS = 50;
155
156static constexpr bool CLUST_HIST_INT_SUM = false;
157
158static constexpr const int32_t COLORCOUNT = 12;
159
160static const constexpr char* EFF_TYPES[6] = {"Rec", "Clone", "Fake", "All", "RecAndClone", "MC"};
161static const constexpr char* FINDABLE_NAMES[2] = {"All", "Findable"};
162static const constexpr char* PRIM_NAMES[2] = {"Prim", "Sec"};
163static const constexpr char* PARAMETER_NAMES[5] = {"Y", "Z", "#Phi", "#lambda", "Relative #it{p}_{T}"};
164static const constexpr char* PARAMETER_NAMES_NATIVE[5] = {"Y", "Z", "sin(#Phi)", "tan(#lambda)", "q/#it{p}_{T} (curvature)"};
165static const constexpr char* VSPARAMETER_NAMES[6] = {"Y", "Z", "Phi", "Eta", "Pt", "Pt_log"};
166static const constexpr char* EFF_NAMES[3] = {"Efficiency", "Clone Rate", "Fake Rate"};
167static const constexpr char* EFFICIENCY_TITLES[4] = {"Efficiency (Primary Tracks, Findable)", "Efficiency (Secondary Tracks, Findable)", "Efficiency (Primary Tracks)", "Efficiency (Secondary Tracks)"};
168static const constexpr double SCALE[5] = {10., 10., 1000., 1000., 100.};
169static const constexpr double SCALE_NATIVE[5] = {10., 10., 1000., 1000., 1.};
170static const constexpr char* XAXIS_TITLES[5] = {"#it{y}_{mc} (cm)", "#it{z}_{mc} (cm)", "#Phi_{mc} (rad)", "#eta_{mc}", "#it{p}_{Tmc} (GeV/#it{c})"};
171static const constexpr char* AXIS_TITLES[5] = {"#it{y}-#it{y}_{mc} (mm) (Resolution)", "#it{z}-#it{z}_{mc} (mm) (Resolution)", "#phi-#phi_{mc} (mrad) (Resolution)", "#lambda-#lambda_{mc} (mrad) (Resolution)", "(#it{p}_{T} - #it{p}_{Tmc}) / #it{p}_{Tmc} (%) (Resolution)"};
172static const constexpr char* AXIS_TITLES_NATIVE[5] = {"#it{y}-#it{y}_{mc} (mm) (Resolution)", "#it{z}-#it{z}_{mc} (mm) (Resolution)", "sin(#phi)-sin(#phi_{mc}) (Resolution)", "tan(#lambda)-tan(#lambda_{mc}) (Resolution)", "q*(q/#it{p}_{T} - q/#it{p}_{Tmc}) (Resolution)"};
173static const constexpr char* AXIS_TITLES_PULL[5] = {"#it{y}-#it{y}_{mc}/#sigma_{y} (Pull)", "#it{z}-#it{z}_{mc}/#sigma_{z} (Pull)", "sin(#phi)-sin(#phi_{mc})/#sigma_{sin(#phi)} (Pull)", "tan(#lambda)-tan(#lambda_{mc})/#sigma_{tan(#lambda)} (Pull)",
174 "q*(q/#it{p}_{T} - q/#it{p}_{Tmc})/#sigma_{q/#it{p}_{T}} (Pull)"};
175static const constexpr char* CLUSTER_NAMES[GPUQA::N_CLS_HIST] = {"Correctly attached clusters", "Fake attached clusters", "Attached + adjacent clusters", "Fake adjacent clusters", "Clusters of reconstructed tracks", "Used in Physics", "Protected", "All clusters"};
176static const constexpr char* CLUSTER_TITLES[GPUQA::N_CLS_TYPE] = {"Clusters Pt Distribution / Attachment", "Clusters Pt Distribution / Attachment (relative to all clusters)", "Clusters Pt Distribution / Attachment (integrated)"};
177static const constexpr char* CLUSTER_NAMES_SHORT[GPUQA::N_CLS_HIST] = {"Attached", "Fake", "AttachAdjacent", "FakeAdjacent", "FoundTracks", "Physics", "Protected", "All"};
178static const constexpr char* CLUSTER_TYPES[GPUQA::N_CLS_TYPE] = {"", "Ratio", "Integral"};
179static const constexpr char* REJECTED_NAMES[3] = {"All", "Rejected", "Fraction"};
180static const constexpr int32_t COLORS_HEX[COLORCOUNT] = {0xB03030, 0x00A000, 0x0000C0, 0x9400D3, 0x19BBBF, 0xF25900, 0x7F7F7F, 0xFFD700, 0x07F707, 0x07F7F7, 0xF08080, 0x000000};
181
182static const constexpr int32_t CONFIG_DASHED_MARKERS = 0;
183
184static const constexpr float AXES_MIN[5] = {-Y_MAX, -Z_MAX, 0.f, -ETA_MAX, PT_MIN};
185static const constexpr float AXES_MAX[5] = {Y_MAX, Z_MAX, 2.f * M_PI, ETA_MAX, PT_MAX};
186static const constexpr int32_t AXIS_BINS[5] = {51, 51, 144, 31, 50};
187static const constexpr int32_t RES_AXIS_BINS[] = {1017, 113}; // Consecutive bin sizes, histograms are binned down until the maximum entry is 50, each bin size should evenly divide its predecessor.
188static const constexpr float RES_AXES[5] = {1., 1., 0.03, 0.03, 1.0};
189static const constexpr float RES_AXES_NATIVE[5] = {1., 1., 0.1, 0.1, 5.0};
190static const constexpr float PULL_AXIS = 10.f;
191
192std::vector<TColor*> GPUQA::mColors;
193int32_t GPUQA::initColors()
194{
195 mColors.reserve(COLORCOUNT);
196 for (int32_t i = 0; i < COLORCOUNT; i++) {
197 float f1 = (float)((COLORS_HEX[i] >> 16) & 0xFF) / (float)0xFF;
198 float f2 = (float)((COLORS_HEX[i] >> 8) & 0xFF) / (float)0xFF;
199 float f3 = (float)((COLORS_HEX[i] >> 0) & 0xFF) / (float)0xFF;
200 mColors.emplace_back(new TColor(10000 + i, f1, f2, f3));
201 }
202 return 0;
203}
204static constexpr Color_t defaultColorNums[COLORCOUNT] = {kRed, kBlue, kGreen, kMagenta, kOrange, kAzure, kBlack, kYellow, kGray, kTeal, kSpring, kPink};
205
206#define TRACK_EXPECTED_REFERENCE_X_DEFAULT 81
207#ifndef GPUCA_RUN2 // Run 3 implementation
208static inline int32_t GPUQA_O2_ConvertFakeLabel(int32_t label) { return label >= 0x7FFFFFFE ? -1 : label; }
209inline uint32_t GPUQA::GetNMCCollissions() const { return mMCInfosCol.size(); }
210inline uint32_t GPUQA::GetNMCTracks(int32_t iCol) const { return mMCInfosCol[iCol].num; }
211inline uint32_t GPUQA::GetNMCTracks(const mcLabelI_t& label) const { return mMCInfosCol[mMCEventOffset[label.getSourceID()] + label.getEventID()].num; }
212inline uint32_t GPUQA::GetNMCLabels() const { return mClNative->clustersMCTruth ? mClNative->clustersMCTruth->getIndexedSize() : 0; }
213inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(uint32_t iTrk, uint32_t iCol) { return mMCInfos[mMCInfosCol[iCol].first + iTrk]; }
214inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(const mcLabel_t& label) { return mMCInfos[mMCInfosCol[mMCEventOffset[label.getSourceID()] + label.getEventID()].first + label.getTrackID()]; }
215inline GPUQA::mcLabels_t GPUQA::GetMCLabel(uint32_t i) { return mClNative->clustersMCTruth->getLabels(i); }
216inline int32_t GPUQA::GetMCLabelNID(const mcLabels_t& label) { return label.size(); }
217inline int32_t GPUQA::GetMCLabelNID(uint32_t i) { return mClNative->clustersMCTruth->getLabels(i).size(); }
218inline GPUQA::mcLabel_t GPUQA::GetMCLabel(uint32_t i, uint32_t j) { return mClNative->clustersMCTruth->getLabels(i)[j]; }
219inline int32_t GPUQA::GetMCLabelID(uint32_t i, uint32_t j) { return GPUQA_O2_ConvertFakeLabel(mClNative->clustersMCTruth->getLabels(i)[j].getTrackID()); }
220inline int32_t GPUQA::GetMCLabelID(const mcLabels_t& label, uint32_t j) { return GPUQA_O2_ConvertFakeLabel(label[j].getTrackID()); }
221inline int32_t GPUQA::GetMCLabelID(const mcLabel_t& label) { return GPUQA_O2_ConvertFakeLabel(label.getTrackID()); }
222inline uint32_t GPUQA::GetMCLabelCol(uint32_t i, uint32_t j) { return mMCEventOffset[mClNative->clustersMCTruth->getLabels(i)[j].getSourceID()] + mClNative->clustersMCTruth->getLabels(i)[j].getEventID(); }
223inline const auto& GPUQA::GetClusterLabels() { return mClNative->clustersMCTruth; }
224inline float GPUQA::GetMCLabelWeight(uint32_t i, uint32_t j) { return 1; }
225inline float GPUQA::GetMCLabelWeight(const mcLabels_t& label, uint32_t j) { return 1; }
226inline float GPUQA::GetMCLabelWeight(const mcLabel_t& label) { return 1; }
227inline bool GPUQA::mcPresent() { return !mConfig.noMC && mTracking && mClNative && mClNative->clustersMCTruth && mMCInfos.size(); }
228uint32_t GPUQA::GetMCLabelCol(const mcLabel_t& label) const { return !label.isValid() ? 0 : (mMCEventOffset[label.getSourceID()] + label.getEventID()); }
229GPUQA::mcLabelI_t GPUQA::GetMCTrackLabel(uint32_t trackId) const { return trackId >= mTrackMCLabels.size() ? MCCompLabel() : mTrackMCLabels[trackId]; }
230bool GPUQA::CompareIgnoreFake(const mcLabelI_t& l1, const mcLabelI_t& l2) { return l1.compare(l2) >= 0; }
231#define TRACK_EXPECTED_REFERENCE_X 78
232#else // Run 2 implementation
233inline GPUQA::mcLabelI_t::mcLabelI_t(const GPUQA::mcLabel_t& l) : track(l.fMCID) {}
234inline bool GPUQA::mcLabelI_t::operator==(const GPUQA::mcLabel_t& l) { return AbsLabelID(track) == l.fMCID; }
235inline uint32_t GPUQA::GetNMCCollissions() const { return 1; }
236inline uint32_t GPUQA::GetNMCTracks(int32_t iCol) const { return mTracking->mIOPtrs.nMCInfosTPC; }
237inline uint32_t GPUQA::GetNMCTracks(const mcLabelI_t& label) const { return mTracking->mIOPtrs.nMCInfosTPC; }
238inline uint32_t GPUQA::GetNMCLabels() const { return mTracking->mIOPtrs.nMCLabelsTPC; }
239inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(uint32_t iTrk, uint32_t iCol) { return mTracking->mIOPtrs.mcInfosTPC[AbsLabelID(iTrk)]; }
240inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(const mcLabel_t& label) { return GetMCTrack(label.fMCID, 0); }
241inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(const mcLabelI_t& label) { return GetMCTrack(label.track, 0); }
242inline const GPUQA::mcLabels_t& GPUQA::GetMCLabel(uint32_t i) { return mTracking->mIOPtrs.mcLabelsTPC[i]; }
243inline const GPUQA::mcLabel_t& GPUQA::GetMCLabel(uint32_t i, uint32_t j) { return mTracking->mIOPtrs.mcLabelsTPC[i].fClusterID[j]; }
244inline int32_t GPUQA::GetMCLabelNID(const mcLabels_t& label) { return 3; }
245inline int32_t GPUQA::GetMCLabelNID(uint32_t i) { return 3; }
246inline int32_t GPUQA::GetMCLabelID(uint32_t i, uint32_t j) { return mTracking->mIOPtrs.mcLabelsTPC[i].fClusterID[j].fMCID; }
247inline int32_t GPUQA::GetMCLabelID(const mcLabels_t& label, uint32_t j) { return label.fClusterID[j].fMCID; }
248inline int32_t GPUQA::GetMCLabelID(const mcLabel_t& label) { return label.fMCID; }
249inline uint32_t GPUQA::GetMCLabelCol(uint32_t i, uint32_t j) { return 0; }
250inline const auto& GPUQA::GetClusterLabels() { return mTracking->mIOPtrs.mcLabelsTPC; }
251inline float GPUQA::GetMCLabelWeight(uint32_t i, uint32_t j) { return mTracking->mIOPtrs.mcLabelsTPC[i].fClusterID[j].fWeight; }
252inline float GPUQA::GetMCLabelWeight(const mcLabels_t& label, uint32_t j) { return label.fClusterID[j].fWeight; }
253inline float GPUQA::GetMCLabelWeight(const mcLabel_t& label) { return label.fWeight; }
254inline int32_t GPUQA::FakeLabelID(int32_t id) { return id < 0 ? id : (-2 - id); }
255inline int32_t GPUQA::AbsLabelID(int32_t id) { return id >= 0 ? id : (-id - 2); }
256inline bool GPUQA::mcPresent() { return !mConfig.noMC && mTracking && GetNMCLabels() && GetNMCTracks(0); }
257uint32_t GPUQA::GetMCLabelCol(const mcLabel_t& label) const { return 0; }
258GPUQA::mcLabelI_t GPUQA::GetMCTrackLabel(uint32_t trackId) const { return trackId >= mTrackMCLabels.size() ? mcLabelI_t() : mTrackMCLabels[trackId]; }
259bool GPUQA::CompareIgnoreFake(const mcLabelI_t& l1, const mcLabelI_t& l2) { return AbsLabelID(l1.track) == AbsLabelID(l2.track); }
260#define TRACK_EXPECTED_REFERENCE_X TRACK_EXPECTED_REFERENCE_X_DEFAULT
261#endif
262template <class T>
263inline auto& GPUQA::GetMCTrackObj(T& obj, const GPUQA::mcLabelI_t& l)
264{
265 return obj[mMCEventOffset[l.getSourceID()] + l.getEventID()][l.getTrackID()];
266}
267
268template <>
269auto GPUQA::getHistArray<TH1F>()
270{
271 return std::make_pair(mHist1D, &mHist1D_pos);
272}
273template <>
274auto GPUQA::getHistArray<TH2F>()
275{
276 return std::make_pair(mHist2D, &mHist2D_pos);
277}
278template <>
279auto GPUQA::getHistArray<TH1D>()
280{
281 return std::make_pair(mHist1Dd, &mHist1Dd_pos);
282}
283template <>
284auto GPUQA::getHistArray<TGraphAsymmErrors>()
285{
286 return std::make_pair(mHistGraph, &mHistGraph_pos);
287}
288template <class T, typename... Args>
289void GPUQA::createHist(T*& h, const char* name, Args... args)
290{
291 const auto& p = getHistArray<T>();
292 if (mHaveExternalHists) {
293 if (p.first->size() <= p.second->size()) {
294 GPUError("Array sizes mismatch: Histograms %lu <= Positions %lu", p.first->size(), p.second->size());
295 throw std::runtime_error("Incoming histogram array incomplete");
296 }
297 if (strcmp((*p.first)[p.second->size()].GetName(), name)) {
298 GPUError("Histogram name mismatch: in array %s, trying to create %s", (*p.first)[p.second->size()].GetName(), name);
299 throw std::runtime_error("Incoming histogram has incorrect name");
300 }
301 } else {
302 if constexpr (std::is_same_v<T, TGraphAsymmErrors>) {
303 p.first->emplace_back();
304 p.first->back().SetName(name);
305 } else {
306 p.first->emplace_back(name, args...);
307 }
308 }
309 h = &((*p.first)[p.second->size()]);
310 p.second->emplace_back(&h);
311}
312
313namespace o2::gpu::internal
314{
316 std::tuple<std::vector<std::unique_ptr<TCanvas>>, std::vector<std::unique_ptr<TLegend>>, std::vector<std::unique_ptr<TPad>>, std::vector<std::unique_ptr<TLatex>>, std::vector<std::unique_ptr<TH1D>>> v;
317};
318} // namespace o2::gpu::internal
319
320template <class T, typename... Args>
321T* GPUQA::createGarbageCollected(Args... args)
322{
323 auto& v = std::get<std::vector<std::unique_ptr<T>>>(mGarbageCollector->v);
324 v.emplace_back(std::make_unique<T>(args...));
325 return v.back().get();
326}
327void GPUQA::clearGarbagageCollector()
328{
329 std::get<std::vector<std::unique_ptr<TPad>>>(mGarbageCollector->v).clear(); // Make sure to delete TPad first due to ROOT ownership (std::tuple has no defined order in its destructor)
330 std::apply([](auto&&... args) { ((args.clear()), ...); }, mGarbageCollector->v);
331}
332
333GPUQA::GPUQA(GPUChainTracking* chain, const GPUSettingsQA* config, const GPUParam* param) : mTracking(chain), mConfig(config ? *config : GPUQA_GetConfig(chain)), mParam(param ? param : &chain->GetParam()), mGarbageCollector(std::make_unique<internal::GPUQAGarbageCollection>())
334{
335 mMCEventOffset.resize(1, 0);
336}
337
339{
340 if (mQAInitialized && !mHaveExternalHists) {
341 delete mHist1D;
342 delete mHist2D;
343 delete mHist1Dd;
344 delete mHistGraph;
345 }
346 clearGarbagageCollector(); // Needed to guarantee correct order for ROOT ownership
347}
348
349bool GPUQA::clusterRemovable(int32_t attach, bool prot) const
350{
351 const auto& r = checkClusterState<false>(attach);
352 if (prot) {
353 return r.protect || r.physics;
354 }
355 return (!r.unattached && !r.physics && !r.protect);
356}
357
358template <class T>
359void GPUQA::SetAxisSize(T* e)
360{
361 e->GetYaxis()->SetTitleOffset(1.0);
362 e->GetYaxis()->SetTitleSize(0.045);
363 e->GetYaxis()->SetLabelSize(0.045);
364 e->GetXaxis()->SetTitleOffset(1.03);
365 e->GetXaxis()->SetTitleSize(0.045);
366 e->GetXaxis()->SetLabelOffset(-0.005);
367 e->GetXaxis()->SetLabelSize(0.045);
368}
369
370void GPUQA::SetLegend(TLegend* l, bool bigText)
371{
372 l->SetTextFont(72);
373 l->SetTextSize(bigText ? 0.03 : 0.016);
374 l->SetFillColor(0);
375}
376
377double* GPUQA::CreateLogAxis(int32_t nbins, float xmin, float xmax)
378{
379 float logxmin = std::log10(xmin);
380 float logxmax = std::log10(xmax);
381 float binwidth = (logxmax - logxmin) / nbins;
382
383 double* xbins = new double[nbins + 1];
384
385 xbins[0] = xmin;
386 for (int32_t i = 1; i <= nbins; i++) {
387 xbins[i] = std::pow(10, logxmin + i * binwidth);
388 }
389 return xbins;
390}
391
392void GPUQA::ChangePadTitleSize(TPad* p, float size)
393{
394 p->Update();
395 TPaveText* pt = (TPaveText*)(p->GetPrimitive("title"));
396 if (pt == nullptr) {
397 GPUError("Error changing title");
398 } else {
399 pt->SetTextSize(size);
400 p->Modified();
401 }
402}
403
404void GPUQA::DrawHisto(TH1* histo, char* filename, char* options)
405{
406 TCanvas tmp;
407 tmp.cd();
408 histo->Draw(options);
409 tmp.Print(filename);
410}
411
412void GPUQA::doPerfFigure(float x, float y, float size)
413{
414 if (mConfig.perfFigure == "") {
415 return;
416 }
417 static constexpr const char* str_perf_figure_1 = "ALICE Performance";
418 static constexpr const char* str_perf_figure_2_mc = "MC, Pb#minusPb, #sqrt{s_{NN}} = 5.36 TeV";
419 static constexpr const char* str_perf_figure_2_data = "Pb#minusPb, #sqrt{s_{NN}} = 5.36 TeV";
420 const char* str_perf_figure_2 = (mConfig.perfFigure == "mc" || mConfig.perfFigure == "MC") ? str_perf_figure_2_mc : (mConfig.perfFigure == "data" ? str_perf_figure_2_data : mConfig.perfFigure.c_str());
421
422 TLatex* t = createGarbageCollected<TLatex>(); // TODO: We could perhaps put everything in a legend, to get a white background if there is a grid
423 t->SetNDC(kTRUE);
424 t->SetTextColor(1);
425 t->SetTextSize(size);
426 t->DrawLatex(x, y, str_perf_figure_1);
427 t->SetTextSize(size * 0.8);
428 t->DrawLatex(x, y - 0.01 - size, str_perf_figure_2);
429}
430
431void GPUQA::SetMCTrackRange(int32_t min, int32_t max)
432{
433 mMCTrackMin = min;
434 mMCTrackMax = max;
435}
436
437int32_t GPUQA::InitQACreateHistograms()
438{
439 char name[2048], fname[1024];
440 if (mQATasks & taskTrackingEff) {
441 // Create Efficiency Histograms
442 for (int32_t i = 0; i < 6; i++) {
443 for (int32_t j = 0; j < 2; j++) {
444 for (int32_t k = 0; k < 2; k++) {
445 for (int32_t l = 0; l < 5; l++) {
446 snprintf(name, 2048, "%s%s%s%sVs%s", "tracks", EFF_TYPES[i], FINDABLE_NAMES[j], PRIM_NAMES[k], VSPARAMETER_NAMES[l]);
447 if (l == 4) {
448 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], k == 0 ? PT_MIN_PRIM : AXES_MIN[4], AXES_MAX[4])};
449 createHist(mEff[i][j][k][l], name, name, AXIS_BINS[l], binsPt.get());
450 } else {
451 createHist(mEff[i][j][k][l], name, name, AXIS_BINS[l], AXES_MIN[l], AXES_MAX[l]);
452 }
453 if (!mHaveExternalHists) {
454 mEff[i][j][k][l]->Sumw2();
455 }
456 strcat(name, "_eff");
457 if (i < 4) {
458 createHist(mEffResult[i][j][k][l], name);
459 }
460 }
461 }
462 }
463 }
464 }
465
466 // Create Resolution Histograms
467 if (mQATasks & taskTrackingRes) {
468 for (int32_t i = 0; i < 5; i++) {
469 for (int32_t j = 0; j < 5; j++) {
470 snprintf(name, 2048, "rms_%s_vs_%s", VSPARAMETER_NAMES[i], VSPARAMETER_NAMES[j]);
471 snprintf(fname, 1024, "mean_%s_vs_%s", VSPARAMETER_NAMES[i], VSPARAMETER_NAMES[j]);
472 if (j == 4) {
473 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], mConfig.resPrimaries == 1 ? PT_MIN_PRIM : AXES_MIN[4], AXES_MAX[4])};
474 createHist(mRes[i][j][0], name, name, AXIS_BINS[j], binsPt.get());
475 createHist(mRes[i][j][1], fname, fname, AXIS_BINS[j], binsPt.get());
476 } else {
477 createHist(mRes[i][j][0], name, name, AXIS_BINS[j], AXES_MIN[j], AXES_MAX[j]);
478 createHist(mRes[i][j][1], fname, fname, AXIS_BINS[j], AXES_MIN[j], AXES_MAX[j]);
479 }
480 snprintf(name, 2048, "res_%s_vs_%s", VSPARAMETER_NAMES[i], VSPARAMETER_NAMES[j]);
481 const float* axis = mConfig.nativeFitResolutions ? RES_AXES_NATIVE : RES_AXES;
482 const int32_t nbins = i == 4 && mConfig.nativeFitResolutions ? (10 * RES_AXIS_BINS[0]) : RES_AXIS_BINS[0];
483 if (j == 4) {
484 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], mConfig.resPrimaries == 1 ? PT_MIN_PRIM : AXES_MIN[4], AXES_MAX[4])};
485 createHist(mRes2[i][j], name, name, nbins, -axis[i], axis[i], AXIS_BINS[j], binsPt.get());
486 } else {
487 createHist(mRes2[i][j], name, name, nbins, -axis[i], axis[i], AXIS_BINS[j], AXES_MIN[j], AXES_MAX[j]);
488 }
489 }
490 }
491 }
492
493 // Create Pull Histograms
494 if (mQATasks & taskTrackingResPull) {
495 for (int32_t i = 0; i < 5; i++) {
496 for (int32_t j = 0; j < 5; j++) {
497 snprintf(name, 2048, "pull_rms_%s_vs_%s", VSPARAMETER_NAMES[i], VSPARAMETER_NAMES[j]);
498 snprintf(fname, 1024, "pull_mean_%s_vs_%s", VSPARAMETER_NAMES[i], VSPARAMETER_NAMES[j]);
499 if (j == 4) {
500 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], AXES_MIN[4], AXES_MAX[4])};
501 createHist(mPull[i][j][0], name, name, AXIS_BINS[j], binsPt.get());
502 createHist(mPull[i][j][1], fname, fname, AXIS_BINS[j], binsPt.get());
503 } else {
504 createHist(mPull[i][j][0], name, name, AXIS_BINS[j], AXES_MIN[j], AXES_MAX[j]);
505 createHist(mPull[i][j][1], fname, fname, AXIS_BINS[j], AXES_MIN[j], AXES_MAX[j]);
506 }
507 snprintf(name, 2048, "pull_%s_vs_%s", VSPARAMETER_NAMES[i], VSPARAMETER_NAMES[j]);
508 if (j == 4) {
509 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], AXES_MIN[4], AXES_MAX[4])};
510 createHist(mPull2[i][j], name, name, RES_AXIS_BINS[0], -PULL_AXIS, PULL_AXIS, AXIS_BINS[j], binsPt.get());
511 } else {
512 createHist(mPull2[i][j], name, name, RES_AXIS_BINS[0], -PULL_AXIS, PULL_AXIS, AXIS_BINS[j], AXES_MIN[j], AXES_MAX[j]);
513 }
514 }
515 }
516 }
517
518 // Create Cluster Histograms
519 if (mQATasks & taskClusterAttach) {
520 for (int32_t i = 0; i < N_CLS_TYPE * N_CLS_HIST - 1; i++) {
521 int32_t ioffset = i >= (2 * N_CLS_HIST - 1) ? (2 * N_CLS_HIST - 1) : (i >= N_CLS_HIST ? N_CLS_HIST : 0);
522 int32_t itype = i >= (2 * N_CLS_HIST - 1) ? 2 : (i >= N_CLS_HIST ? 1 : 0);
523 snprintf(name, 2048, "clusters%s%s", CLUSTER_NAMES_SHORT[i - ioffset], CLUSTER_TYPES[itype]);
524 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], PT_MIN_CLUST, PT_MAX)};
525 createHist(mClusters[i], name, name, AXIS_BINS[4], binsPt.get());
526 }
527
528 const int32_t maxPads = GPUTPCGeometry::NPads(GPUTPCGeometry::NROWS - 1);
529 createHist(mLowestRowDifferenceVsPad, "lowest_row_difference_vs_pad", "Lowest cluster row difference;Pad of lowest reconstructed cluster;Lowest reconstructed row - lowest MC row", maxPads, -0.5, maxPads - 0.5, 2 * GPUTPCGeometry::NROWS - 1, -GPUTPCGeometry::NROWS + 0.5, GPUTPCGeometry::NROWS - 0.5);
530 createHist(mLowestRowTrackMultiplicity, "lowest_row_track_multiplicity", "Track multiplicity for missing inner clusters;Reconstructed tracks with the same MC label;Selected reconstructed tracks", 100, -0.5, 99.5);
531
532 createHist(mPadRow[0], "padrow0", "padrow0", GPUTPCGeometry::NROWS - PADROW_CHECK_MINCLS, 0, GPUTPCGeometry::NROWS - 1 - PADROW_CHECK_MINCLS, GPUTPCGeometry::NROWS - PADROW_CHECK_MINCLS, 0, GPUTPCGeometry::NROWS - 1 - PADROW_CHECK_MINCLS);
533 createHist(mPadRow[1], "padrow1", "padrow1", 100.f, -0.2f, 0.2f, GPUTPCGeometry::NROWS - PADROW_CHECK_MINCLS, 0, GPUTPCGeometry::NROWS - 1 - PADROW_CHECK_MINCLS);
534 createHist(mPadRow[2], "padrow2", "padrow2", 100.f, -0.2f, 0.2f, GPUTPCGeometry::NROWS - PADROW_CHECK_MINCLS, 0, GPUTPCGeometry::NROWS - 1 - PADROW_CHECK_MINCLS);
535 createHist(mPadRow[3], "padrow3", "padrow3", 100.f, 0, 300000, GPUTPCGeometry::NROWS - PADROW_CHECK_MINCLS, 0, GPUTPCGeometry::NROWS - 1 - PADROW_CHECK_MINCLS);
536 }
537
538 if (mQATasks & taskTrackStatistics) {
539 // Create Tracks Histograms
540 for (int32_t i = 0; i < 2; i++) {
541 snprintf(name, 2048, i ? "nrows_with_cluster" : "nclusters");
542 createHist(mNCl[i], name, name, GPUTPCGeometry::NROWS + 1, 0, GPUTPCGeometry::NROWS);
543 }
544 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], PT_MIN_CLUST, PT_MAX)};
545 createHist(mTrackPt, "tracks_pt", "tracks_pt", AXIS_BINS[4], binsPt.get());
546 for (int32_t i = 0; i < 2; i++) {
547 snprintf(name, 2048, i ? "tracks_dedx_max" : "tracks_dedx_tot");
548 createHist(mTrackdEdx[i], name, name, 200, -3, 1, 1000, 0, i ? 1500 : 5000);
549 }
550 const uint32_t maxTime = (mTracking && mTracking->GetParam().continuousMaxTimeBin > 0) ? mTracking->GetParam().continuousMaxTimeBin : constants::TPC_MAX_TIME_BIN_TRIGGERED;
551 createHist(mT0[0], "tracks_t0", "tracks_t0", (maxTime + 1) / 10, 0, maxTime);
552 createHist(mT0[1], "tracks_t0_res", "tracks_t0_res", 1000, -100, 100);
553 createHist(mClXY, "clXY", "clXY", 1000, -250, 250, 1000, -250, 250); // TODO: Pass name only once
554 }
555 if (mQATasks & taskClusterRejection) {
556 const int padCount = GPUTPCGeometry::NPads(GPUTPCGeometry::NROWS - 1);
557 for (int32_t i = 0; i < 3; i++) {
558 snprintf(name, 2048, "clrej_%d", i);
559 createHist(mClRej[i], name, name, 2 * padCount, -padCount / 2 + 0.5f, padCount / 2 - 0.5f, GPUTPCGeometry::NROWS, 0, GPUTPCGeometry::NROWS - 1);
560 }
561 createHist(mClRejP, "clrejp", "clrejp", GPUTPCGeometry::NROWS, 0, GPUTPCGeometry::NROWS - 1);
562 }
563
564 if ((mQATasks & taskClusterCounts) && mConfig.clusterRejectionHistograms) {
565 int32_t num = DoClusterCounts(nullptr, 2);
566 mHistClusterCount.resize(num);
567 DoClusterCounts(nullptr, 1);
568 }
569
570 for (uint32_t i = 0; i < mHist1D->size(); i++) {
571 *mHist1D_pos[i] = &(*mHist1D)[i];
572 }
573 for (uint32_t i = 0; i < mHist2D->size(); i++) {
574 *mHist2D_pos[i] = &(*mHist2D)[i];
575 }
576 for (uint32_t i = 0; i < mHist1Dd->size(); i++) {
577 *mHist1Dd_pos[i] = &(*mHist1Dd)[i];
578 }
579 for (uint32_t i = 0; i < mHistGraph->size(); i++) {
580 *mHistGraph_pos[i] = &(*mHistGraph)[i];
581 }
582
583 return 0;
584}
585
586int32_t GPUQA::loadHistograms(std::vector<TH1F>& i1, std::vector<TH2F>& i2, std::vector<TH1D>& i3, std::vector<TGraphAsymmErrors>& i4, int32_t tasks)
587{
588 if (tasks == tasksAutomatic) {
590 }
591 if (mQAInitialized && (!mHaveExternalHists || tasks != mQATasks)) {
592 throw std::runtime_error("QA not initialized or initialized with different task array");
593 }
594 mHist1D = &i1;
595 mHist2D = &i2;
596 mHist1Dd = &i3;
597 mHistGraph = &i4;
598 mHist1D_pos.clear();
599 mHist2D_pos.clear();
600 mHist1Dd_pos.clear();
601 mHistGraph_pos.clear();
602 mHaveExternalHists = true;
603 if (mConfig.noMC) {
604 tasks &= tasksAllNoQC;
605 }
606 mQATasks = tasks;
607 if (InitQACreateHistograms()) {
608 return 1;
609 }
610 mQAInitialized = true;
611 return 0;
612}
613
614void GPUQA::DumpO2MCData(const char* filename) const
615{
616 FILE* fp = fopen(filename, "w+b");
617 if (fp == nullptr) {
618 return;
619 }
620 uint32_t n = mMCInfos.size();
621 fwrite(&n, sizeof(n), 1, fp);
622 fwrite(mMCInfos.data(), sizeof(mMCInfos[0]), n, fp);
623 n = mMCInfosCol.size();
624 fwrite(&n, sizeof(n), 1, fp);
625 fwrite(mMCInfosCol.data(), sizeof(mMCInfosCol[0]), n, fp);
626 n = mMCEventOffset.size();
627 fwrite(&n, sizeof(n), 1, fp);
628 fwrite(mMCEventOffset.data(), sizeof(mMCEventOffset[0]), n, fp);
629 fclose(fp);
630}
631
633{
634 FILE* fp = fopen(filename, "rb");
635 if (fp == nullptr) {
636 return 1;
637 }
638 uint32_t n;
639 uint32_t x;
640 if ((x = fread(&n, sizeof(n), 1, fp)) != 1) {
641 fclose(fp);
642 return 1;
643 }
644 mMCInfos.resize(n);
645 if (fread(mMCInfos.data(), sizeof(mMCInfos[0]), n, fp) != n) {
646 fclose(fp);
647 return 1;
648 }
649 if ((x = fread(&n, sizeof(n), 1, fp)) != 1) {
650 fclose(fp);
651 return 1;
652 }
653 mMCInfosCol.resize(n);
654 if (fread(mMCInfosCol.data(), sizeof(mMCInfosCol[0]), n, fp) != n) {
655 fclose(fp);
656 return 1;
657 }
658 if ((x = fread(&n, sizeof(n), 1, fp)) != 1) {
659 fclose(fp);
660 return 1;
661 }
662 mMCEventOffset.resize(n);
663 if (fread(mMCEventOffset.data(), sizeof(mMCEventOffset[0]), n, fp) != n) {
664 fclose(fp);
665 return 1;
666 }
667 if (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2) {
668 printf("Read %ld bytes MC Infos\n", ftell(fp));
669 }
670 fclose(fp);
671 if (mTracking) {
672 CopyO2MCtoIOPtr(&mTracking->mIOPtrs);
673 }
674 return 0;
675}
676
677void GPUQA::CopyO2MCtoIOPtr(GPUTrackingInOutPointers* ptr)
678{
679 ptr->mcInfosTPC = mMCInfos.data();
680 ptr->nMCInfosTPC = mMCInfos.size();
681 ptr->mcInfosTPCCol = mMCInfosCol.data();
682 ptr->nMCInfosTPCCol = mMCInfosCol.size();
683}
684
685void GPUQA::InitO2MCData(GPUTrackingInOutPointers* updateIOPtr)
686{
687#ifndef GPUCA_STANDALONE
688 if (!mO2MCDataLoaded) {
689 HighResTimer timer(mTracking && mTracking->GetProcessingSettings().debugLevel);
690 if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
691 GPUInfo("Start reading O2 Track MC information");
692 }
693 static constexpr float PRIM_MAX_T = 0.01f;
694
695 o2::steer::MCKinematicsReader mcReader("collisioncontext.root");
696 std::vector<int32_t> refId;
697
698 auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
699 const auto& evrec = dc->getEventRecords();
700 const auto& evparts = dc->getEventParts();
701 std::vector<std::vector<float>> evTimeBins(mcReader.getNSources());
702 for (uint32_t i = 0; i < evTimeBins.size(); i++) {
703 evTimeBins[i].resize(mcReader.getNEvents(i), -100.f);
704 }
705 for (uint32_t i = 0; i < evrec.size(); i++) {
706 const auto& ir = evrec[i];
707 for (uint32_t j = 0; j < evparts[i].size(); j++) {
708 const int iSim = evparts[i][j].sourceID;
709 const int iEv = evparts[i][j].entryID;
710 if (iSim == o2::steer::QEDSOURCEID || ir.differenceInBC(o2::raw::HBFUtils::Instance().getFirstIR()) >= 0) {
712 float timebin = (float)ir.differenceInBC(ir0) / o2::tpc::constants::LHCBCPERTIMEBIN;
713 if (evTimeBins[iSim][iEv] >= 0) {
714 throw std::runtime_error("Multiple time bins for same MC collision found");
715 }
716 evTimeBins[iSim][iEv] = timebin;
717 }
718 }
719 }
720
721 uint32_t nSimSources = mcReader.getNSources();
722 mMCEventOffset.resize(nSimSources);
723 uint32_t nSimTotalEvents = 0;
724 uint32_t nSimTotalTracks = 0;
725 for (uint32_t i = 0; i < nSimSources; i++) {
726 mMCEventOffset[i] = nSimTotalEvents;
727 nSimTotalEvents += mcReader.getNEvents(i);
728 }
729
730 mMCInfosCol.resize(nSimTotalEvents);
731 for (int32_t iSim = 0; iSim < mcReader.getNSources(); iSim++) {
732 for (int32_t i = 0; i < mcReader.getNEvents(iSim); i++) {
733 const float timebin = evTimeBins[iSim][i];
734
735 const std::vector<o2::MCTrack>& tracks = mcReader.getTracks(iSim, i);
736 const std::vector<o2::TrackReference>& trackRefs = mcReader.getTrackRefsByEvent(iSim, i);
737
738 refId.resize(tracks.size());
739 std::fill(refId.begin(), refId.end(), -1);
740 for (uint32_t j = 0; j < trackRefs.size(); j++) {
741 if (trackRefs[j].getDetectorId() == o2::detectors::DetID::TPC) {
742 int32_t trkId = trackRefs[j].getTrackID();
743 if (refId[trkId] == -1) {
744 refId[trkId] = j;
745 }
746 }
747 }
748 mMCInfosCol[mMCEventOffset[iSim] + i].first = mMCInfos.size();
749 mMCInfosCol[mMCEventOffset[iSim] + i].num = tracks.size();
750 mMCInfos.resize(mMCInfos.size() + tracks.size());
751 for (uint32_t j = 0; j < tracks.size(); j++) {
752 auto& info = mMCInfos[mMCInfosCol[mMCEventOffset[iSim] + i].first + j];
753 const auto& trk = tracks[j];
754 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(trk.GetPdgCode());
755 Int_t pid = -1;
756 if (abs(trk.GetPdgCode()) == kElectron) {
757 pid = 0;
758 }
759 if (abs(trk.GetPdgCode()) == kMuonMinus) {
760 pid = 1;
761 }
762 if (abs(trk.GetPdgCode()) == kPiPlus) {
763 pid = 2;
764 }
765 if (abs(trk.GetPdgCode()) == kKPlus) {
766 pid = 3;
767 }
768 if (abs(trk.GetPdgCode()) == kProton) {
769 pid = 4;
770 }
771
772 info.charge = particle ? particle->Charge() : 0;
773 info.prim = trk.T() < PRIM_MAX_T;
774 info.primDaughters = 0;
775 if (trk.getFirstDaughterTrackId() != -1) {
776 for (int32_t k = trk.getFirstDaughterTrackId(); k <= trk.getLastDaughterTrackId(); k++) {
777 if (tracks[k].T() < PRIM_MAX_T) {
778 info.primDaughters = 1;
779 break;
780 }
781 }
782 }
783 info.pid = pid;
784 info.t0 = timebin;
785 if (refId[j] >= 0) {
786 const auto& trkRef = trackRefs[refId[j]];
787 info.x = trkRef.X();
788 info.y = trkRef.Y();
789 info.z = trkRef.Z();
790 info.pX = trkRef.Px();
791 info.pY = trkRef.Py();
792 info.pZ = trkRef.Pz();
793 info.genRadius = std::sqrt(trk.GetStartVertexCoordinatesX() * trk.GetStartVertexCoordinatesX() + trk.GetStartVertexCoordinatesY() * trk.GetStartVertexCoordinatesY() + trk.GetStartVertexCoordinatesZ() * trk.GetStartVertexCoordinatesZ());
794 } else {
795 info.x = info.y = info.z = info.pX = info.pY = info.pZ = 0;
796 info.genRadius = 0;
797 }
798 }
799 }
800 }
801 if (timer.IsRunning()) {
802 GPUInfo("Finished reading O2 Track MC information (%f seconds)", timer.GetCurrentElapsedTime());
803 }
804 mO2MCDataLoaded = true;
805 }
806 if (updateIOPtr) {
807 CopyO2MCtoIOPtr(updateIOPtr);
808 }
809#endif
810}
811
812int32_t GPUQA::InitQA(int32_t tasks)
813{
814 if (mQAInitialized) {
815 throw std::runtime_error("QA already initialized");
816 }
817 if (tasks == tasksAutomatic) {
818 tasks = tasksDefault;
819 }
820
821 mHist1D = new std::vector<TH1F>;
822 mHist2D = new std::vector<TH2F>;
823 mHist1Dd = new std::vector<TH1D>;
824 mHistGraph = new std::vector<TGraphAsymmErrors>;
825 if (mConfig.noMC) {
826 tasks &= tasksAllNoQC;
827 }
828 mQATasks = tasks;
829
830 if (mTracking->GetProcessingSettings().qcRunFraction != 100.f && mQATasks != taskClusterCounts) {
831 throw std::runtime_error("QA with qcRunFraction only supported for taskClusterCounts");
832 }
833
834 if (mTracking) {
835 mClNative = mTracking->mIOPtrs.clustersNative;
836 }
837
838 if (InitQACreateHistograms()) {
839 return 1;
840 }
841
842 if (mConfig.enableLocalOutput) {
843 mkdir(mConfig.plotsDir.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
844 }
845
846#ifndef GPUCA_STANDALONE
847 if (!mConfig.noMC) {
848 InitO2MCData(mTracking ? &mTracking->mIOPtrs : nullptr);
849 }
850#endif
851
852 if (mConfig.matchMCLabels.size()) {
853 uint32_t nFiles = mConfig.matchMCLabels.size();
854 std::vector<std::unique_ptr<TFile>> files;
855 std::vector<std::vector<std::vector<int32_t>>*> labelsBuffer(nFiles);
856 std::vector<std::vector<std::vector<int32_t>>*> effBuffer(nFiles);
857 for (uint32_t i = 0; i < nFiles; i++) {
858 files.emplace_back(std::make_unique<TFile>(mConfig.matchMCLabels[i].c_str()));
859 labelsBuffer[i] = (std::vector<std::vector<int32_t>>*)files[i]->Get("mcLabelBuffer");
860 effBuffer[i] = (std::vector<std::vector<int32_t>>*)files[i]->Get("mcEffBuffer");
861 if (labelsBuffer[i] == nullptr || effBuffer[i] == nullptr) {
862 GPUError("Error opening / reading from labels file %u/%s: %p %p", i, mConfig.matchMCLabels[i].c_str(), (void*)labelsBuffer[i], (void*)effBuffer[i]);
863 exit(1);
864 }
865 }
866
867 mGoodTracks.resize(labelsBuffer[0]->size());
868 mGoodHits.resize(labelsBuffer[0]->size());
869 for (uint32_t iEvent = 0; iEvent < labelsBuffer[0]->size(); iEvent++) {
870 std::vector<bool> labelsOK((*effBuffer[0])[iEvent].size());
871 for (uint32_t k = 0; k < (*effBuffer[0])[iEvent].size(); k++) {
872 labelsOK[k] = false;
873 for (uint32_t l = 0; l < nFiles; l++) {
874 if ((*effBuffer[0])[iEvent][k] != (*effBuffer[l])[iEvent][k]) {
875 labelsOK[k] = true;
876 break;
877 }
878 }
879 }
880 mGoodTracks[iEvent].resize((*labelsBuffer[0])[iEvent].size());
881 for (uint32_t k = 0; k < (*labelsBuffer[0])[iEvent].size(); k++) {
882 if ((*labelsBuffer[0])[iEvent][k] == MC_LABEL_INVALID) {
883 continue;
884 }
885 mGoodTracks[iEvent][k] = labelsOK[abs((*labelsBuffer[0])[iEvent][k])];
886 }
887 }
888 }
889 mQAInitialized = true;
890 return 0;
891}
892
893void GPUQA::RunQA(bool matchOnly, const std::vector<o2::tpc::TrackTPC>* tracksExternal, const std::vector<o2::MCCompLabel>* tracksExtMC, const o2::tpc::ClusterNativeAccess* clNative)
894{
895 if (!mQAInitialized) {
896 throw std::runtime_error("QA not initialized");
897 }
898 if (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2) {
899 GPUInfo("Running QA - Mask %d, Efficiency %d, Resolution %d, Pulls %d, Cluster Attachment %d, Track Statistics %d, Cluster Counts %d", mQATasks, (int32_t)(mQATasks & taskTrackingEff), (int32_t)(mQATasks & taskTrackingRes), (int32_t)(mQATasks & taskTrackingResPull), (int32_t)(mQATasks & taskClusterAttach), (int32_t)(mQATasks & taskTrackStatistics), (int32_t)(mQATasks & taskClusterCounts));
900 }
901 if (!clNative && mTracking) {
902 clNative = mTracking->mIOPtrs.clustersNative;
903 }
904 mClNative = clNative;
905
906#ifndef GPUCA_RUN2
907 uint32_t nSimEvents = GetNMCCollissions();
908 if (mTrackMCLabelsReverse.size() < nSimEvents) {
909 mTrackMCLabelsReverse.resize(nSimEvents);
910 }
911 if (mRecTracks.size() < nSimEvents) {
912 mRecTracks.resize(nSimEvents);
913 }
914 if (mFakeTracks.size() < nSimEvents) {
915 mFakeTracks.resize(nSimEvents);
916 }
917 if (mMCParam.size() < nSimEvents) {
918 mMCParam.resize(nSimEvents);
919 }
920#endif
921
922 // Initialize Arrays
923 uint32_t nReconstructedTracks = 0;
924 if (tracksExternal) {
925#ifndef GPUCA_STANDALONE
926 nReconstructedTracks = tracksExternal->size();
927#endif
928 } else {
929 nReconstructedTracks = mTracking->mIOPtrs.nMergedTracks;
930 }
931 mTrackMCLabels.resize(nReconstructedTracks);
932 for (uint32_t iCol = 0; iCol < GetNMCCollissions(); iCol++) {
933 mTrackMCLabelsReverse[iCol].resize(GetNMCTracks(iCol));
934 mRecTracks[iCol].resize(GetNMCTracks(iCol));
935 mFakeTracks[iCol].resize(GetNMCTracks(iCol));
936 mMCParam[iCol].resize(GetNMCTracks(iCol));
937 memset(mRecTracks[iCol].data(), 0, mRecTracks[iCol].size() * sizeof(mRecTracks[iCol][0]));
938 memset(mFakeTracks[iCol].data(), 0, mFakeTracks[iCol].size() * sizeof(mFakeTracks[iCol][0]));
939 for (size_t i = 0; i < mTrackMCLabelsReverse[iCol].size(); i++) {
940 mTrackMCLabelsReverse[iCol][i] = -1;
941 }
942 }
943 if (mQATasks & taskClusterAttach && GetNMCLabels()) {
944 mClusterParam.resize(GetNMCLabels());
945 memset(mClusterParam.data(), 0, mClusterParam.size() * sizeof(mClusterParam[0]));
946 }
947 HighResTimer timer(QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2));
948
949 mNEvents++;
950 if (mConfig.writeMCLabels) {
951 mcEffBuffer.resize(mNEvents);
952 mcLabelBuffer.resize(mNEvents);
953 mcEffBuffer[mNEvents - 1].resize(GetNMCTracks(0));
954 mcLabelBuffer[mNEvents - 1].resize(nReconstructedTracks);
955 }
956
957 bool mcAvail = mcPresent() || tracksExtMC;
958
959 if (mcAvail) { // Assign Track MC Labels
960 if (tracksExternal) {
961#ifndef GPUCA_STANDALONE
962 for (uint32_t i = 0; i < tracksExternal->size(); i++) {
963 mTrackMCLabels[i] = (*tracksExtMC)[i];
964 }
965#endif
966 } else {
967 tbb::parallel_for(tbb::blocked_range<uint32_t>(0, nReconstructedTracks, (QA_DEBUG == 0) ? 32 : nReconstructedTracks), [&](const tbb::blocked_range<uint32_t>& range) {
968 auto acc = GPUTPCTrkLbl<true, mcLabelI_t>(GetClusterLabels(), 1.f - mConfig.recThreshold);
969 for (auto i = range.begin(); i < range.end(); i++) {
970 acc.reset();
971 int32_t nClusters = 0;
972 const GPUTPCGMMergedTrack& track = mTracking->mIOPtrs.mergedTracks[i];
973 std::vector<mcLabel_t> labels;
974 for (uint32_t k = 0; k < track.NClusters(); k++) {
975 if (mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].state & GPUTPCGMMergedTrackHit::flagReject) {
976 continue;
977 }
978 nClusters++;
979 uint32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].num;
980 if (hitId >= GetNMCLabels()) {
981 GPUError("Invalid hit id %u > %d (nClusters %d)", hitId, GetNMCLabels(), clNative ? clNative->nClustersTotal : 0);
982 throw std::runtime_error("qa error");
983 }
984 acc.addLabel(hitId);
985 for (int32_t j = 0; j < GetMCLabelNID(hitId); j++) {
986 if (GetMCLabelID(hitId, j) >= (int32_t)GetNMCTracks(GetMCLabelCol(hitId, j))) {
987 GPUError("Invalid label %d > %d (hit %d, label %d, col %d)", GetMCLabelID(hitId, j), GetNMCTracks(GetMCLabelCol(hitId, j)), hitId, j, (int32_t)GetMCLabelCol(hitId, j));
988 throw std::runtime_error("qa error");
989 }
990 if (GetMCLabelID(hitId, j) >= 0) {
991 if (QA_DEBUG >= 3 && track.OK()) {
992 GPUInfo("Track %d Cluster %u Label %d: %d (%f)", i, k, j, GetMCLabelID(hitId, j), GetMCLabelWeight(hitId, j));
993 }
994 }
995 }
996 }
997
998 float maxweight, sumweight;
999 int32_t maxcount;
1000 auto maxLabel = acc.computeLabel(&maxweight, &sumweight, &maxcount);
1001 mTrackMCLabels[i] = maxLabel;
1002 if (QA_DEBUG && track.OK() && GetNMCTracks(maxLabel) > (uint32_t)maxLabel.getTrackID()) {
1003 const mcInfo_t& mc = GetMCTrack(maxLabel);
1004 GPUInfo("Track %d label %d (fake %d) weight %f clusters %d (fitted %d) (%f%% %f%%) Pt %f", i, maxLabel.getTrackID(), (int32_t)(maxLabel.isFake()), maxweight, nClusters, track.NClustersFitted(), 100.f * maxweight / sumweight, 100.f * (float)maxcount / (float)nClusters,
1005 std::sqrt(mc.pX * mc.pX + mc.pY * mc.pY));
1006 }
1007 }
1008 });
1009 }
1010 if (timer.IsRunning()) {
1011 GPUInfo("QA Time: Assign Track Labels:\t\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1012 }
1013
1014 for (uint32_t i = 0; i < nReconstructedTracks; i++) {
1015 const GPUTPCGMMergedTrack* track = mTracking ? &mTracking->mIOPtrs.mergedTracks[i] : nullptr;
1016 mcLabelI_t label = mTrackMCLabels[i];
1017 if (mQATasks & taskClusterAttach) {
1018 // fill cluster attachment status
1019 if (!track->OK()) {
1020 continue;
1021 }
1022 if (!mTrackMCLabels[i].isValid()) {
1023 for (uint32_t k = 0; k < track->NClusters(); k++) {
1024 if (mTracking->mIOPtrs.mergedTrackHits[track->FirstClusterRef() + k].state & GPUTPCGMMergedTrackHit::flagReject) {
1025 continue;
1026 }
1027 mClusterParam[mTracking->mIOPtrs.mergedTrackHits[track->FirstClusterRef() + k].num].fakeAttached++;
1028 }
1029 continue;
1030 }
1031 if (mMCTrackMin == -1 || (label.getTrackID() >= mMCTrackMin && label.getTrackID() < mMCTrackMax)) {
1032 for (uint32_t k = 0; k < track->NClusters(); k++) {
1033 if (mTracking->mIOPtrs.mergedTrackHits[track->FirstClusterRef() + k].state & GPUTPCGMMergedTrackHit::flagReject) {
1034 continue;
1035 }
1036 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track->FirstClusterRef() + k].num;
1037 bool correct = false;
1038 for (int32_t j = 0; j < GetMCLabelNID(hitId); j++) {
1039 if (label == GetMCLabel(hitId, j)) {
1040 correct = true;
1041 break;
1042 }
1043 }
1044 if (correct) {
1045 mClusterParam[hitId].attached++;
1046 } else {
1047 mClusterParam[hitId].fakeAttached++;
1048 }
1049 }
1050 }
1051 }
1052
1053 if (mTrackMCLabels[i].isFake()) {
1054 (GetMCTrackObj(mFakeTracks, label))++;
1055 } else if (tracksExternal || !track->MergedLooper()) {
1056 GetMCTrackObj(mRecTracks, label)++;
1057 if (mMCTrackMin == -1 || (label.getTrackID() >= mMCTrackMin && label.getTrackID() < mMCTrackMax)) {
1058 int32_t& revLabel = GetMCTrackObj(mTrackMCLabelsReverse, label);
1059 if (tracksExternal) {
1060#ifndef GPUCA_STANDALONE
1061 if (revLabel == -1 || fabsf((*tracksExternal)[i].getZ()) < fabsf((*tracksExternal)[revLabel].getZ())) {
1062 revLabel = i;
1063 }
1064#endif
1065 } else {
1066 const auto* trks = mTracking->mIOPtrs.mergedTracks;
1067 bool comp;
1068 if (revLabel == -1) {
1069 comp = true;
1070 } else {
1071 float shift1 = mTracking->GetTPCTransform()->convDeltaTimeToDeltaZinTimeFrame(trks[i].CSide() * GPUChainTracking::NSECTORS / 2, trks[i].GetParam().GetTOffset());
1072 float shift2 = mTracking->GetTPCTransform()->convDeltaTimeToDeltaZinTimeFrame(trks[revLabel].CSide() * GPUChainTracking::NSECTORS / 2, trks[revLabel].GetParam().GetTOffset());
1073 comp = fabsf(trks[i].GetParam().GetZ() + shift1) < fabsf(trks[revLabel].GetParam().GetZ() + shift2);
1074 }
1075 if (revLabel == -1 || !trks[revLabel].OK() || (trks[i].OK() && comp)) {
1076 revLabel = i;
1077 }
1078 }
1079 }
1080 }
1081 }
1082 if ((mQATasks & taskClusterAttach) && !tracksExternal) {
1083 std::vector<uint8_t> lowestPadRow(mTracking->mIOPtrs.nMergedTracks);
1084 // fill cluster adjacent status
1085 if (mTracking->mIOPtrs.mergedTrackHitAttachment) {
1086 for (uint32_t i = 0; i < GetNMCLabels(); i++) {
1087 if (mClusterParam[i].attached == 0 && mClusterParam[i].fakeAttached == 0) {
1088 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[i];
1091 mcLabelI_t trackL = mTrackMCLabels[track];
1092 bool fake = true;
1093 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1094 // GPUInfo("Attach %x Track %d / %d:%d", attach, track, j, GetMCLabelID(i, j));
1095 if (trackL == GetMCLabel(i, j)) {
1096 fake = false;
1097 break;
1098 }
1099 }
1100 if (fake) {
1101 mClusterParam[i].fakeAdjacent++;
1102 } else {
1103 mClusterParam[i].adjacent++;
1104 }
1105 }
1106 }
1107 }
1108 }
1109 if (mTracking->mIOPtrs.nMergedTracks && clNative) {
1110 std::fill(lowestPadRow.begin(), lowestPadRow.end(), 255);
1111 for (uint32_t iSector = 0; iSector < GPUTPCGeometry::NSECTORS; iSector++) {
1112 for (uint32_t iRow = 0; iRow < GPUTPCGeometry::NROWS; iRow++) {
1113 for (uint32_t iCl = 0; iCl < clNative->nClusters[iSector][iRow]; iCl++) {
1114 int32_t i = clNative->clusterOffset[iSector][iRow] + iCl;
1115 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1116 const mcLabelI_t label = GetMCLabel(i, j);
1117 if (!label.isValid() || label.isNoise()) {
1118 continue;
1119 }
1120 uint32_t trackId = GetMCTrackObj(mTrackMCLabelsReverse, label);
1121 if (trackId < lowestPadRow.size() && lowestPadRow[trackId] > iRow) {
1122 lowestPadRow[trackId] = iRow;
1123 }
1124 }
1125 }
1126 }
1127 }
1128 for (uint32_t i = 0; i < mTracking->mIOPtrs.nMergedTracks; i++) {
1129 const auto& trk = mTracking->mIOPtrs.mergedTracks[i];
1130 const auto& label = mTrackMCLabels[i];
1131 if (trk.OK() && trk.NClustersFitted() > 50 && CAMath::Abs(trk.GetParam().GetQPt()) < 1.f && label.isValid() && !label.isNoise() &&
1132 (mMCTrackMin == -1 || label.getTrackID() >= mMCTrackMin) && (mMCTrackMax == -1 || label.getTrackID() < mMCTrackMax)) {
1133 const uint32_t trackId = GetMCTrackObj(mTrackMCLabelsReverse, label);
1134 if (trackId < lowestPadRow.size() && lowestPadRow[trackId] < 10) {
1135 const GPUTPCGMMergedTrackHit* lowestCl = nullptr;
1136 for (uint32_t j = 0; j < trk.NClusters(); j++) {
1137 const auto& cl = mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef() + j];
1138 if (!(cl.state & GPUTPCGMMergedTrackHit::flagReject) && (!lowestCl || cl.row < lowestCl->row)) {
1139 lowestCl = &cl;
1140 }
1141 }
1142 if (lowestCl) {
1143 const float pad = clNative->clustersLinear[lowestCl->num].getPad();
1144 const int32_t difference = (int32_t)lowestCl->row - lowestPadRow[trackId];
1145 mLowestRowDifferenceVsPad->Fill(pad, difference);
1146 if (CAMath::Abs(difference) > 5 && pad > 5.f && pad < GPUTPCGeometry::NPads(lowestCl->row) - 5.f) {
1147 mLowestRowTrackMultiplicity->Fill(GetMCTrackObj(mRecTracks, label));
1148 }
1149 }
1150 }
1151 }
1152 if (trk.OK() && lowestPadRow[i] != 255 && trk.NClustersFitted() >= PADROW_CHECK_MINCLS && CAMath::Abs(trk.GetParam().GetQPt()) < 1.0) {
1153 const auto& lowestCl = mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef()].row < mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef() + trk.NClusters() - 1].row ? mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef()] : mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef() + trk.NClusters() - 1];
1154 const int32_t lowestRow = lowestCl.row;
1155 mPadRow[0]->Fill(lowestPadRow[i], lowestRow, 1.f);
1156 mPadRow[1]->Fill(CAMath::ATan2(trk.GetParam().GetY(), trk.GetParam().GetX()), lowestRow, 1.f);
1157 if (lowestPadRow[i] < 10 && lowestRow > lowestPadRow[i] + 3) {
1158 const auto& cl = clNative->clustersLinear[lowestCl.num];
1159 float x, y, z;
1160 mTracking->GetTPCTransform()->Transform(lowestCl.sector, lowestCl.row, cl.getPad(), cl.getTime(), x, y, z, trk.GetParam().GetTOffset());
1161 float phi = CAMath::ATan2(y, x);
1162 mPadRow[2]->Fill(phi, lowestRow, 1.f);
1163 if (CAMath::Abs(phi) < 0.15) {
1164 const float time = cl.getTime();
1165 mPadRow[3]->Fill(mTracking->GetParam().GetUnscaledMult(time), lowestRow, 1.f);
1166 }
1167 }
1168 }
1169 }
1170 }
1171 }
1172
1173 if (mConfig.matchMCLabels.size()) {
1174 mGoodHits[mNEvents - 1].resize(GetNMCLabels());
1175 std::vector<bool> allowMCLabels(GetNMCTracks(0));
1176 for (uint32_t k = 0; k < GetNMCTracks(0); k++) {
1177 allowMCLabels[k] = false;
1178 }
1179 for (uint32_t i = 0; i < nReconstructedTracks; i++) {
1180 if (!mGoodTracks[mNEvents - 1][i]) {
1181 continue;
1182 }
1183 if (mConfig.matchDisplayMinPt > 0) {
1184 if (!mTrackMCLabels[i].isValid()) {
1185 continue;
1186 }
1187 const mcInfo_t& info = GetMCTrack(mTrackMCLabels[i]);
1188 if (info.pX * info.pX + info.pY * info.pY < mConfig.matchDisplayMinPt * mConfig.matchDisplayMinPt) {
1189 continue;
1190 }
1191 }
1192
1193 const GPUTPCGMMergedTrack& track = mTracking->mIOPtrs.mergedTracks[i];
1194 for (uint32_t j = 0; j < track.NClusters(); j++) {
1195 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + j].num;
1196 if (GetMCLabelNID(hitId)) {
1197 int32_t mcID = GetMCLabelID(hitId, 0);
1198 if (mcID >= 0) {
1199 allowMCLabels[mcID] = true;
1200 }
1201 }
1202 }
1203 }
1204 for (uint32_t i = 0; i < GetNMCLabels(); i++) {
1205 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1206 int32_t mcID = GetMCLabelID(i, j);
1207 if (mcID >= 0 && allowMCLabels[mcID]) {
1208 mGoodHits[mNEvents - 1][i] = true;
1209 }
1210 }
1211 }
1212 }
1213 if (timer.IsRunning()) {
1214 GPUInfo("QA Time: Cluster attach status:\t\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1215 }
1216
1217 if (matchOnly) {
1218 return;
1219 }
1220
1221 // Recompute fNWeightCls (might have changed after merging events into timeframes)
1222 for (uint32_t iCol = 0; iCol < GetNMCCollissions(); iCol++) {
1223 for (uint32_t i = 0; i < GetNMCTracks(iCol); i++) {
1224 mMCParam[iCol][i].nWeightCls = 0.;
1225 }
1226 }
1227 for (uint32_t i = 0; i < GetNMCLabels(); i++) {
1228 float weightTotal = 0.f;
1229 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1230 if (GetMCLabelID(i, j) >= 0) {
1231 weightTotal += GetMCLabelWeight(i, j);
1232 }
1233 }
1234 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1235 if (GetMCLabelID(i, j) >= 0) {
1236 GetMCTrackObj(mMCParam, GetMCLabel(i, j)).nWeightCls += GetMCLabelWeight(i, j) / weightTotal;
1237 }
1238 }
1239 }
1240 if (timer.IsRunning()) {
1241 GPUInfo("QA Time: Compute cluster label weights:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1242 }
1243
1244 // Compute MC Track Parameters for MC Tracks
1245 tbb::parallel_for<uint32_t>(0, GetNMCCollissions(), [&](auto iCol) {
1246 for (uint32_t i = 0; i < GetNMCTracks(iCol); i++) {
1247 const mcInfo_t& info = GetMCTrack(i, iCol);
1248 additionalMCParameters& mc2 = mMCParam[iCol][i];
1249 mc2.pt = std::sqrt(info.pX * info.pX + info.pY * info.pY);
1250 mc2.phi = M_PI + std::atan2(-info.pY, -info.pX);
1251 float p = info.pX * info.pX + info.pY * info.pY + info.pZ * info.pZ;
1252 if (p < 1e-18) {
1253 mc2.theta = mc2.eta = 0.f;
1254 } else {
1255 mc2.theta = info.pZ == 0 ? (M_PI / 2) : (std::acos(info.pZ / std::sqrt(p)));
1256 mc2.eta = -std::log(std::tan(0.5 * mc2.theta));
1257 }
1258 if (mConfig.writeMCLabels) {
1259 std::vector<int32_t>& effBuffer = mcEffBuffer[mNEvents - 1];
1260 effBuffer[i] = mRecTracks[iCol][i] * 1000 + mFakeTracks[iCol][i];
1261 }
1262 } // clang-format off
1263 }, tbb::simple_partitioner()); // clang-format on
1264 if (timer.IsRunning()) {
1265 GPUInfo("QA Time: Compute track mc parameters:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1266 }
1267
1268 // Fill Efficiency Histograms
1269 if (mQATasks & taskTrackingEff) {
1270 for (uint32_t iCol = 0; iCol < GetNMCCollissions(); iCol++) {
1271 for (uint32_t i = 0; i < GetNMCTracks(iCol); i++) {
1272 if ((mMCTrackMin != -1 && (int32_t)i < mMCTrackMin) || (mMCTrackMax != -1 && (int32_t)i >= mMCTrackMax)) {
1273 continue;
1274 }
1275 const mcInfo_t& info = GetMCTrack(i, iCol);
1276 const additionalMCParameters& mc2 = mMCParam[iCol][i];
1277 if (mc2.nWeightCls == 0.f) {
1278 continue;
1279 }
1280 const float& mcpt = mc2.pt;
1281 const float& mcphi = mc2.phi;
1282 const float& mceta = mc2.eta;
1283
1284 if (info.primDaughters) {
1285 continue;
1286 }
1287 if (mc2.nWeightCls < mConfig.minNClEff) {
1288 continue;
1289 }
1290 int32_t findable = mc2.nWeightCls >= mConfig.minNClFindable;
1291 if (info.pid < 0) {
1292 continue;
1293 }
1294 if (info.charge == 0.f) {
1295 continue;
1296 }
1297 if (mConfig.filterCharge && info.charge * mConfig.filterCharge < 0) {
1298 continue;
1299 }
1300 if (mConfig.filterPID >= 0 && info.pid != mConfig.filterPID) {
1301 continue;
1302 }
1303
1304 if (fabsf(mceta) > ETA_MAX || mcpt < PT_MIN || mcpt > PT_MAX) {
1305 continue;
1306 }
1307
1308 float alpha = std::atan2(info.y, info.x);
1309 alpha /= M_PI / 9.f;
1310 alpha = std::floor(alpha);
1311 alpha *= M_PI / 9.f;
1312 alpha += M_PI / 18.f;
1313
1314 float c = std::cos(alpha);
1315 float s = std::sin(alpha);
1316 float localY = -info.x * s + info.y * c;
1317
1318 if (mConfig.dumpToROOTLevel >= 1) {
1319 static auto effdump = GPUROOTDump<TNtuple>::getNew("eff", "alpha:x:y:z:mcphi:mceta:mcpt:rec:fake:findable:prim:ncls");
1320 float localX = info.x * c + info.y * s;
1321 effdump.Fill(alpha, localX, localY, info.z, mcphi, mceta, mcpt, mRecTracks[iCol][i], mFakeTracks[iCol][i], findable, info.prim, mc2.nWeightCls);
1322 }
1323
1324 for (int32_t j = 0; j < 6; j++) {
1325 if (j == 3 || j == 4) {
1326 continue;
1327 }
1328 for (int32_t k = 0; k < 2; k++) {
1329 if (k == 0 && findable == 0) {
1330 continue;
1331 }
1332
1333 int32_t val = (j == 0) ? (mRecTracks[iCol][i] ? 1 : 0) : ((j == 1) ? (mRecTracks[iCol][i] ? mRecTracks[iCol][i] - 1 : 0) : ((j == 2) ? mFakeTracks[iCol][i] : 1));
1334 if (val == 0) {
1335 continue;
1336 }
1337
1338 for (int32_t l = 0; l < 5; l++) {
1339 if (info.prim && mcpt < PT_MIN_PRIM) {
1340 continue;
1341 }
1342 if (l != 3 && fabsf(mceta) > ETA_MAX2) {
1343 continue;
1344 }
1345 if (l < 4 && mcpt < 1.f / mConfig.qpt) {
1346 continue;
1347 }
1348
1349 float pos = l == 0 ? localY : (l == 1 ? info.z : (l == 2 ? mcphi : (l == 3 ? mceta : mcpt)));
1350
1351 mEff[j][k][!info.prim][l]->Fill(pos, val);
1352 }
1353 }
1354 }
1355 }
1356 }
1357 if (timer.IsRunning()) {
1358 GPUInfo("QA Time: Fill efficiency histograms:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1359 }
1360 }
1361
1362 // Fill Resolution Histograms
1363 if (mQATasks & (taskTrackingRes | taskTrackingResPull)) {
1364 GPUTPCGMPropagator prop;
1365 prop.SetMaxSinPhi(.999);
1366 prop.SetMaterialTPC();
1367 prop.SetPolynomialField(&mParam->polynomialField);
1368
1369 for (uint32_t i = 0; i < mTrackMCLabels.size(); i++) {
1370 if (mConfig.writeMCLabels) {
1371 std::vector<int32_t>& labelBuffer = mcLabelBuffer[mNEvents - 1];
1372 labelBuffer[i] = mTrackMCLabels[i].getTrackID();
1373 }
1374 if (mTrackMCLabels[i].isFake()) {
1375 continue;
1376 }
1377 const mcInfo_t& mc1 = GetMCTrack(mTrackMCLabels[i]);
1378 const additionalMCParameters& mc2 = GetMCTrackObj(mMCParam, mTrackMCLabels[i]);
1379
1380 if (mc1.primDaughters) {
1381 continue;
1382 }
1383 if (!tracksExternal) {
1384 if (!mTracking->mIOPtrs.mergedTracks[i].OK()) {
1385 continue;
1386 }
1387 if (mTracking->mIOPtrs.mergedTracks[i].MergedLooper()) {
1388 continue;
1389 }
1390 }
1391 if ((mMCTrackMin != -1 && mTrackMCLabels[i].getTrackID() < mMCTrackMin) || (mMCTrackMax != -1 && mTrackMCLabels[i].getTrackID() >= mMCTrackMax)) {
1392 continue;
1393 }
1394 if (fabsf(mc2.eta) > ETA_MAX || mc2.pt < PT_MIN || mc2.pt > PT_MAX) {
1395 continue;
1396 }
1397 if (mc1.charge == 0.f) {
1398 continue;
1399 }
1400 if (mc1.pid < 0) {
1401 continue;
1402 }
1403#ifndef GPUCA_RUN2
1404 if (mc1.t0 == -100.f) {
1405 continue;
1406 }
1407#endif
1408 if (mConfig.filterCharge && mc1.charge * mConfig.filterCharge < 0) {
1409 continue;
1410 }
1411 if (mConfig.filterPID >= 0 && mc1.pid != mConfig.filterPID) {
1412 continue;
1413 }
1414 if (mc2.nWeightCls < mConfig.minNClRes) {
1415 continue;
1416 }
1417 if (mConfig.resPrimaries == 1 && !mc1.prim) {
1418 continue;
1419 } else if (mConfig.resPrimaries == 2 && mc1.prim) {
1420 continue;
1421 }
1422 if (GetMCTrackObj(mTrackMCLabelsReverse, mTrackMCLabels[i]) != (int32_t)i) {
1423 continue;
1424 }
1425
1427 float alpha = 0.f;
1428 int32_t side;
1429 if (tracksExternal) {
1430#ifndef GPUCA_STANDALONE
1431 for (int32_t k = 0; k < 5; k++) {
1432 param.Par()[k] = (*tracksExternal)[i].getParams()[k];
1433 }
1434 for (int32_t k = 0; k < 15; k++) {
1435 param.Cov()[k] = (*tracksExternal)[i].getCov()[k];
1436 }
1437 param.X() = (*tracksExternal)[i].getX();
1438 param.TOffset() = (*tracksExternal)[i].getTime0();
1439 alpha = (*tracksExternal)[i].getAlpha();
1440 side = (*tracksExternal)[i].hasBothSidesClusters() ? 2 : ((*tracksExternal)[i].hasCSideClusters() ? 1 : 0);
1441#endif
1442 } else {
1443 param = mTracking->mIOPtrs.mergedTracks[i].GetParam();
1444 alpha = mTracking->mIOPtrs.mergedTracks[i].GetAlpha();
1445 side = mTracking->mIOPtrs.mergedTracks[i].CCE() ? 2 : (mTracking->mIOPtrs.mergedTracks[i].CSide() ? 1 : 0);
1446 }
1447
1448 float mclocal[4]; // Rotated x,y,Px,Py mc-coordinates - the MC data should be rotated since the track is propagated best along x
1449 float c = std::cos(alpha);
1450 float s = std::sin(alpha);
1451 float x = mc1.x;
1452 float y = mc1.y;
1453 mclocal[0] = x * c + y * s;
1454 mclocal[1] = -x * s + y * c;
1455 float px = mc1.pX;
1456 float py = mc1.pY;
1457 mclocal[2] = px * c + py * s;
1458 mclocal[3] = -px * s + py * c;
1459
1460 if (mclocal[0] < TRACK_EXPECTED_REFERENCE_X - 3) {
1461 continue;
1462 }
1463 if (mclocal[0] > param.GetX() + 20) {
1464 continue;
1465 }
1466 if (param.GetX() > mConfig.maxResX) {
1467 continue;
1468 }
1469
1470 auto getdz = [this, &param, &mc1, &side, tracksExternal]() {
1471 if (tracksExternal) {
1472 return param.GetZ();
1473 }
1474 if (!mParam->continuousMaxTimeBin) {
1475 return param.GetZ() - mc1.z;
1476 }
1477 float shift = 0;
1478#ifndef GPUCA_RUN2
1479 if (side != 2) {
1480 shift = mTracking->GetTPCTransform()->convDeltaTimeToDeltaZinTimeFrame(side * GPUChainTracking::NSECTORS / 2, param.GetTOffset() - mc1.t0);
1481 }
1482#endif
1483 return param.GetZ() + shift - mc1.z;
1484 };
1485
1486 prop.SetTrack(&param, alpha);
1487 bool inFlyDirection = 0;
1488 if (mConfig.strict) {
1489 const float dx = param.X() - std::max<float>(mclocal[0], TRACK_EXPECTED_REFERENCE_X_DEFAULT); // Limit distance check
1490 const float dy = param.Y() - mclocal[1];
1491 const float dz = getdz();
1492 if (dx * dx + dy * dy + dz * dz > 5.f * 5.f) {
1493 continue;
1494 }
1495 }
1496
1497 if (prop.PropagateToXAlpha(mclocal[0], alpha, inFlyDirection)) {
1498 continue;
1499 }
1500 if (fabsf(param.Y() - mclocal[1]) > (mConfig.strict ? 1.f : 4.f) || fabsf(getdz()) > (mConfig.strict ? 1.f : 4.f)) {
1501 continue;
1502 }
1503 float charge = mc1.charge > 0 ? 1.f : -1.f;
1504
1505 float deltaY = param.GetY() - mclocal[1];
1506 float deltaZ = getdz();
1507 float deltaPhiNative = param.GetSinPhi() - mclocal[3] / mc2.pt;
1508 float deltaPhi = std::asin(param.GetSinPhi()) - std::atan2(mclocal[3], mclocal[2]);
1509 float deltaLambdaNative = param.GetDzDs() - mc1.pZ / mc2.pt;
1510 float deltaLambda = std::atan(param.GetDzDs()) - std::atan2(mc1.pZ, mc2.pt);
1511 float deltaPtNative = (param.GetQPt() - charge / mc2.pt) * charge;
1512 float deltaPt = (fabsf(1.f / param.GetQPt()) - mc2.pt) / mc2.pt;
1513
1514 float paramval[5] = {mclocal[1], mc1.z, mc2.phi, mc2.eta, mc2.pt};
1515 float resval[5] = {deltaY, deltaZ, mConfig.nativeFitResolutions ? deltaPhiNative : deltaPhi, mConfig.nativeFitResolutions ? deltaLambdaNative : deltaLambda, mConfig.nativeFitResolutions ? deltaPtNative : deltaPt};
1516 float pullval[5] = {deltaY / std::sqrt(param.GetErr2Y()), deltaZ / std::sqrt(param.GetErr2Z()), deltaPhiNative / std::sqrt(param.GetErr2SinPhi()), deltaLambdaNative / std::sqrt(param.GetErr2DzDs()), deltaPtNative / std::sqrt(param.GetErr2QPt())};
1517
1518 for (int32_t j = 0; j < 5; j++) {
1519 for (int32_t k = 0; k < 5; k++) {
1520 if (k != 3 && fabsf(mc2.eta) > ETA_MAX2) {
1521 continue;
1522 }
1523 if (k < 4 && mc2.pt < 1.f / mConfig.qpt) {
1524 continue;
1525 }
1526 if (mQATasks & taskTrackingRes) {
1527 mRes2[j][k]->Fill(resval[j], paramval[k]);
1528 }
1529 if (mQATasks & taskTrackingResPull) {
1530 mPull2[j][k]->Fill(pullval[j], paramval[k]);
1531 }
1532 }
1533 }
1534 }
1535 if (timer.IsRunning()) {
1536 GPUInfo("QA Time: Fill resolution histograms:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1537 }
1538 }
1539
1540 if ((mQATasks & taskClusterAttach) && !tracksExternal) {
1541 // Fill cluster histograms
1542 for (uint32_t iTrk = 0; iTrk < nReconstructedTracks; iTrk++) {
1543 const GPUTPCGMMergedTrack& track = mTracking->mIOPtrs.mergedTracks[iTrk];
1544 if (!track.OK()) {
1545 continue;
1546 }
1547 if (!mTrackMCLabels[iTrk].isValid()) {
1548 for (uint32_t k = 0; k < track.NClusters(); k++) {
1549 if (mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].state & GPUTPCGMMergedTrackHit::flagReject) {
1550 continue;
1551 }
1552 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].num;
1553 float totalWeight = 0.;
1554 for (int32_t j = 0; j < GetMCLabelNID(hitId); j++) {
1555 if (GetMCLabelID(hitId, j) >= 0 && GetMCTrackObj(mMCParam, GetMCLabel(hitId, j)).pt > 1.f / mTracking->GetParam().rec.maxTrackQPtB5) {
1556 totalWeight += GetMCLabelWeight(hitId, j);
1557 }
1558 }
1559 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[hitId];
1560 const auto& r = checkClusterState<false>(attach);
1561 if (totalWeight > 0) {
1562 float weight = 1.f / (totalWeight * (mClusterParam[hitId].attached + mClusterParam[hitId].fakeAttached));
1563 for (int32_t j = 0; j < GetMCLabelNID(hitId); j++) {
1564 mcLabelI_t label = GetMCLabel(hitId, j);
1565 if (!label.isFake() && GetMCTrackObj(mMCParam, label).pt > 1.f / mTracking->GetParam().rec.maxTrackQPtB5) {
1566 float pt = GetMCTrackObj(mMCParam, label).pt;
1567 if (pt < PT_MIN_CLUST) {
1568 pt = PT_MIN_CLUST;
1569 }
1570 mClusters[CL_fake]->Fill(pt, GetMCLabelWeight(hitId, j) * weight);
1571 mClusters[CL_att_adj]->Fill(pt, GetMCLabelWeight(hitId, j) * weight);
1572 if (GetMCTrackObj(mRecTracks, label)) {
1573 mClusters[CL_tracks]->Fill(pt, GetMCLabelWeight(hitId, j) * weight);
1574 }
1575 mClusters[CL_all]->Fill(pt, GetMCLabelWeight(hitId, j) * weight);
1576 if (r.protect || r.physics) {
1577 mClusters[CL_prot]->Fill(pt, GetMCLabelWeight(hitId, j) * weight);
1578 }
1579 if (r.physics) {
1580 mClusters[CL_physics]->Fill(pt, GetMCLabelWeight(hitId, j) * weight);
1581 }
1582 }
1583 }
1584 } else {
1585 float weight = 1.f / (mClusterParam[hitId].attached + mClusterParam[hitId].fakeAttached);
1586 mClusters[CL_fake]->Fill(0.f, weight);
1587 mClusters[CL_att_adj]->Fill(0.f, weight);
1588 mClusters[CL_all]->Fill(0.f, weight);
1589 mClusterCounts.nUnaccessible += weight;
1590 if (r.protect || r.physics) {
1591 mClusters[CL_prot]->Fill(0.f, weight);
1592 }
1593 if (r.physics) {
1594 mClusters[CL_physics]->Fill(0.f, weight);
1595 }
1596 }
1597 }
1598 continue;
1599 }
1600 mcLabelI_t label = mTrackMCLabels[iTrk];
1601 if (mMCTrackMin != -1 && (label.getTrackID() < mMCTrackMin || label.getTrackID() >= mMCTrackMax)) {
1602 continue;
1603 }
1604 for (uint32_t k = 0; k < track.NClusters(); k++) {
1605 if (mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].state & GPUTPCGMMergedTrackHit::flagReject) {
1606 continue;
1607 }
1608 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].num;
1609 float pt = GetMCTrackObj(mMCParam, label).pt;
1610 if (pt < PT_MIN_CLUST) {
1611 pt = PT_MIN_CLUST;
1612 }
1613 float weight = 1.f / (mClusterParam[hitId].attached + mClusterParam[hitId].fakeAttached);
1614 bool correct = false;
1615 for (int32_t j = 0; j < GetMCLabelNID(hitId); j++) {
1616 if (label == GetMCLabel(hitId, j)) {
1617 correct = true;
1618 break;
1619 }
1620 }
1621 if (correct) {
1622 mClusters[CL_attached]->Fill(pt, weight);
1623 mClusters[CL_tracks]->Fill(pt, weight);
1624 } else {
1625 mClusters[CL_fake]->Fill(pt, weight);
1626 }
1627 mClusters[CL_att_adj]->Fill(pt, weight);
1628 mClusters[CL_all]->Fill(pt, weight);
1629 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[hitId];
1630 const auto& r = checkClusterState<false>(attach);
1631 if (r.protect || r.physics) {
1632 mClusters[CL_prot]->Fill(pt, weight);
1633 }
1634 if (r.physics) {
1635 mClusters[CL_physics]->Fill(pt, weight);
1636 }
1637 }
1638 }
1639 for (uint32_t i = 0; i < GetNMCLabels(); i++) {
1640 if ((mMCTrackMin != -1 && GetMCLabelID(i, 0) < mMCTrackMin) || (mMCTrackMax != -1 && GetMCLabelID(i, 0) >= mMCTrackMax)) {
1641 continue;
1642 }
1643 if (mClusterParam[i].attached || mClusterParam[i].fakeAttached) {
1644 continue;
1645 }
1646 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[i];
1647 const auto& r = checkClusterState<false>(attach);
1648 if (mClusterParam[i].adjacent) {
1649 int32_t label = mTracking->mIOPtrs.mergedTrackHitAttachment[i] & gputpcgmmergertypes::attachTrackMask;
1650 if (!mTrackMCLabels[label].isValid()) {
1651 float totalWeight = 0.;
1652 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1653 mcLabelI_t labelT = GetMCLabel(i, j);
1654 if (!labelT.isFake() && GetMCTrackObj(mMCParam, labelT).pt > 1.f / mTracking->GetParam().rec.maxTrackQPtB5) {
1655 totalWeight += GetMCLabelWeight(i, j);
1656 }
1657 }
1658 float weight = 1.f / totalWeight;
1659 if (totalWeight > 0) {
1660 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1661 mcLabelI_t labelT = GetMCLabel(i, j);
1662 if (!labelT.isFake() && GetMCTrackObj(mMCParam, labelT).pt > 1.f / mTracking->GetParam().rec.maxTrackQPtB5) {
1663 float pt = GetMCTrackObj(mMCParam, labelT).pt;
1664 if (pt < PT_MIN_CLUST) {
1665 pt = PT_MIN_CLUST;
1666 }
1667 if (GetMCTrackObj(mRecTracks, labelT)) {
1668 mClusters[CL_tracks]->Fill(pt, GetMCLabelWeight(i, j) * weight);
1669 }
1670 mClusters[CL_att_adj]->Fill(pt, GetMCLabelWeight(i, j) * weight);
1671 mClusters[CL_fakeAdj]->Fill(pt, GetMCLabelWeight(i, j) * weight);
1672 mClusters[CL_all]->Fill(pt, GetMCLabelWeight(i, j) * weight);
1673 if (r.protect || r.physics) {
1674 mClusters[CL_prot]->Fill(pt, GetMCLabelWeight(i, j) * weight);
1675 }
1676 if (r.physics) {
1677 mClusters[CL_physics]->Fill(pt, GetMCLabelWeight(i, j) * weight);
1678 }
1679 }
1680 }
1681 } else {
1682 mClusters[CL_att_adj]->Fill(0.f, 1.f);
1683 mClusters[CL_fakeAdj]->Fill(0.f, 1.f);
1684 mClusters[CL_all]->Fill(0.f, 1.f);
1685 mClusterCounts.nUnaccessible++;
1686 if (r.protect || r.physics) {
1687 mClusters[CL_prot]->Fill(0.f, 1.f);
1688 }
1689 if (r.physics) {
1690 mClusters[CL_physics]->Fill(0.f, 1.f);
1691 }
1692 }
1693 } else {
1694 float pt = GetMCTrackObj(mMCParam, mTrackMCLabels[label]).pt;
1695 if (pt < PT_MIN_CLUST) {
1696 pt = PT_MIN_CLUST;
1697 }
1698 mClusters[CL_att_adj]->Fill(pt, 1.f);
1699 mClusters[CL_tracks]->Fill(pt, 1.f);
1700 mClusters[CL_all]->Fill(pt, 1.f);
1701 if (r.protect || r.physics) {
1702 mClusters[CL_prot]->Fill(pt, 1.f);
1703 }
1704 if (r.physics) {
1705 mClusters[CL_physics]->Fill(pt, 1.f);
1706 }
1707 }
1708 } else {
1709 float totalWeight = 0.;
1710 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1711 mcLabelI_t labelT = GetMCLabel(i, j);
1712 if (!labelT.isFake() && GetMCTrackObj(mMCParam, labelT).pt > 1.f / mTracking->GetParam().rec.maxTrackQPtB5) {
1713 totalWeight += GetMCLabelWeight(i, j);
1714 }
1715 }
1716 if (totalWeight > 0) {
1717 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1718 mcLabelI_t label = GetMCLabel(i, j);
1719 if (!label.isFake() && GetMCTrackObj(mMCParam, label).pt > 1.f / mTracking->GetParam().rec.maxTrackQPtB5) {
1720 float pt = GetMCTrackObj(mMCParam, label).pt;
1721 if (pt < PT_MIN_CLUST) {
1722 pt = PT_MIN_CLUST;
1723 }
1724 float weight = GetMCLabelWeight(i, j) / totalWeight;
1725 if (mClusterParam[i].fakeAdjacent) {
1726 mClusters[CL_fakeAdj]->Fill(pt, weight);
1727 }
1728 if (mClusterParam[i].fakeAdjacent) {
1729 mClusters[CL_att_adj]->Fill(pt, weight);
1730 }
1731 if (GetMCTrackObj(mRecTracks, label)) {
1732 mClusters[CL_tracks]->Fill(pt, weight);
1733 }
1734 mClusters[CL_all]->Fill(pt, weight);
1735 if (r.protect || r.physics) {
1736 mClusters[CL_prot]->Fill(pt, weight);
1737 }
1738 if (r.physics) {
1739 mClusters[CL_physics]->Fill(pt, weight);
1740 }
1741 }
1742 }
1743 } else {
1744 if (mClusterParam[i].fakeAdjacent) {
1745 mClusters[CL_fakeAdj]->Fill(0.f, 1.f);
1746 }
1747 if (mClusterParam[i].fakeAdjacent) {
1748 mClusters[CL_att_adj]->Fill(0.f, 1.f);
1749 }
1750 mClusters[CL_all]->Fill(0.f, 1.f);
1751 mClusterCounts.nUnaccessible++;
1752 if (r.protect || r.physics) {
1753 mClusters[CL_prot]->Fill(0.f, 1.f);
1754 }
1755 if (r.physics) {
1756 mClusters[CL_physics]->Fill(0.f, 1.f);
1757 }
1758 }
1759 }
1760 }
1761
1762 if (timer.IsRunning()) {
1763 GPUInfo("QA Time: Fill cluster histograms:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1764 }
1765 }
1766 } else if (!mConfig.inputHistogramsOnly && !mConfig.noMC && (mQATasks & (taskTrackingEff | taskTrackingRes | taskTrackingResPull | taskClusterAttach))) {
1767 GPUWarning("No MC information available, only running partial TPC QA!");
1768 } // mcAvail
1769
1770 if ((mQATasks & taskTrackStatistics) && !tracksExternal) {
1771 // Fill track statistic histograms
1772 std::vector<std::array<float, 3>> clusterAttachCounts;
1773 if (mcAvail) {
1774 clusterAttachCounts.resize(GetNMCLabels(), {0.f, 0.f});
1775 }
1776 for (uint32_t i = 0; i < nReconstructedTracks; i++) {
1777 const GPUTPCGMMergedTrack& track = mTracking->mIOPtrs.mergedTracks[i];
1778 if (!track.OK()) {
1779 continue;
1780 }
1781 mTrackPt->Fill(1.f / fabsf(track.GetParam().GetQPt()));
1782 if (mParam->par.dodEdx && mParam->dodEdxEnabled && track.NClusters() >= 60) {
1783 const GPUdEdxInfo& trackdEdx = mTracking->GetProcessors()->tpcMerger.MergedTracksdEdx()[i];
1784 const float logp = logf(1.f / fabsf(track.GetParam().GetQPt()) * sqrtf(1.f + track.GetParam().GetDzDs() * track.GetParam().GetDzDs()));
1785 mTrackdEdx[0]->Fill(logp, trackdEdx.dEdxTotTPC);
1786 mTrackdEdx[1]->Fill(logp, trackdEdx.dEdxMaxTPC);
1787 }
1788 mNCl[0]->Fill(track.NClustersFitted());
1789 uint32_t nClCorrected = 0;
1790 const auto& trackClusters = mTracking->mIOPtrs.mergedTrackHits;
1791 uint32_t jNext = 0;
1792 for (uint32_t j = 0; j < track.NClusters(); j = jNext) {
1793 uint32_t rowClCount = !(trackClusters[track.FirstClusterRef() + j].state & GPUTPCGMMergedTrackHit::flagReject);
1794 for (jNext = j + 1; j < track.NClusters(); jNext++) {
1795 if (trackClusters[track.FirstClusterRef() + j].sector != trackClusters[track.FirstClusterRef() + jNext].sector || trackClusters[track.FirstClusterRef() + j].row != trackClusters[track.FirstClusterRef() + jNext].row) {
1796 break;
1797 }
1798 rowClCount += !(trackClusters[track.FirstClusterRef() + jNext].state & GPUTPCGMMergedTrackHit::flagReject);
1799 }
1800 if (!track.MergedLooper() && rowClCount) {
1801 nClCorrected++;
1802 }
1803 if (mcAvail && rowClCount) {
1804 for (uint32_t k = j; k < jNext; k++) {
1805 const auto& cl = trackClusters[track.FirstClusterRef() + k];
1806 if (cl.state & GPUTPCGMMergedTrackHit::flagReject) {
1807 continue;
1808 }
1809 bool labelOk = false, labelOkNonFake = false;
1810 const mcLabelI_t& trkLabel = mTrackMCLabels[i];
1811 if (trkLabel.isValid() && !trkLabel.isNoise()) {
1812 for (int32_t l = 0; l < GetMCLabelNID(cl.num); l++) {
1813 const mcLabelI_t& clLabel = GetMCLabel(cl.num, l);
1814 if (clLabel.isValid() && !clLabel.isNoise() && CompareIgnoreFake(trkLabel, clLabel)) {
1815 labelOk = true;
1816 if (!trkLabel.isFake()) {
1817 labelOkNonFake = true;
1818 }
1819 break;
1820 }
1821 }
1822 }
1823 clusterAttachCounts[cl.num][0] += 1.0f;
1824 clusterAttachCounts[cl.num][1] += (float)labelOk / rowClCount;
1825 clusterAttachCounts[cl.num][2] += (float)labelOkNonFake / rowClCount;
1826 }
1827 }
1828 }
1829 if (nClCorrected) {
1830 mNCl[1]->Fill(nClCorrected);
1831 }
1832 mT0[0]->Fill(track.GetParam().GetTOffset());
1833#ifndef GPUCA_RUN2
1834 if (mTrackMCLabels.size() && !mTrackMCLabels[i].isFake() && !track.MergedLooper() && !track.CCE()) {
1835 const auto& info = GetMCTrack(mTrackMCLabels[i]);
1836 if (info.t0 != -100.f) {
1837 mT0[1]->Fill(track.GetParam().GetTOffset() - info.t0);
1838 }
1839 }
1840#endif
1841 }
1842 if (mClNative && mTracking && mTracking->GetTPCTransform()) {
1843 for (uint32_t i = 0; i < GPUChainTracking::NSECTORS; i++) {
1844 for (uint32_t j = 0; j < GPUTPCGeometry::NROWS; j++) {
1845 for (uint32_t k = 0; k < mClNative->nClusters[i][j]; k++) {
1846 const auto& cl = mClNative->clusters[i][j][k];
1847 float x, y, z;
1848 GPUTPCConvertImpl::convert(*mTracking->GetTPCTransform(), mTracking->GetParam(), i, j, cl.getPad(), cl.getTime(), x, y, z);
1849 mTracking->GetParam().Sector2Global(i, x, y, z, &x, &y, &z);
1850 mClXY->Fill(x, y);
1851 }
1852 }
1853 }
1854 }
1855 if (mcAvail) {
1856 double clusterAttachNormalizedCount = 0, clusterAttachNormalizedCountNonFake = 0;
1857 for (uint32_t i = 0; i < clusterAttachCounts.size(); i++) {
1858 if (clusterAttachCounts[i][0]) {
1859 clusterAttachNormalizedCount += clusterAttachCounts[i][1] / clusterAttachCounts[i][0];
1860 clusterAttachNormalizedCountNonFake += clusterAttachCounts[i][2] / clusterAttachCounts[i][0];
1861 }
1862 }
1863 mClusterCounts.nCorrectlyAttachedNormalized = clusterAttachNormalizedCount;
1864 mClusterCounts.nCorrectlyAttachedNormalizedNonFake = clusterAttachNormalizedCountNonFake;
1865 clusterAttachCounts.clear();
1866 }
1867
1868 if (timer.IsRunning()) {
1869 GPUInfo("QA Time: Fill track statistics:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1870 }
1871 }
1872
1873 uint32_t nCl = clNative ? clNative->nClustersTotal : mTracking->GetProcessors()->tpcMerger.NMaxClusters();
1874 mClusterCounts.nTotal += nCl;
1875 if (mQATasks & (taskClusterCounts | taskClusterRejection)) {
1876 for (uint32_t iSector = 0; iSector < GPUTPCGeometry::NSECTORS; iSector++) {
1877 for (uint32_t iRow = 0; iRow < GPUTPCGeometry::NROWS; iRow++) {
1878 for (uint32_t iCl = 0; iCl < clNative->nClusters[iSector][iRow]; iCl++) {
1879 uint32_t i = clNative->clusterOffset[iSector][iRow] + iCl;
1880 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[i];
1881 const auto& r = checkClusterState<true>(attach, &mClusterCounts);
1882
1883 if (mQATasks & taskClusterRejection) {
1884 if (mcAvail) {
1885 float totalWeight = 0, weight400 = 0, weight40 = 0;
1886 for (int32_t j = 0; j < GetMCLabelNID(i); j++) {
1887 const auto& label = GetMCLabel(i, j);
1888 if (GetMCLabelID(label) >= 0) {
1889 totalWeight += GetMCLabelWeight(label);
1890 if (GetMCTrackObj(mMCParam, label).pt >= 0.4) {
1891 weight400 += GetMCLabelWeight(label);
1892 }
1893 if (GetMCTrackObj(mMCParam, label).pt <= 0.04) {
1894 weight40 += GetMCLabelWeight(label);
1895 }
1896 }
1897 }
1898 if (totalWeight > 0 && 10.f * weight400 >= totalWeight) {
1899 if (!r.unattached && !r.protect && !r.physics) {
1900 mClusterCounts.nFakeRemove400++;
1901 int32_t totalFake = weight400 < 0.9f * totalWeight;
1902 if (totalFake) {
1903 mClusterCounts.nFullFakeRemove400++;
1904 }
1905 /*printf("Fake removal (%d): Hit %7d, attached %d lowPt %d looper %d tube200 %d highIncl %d tube %d bad %d recPt %7.2f recLabel %6d", totalFake, i, (int32_t) (mClusterParam[i].attached || mClusterParam[i].fakeAttached),
1906 (int32_t) lowPt, (int32_t) ((attach & gputpcgmmergertypes::attachGoodLeg) == 0), (int32_t) ((attach & gputpcgmmergertypes::attachTube) && mev200),
1907 (int32_t) ((attach & gputpcgmmergertypes::attachHighIncl) != 0), (int32_t) ((attach & gputpcgmmergertypes::attachTube) != 0), (int32_t) ((attach & gputpcgmmergertypes::attachGood) == 0),
1908 fabsf(qpt) > 0 ? 1.f / qpt : 0.f, id);
1909 for (int32_t j = 0;j < GetMCLabelNID(i);j++)
1910 {
1911 //if (GetMCLabelID(i, j) < 0) break;
1912 printf(" - label%d %6d weight %5d", j, GetMCLabelID(i, j), (int32_t) GetMCLabelWeight(i, j));
1913 if (GetMCLabelID(i, j) >= 0) printf(" - pt %7.2f", mMCParam[GetMCLabelID(i, j)].pt);
1914 else printf(" ");
1915 }
1916 printf("\n");*/
1917 }
1918 mClusterCounts.nAbove400++;
1919 }
1920 if (totalWeight > 0 && weight40 >= 0.9 * totalWeight) {
1921 mClusterCounts.nBelow40++;
1922 if (r.protect || r.physics) {
1923 mClusterCounts.nFakeProtect40++;
1924 }
1925 }
1926 }
1927
1928 if (r.physics) {
1929 mClusterCounts.nPhysics++;
1930 }
1931 if (r.protect) {
1932 mClusterCounts.nProt++;
1933 }
1934 if (r.unattached) {
1935 mClusterCounts.nUnattached++;
1936 }
1937 }
1938 if (mQATasks & taskClusterRejection) {
1939 if (mTracking && clNative) {
1940 const auto& cl = clNative->clustersLinear[i];
1941 mClRej[0]->Fill(cl.getPad() - GPUTPCGeometry::NPads(iRow) / 2 + 0.5, iRow, 1.f);
1942 if (!r.unattached && !r.protect) {
1943 mClRej[1]->Fill(cl.getPad() - GPUTPCGeometry::NPads(iRow) / 2 + 0.5, iRow, 1.f);
1944 }
1945 }
1946 }
1947 }
1948 }
1949 }
1950 }
1951
1952 // Process cluster count statistics
1953 if ((mQATasks & taskClusterCounts) && mConfig.clusterRejectionHistograms) {
1954 DoClusterCounts(nullptr);
1955 mClusterCounts = counts_t();
1956 }
1957
1958 if (timer.IsRunning()) {
1959 GPUInfo("QA Time: Cluster Counts:\t%6.0f us", timer.GetCurrentElapsedTime(true) * 1e6);
1960 }
1961
1962 if (mConfig.dumpToROOTLevel >= 1 && !tracksExternal) {
1963 if (!clNative || !mTracking || !mTracking->mIOPtrs.mergedTrackHitAttachment || !mTracking->mIOPtrs.mergedTracks) {
1964 throw std::runtime_error("Cannot dump non o2::tpc::clusterNative clusters, need also hit attachmend and GPU tracks");
1965 }
1966 uint32_t clid = 0;
1967 for (uint32_t i = 0; i < GPUChainTracking::NSECTORS; i++) {
1968 for (uint32_t j = 0; j < GPUTPCGeometry::NROWS; j++) {
1969 for (uint32_t k = 0; k < mClNative->nClusters[i][j]; k++) {
1970 const auto& cl = mClNative->clusters[i][j][k];
1971 uint32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[clid];
1972 float x = 0, y = 0, z = 0;
1975 const auto& trk = mTracking->mIOPtrs.mergedTracks[track];
1976 mTracking->GetTPCTransform()->Transform(i, j, cl.getPad(), cl.getTime(), x, y, z, trk.GetParam().GetTOffset());
1977 mTracking->GetParam().Sector2Global(i, x, y, z, &x, &y, &z);
1978 }
1979 uint32_t extState = mTracking->mIOPtrs.mergedTrackHitStates ? mTracking->mIOPtrs.mergedTrackHitStates[clid] : 0;
1980
1981 if (mConfig.dumpToROOTLevel >= 2) {
1984 memset((void*)&trk, 0, sizeof(trk));
1985 memset((void*)&trkHit, 0, sizeof(trkHit));
1988 trk = mTracking->mIOPtrs.mergedTracks[track];
1989 for (uint32_t l = 0; l < trk.NClusters(); l++) {
1990 const auto& tmp = mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef() + l];
1991 if (tmp.num == clid) {
1992 trkHit = tmp;
1993 break;
1994 }
1995 }
1996 }
1997 static auto cldump = GPUROOTDump<o2::tpc::ClusterNative, GPUTPCGMMergedTrack, GPUTPCGMMergedTrackHit, uint32_t, uint32_t, float, float, float, uint32_t, uint32_t, uint32_t>::getNew("cluster", "track", "trackHit", "attach", "extState", "x", "y", "z", "sector", "row", "nEv", "clusterTree");
1998 cldump.Fill(cl, trk, trkHit, attach, extState, x, y, z, i, j, mNEvents - 1);
1999 } else {
2000 static auto cldump = GPUROOTDump<o2::tpc::ClusterNative, uint32_t, uint32_t, float, float, float, uint32_t, uint32_t, uint32_t>::getNew("cluster", "attach", "extState", "x", "y", "z", "sector", "row", "nEv", "clusterTree");
2001 cldump.Fill(cl, attach, extState, x, y, z, i, j, mNEvents - 1);
2002 }
2003 clid++;
2004 }
2005 }
2006 }
2007
2008 static auto trkdump = GPUROOTDump<uint32_t, GPUTPCGMMergedTrack>::getNew("nEv", "track", "tracksTree");
2009 for (uint32_t i = 0; i < mTracking->mIOPtrs.nMergedTracks; i++) {
2010 if (mTracking->mIOPtrs.mergedTracks[i].OK()) {
2011 trkdump.Fill(mNEvents - 1, mTracking->mIOPtrs.mergedTracks[i]);
2012 }
2013 }
2014
2015 if (mTracking && mTracking->GetProcessingSettings().createO2Output) {
2016 static auto o2trkdump = GPUROOTDump<uint32_t, o2::tpc::TrackTPC>::getNew("nEv", "track", "tracksO2Tree");
2017 for (uint32_t i = 0; i < mTracking->mIOPtrs.nOutputTracksTPCO2; i++) {
2018 o2trkdump.Fill(mNEvents - 1, mTracking->mIOPtrs.outputTracksTPCO2[i]);
2019 }
2020 }
2021 }
2022
2023 if (mConfig.compareTrackStatus) {
2024#ifdef GPUCA_DETERMINISTIC_MODE
2025 if (!mTracking || !mTracking->GetProcessingSettings().deterministicGPUReconstruction)
2026#endif
2027 {
2028 throw std::runtime_error("Need deterministic processing to compare track status");
2029 }
2030 std::vector<uint8_t> status(mTracking->mIOPtrs.nMergedTracks);
2031 for (uint32_t i = 0; i < mTracking->mIOPtrs.nMergedTracks; i++) {
2032 const auto& trk = mTracking->mIOPtrs.mergedTracks[i];
2033 status[i] = trk.OK() && trk.NClusters() && trk.GetParam().GetNDF() > 0 && (mConfig.noMC || (mTrackMCLabels[i].isValid() && !mTrackMCLabels[i].isFake()));
2034 }
2035 if (mConfig.compareTrackStatus == 1) {
2036 std::ofstream("track.status", std::ios::binary).write((char*)status.data(), status.size() * sizeof(status[0]));
2037 } else if (mConfig.compareTrackStatus == 2) {
2038 std::ifstream f("track.status", std::ios::binary | std::ios::ate);
2039 std::vector<uint8_t> comp(f.tellg());
2040 f.seekg(0);
2041 f.read((char*)comp.data(), comp.size());
2042
2043 if (comp.size() != status.size()) {
2044 throw std::runtime_error("Number of tracks candidates in track fit in track.status and in current reconstruction differ");
2045 }
2046 std::vector<uint32_t> missing, missingComp;
2047 for (uint32_t i = 0; i < status.size(); i++) {
2048 if (status[i] && !comp[i]) {
2049 missingComp.emplace_back(i);
2050 }
2051 if (comp[i] && !status[i]) {
2052 missing.emplace_back(i);
2053 }
2054 }
2055 auto printer = [](std::vector<uint32_t> m, const char* name) {
2056 if (m.size()) {
2057 printf("Missing in %s reconstruction: (%zu)\n", name, m.size());
2058 for (uint32_t i = 0; i < m.size(); i++) {
2059 if (i) {
2060 printf(", ");
2061 }
2062 printf("%d", m[i]);
2063 }
2064 printf("\n");
2065 }
2066 };
2067 printer(missing, "current");
2068 printer(missingComp, "comparison");
2069 }
2070 }
2071
2072 mTrackingScratchBuffer.clear();
2073 mTrackingScratchBuffer.shrink_to_fit();
2074}
2075
2076void GPUQA::GetName(char* fname, int32_t k, bool noDash)
2077{
2078 const int32_t nNewInput = mConfig.inputHistogramsOnly ? 0 : 1;
2079 if (k || mConfig.inputHistogramsOnly || mConfig.name.size()) {
2080 if (!(mConfig.inputHistogramsOnly || k)) {
2081 snprintf(fname, 1024, "%s%s", mConfig.name.c_str(), noDash ? "" : " - ");
2082 } else if (mConfig.compareInputNames.size() > (unsigned)(k - nNewInput)) {
2083 snprintf(fname, 1024, "%s%s", mConfig.compareInputNames[k - nNewInput].c_str(), noDash ? "" : " - ");
2084 } else {
2085 strcpy(fname, mConfig.compareInputs[k - nNewInput].c_str());
2086 if (strlen(fname) > 5 && strcmp(fname + strlen(fname) - 5, ".root") == 0) {
2087 fname[strlen(fname) - 5] = 0;
2088 }
2089 if (!noDash) {
2090 strcat(fname, " - ");
2091 }
2092 }
2093 } else {
2094 fname[0] = 0;
2095 }
2096}
2097
2098template <class T>
2099T* GPUQA::GetHist(T*& ee, std::vector<std::unique_ptr<TFile>>& tin, int32_t k, int32_t nNewInput)
2100{
2101 T* e = ee;
2102 if ((mConfig.inputHistogramsOnly || k) && (e = dynamic_cast<T*>(tin[k - nNewInput]->Get(e->GetName()))) == nullptr) {
2103 GPUWarning("Missing histogram in input %s: %s", mConfig.compareInputs[k - nNewInput].c_str(), ee->GetName());
2104 return (nullptr);
2105 }
2106 ee = e;
2107 return (e);
2108}
2109
2110void GPUQA::DrawQAHistogramsCleanup()
2111{
2112 clearGarbagageCollector();
2113}
2114
2115void GPUQA::resetHists()
2116{
2117 if (!mQAInitialized) {
2118 throw std::runtime_error("QA not initialized");
2119 }
2120 if (mHaveExternalHists) {
2121 throw std::runtime_error("Cannot reset external hists");
2122 }
2123 for (auto& h : *mHist1D) {
2124 h.Reset();
2125 }
2126 for (auto& h : *mHist2D) {
2127 h.Reset();
2128 }
2129 for (auto& h : *mHist1Dd) {
2130 h.Reset();
2131 }
2132 for (auto& h : *mHistGraph) {
2133 h = TGraphAsymmErrors();
2134 }
2135 mClusterCounts = counts_t();
2136}
2137
2138int32_t GPUQA::DrawQAHistograms(TObjArray* qcout)
2139{
2140 const auto oldRootIgnoreLevel = gErrorIgnoreLevel;
2141 gErrorIgnoreLevel = kWarning;
2142 if (!mQAInitialized) {
2143 throw std::runtime_error("QA not initialized");
2144 }
2145
2146 if (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2) {
2147 printf("Creating QA Histograms\n");
2148 }
2149
2150 std::vector<Color_t> colorNums(COLORCOUNT);
2151 if (!(qcout || mConfig.writeFileExt == "root" || mConfig.writeFileExt == "C")) {
2152 [[maybe_unused]] static int32_t initColorsInitialized = initColors();
2153 }
2154 for (int32_t i = 0; i < COLORCOUNT; i++) {
2155 colorNums[i] = (qcout || mConfig.writeFileExt == "root" || mConfig.writeFileExt == "C") ? defaultColorNums[i] : mColors[i]->GetNumber();
2156 }
2157
2158 bool mcAvail = mcPresent();
2159 char name[2048], fname[1024];
2160
2161 const int32_t nNewInput = mConfig.inputHistogramsOnly ? 0 : 1;
2162 const int32_t ConfigNumInputs = nNewInput + mConfig.compareInputs.size();
2163
2164 std::vector<std::unique_ptr<TFile>> tin;
2165 for (uint32_t i = 0; i < mConfig.compareInputs.size(); i++) {
2166 tin.emplace_back(std::make_unique<TFile>(mConfig.compareInputs[i].c_str()));
2167 }
2168 std::unique_ptr<TFile> tout = nullptr;
2169 if (mConfig.output.size()) {
2170 tout = std::make_unique<TFile>(mConfig.output.c_str(), "RECREATE");
2171 }
2172
2173 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2174 float legendSpacingString = 0.025;
2175 for (int32_t i = 0; i < ConfigNumInputs; i++) {
2176 GetName(fname, i);
2177 if (strlen(fname) * 0.006 > legendSpacingString) {
2178 legendSpacingString = strlen(fname) * 0.006;
2179 }
2180 }
2181
2182 // Create Canvas / Pads for Efficiency Histograms
2183 if (mQATasks & taskTrackingEff) {
2184 for (int32_t ii = 0; ii < 6; ii++) {
2185 snprintf(name, 1024, "eff_vs_%s_layout", VSPARAMETER_NAMES[ii]);
2186 mCEff[ii] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2187 mCEff[ii]->cd();
2188 float dy = 1. / 2.;
2189 mPEff[ii][0] = createGarbageCollected<TPad>("p0", "", 0.0, dy * 0, 0.5, dy * 1);
2190 mPEff[ii][0]->Draw();
2191 mPEff[ii][0]->SetRightMargin(0.04);
2192 mPEff[ii][1] = createGarbageCollected<TPad>("p1", "", 0.5, dy * 0, 1.0, dy * 1);
2193 mPEff[ii][1]->Draw();
2194 mPEff[ii][1]->SetRightMargin(0.04);
2195 mPEff[ii][2] = createGarbageCollected<TPad>("p2", "", 0.0, dy * 1, 0.5, dy * 2 - .001);
2196 mPEff[ii][2]->Draw();
2197 mPEff[ii][2]->SetRightMargin(0.04);
2198 mPEff[ii][3] = createGarbageCollected<TPad>("p3", "", 0.5, dy * 1, 1.0, dy * 2 - .001);
2199 mPEff[ii][3]->Draw();
2200 mPEff[ii][3]->SetRightMargin(0.04);
2201 mLEff[ii] = createGarbageCollected<TLegend>(0.92 - legendSpacingString * 1.45, 0.83 - (0.93 - 0.82) / 2. * (float)ConfigNumInputs, 0.98, 0.849);
2202 SetLegend(mLEff[ii]);
2203 }
2204 }
2205
2206 // Create Canvas / Pads for Resolution Histograms
2207 if (mQATasks & taskTrackingRes) {
2208 for (int32_t ii = 0; ii < 7; ii++) {
2209 if (ii == 6) {
2210 snprintf(name, 1024, "res_integral_layout");
2211 } else {
2212 snprintf(name, 1024, "res_vs_%s_layout", VSPARAMETER_NAMES[ii]);
2213 }
2214 mCRes[ii] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2215 mCRes[ii]->cd();
2216 gStyle->SetOptFit(1);
2217
2218 float dy = 1. / 2.;
2219 mPRes[ii][3] = createGarbageCollected<TPad>("p0", "", 0.0, dy * 0, 0.5, dy * 1);
2220 mPRes[ii][3]->Draw();
2221 mPRes[ii][3]->SetRightMargin(0.04);
2222 mPRes[ii][4] = createGarbageCollected<TPad>("p1", "", 0.5, dy * 0, 1.0, dy * 1);
2223 mPRes[ii][4]->Draw();
2224 mPRes[ii][4]->SetRightMargin(0.04);
2225 mPRes[ii][0] = createGarbageCollected<TPad>("p2", "", 0.0, dy * 1, 1. / 3., dy * 2 - .001);
2226 mPRes[ii][0]->Draw();
2227 mPRes[ii][0]->SetRightMargin(0.04);
2228 mPRes[ii][0]->SetLeftMargin(0.15);
2229 mPRes[ii][1] = createGarbageCollected<TPad>("p3", "", 1. / 3., dy * 1, 2. / 3., dy * 2 - .001);
2230 mPRes[ii][1]->Draw();
2231 mPRes[ii][1]->SetRightMargin(0.04);
2232 mPRes[ii][1]->SetLeftMargin(0.135);
2233 mPRes[ii][2] = createGarbageCollected<TPad>("p4", "", 2. / 3., dy * 1, 1.0, dy * 2 - .001);
2234 mPRes[ii][2]->Draw();
2235 mPRes[ii][2]->SetRightMargin(0.06);
2236 mPRes[ii][2]->SetLeftMargin(0.135);
2237 if (ii < 6) {
2238 mLRes[ii] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.45, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949);
2239 SetLegend(mLRes[ii]);
2240 }
2241 }
2242 }
2243
2244 // Create Canvas / Pads for Pull Histograms
2245 if (mQATasks & taskTrackingResPull) {
2246 for (int32_t ii = 0; ii < 7; ii++) {
2247 if (ii == 6) {
2248 snprintf(name, 1024, "pull_integral_layout");
2249 } else {
2250 snprintf(name, 1024, "pull_vs_%s_layout", VSPARAMETER_NAMES[ii]);
2251 }
2252 mCPull[ii] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2253 mCPull[ii]->cd();
2254 gStyle->SetOptFit(1);
2255
2256 float dy = 1. / 2.;
2257 mPPull[ii][3] = createGarbageCollected<TPad>("p0", "", 0.0, dy * 0, 0.5, dy * 1);
2258 mPPull[ii][3]->Draw();
2259 mPPull[ii][3]->SetRightMargin(0.04);
2260 mPPull[ii][4] = createGarbageCollected<TPad>("p1", "", 0.5, dy * 0, 1.0, dy * 1);
2261 mPPull[ii][4]->Draw();
2262 mPPull[ii][4]->SetRightMargin(0.04);
2263 mPPull[ii][0] = createGarbageCollected<TPad>("p2", "", 0.0, dy * 1, 1. / 3., dy * 2 - .001);
2264 mPPull[ii][0]->Draw();
2265 mPPull[ii][0]->SetRightMargin(0.04);
2266 mPPull[ii][0]->SetLeftMargin(0.15);
2267 mPPull[ii][1] = createGarbageCollected<TPad>("p3", "", 1. / 3., dy * 1, 2. / 3., dy * 2 - .001);
2268 mPPull[ii][1]->Draw();
2269 mPPull[ii][1]->SetRightMargin(0.04);
2270 mPPull[ii][1]->SetLeftMargin(0.135);
2271 mPPull[ii][2] = createGarbageCollected<TPad>("p4", "", 2. / 3., dy * 1, 1.0, dy * 2 - .001);
2272 mPPull[ii][2]->Draw();
2273 mPPull[ii][2]->SetRightMargin(0.06);
2274 mPPull[ii][2]->SetLeftMargin(0.135);
2275 if (ii < 6) {
2276 mLPull[ii] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.45, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949);
2277 SetLegend(mLPull[ii]);
2278 }
2279 }
2280 }
2281
2282 // Create Canvas for Cluster Histos
2283 if (mQATasks & taskClusterAttach) {
2284 for (int32_t i = 0; i < 3; i++) {
2285 snprintf(name, 1024, "clusters_%s_layout", CLUSTER_TYPES[i]);
2286 mCClust[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2287 mCClust[i]->cd();
2288 mPClust[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2289 mPClust[i]->Draw();
2290 float y1 = i != 1 ? 0.77 : 0.27, y2 = i != 1 ? 0.9 : 0.42;
2291 mLClust[i] = createGarbageCollected<TLegend>(i == 2 ? 0.1 : (0.65 - legendSpacingString * 1.45), y2 - (y2 - y1) * (ConfigNumInputs + (i != 1) / 2.) + 0.005, i == 2 ? (0.3 + legendSpacingString * 1.45) : 0.9, y2);
2292 SetLegend(mLClust[i]);
2293 }
2294 }
2295
2296 // Create Canvas for track statistic histos
2297 if (mQATasks & taskTrackStatistics) {
2298 for (int32_t i = 0; i < 3; i++) {
2299 snprintf(name, 2048, "ctracks%s", i ? (i == 2 ? "dedxmax" : "dedxtot") : "pt");
2300 mCTracks[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2301 mCTracks[i]->cd();
2302 mPTracks[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2303 mPTracks[i]->Draw();
2304 mLTracks[i] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.5, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949);
2305 SetLegend(mLTracks[i], true);
2306 }
2307
2308 for (int32_t i = 0; i < 2; i++) {
2309 snprintf(name, 2048, "ctrackst0%d", i);
2310 mCT0[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2311 mCT0[i]->cd();
2312 mPT0[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2313 mPT0[i]->Draw();
2314 mLT0[i] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.45, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949);
2315 SetLegend(mLT0[i]);
2316
2317 snprintf(name, 2048, "cncl%d", i);
2318 mCNCl[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2319 mCNCl[i]->cd();
2320 mPNCl[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2321 mPNCl[i]->Draw();
2322 mLNCl[i] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.45, 0.93 - (0.93 - 0.86) / 2. * (float)ConfigNumInputs, 0.98, 0.949); // TODO: Fix sizing of legend, and also fix font size
2323 SetLegend(mLNCl[i], true);
2324 }
2325
2326 mCClXY = createGarbageCollected<TCanvas>("clxy", "clxy", 0, 0, 700, 700. * 2. / 3.);
2327 mCClXY->cd();
2328 mPClXY = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2329 mPClXY->Draw();
2330 }
2331
2332 if (mQATasks & taskClusterRejection) {
2333 for (int32_t i = 0; i < 3; i++) {
2334 snprintf(name, 2048, "cnclrej%d", i);
2335 mCClRej[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2336 mCClRej[i]->cd();
2337 mPClRej[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2338 mPClRej[i]->Draw();
2339 }
2340 mCClRejP = createGarbageCollected<TCanvas>("cnclrejp", "cnclrejp", 0, 0, 700, 700. * 2. / 3.);
2341 mCClRejP->cd();
2342 mPClRejP = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2343 mPClRejP->Draw();
2344 }
2345
2346 if (mQATasks & taskClusterAttach) {
2347 for (int32_t i = 0; i < 4; i++) {
2348 snprintf(name, 2048, "cpadrow%d", i);
2349 mCPadRow[i] = createGarbageCollected<TCanvas>(name, name, 0, 0, 700, 700. * 2. / 3.);
2350 mCPadRow[i]->cd();
2351 mPPadRow[i] = createGarbageCollected<TPad>("p0", "", 0.0, 0.0, 1.0, 1.0);
2352 mPPadRow[i]->Draw();
2353 }
2354 }
2355 }
2356
2357 if (mConfig.enableLocalOutput && !mConfig.inputHistogramsOnly && (mQATasks & taskTrackingEff) && mcPresent()) {
2358 GPUInfo("QA Stats: Eff: Tracks Prim %d (Eta %d, Pt %d) %f%% (%f%%) Sec %d (Eta %d, Pt %d) %f%% (%f%%) - Res: Tracks %d (Eta %d, Pt %d)", (int32_t)mEff[3][1][0][0]->GetEntries(), (int32_t)mEff[3][1][0][3]->GetEntries(), (int32_t)mEff[3][1][0][4]->GetEntries(),
2359 mEff[0][0][0][0]->GetSumOfWeights() / std::max(1., mEff[3][0][0][0]->GetSumOfWeights()), mEff[0][1][0][0]->GetSumOfWeights() / std::max(1., mEff[3][1][0][0]->GetSumOfWeights()), (int32_t)mEff[3][1][1][0]->GetEntries(), (int32_t)mEff[3][1][1][3]->GetEntries(),
2360 (int32_t)mEff[3][1][1][4]->GetEntries(), mEff[0][0][1][0]->GetSumOfWeights() / std::max(1., mEff[3][0][1][0]->GetSumOfWeights()), mEff[0][1][1][0]->GetSumOfWeights() / std::max(1., mEff[3][1][1][0]->GetSumOfWeights()), (int32_t)mRes2[0][0]->GetEntries(),
2361 (int32_t)mRes2[0][3]->GetEntries(), (int32_t)mRes2[0][4]->GetEntries());
2362 }
2363
2364 int32_t flagShowVsPtLog = (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) ? 1 : 0;
2365
2366 if (mQATasks & taskTrackingEff) {
2367 // Process / Draw Efficiency Histograms
2368 for (int32_t ii = 0; ii < 5 + flagShowVsPtLog; ii++) {
2369 int32_t i = ii == 5 ? 4 : ii;
2370 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2371 for (int32_t j = 0; j < 4; j++) {
2372 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2373 mPEff[ii][j]->cd();
2374 if (ii == 5) {
2375 mPEff[ii][j]->SetLogx();
2376 }
2377 }
2378 for (int32_t l = 0; l < 3; l++) {
2379 if (k == 0 && mConfig.inputHistogramsOnly == 0 && ii != 5) {
2380 if (l == 0) {
2381 // Divide eff, compute all for fake/clone
2382 auto oldLevel = gErrorIgnoreLevel;
2383 gErrorIgnoreLevel = kError;
2384 mEffResult[0][j / 2][j % 2][i]->Divide(mEff[l][j / 2][j % 2][i], mEff[5][j / 2][j % 2][i], "cl=0.683 b(1,1) mode");
2385 gErrorIgnoreLevel = oldLevel;
2386 mEff[3][j / 2][j % 2][i]->Reset(); // Sum up rec + clone + fake for fake rate
2387 mEff[3][j / 2][j % 2][i]->Add(mEff[0][j / 2][j % 2][i]);
2388 mEff[3][j / 2][j % 2][i]->Add(mEff[1][j / 2][j % 2][i]);
2389 mEff[3][j / 2][j % 2][i]->Add(mEff[2][j / 2][j % 2][i]);
2390 mEff[4][j / 2][j % 2][i]->Reset(); // Sum up rec + clone for clone rate
2391 mEff[4][j / 2][j % 2][i]->Add(mEff[0][j / 2][j % 2][i]);
2392 mEff[4][j / 2][j % 2][i]->Add(mEff[1][j / 2][j % 2][i]);
2393 } else {
2394 // Divide fake/clone
2395 auto oldLevel = gErrorIgnoreLevel;
2396 gErrorIgnoreLevel = kError;
2397 mEffResult[l][j / 2][j % 2][i]->Divide(mEff[l][j / 2][j % 2][i], mEff[l == 1 ? 4 : 3][j / 2][j % 2][i], "cl=0.683 b(1,1) mode");
2398 gErrorIgnoreLevel = oldLevel;
2399 }
2400 }
2401
2402 TGraphAsymmErrors* e = mEffResult[l][j / 2][j % 2][i];
2403
2404 if (!mConfig.inputHistogramsOnly && k == 0) {
2405 if (tout) {
2406 mEff[l][j / 2][j % 2][i]->Write();
2407 e->Write();
2408 if (l == 2) {
2409 mEff[3][j / 2][j % 2][i]->Write(); // Store also all histogram!
2410 mEff[4][j / 2][j % 2][i]->Write(); // Store also all histogram!
2411 }
2412 }
2413 } else if (GetHist(e, tin, k, nNewInput) == nullptr) {
2414 continue;
2415 }
2416 e->SetTitle(EFFICIENCY_TITLES[j]);
2417 e->GetYaxis()->SetTitle("(Efficiency)");
2418 e->GetXaxis()->SetTitle(XAXIS_TITLES[i]);
2419
2420 e->SetLineWidth(1);
2421 e->SetLineStyle(CONFIG_DASHED_MARKERS ? k + 1 : 1);
2422 SetAxisSize(e);
2423 if (qcout && !mConfig.shipToQCAsCanvas) {
2424 qcout->Add(e);
2425 }
2426 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2427 continue;
2428 }
2429 e->SetMarkerColor(kBlack);
2430 e->SetLineColor(colorNums[(k < 3 ? (l * 3 + k) : (k * 3 + l)) % COLORCOUNT]);
2431 e->GetHistogram()->GetYaxis()->SetRangeUser(-0.02, 1.02);
2432 e->Draw(k || l ? "same P" : "AP");
2433 if (j == 0) {
2434 GetName(fname, k);
2435 mLEff[ii]->AddEntry(e, Form("%s%s", fname, EFF_NAMES[l]), "l");
2436 }
2437 }
2438 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2439 continue;
2440 }
2441 mCEff[ii]->cd();
2442 ChangePadTitleSize(mPEff[ii][j], 0.056);
2443 }
2444 }
2445 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2446 continue;
2447 }
2448
2449 mLEff[ii]->Draw();
2450
2451 if (qcout) {
2452 qcout->Add(mCEff[ii]);
2453 }
2454 if (!mConfig.enableLocalOutput) {
2455 continue;
2456 }
2457 doPerfFigure(0.2, 0.295, 0.025);
2458 mCEff[ii]->Print(Form("%s/eff_vs_%s.pdf", mConfig.plotsDir.c_str(), VSPARAMETER_NAMES[ii]));
2459 if (mConfig.writeFileExt != "") {
2460 mCEff[ii]->Print(Form("%s/eff_vs_%s.%s", mConfig.plotsDir.c_str(), VSPARAMETER_NAMES[ii], mConfig.writeFileExt.c_str()));
2461 }
2462 }
2463 }
2464
2465 if (mQATasks & (taskTrackingRes | taskTrackingResPull)) {
2466 // Process / Draw Resolution Histograms
2467 TH1D *resIntegral[5] = {}, *pullIntegral[5] = {};
2468 TCanvas* cfit = nullptr;
2469 std::unique_ptr<TF1> customGaus = std::make_unique<TF1>("G", "[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))");
2470 for (int32_t p = 0; p < 2; p++) {
2471 if ((p == 0 && (mQATasks & taskTrackingRes) == 0) || (p == 1 && (mQATasks & taskTrackingResPull) == 0)) {
2472 continue;
2473 }
2474 for (int32_t ii = 0; ii < 5 + flagShowVsPtLog; ii++) {
2475 TCanvas* can = p ? mCPull[ii] : mCRes[ii];
2476 TLegend* leg = p ? mLPull[ii] : mLRes[ii];
2477 int32_t i = ii == 5 ? 4 : ii;
2478 for (int32_t j = 0; j < 5; j++) {
2479 TH2F* src = p ? mPull2[j][i] : mRes2[j][i];
2480 TH1F** dst = p ? mPull[j][i] : mRes[j][i];
2481 TH1D*& dstIntegral = p ? pullIntegral[j] : resIntegral[j];
2482 TPad* pad = p ? mPPull[ii][j] : mPRes[ii][j];
2483
2484 if (!mConfig.inputHistogramsOnly && ii != 5) {
2485 if (cfit == nullptr) {
2486 cfit = createGarbageCollected<TCanvas>();
2487 }
2488 cfit->cd();
2489
2490 TAxis* axis = src->GetYaxis();
2491 int32_t nBins = axis->GetNbins();
2492 int32_t integ = 1;
2493 for (int32_t bin = 1; bin <= nBins; bin++) {
2494 int32_t bin0 = std::max(bin - integ, 0);
2495 int32_t bin1 = std::min(bin + integ, nBins);
2496 std::unique_ptr<TH1D> proj{src->ProjectionX("proj", bin0, bin1)};
2497 proj->ClearUnderflowAndOverflow();
2498 if (proj->GetEntries()) {
2499 uint32_t rebin = 1;
2500 while (proj->GetMaximum() < 50 && rebin < sizeof(RES_AXIS_BINS) / sizeof(RES_AXIS_BINS[0])) {
2501 proj->Rebin(RES_AXIS_BINS[rebin - 1] / RES_AXIS_BINS[rebin]);
2502 rebin++;
2503 }
2504
2505 if (proj->GetEntries() < 20 || proj->GetRMS() < 0.00001) {
2506 dst[0]->SetBinContent(bin, proj->GetRMS());
2507 dst[0]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2508 dst[1]->SetBinContent(bin, proj->GetMean());
2509 dst[1]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2510 } else {
2511 proj->GetXaxis()->SetRange(0, 0);
2512 proj->GetXaxis()->SetRangeUser(std::max(proj->GetXaxis()->GetXmin(), proj->GetMean() - 3. * proj->GetRMS()), std::min(proj->GetXaxis()->GetXmax(), proj->GetMean() + 3. * proj->GetRMS()));
2513 bool forceLogLike = proj->GetMaximum() < 20;
2514 for (int32_t k = forceLogLike ? 2 : 0; k < 3; k++) {
2515 proj->Fit("gaus", forceLogLike || k == 2 ? "sQl" : (k ? "sQww" : "sQ"));
2516 TF1* fitFunc = proj->GetFunction("gaus");
2517
2518 if (k && !forceLogLike) {
2519 customGaus->SetParameters(fitFunc->GetParameter(0), fitFunc->GetParameter(1), fitFunc->GetParameter(2));
2520 proj->Fit(customGaus.get(), "sQ");
2521 fitFunc = customGaus.get();
2522 }
2523
2524 const float sigma = fabs(fitFunc->GetParameter(2));
2525 dst[0]->SetBinContent(bin, sigma);
2526 dst[1]->SetBinContent(bin, fitFunc->GetParameter(1));
2527 dst[0]->SetBinError(bin, fitFunc->GetParError(2));
2528 dst[1]->SetBinError(bin, fitFunc->GetParError(1));
2529
2530 const bool fail1 = sigma <= 0.f;
2531 const bool fail2 = fabs(proj->GetMean() - dst[1]->GetBinContent(bin)) > std::min<float>(p ? PULL_AXIS : (mConfig.nativeFitResolutions ? RES_AXES_NATIVE[j] : RES_AXES[j]), 3.f * proj->GetRMS());
2532 const bool fail3 = dst[0]->GetBinContent(bin) > 3.f * proj->GetRMS() || dst[0]->GetBinError(bin) > 1 || dst[1]->GetBinError(bin) > 1;
2533 const bool fail4 = fitFunc->GetParameter(0) < proj->GetMaximum() / 5.;
2534 const bool fail = fail1 || fail2 || fail3 || fail4;
2535 // if (p == 0 && ii == 4 && j == 2) DrawHisto(proj, Form("Hist_bin_%d-%d_vs_%d____%d_%d___%f-%f___%f-%f___%d.pdf", p, j, ii, bin, k, dst[0]->GetBinContent(bin), proj->GetRMS(), dst[1]->GetBinContent(bin), proj->GetMean(), (int32_t) fail), "");
2536
2537 if (!fail) {
2538 break;
2539 } else if (k >= 2) {
2540 dst[0]->SetBinContent(bin, proj->GetRMS());
2541 dst[0]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2542 dst[1]->SetBinContent(bin, proj->GetMean());
2543 dst[1]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2544 }
2545 }
2546 }
2547 } else {
2548 dst[0]->SetBinContent(bin, 0.f);
2549 dst[0]->SetBinError(bin, 0.f);
2550 dst[1]->SetBinContent(bin, 0.f);
2551 dst[1]->SetBinError(bin, 0.f);
2552 }
2553 }
2554 if (ii == 0) {
2555 dstIntegral = src->ProjectionX(mConfig.nativeFitResolutions ? PARAMETER_NAMES_NATIVE[j] : PARAMETER_NAMES[j], 0, nBins + 1);
2556 uint32_t rebin = 1;
2557 while (dstIntegral->GetMaximum() < 50 && rebin < sizeof(RES_AXIS_BINS) / sizeof(RES_AXIS_BINS[0])) {
2558 dstIntegral->Rebin(RES_AXIS_BINS[rebin - 1] / RES_AXIS_BINS[rebin]);
2559 rebin++;
2560 }
2561 }
2562 }
2563 if (ii == 0) {
2564 if (mConfig.inputHistogramsOnly) {
2565 dstIntegral = createGarbageCollected<TH1D>();
2566 }
2567 dstIntegral->SetName(Form(p ? "IntPull%s" : "IntRes%s", VSPARAMETER_NAMES[j]));
2568 dstIntegral->SetTitle(Form(p ? "%s Pull" : "%s Resolution", p || mConfig.nativeFitResolutions ? PARAMETER_NAMES_NATIVE[j] : PARAMETER_NAMES[j]));
2569 }
2570 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2571 pad->cd();
2572 }
2573 int32_t numColor = 0;
2574 float tmpMax = -1000.;
2575 float tmpMin = 1000.;
2576
2577 for (int32_t l = 0; l < 2; l++) {
2578 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2579 TH1F* e = dst[l];
2580 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2581 continue;
2582 }
2583 if (nNewInput && k == 0 && ii != 5) {
2584 if (p == 0) {
2585 e->Scale(mConfig.nativeFitResolutions ? SCALE_NATIVE[j] : SCALE[j]);
2586 }
2587 }
2588 if (ii == 4) {
2589 e->GetXaxis()->SetRangeUser(0.2, PT_MAX);
2590 } else if (LOG_PT_MIN > 0 && ii == 5) {
2591 e->GetXaxis()->SetRangeUser(LOG_PT_MIN, PT_MAX);
2592 } else if (ii == 5) {
2593 e->GetXaxis()->SetRange(1, 0);
2594 }
2595 e->SetMinimum(-1111);
2596 e->SetMaximum(-1111);
2597
2598 if (e->GetMaximum() > tmpMax) {
2599 tmpMax = e->GetMaximum();
2600 }
2601 if (e->GetMinimum() < tmpMin) {
2602 tmpMin = e->GetMinimum();
2603 }
2604 }
2605 }
2606
2607 float tmpSpan;
2608 tmpSpan = tmpMax - tmpMin;
2609 tmpMax += tmpSpan * .02;
2610 tmpMin -= tmpSpan * .02;
2611 if (j == 2 && i < 3) {
2612 tmpMax += tmpSpan * 0.13 * ConfigNumInputs;
2613 }
2614
2615 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2616 for (int32_t l = 0; l < 2; l++) {
2617 TH1F* e = dst[l];
2618 if (!mConfig.inputHistogramsOnly && k == 0) {
2619 e->SetTitle(Form(p ? "%s Pull" : "%s Resolution", p || mConfig.nativeFitResolutions ? PARAMETER_NAMES_NATIVE[j] : PARAMETER_NAMES[j]));
2620 e->SetStats(kFALSE);
2621 if (tout) {
2622 if (l == 0) {
2623 mRes2[j][i]->SetOption("colz");
2624 mRes2[j][i]->Write();
2625 }
2626 e->Write();
2627 }
2628 } else if (GetHist(e, tin, k, nNewInput) == nullptr) {
2629 continue;
2630 }
2631 e->SetMaximum(tmpMax);
2632 e->SetMinimum(tmpMin);
2633 e->SetLineWidth(1);
2634 e->SetLineStyle(CONFIG_DASHED_MARKERS ? k + 1 : 1);
2635 SetAxisSize(e);
2636 e->GetYaxis()->SetTitle(p ? AXIS_TITLES_PULL[j] : (mConfig.nativeFitResolutions ? AXIS_TITLES_NATIVE[j] : AXIS_TITLES[j]));
2637 e->GetXaxis()->SetTitle(XAXIS_TITLES[i]);
2638 if (LOG_PT_MIN > 0 && ii == 5) {
2639 e->GetXaxis()->SetRangeUser(LOG_PT_MIN, PT_MAX);
2640 }
2641
2642 if (j == 0) {
2643 e->GetYaxis()->SetTitleOffset(1.5);
2644 } else if (j < 3) {
2645 e->GetYaxis()->SetTitleOffset(1.4);
2646 }
2647 if (qcout && !mConfig.shipToQCAsCanvas) {
2648 qcout->Add(e);
2649 }
2650 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2651 continue;
2652 }
2653
2654 e->SetMarkerColor(kBlack);
2655 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2656 e->Draw(k || l ? "same" : "");
2657 if (j == 0) {
2658 GetName(fname, k);
2659 leg->AddEntry(e, Form("%s%s", fname, l ? "Mean" : (p ? "Pull" : "Resolution")), "l");
2660 }
2661 }
2662 }
2663 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2664 continue;
2665 }
2666
2667 if (ii == 5) {
2668 pad->SetLogx();
2669 }
2670 can->cd();
2671 if (j == 4) {
2672 ChangePadTitleSize(pad, 0.056);
2673 }
2674 }
2675 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2676 continue;
2677 }
2678
2679 leg->Draw();
2680
2681 if (qcout) {
2682 qcout->Add(can);
2683 }
2684 if (!mConfig.enableLocalOutput) {
2685 continue;
2686 }
2687 doPerfFigure(0.2, 0.295, 0.025);
2688 can->Print(Form(p ? "%s/pull_vs_%s.pdf" : "%s/res_vs_%s.pdf", mConfig.plotsDir.c_str(), VSPARAMETER_NAMES[ii]));
2689 if (mConfig.writeFileExt != "") {
2690 can->Print(Form(p ? "%s/pull_vs_%s.%s" : "%s/res_vs_%s.%s", mConfig.plotsDir.c_str(), VSPARAMETER_NAMES[ii], mConfig.writeFileExt.c_str()));
2691 }
2692 }
2693 }
2694
2695 // Process Integral Resolution Histogreams
2696 for (int32_t p = 0; p < 2; p++) {
2697 if ((p == 0 && (mQATasks & taskTrackingRes) == 0) || (p == 1 && (mQATasks & taskTrackingResPull) == 0)) {
2698 continue;
2699 }
2700 TCanvas* can = p ? mCPull[6] : mCRes[6];
2701 for (int32_t i = 0; i < 5; i++) {
2702 TPad* pad = p ? mPPull[6][i] : mPRes[6][i];
2703 TH1D* hist = p ? pullIntegral[i] : resIntegral[i];
2704 int32_t numColor = 0;
2705 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2706 pad->cd();
2707 }
2708 if (!mConfig.inputHistogramsOnly && mcAvail) {
2709 TH1D* e = hist;
2710 if (e && e->GetEntries()) {
2711 e->Fit("gaus", "sQ");
2712 }
2713 }
2714
2715 float tmpMax = 0;
2716 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2717 TH1D* e = hist;
2718 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2719 continue;
2720 }
2721 e->SetMaximum(-1111);
2722 if (e->GetMaximum() > tmpMax) {
2723 tmpMax = e->GetMaximum();
2724 }
2725 }
2726
2727 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2728 TH1D* e = hist;
2729 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2730 continue;
2731 }
2732 e->SetMaximum(tmpMax * 1.02);
2733 e->SetMinimum(tmpMax * -0.02);
2734 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2735 e->Write();
2736 }
2737 if (qcout && !mConfig.shipToQCAsCanvas) {
2738 qcout->Add(e);
2739 }
2740 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2741 continue;
2742 }
2743
2744 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2745 e->Draw(k == 0 ? "" : "same");
2746 }
2747 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2748 continue;
2749 }
2750 can->cd();
2751 }
2752 if (qcout) {
2753 qcout->Add(can);
2754 }
2755 if (!mConfig.enableLocalOutput) {
2756 continue;
2757 }
2758
2759 can->Print(Form(p ? "%s/pull_integral.pdf" : "%s/res_integral.pdf", mConfig.plotsDir.c_str()));
2760 if (mConfig.writeFileExt != "") {
2761 can->Print(Form(p ? "%s/pull_integral.%s" : "%s/res_integral.%s", mConfig.plotsDir.c_str(), mConfig.writeFileExt.c_str()));
2762 }
2763 }
2764 }
2765
2766 uint64_t attachClusterCounts[N_CLS_HIST];
2767 if (mQATasks & taskClusterAttach) {
2768 // Process Cluster Attachment Histograms
2769 if (mConfig.inputHistogramsOnly == 0) {
2770 for (int32_t i = N_CLS_HIST; i < N_CLS_TYPE * N_CLS_HIST - 1; i++) {
2771 mClusters[i]->Sumw2(true);
2772 }
2773 double totalVal = 0;
2774 if (!CLUST_HIST_INT_SUM) {
2775 for (int32_t j = 0; j < mClusters[N_CLS_HIST - 1]->GetXaxis()->GetNbins() + 2; j++) {
2776 totalVal += mClusters[N_CLS_HIST - 1]->GetBinContent(j);
2777 }
2778 }
2779 if (totalVal == 0.) {
2780 totalVal = 1.;
2781 }
2782 for (int32_t i = 0; i < N_CLS_HIST; i++) {
2783 double val = 0;
2784 for (int32_t j = 0; j < mClusters[i]->GetXaxis()->GetNbins() + 2; j++) {
2785 val += mClusters[i]->GetBinContent(j);
2786 mClusters[2 * N_CLS_HIST - 1 + i]->SetBinContent(j, val / totalVal);
2787 }
2788 attachClusterCounts[i] = val;
2789 }
2790
2791 if (!CLUST_HIST_INT_SUM) {
2792 for (int32_t i = 0; i < N_CLS_HIST; i++) {
2793 mClusters[2 * N_CLS_HIST - 1 + i]->SetMaximum(1.02);
2794 mClusters[2 * N_CLS_HIST - 1 + i]->SetMinimum(-0.02);
2795 }
2796 }
2797
2798 for (int32_t i = 0; i < N_CLS_HIST - 1; i++) {
2799 auto oldLevel = gErrorIgnoreLevel;
2800 gErrorIgnoreLevel = kError;
2801 mClusters[N_CLS_HIST + i]->Divide(mClusters[i], mClusters[N_CLS_HIST - 1], 1, 1, "B");
2802 gErrorIgnoreLevel = oldLevel;
2803 mClusters[N_CLS_HIST + i]->SetMinimum(-0.02);
2804 mClusters[N_CLS_HIST + i]->SetMaximum(1.02);
2805 }
2806 }
2807
2808 float tmpMax[2] = {0, 0}, tmpMin[2] = {0, 0};
2809 for (int32_t l = 0; l <= CLUST_HIST_INT_SUM; l++) {
2810 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2811 TH1* e = mClusters[l ? (N_CLS_TYPE * N_CLS_HIST - 2) : (N_CLS_HIST - 1)];
2812 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2813 continue;
2814 }
2815 e->SetMinimum(-1111);
2816 e->SetMaximum(-1111);
2817 if (l == 0) {
2818 e->GetXaxis()->SetRange(2, AXIS_BINS[4]);
2819 }
2820 if (e->GetMaximum() > tmpMax[l]) {
2821 tmpMax[l] = e->GetMaximum();
2822 }
2823 if (e->GetMinimum() < tmpMin[l]) {
2824 tmpMin[l] = e->GetMinimum();
2825 }
2826 }
2827 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2828 for (int32_t i = 0; i < N_CLS_HIST; i++) {
2829 TH1* e = mClusters[l ? (2 * N_CLS_HIST - 1 + i) : i];
2830 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2831 continue;
2832 }
2833 e->SetMaximum(tmpMax[l] * 1.02);
2834 e->SetMinimum(tmpMax[l] * -0.02);
2835 }
2836 }
2837 }
2838
2839 for (int32_t i = 0; i < N_CLS_TYPE; i++) {
2840 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2841 mPClust[i]->cd();
2842 mPClust[i]->SetLogx();
2843 }
2844 int32_t begin = i == 2 ? (2 * N_CLS_HIST - 1) : (i == 1 ? N_CLS_HIST : 0);
2845 int32_t end = i == 2 ? (3 * N_CLS_HIST - 1) : (i == 1 ? (2 * N_CLS_HIST - 1) : N_CLS_HIST);
2846 int32_t numColor = 0;
2847 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2848 for (int32_t j = end - 1; j >= begin; j--) {
2849 TH1* e = mClusters[j];
2850 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2851 continue;
2852 }
2853
2854 e->SetTitle(mConfig.plotsNoTitle ? "" : CLUSTER_TITLES[i]);
2855 e->GetYaxis()->SetTitle(i == 0 ? "Number of TPC clusters" : (i == 1 ? "Fraction of TPC clusters" : (CLUST_HIST_INT_SUM ? "Total TPC clusters (integrated)" : "Fraction of TPC clusters (integrated)")));
2856 e->GetXaxis()->SetTitle("#it{p}_{Tmc} (GeV/#it{c})");
2857 e->GetXaxis()->SetTitleOffset(1.1);
2858 e->GetXaxis()->SetLabelOffset(-0.005);
2859 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2860 e->Write();
2861 }
2862 e->SetStats(kFALSE);
2863 e->SetLineWidth(1);
2864 e->SetLineStyle(CONFIG_DASHED_MARKERS ? j + 1 : 1);
2865 if (i == 0) {
2866 e->GetXaxis()->SetRange(2, AXIS_BINS[4]);
2867 }
2868 if (qcout && !mConfig.shipToQCAsCanvas) {
2869 qcout->Add(e);
2870 }
2871 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2872 continue;
2873 }
2874
2875 e->SetMarkerColor(kBlack);
2876 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2877 e->Draw(j == end - 1 && k == 0 ? "" : "same");
2878 GetName(fname, k);
2879 mLClust[i]->AddEntry(e, Form("%s%s", fname, CLUSTER_NAMES[j - begin]), "l");
2880 }
2881 }
2882 if (ConfigNumInputs == 1) {
2883 TH1* e = reinterpret_cast<TH1F*>(mClusters[begin + CL_att_adj]->Clone());
2884 e->Add(mClusters[begin + CL_prot], -1);
2885 if (qcout && !mConfig.shipToQCAsCanvas) {
2886 qcout->Add(e);
2887 }
2888 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2889 continue;
2890 }
2891
2892 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2893 e->Draw("same");
2894 mLClust[i]->AddEntry(e, "Removed (Strategy A)", "l");
2895 }
2896 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2897 continue;
2898 }
2899
2900 mLClust[i]->Draw();
2901
2902 if (qcout) {
2903 qcout->Add(mCClust[i]);
2904 }
2905 if (!mConfig.enableLocalOutput) {
2906 continue;
2907 }
2908 doPerfFigure(i == 0 ? 0.37 : (i == 1 ? 0.34 : 0.6), 0.295, 0.030);
2909 mCClust[i]->cd();
2910 mCClust[i]->Print(Form(i == 2 ? "%s/clusters_integral.pdf" : (i == 1 ? "%s/clusters_relative.pdf" : "%s/clusters.pdf"), mConfig.plotsDir.c_str()));
2911 if (mConfig.writeFileExt != "") {
2912 mCClust[i]->Print(Form(i == 2 ? "%s/clusters_integral.%s" : (i == 1 ? "%s/clusters_relative.%s" : "%s/clusters.%s"), mConfig.plotsDir.c_str(), mConfig.writeFileExt.c_str()));
2913 }
2914 }
2915
2916 for (int32_t i = 0; i < 4; i++) {
2917 auto* e = mPadRow[i];
2918 if (tout && !mConfig.inputHistogramsOnly) {
2919 e->Write();
2920 }
2921 mPPadRow[i]->cd();
2922 e->SetOption("colz");
2923 std::string title = "First Track Pad Row (p_{T} > 1GeV, N_{Cl} #geq " + std::to_string(PADROW_CHECK_MINCLS);
2924 if (i >= 2) {
2925 title += ", row_{trk} > row_{MC} + 3, row_{MC} < 10";
2926 }
2927 if (i >= 3) {
2928 title += ", #Phi_{Cl} < 0.15";
2929 }
2930 title += ")";
2931
2932 e->SetTitle(mConfig.plotsNoTitle ? "" : title.c_str());
2933 e->GetXaxis()->SetTitle(i == 3 ? "Local Occupancy" : (i ? "#Phi_{Cl} (sector)" : "First MC Pad Row"));
2934 e->GetYaxis()->SetTitle("First Pad Row");
2935 e->Draw();
2936 mCPadRow[i]->cd();
2937 static const constexpr char* PADROW_NAMES[4] = {"MC", "Phi", "Phi1", "Occ"};
2938 mCPadRow[i]->Print(Form("%s/padRow%s.pdf", mConfig.plotsDir.c_str(), PADROW_NAMES[i]));
2939 if (mConfig.writeFileExt != "") {
2940 mCPadRow[i]->Print(Form("%s/padRow%s.%s", mConfig.plotsDir.c_str(), PADROW_NAMES[i], mConfig.writeFileExt.c_str()));
2941 }
2942 }
2943 }
2944
2945 if (mQATasks & taskClusterAttach) {
2946 for (TH1* h : {static_cast<TH1*>(mLowestRowDifferenceVsPad), static_cast<TH1*>(mLowestRowTrackMultiplicity)}) {
2947 if (tout && !mConfig.inputHistogramsOnly) {
2948 h->Write();
2949 }
2950 if (qcout) {
2951 qcout->Add(h);
2952 }
2953 }
2954 }
2955
2956 // Process cluster count statistics
2957 if ((mQATasks & taskClusterCounts) && !mHaveExternalHists && !mConfig.clusterRejectionHistograms && !mConfig.inputHistogramsOnly) {
2958 DoClusterCounts(attachClusterCounts);
2959 }
2960 if ((qcout || tout) && (mQATasks & taskClusterCounts) && mConfig.clusterRejectionHistograms) {
2961 for (uint32_t i = 0; i < mHistClusterCount.size(); i++) {
2962 if (tout) {
2963 mHistClusterCount[i]->Write();
2964 }
2965 if (qcout) {
2966 qcout->Add(mHistClusterCount[i]);
2967 }
2968 }
2969 }
2970
2971 if (mQATasks & taskTrackStatistics) {
2972 // Process track statistic histograms
2973
2974 float tmpMax = 0.;
2975 for (int32_t k = 0; k < ConfigNumInputs; k++) { // TODO: Simplify this drawing, avoid copy&paste
2976 TH1F* e = mTrackPt;
2977 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2978 continue;
2979 }
2980 e->SetMaximum(-1111);
2981 if (e->GetMaximum() > tmpMax) {
2982 tmpMax = e->GetMaximum();
2983 }
2984 }
2985 mPTracks[0]->cd();
2986 mPTracks[0]->SetLogx();
2987 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2988 TH1F* e = mTrackPt;
2989 if (GetHist(e, tin, k, nNewInput) == nullptr) {
2990 continue;
2991 }
2992 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2993 e->Write();
2994 }
2995 e->SetMaximum(tmpMax * 1.02);
2996 e->SetMinimum(tmpMax * -0.02);
2997 e->SetStats(kFALSE);
2998 e->SetLineWidth(1);
2999 e->SetTitle(mConfig.plotsNoTitle ? "" : "Number of Tracks vs #it{p}_{T}");
3000 e->GetYaxis()->SetTitle("Number of Tracks");
3001 e->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
3002 e->GetXaxis()->SetTitleOffset(1.2);
3003 if (qcout) {
3004 qcout->Add(e);
3005 }
3006 e->SetMarkerColor(kBlack);
3007 e->SetLineColor(colorNums[k % COLORCOUNT]);
3008 e->Draw(k == 0 ? "" : "same");
3009 GetName(fname, k, mConfig.inputHistogramsOnly);
3010 mLTracks[0]->AddEntry(e, Form(mConfig.inputHistogramsOnly ? "%s" : "%sTrack #it{p}_{T}", fname), "l");
3011 }
3012 mLTracks[0]->Draw();
3013 doPerfFigure(0.63, 0.7, 0.030);
3014 mCTracks[0]->cd();
3015 mCTracks[0]->Print(Form("%s/tracks%s.pdf", mConfig.plotsDir.c_str(), "pt"));
3016 if (mConfig.writeFileExt != "") {
3017 mCTracks[0]->Print(Form("%s/tracks%s.%s", mConfig.plotsDir.c_str(), "pt", mConfig.writeFileExt.c_str()));
3018 }
3019
3020 for (int32_t i = 0; i < 2; i++) {
3021 mPTracks[1 + i]->cd();
3022 {
3023 TH2F* e = mTrackdEdx[i];
3024 if (tout && !mConfig.inputHistogramsOnly) {
3025 e->Write();
3026 }
3027 // e->SetStats(kFALSE);
3028 e->SetTitle(mConfig.plotsNoTitle ? "" : (i ? "Track dE/dx (Max)" : "Track dE/dx (Tot)"));
3029 e->GetYaxis()->SetTitle(i ? "dE/dx (max)" : "dE/dx (tot)");
3030 e->GetXaxis()->SetTitle("log(#it{p})");
3031 e->GetXaxis()->SetTitleOffset(1.2);
3032 if (qcout) {
3033 qcout->Add(e);
3034 }
3035 e->SetOption("colz");
3036 e->Draw();
3037 }
3038 mCTracks[1 + i]->cd();
3039 mCTracks[1 + i]->Print(Form("%s/tracks%s.pdf", mConfig.plotsDir.c_str(), i ? "dedx_max" : "dedx_tot"));
3040 if (mConfig.writeFileExt != "") {
3041 mCTracks[1 + i]->Print(Form("%s/tracks%s.%s", mConfig.plotsDir.c_str(), i ? "dedx_max" : "dedx_tot", mConfig.writeFileExt.c_str()));
3042 }
3043
3044 tmpMax = 0.;
3045 for (int32_t k = 0; k < ConfigNumInputs; k++) {
3046 TH1F* e = mT0[i];
3047 if (GetHist(e, tin, k, nNewInput) == nullptr) {
3048 continue;
3049 }
3050 e->SetMaximum(-1111);
3051 if (e->GetMaximum() > tmpMax) {
3052 tmpMax = e->GetMaximum();
3053 }
3054 }
3055 mPT0[i]->cd();
3056 for (int32_t k = 0; k < ConfigNumInputs; k++) {
3057 TH1F* e = mT0[i];
3058 if (GetHist(e, tin, k, nNewInput) == nullptr) {
3059 continue;
3060 }
3061 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
3062 e->Write();
3063 }
3064 e->SetMaximum(tmpMax * 1.02);
3065 e->SetMinimum(tmpMax * -0.02);
3066 e->SetStats(kFALSE);
3067 e->SetLineWidth(1);
3068 e->SetTitle(mConfig.plotsNoTitle ? "" : (i ? "Track t_{0} resolution" : "Track t_{0} distribution"));
3069 e->GetYaxis()->SetTitle("a.u.");
3070 e->GetXaxis()->SetTitle(i ? "t_{0} - t_{0, mc}" : "t_{0}");
3071 if (qcout) {
3072 qcout->Add(e);
3073 }
3074 e->SetMarkerColor(kBlack);
3075 e->SetLineColor(colorNums[k % COLORCOUNT]);
3076 e->Draw(k == 0 ? "" : "same");
3077 GetName(fname, k, mConfig.inputHistogramsOnly);
3078 mLT0[i]->AddEntry(e, Form(mConfig.inputHistogramsOnly ? "%s (%s)" : "%sTrack t_{0} %s", fname, i ? "" : "resolution"), "l");
3079 }
3080 mLT0[i]->Draw();
3081 doPerfFigure(0.63, 0.7, 0.030);
3082 mCT0[i]->cd();
3083 mCT0[i]->Print(Form("%s/t0%s.pdf", mConfig.plotsDir.c_str(), i ? "_res" : ""));
3084 if (mConfig.writeFileExt != "") {
3085 mCT0[i]->Print(Form("%s/t0%s.%s", mConfig.plotsDir.c_str(), i ? "_res" : "", mConfig.writeFileExt.c_str()));
3086 }
3087
3088 tmpMax = 0.;
3089 for (int32_t k = 0; k < ConfigNumInputs; k++) {
3090 TH1F* e = mNCl[i];
3091 if (GetHist(e, tin, k, nNewInput) == nullptr) {
3092 continue;
3093 }
3094 e->SetMaximum(-1111);
3095 if (e->GetMaximum() > tmpMax) {
3096 tmpMax = e->GetMaximum();
3097 }
3098 }
3099 mPNCl[i]->cd();
3100 for (int32_t k = 0; k < ConfigNumInputs; k++) {
3101 TH1F* e = mNCl[i];
3102 if (GetHist(e, tin, k, nNewInput) == nullptr) {
3103 continue;
3104 }
3105 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
3106 e->Write();
3107 }
3108 e->SetMaximum(tmpMax * 1.02);
3109 e->SetMinimum(tmpMax * -0.02);
3110 e->SetStats(kFALSE);
3111 e->SetLineWidth(1);
3112 e->SetTitle(mConfig.plotsNoTitle ? "" : (i ? "Number of Rows with attached Cluster" : "Number of Clusters"));
3113 e->GetYaxis()->SetTitle("a.u.");
3114 e->GetXaxis()->SetTitle(i ? "N_{Rows with Clusters}" : "N_{Clusters}");
3115 if (qcout) {
3116 qcout->Add(e);
3117 }
3118 e->SetMarkerColor(kBlack);
3119 e->SetLineColor(colorNums[k % COLORCOUNT]);
3120 e->Draw(k == 0 ? "" : "same");
3121 GetName(fname, k, mConfig.inputHistogramsOnly);
3122 mLNCl[i]->AddEntry(e, Form(mConfig.inputHistogramsOnly ? "%s" : (i ? "%sN_{Clusters}" : "%sN_{Rows with Clusters}"), fname), "l");
3123 }
3124 mLNCl[i]->Draw();
3125 doPerfFigure(0.6, 0.7, 0.030);
3126 mCNCl[i]->cd();
3127 mCNCl[i]->Print(Form("%s/nClusters%s.pdf", mConfig.plotsDir.c_str(), i ? "_corrected" : ""));
3128 if (mConfig.writeFileExt != "") {
3129 mCNCl[i]->Print(Form("%s/nClusters%s.%s", mConfig.plotsDir.c_str(), i ? "_corrected" : "", mConfig.writeFileExt.c_str()));
3130 }
3131 }
3132
3133 mPClXY->cd(); // TODO: This should become a separate task category
3134 mClXY->SetOption("colz");
3135 mClXY->Draw();
3136 mCClXY->cd();
3137 mCClXY->Print(Form("%s/clustersXY.pdf", mConfig.plotsDir.c_str()));
3138 if (mConfig.writeFileExt != "") {
3139 mCClXY->Print(Form("%s/clustersXY.%s", mConfig.plotsDir.c_str(), mConfig.writeFileExt.c_str()));
3140 }
3141 }
3142
3143 if (mQATasks & taskClusterRejection) {
3144 mClRej[2]->Divide(mClRej[1], mClRej[0]);
3145
3146 for (int32_t i = 0; i < 3; i++) {
3147 if (tout && !mConfig.inputHistogramsOnly) {
3148 mClRej[i]->Write();
3149 }
3150 mPClRej[i]->cd();
3151 mClRej[i]->SetTitle(mConfig.plotsNoTitle ? "" : REJECTED_NAMES[i]);
3152 mClRej[i]->SetOption("colz");
3153 mClRej[i]->Draw();
3154 mCClRej[i]->cd();
3155 mCClRej[i]->Print(Form("%s/clustersRej%d%s.pdf", mConfig.plotsDir.c_str(), i, REJECTED_NAMES[i]));
3156 if (mConfig.writeFileExt != "") {
3157 mCClRej[i]->Print(Form("%s/clustersRej%d%s.%s", mConfig.plotsDir.c_str(), i, REJECTED_NAMES[i], mConfig.writeFileExt.c_str()));
3158 }
3159 }
3160
3161 mPClRejP->cd();
3162 for (int32_t k = 0; k < ConfigNumInputs; k++) {
3163 auto* tmp = mClRej[0];
3164 if (GetHist(tmp, tin, k, nNewInput) == nullptr) {
3165 continue;
3166 }
3167 TH1D* proj1 = tmp->ProjectionY(Form("clrejptmp1%d", k)); // TODO: Clean up names
3168 proj1->SetDirectory(nullptr);
3169 tmp = mClRej[1];
3170 if (GetHist(tmp, tin, k, nNewInput) == nullptr) {
3171 continue;
3172 }
3173 TH1D* proj2 = tmp->ProjectionY(Form("clrejptmp2%d", k));
3174 proj2->SetDirectory(nullptr);
3175
3176 auto* e = mClRejP;
3177 if (GetHist(e, tin, k, nNewInput) == nullptr) {
3178 continue;
3179 }
3180 e->Divide(proj2, proj1);
3181 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
3182 e->Write();
3183 }
3184 delete proj1;
3185 delete proj2;
3186 e->SetMinimum(-0.02);
3187 e->SetMaximum(0.22);
3188 e->SetTitle(mConfig.plotsNoTitle ? "" : "Rejected Clusters");
3189 e->GetXaxis()->SetTitle("Pad Row");
3190 e->GetYaxis()->SetTitle("Rejected Clusters (fraction)");
3191 e->Draw(k == 0 ? "" : "same");
3192 }
3193 mPClRejP->Print(Form("%s/clustersRejProjected.pdf", mConfig.plotsDir.c_str()));
3194 if (mConfig.writeFileExt != "") {
3195 mPClRejP->Print(Form("%s/clustersRejProjected.%s", mConfig.plotsDir.c_str(), mConfig.writeFileExt.c_str()));
3196 }
3197 }
3198
3199 if (tout && !mConfig.inputHistogramsOnly && mConfig.writeMCLabels) {
3200 gInterpreter->GenerateDictionary("vector<vector<int32_t>>", "");
3201 tout->WriteObject(&mcEffBuffer, "mcEffBuffer");
3202 tout->WriteObject(&mcLabelBuffer, "mcLabelBuffer");
3203 remove("AutoDict_vector_vector_int__.cxx");
3204 remove("AutoDict_vector_vector_int___cxx_ACLiC_dict_rdict.pcm");
3205 remove("AutoDict_vector_vector_int___cxx.d");
3206 remove("AutoDict_vector_vector_int___cxx.so");
3207 }
3208
3209 if (tout) {
3210 tout->Close();
3211 }
3212 for (uint32_t i = 0; i < mConfig.compareInputs.size(); i++) {
3213 tin[i]->Close();
3214 }
3215 if (!qcout) {
3216 clearGarbagageCollector();
3217 }
3218 GPUInfo("GPU TPC QA histograms have been written to pdf%s%s files", mConfig.writeFileExt == "" ? "" : " and ", mConfig.writeFileExt.c_str());
3219 gErrorIgnoreLevel = oldRootIgnoreLevel;
3220 return (0);
3221}
3222
3223void GPUQA::PrintClusterCount(int32_t mode, int32_t& num, const char* name, uint64_t n, uint64_t normalization)
3224{
3225 if (mode == 2) {
3226 // do nothing, just count num
3227 } else if (mode == 1) {
3228 char name2[128];
3229 snprintf(name2, 128, "clusterCount%d_", num);
3230 char* ptr = name2 + strlen(name2);
3231 for (uint32_t i = 0; i < strlen(name); i++) {
3232 if ((name[i] >= 'a' && name[i] <= 'z') || (name[i] >= 'A' && name[i] <= 'Z') || (name[i] >= '0' && name[i] <= '9')) {
3233 *(ptr++) = name[i];
3234 }
3235 }
3236 *ptr = 0;
3237 createHist(mHistClusterCount[num], name2, name, 1000, 0, mConfig.histMaxNClusters, 1000, 0, 100);
3238 } else if (mode == 0) {
3239 if (normalization && mConfig.enableLocalOutput) {
3240 for (uint32_t i = 0; i < 1 + (mTextDump != nullptr); i++) {
3241 fprintf(i ? mTextDump : stdout, "\t%40s: %'12" PRIu64 " (%6.2f%%)\n", name, n, 100.f * n / normalization);
3242 }
3243 }
3244 if (mConfig.clusterRejectionHistograms) {
3245 float ratio = 100.f * n / std::max<uint64_t>(normalization, 1);
3246 mHistClusterCount[num]->Fill(normalization, ratio, 1);
3247 }
3248 }
3249 num++;
3250}
3251
3252int32_t GPUQA::DoClusterCounts(uint64_t* attachClusterCounts, int32_t mode)
3253{
3254 if (mConfig.enableLocalOutput && !mConfig.inputHistogramsOnly && mConfig.plotsDir != "") {
3255 mTextDump = fopen((mConfig.plotsDir + "/clusterCounts.txt").c_str(), "w+");
3256 }
3257 int32_t num = 0;
3258 if (mcPresent() && (mQATasks & taskClusterAttach) && attachClusterCounts) {
3259 for (int32_t i = 0; i < N_CLS_HIST; i++) { // TODO: Check that these counts are still printed correctly!
3260 PrintClusterCount(mode, num, CLUSTER_NAMES[i], attachClusterCounts[i], mClusterCounts.nTotal);
3261 }
3262 PrintClusterCount(mode, num, "Unattached", attachClusterCounts[N_CLS_HIST - 1] - attachClusterCounts[CL_att_adj], mClusterCounts.nTotal);
3263 PrintClusterCount(mode, num, "Removed (Strategy A)", attachClusterCounts[CL_att_adj] - attachClusterCounts[CL_prot], mClusterCounts.nTotal); // Attached + Adjacent (also fake) - protected
3264 PrintClusterCount(mode, num, "Unaccessible", mClusterCounts.nUnaccessible, mClusterCounts.nTotal); // No contribution from track >= 10 MeV, unattached or fake-attached/adjacent
3265 } else {
3266 PrintClusterCount(mode, num, "All Clusters", mClusterCounts.nTotal, mClusterCounts.nTotal);
3267 PrintClusterCount(mode, num, "Used in Physics", mClusterCounts.nPhysics, mClusterCounts.nTotal);
3268 PrintClusterCount(mode, num, "Protected", mClusterCounts.nProt, mClusterCounts.nTotal);
3269 PrintClusterCount(mode, num, "Unattached", mClusterCounts.nUnattached, mClusterCounts.nTotal);
3270 PrintClusterCount(mode, num, "Removed (Strategy A)", mClusterCounts.nTotal - mClusterCounts.nUnattached - mClusterCounts.nProt, mClusterCounts.nTotal);
3271 PrintClusterCount(mode, num, "Removed (Strategy B)", mClusterCounts.nTotal - mClusterCounts.nProt, mClusterCounts.nTotal);
3272 }
3273
3274 PrintClusterCount(mode, num, "Merged Loopers (Track Merging)", mClusterCounts.nMergedLooperConnected, mClusterCounts.nTotal);
3275 PrintClusterCount(mode, num, "Merged Loopers (Afterburner)", mClusterCounts.nMergedLooperUnconnected, mClusterCounts.nTotal);
3276 PrintClusterCount(mode, num, "Looping Legs (other)", mClusterCounts.nLoopers, mClusterCounts.nTotal);
3277 PrintClusterCount(mode, num, "High Inclination Angle", mClusterCounts.nHighIncl, mClusterCounts.nTotal);
3278 PrintClusterCount(mode, num, "Rejected", mClusterCounts.nRejected, mClusterCounts.nTotal);
3279 PrintClusterCount(mode, num, "Tube (> 200 MeV)", mClusterCounts.nTube, mClusterCounts.nTotal);
3280 PrintClusterCount(mode, num, "Tube (< 200 MeV)", mClusterCounts.nTube200, mClusterCounts.nTotal);
3281 PrintClusterCount(mode, num, "Low Pt < 50 MeV", mClusterCounts.nLowPt, mClusterCounts.nTotal);
3282 PrintClusterCount(mode, num, "Low Pt < 200 MeV", mClusterCounts.n200MeV, mClusterCounts.nTotal);
3283
3284 if (mcPresent() && (mQATasks & taskClusterAttach)) {
3285 PrintClusterCount(mode, num, "Tracks > 400 MeV", mClusterCounts.nAbove400, mClusterCounts.nTotal);
3286 PrintClusterCount(mode, num, "Fake Removed (> 400 MeV)", mClusterCounts.nFakeRemove400, mClusterCounts.nAbove400);
3287 PrintClusterCount(mode, num, "Full Fake Removed (> 400 MeV)", mClusterCounts.nFullFakeRemove400, mClusterCounts.nAbove400);
3288 PrintClusterCount(mode, num, "Tracks < 40 MeV", mClusterCounts.nBelow40, mClusterCounts.nTotal);
3289 PrintClusterCount(mode, num, "Fake Protect (< 40 MeV)", mClusterCounts.nFakeProtect40, mClusterCounts.nBelow40);
3290 }
3291 if (mcPresent() && (mQATasks & taskTrackStatistics)) {
3292 PrintClusterCount(mode, num, "Correctly Attached all-trk normalized", mClusterCounts.nCorrectlyAttachedNormalized, mClusterCounts.nTotal);
3293 PrintClusterCount(mode, num, "Correctly Attached non-fake normalized", mClusterCounts.nCorrectlyAttachedNormalizedNonFake, mClusterCounts.nTotal);
3294 }
3295 if (mTextDump) {
3296 fclose(mTextDump);
3297 mTextDump = nullptr;
3298 }
3299 return num;
3300}
3301
3302void* GPUQA::AllocateScratchBuffer(size_t nBytes)
3303{
3304 mTrackingScratchBuffer.resize((nBytes + sizeof(mTrackingScratchBuffer[0]) - 1) / sizeof(mTrackingScratchBuffer[0]));
3305 return mTrackingScratchBuffer.data();
3306}
A const (ready only) version of MCTruthContainer.
int16_t charge
Definition RawEventData.h:5
int16_t time
Definition RawEventData.h:4
int32_t i
#define TRACK_EXPECTED_REFERENCE_X_DEFAULT
Definition GPUQA.cxx:206
#define TRACK_EXPECTED_REFERENCE_X
Definition GPUQA.cxx:231
#define QA_DEBUG
Definition GPUQA.cxx:15
#define QA_TIMING
Definition GPUQA.cxx:16
int16_t Color_t
Definition GPUQA.h:33
uint32_t iSector
GPUChain * chain
configStandalone Light OpenGL display display_light Apply smoothing to points Enable Z buffer Highlight links Highlight seeds Highlight sector tracks Highlight final tracks Show individual Show related Highlight grid Show only individual collision Color clusters belonging to track objects Set line Enable drawing TPC data Enable drawing TOF data Enable drawing magnetic field Set field line step count Invert colors configStandalone Heavy OpenGL display display_heavy Show tracks with ITS contribution Show tracks with TOF contribution Enable propagation of loopers Show clusters at nominal x of pad row for early transformed data Split CE tracks when they cross the central Mark clusters Mark adjacent clusters Hide rejected tracks Apply filter on tracks to be displayed Separate increase X h configStandalone Camera window renderer OpenGL display display_camera Orient the camera such that the y axis is always upwards Full Screen MultiSample Anti Downsampling Anti Enable Vertical Sync Target frame or if the backend cannot limit to the refresh Use OpenGL indirect draws to reduce number of draw calls Display h configStandalone Vulkan display display_vulkan Have uniform buffers in host accessible device memory h internal
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error s
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char Directory with events to Events directory is and not inside events Show standalone event Run event generator maximum time bin of continuous for triggered for automatic continuous for automatic continuous triggered uint64_t
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char e
configStandalone Reconstruction rec_tpc multiplier for the combined cluster track error during track following area in cm for the search of for z only used if decrease removal by when above this local occupancy rowx Shift of TPC revert cluster flag edge bit distance to edge exceeds this ZOffset for secondary multiple of radius max P of electron K BB bands crossing max P of electron p BB bands crossing max P of electron d BB bands crossing max P of electron t BB bands crossing min accepted dEdxMaxTPC of the track Additive extra cluster error for Y2 if edge flag set Additive extra cluster error for Y2 if single set Additive extra cluster error for Y2 if or single set Additive extra cluster error for Z2 if or single set Extra errors to add to track when crossing depending on addErrorsCECrossing Inner portion of y range in sector that is not used in searching for extrapolated track candidates factor2YS for track merging factor2K for track merging Margin in cm of Y position when rejecting edge clusters based on uncorrected track Y Max fraction of shared clusters for tracklet Except pad from cluster finding if total number of charges in a fragment is above this Except pad from cluster finding if number of consecutive charges in a fragment is above this Threshold where a timebin is considered disabling the noisy pad check for that pad Treat signal above this minimum as saturated Smoothing factor for the exponential Highly Ionising Particle tail filter Number of timebins to use for the averaging Add additional custom track errors when crossing relative as fraction of Do N fit passes in final fit of< 0:reject at max -n hits") AddOptionRTC(rejectIFCLowRadiusCluster, uint8_t, 1, "", 0, "Reject clusters that get the IFC mask error during refit") AddOptionRTC(dEdxTruncLow, uint8_t, 2, "", 0, "Low truncation threshold, fraction of 128") AddOptionRTC(dEdxTruncHigh, uint8_t, 77, "", 0, "High truncation threshold, fraction of 128") AddOptionRTC(extrapolationTracking, int8_t, 1, "", 0, "Enable Extrapolation Tracking(prolong tracks to adjacent sectors to find short segments)") AddOptionRTC(disableRefitAttachment, uint8_t, 0, "", 0, "Bitmask to disable certain attachment steps during refit(1:attachment, 2:propagation, 4:loop following, 8:mirroring)") AddOptionRTC(rejectionStrategy, uint8_t, o2::gpu::GPUSettings::RejectionStrategyA, "", 0, "Enable rejection of TPC clusters for compression(0=no, 1=strategy A, 2=strategy B)") AddOptionRTC(mergeLoopersAfterburner, uint8_t, 1, "", 0, "Run afterburner for additional looper merging") AddOptionRTC(compressionTypeMask, uint8_t, o2::gpu::GPUSettings::CompressionFull, "", 0, "TPC Compression mode bits(1=truncate charge/width LSB, 2=differences, 4=track-model)") AddOptionRTC(compressionSortOrder, uint8_t, o2::gpu::GPUSettings::SortTime, "", 0, "Sort order of TPC compression(0=time, 1=pad, 2=Z-time-pad, 3=Z-pad-time, 4=no sorting(use incoming order))") AddOptionRTC(sigBitsCharge, uint8_t, 4, "", 0, "Number of significant bits for TPC cluster charge in compression mode 1") AddOptionRTC(sigBitsWidth, uint8_t, 3, "", 0, "Number of significant bits for TPC cluster width in compression mode 1") AddOptionRTC(dropLoopers, uint8_t, 0, "", 0, "Drop looping tracks starting from second loop") AddOptionRTC(mergerInterpolateErrors, uint8_t, 1, "", 0, "Use interpolation instead of extrapolation for chi2 based cluster rejection") AddOptionRTC(mergerInterpolateRejectAlsoOnCurrentPosition, uint8_t, 1, "", 0, "When using mergerInterpolateErrors, reject based on chi2 twice computed with interpolated and current track position starting from NDF > mergerNonInterpolateRejectMinNDF Merge tracks accross the central electrode Enable PID response Add sub threshold clusters in TPC dEdx computation OR mask of TPC flags that will reject the cluster in alternative dEdx Reject edge clusters based on uncorrected track Y during track fit Increase decrease removal by when on this number of pads from the edge h configStandalone Reconstruction rec_trd Max chi2 for TRD tracklets to be matched to a track Chi2 cut for strict matching mode Number of sigmas for ITS TPC track time error estimate Addition to search road around track prolongation along Y in cm Optimal tracklet rphi resolution in maximum chi2 per attached tracklet for TRD tracks take into account the full likelihood of the slope including exponential tails Tracks with less attached tracklets are discarded after the tracking Post trigger Pile up integration time in BCs fits all iterations but the last one in projections Use material lookup table for TPC refit Apply workarounds to avoid cyclic merge should not be needed tpc dyn configStandalone Processing proc_rtc Replace constant variables static by constexpr expressions Compile RTC in deterministic with NO_FAST_MATH flags and GPUCA_DETERMINISTIC_MODE define Use RTC to optimize GPU code h configStandalone Processing proc_rtctech Use a file lock to serialize access to the cache folder int32_t
configStandalone Reconstruction rec_tpc float
std::vector< o2::MCCompLabel > labels
std::vector< o2::its::TrackITS > tracks
Definition of the MCTrack class.
Definition of the Names Generator class.
uint16_t pos
Definition RawData.h:3
uint32_t j
Definition RawData.h:0
uint32_t side
Definition RawData.h:0
uint16_t pid
Definition RawData.h:2
uint32_t c
Definition RawData.h:2
POD correction map.
TBranch * ptr
int nClusters
double num
Class for time synchronization of RawReader instances.
static constexpr int32_t NSECTORS
Definition GPUChain.h:58
int32_t ReadO2MCData(const char *filename)
Definition GPUQA.h:55
bool clusterRemovable(int32_t attach, bool prot) const
Definition GPUQA.h:53
~GPUQA()=default
Definition GPUQA.cxx:338
void * AllocateScratchBuffer(size_t nBytes)
Definition GPUQA.h:56
void SetMCTrackRange(int32_t min, int32_t max)
Definition GPUQA.h:48
mcLabelI_t GetMCTrackLabel(uint32_t trackId) const
Definition GPUQA.h:52
int32_t DrawQAHistograms()
Definition GPUQA.h:47
int32_t InitQA(int32_t tasks=0)
Definition GPUQA.h:45
void DumpO2MCData(const char *filename) const
Definition GPUQA.h:54
int32_t mcLabelI_t
Definition GPUQA.h:44
void RunQA(bool matchOnly=false)
Definition GPUQA.h:46
GPUQA(void *chain)
Definition GPUQA.h:42
static DigitizationContext * loadFromFile(std::string_view filename="")
GLdouble n
Definition glcorearb.h:1982
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
const GLfloat * m
Definition glcorearb.h:4066
GLuint GLfloat GLfloat GLfloat GLfloat y1
Definition glcorearb.h:5034
GLsizeiptr size
Definition glcorearb.h:659
GLuint GLuint end
Definition glcorearb.h:469
const GLdouble * v
Definition glcorearb.h:832
GLuint const GLchar * name
Definition glcorearb.h:781
GLdouble f
Definition glcorearb.h:310
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
GLenum GLint * range
Definition glcorearb.h:1899
GLint y
Definition glcorearb.h:270
GLenum GLenum dst
Definition glcorearb.h:1767
GLuint GLsizei const GLchar * label
Definition glcorearb.h:2519
GLuint GLfloat * val
Definition glcorearb.h:1582
GLboolean r
Definition glcorearb.h:1233
GLenum GLfloat param
Definition glcorearb.h:271
GLuint id
Definition glcorearb.h:650
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
float float float y2
Definition MathUtils.h:49
o2::track::TrackParCov int int int float int nCl
const bool const int TrackITSInternal< NLayers > & track
Enum< T >::Iterator begin(Enum< T >)
Definition Defs.h:158
value_T f3
Definition TrackUtils.h:93
value_T f1
Definition TrackUtils.h:91
value_T f2
Definition TrackUtils.h:92
std::string to_string(gsl::span< T, Size > span)
Definition common.h:52
std::string filename()
bool isValid(std::string alias)
int64_t differenceInBC(const InteractionRecord &other) const
std::tuple< std::vector< std::unique_ptr< TCanvas > >, std::vector< std::unique_ptr< TLegend > >, std::vector< std::unique_ptr< TPad > >, std::vector< std::unique_ptr< TLatex > >, std::vector< std::unique_ptr< TH1D > > > v
Definition GPUQA.cxx:316
IR getFirstIRofTF(const IR &rec) const
get 1st IR of TF corresponding to the 1st sampled orbit (in MC)
Definition HBFUtils.h:71
IR getFirstIR() const
Definition HBFUtils.h:47
unsigned int nClusters[constants::MAXSECTOR][constants::MAXGLOBALPADROW]
unsigned int clusterOffset[constants::MAXSECTOR][constants::MAXGLOBALPADROW]
const ClusterNative * clustersLinear
constexpr size_t min
constexpr size_t max
o2::InteractionRecord ir(0, 0)
vec clear()
o2::InteractionRecord ir0(3, 5)