FairRoot/PandaRoot
PndGemFindHitsQA.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndGemFindHitsQA source file -----
3 // ----- Created 02.06.2009 by R. Karabowicz -----
4 // -------------------------------------------------------------------------
5 
6 #include "PndGemFindHitsQA.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 "PndGemDigi.h"
18 #include "PndGemDigiPar.h"
19 #include "PndTrack.h"
20 #include "PndTrackCand.h"
21 #include "PndTrackCandHit.h"
22 #include "PndDetectorList.h"
23 // ROOT includes
24 #include "TClonesArray.h"
25 #include "TGeoManager.h"
26 #include "TMatrixFSym.h"
27 
28 // C++ includes
29 #include <iostream>
30 #include <iomanip>
31 #include <map>
32 #include <cmath>
33 using std::cout;
34 using std::endl;
35 using std::map;
36 using std::pair;
37 
38 // ----- Default constructor ------------------------------------------
40  : FairTask("GEM Find Hits QA", 1)
41  ,fDigiPar (NULL)
42  ,fMCPointArray (NULL)
43  ,fGemHitArray (NULL)
44  ,fNofEvents (0)
45  ,fHistoList (NULL)
46  ,fHistPlaneDivs (0.0)
47  ,fPointEffDist (0.1)
48  ,fhTrueMatchDiXY (NULL)
49  ,fhTrueMatchDist (NULL)
50  ,fhTrueMatchValue (NULL)
51  ,fhTrueMatchDistValue (NULL)
52 
53  ,fhTrueMatchNofPerHit (NULL)
54  ,fhTrueMatchNofPerPoint (NULL)
55 {
56  for ( Int_t istat = 0 ; istat < 4 ; istat++ ) {
57  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
58  fHistWidth[istat][isens]=0.;
59  for ( Int_t iregx = 0 ; iregx < 4 ; iregx++ ) {
60  for ( Int_t iregy = 0 ; iregy < 4 ; iregy++ ) {
61  fhPointToHit[istat][isens][iregx][iregy] = NULL;
62  }
63  }
64  // fhPointClosest [istat][isens] = NULL;
65  fhPointNof [istat][isens] = NULL;
66  fhPointReco [istat][isens] = NULL;
67  fhPointRecoEff [istat][isens] = NULL;
68  fhPointRadNof [istat][isens] = NULL;
69  fhPointRadReco [istat][isens] = NULL;
70  fhPointRadRecoEff [istat][isens] = NULL;
71 
72  fhPointMatch [istat][isens] = NULL;
73  fhPointMatchEff [istat][isens] = NULL;
74  fhPointRadMatch [istat][isens] = NULL;
75  fhPointRadMatchEff [istat][isens] = NULL;
76 
77  fhHitNof [istat][isens] = NULL;
78  fhHitFake [istat][isens] = NULL;
79  fhHitFakeProb [istat][isens] = NULL;
80  fhHitRadNof [istat][isens] = NULL;
81  fhHitRadFake [istat][isens] = NULL;
82  fhHitRadFakeProb [istat][isens] = NULL;
83  fhHitMultipleRate [istat][isens] = NULL;
84  fhCloseHits [istat][isens] = NULL;
85  fhTrueMatchDiXYPerSt[istat][isens] = NULL;
86  fhTrueMatchDistPerSt[istat][isens] = NULL;
87  }
88  }
89 }
90 // -------------------------------------------------------------------------
91 
92 
93 
94 // ----- Standard constructor ------------------------------------------
96  : FairTask("GEM Find Hits QA", iVerbose)
97  ,fDigiPar (NULL)
98  ,fMCPointArray (NULL)
99  ,fGemHitArray (NULL)
100  ,fNofEvents (0)
101  ,fHistoList (NULL)
102  ,fHistPlaneDivs (0.0)
103  ,fPointEffDist (0.1)
104  ,fhTrueMatchDiXY (NULL)
105  ,fhTrueMatchDist (NULL)
106  ,fhTrueMatchValue (NULL)
107  ,fhTrueMatchDistValue (NULL)
108 
109  ,fhTrueMatchNofPerHit (NULL)
110  ,fhTrueMatchNofPerPoint (NULL)
111 {
112  for ( Int_t istat = 0 ; istat < 4 ; istat++ ) {
113  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
114  fHistWidth[istat][isens]=0.;
115  for ( Int_t iregx = 0 ; iregx < 4 ; iregx++ ) {
116  for ( Int_t iregy = 0 ; iregy < 4 ; iregy++ ) {
117  fhPointToHit[istat][isens][iregx][iregy] = NULL;
118  }
119  }
120  // fhPointClosest [istat][isens] = NULL;
121  fhPointNof [istat][isens] = NULL;
122  fhPointReco [istat][isens] = NULL;
123  fhPointRecoEff [istat][isens] = NULL;
124  fhPointRadNof [istat][isens] = NULL;
125  fhPointRadReco [istat][isens] = NULL;
126  fhPointRadRecoEff [istat][isens] = NULL;
127 
128  fhPointMatch [istat][isens] = NULL;
129  fhPointMatchEff [istat][isens] = NULL;
130  fhPointRadMatch [istat][isens] = NULL;
131  fhPointRadMatchEff [istat][isens] = NULL;
132 
133  fhHitNof [istat][isens] = NULL;
134  fhHitFake [istat][isens] = NULL;
135  fhHitFakeProb [istat][isens] = NULL;
136  fhHitRadNof [istat][isens] = NULL;
137  fhHitRadFake [istat][isens] = NULL;
138  fhHitRadFakeProb [istat][isens] = NULL;
139  fhHitMultipleRate [istat][isens] = NULL;
140  fhCloseHits [istat][isens] = NULL;
141  fhTrueMatchDiXYPerSt[istat][isens] = NULL;
142  fhTrueMatchDistPerSt[istat][isens] = NULL;
143  }
144  }
145 }
146 // -------------------------------------------------------------------------
147 
148 // ----- Destructor ----------------------------------------------------
150 
151 // ----- Init -----------------------------------------------------------
153 
154  // Get and check FairRootManager
155  FairRootManager* ioman = FairRootManager::Instance();
156  if( !ioman ) {
157  cout << "-E- "<< GetName() <<"::Init: "
158  << "RootManager not instantised!" << endl;
159  return kERROR;
160  }
161 
162  // Get the pointer to the singleton FairRunAna object
163  FairRunAna* ana = FairRunAna::Instance();
164  if(NULL == ana) {
165  cout << "-E- "<< GetName() <<"::Init :"
166  <<" no FairRunAna object!" << endl;
167  return kERROR;
168  }
169  // Get the pointer to run-time data base
170  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
171  if(NULL == rtdb) {
172  cout << "-E- "<< GetName() <<"::Init :"
173  <<" no runtime database!" << endl;
174  return kERROR;
175  }
176 
177  // Get PndGemPoint (MCPoint) array
178  fMCPointArray = (TClonesArray*) ioman->GetObject("GEMPoint");
179  if( !fMCPointArray ) {
180  cout << "-E- "<< GetName() <<"::Init: No MCPoint array!"
181  << endl;
182  return kERROR;
183  }
184 
185  // Get Gem hit Array
186  fGemHitArray = (TClonesArray*) ioman->GetObject("GEMHit");
187  if ( !fGemHitArray ) {
188  cout << "-E- " << GetName() << "::Init: No PndGemHit array!" << endl;
189  return kERROR;
190  }
191 
192  // Get GEM digitisation parameter container
193  fDigiPar = (PndGemDigiPar*)(rtdb->getContainer("PndGemDetectors"));
194 
195  cout << "-I- " << fName.Data() << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
196  cout << "-I- " << fName.Data() << "::Init(). Initialization succesfull." << endl;
197 
198  TList* branchNames = FairRootManager::Instance()->GetBranchNameList();
199 
200  Int_t nofGemBranches = 0;
201  fGemPointNumber = -1;
202  for (int i = 0; i < branchNames->GetEntries(); i++){
203  TObjString* branchName = (TObjString*)branchNames->At(i);
204  if ( !branchName->GetString().BeginsWith("GEM") ) continue;
205  if ( branchName->GetString() == "GEMPoint" ) fGemPointNumber = i;
206  fGemData[nofGemBranches] = (TClonesArray*) ioman->GetObject(branchName->GetString().Data());
207  fGemDataPointer[i] = nofGemBranches;
208  cout << "GEM Data [" << nofGemBranches << "] = " << branchName->GetString().Data() << " referred to as " << i << (fGemPointNumber==i?" ***":"") << endl;
209  nofGemBranches++;
210  }
211 
212  CreateHistos();
213 
214  return kSUCCESS;
215 }
216 // -------------------------------------------------------------------------
217 
218 // ----- Private method ReInit -----------------------------------------
220 
221  return kSUCCESS;
222 }
223 // -------------------------------------------------------------------------
224 
225 // ----- Private method SetParContainers -------------------------------
227 
228  // Get run and runtime database
229  FairRunAna* run = FairRunAna::Instance();
230  if ( ! run ) Fatal("SetParContainers", "No analysis run");
231 
232  FairRuntimeDb* db = run->GetRuntimeDb();
233  if ( ! db ) Fatal("SetParContainers", "No runtime database");
234 
235  // Get GEM digitisation parameter container
236  fDigiPar = (PndGemDigiPar*)(db->getContainer("PndGemDetectors"));
237 
238  fNofEvents = 0;
239 }
240 // -------------------------------------------------------------------------
241 
242 
243 // ----- Public method Exec --------------------------------------------
244 void PndGemFindHitsQA::Exec(Option_t*) {
245 
246  Bool_t printInfo = kFALSE;
247 
248  fNofEvents++;
249 
250  Int_t nofGemHits = fGemHitArray->GetEntries();
251  Int_t nofGemPnts = fMCPointArray->GetEntries();
252 
253  PndGemHit* gemHit = NULL;
254  PndGemMCPoint* gemPnt = NULL;
255 
256  // cout << "---------------------------------" << endl;
257  // cout << "gem " << fHistWidth[istat][isens] << " wi
258  if ( printInfo ) {
259  cout << "GEM Point Array has " << nofGemPnts << " entries" << endl;
260  cout << "GEM Hit Array has " << nofGemHits << " entries" << endl;
261  }
262 
263  for ( Int_t ipnt = 0 ; ipnt < nofGemPnts ; ipnt++ ) {
264  gemPnt = (PndGemMCPoint*)fMCPointArray->At(ipnt);
265  Int_t sensorId = gemPnt->GetSensorId();
266  Int_t station = fDigiPar->GetStationNr(sensorId)-1;
267  Int_t sensor = fDigiPar->GetSensorNr (sensorId)-1;
268  Double_t pntX = (gemPnt->GetX()+gemPnt->GetXOut())/2.;
269  Double_t pntY = (gemPnt->GetY()+gemPnt->GetYOut())/2.;
270  Int_t regX = fHistPlaneDivs+TMath::Floor(pntX/fHistWidth[station][sensor]);
271  Int_t regY = fHistPlaneDivs+TMath::Floor(pntY/fHistWidth[station][sensor]);
272 
273  if ( printInfo ) {
274  cout << "GEM Point at " << gemPnt->GetX() << " , " << gemPnt->GetY() << " , " << gemPnt->GetZ() << endl;
275  cout << "----> Will go to [ " << station << " ] [ " << sensor << " ] [ " << regX << " ] [ " << regY << " ]" << endl;
276  }
277 
278  Int_t nofCloseHits = 0;
279  Double_t distToClosestHit = 100000.;
280  for ( Int_t ihit = 0 ; ihit < nofGemHits ; ihit++ ) {
281  gemHit = (PndGemHit*)fGemHitArray->At(ihit);
282  if ( gemHit->GetStationNr()!=station+1 || gemHit->GetSensorNr()!=sensor+1 ) continue;
283 
284  fhPointToHit[station][sensor][regX][regY]->Fill(pntX-gemHit->GetX(),
285  pntY-gemHit->GetY());
286 
287  if ( printInfo ) {
288  cout << "**** " << pntX << " : " << gemHit->GetX() << " *** " << pntY << " : " << gemHit->GetY() << " ***" << endl;
289  }
290 
291  Double_t p2hDistSq = (pntX-gemHit->GetX())*(pntX-gemHit->GetX())+(pntY-gemHit->GetY())*(pntY-gemHit->GetY());
292  if ( distToClosestHit > p2hDistSq )
293  distToClosestHit = p2hDistSq;
294  if ( p2hDistSq < fPointEffDist*fPointEffDist ) {
295  if ( printInfo ) {
296  cout << "POINT " << ipnt << " AND HIT " << ihit << " ARE CLOSE!!!" << endl;
297  }
298  nofCloseHits++;
299  }
300  }
301  // fhPointClosest [station][sensor]->Fill(pntX,pntY,distToClosestHit);
302  fhCloseHits [station][sensor]->Fill(nofCloseHits);
303  fhPointNof [station][sensor]->Fill(pntX,pntY);
304  fhPointRadNof [station][sensor]->Fill(TMath::Sqrt(pntX*pntX+pntY*pntY));
305  if ( nofCloseHits > 0 ) {
306  fhPointReco [station][sensor]->Fill(pntX,pntY);
307  fhPointRadReco [station][sensor]->Fill(TMath::Sqrt(pntX*pntX+pntY*pntY));
308  if ( nofCloseHits > 1 ) {
309  fhHitMultipleRate [station][sensor]->Fill(pntX,pntY);
310  }
311  }
312  else {
313  if ( fVerbose ) {
314  cout << "NO MATCHING POINT IN EVENT " << fNofEvents << " FOR POINT ( " << pntX << " , " << pntY << " , " << gemPnt->GetZ() << " ) " << "[" << gemPnt->GetX() << "X" << gemPnt->GetXOut() << "] [" << gemPnt->GetY() << "Y" << gemPnt->GetYOut() << "]" << endl;
315  }
316  }
317  }
318 
319  for ( Int_t ihit = 0 ; ihit < nofGemHits ; ihit++ ) {
320  gemHit = (PndGemHit*)fGemHitArray->At(ihit);
321 
322  Int_t station = gemHit->GetStationNr()-1;
323  Int_t sensor = gemHit->GetSensorNr()-1;
324  Double_t hitX = gemHit->GetX();
325  Double_t hitY = gemHit->GetY();
326 
327  Bool_t hitHasMatchingPoint = kFALSE;
328  // Double_t distToClosestPoint = 100000.;
329  for ( Int_t ipnt = 0 ; ipnt < nofGemPnts ; ipnt++ ) {
330  gemPnt = (PndGemMCPoint*)fMCPointArray->At(ipnt);
331  Int_t sensorId = gemPnt->GetSensorId();
332  if ( fDigiPar->GetStationNr(sensorId)!=station+1 || fDigiPar->GetSensorNr(sensorId)!=sensor+1 ) continue;
333 
334  Double_t pntX = (gemPnt->GetX()+gemPnt->GetXOut())/2.;
335  Double_t pntY = (gemPnt->GetY()+gemPnt->GetYOut())/2.;
336 
337  // Double_t p2hDistSq = (hitX-gemPnt->GetX())*(hitX-gemPnt->GetX())+(hitY-gemPnt->GetY())*(hitY-gemPnt->GetY());
338  // if ( distToClosestPoint > p2hDistSq )
339  // distToClosestPoint = p2hDistSq;
340  if ( TMath::Sqrt((hitX-pntX)*(hitX-pntX)+(hitY-pntY)*(hitY-pntY)) < fPointEffDist )
341  hitHasMatchingPoint = kTRUE;
342  }
343  fhHitNof [station][sensor]->Fill(hitX,hitY);
344  fhHitRadNof [station][sensor]->Fill(TMath::Sqrt(hitX*hitX+hitY*hitY));
345  // if ( distToClosestPoint > fPointEffDist*fPointEffDist )
346  if ( !hitHasMatchingPoint ) {
347  fhHitFake [station][sensor]->Fill(hitX,hitY);
348  fhHitRadFake [station][sensor]->Fill(TMath::Sqrt(hitX*hitX+hitY*hitY));
349  }
350  }
351 
352  std::vector<std::pair<Int_t,Int_t> > mcMatchPointHit;
353 
354  Bool_t printMCMatching = kFALSE;
355  if ( fVerbose >= 1 ) printMCMatching = kTRUE;
356 
357  for ( Int_t ihit = 0 ; ihit < nofGemHits ; ihit++ ) {
358  gemHit = (PndGemHit*)fGemHitArray->At(ihit);
359  if ( printMCMatching ) {
360  cout << "---> hit " << ihit << " has " << gemHit->GetNLinks() << " links" << endl;
361  for ( Int_t ilink = 0 ; ilink < gemHit->GetNLinks() ; ilink++ ) {
362  cout << " " << ilink << " > " << gemHit->GetLink(ilink).GetType() << " . " << gemHit->GetLink(ilink).GetIndex() << endl;
363  }
364  }
365 
366  Int_t quickPointIndex = -1;
367  for ( Int_t ilink = 0 ; ilink < gemHit->GetNLinks() ; ilink++) {
368  if ( gemHit->GetLink(ilink).GetType() == fGemPointNumber ) {
369  if ( quickPointIndex != -1 ) // already have found different point index - so it is a ghost
370  { // but what about hits from clusters, if 98% one point, and only 2% other point???
371  quickPointIndex = -2;
372  break;
373  }
374  quickPointIndex = gemHit->GetLink(ilink).GetIndex();
375  }
376  }
377  if ( quickPointIndex >= 0 ) {
378  pair<Int_t, Int_t> a (quickPointIndex,ihit);
379  mcMatchPointHit.push_back(a);
380  }
381 
382  if ( quickPointIndex != -1 ) continue; // do not have any links to points
383  if ( gemHit->GetNLinks() == 2 ) {
384  Int_t maxPnt0 = -1, maxPnt1 = -1;
385  std::vector<Int_t> pointVector0;
386  std::vector<Int_t> pointVector1;
387  GetPointVector(gemHit->GetLink(0).GetType(),gemHit->GetLink(0).GetIndex(),pointVector0,printMCMatching);
388  GetPointVector(gemHit->GetLink(1).GetType(),gemHit->GetLink(1).GetIndex(),pointVector1,printMCMatching);
389  if ( printMCMatching )
390  cout << "VECT0: (" << gemHit->GetLink(0).GetIndex() << ") " << flush;
391  for ( size_t ipnt = 0 ; ipnt < pointVector0.size() ; ipnt++ ) {
392  if ( printMCMatching )
393  cout << pointVector0[ipnt] << " . " << flush;
394  if ( maxPnt0 < pointVector0[ipnt] ) {
395  maxPnt0 = pointVector0[ipnt];
396  }
397  }
398  if ( printMCMatching )
399  cout << "\b\b" << endl;
400  if ( printMCMatching )
401  cout << "VECT1: (" << gemHit->GetLink(1).GetIndex() << ") " << flush;
402  for ( size_t ipnt = 0 ; ipnt < pointVector1.size() ; ipnt++ ) {
403  if ( printMCMatching )
404  cout << pointVector1[ipnt] << " . " << flush;
405  if ( maxPnt1 < pointVector1[ipnt] ) {
406  maxPnt1 = pointVector1[ipnt];
407  }
408  }
409  if ( printMCMatching )
410  cout << "\b\b" << endl;
411  if ( printMCMatching )
412  cout << "highest points are " << maxPnt0 << " , " << maxPnt1 << endl;
413  std::vector<Int_t> countPointV0(maxPnt0+1,0);
414  std::vector<Int_t> countPointV1(maxPnt1+1,0);
415  for ( size_t ipnt = 0 ; ipnt < pointVector0.size() ; ipnt++ ) {
416  ++countPointV0[pointVector0[ipnt]];
417  }
418  for ( size_t ipnt = 0 ; ipnt < pointVector1.size() ; ipnt++ ) {
419  ++countPointV1[pointVector1[ipnt]];
420  }
421  if ( maxPnt0 > maxPnt1 )
422  maxPnt0 = maxPnt1;
423  for ( Int_t ipnt = 0 ; ipnt < maxPnt0+1 ; ipnt++ ) {
424  if ( printMCMatching )
425  cout << ipnt << " - " << countPointV0[ipnt] << " ? " << countPointV1[ipnt] << endl;
426  if ( countPointV0[ipnt] > 0 ) {
427  if ( countPointV1[ipnt] > 0 ) {
428  Double_t tempMatch = 100.*Double_t(countPointV0[ipnt]*countPointV1[ipnt])/(Double_t(pointVector0.size()*pointVector1.size()));
429  if ( printMCMatching )
430  cout << "HIT " << ihit << " IS " << countPointV0[ipnt]*countPointV1[ipnt] << " / " << pointVector0.size()*pointVector1.size()
431  << " ( " << tempMatch
432  << " % ) POINT " << ipnt << endl;
433  pair<Int_t, Int_t> a (ipnt,ihit);
434  mcMatchPointHit.push_back(a);
435  }
436  }
437  }
438  }
439  }
440 
441  std::vector<std::pair<Int_t,Int_t> >::iterator iter;
442  for ( iter = mcMatchPointHit.begin() ; iter != mcMatchPointHit.end() ; iter++ ) {
443  gemPnt = (PndGemMCPoint*)fMCPointArray->At(iter->first);
444  gemHit = (PndGemHit*) fGemHitArray ->At(iter->second);
445 
446  Double_t pntX = (gemPnt->GetX()+gemPnt->GetXOut())/2.;
447  Double_t pntY = (gemPnt->GetY()+gemPnt->GetYOut())/2.;
448 
449  Int_t station = gemHit->GetStationNr()-1;
450  Int_t sensor = gemHit->GetSensorNr()-1;
451  fhTrueMatchDiXYPerSt[station][sensor]->Fill(pntX-gemHit->GetX(),
452  pntY-gemHit->GetY());
453  fhTrueMatchDiXY ->Fill(pntX-gemHit->GetX(),
454  pntY-gemHit->GetY());
455  Double_t tempDist = TMath::Sqrt((pntX-gemHit->GetX())*(pntX-gemHit->GetX())+
456  (pntY-gemHit->GetY())*(pntY-gemHit->GetY()));
457  if ( tempDist > 5. ) {
458  cout << "Event " << fNofEvents << ": point " << iter->first
459  << " at ( " << gemPnt->GetX() << "," << gemPnt->GetY() << "," << gemPnt->GetZ() << " )"
460  << " to ( " << gemPnt->GetXOut() << "," << gemPnt->GetYOut() << "," << gemPnt->GetZOut() << " )" << endl
461  << " hit " << iter->second << " at ( " << gemHit->GetX() << "," << gemHit->GetY() << "," << gemHit->GetZ() << " )" << endl
462  << " ---> dist " << tempDist << endl;//" with match value = " << tempMatch << endl;
463  }
464  fhTrueMatchDistPerSt[station][sensor]->Fill(tempDist);
465  fhTrueMatchDist ->Fill(tempDist);
466  // fhTrueMatchValue ->Fill(tempMatch);
467  // fhTrueMatchDistValue ->Fill(tempMatch,tempDist);
468  }
469 
470  if ( printMCMatching )
471  cout << "True MC Matches betwen points and hits: " << mcMatchPointHit.size() << endl;
472 
473  std::vector<Int_t> nofMatchesPerHit (nofGemHits,0);
474  std::vector<Int_t> nofMatchesPerPoint(nofGemPnts,0);
475  for ( iter = mcMatchPointHit.begin() ; iter != mcMatchPointHit.end() ; iter++ ) {
476  ++nofMatchesPerPoint[iter->first];
477  ++nofMatchesPerHit [iter->second];
478  }
479  for ( Int_t ihit = 0 ; ihit < nofGemHits ; ihit++ ) {
480  fhTrueMatchNofPerHit->Fill(nofMatchesPerHit[ihit]);
481  }
482  for ( Int_t ipnt = 0 ; ipnt < nofGemPnts ; ipnt++ ) {
483  fhTrueMatchNofPerPoint->Fill(nofMatchesPerPoint[ipnt]);
484  gemPnt = (PndGemMCPoint*)fMCPointArray->At(ipnt);
485  if ( nofMatchesPerPoint[ipnt] > 0 ) {
486  Double_t pntX = (gemPnt->GetX()+gemPnt->GetXOut())/2.;
487  Double_t pntY = (gemPnt->GetY()+gemPnt->GetYOut())/2.;
488  Int_t sensorId = gemPnt->GetSensorId();
489  Int_t station = fDigiPar->GetStationNr(sensorId)-1;
490  Int_t sensor = fDigiPar->GetSensorNr (sensorId)-1;
491  fhPointMatch [station][sensor]->Fill(pntX,pntY);
492  fhPointRadMatch [station][sensor]->Fill(TMath::Sqrt(pntX*pntX+pntY*pntY));
493  }
494 
495  /*
496  if ( nofMatchesPerPoint[ipnt] == 0 ) {
497  // cout << "POINT " << ipnt << " AT "
498  // << gemPnt->GetX() << "-" << gemPnt->GetXOut() << " "
499  // << gemPnt->GetY() << "-" << gemPnt->GetYOut() << " "
500  // << gemPnt->GetZ() << "-" << gemPnt->GetZOut() << " WAS NOT RECONSTRUCTED" << endl;
501  for ( Int_t ihit = 0 ; ihit < nofGemHits ; ihit++ ) {
502  gemHit = (PndGemHit*)fGemHitArray->At(ihit);
503 
504  for ( Int_t ilink = 0 ; ilink < gemHit->GetNLinks() ; ilink++ ) {
505  std::vector<Int_t> pointVector;
506  Bool_t pointFoundHere = kFALSE;
507  GetPointVector(gemHit->GetLink(ilink).GetType(),gemHit->GetLink(ilink).GetIndex(),pointVector,printMCMatching);
508  for ( Int_t imatch = 0 ; imatch < pointVector.size() ; imatch++ ) {
509  if ( pointVector[imatch] == ipnt ) {
510  pointFoundHere = kTRUE;
511  }
512  }
513  // if ( pointFoundHere ) cout << "FOUND POINT " << ipnt << " IN HIT " << ihit << " LINK " << ilink << endl;
514  }
515  }
516  }
517  */
518  }
519 }
520 // ------------------------------------------------------------
521 
522 // ----- Private method GetPointVector, recurrence function --------------------------------------------
523 Int_t PndGemFindHitsQA::GetPointVector(Int_t arrayId, Int_t entryId, std::vector<Int_t>& pointVector, Bool_t printInfo) {
524  if ( printInfo ) {
525  for ( Int_t itemp = 0 ; itemp < 10 - arrayId ; itemp++ ) cout << " " << flush;
526  cout << "GetPointVector(" << arrayId << ", " << entryId << ")" << endl;
527  }
528  if ( arrayId == fGemPointNumber ) {
529  pointVector.push_back(entryId);
530  return pointVector.size();
531  }
532  if ( arrayId == 0 ) { // this is probably MCTrack!
533  return 0;
534  }
535 
536  FairMultiLinkedData* tempData = (FairMultiLinkedData*)fGemData[fGemDataPointer[arrayId]]->At(entryId);
537  for ( Int_t ilink = 0 ; ilink < tempData->GetNLinks() ; ilink++ ) {
538  GetPointVector(tempData->GetLink(ilink).GetType(),tempData->GetLink(ilink).GetIndex(),pointVector,printInfo);
539  }
540  return pointVector.size();
541 }
542 // ------------------------------------------------------------
543 
544 // ----- Private method CreateHistos --------------------------------------------
546  fHistoList = new TList();
547  // number of mc tracks, reco tracks, efficiency as function of MOMENTUM
548 
549  Int_t nStations = fDigiPar->GetNStations();
550 
551  fHistPlaneDivs = 2;
552 
553  for ( Int_t istat = 0 ; istat < nStations ; istat++ ) {
554  PndGemStation* station = (PndGemStation*)fDigiPar->GetStation(istat);
555  Int_t nSensors = station->GetNSensors();
556  for ( Int_t isens = 0 ; isens < nSensors ; isens++ ) {
557  PndGemSensor* sensor = (PndGemSensor*)station->GetSensor(isens);
558  Double_t sensOutRad = sensor->GetOuterRadius();
559  fHistWidth[istat][isens] = TMath::Ceil(sensOutRad/fHistPlaneDivs); // make it an even number
560  Double_t histBeginX = -fHistWidth[istat][isens]*fHistPlaneDivs;
561  for ( Int_t iregx = 0 ; iregx < 2*fHistPlaneDivs ; iregx++ ) {
562  Double_t histBeginY = -fHistWidth[istat][isens]*fHistPlaneDivs;
563  for ( Int_t iregy = 0 ; iregy < 2*fHistPlaneDivs ; iregy++ ) {
564  fhPointToHit[istat][isens][iregx][iregy] = new TH2F(Form("fhPointToHit_s%d_s%d_x%d_y%d",istat,isens,iregx,iregy),
565  Form("Points vs hits, station %d, sensor %d at z=%.1fcm, %.1f<x<%.1f, %.1f<y<%.1f",
566  istat,isens,sensor->GetZ0(),
567  histBeginX,histBeginX+fHistWidth[istat][isens],
568  histBeginY,histBeginY+fHistWidth[istat][isens]),
569  2000,-1.,1.,
570  2000,-1.,1.);
571  histBeginY+=fHistWidth[istat][isens];
572  fHistoList->Add(fhPointToHit[istat][isens][iregx][iregy]);
573  }
574  histBeginX+=fHistWidth[istat][isens];
575  }
576  // fhPointClosest [istat][isens] = new TH3F(Form("fhPointClosest_s%d_s%d",istat,isens),
577  // Form("Distance from point to closest hits, station %d, sensor %d",istat,isens),
578  // 2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
579  // 2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
580  // 1000,0.,1.);
581  fhPointNof [istat][isens] = new TH2F(Form("fhPointNof_s%d_s%d",istat,isens),
582  Form("Number of points, station %d, sensor %d",istat,isens),
583  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
584  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
585  fhPointReco [istat][isens] = new TH2F(Form("fhPointReco_s%d_s%d",istat,isens),
586  Form("Number of reconstructed points (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
587  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
588  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
589  fhPointRecoEff [istat][isens] = new TH2F(Form("fhPointRecoEff_s%d_s%d",istat,isens),
590  Form("Hit finding efficiency (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
591  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
592  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
593  fhPointRadNof [istat][isens] = new TH1F(Form("fhPointRadNof_s%d_s%d",istat,isens),
594  Form("Number of points, station %d, sensor %d",istat,isens),
595  1000,0,100);
596  fhPointRadReco [istat][isens] = new TH1F(Form("fhPointRadReco_s%d_s%d",istat,isens),
597  Form("Number of reconstructed points (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
598  1000,0,100);
599  fhPointRadRecoEff [istat][isens] = new TH1F(Form("fhPointRadRecoEff_s%d_s%d",istat,isens),
600  Form("Hit finding efficiency (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
601  1000,0,100);
602 
603  fhPointMatch [istat][isens] = new TH2F(Form("fhPointMatch_s%d_s%d",istat,isens),
604  Form("Number of matched points (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
605  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
606  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
607  fhPointMatchEff [istat][isens] = new TH2F(Form("fhPointMatchEff_s%d_s%d",istat,isens),
608  Form("Matched points efficiency (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
609  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
610  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
611  fhPointRadMatch [istat][isens] = new TH1F(Form("fhPointRadMatch_s%d_s%d",istat,isens),
612  Form("Number of matched points (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
613  1000,0,100);
614  fhPointRadMatchEff [istat][isens] = new TH1F(Form("fhPointRadMatchEff_s%d_s%d",istat,isens),
615  Form("Matched points efficiency (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
616  1000,0,100);
617 
618  fhHitNof [istat][isens] = new TH2F(Form("fhHitNof_s%d_s%d",istat,isens),
619  Form("Number of all hits (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
620  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
621  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
622  fhHitFake [istat][isens] = new TH2F(Form("fhHitFake_s%d_s%d",istat,isens),
623  Form("Number of fake hits (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
624  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
625  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
626  fhHitFakeProb [istat][isens] = new TH2F(Form("fhHitFakeProb_s%d_s%d",istat,isens),
627  Form("Hit fake probability (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
628  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
629  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
630  fhHitRadNof [istat][isens] = new TH1F(Form("fhHitRadNof_s%d_s%d",istat,isens),
631  Form("Number of all hits (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
632  1000,0,100);
633  fhHitRadFake [istat][isens] = new TH1F(Form("fhHitRadFake_s%d_s%d",istat,isens),
634  Form("Number of fake hits (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
635  1000,0,100);
636  fhHitRadFakeProb [istat][isens] = new TH1F(Form("fhHitRadFakeProb_s%d_s%d",istat,isens),
637  Form("Hit fake probability (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
638  1000,0,100);
639 
640  fhHitMultipleRate [istat][isens] = new TH2F(Form("fhHitMultipleRate_s%d_s%d",istat,isens),
641  Form("Distance from point to closest hits (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
642  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad),
643  2*TMath::Ceil(sensOutRad),-TMath::Ceil(sensOutRad),TMath::Ceil(sensOutRad));
644  fhCloseHits [istat][isens] = new TH1F(Form("fhCloseHits_s%d_s%d",istat,isens),
645  Form("Number of close hits per point (point-hit %.2f cm), station %d, sensor %d",fPointEffDist,istat,isens),
646  101,-0.5,100.5);
647  fhTrueMatchDiXYPerSt[istat][isens] = new TH2F(Form("fhTrueMatchDiXYPerSt_s%d_s%d",istat,isens),
648  Form("True matches: difference in XY between points and hits, station %d, sensor %d",istat,isens),
649  200,-1.,1.,
650  200,-1.,1.);
651  fhTrueMatchDistPerSt[istat][isens] = new TH1F(Form("fhTrueMatchDistPerSt_s%d_s%d",istat,isens),
652  Form("True matches: distance between points and hits, station %d, sensor %d",istat,isens),
653  1000,0,10);
654  // fHistoList->Add(fhPointClosest [istat][isens]);
655  fHistoList->Add(fhPointNof [istat][isens]);
656  fHistoList->Add(fhPointReco [istat][isens]);
657  fHistoList->Add(fhPointRecoEff [istat][isens]);
658  fHistoList->Add(fhPointRadNof [istat][isens]);
659  fHistoList->Add(fhPointRadReco [istat][isens]);
660  fHistoList->Add(fhPointRadRecoEff [istat][isens]);
661  fHistoList->Add(fhPointMatch [istat][isens]);
662  fHistoList->Add(fhPointMatchEff [istat][isens]);
663  fHistoList->Add(fhPointRadMatch [istat][isens]);
664  fHistoList->Add(fhPointRadMatchEff [istat][isens]);
665  fHistoList->Add(fhHitNof [istat][isens]);
666  fHistoList->Add(fhHitFake [istat][isens]);
667  fHistoList->Add(fhHitFakeProb [istat][isens]);
668  fHistoList->Add(fhHitRadNof [istat][isens]);
669  fHistoList->Add(fhHitRadFake [istat][isens]);
670  fHistoList->Add(fhHitRadFakeProb [istat][isens]);
671  fHistoList->Add(fhHitMultipleRate [istat][isens]);
672  fHistoList->Add(fhCloseHits [istat][isens]);
673  fHistoList->Add(fhTrueMatchDiXYPerSt[istat][isens]);
674  fHistoList->Add(fhTrueMatchDistPerSt[istat][isens]);
675  }
676  }
677  fhTrueMatchDiXY = new TH2F("fhTrueMatchDiXY",
678  "True matches: difference in XY between points and hits",
679  2000,-1.,1.,
680  2000,-1.,1.);
681  fhTrueMatchDist = new TH1F("fhTrueMatchDist",
682  "True matches: distance between points and hits",
683  1000,0,10);
684  fhTrueMatchValue = new TH1F("fhTrueMatchValue",
685  "True matches: hit-point match value in pecent",
686  1010,0.,101.);
687  fhTrueMatchDistValue = new TH2F("fhTrueMatchDistValue",
688  "True matches: hit-point match value vs distance",
689  101,0.,101.,
690  200,0.,2.);
691 
692  fhTrueMatchNofPerHit = new TH1F("fhTrueMatchNofPerHit",
693  "True matches: number of matched points per hit",
694  101,-0.5,100.5);
695  fhTrueMatchNofPerPoint = new TH1F("fhTrueMatchNofPerPoint",
696  "True matches: number of matched hits per point",
697  101,-0.5,100.5);
698 
703 
706 }
707 // ------------------------------------------------------------
708 
709 // ----- Private method DivideHistos --------------------------------------------
711  hist1->Sumw2();
712  hist2->Sumw2();
713  hist3->Divide(hist1,hist2,1,1,"B");
714 }
715 // ------------------------------------------------------------
716 
717 // ----- Private method Finish --------------------------------------------
719 
720  cout << "-------------------- PndGemFindHitsQA : Finish -------------------" << endl;
721  // cout << " dividing histos" << endl;
722 
723  Int_t nStations = fDigiPar->GetNStations();
724  for ( Int_t istat = 0 ; istat < nStations ; istat++ ) {
725  for ( Int_t isens = 0 ; isens < 2 ; isens++ ) {
726  // cout << "doing " << istat << " / " << isens << endl;
727  fhPointReco [istat][isens]->Sumw2();
728  fhPointNof [istat][isens]->Sumw2();
729  fhPointRecoEff [istat][isens]->Divide(fhPointReco[istat][isens],fhPointNof[istat][isens],1.,1.,"B");
730  fhPointRecoEff [istat][isens]->Scale(100.);
731 
732  fhPointMatch [istat][isens]->Sumw2();
733  fhPointMatchEff[istat][isens]->Divide(fhPointMatch[istat][isens],fhPointNof[istat][isens],1.,1.,"B");
734  fhPointMatchEff[istat][isens]->Scale(100.);
735 
736  fhPointRadReco [istat][isens]->Sumw2();
737  fhPointRadNof [istat][isens]->Sumw2();
738  fhPointRadRecoEff [istat][isens]->Divide(fhPointRadReco[istat][isens],fhPointRadNof[istat][isens],1.,1.,"B");
739  fhPointRadRecoEff [istat][isens]->Scale(100.);
740 
741  fhPointRadMatch [istat][isens]->Sumw2();
742  fhPointRadMatchEff[istat][isens]->Divide(fhPointRadMatch[istat][isens],fhPointRadNof[istat][isens],1.,1.,"B");
743  fhPointRadMatchEff[istat][isens]->Scale(100.);
744 
745  fhHitFake [istat][isens]->Sumw2();
746  fhHitNof [istat][isens]->Sumw2();
747  fhHitFakeProb [istat][isens]->Divide(fhHitFake [istat][isens],fhHitNof [istat][isens]);
748  fhHitFakeProb [istat][isens]->Scale(100.);
749 
750  fhHitRadFake [istat][isens]->Sumw2();
751  fhHitRadNof [istat][isens]->Sumw2();
752  fhHitRadFakeProb [istat][isens]->Divide(fhHitRadFake [istat][isens],fhHitRadNof [istat][isens]);
753  fhHitRadFakeProb [istat][isens]->Scale(100.);
754 
755  }
756  }
757 
758  // cout << "-------------------- PndGemFindHitsQA : Summary ------------------" << endl;
759 
760  cout << " Events: " << setw(10) << fNofEvents << endl;
761 
762  cout << "-------------------- PndGemFindHitsQA : Efficiency ---------------" << endl;
763  cout << " efficiency = number of points that have hit closer than " << fPointEffDist << " cm / number of points" << endl;
764  Int_t nofPointsAll = 0;
765  Int_t nofPointsReco = 0;
766  for ( Int_t istat = 0 ; istat < nStations ; istat++ ) {
767  PndGemStation* station = (PndGemStation*)fDigiPar->GetStation(istat);
768  Int_t nSensors = station->GetNSensors();
769  for ( Int_t isens = 0 ; isens < nSensors ; isens++ ) {
770  nofPointsAll += fhPointRadNof [istat][isens]->GetEntries();
771  nofPointsReco += fhPointRadReco [istat][isens]->GetEntries();
772  cout << "Station " << istat << " sensor " << isens << ", eff = "
773  << ((Double_t)(fhPointRadReco[istat][isens]->GetEntries()))/((Double_t)(fhPointRadNof [istat][isens]->GetEntries()))*100.
774  << "% (" << fhPointRadReco [istat][isens]->GetEntries()
775  << " / " << fhPointRadNof [istat][isens]->GetEntries()
776  << ")" << endl;
777  }
778  }
779  cout << " OVERALL EFF = " << ((Double_t)(nofPointsReco))/((Double_t)(nofPointsAll))*100.
780  << "% (" << nofPointsReco
781  << " / " << nofPointsAll
782  << ")" << endl;
783  cout << "o o o o o o o o o o o o o o o o o" << endl;
784  Int_t pointsAll = fhTrueMatchNofPerPoint->GetEntries();
785  Int_t pointsMCM = pointsAll-fhTrueMatchNofPerPoint->GetBinContent(fhTrueMatchNofPerPoint->FindBin(0.));
786  cout << " TRUE MATCH EFF = " << ((Double_t)(pointsMCM))/((Double_t)(pointsAll))*100. << "% (" << pointsMCM << " / " << pointsAll << ")" << endl;
787 
788  cout << "-----------------------------------------------------------------" << endl;
789 
790  TFile* temp = gFile;
791  FairRootManager* ioman = FairRootManager::Instance();
792  gFile = ioman->GetOutFile();
793  gDirectory = (TDirectory*)gFile;
794 
795  gDirectory->mkdir("GemFindHitsQA");
796  gDirectory->cd("GemFindHitsQA");
797  TIter next(fHistoList);
798  while ( TH1* histo = ((TH1*)next()) ) histo->Write();
799  gDirectory->cd("..");
800 
801  gFile = temp;
802 
803 }
804 // ------------------------------------------------------------
805 
806 
TH2F * fhTrueMatchDiXYPerSt[4][2]
Int_t fGemDataPointer[1000]
virtual InitStatus Init()
int fVerbose
Definition: poormantracks.C:24
TH2F * fhHitMultipleRate[4][2]
TH2F * fhPointNof[4][2]
Int_t run
Definition: autocutx.C:47
TClonesArray * fGemData[10]
Int_t GetPointVector(Int_t arrayId, Int_t entryId, std::vector< Int_t > &pointVector, Bool_t printInfo=kFALSE)
TClonesArray * fMCPointArray
TH1F * fhTrueMatchDistPerSt[4][2]
Int_t i
Definition: run_full.C:25
TH2F * fhPointToHit[4][2][4][4]
TH2F * fhHitNof[4][2]
TH2F * fhHitFake[4][2]
PndTransMap * map
Definition: sim_emc_apd.C:99
Int_t fNofEvents
event counter
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
virtual void Finish()
TClonesArray * fGemHitArray
Double_t GetZOut() const
Definition: PndGemMCPoint.h:85
virtual void SetParContainers()
TGeoVolume * sensor
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
Int_t GetSensorNr(Int_t sensorId)
Definition: PndGemDigiPar.h:56
TH2F * fhPointReco[4][2]
PndGemSensor * GetSensor(Int_t iSensor)
Definition: PndGemStation.h:63
TH1D * hist3
Definition: hist-t7.C:84
TH1F * fhHitRadFake[4][2]
virtual void Exec(Option_t *opt)
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
TH2F * fhPointMatchEff[4][2]
TH1F * fhPointRadReco[4][2]
Int_t a
Definition: anaLmdDigi.C:126
Int_t GetNSensors() const
Definition: PndGemStation.h:60
Int_t GetSensorId() const
Definition: PndGemMCPoint.h:90
Double_t fHistWidth[4][2]
TH1F * fhHitRadNof[4][2]
TH2F * fhHitFakeProb[4][2]
Double_t
PndGemDigiPar * fDigiPar
TH1D * hist2
Definition: hist-t7.C:82
TH1F * fhPointRadNof[4][2]
TH1F * fhCloseHits[4][2]
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
virtual InitStatus ReInit()
TH2F * fhPointMatch[4][2]
TH1F * fhPointRadMatchEff[4][2]
TH1F * fhPointRadRecoEff[4][2]
track finding quality assesment task
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
Int_t GetStationNr() const
Definition: PndGemHit.h:81
Int_t GetStationNr(Int_t sensorId)
Definition: PndGemDigiPar.h:54
TH2F * fhPointRecoEff[4][2]
void DivideHistos(TH1 *hist1, TH1 *hist2, TH1 *hist3)
cout<<"will loop over "<< t-> GetEntries()
Definition: root2ascii.C:17
ClassImp(PndAnaContFact)
Int_t iVerbose
Double_t GetZ0() const
Definition: PndGemSensor.h:100
Double_t GetOuterRadius() const
Definition: PndGemSensor.h:103
static int next[96]
Definition: ranlxd.cxx:374
Double_t GetYOut() const
Definition: PndGemMCPoint.h:84
Double_t GetXOut() const
Definition: PndGemMCPoint.h:83
cout<<"the Event No is "<< i<< endl;{{if(hit_array->GetEntriesFast()!=mc_array->GetEntriesFast()) continue;PndSdsHit *hit=(PndSdsHit *) hit_array-> At(j)
Definition: anaLmdCluster.C:71
TH1F * fhHitRadFakeProb[4][2]
PndGemStation * GetStation(Int_t iStation)
TH1F * fhPointRadMatch[4][2]
TH1D * hist1
Definition: hist-t7.C:80