CbmRoot
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 class header
7 #include "CbmStsFindTracksQa.h"
8 
9 // Includes from C++
10 #include <cassert>
11 #include <iomanip>
12 
13 // Includes from ROOT
14 #include "TClonesArray.h"
15 #include "TGeoManager.h"
16 #include "TH1F.h"
17 
18 // Includes from FairRoot
19 #include "FairEventHeader.h"
20 #include "FairLogger.h"
21 #include "FairRun.h"
22 
23 // Includes from CbmRoot
24 #include "CbmEvent.h"
25 #include "CbmMCDataArray.h"
26 #include "CbmMCDataManager.h"
27 #include "CbmMCTrack.h"
28 #include "CbmStsHit.h"
29 #include "CbmStsPoint.h"
30 #include "CbmStsSetup.h"
31 #include "CbmStsTrack.h"
32 #include "CbmTrackMatchNew.h"
33 
34 
35 using std::fixed;
36 using std::right;
37 using std::setprecision;
38 using std::setw;
39 
40 
41 // ----- Default constructor -------------------------------------------
42 CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t iVerbose)
43  : FairTask("STSFindTracksQA", iVerbose)
44  , fHitMap()
45  , fMatchMap()
46  , fQualiMap()
47  , fEvents()
48  , fMCTracks(NULL)
49  , fStsPoints(NULL)
50  , fStsHits(NULL)
51  , fStsHitMatch(NULL)
52  , fStsTracks(NULL)
53  , fMatches(NULL)
54  , fLegacy(kFALSE)
55  , fTargetPos(0., 0., 0.)
56  , fSetup(NULL)
57  , fNStations(0)
58  , fMinStations(3)
59  , fQuota(0.7)
60  , fhMomAccAll(new TH1F())
61  , fhMomRecAll(new TH1F())
62  , fhMomEffAll(new TH1F())
63  , fhMomAccPrim(new TH1F())
64  , fhMomRecPrim(new TH1F())
65  , fhMomEffPrim(new TH1F())
66  , fhMomAccSec(new TH1F())
67  , fhMomRecSec(new TH1F())
68  , fhMomEffSec(new TH1F())
69  , fhNpAccAll(new TH1F())
70  , fhNpRecAll(new TH1F())
71  , fhNpEffAll(new TH1F())
72  , fhNpAccPrim(new TH1F())
73  , fhNpRecPrim(new TH1F())
74  , fhNpEffPrim(new TH1F())
75  , fhNpAccSec(new TH1F())
76  , fhNpRecSec(new TH1F())
77  , fhNpEffSec(new TH1F())
78  , fhZAccSec(new TH1F())
79  , fhZRecSec(new TH1F())
80  , fhZEffSec(new TH1F())
81  , fhNhClones(new TH1F())
82  , fhNhGhosts(new TH1F())
83  , fHistoList(new TList())
84  , fNAccAll(0)
85  , fNAccPrim(0)
86  , fNAccRef(0)
87  , fNAccSec(0)
88  , fNRecAll(0)
89  , fNRecPrim(0)
90  , fNRecRef(0)
91  , fNRecSec(0)
92  , fNGhosts(0)
93  , fNClones(0)
94  , fNEvents(0)
95  , fNEventsFailed(0)
96  , fTime(0.)
97  , fTimer() {}
98 
99 // -------------------------------------------------------------------------
100 
101 
102 // ----- Standard constructor ------------------------------------------
103 CbmStsFindTracksQa::CbmStsFindTracksQa(Int_t minStations,
104  Double_t quota,
105  Int_t iVerbose)
106  : FairTask("STSFindTracksQA", iVerbose)
107  , fHitMap()
108  , fMatchMap()
109  , fQualiMap()
110  , fEvents(NULL)
111  , fMCTracks(NULL)
112  , fStsPoints(NULL)
113  , fStsHits(NULL)
114  , fStsHitMatch(NULL)
115  , fStsTracks(NULL)
116  , fMatches(NULL)
117  , fLegacy(kFALSE)
118  , fTargetPos(0., 0., 0.)
119  , fSetup(NULL)
120  , fNStations(0)
121  , fMinStations(minStations)
122  , fQuota(quota)
123  , fhMomAccAll(new TH1F())
124  , fhMomRecAll(new TH1F())
125  , fhMomEffAll(new TH1F())
126  , fhMomAccPrim(new TH1F())
127  , fhMomRecPrim(new TH1F())
128  , fhMomEffPrim(new TH1F())
129  , fhMomAccSec(new TH1F())
130  , fhMomRecSec(new TH1F())
131  , fhMomEffSec(new TH1F())
132  , fhNpAccAll(new TH1F())
133  , fhNpRecAll(new TH1F())
134  , fhNpEffAll(new TH1F())
135  , fhNpAccPrim(new TH1F())
136  , fhNpRecPrim(new TH1F())
137  , fhNpEffPrim(new TH1F())
138  , fhNpAccSec(new TH1F())
139  , fhNpRecSec(new TH1F())
140  , fhNpEffSec(new TH1F())
141  , fhZAccSec(new TH1F())
142  , fhZRecSec(new TH1F())
143  , fhZEffSec(new TH1F())
144  , fhNhClones(new TH1F())
145  , fhNhGhosts(new TH1F())
146  , fHistoList(new TList())
147  , fNAccAll(0)
148  , fNAccPrim(0)
149  , fNAccRef(0)
150  , fNAccSec(0)
151  , fNRecAll(0)
152  , fNRecPrim(0)
153  , fNRecRef(0)
154  , fNRecSec(0)
155  , fNGhosts(0)
156  , fNClones(0)
157  , fNEvents(0)
158  , fNEventsFailed(0)
159  , fTime(0.)
160  , fTimer() {}
161 // -------------------------------------------------------------------------
162 
163 
164 // ----- Destructor ----------------------------------------------------
165 CbmStsFindTracksQa::~CbmStsFindTracksQa() {
166 
167  fHistoList->Delete();
168  delete fHistoList;
169 }
170 // -------------------------------------------------------------------------
171 
172 
173 // ----- Task execution ------------------------------------------------
174 void CbmStsFindTracksQa::Exec(Option_t* /*opt*/) {
175 
176  // If there is an event branch: do the event loop
177  if (fEvents) {
178  Int_t nEvents = fEvents->GetEntriesFast();
179  LOG(debug) << GetName() << ": found time slice with " << nEvents
180  << " events.";
181 
182  for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
183  CbmEvent* event = dynamic_cast<CbmEvent*>(fEvents->At(iEvent));
184  assert(event);
185  ProcessEvent(event);
186  }
187  }
188 
189  // If there is no event branch, process the entire tree
190  else {
191  ProcessEvent();
192  }
193 }
194 // -------------------------------------------------------------------------
195 
196 
197 // ----- Public method SetParContainers --------------------------------
198 void CbmStsFindTracksQa::SetParContainers() {}
199 // -------------------------------------------------------------------------
200 
201 
202 // ----- Public method Init --------------------------------------------
203 InitStatus CbmStsFindTracksQa::Init() {
204 
205  LOG(info) << "\n\n====================================================";
206  LOG(info) << GetName() << ": Initialising...";
207 
208  // Get RootManager
209  FairRootManager* ioman = FairRootManager::Instance();
210  assert(ioman);
211 
212  // Get STS setup
213  fSetup = CbmStsSetup::Instance();
214 
215  // Get MCDataManager
216  CbmMCDataManager* mcManager =
217  dynamic_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
218  assert(mcManager);
219 
220  // Get MCTrack array
221  fMCTracks = mcManager->InitBranch("MCTrack");
222  assert(fMCTracks);
223 
224  // Get StsPoint array
225  fStsPoints = mcManager->InitBranch("StsPoint");
226  assert(fStsPoints);
227 
228  // Get Event array
229  fEvents = dynamic_cast<TClonesArray*>(ioman->GetObject("Event"));
230 
231  // Get StsHit array
232  fStsHits = (TClonesArray*) ioman->GetObject("StsHit");
233  assert(fStsHits);
234 
235  // Get StsHitMatch array
236  fStsHitMatch = (TClonesArray*) ioman->GetObject("StsHitMatch");
237  assert(fStsHitMatch);
238 
239  // Get StsTrack array
240  fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
241  assert(fStsTracks);
242 
243  // Get StsTrackMatch array
244  fMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
245  assert(fMatches);
246 
247 
248  // Get the geometry of target and STS
249  InitStatus geoStatus = GetGeometry();
250  if (geoStatus != kSUCCESS) {
251  LOG(error) << GetName() << "::Init: Error in reading geometry!";
252  return geoStatus;
253  }
254 
255  // Create histograms
256  CreateHistos();
257  Reset();
258 
259  // Output
260  LOG(info) << " Number of STS stations : " << fNStations;
261  LOG(info) << " Target position ( " << fTargetPos.X() << ", "
262  << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm";
263  LOG(info) << " Minimum number of STS stations : " << fMinStations;
264  LOG(info) << " Matching quota : " << fQuota;
265  LOG(info) << "====================================================";
266 
267  return geoStatus;
268 }
269 // -------------------------------------------------------------------------
270 
271 
272 // ----- Public method ReInit ------------------------------------------
273 InitStatus CbmStsFindTracksQa::ReInit() {
274 
275  LOG(info) << "\n\n====================================================";
276  LOG(info) << GetName() << ": Re-initialising...";
277 
278  // Get the geometry of target and STS
279  InitStatus geoStatus = GetGeometry();
280  if (geoStatus != kSUCCESS) {
281  LOG(error) << GetName() << "::Init: Error in reading geometry!";
282  return geoStatus;
283  }
284 
285  // --- Screen log
286  LOG(info) << " Number of STS stations : " << fNStations;
287  LOG(info) << " Target position ( " << fTargetPos.X() << ", "
288  << fTargetPos.Y() << ", " << fTargetPos.Z() << ") cm";
289  LOG(info) << " Minimum number of STS stations : " << fMinStations;
290  LOG(info) << " Matching quota : " << fQuota;
291  LOG(info) << "====================================================";
292 
293  return geoStatus;
294 }
295 // -------------------------------------------------------------------------
296 
297 
298 // ----- Public method Exec --------------------------------------------
299 void CbmStsFindTracksQa::ProcessEvent(CbmEvent* event) {
300 
301  // --- Event number. Note that the FairRun counting start with 1.
302  Int_t eventNumber =
303  (event ? event->GetNumber()
304  : FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() - 1);
305 
306  LOG(debug) << GetName() << ": Process event " << eventNumber;
307 
308  // Timer
309  fTimer.Start();
310 
311  // Eventwise counters
312  // Int_t nMCTracks = 0;
313  Int_t nTracks = 0;
314  Int_t nGhosts = 0;
315  Int_t nClones = 0;
316  Int_t nAll = 0;
317  Int_t nAcc = 0;
318  Int_t nRecAll = 0;
319  Int_t nPrim = 0;
320  Int_t nRecPrim = 0;
321  Int_t nRef = 0;
322  Int_t nRecRef = 0;
323  Int_t nSec = 0;
324  Int_t nRecSec = 0;
325  TVector3 momentum;
326  TVector3 vertex;
327 
328  // Fill hit and track maps
329  FillHitMap(event);
330  FillMatchMap(event, nTracks, nGhosts, nClones);
331 
332  // Loop over MCTracks
333  Int_t nMcTracks = fMCTracks->Size(0, eventNumber);
334  for (Int_t mcTrackId = 0; mcTrackId < nMcTracks; mcTrackId++) {
335  CbmMCTrack* mcTrack =
336  dynamic_cast<CbmMCTrack*>(fMCTracks->Get(0, eventNumber, mcTrackId));
337  assert(mcTrack);
338 
339  // Continue only for reconstructible tracks
340  nAll++;
341  if (fHitMap.find(mcTrackId) == fHitMap.end()) continue; // No hits
342  Int_t nStations = fHitMap[mcTrackId].size();
343  if (nStations < fMinStations) continue; // Too few stations
344  nAcc++;
345 
346  // Check origin of MCTrack
347  // TODO: Track origin should rather be compared to MC event vertex
348  // But that is not available from MCDataManager
349  mcTrack->GetStartVertex(vertex);
350  Bool_t isPrim = kFALSE;
351  if (TMath::Abs(vertex.Z() - fTargetPos.Z()) < 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 reconstructible tracks
367  fhMomAccAll->Fill(mom);
368  fhNpAccAll->Fill(Double_t(nStations));
369  if (isPrim) {
370  fhMomAccPrim->Fill(mom);
371  fhNpAccPrim->Fill(Double_t(nStations));
372  } else {
373  fhMomAccSec->Fill(mom);
374  fhNpAccSec->Fill(Double_t(nStations));
375  fhZAccSec->Fill(vertex.Z());
376  }
377 
378  // Get matched StsTrack
379  Int_t trackId = -1;
380  Double_t quali = 0.;
381  // Bool_t isRec = kFALSE;
382  if (fMatchMap.find(mcTrackId) != fMatchMap.end()) {
383  trackId = fMatchMap[mcTrackId];
384  // isRec = kTRUE;
385  CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(trackId);
386  assert(stsTrack);
387  quali = fQualiMap[mcTrackId];
388  assert(quali >= fQuota);
389  CbmTrackMatchNew* match = (CbmTrackMatchNew*) fMatches->At(trackId);
390  assert(match);
391  Int_t nTrue = match->GetNofTrueHits();
392  Int_t nWrong = match->GetNofWrongHits();
393  //Int_t nFake = match->GetNofFakeHits();
394  Int_t nFake = 0;
395  Int_t nAllHits = stsTrack->GetNofStsHits();
396  assert(nTrue + nWrong + nFake == nAllHits);
397 
398  // Verbose output
399  LOG(debug1) << GetName() << ": MCTrack " << mcTrackId << ", stations "
400  << nStations << ", hits " << nAllHits << ", true hits "
401  << nTrue;
402 
403  // Fill histograms for reconstructed tracks
404  nRecAll++;
405  fhMomRecAll->Fill(mom);
406  fhNpRecAll->Fill(Double_t(nAllHits));
407  if (isPrim) {
408  nRecPrim++;
409  fhMomRecPrim->Fill(mom);
410  fhNpRecPrim->Fill(Double_t(nAllHits));
411  if (isRef) nRecRef++;
412  } else {
413  nRecSec++;
414  fhMomRecSec->Fill(mom);
415  fhNpRecSec->Fill(Double_t(nAllHits));
416  fhZRecSec->Fill(vertex.Z());
417  }
418 
419  } // Match found in map?
420 
421  } // Loop over MCTracks
422 
423 
424  // Calculate efficiencies
425  Double_t effAll = Double_t(nRecAll) / Double_t(nAcc);
426  Double_t effPrim = Double_t(nRecPrim) / Double_t(nPrim);
427  Double_t effRef = Double_t(nRecRef) / Double_t(nRef);
428  Double_t effSec = Double_t(nRecSec) / Double_t(nSec);
429 
430  fTimer.Stop();
431 
432 
433  // Event summary
434  LOG(info) << "+ " << setw(20) << GetName() << ": Event " << setw(6) << right
435  << fNEvents << ", real time " << fixed << setprecision(6)
436  << fTimer.RealTime() << " s, MC tracks: all " << nMcTracks
437  << ", acc. " << nAcc << ", rec. " << nRecAll << ", eff. "
438  << setprecision(2) << 100. * effAll << " %";
439  if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug)) {
440  LOG(debug) << "---------- StsFindTracksQa : Event summary ------------";
441  LOG(debug) << "MCTracks : " << nAll << ", reconstructible: " << nAcc
442  << ", reconstructed: " << nRecAll;
443  LOG(debug) << "Vertex : reconstructible: " << nPrim
444  << ", reconstructed: " << nRecPrim << ", efficiency "
445  << effPrim * 100. << "%";
446  LOG(debug) << "Reference : reconstructible: " << nRef
447  << ", reconstructed: " << nRecRef << ", efficiency "
448  << effRef * 100. << "%";
449  LOG(debug) << "Non-vertex : reconstructible: " << nSec
450  << ", reconstructed: " << nRecSec << ", efficiency "
451  << effSec * 100. << "%";
452  LOG(debug) << "STSTracks " << nTracks << ", ghosts " << nGhosts
453  << ", clones " << nClones;
454  LOG(debug)
455  << "-----------------------------------------------------------\n";
456  }
457 
458 
459  // Increase counters
460  fNAccAll += nAcc;
461  fNAccPrim += nPrim;
462  fNAccRef += nRef;
463  fNAccSec += nSec;
464  fNRecAll += nRecAll;
465  fNRecPrim += nRecPrim;
466  fNRecRef += nRecRef;
467  fNRecSec += nRecSec;
468  fNGhosts += nGhosts;
469  fNClones += nClones;
470  fNEvents++;
471  fTime += fTimer.RealTime();
472 }
473 // -------------------------------------------------------------------------
474 
475 
476 // ----- Private method Finish -----------------------------------------
477 void CbmStsFindTracksQa::Finish() {
478 
479  // Divide histograms for efficiency calculation
480  DivideHistos(fhMomRecAll, fhMomAccAll, fhMomEffAll);
481  DivideHistos(fhMomRecPrim, fhMomAccPrim, fhMomEffPrim);
482  DivideHistos(fhMomRecSec, fhMomAccSec, fhMomEffSec);
483  DivideHistos(fhNpRecAll, fhNpAccAll, fhNpEffAll);
484  DivideHistos(fhNpRecPrim, fhNpAccPrim, fhNpEffPrim);
485  DivideHistos(fhNpRecSec, fhNpAccSec, fhNpEffSec);
486  DivideHistos(fhZRecSec, fhZAccSec, fhZEffSec);
487 
488  // Normalise histos for clones and ghosts to one event
489  if (fNEvents) {
490  fhNhClones->Scale(1. / Double_t(fNEvents));
491  fhNhGhosts->Scale(1. / Double_t(fNEvents));
492  }
493 
494  // Calculate integrated efficiencies and rates
495  Double_t effAll = Double_t(fNRecAll) / Double_t(fNAccAll);
496  Double_t effPrim = Double_t(fNRecPrim) / Double_t(fNAccPrim);
497  Double_t effRef = Double_t(fNRecRef) / Double_t(fNAccRef);
498  Double_t effSec = Double_t(fNRecSec) / Double_t(fNAccSec);
499  Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNEvents);
500  Double_t rateClones = Double_t(fNClones) / Double_t(fNEvents);
501 
502  // Run summary to screen
503  std::cout << std::endl;
504  LOG(info) << "=====================================";
505  LOG(info) << fName << ": Run summary ";
506  LOG(info) << "Events processed : " << fNEvents << setprecision(2);
507  LOG(info) << "Eff. all tracks : " << effAll * 100 << " % (" << fNRecAll
508  << "/" << fNAccAll << ")";
509  LOG(info) << "Eff. vertex tracks : " << effPrim * 100 << " % ("
510  << fNRecPrim << "/" << fNAccPrim << ")";
511  LOG(info) << "Eff. reference tracks : " << effRef * 100 << " % (" << fNRecRef
512  << "/" << fNAccRef << ")";
513  LOG(info) << "Eff. secondary tracks : " << effSec * 100 << " % (" << fNRecSec
514  << "/" << fNAccSec << ")";
515  LOG(info) << "Ghost rate : " << rateGhosts << " per event";
516  LOG(info) << "Clone rate : " << rateClones << " per event";
517  LOG(info) << "Time per event : " << setprecision(6)
518  << fTime / Double_t(fNEvents) << " s";
519  LOG(info) << "=====================================";
520 
521  // Write histos to output
522  gDirectory->mkdir("STSFindTracksQA");
523  gDirectory->cd("STSFindTracksQA");
524  TIter next(fHistoList);
525  while (TH1* histo = ((TH1*) next()))
526  histo->Write();
527  gDirectory->cd("..");
528 }
529 // -------------------------------------------------------------------------
530 
531 
532 // ----- Private method GetGeometry ------------------------------------
533 InitStatus CbmStsFindTracksQa::GetGeometry() {
534 
535  // Get target geometry
536  GetTargetPosition();
537 
538  fNStations = CbmStsSetup::Instance()->GetNofStations();
539 
540 
541  return kSUCCESS;
542 }
543 // -------------------------------------------------------------------------
544 
545 
546 // ----- Get target node -----------------------------------------------
547 void CbmStsFindTracksQa::GetTargetPosition() {
548 
549  TGeoNode* target = NULL;
550 
551  gGeoManager->CdTop();
552  TGeoNode* cave = gGeoManager->GetCurrentNode();
553  for (Int_t iNode1 = 0; iNode1 < cave->GetNdaughters(); iNode1++) {
554  TString name = cave->GetDaughter(iNode1)->GetName();
555  if (name.Contains("pipe", TString::kIgnoreCase)) {
556  LOG(debug) << "Found pipe node " << name;
557  gGeoManager->CdDown(iNode1);
558  break;
559  }
560  }
561  for (Int_t iNode2 = 0;
562  iNode2 < gGeoManager->GetCurrentNode()->GetNdaughters();
563  iNode2++) {
564  TString name =
565  gGeoManager->GetCurrentNode()->GetDaughter(iNode2)->GetName();
566  if (name.Contains("pipevac1", TString::kIgnoreCase)) {
567  LOG(debug) << "Found vacuum node " << name;
568  gGeoManager->CdDown(iNode2);
569  break;
570  }
571  }
572  for (Int_t iNode3 = 0;
573  iNode3 < gGeoManager->GetCurrentNode()->GetNdaughters();
574  iNode3++) {
575  TString name =
576  gGeoManager->GetCurrentNode()->GetDaughter(iNode3)->GetName();
577  if (name.Contains("target", TString::kIgnoreCase)) {
578  LOG(debug) << "Found target node " << name;
579  gGeoManager->CdDown(iNode3);
580  target = gGeoManager->GetCurrentNode();
581  break;
582  }
583  }
584  if (!target) {
585  fTargetPos[0] = 0.;
586  fTargetPos[1] = 0.;
587  fTargetPos[2] = 0.;
588  } else {
589  TGeoHMatrix* glbMatrix = gGeoManager->GetCurrentMatrix();
590  Double_t* pos = glbMatrix->GetTranslation();
591  fTargetPos[0] = pos[0];
592  fTargetPos[1] = pos[1];
593  fTargetPos[2] = pos[2];
594  }
595 
596  gGeoManager->CdTop();
597 }
598 // -------------------------------------------------------------------------
599 
600 
601 // ----- Private method CreateHistos -----------------------------------
602 void CbmStsFindTracksQa::CreateHistos() {
603 
604  // Histogram list
605  fHistoList = new TList();
606 
607  // Momentum distributions
608  Double_t minMom = 0.;
609  Double_t maxMom = 10.;
610  Int_t nBinsMom = 40;
611  fhMomAccAll = new TH1F(
612  "hMomAccAll", "all reconstructable tracks", nBinsMom, minMom, maxMom);
613  fhMomRecAll = new TH1F(
614  "hMomRecAll", "all reconstructed tracks", nBinsMom, minMom, maxMom);
615  fhMomEffAll =
616  new TH1F("hMomEffAll", "efficiency all tracks", nBinsMom, minMom, maxMom);
617  fhMomAccPrim = new TH1F(
618  "hMomAccPrim", "reconstructable vertex tracks", nBinsMom, minMom, maxMom);
619  fhMomRecPrim = new TH1F(
620  "hMomRecPrim", "reconstructed vertex tracks", nBinsMom, minMom, maxMom);
621  fhMomEffPrim = new TH1F(
622  "hMomEffPrim", "efficiency vertex tracks", nBinsMom, minMom, maxMom);
623  fhMomAccSec = new TH1F("hMomAccSec",
624  "reconstructable non-vertex tracks",
625  nBinsMom,
626  minMom,
627  maxMom);
628  fhMomRecSec = new TH1F(
629  "hMomRecSec", "reconstructed non-vertex tracks", nBinsMom, minMom, maxMom);
630  fhMomEffSec = new TH1F(
631  "hMomEffSec", "efficiency non-vertex tracks", nBinsMom, minMom, maxMom);
632  fHistoList->Add(fhMomAccAll);
633  fHistoList->Add(fhMomRecAll);
634  fHistoList->Add(fhMomEffAll);
635  fHistoList->Add(fhMomAccPrim);
636  fHistoList->Add(fhMomRecPrim);
637  fHistoList->Add(fhMomEffPrim);
638  fHistoList->Add(fhMomAccSec);
639  fHistoList->Add(fhMomRecSec);
640  fHistoList->Add(fhMomEffSec);
641 
642  // Number-of-points distributions
643  Double_t minNp = -0.5;
644  Double_t maxNp = 15.5;
645  Int_t nBinsNp = 16;
646  fhNpAccAll =
647  new TH1F("hNpAccAll", "all reconstructable tracks", nBinsNp, minNp, maxNp);
648  fhNpRecAll =
649  new TH1F("hNpRecAll", "all reconstructed tracks", nBinsNp, minNp, maxNp);
650  fhNpEffAll =
651  new TH1F("hNpEffAll", "efficiency all tracks", nBinsNp, minNp, maxNp);
652  fhNpAccPrim = new TH1F(
653  "hNpAccPrim", "reconstructable vertex tracks", nBinsNp, minNp, maxNp);
654  fhNpRecPrim = new TH1F(
655  "hNpRecPrim", "reconstructed vertex tracks", nBinsNp, minNp, maxNp);
656  fhNpEffPrim =
657  new TH1F("hNpEffPrim", "efficiency vertex tracks", nBinsNp, minNp, maxNp);
658  fhNpAccSec = new TH1F(
659  "hNpAccSec", "reconstructable non-vertex tracks", nBinsNp, minNp, maxNp);
660  fhNpRecSec = new TH1F(
661  "hNpRecSec", "reconstructed non-vertex tracks", nBinsNp, minNp, maxNp);
662  fhNpEffSec = new TH1F(
663  "hNpEffSec", "efficiency non-vertex tracks", nBinsNp, minNp, maxNp);
664  fHistoList->Add(fhNpAccAll);
665  fHistoList->Add(fhNpRecAll);
666  fHistoList->Add(fhNpEffAll);
667  fHistoList->Add(fhNpAccPrim);
668  fHistoList->Add(fhNpRecPrim);
669  fHistoList->Add(fhNpEffPrim);
670  fHistoList->Add(fhNpAccSec);
671  fHistoList->Add(fhNpRecSec);
672  fHistoList->Add(fhNpEffSec);
673 
674  // z(vertex) distributions
675  Double_t minZ = 0.;
676  Double_t maxZ = 50.;
677  Int_t nBinsZ = 50;
678  fhZAccSec = new TH1F(
679  "hZAccSec", "reconstructable non-vertex tracks", nBinsZ, minZ, maxZ);
680  fhZRecSec = new TH1F(
681  "hZRecSecl", "reconstructed non-vertex tracks", nBinsZ, minZ, maxZ);
682  fhZEffSec =
683  new TH1F("hZEffRec", "efficiency non-vertex tracks", nBinsZ, minZ, maxZ);
684  fHistoList->Add(fhZAccSec);
685  fHistoList->Add(fhZRecSec);
686  fHistoList->Add(fhZEffSec);
687 
688  // Number-of-hit distributions
689  fhNhClones =
690  new TH1F("hNhClones", "number of hits for clones", nBinsNp, minNp, maxNp);
691  fhNhGhosts =
692  new TH1F("hNhGhosts", "number of hits for ghosts", nBinsNp, minNp, maxNp);
693  fHistoList->Add(fhNhClones);
694  fHistoList->Add(fhNhGhosts);
695 }
696 // -------------------------------------------------------------------------
697 
698 
699 // ----- Private method Reset ------------------------------------------
700 void CbmStsFindTracksQa::Reset() {
701 
702  TIter next(fHistoList);
703  while (TH1* histo = ((TH1*) next()))
704  histo->Reset();
705 
706  fNAccAll = fNAccPrim = fNAccRef = fNAccSec = 0;
707  fNRecAll = fNRecPrim = fNRecRef = fNRecSec = 0;
708  fNGhosts = fNClones = fNEvents = 0;
709 }
710 // -------------------------------------------------------------------------
711 
712 
713 // ----- Private method FillHitMap -------------------------------------
714 void CbmStsFindTracksQa::FillHitMap(CbmEvent* event) {
715 
716  // --- Event number. Note that the FairRun counting starts with 1.
717  Int_t eventNumber =
718  (event ? event->GetNumber()
719  : FairRun::Instance()->GetEventHeader()->GetMCEntryNumber() - 1);
720 
721  // --- Fill hit map ( mcTrack -> ( station -> number of hits ) )
722  fHitMap.clear();
723  Int_t nHits = (event ? event->GetNofData(ECbmDataType::kStsHit)
724  : fStsHits->GetEntriesFast());
725  for (Int_t iHit = 0; iHit < nHits; iHit++) {
726  Int_t hitIndex =
727  (event ? event->GetIndex(ECbmDataType::kStsHit, iHit) : iHit);
728  CbmStsHit* hit = (CbmStsHit*) fStsHits->At(hitIndex);
729  CbmMatch* hitMatch = (CbmMatch*) fStsHitMatch->At(hitIndex);
730  Int_t pointIndex = hitMatch->GetMatchedLink().GetIndex();
731  assert(pointIndex >= 0);
732  CbmStsPoint* stsPoint =
733  dynamic_cast<CbmStsPoint*>(fStsPoints->Get(0, eventNumber, pointIndex));
734  assert(stsPoint);
735  Int_t mcTrackIndex = stsPoint->GetTrackID();
736  Int_t station = fSetup->GetStationNumber(hit->GetAddress());
737  fHitMap[mcTrackIndex][station]++;
738  }
739  LOG(debug) << GetName() << ": Filled hit map from " << nHits
740  << " STS hits for " << fHitMap.size() << " MCTracks.";
741 }
742 // -------------------------------------------------------------------------
743 
744 
745 // ------ Private method FillMatchMap ----------------------------------
746 void CbmStsFindTracksQa::FillMatchMap(CbmEvent* event,
747  Int_t& nRec,
748  Int_t& nGhosts,
749  Int_t& nClones) {
750 
751  // Clear matching maps
752  fMatchMap.clear();
753  fQualiMap.clear();
754 
755  // Loop over StsTracks. Check matched MCtrack and fill maps.
756  nGhosts = 0;
757  nClones = 0;
758  Int_t nTracks = (event ? event->GetNofData(ECbmDataType::kStsTrack)
759  : fStsTracks->GetEntriesFast());
760 
761  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
762 
763  // --- StsTrack
764  Int_t trackIndex =
765  (event ? event->GetIndex(ECbmDataType::kStsTrack, iTrack) : iTrack);
766  CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(trackIndex);
767  assert(stsTrack);
768  Int_t nHits = stsTrack->GetNofStsHits();
769 
770  // --- TrackMatch
771  CbmTrackMatchNew* match = (CbmTrackMatchNew*) fMatches->At(trackIndex);
772  assert(match);
773  Int_t nTrue = match->GetNofTrueHits();
774 
775  // --- Matched MCTrack
776  Int_t mcTrackId = -1;
777  if (nTrue > 0) mcTrackId = match->GetMatchedLink().GetIndex();
778  if (mcTrackId < 0) {
779  fhNhGhosts->Fill(nHits);
780  nGhosts++;
781  continue;
782  }
783 
784  // Check matching criterion (quota)
785  Double_t quali = Double_t(nTrue) / Double_t(nHits);
786  if (quali >= fQuota) {
787 
788  // No previous match for this MCTrack
789  if (fMatchMap.find(mcTrackId) == fMatchMap.end()) {
790  fMatchMap[mcTrackId] = iTrack;
791  fQualiMap[mcTrackId] = quali;
792  } //? no previous match
793 
794  // Previous match; take the better one
795  else {
796  if (fQualiMap[mcTrackId] < quali) {
797  CbmStsTrack* oldTrack =
798  (CbmStsTrack*) fStsTracks->At(fMatchMap[mcTrackId]);
799  fhNhClones->Fill(Double_t(oldTrack->GetNofStsHits()));
800  fMatchMap[mcTrackId] = iTrack;
801  fQualiMap[mcTrackId] = quali;
802  } //? new track matches better to MCTrack
803 
804  else
805  fhNhClones->Fill(nHits);
806  nClones++;
807  } //? previous match found
808 
809  } //? true match ratio > quota
810 
811  // If not matched, it's a ghost
812  else {
813  fhNhGhosts->Fill(nHits);
814  nGhosts++;
815  }
816 
817  } // Loop over StsTracks
818  nRec = nTracks;
819  LOG(debug) << GetName() << ": Filled match map for " << nRec
820  << " STS tracks. Ghosts " << nGhosts << " Clones " << nClones;
821 }
822 // -------------------------------------------------------------------------
823 
824 
825 // ----- Private method DivideHistos -----------------------------------
826 void CbmStsFindTracksQa::DivideHistos(TH1* histo1, TH1* histo2, TH1* histo3) {
827 
828  if (!histo1 || !histo2 || !histo3) {
829  LOG(fatal) << GetName() << "::DivideHistos: "
830  << "NULL histogram pointer";
831  }
832 
833  Int_t nBins = histo1->GetNbinsX();
834  if (histo2->GetNbinsX() != nBins || histo3->GetNbinsX() != nBins) {
835  LOG(error) << GetName() << "::DivideHistos: "
836  << "Different bin numbers in histos";
837  LOG(error) << histo1->GetName() << " " << histo1->GetNbinsX();
838  LOG(error) << histo2->GetName() << " " << histo2->GetNbinsX();
839  LOG(error) << histo3->GetName() << " " << histo3->GetNbinsX();
840  return;
841  }
842 
843  Double_t c1, c2, c3, ce;
844  for (Int_t iBin = 0; iBin < nBins; iBin++) {
845  c1 = histo1->GetBinContent(iBin);
846  c2 = histo2->GetBinContent(iBin);
847  if (c2 != 0.) {
848  c3 = c1 / c2;
849  Double_t c4 = (c3 * (1. - c3) / c2);
850  if (c4 >= 0.) {
851  ce = TMath::Sqrt(c3 * (1. - c3) / c2);
852  } else {
853  ce = 0;
854  }
855  } else {
856  c3 = 0.;
857  ce = 0.;
858  }
859  histo3->SetBinContent(iBin, c3);
860  histo3->SetBinError(iBin, ce);
861  }
862 }
863 // -------------------------------------------------------------------------
864 
865 
866 ClassImp(CbmStsFindTracksQa)
CbmStsFindTracksQa::GetTargetPosition
void GetTargetPosition()
Definition: CbmStsFindTracksQa.cxx:547
CbmStsFindTracksQa::fhMomRecSec
TH1F * fhMomRecSec
Definition: CbmStsFindTracksQa.h:152
CbmStsFindTracksQa.h
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmStsFindTracksQa::GetGeometry
InitStatus GetGeometry()
Definition: CbmStsFindTracksQa.cxx:533
CbmStsFindTracksQa::fNRecAll
Int_t fNRecAll
Definition: CbmStsFindTracksQa.h:166
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition: CbmMCTrack.h:172
CbmStsFindTracksQa::fEvents
TClonesArray * fEvents
Definition: CbmStsFindTracksQa.h:125
CbmMatch
Definition: CbmMatch.h:22
CbmStsFindTracksQa::CbmStsFindTracksQa
CbmStsFindTracksQa(Int_t iVerbose=1)
Definition: CbmStsFindTracksQa.cxx:42
ECbmDataType::kStsHit
@ kStsHit
CbmMCDataManager.h
CbmStsFindTracksQa::fhZAccSec
TH1F * fhZAccSec
Definition: CbmStsFindTracksQa.h:156
CbmTrackMatchNew::GetNofWrongHits
Int_t GetNofWrongHits() const
Definition: CbmTrackMatchNew.h:33
CbmStsFindTracksQa::~CbmStsFindTracksQa
virtual ~CbmStsFindTracksQa()
Definition: CbmStsFindTracksQa.cxx:165
CbmStsFindTracksQa::fHistoList
TList * fHistoList
Definition: CbmStsFindTracksQa.h:161
CbmStsFindTracksQa::Exec
virtual void Exec(Option_t *opt)
Definition: CbmStsFindTracksQa.cxx:174
CbmStsSetup.h
CbmStsFindTracksQa::fMinStations
Int_t fMinStations
Definition: CbmStsFindTracksQa.h:145
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition: CbmMCDataArray.cxx:133
CbmStsFindTracksQa::Init
virtual InitStatus Init()
Definition: CbmStsFindTracksQa.cxx:203
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
CbmStsFindTracksQa::FillMatchMap
void FillMatchMap(CbmEvent *event, Int_t &nRec, Int_t &nGhosts, Int_t &nClones)
Definition: CbmStsFindTracksQa.cxx:746
CbmStsFindTracksQa::fQualiMap
std::map< Int_t, Double_t > fQualiMap
Definition: CbmStsFindTracksQa.h:121
CbmStsFindTracksQa::fTargetPos
TVector3 fTargetPos
Definition: CbmStsFindTracksQa.h:138
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition: CbmStsSetup.cxx:293
CbmStsFindTracksQa::fNAccSec
Int_t fNAccSec
Definition: CbmStsFindTracksQa.h:165
CbmMCDataArray.h
CbmStsFindTracksQa::fMCTracks
CbmMCDataArray * fMCTracks
Event.
Definition: CbmStsFindTracksQa.h:126
CbmStsFindTracksQa::fhMomAccAll
TH1F * fhMomAccAll
Definition: CbmStsFindTracksQa.h:150
CbmStsFindTracksQa::Reset
void Reset()
Definition: CbmStsFindTracksQa.cxx:700
CbmStsPoint
Definition: CbmStsPoint.h:27
CbmStsSetup::GetNofStations
Int_t GetNofStations() const
Definition: CbmStsSetup.h:90
CbmStsFindTracksQa::fSetup
CbmStsSetup * fSetup
Definition: CbmStsFindTracksQa.h:139
CbmStsFindTracksQa::fhMomEffAll
TH1F * fhMomEffAll
Definition: CbmStsFindTracksQa.h:150
CbmStsFindTracksQa::DivideHistos
void DivideHistos(TH1 *histo1, TH1 *histo2, TH1 *histo3)
Definition: CbmStsFindTracksQa.cxx:826
CbmStsFindTracksQa::fhMomEffPrim
TH1F * fhMomEffPrim
Definition: CbmStsFindTracksQa.h:151
CbmStsFindTracksQa::fhNpEffAll
TH1F * fhNpEffAll
Definition: CbmStsFindTracksQa.h:153
CbmStsFindTracksQa::fhMomAccPrim
TH1F * fhMomAccPrim
Definition: CbmStsFindTracksQa.h:151
CbmEvent.h
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition: CbmStsHit.h:31
ECbmDataType::kStsTrack
@ kStsTrack
CbmStsFindTracksQa::fStsTracks
TClonesArray * fStsTracks
StsHitMatch.
Definition: CbmStsFindTracksQa.h:130
CbmStsTrack.h
Data class for STS tracks.
CbmStsFindTracksQa::fhNpEffSec
TH1F * fhNpEffSec
Definition: CbmStsFindTracksQa.h:155
CbmStsFindTracksQa::fMatchMap
std::map< Int_t, Int_t > fMatchMap
Definition: CbmStsFindTracksQa.h:117
CbmStsFindTracksQa::fMatches
TClonesArray * fMatches
StsTrack.
Definition: CbmStsFindTracksQa.h:131
CbmStsFindTracksQa::fhNpRecAll
TH1F * fhNpRecAll
Definition: CbmStsFindTracksQa.h:153
CbmStsFindTracksQa::fNRecSec
Int_t fNRecSec
Definition: CbmStsFindTracksQa.h:166
CbmStsFindTracksQa::fhNpRecSec
TH1F * fhNpRecSec
Definition: CbmStsFindTracksQa.h:155
CbmStsFindTracksQa::fhMomRecPrim
TH1F * fhMomRecPrim
Definition: CbmStsFindTracksQa.h:151
CbmStsFindTracksQa::fNAccPrim
Int_t fNAccPrim
Definition: CbmStsFindTracksQa.h:165
CbmStsFindTracksQa::fhNpAccPrim
TH1F * fhNpAccPrim
Definition: CbmStsFindTracksQa.h:154
CbmStsFindTracksQa::fNRecPrim
Int_t fNRecPrim
Definition: CbmStsFindTracksQa.h:166
CbmStsSetup::GetStationNumber
Int_t GetStationNumber(Int_t address)
Definition: CbmStsSetup.cxx:187
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmStsFindTracksQa::ProcessEvent
void ProcessEvent(CbmEvent *event=NULL)
Definition: CbmStsFindTracksQa.cxx:299
CbmStsFindTracksQa::ReInit
virtual InitStatus ReInit()
Definition: CbmStsFindTracksQa.cxx:273
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
CbmTrackMatchNew.h
CbmStsFindTracksQa::fStsHitMatch
TClonesArray * fStsHitMatch
StsHits.
Definition: CbmStsFindTracksQa.h:129
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmStsFindTracksQa::SetParContainers
virtual void SetParContainers()
Definition: CbmStsFindTracksQa.cxx:198
CbmStsFindTracksQa::fTimer
TStopwatch fTimer
Definition: CbmStsFindTracksQa.h:173
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition: CbmMCTrack.h:182
CbmStsFindTracksQa
Definition: CbmStsFindTracksQa.h:27
CbmStsFindTracksQa::fhZEffSec
TH1F * fhZEffSec
Definition: CbmStsFindTracksQa.h:156
CbmStsFindTracksQa::CreateHistos
void CreateHistos()
Definition: CbmStsFindTracksQa.cxx:602
CbmStsFindTracksQa::fhNpAccSec
TH1F * fhNpAccSec
Definition: CbmStsFindTracksQa.h:155
CbmStsFindTracksQa::fhMomAccSec
TH1F * fhMomAccSec
Definition: CbmStsFindTracksQa.h:152
CbmStsFindTracksQa::fStsHits
TClonesArray * fStsHits
StsPoints.
Definition: CbmStsFindTracksQa.h:128
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
CbmStsFindTracksQa::fhZRecSec
TH1F * fhZRecSec
Definition: CbmStsFindTracksQa.h:156
CbmStsFindTracksQa::fNAccAll
Int_t fNAccAll
Definition: CbmStsFindTracksQa.h:165
CbmStsPoint.h
CbmStsFindTracksQa::fNRecRef
Int_t fNRecRef
Definition: CbmStsFindTracksQa.h:166
CbmMCTrack.h
CbmTrackMatchNew::GetNofTrueHits
Int_t GetNofTrueHits() const
Definition: CbmTrackMatchNew.h:32
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmStsFindTracksQa::fhNpRecPrim
TH1F * fhNpRecPrim
Definition: CbmStsFindTracksQa.h:154
CbmStsFindTracksQa::fStsPoints
CbmMCDataArray * fStsPoints
MCtrack.
Definition: CbmStsFindTracksQa.h:127
CbmStsFindTracksQa::fhNpAccAll
TH1F * fhNpAccAll
Definition: CbmStsFindTracksQa.h:153
CbmStsFindTracksQa::Finish
virtual void Finish()
Definition: CbmStsFindTracksQa.cxx:477
CbmStsFindTracksQa::fHitMap
std::map< Int_t, std::map< Int_t, Int_t > > fHitMap
Definition: CbmStsFindTracksQa.h:113
CbmStsFindTracksQa::fNStations
Int_t fNStations
Definition: CbmStsFindTracksQa.h:140
CbmEvent
Class characterising one event by a collection of links (indices) to data objects,...
Definition: CbmEvent.h:30
pos
TVector3 pos
Definition: CbmMvdSensorDigiToHitTask.cxx:60
CbmStsFindTracksQa::FillHitMap
void FillHitMap(CbmEvent *event)
Definition: CbmStsFindTracksQa.cxx:714
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmStsFindTracksQa::fhNhClones
TH1F * fhNhClones
Definition: CbmStsFindTracksQa.h:157
CbmStsFindTracksQa::fNClones
Int_t fNClones
Definition: CbmStsFindTracksQa.h:167
CbmStsFindTracksQa::fhNpEffPrim
TH1F * fhNpEffPrim
Definition: CbmStsFindTracksQa.h:154
CbmStsFindTracksQa::fNAccRef
Int_t fNAccRef
Definition: CbmStsFindTracksQa.h:165
CbmStsFindTracksQa::fhNhGhosts
TH1F * fhNhGhosts
Definition: CbmStsFindTracksQa.h:157
CbmStsFindTracksQa::fNGhosts
Int_t fNGhosts
Definition: CbmStsFindTracksQa.h:167
CbmStsTrack
Definition: CbmStsTrack.h:37
CbmStsFindTracksQa::fhMomRecAll
TH1F * fhMomRecAll
Definition: CbmStsFindTracksQa.h:150
eventNumber
Int_t eventNumber
Definition: riplet/Lx.cxx:78
CbmStsFindTracksQa::fNEvents
Int_t fNEvents
Definition: CbmStsFindTracksQa.h:168
CbmStsFindTracksQa::fhMomEffSec
TH1F * fhMomEffSec
Definition: CbmStsFindTracksQa.h:152
CbmStsFindTracksQa::fTime
Double_t fTime
Definition: CbmStsFindTracksQa.h:170
CbmStsFindTracksQa::fQuota
Double_t fQuota
Definition: CbmStsFindTracksQa.h:146
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsTrack::GetNofStsHits
Int_t GetNofStsHits() const
Definition: CbmStsTrack.h:90