BmnRoot 26.10.0-19-g6979a6c6
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
CbmStsFindTracksQa.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- CbmStsFindTracksQa source file -----
3// ----- Created 11/01/06 by V. Friese -----
4// -------------------------------------------------------------------------
5
6// Includes from sts
7#include "CbmStsFindTracksQa.h"
8
9#include "CbmGeoStsPar.h"
10#include "CbmStsHit.h"
11#include "CbmStsPoint.h"
12#include "CbmStsTrack.h"
13#include "CbmTrackMatch.h"
14
15// Includes from base
16#include "CbmMCTrack.h"
17#include "FairGeoNode.h"
18#include "FairGeoPassivePar.h"
19#include "FairGeoVector.h"
20#include "FairRootManager.h"
21#include "FairRunAna.h"
22#include "FairRuntimeDb.h"
23
24// Includes from ROOT
25#include "TCanvas.h"
26#include "TClonesArray.h"
27#include "TFile.h"
28#include "TH1F.h"
29#include "TList.h"
30#include "TVector3.h"
31
32// Includes from C++
33#include <iomanip>
34#include <iostream>
35using namespace std;
36
37// ----- Default constructor -------------------------------------------
38CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t iVerbose)
39 : FairTask("STSFindTracksQA", iVerbose)
40 , fHitMap()
41 , fMatchMap()
42 , fQualiMap()
43 , fMCTracks(NULL)
44 , fStsPoints(NULL)
45 , fStsHits(NULL)
46 , fStsTracks(NULL)
47 , fMatches(NULL)
48 , fPassGeo(NULL)
49 , fStsGeo(NULL)
50 , fTargetPos(0., 0., 0.)
51 , fNStations(0)
52 , fMinHits(4)
53 , fQuota(0.7)
54 , fhMomAccAll(new TH1F())
55 , fhMomRecAll(new TH1F())
56 , fhMomEffAll(new TH1F())
57 , fhMomAccPrim(new TH1F())
58 , fhMomRecPrim(new TH1F())
59 , fhMomEffPrim(new TH1F())
60 , fhMomAccSec(new TH1F())
61 , fhMomRecSec(new TH1F())
62 , fhMomEffSec(new TH1F())
63 , fhNpAccAll(new TH1F())
64 , fhNpRecAll(new TH1F())
65 , fhNpEffAll(new TH1F())
66 , fhNpAccPrim(new TH1F())
67 , fhNpRecPrim(new TH1F())
68 , fhNpEffPrim(new TH1F())
69 , fhNpAccSec(new TH1F())
70 , fhNpRecSec(new TH1F())
71 , fhNpEffSec(new TH1F())
72 , fhZAccSec(new TH1F())
73 , fhZRecSec(new TH1F())
74 , fhZEffSec(new TH1F())
75 , fhNhClones(new TH1F())
76 , fhNhGhosts(new TH1F())
77 , fHistoList(new TList())
78 , fNAccAll(0)
79 , fNAccPrim(0)
80 , fNAccRef(0)
81 , fNAccSec(0)
82 , fNRecAll(0)
83 , fNRecPrim(0)
84 , fNRecRef(0)
85 , fNRecSec(0)
86 , fNGhosts(0)
87 , fNClones(0)
88 , fNEvents(0)
89 , fNEventsFailed(0)
90 , fTime(0.)
91 , fTimer()
92{}
93
94// -------------------------------------------------------------------------
95
96// ----- Standard constructor ------------------------------------------
97CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t minHits, Double_t quota, Int_t iVerbose)
98 : FairTask("STSFindTracksQA", iVerbose)
99 , fHitMap()
100 , fMatchMap()
101 , fQualiMap()
102 , fMCTracks(NULL)
103 , fStsPoints(NULL)
104 , fStsHits(NULL)
105 , fStsTracks(NULL)
106 , fMatches(NULL)
107 , fPassGeo(NULL)
108 , fStsGeo(NULL)
109 , fTargetPos(0., 0., 0.)
110 , fNStations(0)
111 , fMinHits(minHits)
112 , fQuota(quota)
113 , fhMomAccAll(new TH1F())
114 , fhMomRecAll(new TH1F())
115 , fhMomEffAll(new TH1F())
116 , fhMomAccPrim(new TH1F())
117 , fhMomRecPrim(new TH1F())
118 , fhMomEffPrim(new TH1F())
119 , fhMomAccSec(new TH1F())
120 , fhMomRecSec(new TH1F())
121 , fhMomEffSec(new TH1F())
122 , fhNpAccAll(new TH1F())
123 , fhNpRecAll(new TH1F())
124 , fhNpEffAll(new TH1F())
125 , fhNpAccPrim(new TH1F())
126 , fhNpRecPrim(new TH1F())
127 , fhNpEffPrim(new TH1F())
128 , fhNpAccSec(new TH1F())
129 , fhNpRecSec(new TH1F())
130 , fhNpEffSec(new TH1F())
131 , fhZAccSec(new TH1F())
132 , fhZRecSec(new TH1F())
133 , fhZEffSec(new TH1F())
134 , fhNhClones(new TH1F())
135 , fhNhGhosts(new TH1F())
136 , fHistoList(new TList())
137 , fNAccAll(0)
138 , fNAccPrim(0)
139 , fNAccRef(0)
140 , fNAccSec(0)
141 , fNRecAll(0)
142 , fNRecPrim(0)
143 , fNRecRef(0)
144 , fNRecSec(0)
145 , fNGhosts(0)
146 , fNClones(0)
147 , fNEvents(0)
148 , fNEventsFailed(0)
149 , fTime(0.)
150 , fTimer()
151{}
152// -------------------------------------------------------------------------
153
154// ----- Destructor ----------------------------------------------------
155CbmStsFindTracksQa::~CbmStsFindTracksQa()
156{
157
158 fHistoList->Delete();
159 delete fHistoList;
160}
161// -------------------------------------------------------------------------
162
163// ----- Public method SetParContainers --------------------------------
164void CbmStsFindTracksQa::SetParContainers()
165{
166
167 // Get Run
168 FairRunAna* run = FairRunAna::Instance();
169 if (!run) {
170 cout << "-E- " << GetName() << "::SetParContainers: No FairRunAna!" << endl;
171 return;
172 }
173
174 // Get Runtime Database
175 FairRuntimeDb* runDb = run->GetRuntimeDb();
176 if (!run) {
177 cout << "-E- " << GetName() << "::SetParContainers: No runtime database!" << endl;
178 return;
179 }
180
181 // Get passive geometry parameters
182 fPassGeo = (FairGeoPassivePar*)runDb->getContainer("FairGeoPassivePar");
183 if (!fPassGeo) {
184 cout << "-E- " << GetName() << "::SetParContainers: "
185 << "No passive geometry parameters!" << endl;
186 return;
187 }
188
189 // Get STS geometry parameters
190 fStsGeo = (CbmGeoStsPar*)runDb->getContainer("CbmGeoStsPar");
191 if (!fStsGeo) {
192 cout << "-E- " << GetName() << "::SetParContainers: "
193 << "No STS geometry parameters!" << endl;
194 return;
195 }
196}
197// -------------------------------------------------------------------------
198
199// ----- Public method Init --------------------------------------------
200InitStatus CbmStsFindTracksQa::Init()
201{
202
203 cout << "===========================================================" << endl;
204 ;
205 cout << GetName() << ": Initialising..." << endl;
206
207 // Get RootManager
208 FairRootManager* ioman = FairRootManager::Instance();
209 if (!ioman) {
210 cout << "-E- " << GetName() << "::Init: "
211 << "RootManager not instantised!" << endl;
212 return kFATAL;
213 }
214
215 // Get MCTrack array
216 fMCTracks = (TClonesArray*)ioman->GetObject("MCTrack");
217 if (!fMCTracks) {
218 cout << "-E- " << GetName() << "::Init: No MCTrack array!" << endl;
219 return kFATAL;
220 }
221
222 // Get StsPoint array
223 fStsPoints = (TClonesArray*)ioman->GetObject("StsPoint");
224 if (!fStsPoints) {
225 // cout << "-E- " << GetName() << "::Init: No StsPoint array!" << endl;
226 return kFATAL;
227 }
228
229 // Get StsHit array
230 fStsHits = (TClonesArray*)ioman->GetObject("StsHit");
231 if (!fStsHits) {
232 cout << "-E- " << GetName() << "::Init: No StsHit array!" << endl;
233 return kFATAL;
234 }
235
236 // Get StsTrack array
237 fStsTracks = (TClonesArray*)ioman->GetObject("StsTrack");
238 if (!fStsTracks) {
239 cout << "-E- " << GetName() << "::Init: No StsTrack array!" << endl;
240 return kERROR;
241 }
242
243 // Get StsTrackMatch array
244 fMatches = (TClonesArray*)ioman->GetObject("StsTrackMatch");
245 if (!fMatches) {
246 cout << "-E- " << GetName() << "::Init: No StsTrackMatch array!" << endl;
247 return kERROR;
248 }
249
250 // Get the geometry of target and STS
251 InitStatus geoStatus = GetGeometry();
252 if (geoStatus != kSUCCESS) {
253 cout << "-E- " << GetName() << "::Init: Error in reading geometry!" << endl;
254 return geoStatus;
255 }
256
257 // Create histograms
258 CreateHistos();
259 Reset();
260
261 // Output
262 cout << " Minimum number of STS hits : " << fMinHits << endl;
263 cout << " Matching quota : " << fQuota << endl;
264 cout << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") "
265 << endl;
266 cout << " Number of STS stations : " << fNStations << endl;
267 if (fActive)
268 cout << " ***** Task is ACTIVE *****" << endl;
269 cout << "===========================================================" << endl << endl;
270
271 return geoStatus;
272}
273// -------------------------------------------------------------------------
274
275// ----- Public method ReInit ------------------------------------------
276InitStatus CbmStsFindTracksQa::ReInit()
277{
278
279 cout << "===========================================================" << endl;
280 ;
281 cout << GetName() << ": Reinitialising..." << endl;
282
283 // Get the geometry of target and STS
284 InitStatus geoStatus = GetGeometry();
285 if (geoStatus != kSUCCESS) {
286 cout << "-E- " << GetName() << "::ReInit: Error in reading geometry!" << endl;
287 return geoStatus;
288 }
289
290 // Output
291 cout << " Target position ( " << fTargetPos.X() << ", " << fTargetPos.Y() << ", " << fTargetPos.Z() << ") "
292 << endl;
293 cout << " Number of STS stations : " << fNStations << endl;
294 if (fActive)
295 cout << " ***** Task is ACTIVE *****" << endl;
296 cout << "===========================================================" << endl << endl;
297
298 return geoStatus;
299}
300// -------------------------------------------------------------------------
301
302// ----- Public method Exec --------------------------------------------
303void CbmStsFindTracksQa::Exec(Option_t* opt)
304{
305
306 // Timer
307 fTimer.Start();
308
309 // Eventwise counters
310 Int_t nRec = 0;
311 Int_t nGhosts = 0;
312 Int_t nClones = 0;
313 Int_t nAll = 0;
314 Int_t nAcc = 0;
315 Int_t nRecAll = 0;
316 Int_t nPrim = 0;
317 Int_t nRecPrim = 0;
318 Int_t nRef = 0;
319 Int_t nRecRef = 0;
320 Int_t nSec = 0;
321 Int_t nRecSec = 0;
322 TVector3 momentum;
323 TVector3 vertex;
324
325 // Fill hit and track maps
326 FillHitMap();
327 FillMatchMap(nRec, nGhosts, nClones);
328
329 // Loop over MCTracks
330 Int_t nMC = fMCTracks->GetEntriesFast();
331 for (Int_t iMC = 0; iMC < nMC; iMC++) {
332 CbmMCTrack* mcTrack = (CbmMCTrack*)fMCTracks->At(iMC);
333 if (!mcTrack) {
334 cout << "-E- " << GetName() << "::Exec: "
335 << "No MCTrack at index " << iMC << endl;
336 Fatal("Exec", "No MCTrack in array");
337 }
338
339 // Check reconstructability; continue only for reconstructable tracks
340 nAll++;
341 if (fHitMap.find(iMC) == fHitMap.end())
342 continue; // No hits
343 Int_t nHits = fHitMap[iMC];
344 if (nHits < fMinHits)
345 continue; // Too few hits
346 nAcc++;
347
348 // Check origin of MCTrack
349 mcTrack->GetStartVertex(vertex);
350 Bool_t isPrim = kFALSE;
351 if ((vertex - fTargetPos).Mag() < 1.) {
352 isPrim = kTRUE;
353 nPrim++;
354 } else
355 nSec++;
356
357 // Get momentum
358 mcTrack->GetMomentum(momentum);
359 Double_t mom = momentum.Mag();
360 Bool_t isRef = kFALSE;
361 if (mom > 1. && isPrim) {
362 isRef = kTRUE;
363 nRef++;
364 }
365
366 // Fill histograms for reconstructable tracks
367 fhMomAccAll->Fill(mom);
368 fhNpAccAll->Fill(Double_t(nHits));
369 if (isPrim) {
370 fhMomAccPrim->Fill(mom);
371 fhNpAccPrim->Fill(Double_t(nHits));
372 } else {
373 fhMomAccSec->Fill(mom);
374 fhNpAccSec->Fill(Double_t(nHits));
375 fhZAccSec->Fill(vertex.Z());
376 }
377
378 // Get matched StsTrack
379 Int_t iRec = -1;
380 Double_t quali = 0.;
381 // Bool_t isRec = kFALSE;
382 if (fMatchMap.find(iMC) != fMatchMap.end()) {
383 iRec = fMatchMap[iMC];
384 // isRec = kTRUE;
385 CbmStsTrack* stsTrack = (CbmStsTrack*)fStsTracks->At(iRec);
386 if (!stsTrack) {
387 cout << "-E- " << GetName() << "::Exec: "
388 << "No StsTrack for matched MCTrack " << iMC << endl;
389 Fatal("Exec", "No StsTrack for matched MCTrack");
390 }
391 quali = fQualiMap[iMC];
392 if (quali < fQuota) {
393 cout << "-E- " << GetName() << "::Exec: "
394 << "Matched StsTrack " << iRec << " is below matching "
395 << "criterion ( " << quali << ")" << endl;
396 Fatal("Exec", "Match below matching quota");
397 }
398 CbmTrackMatch* match = (CbmTrackMatch*)fMatches->At(iRec);
399 if (!match) {
400 cout << "-E- " << GetName() << "::Exec: "
401 << "No StsTrackMatch for matched MCTrack " << iMC << endl;
402 Fatal("Exec", "No StsTrackMatch for matched MCTrack");
403 }
404 Int_t nTrue = match->GetNofTrueHits();
405 Int_t nWrong = match->GetNofWrongHits();
406 Int_t nFake = match->GetNofFakeHits();
407 Int_t nAllHits = stsTrack->GetNStsHits();
408 if (nTrue + nWrong + nFake != nAllHits) {
409 cout << "True " << nTrue << " wrong " << nWrong << " Fake " << nFake << " Hits " << nAllHits << endl;
410 Fatal("Exec", "Wrong number of hits");
411 }
412
413 // Verbose output
414 if (fVerbose > 1)
415 cout << "-I- " << GetName() << ": "
416 << "MCTrack " << iMC << ", hits " << nHits << ", StsTrack " << iRec << ", hits " << nAllHits
417 << ", true hits " << nTrue << endl;
418
419 // Fill histograms for reconstructed tracks
420 nRecAll++;
421 fhMomRecAll->Fill(mom);
422 fhNpRecAll->Fill(Double_t(nAllHits));
423 if (isPrim) {
424 nRecPrim++;
425 fhMomRecPrim->Fill(mom);
426 fhNpRecPrim->Fill(Double_t(nAllHits));
427 if (isRef)
428 nRecRef++;
429 } else {
430 nRecSec++;
431 fhMomRecSec->Fill(mom);
432 fhNpRecSec->Fill(Double_t(nAllHits));
433 fhZRecSec->Fill(vertex.Z());
434 }
435
436 } // Match found in map?
437
438 } // Loop over MCTracks
439
440 // Calculate efficiencies
441 Double_t effAll = Double_t(nRecAll) / Double_t(nAcc);
442 Double_t effPrim = Double_t(nRecPrim) / Double_t(nPrim);
443 Double_t effRef = Double_t(nRecRef) / Double_t(nRef);
444 Double_t effSec = Double_t(nRecSec) / Double_t(nSec);
445
446 // Event summary
447 if (fVerbose > 0) {
448 cout << "---------- StsFindTracksQa : Event summary ------------" << endl;
449 cout << "MCTracks : " << nAll << ", reconstructable: " << nAcc << ", reconstructed: " << nRecAll << endl;
450 cout << "Vertex : reconstructable: " << nPrim << ", reconstructed: " << nRecPrim << ", efficiency "
451 << effPrim * 100. << "%" << endl;
452 cout << "Reference : reconstructable: " << nRef << ", reconstructed: " << nRecRef << ", efficiency "
453 << effRef * 100. << "%" << endl;
454 cout << "Non-vertex : reconstructable: " << nSec << ", reconstructed: " << nRecSec << ", efficiency "
455 << effSec * 100. << "%" << endl;
456 cout << "STSTracks " << nRec << ", ghosts " << nGhosts << ", clones " << nClones << endl;
457 cout << "-----------------------------------------------------------" << endl;
458 cout << endl;
459 } else
460 cout << "+ " << setw(15) << left << fName << ": " << setprecision(4) << setw(8) << fixed << right
461 << fTimer.RealTime() << " s, efficiency all " << effAll * 100. << " %, vertex " << effPrim * 100.
462 << " %, reference " << effRef * 100. << " %" << endl;
463
464 // Increase counters
465 fNAccAll += nAcc;
466 fNAccPrim += nPrim;
467 fNAccRef += nRef;
468 fNAccSec += nSec;
469 fNRecAll += nRecAll;
470 fNRecPrim += nRecPrim;
471 fNRecRef += nRecRef;
472 fNRecSec += nRecSec;
473 fNGhosts += nGhosts;
474 fNClones += nClones;
475 fNEvents++;
476 fTime += fTimer.RealTime();
477}
478// -------------------------------------------------------------------------
479
480// ----- Private method Finish -----------------------------------------
481void CbmStsFindTracksQa::Finish()
482{
483
484 // Divide histograms for efficiency calculation
485 DivideHistos(fhMomRecAll, fhMomAccAll, fhMomEffAll);
486 DivideHistos(fhMomRecPrim, fhMomAccPrim, fhMomEffPrim);
487 DivideHistos(fhMomRecSec, fhMomAccSec, fhMomEffSec);
488 DivideHistos(fhNpRecAll, fhNpAccAll, fhNpEffAll);
489 DivideHistos(fhNpRecPrim, fhNpAccPrim, fhNpEffPrim);
490 DivideHistos(fhNpRecSec, fhNpAccSec, fhNpEffSec);
491 DivideHistos(fhZRecSec, fhZAccSec, fhZEffSec);
492
493 // Normalise histos for clones and ghosts to one event
494 if (fNEvents) {
495 fhNhClones->Scale(1. / Double_t(fNEvents));
496 fhNhGhosts->Scale(1. / Double_t(fNEvents));
497 }
498
499 // Calculate integrated efficiencies and rates
500 Double_t effAll = Double_t(fNRecAll) / Double_t(fNAccAll);
501 Double_t effPrim = Double_t(fNRecPrim) / Double_t(fNAccPrim);
502 Double_t effRef = Double_t(fNRecRef) / Double_t(fNAccRef);
503 Double_t effSec = Double_t(fNRecSec) / Double_t(fNAccSec);
504 Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNEvents);
505 Double_t rateClones = Double_t(fNClones) / Double_t(fNEvents);
506
507 // Run summary to screen
508 cout << endl;
509 cout << "============================================================" << endl;
510 cout << "===== " << fName << ": Run summary " << endl;
511 cout << "===== " << endl;
512 cout << "===== Good events : " << setw(6) << fNEvents << endl;
513 cout << "===== Failed events : " << setw(6) << fNEventsFailed << endl;
514 cout << "===== Average time : " << setprecision(4) << setw(8) << right << fTime / Double_t(fNEvents) << " s"
515 << endl;
516 cout << "===== " << endl;
517 cout << "===== Efficiency all tracks : " << effAll * 100 << " % (" << fNRecAll << "/" << fNAccAll << ")"
518 << endl;
519 cout << "===== Efficiency vertex tracks : " << effPrim * 100 << " % (" << fNRecPrim << "/" << fNAccPrim << ")"
520 << endl;
521 cout << "===== Efficiency reference tracks : " << effRef * 100 << " % (" << fNRecRef << "/" << fNAccRef << ")"
522 << endl;
523 cout << "===== Efficiency secondary tracks : " << effSec * 100 << " % (" << fNRecSec << "/" << fNAccSec << ")"
524 << endl;
525 cout << "===== Ghost rate " << rateGhosts << " per event" << endl;
526 cout << "===== Clone rate " << rateClones << " per event" << endl;
527 cout << "============================================================" << endl;
528
529 // Write histos to output
530 gDirectory->mkdir("STSFindTracksQA");
531 gDirectory->cd("STSFindTracksQA");
532 TIter next(fHistoList);
533 while (TH1* histo = ((TH1*)next()))
534 histo->Write();
535 gDirectory->cd("..");
536}
537// -------------------------------------------------------------------------
538
539// ----- Private method GetGeometry ------------------------------------
540InitStatus CbmStsFindTracksQa::GetGeometry()
541{
542
543 // Get target geometry
544 if (!fPassGeo) {
545 cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!" << endl;
546 fTargetPos.SetXYZ(0., 0., 0.);
547 return kERROR;
548 }
549 TObjArray* passNodes = fPassGeo->GetGeoPassiveNodes();
550 if (!passNodes) {
551 cout << "-W- " << GetName() << "::GetGeometry: No passive node array" << endl;
552 fTargetPos.SetXYZ(0., 0., 0.);
553 return kERROR;
554 }
555 FairGeoNode* target = (FairGeoNode*)passNodes->FindObject("targ");
556 if (!target) {
557 cout << "-E- " << GetName() << "::GetGeometry: No target node" << endl;
558 fTargetPos.SetXYZ(0., 0., 0.);
559 return kERROR;
560 }
561 FairGeoVector targetPos = target->getLabTransform()->getTranslation();
562 FairGeoVector centerPos = target->getCenterPosition().getTranslation();
563 Double_t targetX = targetPos.X() + centerPos.X();
564 Double_t targetY = targetPos.Y() + centerPos.Y();
565 Double_t targetZ = targetPos.Z() + centerPos.Z();
566 fTargetPos.SetXYZ(targetX, targetY, targetZ);
567
568 // Get STS geometry
569 if (!fStsGeo) {
570 cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!" << endl;
571 fNStations = 0;
572 return kERROR;
573 }
574 TObjArray* stsNodes = fStsGeo->GetGeoSensitiveNodes();
575 if (!stsNodes) {
576 cout << "-E- " << GetName() << "::GetGeometry: No STS node array" << endl;
577 fNStations = 0;
578 return kERROR;
579 }
580 fNStations = stsNodes->GetEntries();
581
582 return kSUCCESS;
583}
584// -------------------------------------------------------------------------
585
586// ----- Private method CreateHistos -----------------------------------
587void CbmStsFindTracksQa::CreateHistos()
588{
589
590 // Histogram list
591 fHistoList = new TList();
592
593 // Momentum distributions
594 Double_t minMom = 0.;
595 Double_t maxMom = 10.;
596 Int_t nBinsMom = 40;
597 fhMomAccAll = new TH1F("hMomAccAll", "all reconstructable tracks", nBinsMom, minMom, maxMom);
598 fhMomRecAll = new TH1F("hMomRecAll", "all reconstructed tracks", nBinsMom, minMom, maxMom);
599 fhMomEffAll = new TH1F("hMomEffAll", "efficiency all tracks", nBinsMom, minMom, maxMom);
600 fhMomAccPrim = new TH1F("hMomAccPrim", "reconstructable vertex tracks", nBinsMom, minMom, maxMom);
601 fhMomRecPrim = new TH1F("hMomRecPrim", "reconstructed vertex tracks", nBinsMom, minMom, maxMom);
602 fhMomEffPrim = new TH1F("hMomEffPrim", "efficiency vertex tracks", nBinsMom, minMom, maxMom);
603 fhMomAccSec = new TH1F("hMomAccSec", "reconstructable non-vertex tracks", nBinsMom, minMom, maxMom);
604 fhMomRecSec = new TH1F("hMomRecSec", "reconstructed non-vertex tracks", nBinsMom, minMom, maxMom);
605 fhMomEffSec = new TH1F("hMomEffSec", "efficiency non-vertex tracks", nBinsMom, minMom, maxMom);
606 fHistoList->Add(fhMomAccAll);
607 fHistoList->Add(fhMomRecAll);
608 fHistoList->Add(fhMomEffAll);
609 fHistoList->Add(fhMomAccPrim);
610 fHistoList->Add(fhMomRecPrim);
611 fHistoList->Add(fhMomEffPrim);
612 fHistoList->Add(fhMomAccSec);
613 fHistoList->Add(fhMomRecSec);
614 fHistoList->Add(fhMomEffSec);
615
616 // Number-of-points distributions
617 Double_t minNp = -0.5;
618 Double_t maxNp = 15.5;
619 Int_t nBinsNp = 16;
620 fhNpAccAll = new TH1F("hNpAccAll", "all reconstructable tracks", nBinsNp, minNp, maxNp);
621 fhNpRecAll = new TH1F("hNpRecAll", "all reconstructed tracks", nBinsNp, minNp, maxNp);
622 fhNpEffAll = new TH1F("hNpEffAll", "efficiency all tracks", nBinsNp, minNp, maxNp);
623 fhNpAccPrim = new TH1F("hNpAccPrim", "reconstructable vertex tracks", nBinsNp, minNp, maxNp);
624 fhNpRecPrim = new TH1F("hNpRecPrim", "reconstructed vertex tracks", nBinsNp, minNp, maxNp);
625 fhNpEffPrim = new TH1F("hNpEffPrim", "efficiency vertex tracks", nBinsNp, minNp, maxNp);
626 fhNpAccSec = new TH1F("hNpAccSec", "reconstructable non-vertex tracks", nBinsNp, minNp, maxNp);
627 fhNpRecSec = new TH1F("hNpRecSec", "reconstructed non-vertex tracks", nBinsNp, minNp, maxNp);
628 fhNpEffSec = new TH1F("hNpEffSec", "efficiency non-vertex tracks", nBinsNp, minNp, maxNp);
629 fHistoList->Add(fhNpAccAll);
630 fHistoList->Add(fhNpRecAll);
631 fHistoList->Add(fhNpEffAll);
632 fHistoList->Add(fhNpAccPrim);
633 fHistoList->Add(fhNpRecPrim);
634 fHistoList->Add(fhNpEffPrim);
635 fHistoList->Add(fhNpAccSec);
636 fHistoList->Add(fhNpRecSec);
637 fHistoList->Add(fhNpEffSec);
638
639 // z(vertex) distributions
640 Double_t minZ = 0.;
641 Double_t maxZ = 50.;
642 Int_t nBinsZ = 50;
643 fhZAccSec = new TH1F("hZAccSec", "reconstructable non-vertex tracks", nBinsZ, minZ, maxZ);
644 fhZRecSec = new TH1F("hZRecSecl", "reconstructed non-vertex tracks", nBinsZ, minZ, maxZ);
645 fhZEffSec = new TH1F("hZEffRec", "efficiency non-vertex tracks", nBinsZ, minZ, maxZ);
646 fHistoList->Add(fhZAccSec);
647 fHistoList->Add(fhZRecSec);
648 fHistoList->Add(fhZEffSec);
649
650 // Number-of-hit distributions
651 fhNhClones = new TH1F("hNhClones", "number of hits for clones", nBinsNp, minNp, maxNp);
652 fhNhGhosts = new TH1F("hNhGhosts", "number of hits for ghosts", nBinsNp, minNp, maxNp);
653 fHistoList->Add(fhNhClones);
654 fHistoList->Add(fhNhGhosts);
655}
656// -------------------------------------------------------------------------
657
658// ----- Private method Reset ------------------------------------------
659void CbmStsFindTracksQa::Reset()
660{
661
662 TIter next(fHistoList);
663 while (TH1* histo = ((TH1*)next()))
664 histo->Reset();
665
666 fNAccAll = fNAccPrim = fNAccRef = fNAccSec = 0;
667 fNRecAll = fNRecPrim = fNRecRef = fNRecSec = 0;
668 fNGhosts = fNClones = fNEvents = 0;
669}
670// -------------------------------------------------------------------------
671
672// ----- Private method FillHitMap -------------------------------------
673void CbmStsFindTracksQa::FillHitMap()
674{
675 fHitMap.clear();
676 Int_t nHits = fStsHits->GetEntriesFast();
677 for (Int_t iHit = 0; iHit < nHits; iHit++) {
678 CbmStsHit* hit = (CbmStsHit*)fStsHits->At(iHit);
679 Int_t iMc = hit->GetRefIndex();
680 if (iMc < 0)
681 continue;
682 CbmStsPoint* stsPoint = (CbmStsPoint*)fStsPoints->At(iMc);
683 Int_t iTrack = stsPoint->GetTrackID();
684 fHitMap[iTrack]++;
685 }
686}
687// -------------------------------------------------------------------------
688
689// ------ Private method FillMatchMap ----------------------------------
690void CbmStsFindTracksQa::FillMatchMap(Int_t& nRec, Int_t& nGhosts, Int_t& nClones)
691{
692
693 // Clear matching maps
694 fMatchMap.clear();
695 fQualiMap.clear();
696
697 // Loop over StsTracks. Check matched MCtrack and fill maps.
698 nGhosts = 0;
699 nClones = 0;
700 nRec = fStsTracks->GetEntriesFast();
701 Int_t nMtc = fMatches->GetEntriesFast();
702 if (nMtc != nRec) {
703 cout << "-E- " << GetName() << "::Exec: Number of StsMatches (" << nMtc
704 << ") does not equal number of StsTracks (" << nRec << ")" << endl;
705 Fatal("Exec", "Inequal number of StsTrack and StsTrackMatch");
706 }
707 for (Int_t iRec = 0; iRec < nRec; iRec++) {
708
709 CbmStsTrack* stsTrack = (CbmStsTrack*)fStsTracks->At(iRec);
710 if (!stsTrack) {
711 cout << "-E- " << GetName() << "::Exec: "
712 << "No StsTrack at index " << iRec << endl;
713 Fatal("Exec", "No StsTrack in array");
714 }
715 Int_t nHits = stsTrack->GetNStsHits();
716
717 CbmTrackMatch* match = (CbmTrackMatch*)fMatches->At(iRec);
718 if (!match) {
719 cout << "-E- " << GetName() << "::Exec: "
720 << "No StsTrackMatch at index " << iRec << endl;
721 Fatal("Exec", "No StsTrackMatch in array");
722 }
723 Int_t nTrue = match->GetNofTrueHits();
724
725 Int_t iMC = match->GetMCTrackId();
726 if (iMC == -1) { // no common point with MC, really ghastly!
727 if (fVerbose > 1)
728 cout << "-I- " << GetName() << ":"
729 << "No MC match for StsTrack " << iRec << endl;
730 fhNhGhosts->Fill(nHits);
731 nGhosts++;
732 continue;
733 }
734
735 // Check matching criterion (quota)
736 Double_t quali = Double_t(nTrue) / Double_t(nHits);
737 if (quali >= fQuota) {
738
739 // No previous match for this MCTrack
740 if (fMatchMap.find(iMC) == fMatchMap.end()) {
741 fMatchMap[iMC] = iRec;
742 fQualiMap[iMC] = quali;
743 }
744
745 // Previous match; take the better one
746 else
747 {
748 if (fVerbose > 1)
749 cout << "-I- " << GetName() << ": "
750 << "MCTrack " << iMC << " doubly matched."
751 << "Current match " << iRec << ", previous match " << fMatchMap[iMC] << endl;
752 if (fQualiMap[iMC] < quali) {
753 CbmStsTrack* oldTrack = (CbmStsTrack*)fStsTracks->At(fMatchMap[iMC]);
754 fhNhClones->Fill(Double_t(oldTrack->GetNStsHits()));
755 fMatchMap[iMC] = iRec;
756 fQualiMap[iMC] = quali;
757 } else
758 fhNhClones->Fill(nHits);
759 nClones++;
760 }
761
762 }
763
764 // If not matched, it's a ghost
765 else
766 {
767 if (fVerbose > 1)
768 cout << "-I- " << GetName() << ":"
769 << "StsTrack " << iRec << " below matching criterion "
770 << "(" << quali << ")" << endl;
771 fhNhGhosts->Fill(nHits);
772 nGhosts++;
773 }
774
775 } // Loop over StsTracks
776}
777// -------------------------------------------------------------------------
778
779// ----- Private method DivideHistos -----------------------------------
780void CbmStsFindTracksQa::DivideHistos(TH1* histo1, TH1* histo2, TH1* histo3)
781{
782
783 if (!histo1 || !histo2 || !histo3) {
784 cout << "-E- " << GetName() << "::DivideHistos: "
785 << "NULL histogram pointer" << endl;
786 Fatal("DivideHistos", "Null histo pointer");
787 }
788
789 Int_t nBins = histo1->GetNbinsX();
790 if (histo2->GetNbinsX() != nBins || histo3->GetNbinsX() != nBins) {
791 cout << "-E- " << GetName() << "::DivideHistos: "
792 << "Different bin numbers in histos" << endl;
793 cout << histo1->GetName() << " " << histo1->GetNbinsX() << endl;
794 cout << histo2->GetName() << " " << histo2->GetNbinsX() << endl;
795 cout << histo3->GetName() << " " << histo3->GetNbinsX() << endl;
796 return;
797 }
798
799 Double_t c1, c2, c3, ce;
800 for (Int_t iBin = 0; iBin < nBins; iBin++) {
801 c1 = histo1->GetBinContent(iBin);
802 c2 = histo2->GetBinContent(iBin);
803 if (c2 != 0.) {
804 c3 = c1 / c2;
805 Double_t c4 = (c3 * (1. - c3) / c2);
806 if (c4 >= 0.) {
807 ce = TMath::Sqrt(c3 * (1. - c3) / c2);
808 } else {
809 ce = 0;
810 }
811 } else {
812 c3 = 0.;
813 ce = 0.;
814 }
815 histo3->SetBinContent(iBin, c3);
816 histo3->SetBinError(iBin, ce);
817 }
818}
819// -------------------------------------------------------------------------
Bool_t
bool Bool_t
Int_t
int Int_t
kFALSE
constexpr Bool_t kFALSE
Double_t
double Double_t
kTRUE
constexpr Bool_t kTRUE
Option_t
const char Option_t
gDirectory
#define gDirectory
Fatal
void Fatal(const char *location, const char *msgfmt,...)
CbmMCTrack
CbmMCTrack.h.
Definition CbmMCTrack.h:26
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:139
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:149
CbmStsFindTracksQa::fTimer
TStopwatch fTimer
Total real time used for good events.
Definition CbmStsFindTracksQa.h:136
CbmStsFindTracksQa::Finish
virtual void Finish()
Finish.
Definition CbmStsFindTracksQa.cxx:481
CbmStsFindTracksQa::CbmStsFindTracksQa
CbmStsFindTracksQa(Int_t iVerbose=1)
Default constructor.
Definition CbmStsFindTracksQa.cxx:38
CbmStsFindTracksQa::FillHitMap
void FillHitMap()
Fill a map from MCTrack index to number of corresponding StsHits.
Definition CbmStsFindTracksQa.cxx:673
CbmStsFindTracksQa::CreateHistos
void CreateHistos()
Create histograms.
Definition CbmStsFindTracksQa.cxx:587
CbmStsFindTracksQa::fStsTracks
TClonesArray * fStsTracks
Definition CbmStsFindTracksQa.h:101
CbmStsFindTracksQa::Init
virtual InitStatus Init()
Initialisation.
Definition CbmStsFindTracksQa.cxx:200
CbmStsFindTracksQa::fPassGeo
FairGeoPassivePar * fPassGeo
Geometry parameters.
Definition CbmStsFindTracksQa.h:105
CbmStsFindTracksQa::fNEventsFailed
Int_t fNEventsFailed
Number of events with success.
Definition CbmStsFindTracksQa.h:132
CbmStsFindTracksQa::fMatchMap
std::map< Int_t, Int_t > fMatchMap
Map from MCTrack index to matched StsTrack index.
Definition CbmStsFindTracksQa.h:92
CbmStsFindTracksQa::FillMatchMap
void FillMatchMap(Int_t &nRec, Int_t &nGhosts, Int_t &nClones)
Fill a map from MCTrack index to matched StsTrack index.
Definition CbmStsFindTracksQa.cxx:690
CbmStsFindTracksQa::~CbmStsFindTracksQa
virtual ~CbmStsFindTracksQa()
Destructor.
Definition CbmStsFindTracksQa.cxx:155
CbmStsFindTracksQa::Exec
virtual void Exec(Option_t *opt)
Execution.
Definition CbmStsFindTracksQa.cxx:303
CbmStsFindTracksQa::fMCTracks
TClonesArray * fMCTracks
Pointers to data arrays.
Definition CbmStsFindTracksQa.h:98
CbmStsFindTracksQa::fQualiMap
std::map< Int_t, Double_t > fQualiMap
Map from MCTrack index to percentage of matched hits.
Definition CbmStsFindTracksQa.h:95
CbmStsFindTracksQa::fhMomAccAll
TH1F * fhMomAccAll
Histograms.
Definition CbmStsFindTracksQa.h:115
CbmStsFindTracksQa::GetGeometry
InitStatus GetGeometry()
Read the geometry parameters.
Definition CbmStsFindTracksQa.cxx:540
CbmStsFindTracksQa::fMinHits
Int_t fMinHits
Task parameters.
Definition CbmStsFindTracksQa.h:111
CbmStsFindTracksQa::fStsPoints
TClonesArray * fStsPoints
Definition CbmStsFindTracksQa.h:99
CbmStsFindTracksQa::SetParContainers
virtual void SetParContainers()
Set parameter containers.
Definition CbmStsFindTracksQa.cxx:164
CbmStsFindTracksQa::fHitMap
std::map< Int_t, Int_t > fHitMap
Map from MCTrack index to number of attached StsHits.
Definition CbmStsFindTracksQa.h:89
CbmStsFindTracksQa::fHistoList
TList * fHistoList
List of histograms.
Definition CbmStsFindTracksQa.h:125
CbmStsFindTracksQa::fTime
Double_t fTime
Number of events with failure.
Definition CbmStsFindTracksQa.h:133
CbmStsFindTracksQa::Reset
void Reset()
Reset histograms and counters.
Definition CbmStsFindTracksQa.cxx:659
CbmStsFindTracksQa::DivideHistos
void DivideHistos(TH1 *histo1, TH1 *histo2, TH1 *histo3)
Divide histograms (reco/all) with correct error for the efficiency.
Definition CbmStsFindTracksQa.cxx:780
CbmStsFindTracksQa::ReInit
virtual InitStatus ReInit()
Reinitialisation.
Definition CbmStsFindTracksQa.cxx:276
CbmStsTrack::GetNStsHits
Int_t GetNStsHits() const
Accessors.
Definition CbmStsTrack.h:60
CbmTrackMatch
CbmStsTrackMatch.h.
Definition CbmTrackMatch.h:16
TH1::GetNbinsX
virtual Int_t GetNbinsX() const
TH1::SetBinError
virtual void SetBinError(Int_t bin, Double_t error)
TH1::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
TNamed::GetName
const char * GetName() const override
TObjArray::GetEntries
Int_t GetEntries() const override
TObjArray::FindObject
TObject * FindObject(const char *name) const override
TVector3::Z
Double_t Z() const
TVector3::Mag
Double_t Mag() const
TMath::Sqrt
Double_t Sqrt(Double_t x)