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