FairRoot/PandaRoot
PndGemTrackFinderQA.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndGemTrackFinderQA source file -----
3 // ----- Created 02.06.2009 by R. Karabowicz -----
4 // -------------------------------------------------------------------------
5 
6 #include "PndGemTrackFinderQA.h"
7 
8 // FairRoot includes
9 #include "FairRootManager.h"
10 #include "FairRunAna.h"
11 #include "FairRuntimeDb.h"
12 #include "FairBaseParSet.h"
13 #include "FairTrackParam.h"
14 #include "FairRootManager.h"
15 // Pnd includes
16 #include "PndGemMCPoint.h"
17 #include "PndGemDigiPar.h"
18 #include "PndTrack.h"
19 #include "PndTrackCand.h"
20 #include "PndTrackCandHit.h"
21 #include "PndDetectorList.h"
22 // ROOT includes
23 #include "TClonesArray.h"
24 #include "TGeoManager.h"
25 #include "TMatrixFSym.h"
26 
27 // C++ includes
28 #include <iostream>
29 #include <iomanip>
30 #include <map>
31 #include <cmath>
32 using std::cout;
33 using std::endl;
34 using std::map;
35 
36 // ----- Default constructor ------------------------------------------
37 PndGemTrackFinderQA::PndGemTrackFinderQA() : FairTask("GEM Track Finder QA", 1) {
38  fDigiPar = NULL;
39  fMCTrackArray = NULL;
40  fMCPointArray = NULL;
41  fGemHitArray = NULL;
42  fGemTrackArray = NULL;
43  fNofEvents = 0;
44 
45  fNeededStationsToRecoTrack = -666;
46  fMinQuota = -666;
47 
48  fNofMCAll = 0;
49  fNofMCAcc = 0;
50  fNofMCPrim = 0;
51  fNofMCSec = 0;
52  fNofMCRef = 0;
53 
54  fNofRecoAcc = 0;
55  fNofRecoPrim = 0;
56  fNofRecoSec = 0;
57  fNofRecoRef = 0;
58  fNofRecoGhosts = 0;
59  fNofRecoClones = 0;
60 
61  fHistoList = NULL;
62  fhMCAllVsP = NULL;
63  fhMCAccVsP = NULL;
64  fhMCPrimVsP = NULL;
65  fhMCSecVsP = NULL;
66  fhMCRefVsP = NULL;
67  fhRecoAccVsP = NULL;
68  fhRecoPrimVsP = NULL;
69  fhRecoSecVsP = NULL;
70  fhRecoRefVsP = NULL;
71  fhEffAccVsP = NULL;
72  fhEffPrimVsP = NULL;
73  fhEffSecVsP = NULL;
74  fhEffRefVsP = NULL;
75 
76  fhMCAllVsT = NULL;
77  fhMCAccVsT = NULL;
78  fhMCPrimVsT = NULL;
79  fhMCSecVsT = NULL;
80  fhMCRefVsT = NULL;
81  fhRecoAccVsT = NULL;
82  fhRecoPrimVsT = NULL;
83  fhRecoSecVsT = NULL;
84  fhRecoRefVsT = NULL;
85  fhEffAccVsT = NULL;
86  fhEffPrimVsT = NULL;
87  fhEffSecVsT = NULL;
88  fhEffRefVsT = NULL;
89 
90  fhMCAllVsA = NULL;
91  fhMCAccVsA = NULL;
92  fhMCPrimVsA = NULL;
93  fhMCSecVsA = NULL;
94  fhMCRefVsA = NULL;
95  fhRecoAccVsA = NULL;
96  fhRecoPrimVsA = NULL;
97  fhRecoSecVsA = NULL;
98  fhRecoRefVsA = NULL;
99  fhEffAccVsA = NULL;
100  fhEffPrimVsA = NULL;
101  fhEffSecVsA = NULL;
102  fhEffRefVsA = NULL;
103 
104  fhMCAllVsN = NULL;
105  fhMCAccVsN = NULL;
106  fhMCPrimVsN = NULL;
107  fhMCSecVsN = NULL;
108  fhMCRefVsN = NULL;
109  fhRecoAccVsN = NULL;
110  fhRecoPrimVsN = NULL;
111  fhRecoSecVsN = NULL;
112  fhRecoRefVsN = NULL;
113  fhEffAccVsN = NULL;
114  fhEffPrimVsN = NULL;
115  fhEffSecVsN = NULL;
116  fhEffRefVsN = NULL;
117 
118  fhRecoAllP = NULL;
119  fhRecoPrimP = NULL;
120  fhRecoSecP = NULL;
121  fhRecoAllT = NULL;
122  fhRecoPrimT = NULL;
123  fhRecoSecT = NULL;
124  fhRecoAllA = NULL;
125  fhRecoPrimA = NULL;
126  fhRecoSecA = NULL;
127 
128  fhMomResAccVsP = NULL;
129  fhMomResPrimVsP = NULL;
130  fhMomResSecVsP = NULL;
131  fhMomResRefVsP = NULL;
132  fhMomResAccVsT = NULL;
133  fhMomResPrimVsT = NULL;
134  fhMomResSecVsT = NULL;
135  fhMomResRefVsT = NULL;
136  fhMomResAccVsA = NULL;
137  fhMomResPrimVsA = NULL;
138  fhMomResSecVsA = NULL;
139  fhMomResRefVsA = NULL;
140 
141  fhNofHitsPerTrack = NULL;
142  fhNofHitsPerRecoTrack = NULL;
143  fhNofHitsPerGhost = NULL;
144  fhNofHitsPerClone = NULL;
145 
146  fhNofCorrHitsPerRecoTrack = NULL;
147  fhNofOthTHitsPerRecoTrack = NULL;
148  fhNofNoTrHitsPerRecoTrack = NULL;
149 
150  fhNofMCTracksPerEvent = NULL;
151  fhNofRecoTracksPerEvent = NULL;
152 }
153 // -------------------------------------------------------------------------
154 
155 
156 
157 // ----- Standard constructor ------------------------------------------
158 PndGemTrackFinderQA::PndGemTrackFinderQA(Int_t iVerbose)
159  : FairTask("GEM Track Finder QA", iVerbose) {
160  fDigiPar = NULL;
161  fMCTrackArray = NULL;
162  fMCPointArray = NULL;
163  fGemHitArray = NULL;
164  fGemTrackArray = NULL;
165  fNofEvents = 0;
166 
167  fNeededStationsToRecoTrack = -666;
168  fMinQuota = -666;
169 
170  fNofMCAll = 0;
171  fNofMCAcc = 0;
172  fNofMCPrim = 0;
173  fNofMCSec = 0;
174  fNofMCRef = 0;
175 
176  fNofRecoAcc = 0;
177  fNofRecoPrim = 0;
178  fNofRecoSec = 0;
179  fNofRecoRef = 0;
180  fNofRecoGhosts = 0;
181  fNofRecoClones = 0;
182 
183  fHistoList = NULL;
184  fhMCAllVsP = NULL;
185  fhMCAccVsP = NULL;
186  fhMCPrimVsP = NULL;
187  fhMCSecVsP = NULL;
188  fhMCRefVsP = NULL;
189  fhRecoAccVsP = NULL;
190  fhRecoPrimVsP = NULL;
191  fhRecoSecVsP = NULL;
192  fhRecoRefVsP = NULL;
193  fhEffAccVsP = NULL;
194  fhEffPrimVsP = NULL;
195  fhEffSecVsP = NULL;
196  fhEffRefVsP = NULL;
197 
198  fhMCAllVsT = NULL;
199  fhMCAccVsT = NULL;
200  fhMCPrimVsT = NULL;
201  fhMCSecVsT = NULL;
202  fhMCRefVsT = NULL;
203  fhRecoAccVsT = NULL;
204  fhRecoPrimVsT = NULL;
205  fhRecoSecVsT = NULL;
206  fhRecoRefVsT = NULL;
207  fhEffAccVsT = NULL;
208  fhEffPrimVsT = NULL;
209  fhEffSecVsT = NULL;
210  fhEffRefVsT = NULL;
211 
212  fhMCAllVsA = NULL;
213  fhMCAccVsA = NULL;
214  fhMCPrimVsA = NULL;
215  fhMCSecVsA = NULL;
216  fhMCRefVsA = NULL;
217  fhRecoAccVsA = NULL;
218  fhRecoPrimVsA = NULL;
219  fhRecoSecVsA = NULL;
220  fhRecoRefVsA = NULL;
221  fhEffAccVsA = NULL;
222  fhEffPrimVsA = NULL;
223  fhEffSecVsA = NULL;
224  fhEffRefVsA = NULL;
225 
226  fhMCAllVsN = NULL;
227  fhMCAccVsN = NULL;
228  fhMCPrimVsN = NULL;
229  fhMCSecVsN = NULL;
230  fhMCRefVsN = NULL;
231  fhRecoAccVsN = NULL;
232  fhRecoPrimVsN = NULL;
233  fhRecoSecVsN = NULL;
234  fhRecoRefVsN = NULL;
235  fhEffAccVsN = NULL;
236  fhEffPrimVsN = NULL;
237  fhEffSecVsN = NULL;
238  fhEffRefVsN = NULL;
239 
240  fhRecoAllP = NULL;
241  fhRecoPrimP = NULL;
242  fhRecoSecP = NULL;
243  fhRecoAllT = NULL;
244  fhRecoPrimT = NULL;
245  fhRecoSecT = NULL;
246  fhRecoAllA = NULL;
247  fhRecoPrimA = NULL;
248  fhRecoSecA = NULL;
249 
250  fhMomResAccVsP = NULL;
251  fhMomResPrimVsP = NULL;
252  fhMomResSecVsP = NULL;
253  fhMomResRefVsP = NULL;
254  fhMomResAccVsT = NULL;
255  fhMomResPrimVsT = NULL;
256  fhMomResSecVsT = NULL;
257  fhMomResRefVsT = NULL;
258  fhMomResAccVsA = NULL;
259  fhMomResPrimVsA = NULL;
260  fhMomResSecVsA = NULL;
261  fhMomResRefVsA = NULL;
262 
263  fhNofHitsPerTrack = NULL;
264  fhNofHitsPerRecoTrack = NULL;
265  fhNofHitsPerGhost = NULL;
266  fhNofHitsPerClone = NULL;
267 
268  fhNofCorrHitsPerRecoTrack = NULL;
269  fhNofOthTHitsPerRecoTrack = NULL;
270  fhNofNoTrHitsPerRecoTrack = NULL;
271 
272  fhNofMCTracksPerEvent = NULL;
273  fhNofRecoTracksPerEvent = NULL;
274 }
275 // -------------------------------------------------------------------------
276 
277 // ----- Destructor ----------------------------------------------------
278 PndGemTrackFinderQA::~PndGemTrackFinderQA() { }
279 
280 // ----- Init -----------------------------------------------------------
281 InitStatus PndGemTrackFinderQA::Init() {
282 
283  // Get and check FairRootManager
284  FairRootManager* ioman = FairRootManager::Instance();
285 
286 #if (ROOT_VERSION_CODE >= ROOT_VERSION(5,34,10))
287  ioman->SetUseFairLinks(kTRUE);
288 #endif
289 
290  if( !ioman ) {
291  cout << "-E- "<< GetName() <<"::Init: "
292  << "RootManager not instantised!" << endl;
293  return kERROR;
294  }
295 
296  // Get the pointer to the singleton FairRunAna object
297  FairRunAna* ana = FairRunAna::Instance();
298  if(NULL == ana) {
299  cout << "-E- "<< GetName() <<"::Init :"
300  <<" no FairRunAna object!" << endl;
301  return kERROR;
302  }
303  // Get the pointer to run-time data base
304  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
305  if(NULL == rtdb) {
306  cout << "-E- "<< GetName() <<"::Init :"
307  <<" no runtime database!" << endl;
308  return kERROR;
309  }
310 
311  // Get MCTrack array
312  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
313  if( !fMCTrackArray ) {
314  cout << "-E- "<< GetName() <<"::Init: No MCTrack array!"
315  << endl;
316  return kERROR;
317  }
318 
319  // Get PndGemPoint (MCPoint) array
320  fMCPointArray = (TClonesArray*) ioman->GetObject("GEMPoint");
321  if( !fMCPointArray ) {
322  cout << "-E- "<< GetName() <<"::Init: No MCPoint array!"
323  << endl;
324  return kERROR;
325  }
326 
327  // Get Gem hit Array
328  fGemHitArray = (TClonesArray*) ioman->GetObject("GEMHit");
329  if ( !fGemHitArray ) {
330  cout << "-E- " << GetName() << "::Init: No PndGemHit array!" << endl;
331  return kERROR;
332  }
333 
334  // Get Gem track Array
335  fGemTrackArray = (TClonesArray*) ioman->GetObject("GEMTrack");
336  if ( !fGemTrackArray ) {
337  cout << "-E- " << GetName() << "::Init: No PndGemTrack array!" << endl;
338  return kERROR;
339  }
340 
341 
342  // Get GEM digitisation parameter container
343  fDigiPar = (PndGemDigiPar*)(rtdb->getContainer("PndGemDetectors"));
344 
345  cout << "-I- " << fName.Data() << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
346  cout << "-I- " << fName.Data() << "::Init(). Initialization succesfull." << endl;
347 
348  TList* branchNames = FairRootManager::Instance()->GetBranchNameList();
349 
350  Int_t nofGemBranches = 0;
351  fGemPointNumber = -1;
352  for (int i = 0; i < branchNames->GetEntries(); i++){
353  TObjString* branchName = (TObjString*)branchNames->At(i);
354  if ( !branchName->GetString().BeginsWith("GEM") ) continue;
355  if ( branchName->GetString() == "GEMPoint" ) fGemPointNumber = i;
356  fGemData[nofGemBranches] = (TClonesArray*) ioman->GetObject(branchName->GetString().Data());
357  fGemDataPointer[i] = nofGemBranches;
358  cout << "GEM Data [" << nofGemBranches << "] = " << branchName->GetString().Data() << " referred to as " << i << (fGemPointNumber==i?" ***":"") << endl;
359  nofGemBranches++;
360  }
361 
362 
363  CreateHistos();
364 
365  return kSUCCESS;
366 }
367 // -------------------------------------------------------------------------
368 
369 // ----- Private method ReInit -----------------------------------------
370 InitStatus PndGemTrackFinderQA::ReInit() {
371 
372  return kSUCCESS;
373 }
374 // -------------------------------------------------------------------------
375 
376 // ----- Private method SetParContainers -------------------------------
377 void PndGemTrackFinderQA::SetParContainers() {
378 
379  // Get run and runtime database
380  FairRunAna* run = FairRunAna::Instance();
381  if ( ! run ) Fatal("SetParContainers", "No analysis run");
382 
383  FairRuntimeDb* db = run->GetRuntimeDb();
384  if ( ! db ) Fatal("SetParContainers", "No runtime database");
385 
386  // Get GEM digitisation parameter container
387  fDigiPar = (PndGemDigiPar*)(db->getContainer("PndGemDetectors"));
388 
389  fNofEvents = 0;
390  fTargetPos.SetXYZ(0.,0.,0.);
391  fNeededStationsToRecoTrack = 3;
392  fMinQuota = 0.6;
393 
394  fNofMCAll = 0;
395  fNofMCAcc = 0;
396  fNofMCPrim = 0;
397  fNofMCSec = 0;
398  fNofMCRef = 0; // primary, mom.mag > 0.5GeV/c
399 
400  fNofRecoAcc = 0;
401  fNofRecoPrim = 0;
402  fNofRecoSec = 0;
403  fNofRecoRef = 0; // primary, mom.mag > 0.5GeV/c
404 
405  fNofRecoGhosts = 0;
406  fNofRecoClones = 0;
407 
408 }
409 // -------------------------------------------------------------------------
410 
411 
412 // ----- Public method Exec --------------------------------------------
413 void PndGemTrackFinderQA::Exec(Option_t*) {
414 
415  fNofEvents++;
416 
417  Int_t nofMCAll = 0;
418  Int_t nofMCAcc = 0;
419  Int_t nofMCPrim = 0;
420  Int_t nofMCSec = 0;
421  Int_t nofMCRef = 0; // primary, mom.mag > 0.5GeV/c
422 
423  Int_t nofRecoAcc = 0;
424  Int_t nofRecoPrim = 0;
425  Int_t nofRecoSec = 0;
426  Int_t nofRecoRef = 0; // primary, mom.mag > 0.5GeV/c
427 
428  Int_t nofRecoClones = 0;
429  Int_t nofRecoGhosts = 0;
430 
431  PrepareMCTracks();
432  MatchRecoTracks();
433 
434  Int_t nofMCTracks = fMCTrackArray->GetEntriesFast();
435  Int_t nofRecoTracks = fGemTrackArray->GetEntriesFast();
436 
437  fhNofMCTracksPerEvent ->Fill(fNofEvents-1,nofMCTracks);
438  fhNofRecoTracksPerEvent->Fill(fNofEvents-1,nofRecoTracks);
439 
440  PndMCTrack* mcTrack;
441  PndTrack* gemTrack;
442  //PndGemHit* gemHit; //[R.K. 01/2017] unused variable?
443 
444  for ( Int_t imct = 0 ; imct < nofMCTracks ; imct++ ) {
445 
446  mcTrack = (PndMCTrack*) fMCTrackArray->At(imct);
447 
448  if ( mcTrack->GetNPoints(kGEM) == 0 ) continue;
449 
450  TVector3 mcVertex = mcTrack->GetStartVertex();
451  Bool_t isPrim = kFALSE;
452  TVector3 mcMomVec = mcTrack->GetMomentum();
453  Double_t mcMomMag = mcMomVec.Mag();
454  Double_t mcMomThe = mcMomVec.Theta()*TMath::RadToDeg();
455  Double_t mcMomPhi = mcMomVec.Phi()*TMath::RadToDeg();
456  Int_t mcPoints = fMCTrackNofGemPoints[imct];
457 
458  if ( (mcVertex-fTargetPos).Mag() < 1. )
459  isPrim = kTRUE;
460 
461  fhMCAllVsP->Fill(mcMomMag);
462  fhMCAllVsT->Fill(mcMomThe);
463  fhMCAllVsA->Fill(mcMomPhi);
464  fhMCAllVsN->Fill(mcPoints);
465 
466  // check reconstrubility, continue only for reconstruable tracks
467  nofMCAll ++;
468  if ( fMCTrackNofCrossedGemStations[imct] < fNeededStationsToRecoTrack ) continue;
469  nofMCAcc ++;
470 
471  fhMCAccVsP->Fill(mcMomMag);
472  fhMCAccVsT->Fill(mcMomThe);
473  fhMCAccVsA->Fill(mcMomPhi);
474  fhMCAccVsN->Fill(mcPoints);
475 
476  if ( isPrim ) {
477  nofMCPrim ++;
478  fhMCPrimVsP->Fill(mcMomMag);
479  fhMCPrimVsT->Fill(mcMomThe);
480  fhMCPrimVsA->Fill(mcMomPhi);
481  fhMCPrimVsN->Fill(mcPoints);
482  }
483  else {
484  nofMCSec ++;
485  fhMCSecVsP->Fill(mcMomMag);
486  fhMCSecVsT->Fill(mcMomThe);
487  fhMCSecVsA->Fill(mcMomPhi);
488  fhMCSecVsN->Fill(mcPoints);
489  }
490 
491  Bool_t isRefP = kFALSE;
492  Bool_t isRefT = kFALSE;
493  if ( isPrim && mcMomMag > 0.5 ) {
494  isRefP = kTRUE;
495  fhMCRefVsT->Fill(mcMomThe);
496  }
497  if ( isPrim && mcMomThe > 5. && mcMomThe < 20. ) {
498  isRefT = kTRUE;
499  fhMCRefVsP->Fill(mcMomMag);
500  }
501  if ( isRefP && isRefT ) {
502  nofMCRef ++;
503  fhMCRefVsN->Fill(mcPoints);
504  fhMCRefVsA->Fill(mcMomPhi);
505  }
506 
507  Bool_t mcTrackReconstructed = kFALSE;
508  for ( Int_t irtr = 0 ; irtr < nofRecoTracks ; irtr++ ) {
509  if ( fRecoTrackMCMatch[irtr] != imct ) continue;
510  mcTrackReconstructed = kTRUE;
511  gemTrack = (PndTrack*) fGemTrackArray->At(irtr);
512 
513  TVector3 recoTrackMom = gemTrack->GetParamFirst().GetMomentum();
514  Double_t recoMomMag = recoTrackMom.Mag();
515 
516  nofRecoAcc ++;
517 
518  fhRecoAccVsP->Fill(mcMomMag);
519  fhRecoAccVsT->Fill(mcMomThe);
520  fhRecoAccVsA->Fill(mcMomPhi);
521  fhRecoAccVsN->Fill(mcPoints);
522  fhMomResAccVsP->Fill(mcMomMag,100.*(mcMomMag-recoMomMag)/mcMomMag);
523  fhMomResAccVsT->Fill(mcMomThe,100.*(mcMomMag-recoMomMag)/mcMomMag);
524  fhMomResAccVsA->Fill(mcMomPhi,100.*(mcMomMag-recoMomMag)/mcMomMag);
525 
526  if ( isRefP ) {
527  fhRecoRefVsT->Fill(mcMomThe);
528  fhMomResRefVsT->Fill(mcMomThe,100.*(mcMomMag-recoMomMag)/mcMomMag);
529  }
530  if ( isRefT ) {
531  fhRecoRefVsP->Fill(mcMomMag);
532  fhMomResRefVsP->Fill(mcMomMag,100.*(mcMomMag-recoMomMag)/mcMomMag);
533  }
534  if ( isRefP && isRefT ) {
535  nofRecoRef ++;
536  fhRecoRefVsN->Fill(mcPoints);
537  fhRecoRefVsA->Fill(mcMomPhi);
538  }
539  if ( isPrim ) {
540  nofRecoPrim ++;
541  fhRecoPrimVsP->Fill(mcMomMag);
542  fhRecoPrimVsT->Fill(mcMomThe);
543  fhRecoPrimVsA->Fill(mcMomPhi);
544  fhRecoPrimVsN->Fill(mcPoints);
545  fhMomResPrimVsP->Fill(mcMomMag,100.*(mcMomMag-recoMomMag)/mcMomMag);
546  fhMomResPrimVsT->Fill(mcMomThe,100.*(mcMomMag-recoMomMag)/mcMomMag);
547  fhMomResPrimVsA->Fill(mcMomPhi,100.*(mcMomMag-recoMomMag)/mcMomMag);
548 
549  fhRecoPrimP->Fill(recoTrackMom.Mag());
550  fhRecoPrimT->Fill(recoTrackMom.Theta()*TMath::RadToDeg());
551  fhRecoPrimA->Fill(recoTrackMom.Phi()*TMath::RadToDeg());
552  }
553  else {
554  nofRecoSec ++;
555  fhRecoSecVsP->Fill(mcMomMag);
556  fhRecoSecVsT->Fill(mcMomThe);
557  fhRecoSecVsA->Fill(mcMomPhi);
558  fhRecoSecVsN->Fill(mcPoints);
559  fhMomResSecVsP->Fill(mcMomMag,100.*(mcMomMag-recoMomMag)/mcMomMag);
560  fhMomResSecVsT->Fill(mcMomThe,100.*(mcMomMag-recoMomMag)/mcMomMag);
561  fhMomResSecVsA->Fill(mcMomPhi,100.*(mcMomMag-recoMomMag)/mcMomMag);
562 
563  fhRecoSecP->Fill(recoTrackMom.Mag());
564  fhRecoSecT->Fill(recoTrackMom.Theta()*TMath::RadToDeg());
565  fhRecoSecA->Fill(recoTrackMom.Phi()*TMath::RadToDeg());
566  }
567  break; //do not include clones, the break is good for efficiency plots but bad(?) for momentum resolution
568  }
569  if ( mcTrackReconstructed == kFALSE && fVerbose > 0 ) {
570  cout << "MC TRACK WAS NOT RECONSTRUCTED"
571  << " mom = " << mcMomMag
572  << " the = " << mcMomThe
573  << " phi = " << mcMomPhi
574  << endl;
575  }
576  }
577 
578  for ( Int_t irtr = 0 ; irtr < nofRecoTracks ; irtr++ ) {
579  gemTrack = (PndTrack*) fGemTrackArray->At(irtr);
580 
581  fhNofHitsPerTrack->Fill(gemTrack->GetTrackCand().GetNHits());
582 
583  if ( fRecoTrackMCMatch[irtr] == -1 ) {
584  nofRecoGhosts ++;
585  fhNofHitsPerGhost->Fill(gemTrack->GetTrackCand().GetNHits());
586  continue;
587  }
588  fhNofHitsPerRecoTrack->Fill(gemTrack->GetTrackCand().GetNHits());
589 
590  Int_t nofCorrHits = 0;
591  Int_t nofOthTHits = 0;
592  Int_t nofNoTrHits = 0;
593  for ( size_t ihit = 0 ; ihit < gemTrack->GetTrackCand().GetNHits() ; ihit++ ) {
594  PndTrackCandHit tch = gemTrack->GetTrackCand().GetSortedHit(ihit);
595  Int_t bestPointIndex = FindMatchingPoint(tch.GetHitId());
596 
597  if ( bestPointIndex == -1 ) { nofNoTrHits++; continue; }
598  if ( bestPointIndex != fRecoTrackMCMatch[irtr] ) { nofCorrHits++; continue; }
599  nofCorrHits ++;
600  }
601 
602  fhNofCorrHitsPerRecoTrack->Fill(nofCorrHits);
603  fhNofOthTHitsPerRecoTrack->Fill(nofOthTHits);
604  fhNofNoTrHitsPerRecoTrack->Fill(nofNoTrHits);
605 
606  for ( Int_t irtr2 = 0 ; irtr2 < nofRecoTracks ; irtr2++ ) {
607  if ( irtr == irtr2 ) continue;
608  if ( fRecoTrackMCMatch[irtr2] == fRecoTrackMCMatch[irtr] ) {
609  fhNofHitsPerClone->Fill(gemTrack->GetTrackCand().GetNHits());
610  nofRecoClones ++;
611  break;
612  }
613  }
614  }
615 
616  fNofMCAll += nofMCAll;
617  fNofMCAcc += nofMCAcc;
618  fNofMCPrim += nofMCPrim;
619  fNofMCSec += nofMCSec;
620  fNofMCRef += nofMCRef;
621 
622  fNofRecoAcc += nofRecoAcc;
623  fNofRecoPrim += nofRecoPrim;
624  fNofRecoSec += nofRecoSec;
625  fNofRecoRef += nofRecoRef;
626 
627  fNofRecoGhosts += nofRecoGhosts;
628  fNofRecoClones += nofRecoClones;
629 
630  if ( fVerbose ) {
631  Double_t effAcc = 1.;
632  Double_t effPrim = 1.;
633  Double_t effSec = 1.;
634  Double_t effRef = 1.;
635  if ( nofMCAcc ) effAcc = 100.*(Double_t)nofRecoAcc /((Double_t)nofMCAcc);
636  if ( nofMCPrim ) effPrim = 100.*(Double_t)nofRecoPrim/((Double_t)nofMCPrim);
637  if ( nofMCSec ) effSec = 100.*(Double_t)nofRecoSec /((Double_t)nofMCSec);
638  if ( nofMCRef ) effRef = 100.*(Double_t)nofRecoRef /((Double_t)nofMCRef);
639  cout << "---- PndGemTrackFinderQA : Event " << fNofEvents << " summary -----" << endl;
640  cout << "MC Tracks: " << nofMCTracks << endl;
641  cout << " reconstruable: " << nofMCAcc << " reconstructed " << nofRecoAcc << " >>>> " << effAcc << endl;
642  cout << " primaries : " << nofMCPrim << " reconstructed " << nofRecoPrim << " >>>> " << effPrim << endl;
643  cout << " secondaries : " << nofMCSec << " reconstructed " << nofRecoSec << " >>>> " << effSec << endl;
644  cout << " reference : " << nofMCRef << " reconstructed " << nofRecoRef << " >>>> " << effRef << endl;
645  cout << " ghosts : " << nofRecoGhosts << " clones " << nofRecoClones << endl;
646  cout << "---------------------------------------------------------------" << endl;
647  }
648 
649 }
650 // ------------------------------------------------------------
651 
652 // ----- Private method PrepareMCTrack --------------------------------------------
653 void PndGemTrackFinderQA::PrepareMCTracks() {
654 
655  fMCTrackNofCrossedGemStations.clear();
656  fMCTrackNofCrossedGemStations.resize(fMCTrackArray->GetEntriesFast(),0);
657  fMCTrackNofGemPoints.clear();
658  fMCTrackNofGemPoints.resize(fMCTrackArray->GetEntriesFast(),0);
659 
660  const Int_t nofGemPoints = fMCPointArray->GetEntriesFast();
661 
662  const Int_t nofMCTracks = fMCTrackArray->GetEntriesFast();
663  const Int_t nofGemStations = fDigiPar->GetNStations();
664  Int_t nofPointsPerStation[nofMCTracks][nofGemStations];
665  for ( Int_t imct = 0 ; imct < nofMCTracks ; imct++ )
666  for ( Int_t is = 0 ; is < nofGemStations ; is++ )
667  nofPointsPerStation[imct][is] = 0;
668 
669  PndGemMCPoint* mcPoint;
670  for ( Int_t imcp = 0 ; imcp < nofGemPoints ; imcp++ ) {
671  mcPoint = (PndGemMCPoint*)fMCPointArray->At(imcp);
672 
673  Int_t stationNr = fDigiPar->GetStationNr(mcPoint->GetSensorId());
674 
675  nofPointsPerStation[mcPoint->GetTrackID()][stationNr-1] += 1;
676 
677  fMCTrackNofGemPoints[mcPoint->GetTrackID()] += 1;
678  }
679 
680  Int_t nofCGM = 0;
681  for ( Int_t imct = 0 ; imct < nofMCTracks ; imct++ ) {
682  // fMCTrackNofCrossedGemStations[imct] = nofCGM;
683  nofCGM = 0;
684  for ( Int_t is = 0 ; is < nofGemStations ; is++ )
685  if ( nofPointsPerStation[imct][is] > 0 )
686  nofCGM++;
687  fMCTrackNofCrossedGemStations[imct] = nofCGM;
688  }
689 }
690 // ------------------------------------------------------------
691 
692 // ----- Private method PrepareMCTrack --------------------------------------------
693 void PndGemTrackFinderQA::MatchRecoTracks() {
694 
695  Int_t nofRecoTracks = fGemTrackArray->GetEntriesFast();
696 
697  fRecoTrackMCMatch.clear();
698  fRecoTrackMCMatch.resize(nofRecoTracks,-1);
699 
700  //PndGemHit* gemHit; //[R.K. 01/2017] unused variable?
701  FairMCPoint* mcPoint;
702  PndTrack* gemTrack;
703 
704  //const Int_t nsh = 2*fDigiPar->GetNStations(); //[R.K. 01/2017] unused variable?
705 
706  for ( Int_t irtr = 0 ; irtr < nofRecoTracks ; irtr++ ) {
707  gemTrack = (PndTrack*) fGemTrackArray->At(irtr);
708  vector<Int_t> nofTrMCId(500,0);
709 
710  for ( size_t ihit = 0 ; ihit < gemTrack->GetTrackCand().GetNHits() ; ihit++ ) {
711  PndTrackCandHit tch = gemTrack->GetTrackCand().GetSortedHit(ihit);
712 
713  Int_t bestPointIndex = FindMatchingPoint(tch.GetHitId());
714 
715  if ( bestPointIndex == -1 ) continue;
716  mcPoint = (FairMCPoint*) fMCPointArray->At(bestPointIndex);
717  nofTrMCId[mcPoint->GetTrackID()] += 1;
718  }
719  Int_t bestMCId = -1;
720  Int_t largestNofMCId = 0;
721  for ( Int_t itm = 0 ; itm < 500 ; itm++ ) {
722  if ( nofTrMCId[itm] == largestNofMCId ) { bestMCId = -1; }
723  if ( nofTrMCId[itm] > largestNofMCId ) { bestMCId = itm; largestNofMCId = nofTrMCId[itm]; }
724  }
725  if ( bestMCId == -1 ) {
726  if ( fVerbose ) {
727  cout << "GHOST TRACK WITH"
728  << " mom = " << gemTrack->GetParamFirst().GetMomentum().Mag()
729  << " the = " << gemTrack->GetParamFirst().GetMomentum().Theta()*TMath::RadToDeg()
730  << " phi = " << gemTrack->GetParamFirst().GetMomentum().Phi()*TMath::RadToDeg()
731  << endl;
732  }
733  continue;
734  }
735 
736  if ( largestNofMCId < fMinQuota*gemTrack->GetTrackCand().GetNHits() ) continue;
737 
738  fRecoTrackMCMatch[irtr] = bestMCId;
739 
740  TVector3 recoTrackMom = gemTrack->GetParamFirst().GetMomentum();
741 
742  fhRecoAllP->Fill(recoTrackMom.Mag());
743  fhRecoAllT->Fill(recoTrackMom.Theta()*TMath::RadToDeg());
744  fhRecoAllA->Fill(recoTrackMom.Phi()*TMath::RadToDeg());
745  }
746 
747 }
748 // ------------------------------------------------------------
749 
750 
751 // ----- Private method GetPointVector, recurrence function --------------------------------------------
752 Int_t PndGemTrackFinderQA::FindMatchingPoint(Int_t gemHitIndex) {
753  Bool_t printMCMatching = kFALSE;
754 
755  PndGemHit* gemHit = (PndGemHit*)fGemHitArray->At(gemHitIndex);
756  if ( printMCMatching ) {
757  cout << "---> hit " << gemHitIndex << " has " << gemHit->GetNLinks() << " links" << endl;
758  for ( Int_t ilink = 0 ; ilink < gemHit->GetNLinks() ; ilink++) {
759  std::cout << " --> " << gemHit->GetLink(ilink).GetEntry()
760  << " --> " << gemHit->GetLink(ilink).GetType()
761  << " --> " << gemHit->GetLink(ilink).GetIndex()
762  << std::endl;
763  }
764  }
765 
766 
767  Int_t bestPointIndex = -1;
768  Double_t bestPointValue = 0.;
769 
770  for ( Int_t ilink = 0 ; ilink < gemHit->GetNLinks() ; ilink++) {
771  if ( gemHit->GetLink(ilink).GetType() == fGemPointNumber ) {
772  return gemHit->GetLink(ilink).GetIndex();
773  }
774  }
775 
776  // if ( gemHit->GetNLinks() != 2 ) {
777  // cout << "THERE ARE " << gemHit->GetNLinks() << " FOR HIT " << gemHitIndex << " IN EVENT " << fNofEvents << endl;
778  // }
779  // if ( gemHit->GetNLinks() == 2 ) {
780  // IT SHOULD BE REDONE TO TAKE THE BEST MATCHING MC POINT
781  Int_t maxPnt0 = -1, maxPnt1 = -1;
782  std::vector<Int_t> pointVector0;
783  std::vector<Int_t> pointVector1;
784  GetPointVector(gemHit->GetLink(0).GetType(),gemHit->GetLink(0).GetIndex(),pointVector0,printMCMatching);
785  GetPointVector(gemHit->GetLink(1).GetType(),gemHit->GetLink(1).GetIndex(),pointVector1,printMCMatching);
786  if ( printMCMatching )
787  cout << "VECT0: (" << gemHit->GetLink(0).GetIndex() << ") " << flush;
788  for ( size_t ipnt = 0 ; ipnt < pointVector0.size() ; ipnt++ ) {
789  if ( printMCMatching )
790  cout << pointVector0[ipnt] << " . " << flush;
791  if ( maxPnt0 < pointVector0[ipnt] ) {
792  maxPnt0 = pointVector0[ipnt];
793  }
794  }
795  if ( printMCMatching )
796  cout << "\b\b" << endl;
797  if ( printMCMatching )
798  cout << "VECT1: (" << gemHit->GetLink(1).GetIndex() << ") " << flush;
799  for ( size_t ipnt = 0 ; ipnt < pointVector1.size() ; ipnt++ ) {
800  if ( printMCMatching )
801  cout << pointVector1[ipnt] << " . " << flush;
802  if ( maxPnt1 < pointVector1[ipnt] ) {
803  maxPnt1 = pointVector1[ipnt];
804  }
805  }
806  if ( printMCMatching )
807  cout << "\b\b" << endl;
808  if ( printMCMatching )
809  cout << "highest points are " << maxPnt0 << " , " << maxPnt1 << endl;
810  std::vector<Int_t> countPointV0(maxPnt0+1,0);
811  std::vector<Int_t> countPointV1(maxPnt1+1,0);
812  for ( size_t ipnt = 0 ; ipnt < pointVector0.size() ; ipnt++ ) {
813  ++countPointV0[pointVector0[ipnt]];
814  }
815  for ( size_t ipnt = 0 ; ipnt < pointVector1.size() ; ipnt++ ) {
816  ++countPointV1[pointVector1[ipnt]];
817  }
818  if ( maxPnt0 > maxPnt1 )
819  maxPnt0 = maxPnt1;
820  for ( Int_t ipnt = 0 ; ipnt < maxPnt0+1 ; ipnt++ ) {
821  if ( printMCMatching )
822  cout << ipnt << " - " << countPointV0[ipnt] << " ? " << countPointV1[ipnt] << endl;
823  if ( countPointV0[ipnt] > 0 ) {
824  if ( countPointV1[ipnt] > 0 ) {
825  Double_t tempMatch = 100.*Double_t(countPointV0[ipnt]*countPointV1[ipnt])/(Double_t(pointVector0.size()*pointVector1.size()));
826  if ( bestPointValue < tempMatch ) {
827  bestPointValue = tempMatch;
828  bestPointIndex = ipnt;
829  }
830  }
831  }
832  }
833 
834  if ( printMCMatching ) {
835  cout << "RETURNING " << bestPointIndex << " (with probability " << bestPointValue << " %)" << endl;
836  }
837  return bestPointIndex;
838 }
839 // ------------------------------------------------------------
840 
841 // ----- Private method GetPointVector, recurrence function --------------------------------------------
842 Int_t PndGemTrackFinderQA::GetPointVector(Int_t arrayId, Int_t entryId, std::vector<Int_t>& pointVector, Bool_t printInfo) {
843  if ( printInfo ) {
844  for ( Int_t itemp = 0 ; itemp < 10 - arrayId ; itemp++ ) cout << " " << flush;
845  cout << "GetPointVector(" << arrayId << ", " << entryId << ")" << endl;
846  }
847  if ( arrayId == fGemPointNumber ) {
848  pointVector.push_back(entryId);
849  return pointVector.size();
850  }
851  if ( arrayId == 0 ) { // this is probably MCTrack!
852  return 0;
853  }
854 
855  FairMultiLinkedData* tempData = (FairMultiLinkedData*)fGemData[fGemDataPointer[arrayId]]->At(entryId);
856  for ( Int_t ilink = 0 ; ilink < tempData->GetNLinks() ; ilink++ ) {
857  GetPointVector(tempData->GetLink(ilink).GetType(),tempData->GetLink(ilink).GetIndex(),pointVector,printInfo);
858  }
859  return pointVector.size();
860 }
861 // ------------------------------------------------------------
862 
863 // ----- Private method CreateHistos --------------------------------------------
864 void PndGemTrackFinderQA::CreateHistos() {
865  fHistoList = new TList();
866  // number of mc tracks, reco tracks, efficiency as function of MOMENTUM
867  Double_t minMom = 0.;
868  Double_t maxMom = 10.;
869  Int_t momBins = 40;
870  fhMCAllVsP = new TH1F("hMCAllVsP" ,"all mc tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
871  fhMCAccVsP = new TH1F("hMCAccVsP" ,"reconstruable mc tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
872  fhMCPrimVsP = new TH1F("hMCPrimVsP","primary mc tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
873  fhMCSecVsP = new TH1F("hMCSecVsP" ,"secondary mc tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
874  fhRecoAccVsP = new TH1F("hRecoAccVsP" ,"reconstructed tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
875  fhRecoPrimVsP = new TH1F("hRecoPrimVsP","reco primary tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
876  fhRecoSecVsP = new TH1F("hRecoSecVsP" ,"reco secondary tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
877  fhEffAccVsP = new TH1F("hEffAccVsP" ,"eff all tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
878  fhEffPrimVsP = new TH1F("hEffPrimVsP","eff primary tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
879  fhEffSecVsP = new TH1F("hEffSecVsP" ,"eff secondary tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
880  fhMCRefVsP = new TH1F("hMCRefVsP" ,"reference mc tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
881  fhRecoRefVsP = new TH1F("hRecoRefVsP" ,"reco reference tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
882  fhEffRefVsP = new TH1F("hEffRefVsP" ,"eff reference tracks;p [GeV/c];yield [a.u.]",momBins,minMom,maxMom);
883 
884  fHistoList->Add(fhMCAllVsP);
885  fHistoList->Add(fhMCAccVsP);
886  fHistoList->Add(fhMCPrimVsP);
887  fHistoList->Add(fhMCSecVsP);
888  fHistoList->Add(fhRecoAccVsP);
889  fHistoList->Add(fhRecoPrimVsP);
890  fHistoList->Add(fhRecoSecVsP);
891  fHistoList->Add(fhEffAccVsP);
892  fHistoList->Add(fhEffPrimVsP);
893  fHistoList->Add(fhEffSecVsP);
894  fHistoList->Add(fhMCRefVsP);
895  fHistoList->Add(fhRecoRefVsP);
896  fHistoList->Add(fhEffRefVsP);
897 
898  // number of mc tracks, reco tracks, efficiency as function of THETA
899  Double_t minThe = 0.;
900  Double_t maxThe = 40.;
901  Int_t theBins = 40;
902  fhMCAllVsT = new TH1F("hMCAllVsT" ,"all mc tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
903  fhMCAccVsT = new TH1F("hMCAccVsT" ,"reconstruable mc tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
904  fhMCPrimVsT = new TH1F("hMCPrimVsT","primary mc tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
905  fhMCSecVsT = new TH1F("hMCSecVsT" ,"secondary mc tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
906  fhRecoAccVsT = new TH1F("hRecoAccVsT" ,"reconstructed tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
907  fhRecoPrimVsT = new TH1F("hRecoPrimVsT","reco primary tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
908  fhRecoSecVsT = new TH1F("hRecoSecVsT" ,"reco secondary tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
909  fhEffAccVsT = new TH1F("hEffAccVsT" ,"eff all tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
910  fhEffPrimVsT = new TH1F("hEffPrimVsT","eff primary tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
911  fhEffSecVsT = new TH1F("hEffSecVsT" ,"eff secondary tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
912  fhMCRefVsT = new TH1F("hMCRefVsT" ,"reference mc tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
913  fhRecoRefVsT = new TH1F("hRecoRefVsT" ,"reco reference tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
914  fhEffRefVsT = new TH1F("hEffRefVsT" ,"eff reference tracks;theta [deg];yield [a.u.]",theBins,minThe,maxThe);
915 
916  fHistoList->Add(fhMCAllVsT);
917  fHistoList->Add(fhMCAccVsT);
918  fHistoList->Add(fhMCPrimVsT);
919  fHistoList->Add(fhMCSecVsT);
920  fHistoList->Add(fhRecoAccVsT);
921  fHistoList->Add(fhRecoPrimVsT);
922  fHistoList->Add(fhRecoSecVsT);
923  fHistoList->Add(fhEffAccVsT);
924  fHistoList->Add(fhEffPrimVsT);
925  fHistoList->Add(fhEffSecVsT);
926  fHistoList->Add(fhMCRefVsT);
927  fHistoList->Add(fhRecoRefVsT);
928  fHistoList->Add(fhEffRefVsT);
929 
930  // number of mc tracks, reco tracks, efficiency as function of PHI
931  Double_t minPhi = -180.;
932  Double_t maxPhi = 180.;
933  Int_t phiBins = 72;
934  fhMCAllVsA = new TH1F("hMCAllVsA" ,"all mc tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
935  fhMCAccVsA = new TH1F("hMCAccVsA" ,"reconstruable mc tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
936  fhMCPrimVsA = new TH1F("hMCPrimVsA","primary mc tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
937  fhMCSecVsA = new TH1F("hMCSecVsA" ,"secondary mc tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
938  fhRecoAccVsA = new TH1F("hRecoAccVsA" ,"reconstructed tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
939  fhRecoPrimVsA = new TH1F("hRecoPrimVsA","reco primary tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
940  fhRecoSecVsA = new TH1F("hRecoSecVsA" ,"reco secondary tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
941  fhEffAccVsA = new TH1F("hEffAccVsA" ,"eff all tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
942  fhEffPrimVsA = new TH1F("hEffPrimVsA","eff primary tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
943  fhEffSecVsA = new TH1F("hEffSecVsA" ,"eff secondary tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
944  fhMCRefVsA = new TH1F("hMCRefVsA" ,"reference mc tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
945  fhRecoRefVsA = new TH1F("hRecoRefVsA" ,"reco reference tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
946  fhEffRefVsA = new TH1F("hEffRefVsA" ,"eff reference tracks;phi [deg];yield [a.u.]",phiBins,minPhi,maxPhi);
947 
948  fHistoList->Add(fhMCAllVsA);
949  fHistoList->Add(fhMCAccVsA);
950  fHistoList->Add(fhMCPrimVsA);
951  fHistoList->Add(fhMCSecVsA);
952  fHistoList->Add(fhRecoAccVsA);
953  fHistoList->Add(fhRecoPrimVsA);
954  fHistoList->Add(fhRecoSecVsA);
955  fHistoList->Add(fhEffAccVsA);
956  fHistoList->Add(fhEffPrimVsA);
957  fHistoList->Add(fhEffSecVsA);
958  fHistoList->Add(fhMCRefVsA);
959  fHistoList->Add(fhRecoRefVsA);
960  fHistoList->Add(fhEffRefVsA);
961 
962  // number of mc tracks, reco tracks, efficiency as function of NUMBER OF POINTS
963  Double_t minPnt = -0.5;
964  Double_t maxPnt = 40.5;
965  Int_t pntBins = 41;
966  fhMCAllVsN = new TH1F("hMCAllVsN" ,"all mc tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
967  fhMCAccVsN = new TH1F("hMCAccVsN" ,"reconstruable mc tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
968  fhMCPrimVsN = new TH1F("hMCPrimVsN","primary mc tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
969  fhMCSecVsN = new TH1F("hMCSecVsN" ,"secondary mc tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
970  fhRecoAccVsN = new TH1F("hRecoAccVsN" ,"reconstructed tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
971  fhRecoPrimVsN = new TH1F("hRecoPrimVsN","reco primary tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
972  fhRecoSecVsN = new TH1F("hRecoSecVsN" ,"reco secondary tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
973  fhEffAccVsN = new TH1F("hEffAccVsN" ,"eff all tracks;# points;efficiency [%]",pntBins,minPnt,maxPnt);
974  fhEffPrimVsN = new TH1F("hEffPrimVsN","eff primary tracks;# points;efficiency [%]",pntBins,minPnt,maxPnt);
975  fhEffSecVsN = new TH1F("hEffSecVsN" ,"eff secondary tracks;# points;efficiency [%]",pntBins,minPnt,maxPnt);
976  fhMCRefVsN = new TH1F("hMCRefVsN" ,"reference mc tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
977  fhRecoRefVsN = new TH1F("hRecoRefVsN" ,"reco reference tracks;# points;yield [a.u.]",pntBins,minPnt,maxPnt);
978  fhEffRefVsN = new TH1F("hEffRefVsN" ,"eff reference tracks;# points;efficiency [%]",pntBins,minPnt,maxPnt);
979 
980  fHistoList->Add(fhMCAllVsN);
981  fHistoList->Add(fhMCAccVsN);
982  fHistoList->Add(fhMCPrimVsN);
983  fHistoList->Add(fhMCSecVsN);
984  fHistoList->Add(fhRecoAccVsN);
985  fHistoList->Add(fhRecoPrimVsN);
986  fHistoList->Add(fhRecoSecVsN);
987  fHistoList->Add(fhEffAccVsN);
988  fHistoList->Add(fhEffPrimVsN);
989  fHistoList->Add(fhEffSecVsN);
990  fHistoList->Add(fhMCRefVsN);
991  fHistoList->Add(fhRecoRefVsN);
992  fHistoList->Add(fhEffRefVsN);
993 
994  // momentum resolution vs MOMENTUM
995  fhMomResAccVsP = new TH2F("hMomResAccVsP" ,"momentum resolution for all tracks;p [GeV/c];#delta{p}/p [%]",momBins,minMom,maxMom,400,-20.,20.);
996  fhMomResPrimVsP = new TH2F("hMomResPrimVsP","momentum resolution for primary tracks;p [GeV/c];#delta{p}/p [%]",momBins,minMom,maxMom,400,-20.,20.);
997  fhMomResSecVsP = new TH2F("hMomResSecVsP" ,"momentum resolution for secondary tracks;p [GeV/c];#delta{p}/p [%]",momBins,minMom,maxMom,400,-20.,20.);
998  fhMomResRefVsP = new TH2F("hMomResRefVsP" ,"momentum resolution for reference tracks;p [GeV/c];#delta{p}/p [%]",momBins,minMom,maxMom,400,-20.,20.);
999  fHistoList->Add(fhMomResAccVsP);
1000  fHistoList->Add(fhMomResPrimVsP);
1001  fHistoList->Add(fhMomResSecVsP);
1002  fHistoList->Add(fhMomResRefVsP);
1003  // momentum resolution vs THETA
1004  fhMomResAccVsT = new TH2F("hMomResAccVsT" ,"momentum resolution for all tracks;theta [deg];#delta{p}/p [%]",theBins,minThe,maxThe,400,-20.,20.);
1005  fhMomResPrimVsT = new TH2F("hMomResPrimVsT","momentum resolution for primary tracks;theta [deg];#delta{p}/p [%]",theBins,minThe,maxThe,400,-20.,20.);
1006  fhMomResSecVsT = new TH2F("hMomResSecVsT" ,"momentum resolution for secondary tracks;theta [deg];#delta{p}/p [%]",theBins,minThe,maxThe,400,-20.,20.);
1007  fhMomResRefVsT = new TH2F("hMomResRefVsT" ,"momentum resolution for reference tracks;theta [deg];#delta{p}/p [%]",theBins,minThe,maxThe,400,-20.,20.);
1008  fHistoList->Add(fhMomResAccVsT);
1009  fHistoList->Add(fhMomResPrimVsT);
1010  fHistoList->Add(fhMomResSecVsT);
1011  fHistoList->Add(fhMomResRefVsT);
1012  // momentum resolution vs PHI
1013  fhMomResAccVsA = new TH2F("hMomResAccVsA" ,"momentum resolution for all tracks;phi [deg];#delta{p}/p [%]",phiBins,minPhi,maxPhi,400,-20.,20.);
1014  fhMomResPrimVsA = new TH2F("hMomResPrimVsA","momentum resolution for primary tracks;phi [deg];#delta{p}/p [%]",phiBins,minPhi,maxPhi,400,-20.,20.);
1015  fhMomResSecVsA = new TH2F("hMomResSecVsA" ,"momentum resolution for secondary tracks;phi [deg];#delta{p}/p [%]",phiBins,minPhi,maxPhi,400,-20.,20.);
1016  fhMomResRefVsA = new TH2F("hMomResRefVsA" ,"momentum resolution for reference tracks;phi [deg];#delta{p}/p [%]",phiBins,minPhi,maxPhi,400,-20.,20.);
1017  fHistoList->Add(fhMomResAccVsA);
1018  fHistoList->Add(fhMomResPrimVsA);
1019  fHistoList->Add(fhMomResSecVsA);
1020  fHistoList->Add(fhMomResRefVsA);
1021 
1022  fhRecoAllP = new TH1F("hRecoAllP" ,"reconstructed track momentum for all",2000,0.,20.);
1023  fhRecoPrimP = new TH1F("hRecoPrimP","reconstructed track momentum for primaries",2000,0.,20.);
1024  fhRecoSecP = new TH1F("hRecoSecP" ,"reconstructed track momentum for secondaries",2000,0.,20.);
1025  fhRecoAllT = new TH1F("hRecoAllT" ,"reconstructed track theta for all",400,0.,40.);
1026  fhRecoPrimT = new TH1F("hRecoPrimT","reconstructed track theta for primaries",400,0.,40.);
1027  fhRecoSecT = new TH1F("hRecoSecT" ,"reconstructed track theta for secondaries",400,0.,40.);
1028  fhRecoAllA = new TH1F("hRecoAllA" ,"reconstructed track phi for all",360,-180.,180.);
1029  fhRecoPrimA = new TH1F("hRecoPrimA","reconstructed track phi for primaries",360,-180.,180.);
1030  fhRecoSecA = new TH1F("hRecoSecA" ,"reconstructed track phi for secondaries",360,-180.,180.);
1031  fHistoList->Add(fhRecoAllP);
1032  fHistoList->Add(fhRecoPrimP);
1033  fHistoList->Add(fhRecoSecP);
1034  fHistoList->Add(fhRecoAllT);
1035  fHistoList->Add(fhRecoPrimT);
1036  fHistoList->Add(fhRecoSecT);
1037  fHistoList->Add(fhRecoAllA);
1038  fHistoList->Add(fhRecoPrimA);
1039  fHistoList->Add(fhRecoSecA);
1040 
1041  fhNofHitsPerTrack = new TH1F("hNofHitsPerTrack","nof hits per track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1042  fhNofHitsPerRecoTrack = new TH1F("hNofHitsPerRecoTrack","nof hits per reco track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1043  fhNofHitsPerGhost = new TH1F("hNofHitsPerGhost","nof hits per ghost track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1044  fhNofHitsPerClone = new TH1F("hNofHitsPerClone","nof hits per clone track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1045  fHistoList->Add(fhNofHitsPerTrack);
1046  fHistoList->Add(fhNofHitsPerRecoTrack);
1047  fHistoList->Add(fhNofHitsPerGhost);
1048  fHistoList->Add(fhNofHitsPerClone);
1049 
1050  fhNofCorrHitsPerRecoTrack = new TH1F("hNofCorrHitsPerRecoTrack","nof hits with correct MCId per track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1051  fhNofOthTHitsPerRecoTrack = new TH1F("hNofOthTHitsPerRecoTrack","nof hits with different MCId per track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1052  fhNofNoTrHitsPerRecoTrack = new TH1F("hNofNoTrHitsPerRecoTrack","nof hits without any MCId per track;# hits;yield [a.u.]",pntBins,minPnt,maxPnt);
1053  fHistoList->Add(fhNofCorrHitsPerRecoTrack);
1054  fHistoList->Add(fhNofOthTHitsPerRecoTrack);
1055  fHistoList->Add(fhNofNoTrHitsPerRecoTrack);
1056 
1057 
1058  fhNofMCTracksPerEvent = new TH1F("hNofMCTracksPerEvent","nof MC tracks per event;event no;nof mc tracks",100001,-0.5,100000.5);
1059  fhNofRecoTracksPerEvent = new TH1F("hNofRecoTracksPerEvent","nof reco tracks per event;event no;nof mc tracks",100001,-0.5,100000.5);
1060  fHistoList->Add(fhNofMCTracksPerEvent);
1061  fHistoList->Add(fhNofRecoTracksPerEvent);
1062 }
1063 // ------------------------------------------------------------
1064 
1065 // ----- Private method DivideHistos --------------------------------------------
1066 void PndGemTrackFinderQA::DivideHistos(TH1* hist1, TH1* hist2, TH1* hist3) {
1067  hist1->Sumw2();
1068  hist2->Sumw2();
1069  hist3->Divide(hist1,hist2,1,1,"B");
1070 }
1071 // ------------------------------------------------------------
1072 
1073 // ----- Private method Finish --------------------------------------------
1074 void PndGemTrackFinderQA::Finish() {
1075  Double_t effAcc = 0.;
1076  Double_t effPrim = 0.;
1077  Double_t effSec = 0.;
1078  Double_t effRef = 0.;
1079  if ( fNofMCAcc ) effAcc = 100.*(Double_t)fNofRecoAcc /((Double_t)fNofMCAcc);
1080  if ( fNofMCPrim ) effPrim = 100.*(Double_t)fNofRecoPrim/((Double_t)fNofMCPrim);
1081  if ( fNofMCSec ) effSec = 100.*(Double_t)fNofRecoSec /((Double_t)fNofMCSec);
1082  if ( fNofMCRef ) effRef = 100.*(Double_t)fNofRecoRef /((Double_t)fNofMCRef);
1083 
1084  Double_t ghPerEv = 0.;
1085  Double_t ghPerTr = 0.;
1086  Double_t clPerEv = 0.;
1087  Double_t clPerTr = 0.;
1088  if ( fNofEvents ) ghPerEv = (Double_t)fNofRecoGhosts/((Double_t)fNofEvents);
1089  if ( fNofMCAll ) ghPerTr = (Double_t)fNofRecoGhosts/((Double_t)fNofMCAll);
1090  if ( fNofEvents ) clPerEv = (Double_t)fNofRecoClones/((Double_t)fNofEvents);
1091  if ( fNofMCAll ) clPerTr = (Double_t)fNofRecoClones/((Double_t)fNofMCAll);
1092 
1093  cout << "-------------------- PndGemTrackFinderQA : Summary ------------------" << endl;
1094  cout << " Events: " << setw(10) << fNofEvents << endl;
1095  cout << " MC Tracks: " << setw(10) << fNofMCAll << endl;
1096  cout << " reconstruable: " << setw(10) << fNofMCAcc << " reconstructed: " << setw(10) << fNofRecoAcc << " >>>> " << effAcc << "%" << endl;
1097  cout << " primaries : " << setw(10) << fNofMCPrim << " reconstructed: " << setw(10) << fNofRecoPrim << " >>>> " << effPrim << "%" << endl;
1098  cout << " reference : " << setw(10) << fNofMCRef << " reconstructed: " << setw(10) << fNofRecoRef << " >>>> " << effRef << "%" << endl;
1099  cout << " secondaries : " << setw(10) << fNofMCSec << " reconstructed: " << setw(10) << fNofRecoSec << " >>>> " << effSec << "%" << endl;
1100  cout << " ghosts : " << setw(10) << fNofRecoGhosts << " >>> " << setw(10) << ghPerEv << " per event >>> " << setw(10) << ghPerTr << " per MC Track" << endl;
1101  cout << " clones : " << setw(10) << fNofRecoClones << " >>> " << setw(10) << clPerEv << " per event >>> " << setw(10) << clPerTr << " per MC Track" << endl;
1102  cout << "---------------------------------------------------------------------" << endl;
1103 
1104 
1105  DivideHistos(fhRecoAccVsP, fhMCAccVsP, fhEffAccVsP);
1106  DivideHistos(fhRecoPrimVsP,fhMCPrimVsP,fhEffPrimVsP);
1107  DivideHistos(fhRecoSecVsP, fhMCSecVsP, fhEffSecVsP);
1108  DivideHistos(fhRecoRefVsP, fhMCRefVsP, fhEffRefVsP);
1109  fhEffAccVsP ->Scale(100.);
1110  fhEffPrimVsP->Scale(100.);
1111  fhEffSecVsP ->Scale(100.);
1112  fhEffRefVsP ->Scale(100.);
1113  DivideHistos(fhRecoAccVsT, fhMCAccVsT, fhEffAccVsT);
1114  DivideHistos(fhRecoPrimVsT,fhMCPrimVsT,fhEffPrimVsT);
1115  DivideHistos(fhRecoSecVsT, fhMCSecVsT, fhEffSecVsT);
1116  DivideHistos(fhRecoRefVsT, fhMCRefVsT, fhEffRefVsT);
1117  fhEffAccVsT ->Scale(100.);
1118  fhEffPrimVsT->Scale(100.);
1119  fhEffSecVsT ->Scale(100.);
1120  fhEffRefVsT ->Scale(100.);
1121  DivideHistos(fhRecoAccVsA, fhMCAccVsA, fhEffAccVsA);
1122  DivideHistos(fhRecoPrimVsA,fhMCPrimVsA,fhEffPrimVsA);
1123  DivideHistos(fhRecoSecVsA, fhMCSecVsA, fhEffSecVsA);
1124  DivideHistos(fhRecoRefVsA, fhMCRefVsA, fhEffRefVsA);
1125  fhEffAccVsA ->Scale(100.);
1126  fhEffPrimVsA->Scale(100.);
1127  fhEffSecVsA ->Scale(100.);
1128  fhEffRefVsA ->Scale(100.);
1129  DivideHistos(fhRecoAccVsN, fhMCAccVsN, fhEffAccVsN);
1130  DivideHistos(fhRecoPrimVsN,fhMCPrimVsN,fhEffPrimVsN);
1131  DivideHistos(fhRecoSecVsN, fhMCSecVsN, fhEffSecVsN);
1132  DivideHistos(fhRecoRefVsN, fhMCRefVsN, fhEffRefVsN);
1133  fhEffAccVsN ->Scale(100.);
1134  fhEffPrimVsN->Scale(100.);
1135  fhEffSecVsN ->Scale(100.);
1136  fhEffRefVsN ->Scale(100.);
1137 
1138 
1139  TFile* temp = gFile;
1140  FairRootManager* ioman = FairRootManager::Instance();
1141  gFile = ioman->GetOutFile();
1142  gDirectory = (TDirectory*)gFile;
1143 
1144  gDirectory->mkdir("GemTrackFinderQA");
1145  gDirectory->cd("GemTrackFinderQA");
1146  TIter next(fHistoList);
1147  while ( TH1* histo = ((TH1*)next()) ) histo->Write();
1148  gDirectory->cd("..");
1149 
1150  gFile = temp;
1151 
1152 }
1153 // ------------------------------------------------------------
1154 
1155 
1156 ClassImp(PndGemTrackFinderQA)
PndGemTrackFinderQA::fMCTrackNofGemPoints
std::vector< Int_t > fMCTrackNofGemPoints
Definition: PndGemTrackFinderQA.h:72
fVerbose
int fVerbose
Definition: poormantracks.C:24
run
Int_t run
Definition: autocutx.C:47
PndGemTrackFinderQA::SetParContainers
virtual void SetParContainers()
Definition: PndGemTrackFinderQA.cxx:377
i
Int_t i
Definition: run_full.C:25
map
PndTransMap * map
Definition: sim_emc_apd.C:99
PndGemTrackFinderQA::fMCPointArray
TClonesArray * fMCPointArray
Definition: PndGemTrackFinderQA.h:57
PndMCTrack::GetNPoints
Int_t GetNPoints(DetectorId detId) const
Definition: PndMCTrack.cxx:120
PndTrackCand::GetSortedHit
PndTrackCandHit GetSortedHit(UInt_t i)
Definition: PndTrackCand.h:54
PndGemTrackFinderQA::DivideHistos
void DivideHistos(TH1 *hist1, TH1 *hist2, TH1 *hist3)
Definition: PndGemTrackFinderQA.cxx:1066
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndGemDigiPar
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
PndGemTrackFinderQA::fDigiPar
PndGemDigiPar * fDigiPar
Definition: PndGemTrackFinderQA.h:48
hist3
TH1D * hist3
Definition: hist-t7.C:84
PndGemTrackFinderQA::fGemTrackArray
TClonesArray * fGemTrackArray
Output array of PndGemTracks.
Definition: PndGemTrackFinderQA.h:60
PndGemTrackFinderQA::fMCTrackNofCrossedGemStations
std::vector< Int_t > fMCTrackNofCrossedGemStations
Definition: PndGemTrackFinderQA.h:71
At
cout<< "blue = Monte Carlo "<< endl;cout<< "red = Helix Hit "<< endl;cout<< "green = Center Of Tubes "<< endl;for(Int_t k=0;k< track->GetEntriesFast();k++){PndSttTrack *stttrack=(PndSttTrack *) track-> At(k)
Definition: checkhelixhit.C:56
PndGemTrackFinderQA::fNofEvents
Int_t fNofEvents
event counter
Definition: PndGemTrackFinderQA.h:68
PndGemMCPoint::GetSensorId
Int_t GetSensorId() const
Definition: PndGemMCPoint.h:90
PndGemTrackFinderQA::ReInit
virtual InitStatus ReInit()
Definition: PndGemTrackFinderQA.cxx:370
PndTrack::GetTrackCand
PndTrackCand GetTrackCand()
Definition: PndTrack.h:47
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndGemTrackFinderQA::fGemHitArray
TClonesArray * fGemHitArray
Definition: PndGemTrackFinderQA.h:59
hist2
TH1D * hist2
Definition: hist-t7.C:82
PndGemTrackFinderQA::FindMatchingPoint
Int_t FindMatchingPoint(Int_t gemHitIndex)
Definition: PndGemTrackFinderQA.cxx:752
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
is
static int is
Definition: ranlxd.cxx:374
PndGemDigiPar::GetNStations
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
PndGemTrackFinderQA::fRecoTrackMCMatch
std::vector< Int_t > fRecoTrackMCMatch
Definition: PndGemTrackFinderQA.h:73
PndGemTrackFinderQA::Exec
virtual void Exec(Option_t *opt)
Definition: PndGemTrackFinderQA.cxx:413
PndGemDigiPar::GetStationNr
Int_t GetStationNr(Int_t sensorId)
Definition: PndGemDigiPar.h:54
PndTrackCand::GetNHits
UInt_t GetNHits() const
Definition: PndTrackCand.h:59
ClassImp
ClassImp(PndAnaContFact)
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
kGEM
Definition: production/tdrct/eta_c/track_check.C:4
PndGemTrackFinderQA
track finding quality assesment task
Definition: PndGemTrackFinderQA.h:30
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndGemTrackFinderQA::fMCTrackArray
TClonesArray * fMCTrackArray
Definition: PndGemTrackFinderQA.h:56
next
static int next[96]
Definition: ranlxd.cxx:374
PndTrackCandHit::GetHitId
Int_t GetHitId() const
Definition: PndTrackCandHit.h:55
PndGemTrackFinderQA::GetPointVector
Int_t GetPointVector(Int_t arrayId, Int_t entryId, std::vector< Int_t > &pointVector, Bool_t printInfo=kFALSE)
Definition: PndGemTrackFinderQA.cxx:842
PndTrack::GetParamFirst
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
PndGemTrackFinderQA::fGemData
TClonesArray * fGemData[10]
Definition: PndGemTrackFinderQA.h:63
PndGemTrackFinderQA::Init
virtual InitStatus Init()
Definition: PndGemTrackFinderQA.cxx:281
hist1
TH1D * hist1
Definition: hist-t7.C:80
Generated on Sun Aug 18 2019 03:01:29 for FairRoot/PandaRoot by   doxygen 1.8.8