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