CbmRoot
Loading...
Searching...
No Matches
CbmRichRecoTbQa.cxx
Go to the documentation of this file.
1/* Copyright (C) 2018-2021 GSI/JINR-LIT, Darmstadt/Dubna
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Andrey Lebedev [committer] */
4
5#include "CbmRichRecoTbQa.h"
6
7#include "CbmDigiManager.h"
8#include "CbmDrawHist.h"
9#include "CbmGlobalTrack.h"
10#include "CbmHistManager.h"
11#include "CbmMCDataArray.h"
12#include "CbmMCDataManager.h"
13#include "CbmMCEventList.h"
14#include "CbmMCTrack.h"
15#include "CbmMatchRecoToMC.h"
16#include "CbmRichDigi.h"
17#include "CbmRichDraw.h"
18#include "CbmRichGeoManager.h"
19#include "CbmRichHit.h"
20#include "CbmRichPoint.h"
21#include "CbmRichRing.h"
22#include "CbmRichUtil.h"
23#include "CbmStsPoint.h"
24#include "CbmTrackMatchNew.h"
25#include "CbmUtils.h"
26#include "FairEventHeader.h"
27#include "FairMCPoint.h"
28#include "FairRunAna.h"
29#include "FairRunSim.h"
30#include "FairTrackParam.h"
31#include "TCanvas.h"
32#include "TClonesArray.h"
33#include "TEllipse.h"
34#include "TF1.h"
35#include "TH1.h"
36#include "TH1D.h"
37#include "TMCProcess.h"
38#include "TMarker.h"
39#include "TStyle.h"
40#include "elid/CbmLitGlobalElectronId.h"
41
42#include <Logger.h>
43
44#include <TFile.h>
45
46#include <iostream>
47#include <sstream>
48#include <string>
49
50using namespace std;
51
52CbmRichRecoTbQa::CbmRichRecoTbQa()
53 : FairTask("CbmRichRecoTbQa")
54 , fHM(nullptr)
55 , fTimeSliceNum(0.)
56 , fNofLogEvents(150)
57 , fOutputDir("")
58 , fMCTracks(nullptr)
59 , fRichPoints(nullptr)
60 , fStsPoints(nullptr)
61 , fDigiMan(nullptr)
62 , fRichHits(nullptr)
63 , fRichRings(nullptr)
64 , fRichRingMatches(nullptr)
65 , fEventList(nullptr)
66 , fRecRings()
67{
68}
69
70
71InitStatus CbmRichRecoTbQa::Init()
72{
73 cout << "CbmRichRecoTbQa::Init" << endl;
74 FairRootManager* ioman = FairRootManager::Instance();
75 if (nullptr == ioman) {
76 LOG(fatal) << GetName() << "::Init: RootManager not instantised!";
77 }
78
79 CbmMCDataManager* mcManager = (CbmMCDataManager*) ioman->GetObject("MCDataManager");
80 if (mcManager == nullptr) {
81 LOG(fatal) << GetName() << "::Init: No MCDataManager!";
82 }
83
84 fMCTracks = mcManager->InitBranch("MCTrack");
85 fRichPoints = mcManager->InitBranch("RichPoint");
86 fStsPoints = mcManager->InitBranch("StsPoint");
87
88 fDigiMan = CbmDigiManager::Instance();
89 fDigiMan->Init();
90
91 fRichHits = (TClonesArray*) ioman->GetObject("RichHit");
92 if (nullptr == fRichHits) {
93 LOG(fatal) << GetName() << "::Init: No Rich hits!";
94 }
95
96 fRichRings = (TClonesArray*) ioman->GetObject("RichRing");
97 if (nullptr == fRichRings) {
98 LOG(fatal) << GetName() << "::Init: No Rich rings!";
99 }
100
101 fRichRingMatches = (TClonesArray*) ioman->GetObject("RichRingMatch");
102 if (nullptr == fRichRingMatches) {
103 LOG(fatal) << GetName() << "::Init: No RichRingMatch!";
104 }
105
106 fEventList = (CbmMCEventList*) ioman->GetObject("MCEventList.");
107 if (nullptr == fEventList) {
108 LOG(fatal) << GetName() << "::Init: No MCEventList!";
109 }
110
111 InitHistograms();
112
113 return kSUCCESS;
114}
115
116void CbmRichRecoTbQa::InitHistograms()
117{
118 fHM = new CbmHistManager();
119
120 Int_t nBinsLog = 20000;
121 Double_t minLog = 0.;
122 Double_t maxLog = 10000.;
123
124
125 fHM->Create1<TH1D>("fhNofRichDigisPerTS", "fhNofRichDigisPerTS;RICH digis per time slice;Yield", 100, -1, -1.);
126 fHM->Create1<TH1D>("fhNofRichHitsPerTS", "fhNofRichHitsPerTS;RICH hits per time slice;Yield", 100, -1., -1.);
127 fHM->Create1<TH1D>("fhNofRichRingsPerTS", "fhNofRichRingsPerTS;RICH rings per time slice;Yield", 100, -1., -1.);
128
129 fHM->Create1<TH1D>("fhRichDigiTimeLog", "fhRichDigiTimeLog;RICH digi time [ns];Yield", nBinsLog, minLog, maxLog);
130
131 fHM->Create1<TH1D>("fhRichRingTimeLog", "fhRichRingTimeLog;RICH ring time [ns];Yield", nBinsLog, minLog, maxLog);
132
133 fHM->Create1<TH1D>("fhRichPointTime", "fhRichPointTime;RICH MC point time [ns];Yield", 100, 0., 100.);
134 fHM->Create1<TH1D>("fhRichPointTimeChPhoton", "fhRichPointTimeChPhoton;RICH MC point time [ns];Yield", 100, 0., 100.);
135 fHM->Create1<TH1D>("fhRichPointTimeNotChPhoton", "fhRichPointTimeNotChPhoton;RICH MC point time [ns];Yield", 100, 0.,
136 100.);
137 fHM->Create1<TH1D>("fhRichPointTimePrimEl", "fhRichPointTimePrimEl;RICH MC point time [ns];Yield", 100, 0., 100.);
138 fHM->Create1<TH1D>("fhRichPointTimeSecEl", "fhRichPointTimeSecEl;RICH MC point time [ns];Yield", 100, 0., 100.);
139 fHM->Create1<TH1D>("fhRichPointTimeOther", "fhRichPointTimeOther;RICH MC point time [ns];Yield", 100, 0., 100.);
140 fHM->Create1<TH1D>("fhRichPointTimePion", "fhRichPointTimePion;RICH MC point time [ns];Yield", 100, 0., 100.);
141
142 fHM->Create1<TH1D>("fhStsPointTime", "fhStsPointTime;STS MC point time [ns];Yield", 400, 0., 400.);
143
144 // -1 for noise hits
145 for (int iEv = -1; iEv < fNofLogEvents; iEv++) {
146 string richDigiStr = "fhRichDigiTimeLog_" + to_string(iEv);
147 fHM->Create1<TH1D>(richDigiStr, richDigiStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
148
149 string richRingStr = "fhRichRingTimeLog_" + to_string(iEv);
150 fHM->Create1<TH1D>(richRingStr, richRingStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
151
152 string richStr = "fhRichPointTimeLog_" + to_string(iEv);
153 fHM->Create1<TH1D>(richStr, richStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
154
155 string stsStr = "fhStsPointTimeLog_" + to_string(iEv);
156 fHM->Create1<TH1D>(stsStr, stsStr + ";Time [ns];Yield", nBinsLog, minLog, maxLog);
157 }
158
159 fHM->Create1<TH1D>("fhTimeBetweenEvents", "fhTimeBetweenEvents;Time between events [ns];Yield", 100, 0., 500.);
160
161
162 // efficiency histograms
163 fHM->Create1<TH1D>("fhMomElAcc", "fhMomElAcc;P [GeV/c];Yield", 10, 0., 10.);
164 fHM->Create1<TH1D>("fhMomElRec", "fhMomElRec;P [GeV/c];Yield", 10, 0., 10.);
165 fHM->Create1<TH1D>("fhMomRefElAcc", "fhMomRefElAcc;P [GeV/c];Yield", 10, 0., 10.);
166 fHM->Create1<TH1D>("fhMomRefElRec", "fhMomRefElRec;P [GeV/c];Yield", 10, 0., 10.);
167
168
169 fHM->Create1<TH1D>("fhNofHitsElAcc", "fhNofHitsElAcc;Nof hits in ring;Yield", 20, -.5, 39.5);
170 fHM->Create1<TH1D>("fhNofHitsElRec", "fhNofHitsElRec;Nof hits in ring;Yield", 20, -.5, 39.5);
171 fHM->Create1<TH1D>("fhNofHitsRefElAcc", "fhNofHitsRefElAcc;Nof hits in ring;Yield", 20, -.5, 39.5);
172 fHM->Create1<TH1D>("fhNofHitsRefElRec", "fhNofHitsRefElRec;Nof hits in ring;Yield", 20, -.5, 39.5);
173
174 fHM->Create1<TH1D>("fhEventMultElAcc", "fhEventMultElAcc;Nof MC tracks;Yield", 10, 0., 2000);
175 fHM->Create1<TH1D>("fhEventMultElRec", "fhEventMultElRec;Nof MC tracks;Yield", 10, 0, 2000);
176}
177
178void CbmRichRecoTbQa::Exec(Option_t* /*option*/)
179{
180 fTimeSliceNum++;
181 int nofMcEvents = fEventList->GetNofEvents();
182
183 cout << "CbmRichRecoTbQa, time slice:" << fTimeSliceNum << " nofMcEvents:" << nofMcEvents << endl;
184
185 Process();
186
187 RingRecoEfficiency();
188}
189
190void CbmRichRecoTbQa::Process()
191{
192 Int_t nofRichDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
193 Int_t nofRichHits = fRichHits->GetEntriesFast();
194 Int_t nofRichRings = fRichRings->GetEntriesFast();
195 fHM->H1("fhNofRichDigisPerTS")->Fill(nofRichDigis);
196 fHM->H1("fhNofRichHitsPerTS")->Fill(nofRichHits);
197 fHM->H1("fhNofRichRingsPerTS")->Fill(nofRichRings);
198
199 LOG(debug) << "nofRichDigis:" << nofRichDigis;
200 for (Int_t iDigi = 0; iDigi < nofRichDigis; iDigi++) {
201 const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(iDigi);
202 const CbmMatch* digiMatch = fDigiMan->GetMatch(ECbmModuleId::kRich, iDigi);
203 fHM->H1("fhRichDigiTimeLog")->Fill(digi->GetTime());
204
205 Int_t eventNum = digiMatch->GetMatchedLink().GetEntry();
206 Int_t index = digiMatch->GetMatchedLink().GetIndex();
207
208 if (eventNum < 0 || index < 0) {
209 fHM->H1("fhRichDigiTimeLog_-1")->Fill(digi->GetTime());
210 }
211 else {
212 if (eventNum < fNofLogEvents) {
213 string hName = "fhRichDigiTimeLog_" + to_string(eventNum);
214 fHM->H1(hName)->Fill(digi->GetTime());
215 }
216 }
217 }
218
219 for (Int_t iR = 0; iR < nofRichRings; iR++) {
220 CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
221 CbmTrackMatchNew* ringMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(iR));
222 if (ring == nullptr || ringMatch == nullptr) continue;
223
224 fHM->H1("fhRichRingTimeLog")->Fill(ring->GetTime());
225
226 Int_t eventNum = ringMatch->GetMatchedLink().GetEntry();
227 Int_t index = ringMatch->GetMatchedLink().GetIndex();
228
229 if (eventNum < 0 || index < 0) {
230 fHM->H1("fhRichRingTimeLog_-1")->Fill(ring->GetTime());
231 }
232 else {
233 if (eventNum < fNofLogEvents) {
234 string hName = "fhRichRingTimeLog_" + to_string(eventNum);
235 fHM->H1(hName)->Fill(ring->GetTime());
236 }
237 }
238 }
239}
240
241void CbmRichRecoTbQa::ProcessMc()
242{
243 int fileId = 0;
244
245 // int nMCEvents = fEventList->GetNofEvents();
246 for (int iEv = 0; fEventList->GetEventTime(iEv, fileId) > 0.; iEv++) {
247 double dT = (iEv == 0) ? fEventList->GetEventTime(iEv, fileId)
248 : fEventList->GetEventTime(iEv, fileId) - fEventList->GetEventTime(iEv - 1, fileId);
249 fHM->H1("fhTimeBetweenEvents")->Fill(dT);
250 }
251
252 for (Int_t iEv = 0; fRichPoints->Size(fileId, iEv) >= 0; iEv++) {
253 Int_t nofRichPoints = fRichPoints->Size(fileId, iEv);
254 for (Int_t iP = 0; iP < nofRichPoints; iP++) {
255 CbmRichPoint* point = (CbmRichPoint*) fRichPoints->Get(fileId, iEv, iP);
256 fHM->H1("fhRichPointTime")->Fill(point->GetTime());
257 if (IsCherenkovPhoton(point, fileId, iEv)) {
258 fHM->H1("fhRichPointTimeChPhoton")->Fill(point->GetTime());
259 }
260 else {
261 fHM->H1("fhRichPointTimeNotChPhoton")->Fill(point->GetTime());
262 }
263
264 if (IsCherenkovPhotonFromPrimaryElectron(point, fileId, iEv)) {
265 fHM->H1("fhRichPointTimePrimEl")->Fill(point->GetTime());
266 }
267 else if (IsCherenkovPhotonFromPion(point, fileId, iEv)) {
268 fHM->H1("fhRichPointTimePion")->Fill(point->GetTime());
269 }
270 if (IsCherenkovPhotonFromSecondaryElectron(point, fileId, iEv)) {
271 fHM->H1("fhRichPointTimeSecEl")->Fill(point->GetTime());
272 }
273 else {
274 fHM->H1("fhRichPointTimeOther")->Fill(point->GetTime());
275 }
276
277 if (iEv < fNofLogEvents) {
278 string richStr = "fhRichPointTimeLog_" + to_string(iEv);
279 double eventTime = fEventList->GetEventTime(iEv, fileId);
280 fHM->H1(richStr)->Fill(eventTime + point->GetTime());
281 }
282 }
283 }
284
285
286 for (Int_t iEv = 0; fStsPoints->Size(fileId, iEv) >= 0; iEv++) {
287 Int_t nofStsPoints = fStsPoints->Size(fileId, iEv);
288 for (Int_t j = 0; j < nofStsPoints; j++) {
289 const CbmStsPoint* point = (CbmStsPoint*) fStsPoints->Get(fileId, iEv, j);
290 fHM->H1("fhStsPointTime")->Fill(point->GetTime());
291 if (iEv < fNofLogEvents) {
292 string stsStr = "fhStsPointTimeLog_" + to_string(iEv);
293 double eventTime = fEventList->GetEventTime(iEv, fileId);
294 fHM->H1(stsStr)->Fill(eventTime + point->GetTime());
295 }
296 }
297 }
298}
299
300Int_t CbmRichRecoTbQa::GetNofPrimaryMcTracks(Int_t iEv)
301{
302 Int_t fileId = 0;
303 Int_t counter = 0;
304 Int_t nofMcTracks = fMCTracks->Size(fileId, iEv);
305 for (Int_t j = 0; j < nofMcTracks; j++) {
306 const CbmMCTrack* track = static_cast<CbmMCTrack*>(fMCTracks->Get(fileId, iEv, j));
307 if (track->GetGeantProcessId() == kPPrimary) counter++;
308 }
309 return counter;
310}
311
312void CbmRichRecoTbQa::RingRecoEfficiency()
313{
314 map<CbmLink, Int_t> nofHitsInRing;
315 Int_t nofRichHits = fRichHits->GetEntriesFast();
316 for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
317 const CbmRichHit* hit = static_cast<const CbmRichHit*>(fRichHits->At(iHit));
318 if (nullptr == hit) continue;
319 vector<CbmLink> motherIds = CbmMatchRecoToMC::GetMcTrackMotherIdsForRichHit(fDigiMan, hit, fRichPoints, fMCTracks);
320 for (const auto& motherId : motherIds) {
321 nofHitsInRing[motherId]++;
322 }
323 }
324
325 Int_t fileId = 0;
326 for (Int_t iEv = 0; fMCTracks->Size(fileId, iEv) >= 0; iEv++) {
327 Int_t nofMcTracks = fMCTracks->Size(fileId, iEv);
328 Int_t nofPrimMcTracks = GetNofPrimaryMcTracks(iEv);
329 for (Int_t j = 0; j < nofMcTracks; j++) {
330 CbmLink val(1., j, iEv, fileId);
331 const CbmMCTrack* track = static_cast<CbmMCTrack*>(fMCTracks->Get(val));
332 Int_t nofRichHitsInRing = nofHitsInRing[val];
333 Double_t mom = track->GetP();
334 if (nofRichHitsInRing >= 7 && IsMcPrimaryElectron(track)) {
335 fHM->H1("fhMomElAcc")->Fill(mom);
336 fHM->H1("fhNofHitsElAcc")->Fill(nofRichHitsInRing);
337 fHM->H1("fhEventMultElAcc")->Fill(nofPrimMcTracks);
338 if (nofRichHitsInRing >= 15) {
339 fHM->H1("fhMomRefElAcc")->Fill(mom);
340 fHM->H1("fhNofHitsRefElAcc")->Fill(nofRichHitsInRing);
341 }
342 }
343 }
344 }
345
346 Int_t nofRings = fRichRings->GetEntriesFast();
347 for (Int_t iRing = 0; iRing < nofRings; iRing++) {
348 const CbmTrackMatchNew* richRingMatch = static_cast<const CbmTrackMatchNew*>(fRichRingMatches->At(iRing));
349 Bool_t isRichOk = richRingMatch->GetTrueOverAllHitsRatio() >= 0.6;
350
351 const CbmMCTrack* track = static_cast<CbmMCTrack*>(
352 fMCTracks->Get(fileId, richRingMatch->GetMatchedLink().GetEntry(), richRingMatch->GetMatchedLink().GetIndex()));
353 Int_t nofPrimMcTracks = GetNofPrimaryMcTracks(richRingMatch->GetMatchedLink().GetEntry());
354 auto val = richRingMatch->GetMatchedLink();
355 Int_t nofRichHitsInRing = nofHitsInRing[val];
356 Double_t mom = track->GetP();
357 Bool_t isClone = (std::find(fRecRings.begin(), fRecRings.end(), val) != fRecRings.end());
358 if (nofRichHitsInRing >= 7 && isRichOk && IsMcPrimaryElectron(track) && !isClone) {
359 fHM->H1("fhMomElRec")->Fill(mom);
360 fHM->H1("fhNofHitsElRec")->Fill(nofRichHitsInRing);
361 fHM->H1("fhEventMultElRec")->Fill(nofPrimMcTracks);
362 fRecRings.push_back(val);
363 if (nofRichHitsInRing >= 15) {
364 fHM->H1("fhMomRefElRec")->Fill(mom);
365 fHM->H1("fhNofHitsRefElRec")->Fill(nofRichHitsInRing);
366 }
367 }
368 }
369}
370
371
372Bool_t CbmRichRecoTbQa::IsCherenkovPhoton(const CbmRichPoint* point, Int_t fileId, Int_t eventId)
373{
374 if (point == nullptr) return false;
375 Int_t trackId = point->GetTrackID();
376 if (trackId < 0) return false;
377 CbmMCTrack* p = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
378 return (p != nullptr && TMath::Abs(p->GetPdgCode()) == 50000050);
379}
380
381Bool_t CbmRichRecoTbQa::IsCherenkovPhotonFromPrimaryElectron(const CbmRichPoint* point, Int_t fileId, Int_t eventId)
382{
383 if (point == nullptr) return false;
384 Int_t trackId = point->GetTrackID();
385 if (trackId < 0) return false;
386 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
387 if (mcTrack == nullptr || TMath::Abs(mcTrack->GetPdgCode()) != 50000050) return false;
388
389 Int_t motherId = mcTrack->GetMotherId();
390 if (motherId < 0) return false;
391 CbmMCTrack* mcMotherTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, motherId);
392 return IsMcPrimaryElectron(mcMotherTrack);
393}
394
395Bool_t CbmRichRecoTbQa::IsCherenkovPhotonFromSecondaryElectron(const CbmRichPoint* point, Int_t fileId, Int_t eventId)
396{
397 if (point == nullptr) return false;
398 Int_t trackId = point->GetTrackID();
399 if (trackId < 0) return false;
400 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
401 if (mcTrack == nullptr || TMath::Abs(mcTrack->GetPdgCode()) != 50000050) return false;
402
403 Int_t motherId = mcTrack->GetMotherId();
404 if (motherId < 0) return false;
405 CbmMCTrack* mcMotherTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, motherId);
406 if (mcMotherTrack == nullptr) return false;
407 int pdg = TMath::Abs(mcMotherTrack->GetPdgCode());
408 if (mcMotherTrack->GetGeantProcessId() != kPPrimary && pdg == 11) return true;
409
410 return false;
411}
412
413Bool_t CbmRichRecoTbQa::IsMcPrimaryElectron(const CbmMCTrack* mctrack)
414{
415 if (mctrack == nullptr) return false;
416 int pdg = TMath::Abs(mctrack->GetPdgCode());
417 if (mctrack->GetGeantProcessId() == kPPrimary && pdg == 11) return true;
418 return false;
419}
420
421Bool_t CbmRichRecoTbQa::IsMcPion(const CbmMCTrack* mctrack)
422{
423 if (mctrack == nullptr) return false;
424 int pdg = TMath::Abs(mctrack->GetPdgCode());
425 if (pdg == 211) return true;
426 return false;
427}
428
429Bool_t CbmRichRecoTbQa::IsCherenkovPhotonFromPion(const CbmRichPoint* point, Int_t fileId, Int_t eventId)
430{
431 if (point == nullptr) return false;
432 Int_t trackId = point->GetTrackID();
433 if (trackId < 0) return false;
434 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, trackId);
435 if (mcTrack == nullptr || TMath::Abs(mcTrack->GetPdgCode()) != 50000050) return false;
436
437 Int_t motherId = mcTrack->GetMotherId();
438 if (motherId < 0) return false;
439 CbmMCTrack* mcMotherTrack = (CbmMCTrack*) fMCTracks->Get(fileId, eventId, motherId);
440 return IsMcPion(mcMotherTrack);
441}
442
443void CbmRichRecoTbQa::DrawHist()
444{
445 SetDefaultDrawStyle();
446
447
448 {
449 fHM->CreateCanvas("rich_tb_fhTimeBetweenEvents", "rich_tb_fhTimeBetweenEvents", 800, 800);
450 DrawH1(fHM->H1("fhTimeBetweenEvents"), kLinear, kLog);
451 Double_t min = 0.;
452 Double_t max = 10.;
453 Double_t partIntegral =
454 fHM->H1("fhTimeBetweenEvents")
455 ->Integral(fHM->H1("fhTimeBetweenEvents")->FindBin(min), fHM->H1("fhTimeBetweenEvents")->FindBin(max));
456 Double_t allIntegral = fHM->H1("fhTimeBetweenEvents")->Integral();
457 Double_t ratio = 100. * partIntegral / allIntegral;
458 cout << ratio << "% of all time between events are in range (" << min << " ," << max << ") ns" << endl;
459 }
460
461
462 {
463 fHM->CreateCanvas("rich_tb_fhRichPointTime", "rich_tb_fhRichPointTime", 800, 800);
464 DrawH1(fHM->H1("fhRichPointTime"), kLinear, kLog);
465 }
466
467 {
468 TCanvas* c = fHM->CreateCanvas("rich_tb_fhRichPointTimeSources", "rich_tb_fhRichPointTimeSources", 1600, 800);
469 c->Divide(2, 1);
470 c->cd(1);
471 DrawH1({fHM->H1("fhRichPointTimeChPhoton"), fHM->H1("fhRichPointTimeNotChPhoton")},
472 {"Ch. Photons", "Not Ch. Photons"}, kLinear, kLog);
473
474 c->cd(2);
475 DrawH1({fHM->H1("fhRichPointTimeSecEl"), fHM->H1("fhRichPointTimePrimEl"), fHM->H1("fhRichPointTimePion"),
476 fHM->H1("fhRichPointTimeOther")},
477 {"e^{#pm}_{sec}", "e^{#pm}_{prim}", "#pi^{#pm}", "Other"}, kLinear, kLog);
478 }
479
480
481 {
482 fHM->CreateCanvas("rich_tb_fhStsPointTime", "rich_tb_fhStsPointTime", 800, 800);
483 DrawH1(fHM->H1("fhStsPointTime"), kLinear, kLog);
484 }
485
486 {
487 fHM->CreateCanvas("rich_tb_fhRichPointTimeLog", "rich_tb_fhRichPointTimeLog", 1600, 600);
488 DrawTimeLog("fhRichPointTimeLog", fNofLogEvents);
489 }
490
491 {
492 fHM->CreateCanvas("rich_tb_fhStsPointTimeLog", "rich_tb_fhStsPointTimeLog", 1600, 600);
493 DrawTimeLog("fhStsPointTimeLog", fNofLogEvents);
494 }
495
496 {
497 fHM->CreateCanvas("rich_tb_reco_fhRichDigiTimeLogAll", "rich_tb_reco_fhRichDigiTimeLogAll", 1600, 600);
498 DrawH1(fHM->H1("fhRichDigiTimeLog"), kLinear, kLog);
499 gPad->SetLeftMargin(0.07);
500 gPad->SetRightMargin(0.05);
501 }
502
503 {
504 fHM->CreateCanvas("rich_tb_reco_fhRichDigiTimeLog", "rich_tb_reco_fhRichDigiTimeLog", 1600, 600);
505 DrawTimeLog("fhRichDigiTimeLog", fNofLogEvents, true);
506 }
507
508 {
509 fHM->CreateCanvas("rich_tb_reco_fhRichRingTimeLogAll", "rich_tb_reco_fhRichRingTimeLogAll", 1600, 600);
510 DrawH1(fHM->H1("fhRichRingTimeLog"), kLinear, kLog);
511 gPad->SetLeftMargin(0.07);
512 gPad->SetRightMargin(0.05);
513 }
514
515 {
516 fHM->CreateCanvas("rich_tb_reco_timelog_hits_rings", "rich_tb_reco_timelog_hits_rings", 1600, 600);
517 DrawH1({fHM->H1("fhRichDigiTimeLog"), fHM->H1("fhRichRingTimeLog")}, {"Hits", "Rings"});
518 gPad->SetLeftMargin(0.07);
519 gPad->SetRightMargin(0.05);
520 }
521
522 {
523 fHM->CreateCanvas("rich_tb_reco_fhRichRingTimeLog", "rich_tb_reco_fhRichRingTimeLog", 1600, 600);
524 DrawTimeLog("fhRichRingTimeLog", fNofLogEvents, true);
525 }
526
527 {
528 TCanvas* c =
529 fHM->CreateCanvas("rich_tb_reco_fhRichObectsPerTimeSlice", "rich_tb_reco_fhRichObectsPerTimeSlice", 1500, 500);
530 c->Divide(3, 1);
531 c->cd(1);
532 DrawH1(fHM->H1("fhNofRichDigisPerTS"), kLinear, kLog);
533 c->cd(2);
534 DrawH1(fHM->H1("fhNofRichHitsPerTS"), kLinear, kLog);
535 c->cd(3);
536 DrawH1(fHM->H1("fhNofRichRingsPerTS"), kLinear, kLog);
537 }
538
539 {
540 TCanvas* c = fHM->CreateCanvas("rich_tb_reco_acc_rec", "rich_tb_reco_acc_rec", 1600, 600);
541 c->Divide(2, 1);
542 c->cd(1);
543 Int_t nofAccMom = fHM->H1("fhMomElAcc")->GetEntries();
544 Int_t nofRecMom = fHM->H1("fhMomElRec")->GetEntries();
545 Int_t nofRefAccMom = fHM->H1("fhMomRefElAcc")->GetEntries();
546 Int_t nofRefRecMom = fHM->H1("fhMomRefElRec")->GetEntries();
547 Double_t effMom = 100. * (Double_t) nofRecMom / (Double_t) nofAccMom;
548 Double_t effRefMom = 100. * (Double_t) nofRefRecMom / (Double_t) nofRefAccMom;
549 DrawH1({fHM->H1("fhMomElAcc"), fHM->H1("fhMomElRec")},
550 {"Acc (" + to_string(nofAccMom) + ")", "Rec (" + to_string(nofRecMom) + ")"});
551 c->cd(2);
552 Int_t nofAccNofHits = fHM->H1("fhNofHitsElAcc")->GetEntries();
553 Int_t nofRecNofHits = fHM->H1("fhNofHitsElRec")->GetEntries();
554 Int_t nofRefAccNofHits = fHM->H1("fhNofHitsRefElAcc")->GetEntries();
555 Int_t nofRefRecNofHits = fHM->H1("fhNofHitsRefElRec")->GetEntries();
556 Double_t effNofHits = 100. * (Double_t) nofRecNofHits / (Double_t) nofAccNofHits;
557 Double_t effRefNofHits = 100. * (Double_t) nofRefRecNofHits / (Double_t) nofRefAccNofHits;
558 DrawH1({fHM->H1("fhNofHitsElAcc"), fHM->H1("fhNofHitsElRec")},
559 {"Acc (" + to_string(nofAccNofHits) + ")", "Rec" + to_string(nofRecNofHits) + ")"});
560
561 fHM->CreateCanvas("rich_tb_reco_eff_mom", "rich_tb_reco_eff_mom", 800, 800);
562 TH1D* hEffMom = Cbm::DivideH1(fHM->H1("fhMomElRec"), fHM->H1("fhMomElAcc"));
563 TH1D* hRefEffMom = Cbm::DivideH1(fHM->H1("fhMomRefElRec"), fHM->H1("fhMomRefElAcc"));
564 DrawH1({hEffMom, hRefEffMom},
565 {"Electrons (" + Cbm::NumberToString(effMom, 2) + "%)",
566 "Reference electrons (" + Cbm::NumberToString(effRefMom, 2) + "%)"},
567 kLinear, kLinear, true, 0.5, 0.85, 0.99, 0.99);
568 hEffMom->SetMinimum(0.);
569 hEffMom->SetMaximum(105.);
570
571 fHM->CreateCanvas("rich_tb_reco_eff_nofHits", "rich_tb_reco_eff_nofHits", 800, 800);
572 TH1D* hEffHits = Cbm::DivideH1(fHM->H1("fhNofHitsElRec"), fHM->H1("fhNofHitsElAcc"));
573 TH1D* hRefEffHits = Cbm::DivideH1(fHM->H1("fhNofHitsRefElRec"), fHM->H1("fhNofHitsRefElAcc"));
574 DrawH1({hEffHits, hRefEffHits},
575 {"Electron (" + Cbm::NumberToString(effNofHits, 2) + "%)",
576 "ElectronReference (" + Cbm::NumberToString(effRefNofHits, 2) + "%)"},
577 kLinear, kLinear, true, 0.5, 0.9, 0.99, 0.99);
578 hEffHits->SetMinimum(0.);
579 hEffHits->SetMaximum(105.);
580
581 fHM->CreateCanvas("rich_tb_reco_eff_eventMult", "rich_tb_reco_eff_eventMult", 800, 800);
582 TH1D* hEffEventMult = Cbm::DivideH1(fHM->H1("fhEventMultElRec"), fHM->H1("fhEventMultElAcc"));
583 DrawH1({hEffEventMult}, {"Electron (" + Cbm::NumberToString(effNofHits, 2) + "%)"}, kLinear, kLinear, true, 0.5,
584 0.9, 0.99, 0.99);
585 hEffHits->SetMinimum(0.);
586 hEffHits->SetMaximum(105.);
587 }
588}
589
590void CbmRichRecoTbQa::DrawTimeLog(const string& hMainName, Int_t nofLogEvents, bool withNoise)
591{
592 Double_t max = std::numeric_limits<Double_t>::min();
593 Int_t startInd = (withNoise) ? -1 : 0;
594 for (int iEv = startInd; iEv < nofLogEvents; iEv++) {
595 string hName = hMainName + "_" + to_string(iEv);
596 string option = ((iEv == 0 && !withNoise) || (iEv == -1 && withNoise)) ? "" : "same";
597 int ind = iEv % 6;
598 int color = (ind == 0) ? kBlue
599 : (ind == 1) ? kRed + 1
600 : (ind == 2) ? kGreen + 3
601 : (ind == 3) ? kMagenta + 2
602 : (ind == 4) ? kYellow + 2
603 : kOrange + 8;
604 if (ind == -1) color = kAzure + 6;
605 max = std::max(max, fHM->H1(hName)->GetMaximum());
606 DrawH1(fHM->H1(hName), kLinear, kLog, option, color);
607 }
608 fHM->H1(hMainName + "_" + to_string(startInd))->SetMaximum(max * 1.10);
609 gPad->SetLeftMargin(0.07);
610 gPad->SetRightMargin(0.05);
611}
612
613
614void CbmRichRecoTbQa::Finish()
615{
616 ProcessMc();
617
618 TDirectory* oldir = gDirectory;
619 TFile* outFile = FairRootManager::Instance()->GetOutFile();
620 if (outFile != NULL) {
621 outFile->mkdir(GetName());
622 outFile->cd(GetName());
623 fHM->WriteToFile();
624 }
625
626 DrawHist();
627 fHM->SaveCanvasToImage(fOutputDir);
628
629 gDirectory->cd(oldir->GetPath());
630}
631
632
633ClassImp(CbmRichRecoTbQa)
ClassImp
ClassImp(CbmConverterManager)
ECbmModuleId::kRich
@ kRich
Ring-Imaging Cherenkov Detector.
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
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.
fDigiMan
CbmDigiManager * fDigiMan
Definition CbmTofAnaTestbeam.cxx:91
CbmDigiManager::GetNofDigis
static Int_t GetNofDigis(ECbmModuleId systemId)
Definition CbmDigiManager.cxx:67
CbmDigiManager::Init
InitStatus Init()
Initialisation.
Definition CbmDigiManager.cxx:77
CbmDigiManager::Get
const Digi * Get(Int_t index) const
Get a digi object.
Definition CbmDigiManager.h:57
CbmDigiManager::Instance
static CbmDigiManager * Instance()
Static instance.
Definition CbmDigiManager.h:126
CbmDigiManager::GetMatch
const CbmMatch * GetMatch(ECbmModuleId systemId, UInt_t index) const
Get a match object.
Definition CbmDigiManager.cxx:59
CbmHistManager
Histogram manager.
Definition CbmHistManager.h:48
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::SaveCanvasToImage
void SaveCanvasToImage(const std::string &outputDir, const std::string &options="png,eps")
Save all stored canvases to images.
Definition core/base/CbmHistManager.cxx:310
CbmHistManager::WriteToFile
void WriteToFile()
Write all objects to current opened file.
Definition core/base/CbmHistManager.cxx:116
CbmHistManager::Create1
void Create1(const std::string &name, const std::string &title, Int_t nofBins, Double_t minBin, Double_t maxBin)
Helper function for creation of 1-dimensional histograms and profiles. Template argument is a real ob...
Definition CbmHistManager.h:88
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition CbmHistManager.h:207
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition CbmMCDataArray.h:51
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition CbmMCDataArray.cxx:137
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition CbmMCDataManager.h:31
CbmMCDataManager::GetObject
CbmMCDataObject * GetObject(const char *name)
Definition CbmMCDataManager.cxx:147
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition CbmMCDataManager.cxx:115
CbmMCEventList
Container class for MC events with number, file and start time.
Definition CbmMCEventList.h:43
CbmMCEventList::GetNofEvents
std::size_t GetNofEvents() const
Number of events in the list @value Number of events.
Definition CbmMCEventList.h:116
CbmMCEventList::GetEventTime
double GetEventTime(uint32_t event, uint32_t file)
Event start time.
Definition CbmMCEventList.cxx:89
CbmMCTrack::GetP
double GetP() const
Definition CbmMCTrack.h:98
CbmMCTrack::GetGeantProcessId
uint32_t GetGeantProcessId() const
Definition CbmMCTrack.h:67
CbmMCTrack::GetMotherId
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
CbmMCTrack::GetPdgCode
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
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::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
CbmRichDigi::GetTime
double GetTime() const
Definition CbmRichDigi.h:72
CbmRichRecoTbQa::IsCherenkovPhotonFromPrimaryElectron
Bool_t IsCherenkovPhotonFromPrimaryElectron(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition CbmRichRecoTbQa.cxx:381
CbmRichRecoTbQa::IsMcPrimaryElectron
Bool_t IsMcPrimaryElectron(const CbmMCTrack *mctrack)
Definition CbmRichRecoTbQa.cxx:413
CbmRichRecoTbQa::fRichRings
TClonesArray * fRichRings
Definition CbmRichRecoTbQa.h:116
CbmRichRecoTbQa::fRichRingMatches
TClonesArray * fRichRingMatches
Definition CbmRichRecoTbQa.h:117
CbmRichRecoTbQa::fStsPoints
CbmMCDataArray * fStsPoints
Definition CbmRichRecoTbQa.h:113
CbmRichRecoTbQa::InitHistograms
void InitHistograms()
Initialize histograms.
Definition CbmRichRecoTbQa.cxx:116
CbmRichRecoTbQa::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition CbmRichRecoTbQa.cxx:71
CbmRichRecoTbQa::fRichPoints
CbmMCDataArray * fRichPoints
Definition CbmRichRecoTbQa.h:112
CbmRichRecoTbQa::DrawTimeLog
void DrawTimeLog(const string &hMainName, Int_t nofLogEvents, bool withNoise=false)
Definition CbmRichRecoTbQa.cxx:590
CbmRichRecoTbQa::GetNofPrimaryMcTracks
Int_t GetNofPrimaryMcTracks(Int_t iEv)
Definition CbmRichRecoTbQa.cxx:300
CbmRichRecoTbQa::IsCherenkovPhoton
Bool_t IsCherenkovPhoton(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition CbmRichRecoTbQa.cxx:372
CbmRichRecoTbQa::CbmRichRecoTbQa
CbmRichRecoTbQa()
Standard constructor.
Definition CbmRichRecoTbQa.cxx:52
CbmRichRecoTbQa::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition CbmRichRecoTbQa.cxx:178
CbmRichRecoTbQa::fHM
CbmHistManager * fHM
Definition CbmRichRecoTbQa.h:104
CbmRichRecoTbQa::IsCherenkovPhotonFromPion
Bool_t IsCherenkovPhotonFromPion(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition CbmRichRecoTbQa.cxx:429
CbmRichRecoTbQa::fRichHits
TClonesArray * fRichHits
Definition CbmRichRecoTbQa.h:115
CbmRichRecoTbQa::fDigiMan
CbmDigiManager * fDigiMan
Definition CbmRichRecoTbQa.h:114
CbmRichRecoTbQa::fMCTracks
CbmMCDataArray * fMCTracks
Definition CbmRichRecoTbQa.h:111
CbmRichRecoTbQa::fRecRings
vector< CbmLink > fRecRings
Definition CbmRichRecoTbQa.h:120
CbmRichRecoTbQa::IsCherenkovPhotonFromSecondaryElectron
Bool_t IsCherenkovPhotonFromSecondaryElectron(const CbmRichPoint *point, Int_t fileId, Int_t eventId)
Definition CbmRichRecoTbQa.cxx:395
CbmRichRecoTbQa::fEventList
CbmMCEventList * fEventList
Definition CbmRichRecoTbQa.h:118
CbmRichRecoTbQa::IsMcPion
Bool_t IsMcPion(const CbmMCTrack *mctrack)
Definition CbmRichRecoTbQa.cxx:421
CbmRichRecoTbQa::Finish
virtual void Finish()
Inherited from FairTask.
Definition CbmRichRecoTbQa.cxx:614
CbmRichRing::GetTime
double GetTime() const
Definition CbmRichRing.h:99
CbmTrackMatchNew::GetTrueOverAllHitsRatio
double GetTrueOverAllHitsRatio() const
Definition CbmTrackMatchNew.h:39
Cbm::DivideH1
TH1D * DivideH1(TH1 *h1, TH1 *h2, const string &histName, double scale, const string &titleYaxis)
Definition CbmUtils.cxx:84
Cbm::NumberToString
std::string NumberToString(const T &value, int precision=1)
Definition CbmUtils.h:34
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21