CbmRoot
Loading...
Searching...
No Matches
qa/CbmRichRecoQa.cxx
Go to the documentation of this file.
1/* Copyright (C) 2016-2024 GSI/JINR-LIT, Darmstadt/Dubna
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Semen Lebedev, Andrey Lebedev [committer], Martin Beyer */
4
5#include "CbmRichRecoQa.h"
6
7#include "CbmDigiManager.h"
8#include "CbmDrawHist.h"
9#include "CbmGlobalTrack.h"
10#include "CbmHistManager.h"
11#include "CbmMCDataArray.h"
12#include "CbmMCEventList.h"
13#include "CbmMCTrack.h"
14#include "CbmMatchRecoToMC.h"
15#include "CbmRichDraw.h"
16#include "CbmRichGeoManager.h"
17#include "CbmRichHit.h"
18#include "CbmRichPoint.h"
19#include "CbmRichRing.h"
20#include "CbmRichUtil.h"
21#include "CbmTrackMatchNew.h"
22#include "CbmUtils.h"
23#include "FairMCPoint.h"
24#include "FairTrackParam.h"
25#include "TCanvas.h"
26#include "TClonesArray.h"
27#include "TEllipse.h"
28#include "TF1.h"
29#include "TH1.h"
30#include "TH1D.h"
31#include "TMCProcess.h"
32#include "TMarker.h"
33#include "TStyle.h"
34#include "elid/CbmLitGlobalElectronId.h"
35
36#include <TFile.h>
37
38#include <iostream>
39#include <sstream>
40#include <string>
41
42using namespace std;
43using namespace Cbm;
44
45CbmRichRecoQa::CbmRichRecoQa() : FairTask("CbmRichRecoQa") {}
46
47
48InitStatus CbmRichRecoQa::Init()
49{
50 fMcTracks = InitOrFatalMc("MCTrack", "CbmRichRecoQa::Init");
51 fRichPoints = InitOrFatalMc("RichPoint", "CbmRichRecoQa::Init");
52 fRichHits = GetOrFatal<TClonesArray>("RichHit", "CbmRichRecoQa::Init");
53 fRichRings = GetOrFatal<TClonesArray>("RichRing", "CbmRichRecoQa::Init");
54 fRichRingMatches = GetOrFatal<TClonesArray>("RichRingMatch", "CbmRichRecoQa::Init");
55 fRichProjections = GetOrFatal<TClonesArray>("RichProjection", "CbmRichRecoQa::Init");
56 fGlobalTracks = GetOrFatal<TClonesArray>("GlobalTrack", "CbmRichRecoQa::Init");
57 fStsTracks = GetOrFatal<TClonesArray>("StsTrack", "CbmRichRecoQa::Init");
58 fStsTrackMatches = GetOrFatal<TClonesArray>("StsTrackMatch", "CbmRichRecoQa::Init");
59 fEventList = GetOrFatal<CbmMCEventList>("MCEventList.", "CbmRichUrqmdTest::Init");
60
61 fDigiMan = CbmDigiManager::Instance();
62 fDigiMan->Init();
63
64 InitHistograms();
65
66 // CbmLitGlobalElectronId::GetInstance();
67
68 return kSUCCESS;
69}
70
71void CbmRichRecoQa::InitHistograms()
72{
73 fHM = new CbmHistManager();
74
75 // 2D histogramm XY
76 double xMin = -120.;
77 double xMax = 120.;
78 int nBinsX = 30;
79 double yMin = -208;
80 double yMax = 208.;
81 int nBinsY = 52;
82
83 // 1D histogram X or Y
84 int nBinsX1 = 60;
85 int xMin1 = -120;
86 int xMax1 = 120;
87 int nBinsY1 = 25;
88 int yMin1 = 100;
89 int yMax1 = 200;
90
91 vector<string> matchTypes{"Primel", "Pi", "PrimelPlus", "PrimelMinus"};
92 for (const string& t : matchTypes) {
93 fHM->Create2<TH2D>("fhRTDistVsMomTruematch" + t,
94 "fhRTDistVsMomTruematch" + t + ";P [GeV/c];Ring-track distance [cm];Yield", 20, 0., 10., 100, 0.,
95 5.);
96 fHM->Create2<TH2D>("fhRTDistVsMomWrongmatch" + t,
97 "fhRTDistVsMomWrongmatch" + t + ";P [GeV/c];Ring-track distance [cm];Yield", 20, 0., 10., 100,
98 0., 5.);
99
100 fHM->Create2<TH2D>("fhRTDistVsNofHitsTruematch" + t,
101 "fhRTDistVsNofHitsTruematch" + t + ";# hits/ring;Ring-track distance [cm];Yield", 40, -.5, 39.5,
102 100, 0., 5.);
103 fHM->Create3<TH3D>("fhRTDistVsXYTruematch" + t,
104 "fhRTDistVsXYTruematch" + t + ";X [cm];Y [cm];Ring-track distance [cm]", nBinsX, xMin, xMax,
105 nBinsY, yMin, yMax, 100, 0., 5.);
106 fHM->Create2<TH2D>("fhRTDistVsXTruematch" + t, "fhRTDistVsXTruematch" + t + ";X [cm];Ring-track distance [cm]",
107 nBinsX1, xMin1, xMax1, 100, 0., 5.);
108 fHM->Create2<TH2D>("fhRTDistVsYTruematch" + t, "fhRTDistVsYTruematch" + t + ";Abs(Y) [cm];Ring-track distance [cm]",
109 nBinsY1, yMin1, yMax1, 100, 0., 5.);
110 }
111
112 // after electron identification
113 fHM->Create2<TH2D>("fhRTDistVsMomTruematchElId",
114 "fhRTDistVsMomTruematchElId;P [GeV/c];Ring-track distance [cm];Yield", 20, 0., 10., 100, 0., 5.);
115
116 fHM->Create1<TH1D>("fhMismatchSrc", "fhMismatchSrc;Global track category;% from MC", 13, -0.5, 12.5);
117
118 vector<string> mismatchTypes{
119 "Mc", "Sts", "StsAccRich", "StsRich", "StsRichTrue", "StsNoRich", "StsNoRichRF",
120 "StsNoRichRM", "StsNoRichNoRF", "StsNoRichNoProj", "StsRichWrong", "StsRichWrongRF", "StsRichWrongRM"};
121 for (const string& t : mismatchTypes) {
122 fHM->Create1<TH1D>("fhMismatchSrcMom" + t, "fhMismatchSrcMom" + t + ";P [GeV/c];Yield", 40, 0., 10.);
123 }
124}
125
126void CbmRichRecoQa::Exec(Option_t* /*option*/)
127{
128 fEventNum++;
129 LOG(info) << "CbmRichRecoQa, event No. " << fEventNum;
130 FillRichRingNofHits();
131 FillRingTrackDistance();
132 RingTrackMismatchSource();
133}
134
135void CbmRichRecoQa::FillRichRingNofHits()
136{
137 fNofHitsInRingMap.clear();
138 int nofRichHits = fRichHits->GetEntriesFast();
139 for (int iHit = 0; iHit < nofRichHits; iHit++) {
140 const CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(iHit));
141 if (nullptr == hit) continue;
142 vector<CbmLink> motherIds = CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(fDigiMan, hit, fRichPoints, fMcTracks);
143
144 for (const auto& motherId : motherIds) {
145 fNofHitsInRingMap[motherId]++;
146 }
147 }
148}
149
150void CbmRichRecoQa::RingTrackMismatchSource()
151{
152
153 int nofEvents = fEventList->GetNofEvents();
154 for (int iE = 0; iE < nofEvents; iE++) {
155 int fileId = fEventList->GetFileIdByIndex(iE);
156 int eventId = fEventList->GetEventIdByIndex(iE);
157
158 int nMcTracks = fMcTracks->Size(fileId, eventId);
159 for (int i = 0; i < nMcTracks; i++) {
160 //At least one hit in RICH
161 CbmLink val(1., i, eventId, fileId);
162 if (fNofHitsInRingMap[val] < 1) continue;
163 const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(val));
164 if (IsMcPrimaryElectron(mcTrack)) {
165 fHM->H1("fhMismatchSrc")->Fill(0); // MC
166 fHM->H1("fhMismatchSrcMomMc")->Fill(mcTrack->GetP());
167 }
168 }
169 }
170
171 int nofGlobalTracks = fGlobalTracks->GetEntriesFast();
172 for (int iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
173 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
174
175 int stsId = globalTrack->GetStsTrackIndex();
176 if (stsId < 0) continue;
177 const CbmTrackMatchNew* stsTrackMatch = static_cast<const CbmTrackMatchNew*>(fStsTrackMatches->At(stsId));
178 if (stsTrackMatch == nullptr || stsTrackMatch->GetNofLinks() < 1) continue;
179 auto stsMatchedLink = stsTrackMatch->GetMatchedLink();
180 const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(stsMatchedLink));
181 if (!IsMcPrimaryElectron(mcTrack)) continue;
182 double mom = mcTrack->GetP();
183 //STS
184 fHM->H1("fhMismatchSrc")->Fill(1);
185 fHM->H1("fhMismatchSrcMomSts")->Fill(mom);
186
187 if (fNofHitsInRingMap[stsMatchedLink] >= 7) {
188 //Sts-AccRich
189 fHM->H1("fhMismatchSrc")->Fill(2);
190 fHM->H1("fhMismatchSrcMomStsAccRich")->Fill(mcTrack->GetP());
191 }
192
193 int richId = globalTrack->GetRichRingIndex();
194 // No RICH segment
195 if (richId < 0) {
196 fHM->H1("fhMismatchSrc")->Fill(5); //STS-noRICH
197 fHM->H1("fhMismatchSrcMomStsNoRich")->Fill(mom);
198 bool ringFound = WasRingFound(stsMatchedLink);
199 bool ringMatched = WasRingMatched(stsMatchedLink);
200 bool hasProj = HasRichProjection(stsId);
201 if (ringFound) {
202 fHM->H1("fhMismatchSrc")->Fill(6); //STS-NoRich RF
203 fHM->H1("fhMismatchSrcMomStsNoRichRF")->Fill(mom);
204 }
205 else {
206 fHM->H1("fhMismatchSrc")->Fill(8); //STS-NoRich NoRF
207 fHM->H1("fhMismatchSrcMomStsNoRichNoRF")->Fill(mom);
208 }
209
210 if (ringMatched) {
211 fHM->H1("fhMismatchSrc")->Fill(7); //STS-NoRich RM
212 fHM->H1("fhMismatchSrcMomStsNoRichRM")->Fill(mom);
213 }
214
215 if (!hasProj) {
216 fHM->H1("fhMismatchSrc")->Fill(9); //STS-NoRich NoProj
217 fHM->H1("fhMismatchSrcMomStsNoRichNoProj")->Fill(mom);
218 }
219 }
220 else {
221
222 //STS-RICH
223 fHM->H1("fhMismatchSrc")->Fill(3);
224 fHM->H1("fhMismatchSrcMomStsRich")->Fill(mom);
225 const CbmTrackMatchNew* richRingMatch = static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
226 if (richRingMatch == nullptr || richRingMatch->GetNofLinks() < 1) continue;
227 auto richMcTrackLink = richRingMatch->GetMatchedLink();
228 const CbmRichRing* ring = static_cast<const CbmRichRing*>(fRichRings->At(richId));
229 if (nullptr == ring) continue;
230
231 if (stsMatchedLink == richMcTrackLink) {
232 fHM->H1("fhMismatchSrc")->Fill(4); //STS-RICH true
233 fHM->H1("fhMismatchSrcMomStsRichTrue")->Fill(mom);
234 }
235 else {
236 fHM->H1("fhMismatchSrc")->Fill(10); //STS-RICH wrong
237 fHM->H1("fhMismatchSrcMomStsRichWrong")->Fill(mom);
238 if (WasRingFound(stsMatchedLink)) {
239 fHM->H1("fhMismatchSrc")->Fill(11); //STS-RICH wrong RF
240 fHM->H1("fhMismatchSrcMomStsRichWrongRF")->Fill(mom);
241 }
242
243 if (WasRingFound(stsMatchedLink)) {
244 fHM->H1("fhMismatchSrc")->Fill(12); //STS-RICH wrong RM
245 fHM->H1("fhMismatchSrcMomStsRichWrongRM")->Fill(mom);
246 }
247 }
248 }
249 }
250}
251
252bool CbmRichRecoQa::WasRingFound(const CbmLink& mcTrackLink)
253{
254 int nofRings = fRichRings->GetEntriesFast();
255 for (int iR = 0; iR < nofRings; iR++) {
256 const CbmRichRing* ring = static_cast<const CbmRichRing*>(fRichRings->At(iR));
257 if (ring == nullptr) continue;
258 const CbmTrackMatchNew* richRingMatch = static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(iR));
259 if (richRingMatch == nullptr || richRingMatch->GetNofLinks() < 1) continue;
260 auto richMcTrackLink = richRingMatch->GetMatchedLink();
261 if (richMcTrackLink == mcTrackLink) return true;
262 }
263 return false;
264}
265
266bool CbmRichRecoQa::WasRingMatched(const CbmLink& mcTrackLink)
267{
268 int nofGlobalTracks = fGlobalTracks->GetEntriesFast();
269 for (int iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
270 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
271
272 int richId = globalTrack->GetRichRingIndex();
273 if (richId < 0) continue;
274 const CbmTrackMatchNew* richRingMatch = static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
275 if (richRingMatch == nullptr || richRingMatch->GetNofLinks() < 1) continue;
276 auto richMcTrackLink = richRingMatch->GetMatchedLink();
277 if (richMcTrackLink == mcTrackLink) {
278 return true;
279 }
280 }
281
282 return false;
283}
284
285
286bool CbmRichRecoQa::HasRichProjection(int stsTrackId)
287{
288 if (stsTrackId < 0) return false;
289 const FairTrackParam* pTrack = static_cast<FairTrackParam*>(fRichProjections->At(stsTrackId));
290 if (pTrack == nullptr) return false;
291
292 if (pTrack->GetX() == 0. && pTrack->GetY() == 0.) return false;
293
294 return true;
295}
296
297void CbmRichRecoQa::FillRingTrackDistance()
298{
299 int nofGlobalTracks = fGlobalTracks->GetEntriesFast();
300 for (int iTrack = 0; iTrack < nofGlobalTracks; iTrack++) {
301 const CbmGlobalTrack* globalTrack = static_cast<const CbmGlobalTrack*>(fGlobalTracks->At(iTrack));
302
303 int stsId = globalTrack->GetStsTrackIndex();
304 int richId = globalTrack->GetRichRingIndex();
305 if (stsId < 0 || richId < 0) continue;
306
307 const CbmTrackMatchNew* stsTrackMatch = static_cast<const CbmTrackMatchNew*>(fStsTrackMatches->At(stsId));
308 if (stsTrackMatch == nullptr || stsTrackMatch->GetNofLinks() < 1) continue;
309 auto stsMcMatchedLink = stsTrackMatch->GetMatchedLink();
310
311 const CbmTrackMatchNew* richRingMatch = static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(richId));
312 if (richRingMatch == nullptr || richRingMatch->GetNofLinks() < 1) continue;
313 auto richMcTrackLink = richRingMatch->GetMatchedLink();
314 const CbmRichRing* ring = static_cast<const CbmRichRing*>(fRichRings->At(richId));
315 if (nullptr == ring) continue;
316
317 double rtDistance = CbmRichUtil::GetRingTrackDistance(iTrack);
318 double xc = ring->GetCenterX();
319 double yc = ring->GetCenterY();
320 int nofHits = ring->GetNofHits();
321
322 const CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->Get(stsMcMatchedLink));
323 if (mcTrack == nullptr) continue;
324 double mom = mcTrack->GetP();
325 int charge = mcTrack->GetCharge();
326
327 bool isEl = IsMcPrimaryElectron(mcTrack);
328 bool isPi = IsMcPion(mcTrack);
329
330 if (!isEl && !isPi) continue;
331
332 vector<string> matchTypes{"Primel", "Pi", "PrimelPlus", "PrimelMinus"};
333 for (const string& t : matchTypes) {
334 //cout << "t:" << t << " iT:" << iTrack << " dist:" << rtDistance << " if:" << (!(t == "Primel" && isEl)?"true":"false")<< endl;
335 bool isGood = false;
336 if (t == "Primel" && isEl) isGood = true;
337 if (t == "Pi" && isPi) isGood = true;
338 if (t == "PrimelPlus" && isEl && charge > 0) isGood = true;
339 if (t == "PrimelMinus" && isEl && charge < 0) isGood = true;
340 if (!isGood) continue;
341 if (stsMcMatchedLink == richMcTrackLink) {
342 fHM->H2("fhRTDistVsMomTruematch" + t)->Fill(mom, rtDistance);
343 fHM->H3("fhRTDistVsXYTruematch" + t)->Fill(xc, yc, rtDistance);
344 fHM->H2("fhRTDistVsXTruematch" + t)->Fill(xc, rtDistance);
345 fHM->H2("fhRTDistVsYTruematch" + t)->Fill(abs(yc), rtDistance);
346 fHM->H2("fhRTDistVsNofHitsTruematch" + t)->Fill(nofHits, rtDistance);
347
348 if (t == "Primel" && isEl) {
349 //if (CbmLitGlobalElectronId::GetInstance().IsRichElectron(iTrack, mom)){
350 fHM->H2("fhRTDistVsMomTruematchElId")->Fill(mom, rtDistance);
351 //}
352 }
353 }
354 else {
355 fHM->H2("fhRTDistVsMomWrongmatch" + t)->Fill(mom, rtDistance);
356 }
357 }
358 }
359}
360
361void CbmRichRecoQa::DrawMismatchSrc()
362{
363 gStyle->SetPaintTextFormat("4.1f");
364 fHM->CreateCanvas("richqa_mismatch_src", "richqa_mismatch_src", 1000, 800);
365 double nofMcEl = fHM->H1("fhMismatchSrc")->GetBinContent(1);
366 fHM->Scale("fhMismatchSrc", 100. / nofMcEl);
367 DrawH1(fHM->H1("fhMismatchSrc"), kLinear, kLog, "hist text");
368 fHM->H1("fhMismatchSrc")->SetMarkerSize(1.9);
369
370 vector<string> labels{"MC",
371 "STS",
372 "STS-Acc RICH",
373 "STS-RICH",
374 "STS-RICH true",
375 "STS-NoRICH",
376 "STS-NoRICH RF",
377 "STS-NoRICH RM",
378 "STS-NoRICH NoRF",
379 "STS-NoRICH NoPrj",
380 "STS-RICH wrong",
381 "STS-RICH wrong RF",
382 "STS-RICH wrong RM"};
383
384 for (size_t i = 0; i < labels.size(); i++) {
385 fHM->H1("fhMismatchSrc")->GetXaxis()->SetBinLabel(i + 1, labels[i].c_str());
386 }
387 fHM->H1("fhMismatchSrc")->GetXaxis()->SetLabelSize(0.03);
388}
389
390void CbmRichRecoQa::DrawHist()
391{
392 SetDefaultDrawStyle();
393
394 {
395 fHM->CreateCanvas("richqa_RTDist_truematch_epem", "richqa_RTDist_truematch_epem", 1000, 1000);
396 TH1D* py_minus = (TH1D*) (fHM->H2("fhRTDistVsMomTruematchPrimelMinus")
397 ->ProjectionY("fhRTDistVsMomTruematchPrimelMinus_py")
398 ->Clone());
399 py_minus->Scale(1. / py_minus->Integral());
400 TH1D* py_plus = (TH1D*) (fHM->H2("fhRTDistVsMomTruematchPrimelPlus")
401 ->ProjectionY("fhRTDistVsMomTruematchPrimelPlus_py")
402 ->Clone());
403 py_plus->Scale(1. / py_plus->Integral());
404 DrawH1({py_minus, py_plus},
405 {"e_{prim}^{-} (" + GetMeanRmsOverflowString(py_minus, false) + ")",
406 "e_{prim}^{+} (" + GetMeanRmsOverflowString(py_plus, false) + ")"},
407 kLinear, kLog, true, 0.5, 0.85, 0.99, 0.99, "hist");
408 }
409
410 {
411 fHM->CreateCanvas("richqa_RTDist_truematch_elpi", "richqa_RTDist_truematch_elpi", 800, 800);
412 TH1D* py_el =
413 (TH1D*) (fHM->H2("fhRTDistVsMomTruematchPrimel")->ProjectionY("fhRTDistVsMomTruematchPrimel_py")->Clone());
414 py_el->Scale(1. / py_el->Integral());
415 TH1D* py_pi = (TH1D*) (fHM->H2("fhRTDistVsMomTruematchPi")->ProjectionY("fhRTDistVsMomTruematchPi_py")->Clone());
416 py_pi->Scale(1. / py_pi->Integral());
417 DrawH1({py_el, py_pi},
418 {"e_{prim}^{#pm} (" + GetMeanRmsOverflowString(py_el, false) + ")",
419 "#pi^{#pm} (" + GetMeanRmsOverflowString(py_pi, false) + ")"},
420 kLinear, kLog, true, 0.5, 0.85, 0.99, 0.99, "hist");
421 }
422
423 {
424 fHM->CreateCanvas("richqa_mismatch_src_mom", "richqa_mismatch_src_mom", 1000, 800);
425 vector<string> labels{"MC",
426 "STS",
427 "STS-Acc RICH",
428 "STS-RICH",
429 "STS-RICH true",
430 "STS-NoRICH",
431 "STS-NoRICH RF",
432 "STS-NoRICH RM",
433 "STS-NoRICH NoRF",
434 "STS-NoRICH NoPrj",
435 "STS-RICH wrong",
436 "STS-RICH wrong RF",
437 "STS-RICH wrong RM"};
438 vector<TH1*> hists = fHM->H1Vector(
439 {"fhMismatchSrcMomMc", "fhMismatchSrcMomSts", "fhMismatchSrcMomStsAccRich", "fhMismatchSrcMomStsRich",
440 "fhMismatchSrcMomStsRichTrue", "fhMismatchSrcMomStsNoRich", "fhMismatchSrcMomStsNoRichRF",
441 "fhMismatchSrcMomStsNoRichRM", "fhMismatchSrcMomStsNoRichNoRF", "fhMismatchSrcMomStsNoRichNoProj",
442 "fhMismatchSrcMomStsRichWrong", "fhMismatchSrcMomStsRichWrongRF", "fhMismatchSrcMomStsRichWrongRM"});
443
444 DrawH1(hists, labels, kLinear, kLog, true, 0.8, 0.35, 0.99, 0.99);
445 fHM->H1("fhMismatchSrcMomMc")->SetMinimum(0.9);
446 }
447
448 DrawRingTrackDist("Primel");
449 DrawRingTrackDist("PrimelPlus");
450 DrawRingTrackDist("PrimelMinus");
451 DrawRingTrackDist("Pi");
452
453 // before and after electron identification
454 {
455 TCanvas* c = fHM->CreateCanvas("richqa_RTDist_truematch_elid", "richqa_RTDist_truematch_elid", 1800, 600);
456 c->Divide(3, 1);
457 c->cd(1);
458 DrawH2WithProfile(fHM->H2("fhRTDistVsMomTruematchPrimel"), false, true);
459 c->cd(2);
460 DrawH2WithProfile(fHM->H2("fhRTDistVsMomTruematchElId"), false, true);
461 c->cd(3);
462 TH1D* py = (TH1D*) (fHM->H2("fhRTDistVsMomTruematchPrimel")
463 ->ProjectionY(string("fhRTDistVsMomTruematchPrimel_py2").c_str())
464 ->Clone());
465 TH1D* pyElId = (TH1D*) (fHM->H2("fhRTDistVsMomTruematchElId")
466 ->ProjectionY(string("fhRTDistVsMomTruematchElId_py").c_str())
467 ->Clone());
468 DrawH1(
469 {py, pyElId},
470 {"before ElId (" + GetMeanRmsOverflowString(py) + ")", "after ElId (" + GetMeanRmsOverflowString(pyElId) + ")"},
471 kLinear, kLog, true, 0.3, 0.75, 0.99, 0.99);
472 fHM->H2("fhRTDistVsMomTruematchPrimel")->GetYaxis()->SetRangeUser(0., 2.);
473 fHM->H2("fhRTDistVsMomTruematchElId")->GetYaxis()->SetRangeUser(0., 2.);
474 }
475}
476
477string CbmRichRecoQa::GetMeanRmsOverflowString(TH1* h, bool withOverflow)
478{
479 if (withOverflow) {
480 double overflow = h->GetBinContent(h->GetNbinsX() + 1);
481 return Cbm::NumberToString<double>(h->GetMean(), 2) + " / " + Cbm::NumberToString<double>(h->GetRMS(), 2) + " / "
482 + Cbm::NumberToString<double>(100. * overflow / h->Integral(0, h->GetNbinsX() + 1), 2) + "%";
483 }
484 else {
485 return Cbm::NumberToString<double>(h->GetMean(), 2) + " / " + Cbm::NumberToString<double>(h->GetRMS(), 2);
486 }
487}
488
489
490void CbmRichRecoQa::DrawRingTrackDist(const string& opt)
491{
492 {
493 TCanvas* c = fHM->CreateCanvas("richqa_RTDist_truematch_" + opt, "richqa_RTDist_truematch_" + opt, 1800, 600);
494 c->Divide(3, 1);
495 c->cd(1);
496 DrawH2WithProfile(fHM->H2("fhRTDistVsMomTruematch" + opt), false, true);
497 c->cd(2);
498 DrawH2WithProfile(fHM->H2("fhRTDistVsNofHitsTruematch" + opt), false, true);
499 c->cd(3);
500 TH1D* py = (TH1D*) (fHM->H2("fhRTDistVsMomTruematch" + opt)
501 ->ProjectionY(("fhRTDistVsMomTruematch_py" + opt).c_str())
502 ->Clone());
503 DrawH1(py);
504 DrawTextOnPad(this->GetMeanRmsOverflowString(py), 0.2, 0.9, 0.8, 0.99);
505 gPad->SetLogy(true);
506
507 fHM->H2("fhRTDistVsMomTruematch" + opt)->GetYaxis()->SetRangeUser(0., 2.);
508 fHM->H2("fhRTDistVsNofHitsTruematch" + opt)->GetYaxis()->SetRangeUser(0., 2.);
509 }
510
511 {
512 fHM->CreateCanvas("richqa_RTDist_wrongmatch" + opt, "richqa_RTDist_wrongmatch" + opt, 600, 600);
513 TH1D* py = (TH1D*) (fHM->H2("fhRTDistVsMomWrongmatch" + opt)
514 ->ProjectionY(("fhRTDistVsMomWrongmatch_py" + opt).c_str())
515 ->Clone());
516 DrawH1(py);
517 DrawTextOnPad(this->GetMeanRmsOverflowString(py), 0.2, 0.9, 0.8, 0.99);
518 gPad->SetLogy(true);
519 }
520
521 {
522 TCanvas* c = fHM->CreateCanvas("richqa_RTDist_xy_truematch" + opt, "richqa_RTDist_xy_truematch" + opt, 1800, 600);
523 c->Divide(3, 1);
524 c->cd(1);
525 DrawH3Profile(fHM->H3("fhRTDistVsXYTruematch" + opt), true, false, 0., .4);
526 c->cd(2);
527 DrawH2WithProfile(fHM->H2("fhRTDistVsXTruematch" + opt), false, true);
528 fHM->H2("fhRTDistVsXTruematch" + opt)->GetYaxis()->SetRangeUser(0., 2.);
529 c->cd(3);
530 DrawH2WithProfile(fHM->H2("fhRTDistVsYTruematch" + opt), false, true);
531 fHM->H2("fhRTDistVsYTruematch" + opt)->GetYaxis()->SetRangeUser(0., 2.);
532 }
533}
534
535bool CbmRichRecoQa::IsMcPrimaryElectron(const CbmMCTrack* mctrack)
536{
537 if (mctrack == nullptr) return false;
538 if (mctrack->GetGeantProcessId() == kPPrimary && std::abs(mctrack->GetPdgCode()) == 11) return true;
539 return false;
540}
541
542bool CbmRichRecoQa::IsMcPion(const CbmMCTrack* mctrack)
543{
544 if (mctrack == nullptr || std::abs(mctrack->GetPdgCode()) != 211) return false;
545 return true;
546}
547
548void CbmRichRecoQa::Finish()
549{
550 DrawMismatchSrc();
551
552 TDirectory* oldir = gDirectory;
553 TFile* outFile = FairRootManager::Instance()->GetOutFile();
554 if (outFile != nullptr) {
555 outFile->mkdir(GetName());
556 outFile->cd(GetName());
557 fHM->WriteToFile();
558 }
559
560 DrawHist();
561 fHM->SaveCanvasToImage(fOutputDir);
562 fHM->Clear();
563 gDirectory->cd(oldir->GetPath());
564}
565
566void CbmRichRecoQa::DrawFromFile(const string& fileName, const string& outputDir)
567{
568 fOutputDir = outputDir;
569
571 TFile* oldFile = gFile;
572 TDirectory* oldDir = gDirectory;
573
574 if (fHM != nullptr) delete fHM;
575
576 fHM = new CbmHistManager();
577 TFile* file = new TFile(fileName.c_str());
578 fHM->ReadFromFile(file);
579
580 DrawMismatchSrc();
581 DrawHist();
582
583 fHM->SaveCanvasToImage(fOutputDir);
584
586 gFile = oldFile;
587 gDirectory = oldDir;
588}
589
590ClassImp(CbmRichRecoQa)
ClassImp
ClassImp(CbmConverterManager)
DrawH2WithProfile
void DrawH2WithProfile(TH2 *hist, Bool_t doGaussFit, Bool_t drawOnlyMean, const string &drawOpt2D, Int_t profileColor, Int_t profileLineWidth)
Definition CbmDrawHist.cxx:247
DrawTextOnPad
void DrawTextOnPad(const string &text, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition CbmDrawHist.cxx:212
SetDefaultDrawStyle
void SetDefaultDrawStyle()
Definition CbmDrawHist.cxx:37
DrawH1
void DrawH1(TH1 *hist, HistScale logx, HistScale logy, const string &drawOpt, Int_t color, Int_t lineWidth, Int_t lineStyle, Int_t markerSize, Int_t markerStyle)
Definition CbmDrawHist.cxx:54
DrawH3Profile
TH2D * DrawH3Profile(TH3 *h, Bool_t drawMean, Bool_t doGaussFit, Double_t zMin, Double_t zMax)
Definition CbmDrawHist.cxx:297
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
kLinear
@ kLinear
Definition CbmDrawHist.h:69
kLog
@ kLog
Definition CbmDrawHist.h:68
CbmHistManager.h
Histogram manager.
CbmMatchRecoToMC.h
FairTask for matching RECO data to MC.
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition CbmDigiManager.h:126
CbmGlobalTrack::GetStsTrackIndex
int32_t GetStsTrackIndex() const
Definition CbmGlobalTrack.h:44
CbmGlobalTrack::GetRichRingIndex
int32_t GetRichRingIndex() const
Definition CbmGlobalTrack.h:47
CbmHistManager::CreateCanvas
TCanvas * CreateCanvas(const std::string &name, const std::string &title, Int_t width, Int_t height)
Create and draw TCanvas and store pointer to it.
Definition core/base/CbmHistManager.cxx:303
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition CbmHistManager.h:243
CbmMCTrack::GetCharge
double GetCharge() const
Charge of the associated particle.
Definition CbmMCTrack.cxx:150
CbmMCTrack::GetP
double GetP() const
Definition CbmMCTrack.h:97
CbmMCTrack::GetGeantProcessId
uint32_t GetGeantProcessId() const
Definition CbmMCTrack.h:66
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:67
CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit
static std::vector< CbmLink > GetMcTrackMotherIdsForRichHit(CbmDigiManager *digiMan, const CbmRichHit *hit, CbmMCDataArray *richPoints, CbmMCDataArray *mcTracks)
Get CbmLinks of Mother MC Tracks for RICH hit.
Definition CbmMatchRecoToMC.cxx:901
CbmMatch::GetNofLinks
int32_t GetNofLinks() const
Definition CbmMatch.h:42
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
CbmRichRecoQa::fRichProjections
TClonesArray * fRichProjections
Definition alignment/CbmRichRecoQa.h:154
CbmRichRecoQa::RingTrackMismatchSource
void RingTrackMismatchSource()
Fill histograms related to study of the source of ring-track mismatch.
Definition alignment/CbmRichRecoQa.cxx:281
CbmRichRecoQa::fGlobalTracks
TClonesArray * fGlobalTracks
Definition alignment/CbmRichRecoQa.h:151
CbmRichRecoQa::fNofHitsInRingMap
std::map< Int_t, Int_t > fNofHitsInRingMap
Definition alignment/CbmRichRecoQa.h:157
CbmRichRecoQa::DrawRingTrackDist
void DrawRingTrackDist(const std::string &opt)
Draw histograms related to ring-track distance for pions or electrons (+/-).
Definition qa/CbmRichRecoQa.cxx:490
CbmRichRecoQa::fMcTracks
CbmMCDataArray * fMcTracks
Definition qa/CbmRichRecoQa.h:148
CbmRichRecoQa::FillRingTrackDistance
void FillRingTrackDistance()
Fill histogramms related to ring track distance.
Definition alignment/CbmRichRecoQa.cxx:444
CbmRichRecoQa::InitHistograms
void InitHistograms()
Initialize histograms.
Definition alignment/CbmRichRecoQa.cxx:128
CbmRichRecoQa::GetMeanRmsOverflowString
string GetMeanRmsOverflowString(TH1 *h, Bool_t withOverflow=true)
Return string with mean, RMS and overflow percent for input TH1.
Definition alignment/CbmRichRecoQa.cxx:707
CbmRichRecoQa::DrawMismatchSrc
void DrawMismatchSrc()
Draw MismatchSrc histogram and canvas.
Definition qa/CbmRichRecoQa.cxx:361
CbmRichRecoQa::IsMcPrimaryElectron
static Bool_t IsMcPrimaryElectron(const CbmMCTrack *mctrack)
Definition alignment/CbmRichRecoQa.cxx:865
CbmRichRecoQa::fRichRings
TClonesArray * fRichRings
Definition alignment/CbmRichRecoQa.h:149
CbmRichRecoQa::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition alignment/CbmRichRecoQa.cxx:63
CbmRichRecoQa::HasRichProjection
bool HasRichProjection(Int_t stsTrackId)
Definition alignment/CbmRichRecoQa.cxx:426
CbmRichRecoQa::WasRingMatched
bool WasRingMatched(Int_t mcTrackId)
Definition alignment/CbmRichRecoQa.cxx:407
CbmRichRecoQa::fRichRingMatches
TClonesArray * fRichRingMatches
Definition alignment/CbmRichRecoQa.h:150
CbmRichRecoQa::DrawFromFile
void DrawFromFile(const string &fileName, const string &outputDir)
Draw histogram from file.
Definition alignment/CbmRichRecoQa.cxx:897
CbmRichRecoQa::fStsTracks
TClonesArray * fStsTracks
Definition alignment/CbmRichRecoQa.h:152
CbmRichRecoQa::fDigiMan
CbmDigiManager * fDigiMan
Definition qa/CbmRichRecoQa.h:157
CbmRichRecoQa::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition alignment/CbmRichRecoQa.h:153
CbmRichRecoQa::fEventList
CbmMCEventList * fEventList
Definition qa/CbmRichRecoQa.h:158
CbmRichRecoQa::WasRingFound
bool WasRingFound(Int_t mcTrackId)
Definition alignment/CbmRichRecoQa.cxx:393
CbmRichRecoQa::IsMcPion
static Bool_t IsMcPion(const CbmMCTrack *mctrack)
Definition alignment/CbmRichRecoQa.cxx:873
CbmRichRecoQa::fRichPoints
TClonesArray * fRichPoints
Definition alignment/CbmRichRecoQa.h:146
CbmRichRecoQa::DrawHist
void DrawHist()
Draw histograms.
Definition alignment/CbmRichRecoQa.cxx:555
CbmRichRecoQa::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition alignment/CbmRichRecoQa.cxx:256
CbmRichRecoQa::CbmRichRecoQa
CbmRichRecoQa()
Standard constructor.
Definition alignment/CbmRichRecoQa.cxx:42
CbmRichRecoQa::Finish
virtual void Finish()
Inherited from FairTask.
Definition alignment/CbmRichRecoQa.cxx:881
CbmRichRecoQa::FillRichRingNofHits
void FillRichRingNofHits()
Fill map mcTrackId -> nof RICH hits.
Definition alignment/CbmRichRecoQa.cxx:266
CbmRichRing::GetNofHits
int32_t GetNofHits() const
Definition CbmRichRing.h:37
CbmRichRing::GetCenterX
float GetCenterX() const
Definition CbmRichRing.h:77
CbmRichRing::GetCenterY
float GetCenterY() const
Definition CbmRichRing.h:78
CbmRichUtil::GetRingTrackDistance
static double GetRingTrackDistance(int globalTrackId)
Definition CbmRichUtil.cxx:24
h
Data class with information on a STS local track.
Cbm::GetOrFatal
T * GetOrFatal(const std::string &objName, const std::string &description="")
Definition CbmUtils.h:60
Cbm::InitOrFatalMc
CbmMCDataArray * InitOrFatalMc(const std::string &objName, const std::string &description)
Definition CbmUtils.cxx:27
Cbm::NumberToString
std::string NumberToString(const T &value, int precision=1)
Definition CbmUtils.h:34
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:22