23#include "TGraphAsymmErrors.h"
41#include "GPUChainTrackingGetters.inc"
50#include "GPUParam.inc"
70#include "TParticlePDG.h"
71#include "TDatabasePDG.h"
80#include <oneapi/tbb.h>
84#ifdef GPUCA_MERGER_BY_MC_LABEL
85#define CHECK_CLUSTER_STATE_INIT_LEG_BY_MC() \
86 if (!unattached && mTrackMCLabels[id].isValid()) { \
87 int32_t mcLabel = mTrackMCLabels[id].getTrackID(); \
88 int32_t mcEvent = mTrackMCLabels[id].getEventID(); \
89 int32_t mcSource = mTrackMCLabels[id].getSourceID(); \
90 if (mTrackMCLabelsReverse[mMCEventOffset[mcSource] + mcEvent][mcLabel] != id) { \
91 attach &= (~gputpcgmmergertypes::attachGoodLeg); \
95#define CHECK_CLUSTER_STATE_INIT_LEG_BY_MC()
98#define CHECK_CLUSTER_STATE_INIT() \
99 bool unattached = attach == 0; \
101 bool lowPt = false; \
102 bool mev200 = false; \
103 bool mergedLooper = false; \
104 int32_t id = attach & gputpcgmmergertypes::attachTrackMask; \
106 qpt = fabsf(mTracking->mIOPtrs.mergedTracks[id].GetParam().GetQPt()); \
107 lowPt = qpt * mTracking->GetParam().qptB5Scaler > mTracking->GetParam().rec.tpc.rejectQPtB5; \
109 mergedLooper = mTracking->mIOPtrs.mergedTracks[id].MergedLooper(); \
111 bool physics = false, protect = false; \
112 CHECK_CLUSTER_STATE_INIT_LEG_BY_MC();
114#define CHECK_CLUSTER_STATE() \
115 CHECK_CLUSTER_STATE_INIT() \
117 mClusterCounts.n200MeV++; \
120 mClusterCounts.nLowPt++; \
121 } else if (mergedLooper) { \
122 mClusterCounts.nMergedLooper++; \
124 GPUTPCClusterRejection::GetProtectionStatus<true>(attach, physics, protect, &mClusterCounts, &mev200); \
127#define CHECK_CLUSTER_STATE_NOCOUNT() \
128 CHECK_CLUSTER_STATE_INIT() \
130 if (!lowPt && !mergedLooper) { \
131 GPUTPCClusterRejection::GetProtectionStatus<false>(attach, physics, protect); \
136 static GPUSettingsQA defaultConfig;
138 return *
chain->mConfigQA;
140 return defaultConfig;
146static const constexpr bool PERF_FIGURE = 0;
149static const constexpr float LOG_PT_MIN = -1.;
151static constexpr float Y_MAX = 40;
152static constexpr float Z_MAX = 100;
155static constexpr float PT_MIN_PRIM = 0.1;
157static constexpr float PT_MAX = 20;
158static constexpr float ETA_MAX = 1.5;
159static constexpr float ETA_MAX2 = 0.9;
161static constexpr bool CLUST_HIST_INT_SUM =
false;
163static constexpr const int32_t COLORCOUNT = 12;
165static const constexpr char* EFF_TYPES[5] = {
"Rec",
"Clone",
"Fake",
"All",
"RecAndClone"};
166static const constexpr char* FINDABLE_NAMES[2] = {
"All",
"Findable"};
167static const constexpr char* PRIM_NAMES[2] = {
"Prim",
"Sec"};
168static const constexpr char* PARAMETER_NAMES[5] = {
"Y",
"Z",
"#Phi",
"#lambda",
"Relative #it{p}_{T}"};
169static const constexpr char* PARAMETER_NAMES_NATIVE[5] = {
"Y",
"Z",
"sin(#Phi)",
"tan(#lambda)",
"q/#it{p}_{T} (curvature)"};
170static const constexpr char* VSPARAMETER_NAMES[6] = {
"Y",
"Z",
"Phi",
"Eta",
"Pt",
"Pt_log"};
171static const constexpr char* EFF_NAMES[3] = {
"Efficiency",
"Clone Rate",
"Fake Rate"};
172static const constexpr char* EFFICIENCY_TITLES[4] = {
"Efficiency (Primary Tracks, Findable)",
"Efficiency (Secondary Tracks, Findable)",
"Efficiency (Primary Tracks)",
"Efficiency (Secondary Tracks)"};
173static const constexpr double SCALE[5] = {10., 10., 1000., 1000., 100.};
174static const constexpr double SCALE_NATIVE[5] = {10., 10., 1000., 1000., 1.};
175static const constexpr char* XAXIS_TITLES[5] = {
"#it{y}_{mc} (cm)",
"#it{z}_{mc} (cm)",
"#Phi_{mc} (rad)",
"#eta_{mc}",
"#it{p}_{Tmc} (GeV/#it{c})"};
176static const constexpr char* AXIS_TITLES[5] = {
"#it{y}-#it{y}_{mc} (mm) (Resolution)",
"#it{z}-#it{z}_{mc} (mm) (Resolution)",
"#phi-#phi_{mc} (mrad) (Resolution)",
"#lambda-#lambda_{mc} (mrad) (Resolution)",
"(#it{p}_{T} - #it{p}_{Tmc}) / #it{p}_{Tmc} (%) (Resolution)"};
177static const constexpr char* AXIS_TITLES_NATIVE[5] = {
"#it{y}-#it{y}_{mc} (mm) (Resolution)",
"#it{z}-#it{z}_{mc} (mm) (Resolution)",
"sin(#phi)-sin(#phi_{mc}) (Resolution)",
"tan(#lambda)-tan(#lambda_{mc}) (Resolution)",
"q*(q/#it{p}_{T} - q/#it{p}_{Tmc}) (Resolution)"};
178static const constexpr char* AXIS_TITLES_PULL[5] = {
"#it{y}-#it{y}_{mc}/#sigma_{y} (Pull)",
"#it{z}-#it{z}_{mc}/#sigma_{z} (Pull)",
"sin(#phi)-sin(#phi_{mc})/#sigma_{sin(#phi)} (Pull)",
"tan(#lambda)-tan(#lambda_{mc})/#sigma_{tan(#lambda)} (Pull)",
179 "q*(q/#it{p}_{T} - q/#it{p}_{Tmc})/#sigma_{q/#it{p}_{T}} (Pull)"};
180static const constexpr char* CLUSTER_NAMES[GPUQA::N_CLS_HIST] = {
"Correctly attached clusters",
"Fake attached clusters",
"Attached + adjacent clusters",
"Fake adjacent clusters",
"Clusters of reconstructed tracks",
"Used in Physics",
"Protected",
"All clusters"};
181static const constexpr char* CLUSTER_TITLES[GPUQA::N_CLS_TYPE] = {
"Clusters Pt Distribution / Attachment",
"Clusters Pt Distribution / Attachment (relative to all clusters)",
"Clusters Pt Distribution / Attachment (integrated)"};
182static const constexpr char* CLUSTER_NAMES_SHORT[GPUQA::N_CLS_HIST] = {
"Attached",
"Fake",
"AttachAdjacent",
"FakeAdjacent",
"FoundTracks",
"Physics",
"Protected",
"All"};
183static const constexpr char* CLUSTER_TYPES[GPUQA::N_CLS_TYPE] = {
"",
"Ratio",
"Integral"};
184static const constexpr int32_t COLORS_HEX[COLORCOUNT] = {0xB03030, 0x00A000, 0x0000C0, 0x9400D3, 0x19BBBF, 0xF25900, 0x7F7F7F, 0xFFD700, 0x07F707, 0x07F7F7, 0xF08080, 0x000000};
186static const constexpr int32_t CONFIG_DASHED_MARKERS = 0;
188static const constexpr float AXES_MIN[5] = {-Y_MAX, -Z_MAX, 0.f, -ETA_MAX, PT_MIN};
189static const constexpr float AXES_MAX[5] = {Y_MAX, Z_MAX, 2.f * M_PI, ETA_MAX, PT_MAX};
190static const constexpr int32_t AXIS_BINS[5] = {51, 51, 144, 31, 50};
191static const constexpr int32_t RES_AXIS_BINS[] = {1017, 113};
192static const constexpr float RES_AXES[5] = {1., 1., 0.03, 0.03, 1.0};
193static const constexpr float RES_AXES_NATIVE[5] = {1., 1., 0.1, 0.1, 5.0};
194static const constexpr float PULL_AXIS = 10.f;
196std::vector<TColor*> GPUQA::mColors;
197int32_t GPUQA::initColors()
199 mColors.reserve(COLORCOUNT);
200 for (int32_t
i = 0;
i < COLORCOUNT;
i++) {
201 float f1 = (float)((COLORS_HEX[
i] >> 16) & 0xFF) / (
float)0xFF;
202 float f2 = (float)((COLORS_HEX[
i] >> 8) & 0xFF) / (
float)0xFF;
203 float f3 = (float)((COLORS_HEX[
i] >> 0) & 0xFF) / (
float)0xFF;
204 mColors.emplace_back(
new TColor(10000 +
i, f1, f2, f3));
208static constexpr Color_t defaultColorNums[COLORCOUNT] = {kRed, kBlue, kGreen, kMagenta, kOrange, kAzure, kBlack, kYellow, kGray, kTeal, kSpring, kPink};
210#define TRACK_EXPECTED_REFERENCE_X_DEFAULT 81
211#ifdef GPUCA_TPC_GEOMETRY_O2
212static inline int32_t GPUQA_O2_ConvertFakeLabel(int32_t
label) {
return label >= 0x7FFFFFFE ? -1 :
label; }
213inline uint32_t GPUQA::GetNMCCollissions()
const {
return mMCInfosCol.size(); }
214inline uint32_t GPUQA::GetNMCTracks(int32_t iCol)
const {
return mMCInfosCol[iCol].num; }
215inline uint32_t GPUQA::GetNMCTracks(
const mcLabelI_t&
label)
const {
return mMCInfosCol[mMCEventOffset[
label.getSourceID()] +
label.getEventID()].num; }
216inline uint32_t GPUQA::GetNMCLabels()
const {
return mClNative->clustersMCTruth ? mClNative->clustersMCTruth->getIndexedSize() : 0; }
217inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(uint32_t iTrk, uint32_t iCol) {
return mMCInfos[mMCInfosCol[iCol].first + iTrk]; }
218inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(
const mcLabel_t&
label) {
return mMCInfos[mMCInfosCol[mMCEventOffset[
label.getSourceID()] +
label.getEventID()].first +
label.getTrackID()]; }
219inline GPUQA::mcLabels_t GPUQA::GetMCLabel(uint32_t
i) {
return mClNative->clustersMCTruth->getLabels(
i); }
220inline int32_t GPUQA::GetMCLabelNID(
const mcLabels_t&
label) {
return label.size(); }
221inline int32_t GPUQA::GetMCLabelNID(uint32_t
i) {
return mClNative->clustersMCTruth->getLabels(
i).size(); }
222inline GPUQA::mcLabel_t GPUQA::GetMCLabel(uint32_t
i, uint32_t
j) {
return mClNative->clustersMCTruth->getLabels(
i)[
j]; }
223inline int32_t GPUQA::GetMCLabelID(uint32_t
i, uint32_t
j) {
return GPUQA_O2_ConvertFakeLabel(mClNative->clustersMCTruth->getLabels(
i)[
j].getTrackID()); }
224inline int32_t GPUQA::GetMCLabelID(
const mcLabels_t&
label, uint32_t
j) {
return GPUQA_O2_ConvertFakeLabel(
label[
j].getTrackID()); }
225inline int32_t GPUQA::GetMCLabelID(
const mcLabel_t&
label) {
return GPUQA_O2_ConvertFakeLabel(
label.getTrackID()); }
226inline uint32_t GPUQA::GetMCLabelCol(uint32_t
i, uint32_t
j) {
return mMCEventOffset[mClNative->clustersMCTruth->getLabels(
i)[
j].getSourceID()] + mClNative->clustersMCTruth->getLabels(
i)[
j].getEventID(); }
227inline const auto& GPUQA::GetClusterLabels() {
return mClNative->clustersMCTruth; }
228inline float GPUQA::GetMCLabelWeight(uint32_t
i, uint32_t
j) {
return 1; }
229inline float GPUQA::GetMCLabelWeight(
const mcLabels_t&
label, uint32_t
j) {
return 1; }
230inline float GPUQA::GetMCLabelWeight(
const mcLabel_t&
label) {
return 1; }
231inline bool GPUQA::mcPresent() {
return !mConfig.noMC && mTracking && mClNative && mClNative->clustersMCTruth && mMCInfos.size(); }
232uint32_t GPUQA::GetMCLabelCol(
const mcLabel_t&
label)
const {
return !
label.isValid() ? 0 : (mMCEventOffset[
label.getSourceID()] +
label.getEventID()); }
234bool GPUQA::CompareIgnoreFake(
const mcLabelI_t& l1,
const mcLabelI_t& l2) {
return l1.compare(l2) >= 0; }
235#define TRACK_EXPECTED_REFERENCE_X 78
237inline GPUQA::mcLabelI_t::mcLabelI_t(
const GPUQA::mcLabel_t& l) : track(l.fMCID) {}
238inline bool GPUQA::mcLabelI_t::operator==(
const GPUQA::mcLabel_t& l) {
return AbsLabelID(track) == l.fMCID; }
239inline uint32_t GPUQA::GetNMCCollissions()
const {
return 1; }
240inline uint32_t GPUQA::GetNMCTracks(int32_t iCol)
const {
return mTracking->mIOPtrs.nMCInfosTPC; }
241inline uint32_t GPUQA::GetNMCTracks(
const mcLabelI_t&
label)
const {
return mTracking->mIOPtrs.nMCInfosTPC; }
242inline uint32_t GPUQA::GetNMCLabels()
const {
return mTracking->mIOPtrs.nMCLabelsTPC; }
243inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(uint32_t iTrk, uint32_t iCol) {
return mTracking->mIOPtrs.mcInfosTPC[AbsLabelID(iTrk)]; }
244inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(
const mcLabel_t&
label) {
return GetMCTrack(
label.fMCID, 0); }
245inline const GPUQA::mcInfo_t& GPUQA::GetMCTrack(
const mcLabelI_t&
label) {
return GetMCTrack(
label.track, 0); }
246inline const GPUQA::mcLabels_t& GPUQA::GetMCLabel(uint32_t
i) {
return mTracking->mIOPtrs.mcLabelsTPC[
i]; }
247inline const GPUQA::mcLabel_t& GPUQA::GetMCLabel(uint32_t
i, uint32_t
j) {
return mTracking->mIOPtrs.mcLabelsTPC[
i].fClusterID[
j]; }
248inline int32_t GPUQA::GetMCLabelNID(
const mcLabels_t&
label) {
return 3; }
249inline int32_t GPUQA::GetMCLabelNID(uint32_t
i) {
return 3; }
250inline int32_t GPUQA::GetMCLabelID(uint32_t
i, uint32_t
j) {
return mTracking->mIOPtrs.mcLabelsTPC[
i].fClusterID[
j].fMCID; }
251inline int32_t GPUQA::GetMCLabelID(
const mcLabels_t&
label, uint32_t
j) {
return label.fClusterID[
j].fMCID; }
252inline int32_t GPUQA::GetMCLabelID(
const mcLabel_t&
label) {
return label.fMCID; }
253inline uint32_t GPUQA::GetMCLabelCol(uint32_t
i, uint32_t
j) {
return 0; }
255inline const auto& GPUQA::GetClusterLabels() {
return mTracking->mIOPtrs.mcLabelsTPC; }
256inline float GPUQA::GetMCLabelWeight(uint32_t
i, uint32_t
j) {
return mTracking->mIOPtrs.mcLabelsTPC[
i].fClusterID[
j].fWeight; }
257inline float GPUQA::GetMCLabelWeight(
const mcLabels_t&
label, uint32_t
j) {
return label.fClusterID[
j].fWeight; }
258inline float GPUQA::GetMCLabelWeight(
const mcLabel_t&
label) {
return label.fWeight; }
259inline int32_t GPUQA::FakeLabelID(int32_t
id) {
return id < 0 ?
id : (-2 -
id); }
260inline int32_t GPUQA::AbsLabelID(int32_t
id) {
return id >= 0 ?
id : (-
id - 2); }
261inline bool GPUQA::mcPresent() {
return !mConfig.noMC && mTracking && GetNMCLabels() && GetNMCTracks(0); }
262uint32_t GPUQA::GetMCLabelCol(
const mcLabel_t&
label)
const {
return 0; }
264bool GPUQA::CompareIgnoreFake(
const mcLabelI_t& l1,
const mcLabelI_t& l2) {
return AbsLabelID(l1) == AbsLabelID(l2); }
265#define TRACK_EXPECTED_REFERENCE_X TRACK_EXPECTED_REFERENCE_X_DEFAULT
270 return obj[mMCEventOffset[l.getSourceID()] + l.getEventID()][l.getTrackID()];
274auto GPUQA::getHistArray<TH1F>()
276 return std::make_pair(mHist1D, &mHist1D_pos);
279auto GPUQA::getHistArray<TH2F>()
281 return std::make_pair(mHist2D, &mHist2D_pos);
284auto GPUQA::getHistArray<TH1D>()
286 return std::make_pair(mHist1Dd, &mHist1Dd_pos);
289auto GPUQA::getHistArray<TGraphAsymmErrors>()
291 return std::make_pair(mHistGraph, &mHistGraph_pos);
293template <
class T,
typename... Args>
294void GPUQA::createHist(T*&
h,
const char*
name, Args... args)
296 const auto&
p = getHistArray<T>();
297 if (mHaveExternalHists) {
298 if (
p.first->size() <=
p.second->size()) {
299 GPUError(
"Array sizes mismatch: Histograms %lu <= Positions %lu",
p.first->size(),
p.second->size());
300 throw std::runtime_error(
"Incoming histogram array incomplete");
302 if (strcmp((*
p.first)[
p.second->size()].GetName(),
name)) {
303 GPUError(
"Histogram name mismatch: in array %s, trying to create %s", (*
p.first)[
p.second->size()].GetName(),
name);
304 throw std::runtime_error(
"Incoming histogram has incorrect name");
307 if constexpr (std::is_same_v<T, TGraphAsymmErrors>) {
308 p.first->emplace_back();
309 p.first->back().SetName(
name);
311 p.first->emplace_back(
name, args...);
314 h = &((*
p.first)[
p.second->size()]);
315 p.second->emplace_back(&
h);
321 std::tuple<std::vector<std::unique_ptr<TCanvas>>, std::vector<std::unique_ptr<TLegend>>, std::vector<std::unique_ptr<TPad>>, std::vector<std::unique_ptr<TLatex>>, std::vector<std::unique_ptr<TH1D>>>
v;
325template <
class T,
typename... Args>
326T* GPUQA::createGarbageCollected(Args... args)
328 auto&
v = std::get<std::vector<std::unique_ptr<T>>>(mGarbageCollector->v);
329 v.emplace_back(std::make_unique<T>(args...));
330 return v.back().get();
332void GPUQA::clearGarbagageCollector()
334 std::get<std::vector<std::unique_ptr<TPad>>>(mGarbageCollector->v).
clear();
335 std::apply([](
auto&&... args) { ((args.clear()), ...); }, mGarbageCollector->v);
340 mMCEventOffset.resize(1, 0);
345 if (mQAInitialized && !mHaveExternalHists) {
351 clearGarbagageCollector();
358 return protect || physics;
360 return (!unattached && !physics && !protect);
364void GPUQA::SetAxisSize(T* e)
366 e->GetYaxis()->SetTitleOffset(1.0);
367 e->GetYaxis()->SetTitleSize(0.045);
368 e->GetYaxis()->SetLabelSize(0.045);
369 e->GetXaxis()->SetTitleOffset(1.03);
370 e->GetXaxis()->SetTitleSize(0.045);
371 e->GetXaxis()->SetLabelOffset(-0.005);
372 e->GetXaxis()->SetLabelSize(0.045);
375void GPUQA::SetLegend(TLegend* l)
378 l->SetTextSize(0.016);
382double* GPUQA::CreateLogAxis(int32_t nbins,
float xmin,
float xmax)
384 float logxmin = std::log10(xmin);
385 float logxmax = std::log10(xmax);
386 float binwidth = (logxmax - logxmin) / nbins;
388 double* xbins =
new double[nbins + 1];
391 for (int32_t
i = 1;
i <= nbins;
i++) {
392 xbins[
i] = std::pow(10, logxmin +
i * binwidth);
397void GPUQA::ChangePadTitleSize(TPad* p,
float size)
400 TPaveText* pt = (TPaveText*)(
p->GetPrimitive(
"title"));
402 GPUError(
"Error changing title");
404 pt->SetTextSize(
size);
409void GPUQA::DrawHisto(TH1* histo,
char*
filename,
char* options)
413 histo->Draw(options);
417void GPUQA::doPerfFigure(
float x,
float y,
float size)
422 TLatex* t = createGarbageCollected<TLatex>();
425 t->SetTextSize(
size);
426 t->DrawLatex(
x,
y, str_perf_figure_1);
427 t->DrawLatex(
x,
y - 0.01 -
size, str_perf_figure_2);
436int32_t GPUQA::InitQACreateHistograms()
438 char name[2048], fname[1024];
439 if (mQATasks & taskTrackingEff) {
441 for (int32_t
i = 0;
i < 5;
i++) {
442 for (int32_t
j = 0;
j < 2;
j++) {
443 for (int32_t k = 0; k < 2; k++) {
444 for (int32_t l = 0; l < 5; l++) {
445 snprintf(
name, 2048,
"%s%s%s%sVs%s",
"tracks", EFF_TYPES[
i], FINDABLE_NAMES[
j], PRIM_NAMES[k], VSPARAMETER_NAMES[l]);
447 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], k == 0 ? PT_MIN_PRIM : AXES_MIN[4], AXES_MAX[4])};
448 createHist(mEff[
i][
j][k][l],
name,
name, AXIS_BINS[l], binsPt.get());
450 createHist(mEff[
i][
j][k][l],
name,
name, AXIS_BINS[l], AXES_MIN[l], AXES_MAX[l]);
452 if (!mHaveExternalHists) {
453 mEff[
i][
j][k][l]->Sumw2();
455 strcat(
name,
"_eff");
457 createHist(mEffResult[
i][
j][k][l],
name);
466 if (mQATasks & taskTrackingRes) {
467 for (int32_t
i = 0;
i < 5;
i++) {
468 for (int32_t
j = 0;
j < 5;
j++) {
469 snprintf(
name, 2048,
"rms_%s_vs_%s", VSPARAMETER_NAMES[
i], VSPARAMETER_NAMES[
j]);
470 snprintf(fname, 1024,
"mean_%s_vs_%s", VSPARAMETER_NAMES[
i], VSPARAMETER_NAMES[
j]);
472 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], mConfig.resPrimaries == 1 ? PT_MIN_PRIM : AXES_MIN[4], AXES_MAX[4])};
473 createHist(mRes[
i][
j][0],
name,
name, AXIS_BINS[
j], binsPt.get());
474 createHist(mRes[
i][
j][1], fname, fname, AXIS_BINS[
j], binsPt.get());
476 createHist(mRes[
i][
j][0],
name,
name, AXIS_BINS[
j], AXES_MIN[
j], AXES_MAX[
j]);
477 createHist(mRes[
i][
j][1], fname, fname, AXIS_BINS[
j], AXES_MIN[
j], AXES_MAX[
j]);
479 snprintf(
name, 2048,
"res_%s_vs_%s", VSPARAMETER_NAMES[
i], VSPARAMETER_NAMES[
j]);
480 const float* axis = mConfig.nativeFitResolutions ? RES_AXES_NATIVE : RES_AXES;
481 const int32_t nbins =
i == 4 && mConfig.nativeFitResolutions ? (10 * RES_AXIS_BINS[0]) : RES_AXIS_BINS[0];
483 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], mConfig.resPrimaries == 1 ? PT_MIN_PRIM : AXES_MIN[4], AXES_MAX[4])};
484 createHist(mRes2[
i][
j],
name,
name, nbins, -axis[
i], axis[
i], AXIS_BINS[
j], binsPt.get());
486 createHist(mRes2[
i][
j],
name,
name, nbins, -axis[
i], axis[
i], AXIS_BINS[
j], AXES_MIN[
j], AXES_MAX[
j]);
493 if (mQATasks & taskTrackingResPull) {
494 for (int32_t
i = 0;
i < 5;
i++) {
495 for (int32_t
j = 0;
j < 5;
j++) {
496 snprintf(
name, 2048,
"pull_rms_%s_vs_%s", VSPARAMETER_NAMES[
i], VSPARAMETER_NAMES[
j]);
497 snprintf(fname, 1024,
"pull_mean_%s_vs_%s", VSPARAMETER_NAMES[
i], VSPARAMETER_NAMES[
j]);
499 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], AXES_MIN[4], AXES_MAX[4])};
500 createHist(mPull[
i][
j][0],
name,
name, AXIS_BINS[
j], binsPt.get());
501 createHist(mPull[
i][
j][1], fname, fname, AXIS_BINS[
j], binsPt.get());
503 createHist(mPull[
i][
j][0],
name,
name, AXIS_BINS[
j], AXES_MIN[
j], AXES_MAX[
j]);
504 createHist(mPull[
i][
j][1], fname, fname, AXIS_BINS[
j], AXES_MIN[
j], AXES_MAX[
j]);
506 snprintf(
name, 2048,
"pull_%s_vs_%s", VSPARAMETER_NAMES[
i], VSPARAMETER_NAMES[
j]);
508 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], AXES_MIN[4], AXES_MAX[4])};
509 createHist(mPull2[
i][
j],
name,
name, RES_AXIS_BINS[0], -PULL_AXIS, PULL_AXIS, AXIS_BINS[
j], binsPt.get());
511 createHist(mPull2[
i][
j],
name,
name, RES_AXIS_BINS[0], -PULL_AXIS, PULL_AXIS, AXIS_BINS[
j], AXES_MIN[
j], AXES_MAX[
j]);
518 if (mQATasks & taskClusterAttach) {
519 for (int32_t
i = 0;
i < N_CLS_TYPE * N_CLS_HIST - 1;
i++) {
520 int32_t ioffset =
i >= (2 * N_CLS_HIST - 1) ? (2 * N_CLS_HIST - 1) :
i >= N_CLS_HIST ? N_CLS_HIST : 0;
521 int32_t itype =
i >= (2 * N_CLS_HIST - 1) ? 2 :
i >= N_CLS_HIST ? 1 : 0;
522 snprintf(
name, 2048,
"clusters%s%s", CLUSTER_NAMES_SHORT[
i - ioffset], CLUSTER_TYPES[itype]);
523 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], PT_MIN_CLUST, PT_MAX)};
524 createHist(mClusters[
i],
name,
name, AXIS_BINS[4], binsPt.get());
528 if (mQATasks & taskTrackStatistics) {
530 for (int32_t
i = 0;
i < 2;
i++) {
531 snprintf(
name, 2048,
i ?
"nrows_with_cluster" :
"nclusters");
532 createHist(mNCl[
i],
name,
name, 160, 0, 159);
534 snprintf(
name, 2048,
"tracks");
535 std::unique_ptr<double[]> binsPt{CreateLogAxis(AXIS_BINS[4], PT_MIN_CLUST, PT_MAX)};
536 createHist(mTracks,
name,
name, AXIS_BINS[4], binsPt.get());
537 createHist(mClXY,
"clXY",
"clXY", 1000, -250, 250, 1000, -250, 250);
540 if ((mQATasks & taskClusterCounts) && mConfig.clusterRejectionHistograms) {
541 int32_t
num = DoClusterCounts(
nullptr, 2);
542 mHistClusterCount.resize(
num);
543 DoClusterCounts(
nullptr, 1);
546 for (uint32_t
i = 0;
i < mHist1D->size();
i++) {
547 *mHist1D_pos[
i] = &(*mHist1D)[
i];
549 for (uint32_t
i = 0;
i < mHist2D->size();
i++) {
550 *mHist2D_pos[
i] = &(*mHist2D)[
i];
552 for (uint32_t
i = 0;
i < mHist1Dd->size();
i++) {
553 *mHist1Dd_pos[
i] = &(*mHist1Dd)[
i];
555 for (uint32_t
i = 0;
i < mHistGraph->size();
i++) {
556 *mHistGraph_pos[
i] = &(*mHistGraph)[
i];
562int32_t GPUQA::loadHistograms(std::vector<TH1F>& i1, std::vector<TH2F>& i2, std::vector<TH1D>& i3, std::vector<TGraphAsymmErrors>& i4, int32_t tasks)
565 tasks = taskDefaultPostprocess;
567 if (mQAInitialized && (!mHaveExternalHists || tasks != mQATasks)) {
568 throw std::runtime_error(
"QA not initialized or initialized with different task array");
576 mHist1Dd_pos.clear();
577 mHistGraph_pos.clear();
578 mHaveExternalHists =
true;
583 if (InitQACreateHistograms()) {
586 mQAInitialized =
true;
596 uint32_t
n = mMCInfos.size();
597 fwrite(&
n,
sizeof(
n), 1, fp);
598 fwrite(mMCInfos.data(),
sizeof(mMCInfos[0]),
n, fp);
599 n = mMCInfosCol.size();
600 fwrite(&
n,
sizeof(
n), 1, fp);
601 fwrite(mMCInfosCol.data(),
sizeof(mMCInfosCol[0]),
n, fp);
602 n = mMCEventOffset.size();
603 fwrite(&
n,
sizeof(
n), 1, fp);
604 fwrite(mMCEventOffset.data(),
sizeof(mMCEventOffset[0]),
n, fp);
616 if ((
x = fread(&
n,
sizeof(
n), 1, fp)) != 1) {
621 if (fread(mMCInfos.data(),
sizeof(mMCInfos[0]),
n, fp) !=
n) {
625 if ((
x = fread(&
n,
sizeof(
n), 1, fp)) != 1) {
629 mMCInfosCol.resize(
n);
630 if (fread(mMCInfosCol.data(),
sizeof(mMCInfosCol[0]),
n, fp) !=
n) {
634 if ((
x = fread(&
n,
sizeof(
n), 1, fp)) != 1) {
638 mMCEventOffset.resize(
n);
639 if (fread(mMCEventOffset.data(),
sizeof(mMCEventOffset[0]),
n, fp) !=
n) {
643 if (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2) {
644 printf(
"Read %ld bytes MC Infos\n", ftell(fp));
648 CopyO2MCtoIOPtr(&mTracking->mIOPtrs);
655 ptr->mcInfosTPC = mMCInfos.data();
656 ptr->nMCInfosTPC = mMCInfos.size();
657 ptr->mcInfosTPCCol = mMCInfosCol.data();
658 ptr->nMCInfosTPCCol = mMCInfosCol.size();
664 if (!mO2MCDataLoaded) {
666 if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
667 GPUInfo(
"Start reading O2 Track MC information");
670 static constexpr float PRIM_MAX_T = 0.01f;
673 std::vector<int32_t> refId;
676 const auto& evrec = dc->getEventRecords();
677 const auto& evparts = dc->getEventParts();
678 std::vector<std::vector<float>> evTimeBins(mcReader.getNSources());
679 for (uint32_t
i = 0;
i < evTimeBins.size();
i++) {
680 evTimeBins[
i].resize(mcReader.getNEvents(
i), -100.f);
682 for (uint32_t
i = 0;
i < evrec.size();
i++) {
683 const auto&
ir = evrec[
i];
684 for (uint32_t
j = 0;
j < evparts[
i].size();
j++) {
685 const int iSim = evparts[
i][
j].sourceID;
686 const int iEv = evparts[
i][
j].entryID;
690 if (evTimeBins[iSim][iEv] >= 0) {
691 throw std::runtime_error(
"Multiple time bins for same MC collision found");
693 evTimeBins[iSim][iEv] = timebin;
698 uint32_t nSimSources = mcReader.getNSources();
699 mMCEventOffset.resize(nSimSources);
700 uint32_t nSimTotalEvents = 0;
701 uint32_t nSimTotalTracks = 0;
702 for (uint32_t
i = 0;
i < nSimSources;
i++) {
703 mMCEventOffset[
i] = nSimTotalEvents;
704 nSimTotalEvents += mcReader.getNEvents(
i);
707 mMCInfosCol.resize(nSimTotalEvents);
708 for (int32_t iSim = 0; iSim < mcReader.getNSources(); iSim++) {
709 for (int32_t
i = 0;
i < mcReader.getNEvents(iSim);
i++) {
710 const float timebin = evTimeBins[iSim][
i];
712 const std::vector<o2::MCTrack>&
tracks = mcReader.getTracks(iSim,
i);
713 const std::vector<o2::TrackReference>& trackRefs = mcReader.getTrackRefsByEvent(iSim,
i);
715 refId.resize(
tracks.size());
716 std::fill(refId.begin(), refId.end(), -1);
717 for (uint32_t
j = 0;
j < trackRefs.size();
j++) {
719 int32_t trkId = trackRefs[
j].getTrackID();
720 if (refId[trkId] == -1) {
725 mMCInfosCol[mMCEventOffset[iSim] +
i].first = mMCInfos.size();
726 mMCInfosCol[mMCEventOffset[iSim] +
i].num =
tracks.size();
727 mMCInfos.resize(mMCInfos.size() +
tracks.size());
728 for (uint32_t
j = 0;
j <
tracks.size();
j++) {
729 auto& info = mMCInfos[mMCInfosCol[mMCEventOffset[iSim] +
i].first +
j];
731 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(trk.GetPdgCode());
733 if (abs(trk.GetPdgCode()) == kElectron) {
736 if (abs(trk.GetPdgCode()) == kMuonMinus) {
739 if (abs(trk.GetPdgCode()) == kPiPlus) {
742 if (abs(trk.GetPdgCode()) == kKPlus) {
745 if (abs(trk.GetPdgCode()) == kProton) {
749 info.charge = particle ? particle->Charge() : 0;
750 info.prim = trk.T() < PRIM_MAX_T;
751 info.primDaughters = 0;
752 if (trk.getFirstDaughterTrackId() != -1) {
753 for (int32_t k = trk.getFirstDaughterTrackId(); k <= trk.getLastDaughterTrackId(); k++) {
754 if (tracks[k].
T() < PRIM_MAX_T) {
755 info.primDaughters = 1;
763 const auto& trkRef = trackRefs[refId[
j]];
767 info.pX = trkRef.Px();
768 info.pY = trkRef.Py();
769 info.pZ = trkRef.Pz();
770 info.genRadius = std::sqrt(trk.GetStartVertexCoordinatesX() * trk.GetStartVertexCoordinatesX() + trk.GetStartVertexCoordinatesY() * trk.GetStartVertexCoordinatesY() + trk.GetStartVertexCoordinatesZ() * trk.GetStartVertexCoordinatesZ());
772 info.x = info.y = info.z = info.pX = info.pY = info.pZ = 0;
778 if (mTracking && mTracking->GetProcessingSettings().debugLevel) {
781 mO2MCDataLoaded =
true;
784 CopyO2MCtoIOPtr(updateIOPtr);
791 if (mQAInitialized) {
792 throw std::runtime_error(
"QA already initialized");
798 mHist1D =
new std::vector<TH1F>;
799 mHist2D =
new std::vector<TH2F>;
800 mHist1Dd =
new std::vector<TH1D>;
801 mHistGraph =
new std::vector<TGraphAsymmErrors>;
807 if (mTracking->GetProcessingSettings().qcRunFraction != 100.f && mQATasks != taskClusterCounts) {
808 throw std::runtime_error(
"QA with qcRunFraction only supported for taskClusterCounts");
812 mClNative = mTracking->mIOPtrs.clustersNative;
815 if (InitQACreateHistograms()) {
819 if (mConfig.enableLocalOutput) {
820 mkdir(
"plots", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
825 InitO2MCData(mTracking ? &mTracking->mIOPtrs : nullptr);
829 if (mConfig.matchMCLabels.size()) {
830 uint32_t nFiles = mConfig.matchMCLabels.size();
831 std::vector<std::unique_ptr<TFile>> files;
832 std::vector<std::vector<std::vector<int32_t>>*> labelsBuffer(nFiles);
833 std::vector<std::vector<std::vector<int32_t>>*> effBuffer(nFiles);
834 for (uint32_t
i = 0;
i < nFiles;
i++) {
835 files.emplace_back(std::make_unique<TFile>(mConfig.matchMCLabels[
i].c_str()));
836 labelsBuffer[
i] = (std::vector<std::vector<int32_t>>*)files[
i]->Get(
"mcLabelBuffer");
837 effBuffer[
i] = (std::vector<std::vector<int32_t>>*)files[
i]->Get(
"mcEffBuffer");
838 if (labelsBuffer[
i] ==
nullptr || effBuffer[
i] ==
nullptr) {
839 GPUError(
"Error opening / reading from labels file %u/%s: %p %p",
i, mConfig.matchMCLabels[
i].c_str(), (
void*)labelsBuffer[
i], (
void*)effBuffer[
i]);
844 mGoodTracks.resize(labelsBuffer[0]->
size());
845 mGoodHits.resize(labelsBuffer[0]->
size());
846 for (uint32_t iEvent = 0; iEvent < labelsBuffer[0]->size(); iEvent++) {
847 std::vector<bool> labelsOK((*effBuffer[0])[iEvent].
size());
848 for (uint32_t k = 0; k < (*effBuffer[0])[iEvent].
size(); k++) {
850 for (uint32_t l = 0; l < nFiles; l++) {
851 if ((*effBuffer[0])[iEvent][k] != (*effBuffer[l])[iEvent][k]) {
857 mGoodTracks[iEvent].resize((*labelsBuffer[0])[iEvent].size());
858 for (uint32_t k = 0; k < (*labelsBuffer[0])[iEvent].
size(); k++) {
859 if ((*labelsBuffer[0])[iEvent][k] == MC_LABEL_INVALID) {
862 mGoodTracks[iEvent][k] = labelsOK[abs((*labelsBuffer[0])[iEvent][k])];
866 mQAInitialized =
true;
872 if (!mQAInitialized) {
873 throw std::runtime_error(
"QA not initialized");
875 if (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2) {
876 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));
878 if (!clNative && mTracking) {
879 clNative = mTracking->mIOPtrs.clustersNative;
881 mClNative = clNative;
883#ifdef GPUCA_TPC_GEOMETRY_O2
884 uint32_t nSimEvents = GetNMCCollissions();
885 if (mTrackMCLabelsReverse.size() < nSimEvents) {
886 mTrackMCLabelsReverse.resize(nSimEvents);
888 if (mRecTracks.size() < nSimEvents) {
889 mRecTracks.resize(nSimEvents);
891 if (mFakeTracks.size() < nSimEvents) {
892 mFakeTracks.resize(nSimEvents);
894 if (mMCParam.size() < nSimEvents) {
895 mMCParam.resize(nSimEvents);
900 uint32_t nReconstructedTracks = 0;
901 if (tracksExternal) {
903 nReconstructedTracks = tracksExternal->size();
906 nReconstructedTracks = mTracking->mIOPtrs.nMergedTracks;
908 mTrackMCLabels.resize(nReconstructedTracks);
909 for (uint32_t iCol = 0; iCol < GetNMCCollissions(); iCol++) {
910 mTrackMCLabelsReverse[iCol].resize(GetNMCTracks(iCol));
911 mRecTracks[iCol].resize(GetNMCTracks(iCol));
912 mFakeTracks[iCol].resize(GetNMCTracks(iCol));
913 mMCParam[iCol].resize(GetNMCTracks(iCol));
914 memset(mRecTracks[iCol].
data(), 0, mRecTracks[iCol].
size() *
sizeof(mRecTracks[iCol][0]));
915 memset(mFakeTracks[iCol].
data(), 0, mFakeTracks[iCol].
size() *
sizeof(mFakeTracks[iCol][0]));
916 for (
size_t i = 0;
i < mTrackMCLabelsReverse[iCol].size();
i++) {
917 mTrackMCLabelsReverse[iCol][
i] = -1;
920 if (mQATasks & taskClusterAttach && GetNMCLabels()) {
921 mClusterParam.resize(GetNMCLabels());
922 memset(mClusterParam.data(), 0, mClusterParam.size() *
sizeof(mClusterParam[0]));
927 if (mConfig.writeMCLabels) {
928 mcEffBuffer.resize(mNEvents);
929 mcLabelBuffer.resize(mNEvents);
930 mcEffBuffer[mNEvents - 1].resize(GetNMCTracks(0));
931 mcLabelBuffer[mNEvents - 1].resize(nReconstructedTracks);
934 bool mcAvail = mcPresent() || tracksExtMC;
939 if (tracksExternal) {
941 for (uint32_t
i = 0;
i < tracksExternal->size();
i++) {
942 mTrackMCLabels[
i] = (*tracksExtMC)[
i];
946 tbb::parallel_for(tbb::blocked_range<uint32_t>(0, nReconstructedTracks, (
QA_DEBUG == 0) ? 32 : nReconstructedTracks), [&](const tbb::blocked_range<uint32_t>&
range) {
947 auto acc = GPUTPCTrkLbl<true, mcLabelI_t>(GetClusterLabels(), 1.f - mConfig.recThreshold);
952 std::vector<mcLabel_t> labels;
953 for (uint32_t k = 0; k < track.NClusters(); k++) {
958 uint32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].num;
959 if (hitId >= GetNMCLabels()) {
960 GPUError(
"Invalid hit id %u > %d (nClusters %d)", hitId, GetNMCLabels(), mTracking->mIOPtrs.clustersNative ? mTracking->mIOPtrs.clustersNative->nClustersTotal : 0);
961 throw std::runtime_error(
"qa error");
964 for (int32_t
j = 0;
j < GetMCLabelNID(hitId);
j++) {
965 if (GetMCLabelID(hitId,
j) >= (int32_t)GetNMCTracks(GetMCLabelCol(hitId,
j))) {
966 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));
967 throw std::runtime_error(
"qa error");
969 if (GetMCLabelID(hitId,
j) >= 0) {
971 GPUInfo(
"Track %d Cluster %u Label %d: %d (%f)",
i, k,
j, GetMCLabelID(hitId,
j), GetMCLabelWeight(hitId,
j));
977 float maxweight, sumweight;
979 auto maxLabel = acc.computeLabel(&maxweight, &sumweight, &maxcount);
980 mTrackMCLabels[
i] = maxLabel;
981 if (
QA_DEBUG && track.OK() && GetNMCTracks(maxLabel) > (uint32_t)maxLabel.getTrackID()) {
982 const mcInfo_t& mc = GetMCTrack(maxLabel);
983 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,
984 std::sqrt(mc.pX * mc.pX + mc.pY * mc.pY));
989 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
993 for (uint32_t
i = 0;
i < nReconstructedTracks;
i++) {
995 mcLabelI_t
label = mTrackMCLabels[
i];
996 if (mQATasks & taskClusterAttach) {
1001 if (!mTrackMCLabels[
i].
isValid()) {
1002 for (uint32_t k = 0; k < track->NClusters(); k++) {
1006 mClusterParam[mTracking->mIOPtrs.mergedTrackHits[track->FirstClusterRef() + k].num].fakeAttached++;
1010 if (mMCTrackMin == -1 || (
label.getTrackID() >= mMCTrackMin &&
label.getTrackID() < mMCTrackMax)) {
1011 for (uint32_t k = 0; k < track->NClusters(); k++) {
1015 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track->FirstClusterRef() + k].num;
1016 bool correct =
false;
1017 for (int32_t
j = 0;
j < GetMCLabelNID(hitId);
j++) {
1018 if (
label == GetMCLabel(hitId,
j)) {
1024 mClusterParam[hitId].attached++;
1026 mClusterParam[hitId].fakeAttached++;
1032 if (mTrackMCLabels[
i].isFake()) {
1033 (GetMCTrackObj(mFakeTracks,
label))++;
1034 }
else if (tracksExternal || !track->MergedLooper()) {
1035 GetMCTrackObj(mRecTracks,
label)++;
1036 if (mMCTrackMin == -1 || (
label.getTrackID() >= mMCTrackMin &&
label.getTrackID() < mMCTrackMax)) {
1037 int32_t& revLabel = GetMCTrackObj(mTrackMCLabelsReverse,
label);
1038 if (tracksExternal) {
1040 if (revLabel == -1 || fabsf((*tracksExternal)[
i].getZ()) < fabsf((*tracksExternal)[revLabel].getZ())) {
1045 const auto* trks = mTracking->mIOPtrs.mergedTracks;
1047 if (revLabel == -1) {
1049 }
else if (mTracking->GetParam().par.earlyTpcTransform) {
1050 comp = fabsf(trks[
i].GetParam().GetZ() + trks[
i].GetParam().GetTZOffset()) < fabsf(trks[revLabel].GetParam().GetZ() + trks[revLabel].GetParam().GetTZOffset());
1052 float shift1 = mTracking->GetTPCTransformHelper()->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(trks[
i].CSide() *
GPUChainTracking::NSECTORS / 2, trks[
i].GetParam().GetTZOffset());
1053 float shift2 = mTracking->GetTPCTransformHelper()->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(trks[revLabel].CSide() *
GPUChainTracking::NSECTORS / 2, trks[revLabel].GetParam().GetTZOffset());
1054 comp = fabsf(trks[
i].GetParam().GetZ() + shift1) < fabsf(trks[revLabel].GetParam().GetZ() + shift2);
1056 if (revLabel == -1 || !trks[revLabel].OK() || (trks[
i].OK() && comp)) {
1063 if ((mQATasks & taskClusterAttach) && mTracking->mIOPtrs.mergedTrackHitAttachment) {
1065 for (uint32_t
i = 0;
i < GetNMCLabels();
i++) {
1066 if (mClusterParam[
i].attached == 0 && mClusterParam[
i].fakeAttached == 0) {
1067 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[
i];
1070 mcLabelI_t trackL = mTrackMCLabels[track];
1072 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1074 if (trackL == GetMCLabel(
i,
j)) {
1080 mClusterParam[
i].fakeAdjacent++;
1082 mClusterParam[
i].adjacent++;
1089 if (mConfig.matchMCLabels.size()) {
1090 mGoodHits[mNEvents - 1].resize(GetNMCLabels());
1091 std::vector<bool> allowMCLabels(GetNMCTracks(0));
1092 for (uint32_t k = 0; k < GetNMCTracks(0); k++) {
1093 allowMCLabels[k] =
false;
1095 for (uint32_t
i = 0;
i < nReconstructedTracks;
i++) {
1096 if (!mGoodTracks[mNEvents - 1][
i]) {
1099 if (mConfig.matchDisplayMinPt > 0) {
1100 if (!mTrackMCLabels[
i].
isValid()) {
1103 const mcInfo_t& info = GetMCTrack(mTrackMCLabels[
i]);
1104 if (info.pX * info.pX + info.pY * info.pY < mConfig.matchDisplayMinPt * mConfig.matchDisplayMinPt) {
1110 for (uint32_t
j = 0;
j < track.NClusters();
j++) {
1111 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() +
j].num;
1112 if (GetMCLabelNID(hitId)) {
1113 int32_t mcID = GetMCLabelID(hitId, 0);
1115 allowMCLabels[mcID] =
true;
1120 for (uint32_t
i = 0;
i < GetNMCLabels();
i++) {
1121 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1122 int32_t mcID = GetMCLabelID(
i,
j);
1123 if (mcID >= 0 && allowMCLabels[mcID]) {
1124 mGoodHits[mNEvents - 1][
i] =
true;
1129 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1138 for (uint32_t iCol = 0; iCol < GetNMCCollissions(); iCol++) {
1139 for (uint32_t
i = 0;
i < GetNMCTracks(iCol);
i++) {
1140 mMCParam[iCol][
i].nWeightCls = 0.;
1143 for (uint32_t
i = 0;
i < GetNMCLabels();
i++) {
1144 float weightTotal = 0.f;
1145 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1146 if (GetMCLabelID(
i,
j) >= 0) {
1147 weightTotal += GetMCLabelWeight(
i,
j);
1150 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1151 if (GetMCLabelID(
i,
j) >= 0) {
1152 GetMCTrackObj(mMCParam, GetMCLabel(
i,
j)).nWeightCls += GetMCLabelWeight(
i,
j) / weightTotal;
1156 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1157 GPUInfo(
"QA Time: Compute cluster label weights:\t%6.0f us", timer.
GetCurrentElapsedTime(
true) * 1e6);
1161 tbb::parallel_for<uint32_t>(0, GetNMCCollissions(), [&](
auto iCol) {
1162 for (uint32_t
i = 0;
i < GetNMCTracks(iCol);
i++) {
1163 const mcInfo_t& info = GetMCTrack(
i, iCol);
1164 additionalMCParameters& mc2 = mMCParam[iCol][
i];
1165 mc2.pt = std::sqrt(info.pX * info.pX + info.pY * info.pY);
1166 mc2.phi = M_PI + std::atan2(-info.pY, -info.pX);
1167 float p = info.pX * info.pX + info.pY * info.pY + info.pZ * info.pZ;
1169 mc2.theta = mc2.eta = 0.f;
1171 mc2.theta = info.pZ == 0 ? (M_PI / 2) : (
std::acos(info.pZ /
std::sqrt(
p)));
1172 mc2.eta = -std::log(std::tan(0.5 * mc2.theta));
1174 if (mConfig.writeMCLabels) {
1175 std::vector<int32_t>& effBuffer = mcEffBuffer[mNEvents - 1];
1176 effBuffer[
i] = mRecTracks[iCol][
i] * 1000 + mFakeTracks[iCol][
i];
1179 }, tbb::simple_partitioner());
1180 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1181 GPUInfo(
"QA Time: Compute track mc parameters:\t%6.0f us", timer.
GetCurrentElapsedTime(
true) * 1e6);
1185 if (mQATasks & taskTrackingEff) {
1186 for (uint32_t iCol = 0; iCol < GetNMCCollissions(); iCol++) {
1187 for (uint32_t
i = 0;
i < GetNMCTracks(iCol);
i++) {
1188 if ((mMCTrackMin != -1 && (int32_t)
i < mMCTrackMin) || (mMCTrackMax != -1 && (int32_t)
i >= mMCTrackMax)) {
1191 const mcInfo_t& info = GetMCTrack(
i, iCol);
1192 const additionalMCParameters& mc2 = mMCParam[iCol][
i];
1193 if (mc2.nWeightCls == 0.f) {
1196 const float& mcpt = mc2.pt;
1197 const float& mcphi = mc2.phi;
1198 const float& mceta = mc2.eta;
1200 if (info.primDaughters) {
1203 if (mc2.nWeightCls < mConfig.minNClEff) {
1206 int32_t findable = mc2.nWeightCls >= mConfig.minNClFindable;
1210 if (info.charge == 0.f) {
1213 if (mConfig.filterCharge && info.charge * mConfig.filterCharge < 0) {
1216 if (mConfig.filterPID >= 0 && info.pid != mConfig.filterPID) {
1220 if (fabsf(mceta) > ETA_MAX || mcpt < PT_MIN || mcpt > PT_MAX) {
1224 float alpha = std::atan2(info.y, info.x);
1225 alpha /= M_PI / 9.f;
1227 alpha *= M_PI / 9.f;
1228 alpha += M_PI / 18.f;
1230 float c = std::cos(
alpha);
1231 float s = std::sin(
alpha);
1232 float localY = -info.x *
s + info.y *
c;
1234 if (mConfig.dumpToROOT) {
1236 float localX = info.x *
c + info.y *
s;
1237 effdump.Fill(
alpha, localX, localY, info.z, mcphi, mceta, mcpt, mRecTracks[iCol][
i], mFakeTracks[iCol][
i], findable, info.prim, mc2.nWeightCls);
1240 for (int32_t
j = 0;
j < 4;
j++) {
1241 for (int32_t k = 0; k < 2; k++) {
1242 if (k == 0 && findable == 0) {
1246 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;
1251 for (int32_t l = 0; l < 5; l++) {
1252 if (info.prim && mcpt < PT_MIN_PRIM) {
1255 if (l != 3 && fabsf(mceta) > ETA_MAX2) {
1258 if (l < 4 && mcpt < 1.f / mConfig.qpt) {
1262 float pos = l == 0 ? localY : l == 1 ? info.z : l == 2 ? mcphi : l == 3 ? mceta : mcpt;
1264 mEff[
j][k][!info.prim][l]->Fill(
pos,
val);
1270 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1276 if (mQATasks & (taskTrackingRes | taskTrackingResPull)) {
1278 prop.SetMaxSinPhi(.999);
1279 prop.SetMaterialTPC();
1280 prop.SetPolynomialField(&mParam->polynomialField);
1282 for (uint32_t
i = 0;
i < mTrackMCLabels.size();
i++) {
1283 if (mConfig.writeMCLabels) {
1284 std::vector<int32_t>& labelBuffer = mcLabelBuffer[mNEvents - 1];
1285 labelBuffer[
i] = mTrackMCLabels[
i].getTrackID();
1287 if (mTrackMCLabels[
i].isFake()) {
1290 const mcInfo_t& mc1 = GetMCTrack(mTrackMCLabels[
i]);
1291 const additionalMCParameters& mc2 = GetMCTrackObj(mMCParam, mTrackMCLabels[
i]);
1293 if (mc1.primDaughters) {
1296 if (!tracksExternal) {
1297 if (!mTracking->mIOPtrs.mergedTracks[
i].OK()) {
1300 if (mTracking->mIOPtrs.mergedTracks[
i].MergedLooper()) {
1304 if ((mMCTrackMin != -1 && mTrackMCLabels[
i].getTrackID() < mMCTrackMin) || (mMCTrackMax != -1 && mTrackMCLabels[
i].getTrackID() >= mMCTrackMax)) {
1307 if (fabsf(mc2.eta) > ETA_MAX || mc2.pt < PT_MIN || mc2.pt > PT_MAX) {
1310 if (mc1.charge == 0.f) {
1316 if (mc1.t0 == -100.f) {
1319 if (mConfig.filterCharge && mc1.charge * mConfig.filterCharge < 0) {
1322 if (mConfig.filterPID >= 0 && mc1.pid != mConfig.filterPID) {
1325 if (mc2.nWeightCls < mConfig.minNClRes) {
1328 if (mConfig.resPrimaries == 1 && !mc1.prim) {
1330 }
else if (mConfig.resPrimaries == 2 && mc1.prim) {
1333 if (GetMCTrackObj(mTrackMCLabelsReverse, mTrackMCLabels[
i]) != (int32_t)
i) {
1340 if (tracksExternal) {
1342 for (int32_t k = 0; k < 5; k++) {
1343 param.Par()[k] = (*tracksExternal)[
i].getParams()[k];
1345 for (int32_t k = 0; k < 15; k++) {
1346 param.Cov()[k] = (*tracksExternal)[
i].getCov()[k];
1348 param.X() = (*tracksExternal)[
i].getX();
1349 param.TZOffset() = (*tracksExternal)[
i].getTime0();
1350 alpha = (*tracksExternal)[
i].getAlpha();
1351 side = (*tracksExternal)[
i].hasBothSidesClusters() ? 2 : ((*tracksExternal)[
i].hasCSideClusters() ? 1 : 0);
1354 param = mTracking->mIOPtrs.mergedTracks[
i].GetParam();
1355 alpha = mTracking->mIOPtrs.mergedTracks[
i].GetAlpha();
1356 side = mTracking->mIOPtrs.mergedTracks[
i].CCE() ? 2 : (mTracking->mIOPtrs.mergedTracks[
i].CSide() ? 1 : 0);
1360 float c = std::cos(
alpha);
1361 float s = std::sin(
alpha);
1364 mclocal[0] =
x *
c +
y *
s;
1365 mclocal[1] = -
x *
s +
y *
c;
1368 mclocal[2] = px *
c + py *
s;
1369 mclocal[3] = -px *
s + py *
c;
1374 if (mclocal[0] >
param.GetX() + 20) {
1377 if (
param.GetX() > mConfig.maxResX) {
1381 auto getdz = [
this, &
param, &mc1, &
side, tracksExternal]() {
1382 if (tracksExternal) {
1383 return param.GetZ();
1385 if (!mParam->continuousMaxTimeBin) {
1386 return param.GetZ() - mc1.z;
1388#ifdef GPUCA_TPC_GEOMETRY_O2
1389 if (!mParam->par.earlyTpcTransform) {
1391 return param.GetZ() + shift - mc1.z;
1398 bool inFlyDirection = 0;
1399 if (mConfig.strict) {
1401 const float dy =
param.Y() - mclocal[1];
1402 const float dz = getdz();
1403 if (dx * dx + dy * dy + dz * dz > 5.f * 5.f) {
1408 if (prop.PropagateToXAlpha(mclocal[0],
alpha, inFlyDirection)) {
1411 if (fabsf(
param.Y() - mclocal[1]) > (mConfig.strict ? 1.f : 4.f) || fabsf(getdz()) > (mConfig.strict ? 1.f : 4.f)) {
1414 float charge = mc1.charge > 0 ? 1.f : -1.f;
1416 float deltaY =
param.GetY() - mclocal[1];
1417 float deltaZ = getdz();
1418 float deltaPhiNative =
param.GetSinPhi() - mclocal[3] / mc2.pt;
1419 float deltaPhi = std::asin(
param.GetSinPhi()) - std::atan2(mclocal[3], mclocal[2]);
1420 float deltaLambdaNative =
param.GetDzDs() - mc1.pZ / mc2.pt;
1421 float deltaLambda = std::atan(
param.GetDzDs()) - std::atan2(mc1.pZ, mc2.pt);
1423 float deltaPt = (fabsf(1.f /
param.GetQPt()) - mc2.pt) / mc2.pt;
1425 float paramval[5] = {mclocal[1], mc1.z, mc2.phi, mc2.eta, mc2.pt};
1426 float resval[5] = {deltaY, deltaZ, mConfig.nativeFitResolutions ? deltaPhiNative : deltaPhi, mConfig.nativeFitResolutions ? deltaLambdaNative : deltaLambda, mConfig.nativeFitResolutions ? deltaPtNative : deltaPt};
1427 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())};
1429 for (int32_t
j = 0;
j < 5;
j++) {
1430 for (int32_t k = 0; k < 5; k++) {
1431 if (k != 3 && fabsf(mc2.eta) > ETA_MAX2) {
1434 if (k < 4 && mc2.pt < 1.f / mConfig.qpt) {
1437 if (mQATasks & taskTrackingRes) {
1438 mRes2[
j][k]->Fill(resval[
j], paramval[k]);
1440 if (mQATasks & taskTrackingResPull) {
1441 mPull2[
j][k]->Fill(pullval[
j], paramval[k]);
1446 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1451 if (mQATasks & taskClusterAttach) {
1453 for (uint32_t iTrk = 0; iTrk < nReconstructedTracks; iTrk++) {
1458 if (!mTrackMCLabels[iTrk].
isValid()) {
1459 for (uint32_t k = 0; k < track.NClusters(); k++) {
1463 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].num;
1464 float totalWeight = 0.;
1465 for (int32_t
j = 0;
j < GetMCLabelNID(hitId);
j++) {
1467 totalWeight += GetMCLabelWeight(hitId,
j);
1470 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[hitId];
1472 if (totalWeight > 0) {
1473 float weight = 1.f / (totalWeight * (mClusterParam[hitId].attached + mClusterParam[hitId].fakeAttached));
1474 for (int32_t
j = 0;
j < GetMCLabelNID(hitId);
j++) {
1475 mcLabelI_t
label = GetMCLabel(hitId,
j);
1477 float pt = GetMCTrackObj(mMCParam,
label).pt;
1478 if (pt < PT_MIN_CLUST) {
1481 mClusters[CL_fake]->Fill(pt, GetMCLabelWeight(hitId,
j) *
weight);
1482 mClusters[CL_att_adj]->Fill(pt, GetMCLabelWeight(hitId,
j) *
weight);
1483 if (GetMCTrackObj(mRecTracks,
label)) {
1484 mClusters[CL_tracks]->Fill(pt, GetMCLabelWeight(hitId,
j) *
weight);
1486 mClusters[CL_all]->Fill(pt, GetMCLabelWeight(hitId,
j) *
weight);
1487 if (protect || physics) {
1488 mClusters[CL_prot]->Fill(pt, GetMCLabelWeight(hitId,
j) *
weight);
1491 mClusters[CL_physics]->Fill(pt, GetMCLabelWeight(hitId,
j) *
weight);
1496 float weight = 1.f / (mClusterParam[hitId].attached + mClusterParam[hitId].fakeAttached);
1497 mClusters[CL_fake]->Fill(0.f,
weight);
1498 mClusters[CL_att_adj]->Fill(0.f,
weight);
1499 mClusters[CL_all]->Fill(0.f,
weight);
1500 mClusterCounts.nUnaccessible +=
weight;
1501 if (protect || physics) {
1502 mClusters[CL_prot]->Fill(0.f,
weight);
1505 mClusters[CL_physics]->Fill(0.f,
weight);
1511 mcLabelI_t
label = mTrackMCLabels[iTrk];
1512 if (mMCTrackMin != -1 && (
label.getTrackID() < mMCTrackMin ||
label.getTrackID() >= mMCTrackMax)) {
1515 for (uint32_t k = 0; k < track.NClusters(); k++) {
1519 int32_t hitId = mTracking->mIOPtrs.mergedTrackHits[track.FirstClusterRef() + k].num;
1520 float pt = GetMCTrackObj(mMCParam,
label).pt;
1521 if (pt < PT_MIN_CLUST) {
1524 float weight = 1.f / (mClusterParam[hitId].attached + mClusterParam[hitId].fakeAttached);
1525 bool correct =
false;
1526 for (int32_t
j = 0;
j < GetMCLabelNID(hitId);
j++) {
1527 if (
label == GetMCLabel(hitId,
j)) {
1533 mClusters[CL_attached]->Fill(pt,
weight);
1534 mClusters[CL_tracks]->Fill(pt,
weight);
1536 mClusters[CL_fake]->Fill(pt,
weight);
1538 mClusters[CL_att_adj]->Fill(pt,
weight);
1539 mClusters[CL_all]->Fill(pt,
weight);
1540 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[hitId];
1542 if (protect || physics) {
1543 mClusters[CL_prot]->Fill(pt,
weight);
1546 mClusters[CL_physics]->Fill(pt,
weight);
1550 for (uint32_t
i = 0;
i < GetNMCLabels();
i++) {
1551 if ((mMCTrackMin != -1 && GetMCLabelID(
i, 0) < mMCTrackMin) || (mMCTrackMax != -1 && GetMCLabelID(
i, 0) >= mMCTrackMax)) {
1554 if (mClusterParam[
i].attached || mClusterParam[
i].fakeAttached) {
1557 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[
i];
1559 if (mClusterParam[
i].adjacent) {
1562 float totalWeight = 0.;
1563 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1564 mcLabelI_t labelT = GetMCLabel(
i,
j);
1566 totalWeight += GetMCLabelWeight(
i,
j);
1569 float weight = 1.f / totalWeight;
1570 if (totalWeight > 0) {
1571 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1572 mcLabelI_t labelT = GetMCLabel(
i,
j);
1574 float pt = GetMCTrackObj(mMCParam, labelT).pt;
1575 if (pt < PT_MIN_CLUST) {
1578 if (GetMCTrackObj(mRecTracks, labelT)) {
1579 mClusters[CL_tracks]->Fill(pt, GetMCLabelWeight(
i,
j) *
weight);
1581 mClusters[CL_att_adj]->Fill(pt, GetMCLabelWeight(
i,
j) *
weight);
1582 mClusters[CL_fakeAdj]->Fill(pt, GetMCLabelWeight(
i,
j) *
weight);
1583 mClusters[CL_all]->Fill(pt, GetMCLabelWeight(
i,
j) *
weight);
1584 if (protect || physics) {
1585 mClusters[CL_prot]->Fill(pt, GetMCLabelWeight(
i,
j) *
weight);
1588 mClusters[CL_physics]->Fill(pt, GetMCLabelWeight(
i,
j) *
weight);
1593 mClusters[CL_att_adj]->Fill(0.f, 1.f);
1594 mClusters[CL_fakeAdj]->Fill(0.f, 1.f);
1595 mClusters[CL_all]->Fill(0.f, 1.f);
1596 mClusterCounts.nUnaccessible++;
1597 if (protect || physics) {
1598 mClusters[CL_prot]->Fill(0.f, 1.f);
1601 mClusters[CL_physics]->Fill(0.f, 1.f);
1605 float pt = GetMCTrackObj(mMCParam, mTrackMCLabels[
label]).pt;
1606 if (pt < PT_MIN_CLUST) {
1609 mClusters[CL_att_adj]->Fill(pt, 1.f);
1610 mClusters[CL_tracks]->Fill(pt, 1.f);
1611 mClusters[CL_all]->Fill(pt, 1.f);
1612 if (protect || physics) {
1613 mClusters[CL_prot]->Fill(pt, 1.f);
1616 mClusters[CL_physics]->Fill(pt, 1.f);
1620 float totalWeight = 0.;
1621 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1622 mcLabelI_t labelT = GetMCLabel(
i,
j);
1624 totalWeight += GetMCLabelWeight(
i,
j);
1627 if (totalWeight > 0) {
1628 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1629 mcLabelI_t
label = GetMCLabel(
i,
j);
1631 float pt = GetMCTrackObj(mMCParam,
label).pt;
1632 if (pt < PT_MIN_CLUST) {
1635 float weight = GetMCLabelWeight(
i,
j) / totalWeight;
1636 if (mClusterParam[
i].fakeAdjacent) {
1637 mClusters[CL_fakeAdj]->Fill(pt,
weight);
1639 if (mClusterParam[
i].fakeAdjacent) {
1640 mClusters[CL_att_adj]->Fill(pt,
weight);
1642 if (GetMCTrackObj(mRecTracks,
label)) {
1643 mClusters[CL_tracks]->Fill(pt,
weight);
1645 mClusters[CL_all]->Fill(pt,
weight);
1646 if (protect || physics) {
1647 mClusters[CL_prot]->Fill(pt,
weight);
1650 mClusters[CL_physics]->Fill(pt,
weight);
1655 if (mClusterParam[
i].fakeAdjacent) {
1656 mClusters[CL_fakeAdj]->Fill(0.f, 1.f);
1658 if (mClusterParam[
i].fakeAdjacent) {
1659 mClusters[CL_att_adj]->Fill(0.f, 1.f);
1661 mClusters[CL_all]->Fill(0.f, 1.f);
1662 mClusterCounts.nUnaccessible++;
1663 if (protect || physics) {
1664 mClusters[CL_prot]->Fill(0.f, 1.f);
1667 mClusters[CL_physics]->Fill(0.f, 1.f);
1673 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1677 }
else if (!mConfig.inputHistogramsOnly && !mConfig.noMC && (mQATasks & (taskTrackingEff | taskTrackingRes | taskTrackingResPull | taskClusterAttach))) {
1678 GPUWarning(
"No MC information available, only running partial TPC QA!");
1681 if (mQATasks & taskTrackStatistics) {
1683 std::vector<std::array<float, 3>> clusterAttachCounts;
1685 clusterAttachCounts.resize(GetNMCLabels(), {0.f, 0.f});
1687 for (uint32_t
i = 0;
i < nReconstructedTracks;
i++) {
1692 mTracks->Fill(1.f / fabsf(track.GetParam().GetQPt()));
1693 mNCl[0]->Fill(track.NClustersFitted());
1694 uint32_t nClCorrected = 0;
1695 const auto& trackClusters = mTracking->mIOPtrs.mergedTrackHits;
1697 for (uint32_t
j = 0;
j < track.NClusters();
j = jNext) {
1699 for (jNext =
j + 1;
j < track.NClusters(); jNext++) {
1700 if (trackClusters[track.FirstClusterRef() +
j].sector != trackClusters[track.FirstClusterRef() + jNext].sector || trackClusters[track.FirstClusterRef() +
j].row != trackClusters[track.FirstClusterRef() + jNext].row) {
1705 if (trackClusters[track.FirstClusterRef() +
j].leg == trackClusters[track.FirstClusterRef() + track.NClusters() - 1].leg && rowClCount) {
1708 if (mcAvail && rowClCount) {
1709 for (uint32_t k =
j; k < jNext; k++) {
1710 const auto& cl = trackClusters[track.FirstClusterRef() + k];
1714 bool labelOk =
false, labelOkNonFake =
false;
1715 const mcLabelI_t& trkLabel = mTrackMCLabels[
i];
1716 if (trkLabel.isValid() && !trkLabel.isNoise()) {
1717 for (int32_t l = 0; l < GetMCLabelNID(cl.num); l++) {
1718 const mcLabelI_t& clLabel = GetMCLabel(cl.num, l);
1719 if (clLabel.isValid() && !clLabel.isNoise() && CompareIgnoreFake(trkLabel, clLabel)) {
1721 if (!trkLabel.isFake()) {
1722 labelOkNonFake =
true;
1728 clusterAttachCounts[cl.num][0] += 1.0f;
1729 clusterAttachCounts[cl.num][1] += (float)labelOk / rowClCount;
1730 clusterAttachCounts[cl.num][2] += (float)labelOkNonFake / rowClCount;
1734 mNCl[1]->Fill(nClCorrected);
1736 if (mClNative && mTracking && mTracking->GetTPCTransformHelper()) {
1739 for (uint32_t k = 0; k < mClNative->nClusters[
i][
j]; k++) {
1740 const auto& cl = mClNative->clusters[
i][
j][k];
1742 GPUTPCConvertImpl::convert(*mTracking->GetTPCTransformHelper()->getCorrMap(), mTracking->GetParam(),
i,
j, cl.getPad(), cl.getTime(),
x,
y,
z);
1743 mTracking->GetParam().Sector2Global(
i,
x,
y,
z, &
x, &
y, &
z);
1750 double clusterAttachNormalizedCount = 0, clusterAttachNormalizedCountNonFake = 0;
1751 for (uint32_t
i = 0;
i < clusterAttachCounts.size();
i++) {
1752 if (clusterAttachCounts[
i][0]) {
1753 clusterAttachNormalizedCount += clusterAttachCounts[
i][1] / clusterAttachCounts[
i][0];
1754 clusterAttachNormalizedCountNonFake += clusterAttachCounts[
i][2] / clusterAttachCounts[
i][0];
1757 mClusterCounts.nCorrectlyAttachedNormalized = clusterAttachNormalizedCount;
1758 mClusterCounts.nCorrectlyAttachedNormalizedNonFake = clusterAttachNormalizedCountNonFake;
1759 clusterAttachCounts.clear();
1762 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1767 uint32_t nCl = clNative ? clNative->
nClustersTotal : mTracking->GetProcessors()->tpcMerger.NMaxClusters();
1768 mClusterCounts.nTotal += nCl;
1769 if (mQATasks & taskClusterCounts) {
1770 for (uint32_t
i = 0;
i < nCl;
i++) {
1771 int32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[
i];
1775 float totalWeight = 0, weight400 = 0, weight40 = 0;
1776 for (int32_t
j = 0;
j < GetMCLabelNID(
i);
j++) {
1777 const auto&
label = GetMCLabel(
i,
j);
1778 if (GetMCLabelID(
label) >= 0) {
1779 totalWeight += GetMCLabelWeight(
label);
1780 if (GetMCTrackObj(mMCParam,
label).pt >= 0.4) {
1781 weight400 += GetMCLabelWeight(
label);
1783 if (GetMCTrackObj(mMCParam,
label).pt <= 0.04) {
1784 weight40 += GetMCLabelWeight(
label);
1788 if (totalWeight > 0 && 10.f * weight400 >= totalWeight) {
1789 if (!unattached && !protect && !physics) {
1790 mClusterCounts.nFakeRemove400++;
1791 int32_t totalFake = weight400 < 0.9f * totalWeight;
1793 mClusterCounts.nFullFakeRemove400++;
1808 mClusterCounts.nAbove400++;
1810 if (totalWeight > 0 && weight40 >= 0.9 * totalWeight) {
1811 mClusterCounts.nBelow40++;
1812 if (protect || physics) {
1813 mClusterCounts.nFakeProtect40++;
1818 mClusterCounts.nPhysics++;
1820 if (physics || protect) {
1821 mClusterCounts.nProt++;
1824 mClusterCounts.nUnattached++;
1830 if ((mQATasks & taskClusterCounts) && mConfig.clusterRejectionHistograms) {
1831 DoClusterCounts(
nullptr);
1832 mClusterCounts = counts_t();
1835 if (
QA_TIMING || (mTracking && mTracking->GetProcessingSettings().debugLevel >= 3)) {
1839 if (mConfig.dumpToROOT) {
1840 if (!clNative || !mTracking || !mTracking->mIOPtrs.mergedTrackHitAttachment || !mTracking->mIOPtrs.mergedTracks) {
1841 throw std::runtime_error(
"Cannot dump non o2::tpc::clusterNative clusters, need also hit attachmend and GPU tracks");
1846 for (uint32_t k = 0; k < mClNative->nClusters[
i][
j]; k++) {
1847 const auto& cl = mClNative->clusters[
i][
j][k];
1848 uint32_t attach = mTracking->mIOPtrs.mergedTrackHitAttachment[clid];
1849 float x = 0,
y = 0,
z = 0;
1852 const auto& trk = mTracking->mIOPtrs.mergedTracks[track];
1853 mTracking->GetTPCTransformHelper()->Transform(
i,
j, cl.getPad(), cl.getTime(),
x,
y,
z, trk.GetParam().GetTZOffset());
1854 mTracking->GetParam().Sector2Global(
i,
x,
y,
z, &
x, &
y, &
z);
1856 uint32_t extState = mTracking->mIOPtrs.mergedTrackHitStates ? mTracking->mIOPtrs.mergedTrackHitStates[clid] : 0;
1858 if (mConfig.dumpToROOT >= 2) {
1861 memset((
void*)&trk, 0,
sizeof(trk));
1862 memset((
void*)&trkHit, 0,
sizeof(trkHit));
1865 trk = mTracking->mIOPtrs.mergedTracks[track];
1866 for (uint32_t l = 0; l < trk.NClusters(); l++) {
1867 const auto& tmp = mTracking->mIOPtrs.mergedTrackHits[trk.FirstClusterRef() + l];
1868 if (tmp.num == clid) {
1874 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");
1875 cldump.Fill(cl, trk, trkHit, attach, extState,
x,
y,
z,
i,
j, mNEvents - 1);
1877 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");
1878 cldump.Fill(cl, attach, extState,
x,
y,
z,
i,
j, mNEvents - 1);
1886 for (uint32_t
i = 0;
i < mTracking->mIOPtrs.nMergedTracks;
i++) {
1887 if (mTracking->mIOPtrs.mergedTracks[
i].OK()) {
1888 trkdump.Fill(mNEvents - 1, mTracking->mIOPtrs.mergedTracks[
i]);
1892 if (mTracking && mTracking->GetProcessingSettings().createO2Output) {
1894 for (uint32_t
i = 0;
i < mTracking->mIOPtrs.nOutputTracksTPCO2;
i++) {
1895 o2trkdump.Fill(mNEvents - 1, mTracking->mIOPtrs.outputTracksTPCO2[
i]);
1899 mTrackingScratchBuffer.clear();
1900 mTrackingScratchBuffer.shrink_to_fit();
1903void GPUQA::GetName(
char* fname, int32_t k)
1905 const int32_t nNewInput = mConfig.inputHistogramsOnly ? 0 : 1;
1906 if (k || mConfig.inputHistogramsOnly || mConfig.name.size()) {
1907 if (!(mConfig.inputHistogramsOnly || k)) {
1908 snprintf(fname, 1024,
"%s - ", mConfig.name.c_str());
1909 }
else if (mConfig.compareInputNames.size() > (
unsigned)(k - nNewInput)) {
1910 snprintf(fname, 1024,
"%s - ", mConfig.compareInputNames[k - nNewInput].c_str());
1912 strcpy(fname, mConfig.compareInputs[k - nNewInput].c_str());
1913 if (strlen(fname) > 5 && strcmp(fname + strlen(fname) - 5,
".root") == 0) {
1914 fname[strlen(fname) - 5] = 0;
1916 strcat(fname,
" - ");
1924T* GPUQA::GetHist(T*& ee, std::vector<std::unique_ptr<TFile>>& tin, int32_t k, int32_t nNewInput)
1927 if ((mConfig.inputHistogramsOnly || k) && (e =
dynamic_cast<T*
>(tin[k - nNewInput]->Get(e->GetName()))) ==
nullptr) {
1928 GPUWarning(
"Missing histogram in input %s: %s", mConfig.compareInputs[k - nNewInput].c_str(), ee->GetName());
1935void GPUQA::DrawQAHistogramsCleanup()
1937 clearGarbagageCollector();
1940void GPUQA::resetHists()
1942 if (!mQAInitialized) {
1943 throw std::runtime_error(
"QA not initialized");
1945 if (mHaveExternalHists) {
1946 throw std::runtime_error(
"Cannot reset external hists");
1948 for (
auto&
h : *mHist1D) {
1951 for (
auto&
h : *mHist2D) {
1954 for (
auto&
h : *mHist1Dd) {
1957 for (
auto&
h : *mHistGraph) {
1958 h = TGraphAsymmErrors();
1960 mClusterCounts = counts_t();
1965 const auto oldRootIgnoreLevel = gErrorIgnoreLevel;
1966 gErrorIgnoreLevel = kWarning;
1967 if (!mQAInitialized) {
1968 throw std::runtime_error(
"QA not initialized");
1971 if (mTracking && mTracking->GetProcessingSettings().debugLevel >= 2) {
1972 printf(
"Creating QA Histograms\n");
1975 std::vector<Color_t> colorNums(COLORCOUNT);
1977 static int32_t initColorsInitialized = initColors();
1978 (
void)initColorsInitialized;
1980 for (int32_t
i = 0;
i < COLORCOUNT;
i++) {
1981 colorNums[
i] = qcout ? defaultColorNums[
i] : mColors[
i]->GetNumber();
1984 bool mcAvail = mcPresent();
1985 char name[2048], fname[1024];
1987 const int32_t nNewInput = mConfig.inputHistogramsOnly ? 0 : 1;
1988 const int32_t ConfigNumInputs = nNewInput + mConfig.compareInputs.size();
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()));
1994 std::unique_ptr<TFile> tout =
nullptr;
1995 if (mConfig.output.size()) {
1996 tout = std::make_unique<TFile>(mConfig.output.c_str(),
"RECREATE");
1999 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2000 float legendSpacingString = 0.025;
2001 for (int32_t
i = 0;
i < ConfigNumInputs;
i++) {
2003 if (strlen(fname) * 0.006 > legendSpacingString) {
2004 legendSpacingString = strlen(fname) * 0.006;
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.);
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]);
2035 if (mQATasks & taskTrackingRes) {
2036 for (int32_t ii = 0; ii < 7; ii++) {
2037 int32_t
i = ii == 5 ? 4 : ii;
2039 snprintf(fname, 1024,
"res_integral_layout");
2040 snprintf(
name, 2048,
"Integral Resolution");
2042 snprintf(fname, 1024,
"res_vs_%s_layout", VSPARAMETER_NAMES[ii]);
2043 snprintf(
name, 2048,
"Resolution versus %s", VSPARAMETER_NAMES[
i]);
2045 mCRes[ii] = createGarbageCollected<TCanvas>(fname,
name, 0, 0, 700, 700. * 2. / 3.);
2047 gStyle->SetOptFit(1);
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);
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]);
2076 if (mQATasks & taskTrackingResPull) {
2077 for (int32_t ii = 0; ii < 7; ii++) {
2078 int32_t
i = ii == 5 ? 4 : ii;
2081 snprintf(fname, 1024,
"pull_integral_layout");
2082 snprintf(
name, 2048,
"Integral Pull");
2084 snprintf(fname, 1024,
"pull_vs_%s_layout", VSPARAMETER_NAMES[ii]);
2085 snprintf(
name, 2048,
"Pull versus %s", VSPARAMETER_NAMES[
i]);
2087 mCPull[ii] = createGarbageCollected<TCanvas>(fname,
name, 0, 0, 700, 700. * 2. / 3.);
2089 gStyle->SetOptFit(1);
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);
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]);
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.);
2123 mPClust[
i] = createGarbageCollected<TPad>(
"p0",
"", 0.0, 0.0, 1.0, 1.0);
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]);
2132 if (mQATasks & taskTrackStatistics) {
2133 mCTracks = createGarbageCollected<TCanvas>(
"ctracks",
"Track Pt", 0, 0, 700, 700. * 2. / 3.);
2135 mPTracks = createGarbageCollected<TPad>(
"p0",
"", 0.0, 0.0, 1.0, 1.0);
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);
2140 for (int32_t
i = 0;
i < 2;
i++) {
2141 snprintf(
name, 2048,
"cncl%d Pull",
i);
2142 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.);
2144 mPNCl[
i] = createGarbageCollected<TPad>(
"p0",
"", 0.0, 0.0, 1.0, 1.0);
2146 mLNCl[
i] = createGarbageCollected<TLegend>(0.9 - legendSpacingString * 1.45, 0.93 - (0.93 - 0.86) / 2. * (
float)ConfigNumInputs, 0.98, 0.949);
2147 SetLegend(mLNCl[
i]);
2150 mCClXY = createGarbageCollected<TCanvas>(
"clxy",
"Number of clusters per X / Y", 0, 0, 700, 700. * 2. / 3.);
2152 mPClXY = createGarbageCollected<TPad>(
"p0",
"", 0.0, 0.0, 1.0, 1.0);
2157 if (mConfig.enableLocalOutput && !mConfig.inputHistogramsOnly && (mQATasks & taskTrackingEff) && mcPresent()) {
2158 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(),
2159 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(),
2160 (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(),
2161 (int32_t)mRes2[0][3]->GetEntries(), (int32_t)mRes2[0][4]->GetEntries());
2164 int32_t flagShowVsPtLog = (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) ? 1 : 0;
2166 if (mQATasks & taskTrackingEff) {
2168 for (int32_t ii = 0; ii < 5 + flagShowVsPtLog; ii++) {
2169 int32_t
i = ii == 5 ? 4 : ii;
2170 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2171 for (int32_t
j = 0;
j < 4;
j++) {
2172 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2175 mPEff[ii][
j]->SetLogx();
2178 for (int32_t l = 0; l < 3; l++) {
2179 if (k == 0 && mConfig.inputHistogramsOnly == 0 && ii != 5) {
2182 auto oldLevel = gErrorIgnoreLevel;
2183 gErrorIgnoreLevel = kError;
2184 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");
2185 gErrorIgnoreLevel = oldLevel;
2186 mEff[3][
j / 2][
j % 2][
i]->Reset();
2187 mEff[3][
j / 2][
j % 2][
i]->Add(mEff[0][
j / 2][
j % 2][
i]);
2188 mEff[3][
j / 2][
j % 2][
i]->Add(mEff[1][
j / 2][
j % 2][
i]);
2189 mEff[3][
j / 2][
j % 2][
i]->Add(mEff[2][
j / 2][
j % 2][
i]);
2190 mEff[4][
j / 2][
j % 2][
i]->Reset();
2191 mEff[4][
j / 2][
j % 2][
i]->Add(mEff[0][
j / 2][
j % 2][
i]);
2192 mEff[4][
j / 2][
j % 2][
i]->Add(mEff[1][
j / 2][
j % 2][
i]);
2195 auto oldLevel = gErrorIgnoreLevel;
2196 gErrorIgnoreLevel = kError;
2197 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");
2198 gErrorIgnoreLevel = oldLevel;
2202 TGraphAsymmErrors* e = mEffResult[l][
j / 2][
j % 2][
i];
2204 if (!mConfig.inputHistogramsOnly && k == 0) {
2206 mEff[l][
j / 2][
j % 2][
i]->Write();
2209 mEff[3][
j / 2][
j % 2][
i]->Write();
2210 mEff[4][
j / 2][
j % 2][
i]->Write();
2213 }
else if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2216 e->SetTitle(EFFICIENCY_TITLES[
j]);
2217 e->GetYaxis()->SetTitle(
"(Efficiency)");
2218 e->GetXaxis()->SetTitle(XAXIS_TITLES[
i]);
2221 e->SetLineStyle(CONFIG_DASHED_MARKERS ? k + 1 : 1);
2223 if (qcout && !mConfig.shipToQCAsCanvas) {
2226 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2229 e->SetMarkerColor(kBlack);
2230 e->SetLineColor(colorNums[(l == 2 ? (ConfigNumInputs * 2 + k) : (k * 2 + l)) % COLORCOUNT]);
2231 e->GetHistogram()->GetYaxis()->SetRangeUser(-0.02, 1.02);
2232 e->Draw(k || l ?
"same P" :
"AP");
2235 snprintf(
name, 2048,
"%s%s", fname, EFF_NAMES[l]);
2236 mLEff[ii]->AddEntry(e,
name,
"l");
2239 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2243 ChangePadTitleSize(mPEff[ii][
j], 0.056);
2246 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2253 qcout->Add(mCEff[ii]);
2255 if (!mConfig.enableLocalOutput) {
2258 doPerfFigure(0.2, 0.295, 0.025);
2259 mCEff[ii]->Print(Form(
"plots/eff_vs_%s.pdf", VSPARAMETER_NAMES[ii]));
2260 if (mConfig.writeRootFiles) {
2261 mCEff[ii]->Print(Form(
"plots/eff_vs_%s.root", VSPARAMETER_NAMES[ii]));
2266 if (mQATasks & (taskTrackingRes | taskTrackingResPull)) {
2268 TH1D *resIntegral[5] = {}, *pullIntegral[5] = {};
2269 TCanvas* cfit =
nullptr;
2270 std::unique_ptr<TF1> customGaus = std::make_unique<TF1>(
"G",
"[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))");
2271 for (int32_t p = 0;
p < 2;
p++) {
2272 if ((p == 0 && (mQATasks & taskTrackingRes) == 0) || (p == 1 && (mQATasks & taskTrackingResPull) == 0)) {
2275 for (int32_t ii = 0; ii < 5 + flagShowVsPtLog; ii++) {
2276 TCanvas* can =
p ? mCPull[ii] : mCRes[ii];
2277 TLegend*
leg =
p ? mLPull[ii] : mLRes[ii];
2278 int32_t
i = ii == 5 ? 4 : ii;
2279 for (int32_t
j = 0;
j < 5;
j++) {
2280 TH2F*
src =
p ? mPull2[
j][
i] : mRes2[
j][
i];
2281 TH1F**
dst =
p ? mPull[
j][
i] : mRes[
j][
i];
2282 TH1D*& dstIntegral =
p ? pullIntegral[
j] : resIntegral[
j];
2283 TPad* pad =
p ? mPPull[ii][
j] : mPRes[ii][
j];
2285 if (!mConfig.inputHistogramsOnly && ii != 5) {
2286 if (cfit ==
nullptr) {
2287 cfit = createGarbageCollected<TCanvas>();
2291 TAxis* axis =
src->GetYaxis();
2292 int32_t nBins = axis->GetNbins();
2294 for (int32_t bin = 1; bin <= nBins; bin++) {
2295 int32_t bin0 = std::max(bin - integ, 0);
2296 int32_t bin1 = std::min(bin + integ, nBins);
2297 std::unique_ptr<TH1D> proj{
src->ProjectionX(
"proj", bin0, bin1)};
2298 proj->ClearUnderflowAndOverflow();
2299 if (proj->GetEntries()) {
2301 while (proj->GetMaximum() < 50 && rebin <
sizeof(RES_AXIS_BINS) /
sizeof(RES_AXIS_BINS[0])) {
2302 proj->Rebin(RES_AXIS_BINS[rebin - 1] / RES_AXIS_BINS[rebin]);
2306 if (proj->GetEntries() < 20 || proj->GetRMS() < 0.00001) {
2307 dst[0]->SetBinContent(bin, proj->GetRMS());
2308 dst[0]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2309 dst[1]->SetBinContent(bin, proj->GetMean());
2310 dst[1]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2312 proj->GetXaxis()->SetRange(0, 0);
2313 proj->GetXaxis()->SetRangeUser(std::max(proj->GetXaxis()->GetXmin(), proj->GetMean() - 3. * proj->GetRMS()), std::min(proj->GetXaxis()->GetXmax(), proj->GetMean() + 3. * proj->GetRMS()));
2314 bool forceLogLike = proj->GetMaximum() < 20;
2315 for (int32_t k = forceLogLike ? 2 : 0; k < 3; k++) {
2316 proj->Fit(
"gaus", forceLogLike || k == 2 ?
"sQl" : k ?
"sQww" :
"sQ");
2317 TF1* fitFunc = proj->GetFunction(
"gaus");
2319 if (k && !forceLogLike) {
2320 customGaus->SetParameters(fitFunc->GetParameter(0), fitFunc->GetParameter(1), fitFunc->GetParameter(2));
2321 proj->Fit(customGaus.get(),
"sQ");
2322 fitFunc = customGaus.get();
2325 const float sigma = fabs(fitFunc->GetParameter(2));
2326 dst[0]->SetBinContent(bin, sigma);
2327 dst[1]->SetBinContent(bin, fitFunc->GetParameter(1));
2328 dst[0]->SetBinError(bin, fitFunc->GetParError(2));
2329 dst[1]->SetBinError(bin, fitFunc->GetParError(1));
2331 const bool fail1 = sigma <= 0.f;
2332 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());
2333 const bool fail3 =
dst[0]->GetBinContent(bin) > 3.f * proj->GetRMS() ||
dst[0]->GetBinError(bin) > 1 ||
dst[1]->GetBinError(bin) > 1;
2334 const bool fail4 = fitFunc->GetParameter(0) < proj->GetMaximum() / 5.;
2335 const bool fail = fail1 || fail2 || fail3 || fail4;
2340 }
else if (k >= 2) {
2341 dst[0]->SetBinContent(bin, proj->GetRMS());
2342 dst[0]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2343 dst[1]->SetBinContent(bin, proj->GetMean());
2344 dst[1]->SetBinError(bin, std::sqrt(proj->GetRMS()));
2349 dst[0]->SetBinContent(bin, 0.f);
2350 dst[0]->SetBinError(bin, 0.f);
2351 dst[1]->SetBinContent(bin, 0.f);
2352 dst[1]->SetBinError(bin, 0.f);
2356 dstIntegral =
src->ProjectionX(mConfig.nativeFitResolutions ? PARAMETER_NAMES_NATIVE[
j] : PARAMETER_NAMES[
j], 0, nBins + 1);
2358 while (dstIntegral->GetMaximum() < 50 && rebin <
sizeof(RES_AXIS_BINS) /
sizeof(RES_AXIS_BINS[0])) {
2359 dstIntegral->Rebin(RES_AXIS_BINS[rebin - 1] / RES_AXIS_BINS[rebin]);
2365 if (mConfig.inputHistogramsOnly) {
2366 dstIntegral = createGarbageCollected<TH1D>();
2368 snprintf(fname, 1024, p ?
"IntPull%s" :
"IntRes%s", VSPARAMETER_NAMES[
j]);
2369 snprintf(
name, 2048, p ?
"%s Pull" :
"%s Resolution",
p || mConfig.nativeFitResolutions ? PARAMETER_NAMES_NATIVE[
j] : PARAMETER_NAMES[
j]);
2370 dstIntegral->SetName(fname);
2371 dstIntegral->SetTitle(
name);
2373 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2376 int32_t numColor = 0;
2377 float tmpMax = -1000.;
2378 float tmpMin = 1000.;
2380 for (int32_t l = 0; l < 2; l++) {
2381 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2383 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2386 if (nNewInput && k == 0 && ii != 5) {
2388 e->Scale(mConfig.nativeFitResolutions ? SCALE_NATIVE[
j] : SCALE[
j]);
2392 e->GetXaxis()->SetRangeUser(0.2, PT_MAX);
2393 }
else if (LOG_PT_MIN > 0 && ii == 5) {
2394 e->GetXaxis()->SetRangeUser(LOG_PT_MIN, PT_MAX);
2395 }
else if (ii == 5) {
2396 e->GetXaxis()->SetRange(1, 0);
2398 e->SetMinimum(-1111);
2399 e->SetMaximum(-1111);
2401 if (e->GetMaximum() > tmpMax) {
2402 tmpMax = e->GetMaximum();
2404 if (e->GetMinimum() < tmpMin) {
2405 tmpMin = e->GetMinimum();
2411 tmpSpan = tmpMax - tmpMin;
2412 tmpMax += tmpSpan * .02;
2413 tmpMin -= tmpSpan * .02;
2414 if (
j == 2 &&
i < 3) {
2415 tmpMax += tmpSpan * 0.13 * ConfigNumInputs;
2418 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2419 for (int32_t l = 0; l < 2; l++) {
2421 if (!mConfig.inputHistogramsOnly && k == 0) {
2422 snprintf(
name, 2048, p ?
"%s Pull" :
"%s Resolution",
p || mConfig.nativeFitResolutions ? PARAMETER_NAMES_NATIVE[
j] : PARAMETER_NAMES[
j]);
2424 e->SetStats(kFALSE);
2427 mRes2[
j][
i]->SetOption(
"colz");
2428 mRes2[
j][
i]->Write();
2432 }
else if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2435 e->SetMaximum(tmpMax);
2436 e->SetMinimum(tmpMin);
2438 e->SetLineStyle(CONFIG_DASHED_MARKERS ? k + 1 : 1);
2440 e->GetYaxis()->SetTitle(p ? AXIS_TITLES_PULL[
j] : mConfig.nativeFitResolutions ? AXIS_TITLES_NATIVE[
j] : AXIS_TITLES[
j]);
2441 e->GetXaxis()->SetTitle(XAXIS_TITLES[
i]);
2442 if (LOG_PT_MIN > 0 && ii == 5) {
2443 e->GetXaxis()->SetRangeUser(LOG_PT_MIN, PT_MAX);
2447 e->GetYaxis()->SetTitleOffset(1.5);
2449 e->GetYaxis()->SetTitleOffset(1.4);
2451 if (qcout && !mConfig.shipToQCAsCanvas) {
2454 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2458 e->SetMarkerColor(kBlack);
2459 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2460 e->Draw(k || l ?
"same" :
"");
2464 snprintf(
name, 2048,
"%s%s", fname, l ?
"Mean" :
"Pull");
2466 snprintf(
name, 2048,
"%s%s", fname, l ?
"Mean" :
"Resolution");
2472 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2481 ChangePadTitleSize(pad, 0.056);
2484 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2493 if (!mConfig.enableLocalOutput) {
2496 doPerfFigure(0.2, 0.295, 0.025);
2497 can->Print(Form(p ?
"plots/pull_vs_%s.pdf" :
"plots/res_vs_%s.pdf", VSPARAMETER_NAMES[ii]));
2498 if (mConfig.writeRootFiles) {
2499 can->Print(Form(p ?
"plots/pull_vs_%s.root" :
"plots/res_vs_%s.root", VSPARAMETER_NAMES[ii]));
2505 for (int32_t p = 0;
p < 2;
p++) {
2506 if ((p == 0 && (mQATasks & taskTrackingRes) == 0) || (p == 1 && (mQATasks & taskTrackingResPull) == 0)) {
2509 TCanvas* can =
p ? mCPull[6] : mCRes[6];
2510 for (int32_t
i = 0;
i < 5;
i++) {
2511 TPad* pad =
p ? mPPull[6][
i] : mPRes[6][
i];
2512 TH1D* hist =
p ? pullIntegral[
i] : resIntegral[
i];
2513 int32_t numColor = 0;
2514 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2517 if (!mConfig.inputHistogramsOnly && mcAvail) {
2519 if (e && e->GetEntries()) {
2520 e->Fit(
"gaus",
"sQ");
2525 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2527 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2530 e->SetMaximum(-1111);
2531 if (e->GetMaximum() > tmpMax) {
2532 tmpMax = e->GetMaximum();
2536 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2538 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2541 e->SetMaximum(tmpMax * 1.02);
2542 e->SetMinimum(tmpMax * -0.02);
2543 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2546 if (qcout && !mConfig.shipToQCAsCanvas) {
2549 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2553 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2554 e->Draw(k == 0 ?
"" :
"same");
2556 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2564 if (!mConfig.enableLocalOutput) {
2568 can->Print(p ?
"plots/pull_integral.pdf" :
"plots/res_integral.pdf");
2569 if (mConfig.writeRootFiles) {
2570 can->Print(p ?
"plots/pull_integral.root" :
"plots/res_integral.root");
2575 uint64_t attachClusterCounts[N_CLS_HIST];
2576 if (mQATasks & taskClusterAttach) {
2578 if (mConfig.inputHistogramsOnly == 0) {
2579 for (int32_t
i = N_CLS_HIST;
i < N_CLS_TYPE * N_CLS_HIST - 1;
i++) {
2580 mClusters[
i]->Sumw2(
true);
2582 double totalVal = 0;
2583 if (!CLUST_HIST_INT_SUM) {
2584 for (int32_t
j = 0;
j < mClusters[N_CLS_HIST - 1]->GetXaxis()->GetNbins() + 2;
j++) {
2585 totalVal += mClusters[N_CLS_HIST - 1]->GetBinContent(
j);
2588 if (totalVal == 0.) {
2591 for (int32_t
i = 0;
i < N_CLS_HIST;
i++) {
2593 for (int32_t
j = 0;
j < mClusters[
i]->GetXaxis()->GetNbins() + 2;
j++) {
2594 val += mClusters[
i]->GetBinContent(
j);
2595 mClusters[2 * N_CLS_HIST - 1 +
i]->SetBinContent(
j,
val / totalVal);
2597 attachClusterCounts[
i] =
val;
2600 if (!CLUST_HIST_INT_SUM) {
2601 for (int32_t
i = 0;
i < N_CLS_HIST;
i++) {
2602 mClusters[2 * N_CLS_HIST - 1 +
i]->SetMaximum(1.02);
2603 mClusters[2 * N_CLS_HIST - 1 +
i]->SetMinimum(-0.02);
2607 for (int32_t
i = 0;
i < N_CLS_HIST - 1;
i++) {
2608 auto oldLevel = gErrorIgnoreLevel;
2609 gErrorIgnoreLevel = kError;
2610 mClusters[N_CLS_HIST +
i]->Divide(mClusters[
i], mClusters[N_CLS_HIST - 1], 1, 1,
"B");
2611 gErrorIgnoreLevel = oldLevel;
2612 mClusters[N_CLS_HIST +
i]->SetMinimum(-0.02);
2613 mClusters[N_CLS_HIST +
i]->SetMaximum(1.02);
2617 float tmpMax[2] = {0, 0}, tmpMin[2] = {0, 0};
2618 for (int32_t l = 0; l <= CLUST_HIST_INT_SUM; l++) {
2619 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2620 TH1* e = mClusters[l ? (N_CLS_TYPE * N_CLS_HIST - 2) : (N_CLS_HIST - 1)];
2621 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2624 e->SetMinimum(-1111);
2625 e->SetMaximum(-1111);
2627 e->GetXaxis()->SetRange(2, AXIS_BINS[4]);
2629 if (e->GetMaximum() > tmpMax[l]) {
2630 tmpMax[l] = e->GetMaximum();
2632 if (e->GetMinimum() < tmpMin[l]) {
2633 tmpMin[l] = e->GetMinimum();
2636 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2637 for (int32_t
i = 0;
i < N_CLS_HIST;
i++) {
2638 TH1* e = mClusters[l ? (2 * N_CLS_HIST - 1 +
i) :
i];
2639 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2642 e->SetMaximum(tmpMax[l] * 1.02);
2643 e->SetMinimum(tmpMax[l] * -0.02);
2648 for (int32_t
i = 0;
i < N_CLS_TYPE;
i++) {
2649 if (mConfig.enableLocalOutput || mConfig.shipToQCAsCanvas) {
2651 mPClust[
i]->SetLogx();
2653 int32_t
begin =
i == 2 ? (2 * N_CLS_HIST - 1) :
i == 1 ? N_CLS_HIST : 0;
2654 int32_t
end =
i == 2 ? (3 * N_CLS_HIST - 1) :
i == 1 ? (2 * N_CLS_HIST - 1) : N_CLS_HIST;
2655 int32_t numColor = 0;
2656 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2658 TH1* e = mClusters[
j];
2659 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2663 e->SetTitle(CLUSTER_TITLES[
i]);
2664 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)");
2665 e->GetXaxis()->SetTitle(
"#it{p}_{Tmc} (GeV/#it{c})");
2666 e->GetXaxis()->SetTitleOffset(1.1);
2667 e->GetXaxis()->SetLabelOffset(-0.005);
2668 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2671 e->SetStats(kFALSE);
2673 e->SetLineStyle(CONFIG_DASHED_MARKERS ?
j + 1 : 1);
2675 e->GetXaxis()->SetRange(2, AXIS_BINS[4]);
2677 if (qcout && !mConfig.shipToQCAsCanvas) {
2680 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2684 e->SetMarkerColor(kBlack);
2685 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2686 e->Draw(
j ==
end - 1 && k == 0 ?
"" :
"same");
2688 snprintf(
name, 2048,
"%s%s", fname, CLUSTER_NAMES[
j - begin]);
2689 mLClust[
i]->AddEntry(e,
name,
"l");
2692 if (ConfigNumInputs == 1) {
2693 TH1* e =
reinterpret_cast<TH1F*
>(mClusters[
begin + CL_att_adj]->Clone());
2694 e->Add(mClusters[begin + CL_prot], -1);
2695 if (qcout && !mConfig.shipToQCAsCanvas) {
2698 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2702 e->SetLineColor(colorNums[numColor++ % COLORCOUNT]);
2704 mLClust[
i]->AddEntry(e,
"Removed (Strategy A)",
"l");
2706 if (!mConfig.enableLocalOutput && !mConfig.shipToQCAsCanvas) {
2713 qcout->Add(mCClust[
i]);
2715 if (!mConfig.enableLocalOutput) {
2718 doPerfFigure(
i != 2 ? 0.37 : 0.6, 0.295, 0.030);
2720 mCClust[
i]->Print(
i == 2 ?
"plots/clusters_integral.pdf" :
i == 1 ?
"plots/clusters_relative.pdf" :
"plots/clusters.pdf");
2721 if (mConfig.writeRootFiles) {
2722 mCClust[
i]->Print(
i == 2 ?
"plots/clusters_integral.root" :
i == 1 ?
"plots/clusters_relative.root" :
"plots/clusters.root");
2728 if ((mQATasks & taskClusterCounts) && !mHaveExternalHists && !mConfig.clusterRejectionHistograms && !mConfig.inputHistogramsOnly) {
2729 DoClusterCounts(attachClusterCounts);
2731 if ((qcout || tout) && (mQATasks & taskClusterCounts) && mConfig.clusterRejectionHistograms) {
2732 for (uint32_t
i = 0;
i < mHistClusterCount.size();
i++) {
2734 mHistClusterCount[
i]->Write();
2737 qcout->Add(mHistClusterCount[
i]);
2742 if (mQATasks & taskTrackStatistics) {
2745 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2747 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2750 e->SetMaximum(-1111);
2751 if (e->GetMaximum() > tmpMax) {
2752 tmpMax = e->GetMaximum();
2756 mPTracks->SetLogx();
2757 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2759 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2762 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2765 e->SetMaximum(tmpMax * 1.02);
2766 e->SetMinimum(tmpMax * -0.02);
2767 e->SetStats(kFALSE);
2769 e->GetYaxis()->SetTitle(
"a.u.");
2770 e->GetXaxis()->SetTitle(
"#it{p}_{Tmc} (GeV/#it{c})");
2774 e->SetMarkerColor(kBlack);
2775 e->SetLineColor(colorNums[k % COLORCOUNT]);
2776 e->Draw(k == 0 ?
"" :
"same");
2778 snprintf(
name, 2048,
"%sTrack Pt", fname);
2779 mLTracks->AddEntry(e,
name,
"l");
2783 mCTracks->Print(
"plots/tracks.pdf");
2784 if (mConfig.writeRootFiles) {
2785 mCTracks->Print(
"plots/tracks.root");
2788 for (int32_t
i = 0;
i < 2;
i++) {
2790 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2792 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2795 e->SetMaximum(-1111);
2796 if (e->GetMaximum() > tmpMax) {
2797 tmpMax = e->GetMaximum();
2801 for (int32_t k = 0; k < ConfigNumInputs; k++) {
2803 if (GetHist(e, tin, k, nNewInput) ==
nullptr) {
2806 if (tout && !mConfig.inputHistogramsOnly && k == 0) {
2809 e->SetMaximum(tmpMax * 1.02);
2810 e->SetMinimum(tmpMax * -0.02);
2811 e->SetStats(kFALSE);
2813 e->GetYaxis()->SetTitle(
"a.u.");
2814 e->GetXaxis()->SetTitle(
"NClusters");
2818 e->SetMarkerColor(kBlack);
2819 e->SetLineColor(colorNums[k % COLORCOUNT]);
2820 e->Draw(k == 0 ?
"" :
"same");
2822 snprintf(
name, 2048,
"%sNClusters%d", fname,
i);
2823 mLNCl[
i]->AddEntry(e,
name,
"l");
2827 snprintf(
name, 2048,
"plots/nClusters%s.pdf",
i ?
"_corrected" :
"");
2828 mCNCl[
i]->Print(
name);
2829 if (mConfig.writeRootFiles) {
2830 snprintf(
name, 2048,
"plots/nClusters%s.root",
i ?
"_corrected" :
"");
2831 mCNCl[
i]->Print(
name);
2836 mClXY->SetOption(
"colz");
2839 mCClXY->Print(
"plots/clustersXY.pdf");
2840 if (mConfig.writeRootFiles) {
2841 mCClXY->Print(
"plots/clustersXY.root");
2845 if (tout && !mConfig.inputHistogramsOnly && mConfig.writeMCLabels) {
2846 gInterpreter->GenerateDictionary(
"vector<vector<int32_t>>",
"");
2847 tout->WriteObject(&mcEffBuffer,
"mcEffBuffer");
2848 tout->WriteObject(&mcLabelBuffer,
"mcLabelBuffer");
2849 remove(
"AutoDict_vector_vector_int__.cxx");
2850 remove(
"AutoDict_vector_vector_int___cxx_ACLiC_dict_rdict.pcm");
2851 remove(
"AutoDict_vector_vector_int___cxx.d");
2852 remove(
"AutoDict_vector_vector_int___cxx.so");
2858 for (uint32_t
i = 0;
i < mConfig.compareInputs.size();
i++) {
2862 clearGarbagageCollector();
2864 GPUInfo(
"GPU TPC QA histograms have been written to %s files", mConfig.writeRootFiles ?
".pdf and .root" :
".pdf");
2865 gErrorIgnoreLevel = oldRootIgnoreLevel;
2869void GPUQA::PrintClusterCount(int32_t
mode, int32_t&
num,
const char*
name, uint64_t
n, uint64_t normalization)
2873 }
else if (
mode == 1) {
2875 snprintf(name2, 128,
"clusterCount%d_",
num);
2876 char*
ptr = name2 + strlen(name2);
2877 for (uint32_t
i = 0;
i < strlen(
name);
i++) {
2883 createHist(mHistClusterCount[
num], name2,
name, 1000, 0, mConfig.histMaxNClusters, 1000, 0, 100);
2884 }
else if (
mode == 0) {
2885 if (normalization && mConfig.enableLocalOutput) {
2886 printf(
"\t%40s: %'12" PRIu64
" (%6.2f%%)\n",
name,
n, 100.f *
n / normalization);
2888 if (mConfig.clusterRejectionHistograms) {
2889 float ratio = 100.f *
n / std::max<uint64_t>(normalization, 1);
2890 mHistClusterCount[
num]->Fill(normalization, ratio, 1);
2896int32_t GPUQA::DoClusterCounts(uint64_t* attachClusterCounts, int32_t
mode)
2899 if (mcPresent() && (mQATasks & taskClusterAttach) && attachClusterCounts) {
2900 for (int32_t
i = 0;
i < N_CLS_HIST;
i++) {
2901 PrintClusterCount(
mode,
num, CLUSTER_NAMES[
i], attachClusterCounts[
i], mClusterCounts.nTotal);
2903 PrintClusterCount(
mode,
num,
"Unattached", attachClusterCounts[N_CLS_HIST - 1] - attachClusterCounts[CL_att_adj], mClusterCounts.nTotal);
2904 PrintClusterCount(
mode,
num,
"Removed (Strategy A)", attachClusterCounts[CL_att_adj] - attachClusterCounts[CL_prot], mClusterCounts.nTotal);
2905 PrintClusterCount(
mode,
num,
"Unaccessible", mClusterCounts.nUnaccessible, mClusterCounts.nTotal);
2907 PrintClusterCount(
mode,
num,
"All Clusters", mClusterCounts.nTotal, mClusterCounts.nTotal);
2908 PrintClusterCount(
mode,
num,
"Used in Physics", mClusterCounts.nPhysics, mClusterCounts.nTotal);
2909 PrintClusterCount(
mode,
num,
"Protected", mClusterCounts.nProt, mClusterCounts.nTotal);
2910 PrintClusterCount(
mode,
num,
"Unattached", mClusterCounts.nUnattached, mClusterCounts.nTotal);
2911 PrintClusterCount(
mode,
num,
"Removed (Strategy A)", mClusterCounts.nTotal - mClusterCounts.nUnattached - mClusterCounts.nProt, mClusterCounts.nTotal);
2912 PrintClusterCount(
mode,
num,
"Removed (Strategy B)", mClusterCounts.nTotal - mClusterCounts.nProt, mClusterCounts.nTotal);
2915 PrintClusterCount(
mode,
num,
"Merged Loopers (Afterburner)", mClusterCounts.nMergedLooper, mClusterCounts.nTotal);
2916 PrintClusterCount(
mode,
num,
"High Inclination Angle", mClusterCounts.nHighIncl, mClusterCounts.nTotal);
2917 PrintClusterCount(
mode,
num,
"Rejected", mClusterCounts.nRejected, mClusterCounts.nTotal);
2918 PrintClusterCount(
mode,
num,
"Tube (> 200 MeV)", mClusterCounts.nTube, mClusterCounts.nTotal);
2919 PrintClusterCount(
mode,
num,
"Tube (< 200 MeV)", mClusterCounts.nTube200, mClusterCounts.nTotal);
2920 PrintClusterCount(
mode,
num,
"Looping Legs", mClusterCounts.nLoopers, mClusterCounts.nTotal);
2921 PrintClusterCount(
mode,
num,
"Low Pt < 50 MeV", mClusterCounts.nLowPt, mClusterCounts.nTotal);
2922 PrintClusterCount(
mode,
num,
"Low Pt < 200 MeV", mClusterCounts.n200MeV, mClusterCounts.nTotal);
2924 if (mcPresent() && (mQATasks & taskClusterAttach)) {
2925 PrintClusterCount(
mode,
num,
"Tracks > 400 MeV", mClusterCounts.nAbove400, mClusterCounts.nTotal);
2926 PrintClusterCount(
mode,
num,
"Fake Removed (> 400 MeV)", mClusterCounts.nFakeRemove400, mClusterCounts.nAbove400);
2927 PrintClusterCount(
mode,
num,
"Full Fake Removed (> 400 MeV)", mClusterCounts.nFullFakeRemove400, mClusterCounts.nAbove400);
2928 PrintClusterCount(
mode,
num,
"Tracks < 40 MeV", mClusterCounts.nBelow40, mClusterCounts.nTotal);
2929 PrintClusterCount(
mode,
num,
"Fake Protect (< 40 MeV)", mClusterCounts.nFakeProtect40, mClusterCounts.nBelow40);
2931 if (mcPresent() && (mQATasks & taskTrackStatistics)) {
2932 PrintClusterCount(
mode,
num,
"Correctly Attached all-trk normalized", mClusterCounts.nCorrectlyAttachedNormalized, mClusterCounts.nTotal);
2933 PrintClusterCount(
mode,
num,
"Correctly Attached non-fake normalized", mClusterCounts.nCorrectlyAttachedNormalizedNonFake, mClusterCounts.nTotal);
2940 mTrackingScratchBuffer.resize((nBytes +
sizeof(mTrackingScratchBuffer[0]) - 1) /
sizeof(mTrackingScratchBuffer[0]));
2941 return mTrackingScratchBuffer.data();
A const (ready only) version of MCTruthContainer.
Helper class to access correction maps.
#define GPUCA_MIN_TRACK_PTB5_DEFAULT
#define TRACK_EXPECTED_REFERENCE_X_DEFAULT
#define TRACK_EXPECTED_REFERENCE_X
#define CHECK_CLUSTER_STATE_NOCOUNT()
#define CHECK_CLUSTER_STATE()
Definition of the MCTrack class.
Definition of the Names Generator class.
double GetCurrentElapsedTime(bool reset=false)
Class for time synchronization of RawReader instances.
static const HBFUtils & Instance()
static constexpr int32_t NSECTORS
int32_t ReadO2MCData(const char *filename)
bool clusterRemovable(int32_t attach, bool prot) const
void * AllocateScratchBuffer(size_t nBytes)
void SetMCTrackRange(int32_t min, int32_t max)
mcLabelI_t GetMCTrackLabel(uint32_t trackId) const
int32_t DrawQAHistograms()
int32_t InitQA(int32_t tasks=0)
void DumpO2MCData(const char *filename) const
void RunQA(bool matchOnly=false)
static GPUROOTDump< T, Args... > getNew(const char *name1, Names... names)
static DigitizationContext * loadFromFile(std::string_view filename="")
GLfloat GLfloat GLfloat alpha
GLuint GLfloat GLfloat GLfloat GLfloat y1
GLuint const GLchar * name
GLuint GLuint GLfloat weight
GLuint GLsizei const GLchar * label
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLdouble GLdouble GLdouble z
constexpr int LHCBCPERTIMEBIN
Enum< T >::Iterator begin(Enum< T >)
struct o2::upgrades_utils::@454 tracks
structure to keep trigger-related info
Defining DataPointCompositeObject explicitly as copiable.
bool isValid(std::string alias)
int64_t differenceInBC(const InteractionRecord &other) const
std::tuple< std::vector< std::unique_ptr< TCanvas > >, std::vector< std::unique_ptr< TLegend > >, std::vector< std::unique_ptr< TPad > >, std::vector< std::unique_ptr< TLatex > >, std::vector< std::unique_ptr< TH1D > > > v
IR getFirstIRofTF(const IR &rec) const
get 1st IR of TF corresponding to the 1st sampled orbit (in MC)
unsigned int nClustersTotal
o2::InteractionRecord ir(0, 0)
o2::InteractionRecord ir0(3, 5)