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