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