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