BmnRoot
Loading...
Searching...
No Matches
CbmStsReconstructionQa.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- CbmStsReconstructionQa source file -----
3// ----- Created 06/02/07 by R. Karabowicz -----
4// -------------------------------------------------------------------------
5#include "CbmStsReconstructionQa.h"
6
7#include "CbmStsHit.h"
8#include "CbmStsDigi.h"
9#include "CbmStsPoint.h"
10#include "CbmStsTrack.h"
11#include "CbmTrackMatch.h"
12#include "CbmGeoStsPar.h"
13
14// Includes from base
15#include "FairGeoNode.h"
16#include "FairGeoPassivePar.h"
17#include "FairGeoVector.h"
18#include "CbmMCTrack.h"
19#include "FairRootManager.h"
20#include "FairRunAna.h"
21#include "FairRuntimeDb.h"
22
23// Includes from ROOT
24#include "TCanvas.h"
25#include "TLine.h"
26#include "TPad.h"
27#include "TLegend.h"
28#include "TPaveText.h"
29#include "TClonesArray.h"
30#include "TFile.h"
31#include "TF1.h"
32#include "TH1F.h"
33#include "TH2F.h"
34#include "TH3F.h"
35#include "TList.h"
36#include "TVector3.h"
37
38#include <iostream>
39#include <iomanip>
40using namespace std;
41
42
43// ----- Default constructor -------------------------------------------
44CbmStsReconstructionQa::CbmStsReconstructionQa(Int_t iVerbose)
45: FairTask("STSReconstructionQA", iVerbose),
46 fHitMap(),
47 fHitTrackMap(),
48 fMatchMap(),
49 fQualiMap(),
50 fMCTracks(NULL), // MCtrack
51 fStsPoints(NULL), // StsPoints
52 fStsHits(NULL), // StsHits
53 fStsTracks(NULL), // StsTrack
54 fMatches(NULL), // StsTrackMatch
55 fStsDigis(NULL), // StsDigi
56 fPassGeo(),
57 fStsGeo(),
58 fTargetPos(),
59 fNStations(0),
60 fMinHits(4),
61 fQuota(0.7),
62 fhMomAccAll(),
63 fhMomRecAll(),
64 fhMomEffAll(),
65 fhMomAccPrim(),
66 fhMomRecPrim(),
67 fhMomEffPrim(),
68 fhMomAccSec(),
69 fhMomRecSec(),
70 fhMomEffSec(),
71 fhNpAccAll(),
72 fhNpRecAll(),
73 fhNpEffAll(),
74 fhNpAccPrim(),
75 fhNpRecPrim(),
76 fhNpEffPrim(),
77 fhNpAccSec(),
78 fhNpRecSec(),
79 fhNpEffSec(),
80 fhZAccSec(),
81 fhZRecSec(),
82 fhZEffSec(),
83 fhNhClones(),
84 fhNhGhosts(),
85 fhMomClones(),
86 fhMomGhosts(),
87 fhMomResAll(),
88 fhMomResPrim(),
89 fhMomResSec(),
90 fhLowBand(),
91 fhHigBand(),
92 fhPrimaryVertex(),
93 fhRefTracks(),
94 fhRecRefTracks(),
95 fhStsTrackFMom(),
96 fhStsTrackLMom(),
97 fhStsTrackChiSq(),
98 fHistoList(),
99 fOccupHList(),
100 fNAccAll(0),
101 fNAccPrim(0),
102 fNAccRef(0),
103 fNAccSec(0),
104 fNRecAll(0),
105 fNRecPrim(0),
106 fNRecRef(0),
107 fNRecSec(0),
108 fNGhosts(0),
109 fNClones(0),
110 fNStsTracks(0),
111 fNEvents(0),
112 fNEventsFailed(0),
113 fTime(0.),
114 fOnlineAnalysis(kFALSE),
115 fOnlineCanvas(),
116 fShowStation1(2),
117 fShowStation2(5),
118 fTimer()
119{}
120
121// ----- Standard constructor ------------------------------------------
122CbmStsReconstructionQa::CbmStsReconstructionQa(Bool_t visualizeBool, Int_t minHits, Double_t quota, Int_t iVerbose)
123: FairTask("STSReconstructionQA", iVerbose),
124 fHitMap(),
125 fHitTrackMap(),
126 fMatchMap(),
127 fQualiMap(),
128 fMCTracks(NULL), // MCtrack
129 fStsPoints(NULL), // StsPoints
130 fStsHits(NULL), // StsHits
131 fStsTracks(NULL), // StsTrack
132 fMatches(NULL), // StsTrackMatch
133 fStsDigis(NULL), // StsDigi
134 fPassGeo(),
135 fStsGeo(),
136 fTargetPos(),
137 fNStations(0),
138 fMinHits(minHits),
139 fQuota(quota),
140 fhMomAccAll(),
141 fhMomRecAll(),
142 fhMomEffAll(),
143 fhMomAccPrim(),
144 fhMomRecPrim(),
145 fhMomEffPrim(),
146 fhMomAccSec(),
147 fhMomRecSec(),
148 fhMomEffSec(),
149 fhNpAccAll(),
150 fhNpRecAll(),
151 fhNpEffAll(),
152 fhNpAccPrim(),
153 fhNpRecPrim(),
154 fhNpEffPrim(),
155 fhNpAccSec(),
156 fhNpRecSec(),
157 fhNpEffSec(),
158 fhZAccSec(),
159 fhZRecSec(),
160 fhZEffSec(),
161 fhNhClones(),
162 fhNhGhosts(),
163 fhMomClones(),
164 fhMomGhosts(),
165 fhMomResAll(),
166 fhMomResPrim(),
167 fhMomResSec(),
168 fhLowBand(),
169 fhHigBand(),
170 fhPrimaryVertex(),
171 fhRefTracks(),
172 fhRecRefTracks(),
173 fhStsTrackFMom(),
174 fhStsTrackLMom(),
175 fhStsTrackChiSq(),
176 fHistoList(),
177 fOccupHList(),
178 fNAccAll(0),
179 fNAccPrim(0),
180 fNAccRef(0),
181 fNAccSec(0),
182 fNRecAll(0),
183 fNRecPrim(0),
184 fNRecRef(0),
185 fNRecSec(0),
186 fNGhosts(0),
187 fNClones(0),
188 fNStsTracks(0),
189 fNEvents(0),
190 fNEventsFailed(0),
191 fTime(0.),
192 fOnlineAnalysis(visualizeBool),
193 fOnlineCanvas(),
194 fShowStation1(2),
195 fShowStation2(5),
196 fTimer()
197{
198 fPartPdgTable[0] = 11; // electron
199 fPartPdgTable[1] = - 11; // positron
200 fPartPdgTable[2] = 211; // pi +
201 fPartPdgTable[3] = - 211; // pi -
202 fPartPdgTable[4] = 321; // kaon +
203 fPartPdgTable[5] = - 321; // kaon -
204 fPartPdgTable[6] = 2212; // proton
205 fPartPdgTable[7] = -2212; // antiproton
206 fPartPdgTable[8] = -7777; // don't use
207 fPartPdgTable[9] = -7777; // don't use
208
209}
210
211// ----- Destructor ----------------------------------------------------
212CbmStsReconstructionQa::~CbmStsReconstructionQa()
213{
214 fHistoList->Delete();
215 delete fHistoList;
216 fOccupHList->Delete();
217 delete fOccupHList;
218}
219
220// ----- Public method SetParContainers --------------------------------
221void CbmStsReconstructionQa::SetParContainers() {
222
223 // Get Run
224 FairRunAna* run = FairRunAna::Instance();
225 if ( ! run ) {
226 cout << "-E- " << GetName() << "::SetParContainers: No FairRunAna!"
227 << endl;
228 return;
229 }
230
231 // Get Runtime Database
232 FairRuntimeDb* runDb = run->GetRuntimeDb();
233 if ( ! run ) {
234 cout << "-E- " << GetName() << "::SetParContainers: No runtime database!"
235 << endl;
236 return;
237 }
238
239 // Get passive geometry parameters
240 fPassGeo = (FairGeoPassivePar*) runDb->getContainer("FairGeoPassivePar");
241 if ( ! fPassGeo ) {
242 cout << "-E- " << GetName() << "::SetParContainers: "
243 << "No passive geometry parameters!" << endl;
244 return;
245 }
246
247 // Get STS geometry parameters
248 fStsGeo = (CbmGeoStsPar*) runDb->getContainer("CbmGeoStsPar");
249 if ( ! fStsGeo ) {
250 cout << "-E- " << GetName() << "::SetParContainers: "
251 << "No STS geometry parameters!" << endl;
252 return;
253 }
254
255}
256// -------------------------------------------------------------------------
257
258
259// ----- Public method Init --------------------------------------------
260InitStatus CbmStsReconstructionQa::Init() {
261
262 cout << "==========================================================="
263 << endl;;
264 cout << GetName() << ": Initialising..." << endl;
265
266 // Get RootManager
267 FairRootManager* ioman = FairRootManager::Instance();
268 if (! ioman) {
269 cout << "-E- " << GetName() << "::Init: "
270 << "RootManager not instantised!" << endl;
271 return kFATAL;
272 }
273
274 // Get MCTrack array
275 fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
276 if ( ! fMCTracks ) {
277 cout << "-E- " << GetName() << "::Init: No MCTrack array!" << endl;
278 return kFATAL;
279 }
280
281 // Get StsPoint array
282 fStsPoints = (TClonesArray*) ioman->GetObject("StsPoint");
283 if ( ! fStsPoints ) {
284 //cout << "-E- " << GetName() << "::Init: No StsPoint array!" << endl;
285 return kFATAL;
286 }
287
288 // Get StsHit array
289 fStsHits = (TClonesArray*) ioman->GetObject("StsHit");
290 if ( ! fStsHits ) {
291 cout << "-E- " << GetName() << "::Init: No StsHit array!" << endl;
292 return kFATAL;
293 }
294
295 // Get StsTrack array
296 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
297 if ( ! fStsTracks ) {
298 cout << "-E- " << GetName() << "::Init: No StsTrack array!" << endl;
299 return kERROR;
300 }
301
302 // Get StsTrackMatch array
303 fMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
304 if ( ! fMatches ) {
305 cout << "-E- " << GetName() << "::Init: No StsTrackMatch array!"
306 << endl;
307 return kERROR;
308 }
309
310 // Get StsDigis array
311 fStsDigis = (TClonesArray*) ioman->GetObject("StsDigi");
312 if ( ! fMatches ) {
313 cout << "-E- " << GetName() << "::Init: No StsDigi array!"
314 << endl;
315 return kERROR;
316 }
317
318 // Get PrimaryVertex array
319 // fPVertex = (CbmVertex*) ioman->GetObject("PrimaryVertex");
320 // CbmVertex* hohoho = (CbmVertex*) ioman->GetObject("PrimaryVertex");
321 // if ( ! hohoho ) {
322 // cout << "-E- " << GetName() << "::Init: No PrimaryVertex array!"
323 // << endl;
324 // return kERROR;
325 // }
326
327 // Get the geometry of target and STS
328 InitStatus geoStatus = GetGeometry();
329 if ( geoStatus != kSUCCESS ) {
330 cout << "-E- " << GetName() << "::Init: Error in reading geometry!"
331 << endl;
332 return geoStatus;
333 }
334
335 // Create histograms
336 CreateHistos();
337 Reset();
338
339 if ( fOnlineAnalysis ) {
340 fOnlineCanvas = new TCanvas("StsRecoOnline","Sts reconstruction online" ,10,10,600,900);
341 fOnlinePad[0] = new TPad("titlePad", "Title pad" ,0.01,0.91,0.99,0.99);
342 fOnlinePad[1] = new TPad("efficiencyPad","Efficiency pad" ,0.01,0.61,0.39,0.89);
343 fOnlinePad[2] = new TPad("resolutionPad","Momentum resolution pad" ,0.01,0.31,0.39,0.59);
344 fOnlinePad[3] = new TPad("hNpAccAll","Nof points reconstructuble tracks" ,0.41,0.66,0.69,0.89);
345 fOnlinePad[4] = new TPad("hNpRecAll","Nof points reconstructed track" ,0.71,0.66,0.99,0.89);
346 fOnlinePad[5] = new TPad("hStsTrackFPosZ","Param First pos Z" ,0.41,0.41,0.69,0.64);
347 fOnlinePad[6] = new TPad("hStsTrackLPosZ","Param Last pos Z" ,0.71,0.41,0.99,0.64);
348 fOnlinePad[7] = new TPad("hMomPrim","Momentum of primary tracks" ,0.41,0.16,0.69,0.41);
349 fOnlinePad[8] = new TPad("hMomSec","Momentum of secondary tracks" ,0.71,0.16,0.99,0.41);
350 fOnlinePad[9] = new TPad("printoutPad","Print information pad" ,0.01,0.01,0.39,0.29);
351 fOnlinePad[7]->SetLogy();
352 fOnlinePad[8]->SetLogy();
353 for ( Int_t ipad = 0 ; ipad < 10 ; ipad++ ) {
354 fOnlinePad[ipad]->SetFillColor(0);
355 fOnlinePad[ipad]->SetBorderMode(0);
356 fOnlinePad[ipad]->Draw();
357 }
358
359 fOnlinePad[0]->cd();
360 TLegend* brp = new TLegend(0.1,0.1,0.9,0.9,"Online STS reconstruction");
361 brp->SetTextAlign(22);
362 brp->SetTextSize(0.6);
363 brp->SetTextColor(1);
364 brp->SetBorderSize(0);
365 brp->SetFillColor(0);
366 brp->Draw();
367 fOnlinePad[0]->Update();
368 }
369
370 // Output
371 cout << " Minimum number of STS hits : " << fMinHits << endl;
372 cout << " Matching quota : " << fQuota << endl;
373 cout << " Target position ( " << fTargetPos.X() << ", "
374 << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl;
375 cout << " Number of STS stations : " << fNStations << endl;
376 if (fActive) cout << " ***** Task is ACTIVE *****" << endl;
377 cout << "==========================================================="
378 << endl << endl;
379
380 return geoStatus;
381
382}
383// -------------------------------------------------------------------------
384
385// ----- Public method ReInit ------------------------------------------
386InitStatus CbmStsReconstructionQa::ReInit() {
387
388 cout << "==========================================================="
389 << endl;;
390 cout << GetName() << ": Reinitialising..." << endl;
391
392 // Get the geometry of target and STS
393 InitStatus geoStatus = GetGeometry();
394 if ( geoStatus != kSUCCESS ) {
395 cout << "-E- " << GetName() << "::ReInit: Error in reading geometry!"
396 << endl;
397 return geoStatus;
398 }
399
400 // Output
401 cout << " Target position ( " << fTargetPos.X() << ", "
402 << fTargetPos.Y() << ", " << fTargetPos.Z() << ") " << endl;
403 cout << " Number of STS stations : " << fNStations << endl;
404 if (fActive) cout << " ***** Task is ACTIVE *****" << endl;
405 cout << "==========================================================="
406 << endl << endl;
407
408 return geoStatus;
409
410}
411// -------------------------------------------------------------------------
412
413
414// ----- Public method Exec --------------------------------------------
415void CbmStsReconstructionQa::Exec(Option_t* opt) {
416 // Timer
417 fTimer.Start();
418
419 // Eventwise counters
420 Int_t nRec = 0;
421 Int_t nGhosts = 0;
422 Int_t nClones = 0;
423 Int_t nAll = 0;
424 Int_t nAcc = 0;
425 Int_t nRecAll = 0;
426 Int_t nPrim = 0;
427 Int_t nRecPrim = 0;
428 Int_t nRef = 0;
429 Int_t nRecRef = 0;
430 Int_t nSec = 0;
431 Int_t nRecSec = 0;
432 TVector3 vertex;
433 TVector3 momentum;
434
435 // Fill hit and track maps
436
437
438 for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
439 for ( Int_t isect = 0 ; isect < fNSectors[istat] ; isect++ ) {
440 fNofHits[istat][isect] = 0;
441 for ( Int_t iside = 0 ; iside < 2 ; iside++ ) {
442 fNofFiredDigis[istat][isect][iside] = 0;
443 for ( Int_t ichip = 0 ; ichip < 8 ; ichip++ )
444 fNofDigisPChip[istat][isect][iside][ichip] = 0;
445 }
446 }
447 }
448
449 // Loop over MCTracks
450 Int_t nMC = fMCTracks->GetEntriesFast();
451
452 FillHitMap();
453 FillMatchMap(nRec, nGhosts, nClones);
454 for (Int_t iMC=0; iMC<nMC; iMC++) {
455 CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(iMC);
456 if ( ! mcTrack ) {
457 cout << "-E- " << GetName() << "::Exec: "
458 << "No MCTrack at index " << iMC
459 << endl;
460 Fatal("Exec", "No MCTrack in array");
461 }
462
463 mcTrack->GetStartVertex(vertex);
464 Bool_t isPrim = kFALSE;
465 mcTrack->GetMomentum(momentum);
466 Double_t mom = momentum.Mag();
467
468 if ( (vertex-fTargetPos).Mag() < 1. )
469 isPrim = kTRUE;
470
471 //Int_t pdgC = mcTrack->GetPdgCode();
472
473 /*Int_t toSave = -1;
474 if ( momentum.Z()>0.1 && momentum.Z()<2.6 )
475 toSave = (Int_t)((momentum.Z()-.1)*10.);
476 */
477 // if ( TMath::Abs(momentum.Z()-1.)<0.5 )
478
479 // if ( isPrim )
480 // if ( toSave >= 0 ) {
481 // if ( pdgC > 0 )
482 // fhDirEmiPrimP[toSave]->Fill(momentum.X(),momentum.Y());
483 // else
484 // fhDirEmiPrimM[toSave]->Fill(momentum.X(),momentum.Y());
485 // }
486
487 // Check reconstructability; continue only for reconstructable tracks
488 nAll++;
489 if ( fHitMap.find(iMC) == fHitMap.end() ) continue; // No hits
490 Int_t nHits = fHitMap[iMC];
491 Int_t nHitsSt = 0;
492 Int_t minHits = fMinHits;
493 for (Int_t i=0; i<9; i++) {
494 nHitsSt = 0;
495 nHitsSt = HitSt[iMC][i];
496 if (nHitsSt>1) {
497 minHits += nHitsSt;
498 }
499 }
500 if (minHits>4) minHits=minHits-1;
501 // if (nHits>3&&nHits<minHits) cout << minHits <<" "<<nHits<< endl;
502 if ( nHits < (minHits) ) continue; // Too few hits
503 nAcc++;
504
505 // Check origin of MCTrack
506 if ( isPrim ) nPrim++;
507 else nSec++;
508
509 // cout << ( isPrim ? "p" : "s" ) << nHits << flush;
510
511 // Get momentum
512 Bool_t isRef = kFALSE;
513 if ( mom > 1. && isPrim) {
514 isRef = kTRUE;
515 nRef++;
516 }
517
518 // Fill histograms for reconstructable tracks
519 fhMomAccAll->Fill(mom);
520 fhNpAccAll->Fill(Double_t(nHits));
521 if ( isPrim) {
522 fhMomAccPrim->Fill(mom);
523 fhNpAccPrim->Fill(Double_t(nHits));
524 // if ( toSave >= 0 ) {
525 // if ( pdgC > 0 )
526 // fhDirAccPrimP[toSave]->Fill(momentum.X(),momentum.Y());
527 // else
528 // fhDirAccPrimM[toSave]->Fill(momentum.X(),momentum.Y());
529 // }
530 }
531 else {
532 fhMomAccSec->Fill(mom);
533 fhNpAccSec->Fill(Double_t(nHits));
534 fhZAccSec->Fill(vertex.Z());
535 }
536
537 // Get matched StsTrack
538 Int_t iRec = -1;
539 Double_t quali = 0.;
540 Bool_t isRec = kFALSE;
541 if (fMatchMap.find(iMC) != fMatchMap.end() ) {
542 /* if ( isPrim ) {
543 cout << "mc mom = " << mom << " ("
544 << mcTrack->GetMomentum().X() << ", "
545 << mcTrack->GetMomentum().Y() << ", "
546 << mcTrack->GetMomentum().Z() << ") " << endl;
547 }*/
548 iRec = fMatchMap[iMC];
549 isRec = kTRUE;
550 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(iRec);
551 if ( ! stsTrack ) {
552 cout << "-E- " << GetName() << "::Exec: "
553 << "No StsTrack for matched MCTrack " << iMC << endl;
554 Fatal("Exec", "No StsTrack for matched MCTrack");
555 }
556 quali = fQualiMap[iMC];
557 if ( quali < fQuota ) {
558 cout << "-E- " << GetName() << "::Exec: "
559 << "Matched StsTrack " << iRec << " is below matching "
560 << "criterion ( " << quali << ")" << endl;
561 Fatal("Exec", "Match below matching quota");
562 }
563 CbmTrackMatch* match = (CbmTrackMatch*) fMatches->At(iRec);
564 if ( ! match ) {
565 cout << "-E- " << GetName() << "::Exec: "
566 << "No StsTrackMatch for matched MCTrack " << iMC << endl;
567 Fatal("Exec", "No StsTrackMatch for matched MCTrack");
568 }
569 Int_t nTrue = match->GetNofTrueHits();
570 Int_t nWrong = match->GetNofWrongHits();
571 Int_t nFake = match->GetNofFakeHits();
572 Int_t nAllHits = stsTrack->GetNStsHits();
573 if ( nTrue + nWrong + nFake != nAllHits ) {
574 cout << "True " << nTrue << " wrong " << nWrong << " Fake "
575 << nFake << " Hits " << nAllHits << endl;
576 Fatal("Exec", "Wrong number of hits");
577 }
578
579 // Verbose output
580 if ( fVerbose > 4 )
581 cout << "-I- " << GetName() << ": "
582 << "MCTrack " << iMC << ", hits "
583 << nAllHits << ", StsTrack " << iRec << ", hits " << nHits
584 << ", true hits " << nTrue << endl;
585
586 // Fill histograms for reconstructed tracks
587 if ( stsTrack->GetParamFirst()->GetQp() )
588 fhMomResAll->Fill(mom,100.*(mom-1./TMath::Abs(stsTrack->GetParamFirst()->GetQp()))/mom);
589 nRecAll++;
590 fhMomRecAll->Fill(mom);
591 fhNpRecAll->Fill(Double_t(nAllHits));
592 if ( isPrim ) {
593 nRecPrim++;
594 fhMomRecPrim->Fill(mom);
595 fhNpRecPrim->Fill(Double_t(nAllHits));
596
597 // if ( toSave >= 0 ) {
598 // if ( pdgC > 0 )
599 // fhDirRecPrimP[toSave]->Fill(momentum.X(),momentum.Y());
600 // else
601 // fhDirRecPrimM[toSave]->Fill(momentum.X(),momentum.Y());
602 // }
603
604 if ( isRef ) nRecRef++;
605 if ( stsTrack->GetParamFirst()->GetQp() )
606 fhMomResPrim->Fill(mom,100.*(mom-1./TMath::Abs(stsTrack->GetParamFirst()->GetQp()))/mom);
607 }
608 else {
609 nRecSec++;
610 fhMomRecSec->Fill(mom);
611 fhNpRecSec->Fill(Double_t(nHits));
612 fhZRecSec->Fill(vertex.Z());
613 if ( stsTrack->GetParamFirst()->GetQp() )
614 fhMomResSec->Fill(mom,100.*(mom-1./TMath::Abs(stsTrack->GetParamFirst()->GetQp()))/mom);
615 }
616
617 } // Match found in map?
618
619 Int_t partPdgCode = mcTrack->GetPdgCode();
620 for ( Int_t itemp = 0 ; itemp < 10 ; itemp++ ) {
621 if ( fPartPdgTable[itemp] == -7777 ) break;
622 if ( fPartPdgTable[itemp] != partPdgCode ) continue;
623 fhMomAccPart[itemp]->Fill(mom);
624 if ( isRec )
625 fhMomRecPart[itemp]->Fill(mom);
626 }
627
628 } // Loop over MCTracks
629
630 Int_t nofStsHits = fStsHits->GetEntriesFast();
631 //Int_t nofStsPoints = fStsPoints->GetEntriesFast();
632 // cout << "there are " << nofStsPoints << " points and " << nofStsHits << " hits." << endl;
633 Int_t hitStationLimits[2][100];
634
635 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
636 hitStationLimits[0][ist] = -1;
637 hitStationLimits[1][ist] = -1;
638 }
639
640 // nof digis
641 Int_t nofStsDigis = fStsDigis->GetEntriesFast();
642 for ( Int_t idigi = 0 ; idigi < nofStsDigis ; idigi++ ) {
643 CbmStsDigi *stsDigi = (CbmStsDigi*)fStsDigis->At(idigi);
644 fNofFiredDigis[stsDigi->GetStationNr()-1][stsDigi->GetSectorNr()][stsDigi->GetSide()] += 1;
645
646 fNofDigisPChip[stsDigi->GetStationNr()-1][stsDigi->GetSectorNr()][stsDigi->GetSide()][(Int_t)(stsDigi->GetChannelNr()/125)] += 1;
647 }
648
649
650 // check for limits of hit indices on different stations...
651 for ( Int_t ihit = 0 ; ihit < nofStsHits ; ihit++ ) {
652 CbmStsHit *stsHit = (CbmStsHit*)fStsHits->At(ihit);
653 fNofHits[stsHit->GetStationNr()-1][stsHit->GetSectorNr()] += 1; // count hits per sector
654 if ( hitStationLimits[0][stsHit->GetStationNr()-1] == -1 )
655 hitStationLimits[0][stsHit->GetStationNr()-1] = ihit;
656 CbmStsHit *stsHitBack = (CbmStsHit*)fStsHits->At(nofStsHits-ihit-1);
657 if ( hitStationLimits[1][stsHitBack->GetStationNr()-1] == -1 ) {
658 hitStationLimits[1][stsHitBack->GetStationNr()-1] = nofStsHits-ihit;
659 }
660 }
661
662 // for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
663 // for ( Int_t isect = 0 ; isect < fNSectors[istat] ; isect++ ) {
664 // fhNofHits[istat][isect] -> Fill(fNofHits[istat][isect]);
665 // for ( Int_t iside = 0 ; iside < 2 ; iside++ ) {
666 // fhNofFiredDigis[istat][isect][iside]->Fill(fNofFiredDigis[istat][isect][iside]);
667 // for ( Int_t ichip = 0 ; ichip < 8 ; ichip++ ) {
668 // if ( fhNofDigisPChip[istat][isect][iside][ichip] )
669 // fhNofDigisPChip[istat][isect][iside][ichip]->Fill(fNofDigisPChip[istat][isect][iside][ichip]);
670 // }
671 // }
672
673 // }
674 // }
675 /*
676 for ( Int_t ipnt = 0 ; ipnt < nofStsPoints ; ipnt++ ) {
677 CbmStsPoint *stsPoint = (CbmStsPoint*)fStsPoints->At(ipnt);
678
679 Int_t startHit = hitStationLimits[0][fStationNrFromMcId[stsPoint->GetDetectorID()]];
680 Int_t finalHit = hitStationLimits[1][fStationNrFromMcId[stsPoint->GetDetectorID()]];
681
682 // fhEnergyLoss[fStationNrFromMcId[stsPoint->GetDetectorID()]]->Fill(stsPoint->GetXIn(),stsPoint->GetYIn(),stsPoint->GetEnergyLoss());
683
684// Float_t zP = stsPoint->GetZ();
685// Float_t xP = stsPoint->GetX(zP);
686// Float_t yP = stsPoint->GetY(zP);
687
688 if ( startHit == -1 && finalHit == -1 ) continue;
689
690 for ( Int_t ihit = startHit ; ihit < finalHit ; ihit++ ) {
691 CbmStsHit *stsHit= (CbmStsHit*)fStsHits->At(ihit);
692 if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ())) < .1 ) &&
693 ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsHit->GetZ())) < .1 ) )
694 // fhHitPointCorrelation[fStationNrFromMcId[stsPoint->GetDetectorID()]]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
695 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
696// if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsPoint->GetZ())) < 1. ) &&
697// ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsPoint->GetZ())) < 1. ) )
698// fhHitPointCorrelation[fStationNrFromMcId[stsPoint->GetDetectorID()]]->Fill(stsHit->GetX()-stsPoint->GetX(stsPoint->GetZ()),
699// stsHit->GetY()-stsPoint->GetY(stsPoint->GetZ()));
700 }
701 }*/
702
703 // Calculate efficiencies
704 Double_t effAll = 1.;
705 if ( nAcc ) effAll = Double_t(nRecAll) / Double_t(nAcc);
706 Double_t effPrim = 1.;
707 if ( nPrim ) effPrim = Double_t(nRecPrim) / Double_t(nPrim);
708 Double_t effRef = 1.;
709 if ( nRef ) effRef = Double_t(nRecRef) / Double_t(nRef);
710 Double_t effSec = 1.;
711 if ( nSec ) effSec = Double_t(nRecSec) / Double_t(nSec);
712
713 fhRefTracks ->SetBinContent(fNEvents+1,nRef);
714 fhRecRefTracks->SetBinContent(fNEvents+1,nRecRef);
715
716 // Event summary
717 if ( fVerbose > 1 ) {
718 cout << "---------- StsReconstructionQa : Event " << fNEvents+1 << " summary ------------"
719 << endl;
720 cout << "MCTracks : " << nAll << ", reconstructable: " << nAcc
721 << ", reconstructed: " << nRecAll << endl;
722 cout << "Vertex : reconstructable: " << nPrim << ", reconstructed: "
723 << nRecPrim << ", efficiency " << effPrim*100. << "%" << endl;
724 cout << "Reference : reconstructable: " << nRef << ", reconstructed: "
725 << nRecRef << ", efficiency " << effRef*100. << "%" << endl;
726 cout << "Non-vertex : reconstructable: " << nSec << ", reconstructed: "
727 << nRecSec << ", efficiency " << effSec*100. << "%" << endl;
728 cout << "STSTracks " << nRec << ", ghosts " << nGhosts
729 << ", clones " << nClones << endl;
730 cout << "-----------------------------------------------------------"
731 << endl;
732 cout << endl;
733 }
734 if ( fVerbose == 1 ) {
735 cout << "\r+ " << setw(15) << left << fName << ": event " << fNEvents+1 << " " << setprecision(4)
736 << setw(8) << fixed << right << fTimer.RealTime()
737 << " s, efficiency all " << effAll*100. << " %, vertex "
738 << effPrim*100. << " %, reference " << effRef*100. << " %" << endl;
739 }
740
741 // Increase counters
742 fNAccAll += nAcc;
743 fNAccPrim += nPrim;
744 fNAccRef += nRef;
745 fNAccSec += nSec;
746 fNRecAll += nRecAll;
747 fNRecPrim += nRecPrim;
748 fNRecRef += nRecRef;
749 fNRecSec += nRecSec;
750 fNGhosts += nGhosts;
751 fNClones += nClones;
752 fNStsTracks += nRec;
753 fNEvents++;
754 fTime += fTimer.RealTime();
755
756 effRef = 1.;
757 if ( fNAccRef )
758 effRef = Double_t(fNRecRef) / Double_t(fNAccRef);
759
760 if ( fOnlineAnalysis ) {
761 fOnlinePad[1]->cd();
762 DivideHistos(fhMomRecAll, fhMomAccAll, fhMomEffAll );
763 DivideHistos(fhMomRecPrim, fhMomAccPrim, fhMomEffPrim);
764 DivideHistos(fhMomRecSec, fhMomAccSec, fhMomEffSec );
765 fhMomEffAll->SetAxisRange(0.,1.1,"Y");
766 fhMomEffAll->SetLineWidth(2);
767 fhMomEffAll->SetLineColor(1);
768 fhMomEffAll->SetTitle("Efficiency");
769 fhMomEffAll->Draw();
770 fhMomEffPrim->SetLineWidth(2);
771 fhMomEffPrim->SetLineColor(2);
772 fhMomEffPrim->Draw("same");
773 fhMomEffSec->SetLineWidth(2);
774 fhMomEffSec->SetLineColor(3);
775 fhMomEffSec->Draw("same");
776 TLegend* effLeg = new TLegend(0.3,0.15,0.48,0.4);
777 effLeg->SetBorderSize(0);
778 effLeg->SetFillColor(0);
779 effLeg->AddEntry(fhMomEffAll, "all" ,"pl");
780 effLeg->AddEntry(fhMomEffPrim,"prim","pl");
781 effLeg->AddEntry(fhMomEffSec, "sec" ,"pl");
782 effLeg->Draw();
783 TLine* oneLine = new TLine(0.0,1.0,10.0,1.0);
784 oneLine->SetLineStyle(2);
785 oneLine->Draw();
786 fOnlinePad[1]->Update();
787
788 fOnlinePad[2]->cd();
789 if ( fhMomResPrim->Integral() )
790 fOnlinePad[2]->SetLogz();
791 fhMomResPrim->SetAxisRange(0.,3.,"Y");
792 fhMomResPrim->Draw("cont0");
793 for ( Int_t ibin = fhMomResPrim->GetXaxis()->GetNbins() ; ibin > 1 ; ibin-- ) {
794 TF1* gausFit = new TF1("gausFit","gaus");
795 TH1F* tempProjY = (TH1F*)fhMomResPrim->ProjectionY("tempProjY",ibin,ibin);
796 tempProjY->Fit("gausFit","QN","",-5.,5.);
797 fhLowBand->SetBinContent(ibin,gausFit->GetParameter(1)-gausFit->GetParameter(2));
798 fhLowBand->SetBinError(ibin,0.01);
799 fhHigBand->SetBinContent(ibin,gausFit->GetParameter(2));
800 fhHigBand->SetBinError(ibin,gausFit->GetParError(2));
801 // fhHigBand->SetBinContent(ibin,gausFit->GetParameter(1)+gausFit->GetParameter(2));
802 // fhHigBand->SetBinError(ibin,0.01);
803 }
804 fhLowBand->SetMarkerSize(0.2);
805 fhLowBand->SetLineWidth(2);
806 fhHigBand->SetMarkerSize(0.1);
807 fhHigBand->SetLineWidth(2);
808 fhLowBand->Draw("Psame");
809 fhHigBand->Draw("Psame");
810 fOnlinePad[2]->Update();
811
812 fOnlinePad[3]->cd();
813 fhNpAccAll->Draw();
814 fOnlinePad[3]->Update();
815
816 fOnlinePad[4]->cd();
817 fhNpRecAll->Draw();
818 fOnlinePad[4]->Update();
819
820 fOnlinePad[5]->cd();
821 fhStsTrackFPos[2]->Draw();
822 fOnlinePad[5]->Update();
823
824 fOnlinePad[6]->cd();
825 fhStsTrackLPos[2]->Draw();
826 fOnlinePad[6]->Update();
827
828 fOnlinePad[7]->cd();
829 fhMomAccPrim->SetLineWidth(2);
830 fhMomAccPrim->SetLineColor(3);
831 fhMomAccPrim->Draw();
832 fhMomRecPrim->SetLineColor(2);
833
834 fhMomRecPrim->Draw("same");
835 TLegend* momLeg = new TLegend(0.55,0.45,0.72,0.8);
836 momLeg->SetBorderSize(0);
837 momLeg->SetFillColor(0);
838 momLeg->SetTextSize(0.07);
839 momLeg->AddEntry(fhMomAccPrim, "acc prim" ,"pl");
840 momLeg->AddEntry(fhMomRecPrim, "rec prim" ,"pl");
841 momLeg->Draw();
842 fOnlinePad[7]->Update();
843
844 fOnlinePad[8]->cd();
845 fhMomAccSec->SetLineWidth(2);
846 fhMomAccSec->SetLineColor(3);
847 fhMomAccSec->Draw();
848 fhMomRecSec->SetLineColor(2);
849
850 fhMomRecSec->Draw("same");
851 TLegend* momsLeg = new TLegend(0.55,0.45,0.72,0.8);
852 momsLeg->SetBorderSize(0);
853 momsLeg->SetFillColor(0);
854 momsLeg->SetTextSize(0.07);
855 momsLeg->AddEntry(fhMomAccSec, "acc sec" ,"pl");
856 momsLeg->AddEntry(fhMomRecSec, "rec sec" ,"pl");
857 momsLeg->Draw();
858 fOnlinePad[8]->Update();
859
860 TF1* allEffFit = new TF1 ("allEffFit","pol0",1.,10.);
861 fhMomEffAll->Fit(allEffFit,"QN","",1,10);
862 Double_t allEff = allEffFit->GetParameter(0);
863 effAll = 1.;
864 if ( fhMomAccAll->Integral() )
865 effAll = fhMomRecAll->Integral()/fhMomAccAll->Integral();
866 TF1* primEffFit = new TF1 ("primEffFit","pol0",1.,10.);
867 fhMomEffPrim->Fit(primEffFit,"QN","",1,10);
868 Double_t primEff = primEffFit->GetParameter(0);
869 effPrim = 1.;
870 if ( fhMomAccPrim->Integral() )
871 effPrim = fhMomRecPrim->Integral()/fhMomAccPrim->Integral();
872 TF1* secEffFit = new TF1 ("secEffFit","pol0",1.,10.);
873 fhMomEffSec->Fit(secEffFit,"QN","",1,10);
874 Double_t secEff = secEffFit->GetParameter(0);
875 effSec = 1.;
876 if ( fhMomAccSec->Integral() )
877 effSec = fhMomRecSec->Integral()/fhMomAccSec->Integral();
878
879 TF1* momentumResFuncPrim = new TF1("momentumResFuncPrim","gaus",-10.,10.);
880 TH1F* momentumResHistPrim = (TH1F*)fhMomResPrim->ProjectionY("momentumResHistPrim");
881 momentumResHistPrim->Fit(momentumResFuncPrim,"QN","",-10.,10.);
882 Double_t momentumResolutionPrim = momentumResFuncPrim->GetParameter(2);
883 TF1* momentumResFuncAll = new TF1("momentumResFuncAll","gaus",-10.,10.);
884 TH1F* momentumResHistAll = (TH1F*)fhMomResAll->ProjectionY("momentumResHistAll");
885 momentumResHistAll->Fit(momentumResFuncAll,"QN","",-10.,10.);
886 Double_t momentumResolutionAll = momentumResFuncAll->GetParameter(2);
887
888 fOnlinePad[9]->cd();
889 TPaveText* printoutPave = new TPaveText(0.0,0.0,1.0,1.0);
890 printoutPave->SetTextAlign(23);
891 printoutPave->SetTextSize(0.05);
892 printoutPave->SetTextColor(1);
893 printoutPave->SetBorderSize(0);
894 printoutPave->SetFillColor(0);
895 printoutPave->AddText(Form("%i events",fNEvents));
896 printoutPave->AddText(Form("%3.2f prim, %3.2f sec, %3.2f gh, %3.2f cl",
897 Double_t (fNRecPrim)/Double_t (fNEvents),
898 Double_t (fNRecSec) /Double_t (fNEvents),
899 Double_t (fNGhosts) /Double_t (fNEvents),
900 Double_t (fNClones) /Double_t (fNEvents)));
901 // printoutPave->AddText("Single Hit Resolutions:");
902 // for ( Int_t ist = 0 ; ist < fNStations ; ist++ )
903 // if ( resolution[0][ist] > 0.01 )
904 // printoutPave->AddText(Form("st#%i,#sigma_{x}=%3.2f#mum,#sigma_{y}=%3.2f#mum",
905 // ist+1,resolution[0][ist],resolution[1][ist]));
906 printoutPave->AddText("Tracking efficiencies (p>1.0 GeV/c):");
907 printoutPave->AddText(Form("all = %2.2f%%(%2.2f%%)",100.*effAll,100.*allEff));
908 printoutPave->AddText(Form("vertex = %2.2f%%(%2.2f%%)",100.*effPrim,100.*primEff));
909 printoutPave->AddText(Form("reference = %2.2f%%",100.*effRef));
910 printoutPave->AddText(Form("non-vertex = %2.2f%%(%2.2f%%)",100.*effSec,100.*secEff));
911 printoutPave->AddText(Form("Momentum resolution = %3.2f%%(%3.2f%%)",momentumResolutionAll,momentumResolutionPrim));
912 fOnlinePad[9]->Clear();
913 printoutPave->Draw();
914 fOnlinePad[9]->Update();
915 }
916}
917// -------------------------------------------------------------------------
918
919
920// ----- Private method Finish -----------------------------------------
921void CbmStsReconstructionQa::Finish() {
922
923 // Divide histograms for efficiency calculation
924 DivideHistos(fhMomRecAll, fhMomAccAll, fhMomEffAll);
925 DivideHistos(fhMomRecPrim, fhMomAccPrim, fhMomEffPrim);
926 DivideHistos(fhMomRecSec, fhMomAccSec, fhMomEffSec);
927 DivideHistos(fhNpRecAll, fhNpAccAll, fhNpEffAll);
928 DivideHistos(fhNpRecPrim, fhNpAccPrim, fhNpEffPrim);
929 DivideHistos(fhNpRecSec, fhNpAccSec, fhNpEffSec);
930 DivideHistos(fhZRecSec, fhZAccSec, fhZEffSec);
931
932 // for ( Int_t itemp = 0 ; itemp < 25 ; itemp++ ) {
933 // fhDirAcMPrimM[itemp]->Divide(fhDirAccPrimM[itemp],fhDirEmiPrimM[itemp]);
934 // fhDirEffPrimM[itemp]->Divide(fhDirRecPrimM[itemp],fhDirAccPrimM[itemp]);
935 // fhDirAcMPrimP[itemp]->Divide(fhDirAccPrimP[itemp],fhDirEmiPrimP[itemp]);
936 // fhDirEffPrimP[itemp]->Divide(fhDirRecPrimP[itemp],fhDirAccPrimP[itemp]);
937 // }
938
939 for ( Int_t itemp = 0 ; itemp < 10 ; itemp++ ) {
940 if ( fPartPdgTable[itemp] == -7777 ) break;
941 DivideHistos(fhMomRecPart[itemp], fhMomAccPart[itemp], fhMomEffPart[itemp]);
942 }
943
944 // Normalise histos for clones and ghosts to one event
945 if ( fNEvents ) {
946 fhNhClones->Scale(1./Double_t(fNEvents));
947 fhNhGhosts->Scale(1./Double_t(fNEvents));
948 }
949
950 // Calculate integrated efficiencies and rates
951 Double_t effAll = 1.;
952 if ( fNAccAll ) effAll = Double_t(fNRecAll) / Double_t(fNAccAll);
953 Double_t effPrim = 1.;
954 if ( fNAccPrim ) effPrim = Double_t(fNRecPrim) / Double_t(fNAccPrim);
955 Double_t effRef = 1.;
956 if ( fNAccRef ) effRef = Double_t(fNRecRef) / Double_t(fNAccRef);
957 Double_t effSec = 1.;
958 if ( fNAccSec ) effSec = Double_t(fNRecSec) / Double_t(fNAccSec);
959 Double_t rateGhosts = Double_t(fNGhosts) / Double_t(fNEvents);
960 Double_t rateClones = Double_t(fNClones) / Double_t(fNEvents);
961
962 // Run summary to screen
963 cout << endl;
964 cout << "============================================================"
965 << endl;
966 cout << "===== " << fName << ": Run summary " << endl;
967 cout << "===== " << endl;
968 cout << "===== Good events : " << setw(6) << fNEvents << endl;
969 cout << "===== Failed events : " << setw(6) << fNEventsFailed << endl;
970 cout << "===== Average time : " << setprecision(4) << setw(8) << right
971 << fTime / Double_t(fNEvents) << " s" << endl;
972 cout << "===== " << endl;
973 cout << "===== Efficiency all tracks : " << effAll*100 << " % ("
974 << fNRecAll << "/" << fNAccAll <<")" << endl;
975 cout << "===== Efficiency vertex tracks : " << effPrim*100 << " % ("
976 << fNRecPrim << "/" << fNAccPrim <<")" << endl;
977 cout << "===== Efficiency reference tracks : " << effRef*100 << " % ("
978 << fNRecRef << "/" << fNAccRef <<")" << endl;
979 cout << "===== Efficiency secondary tracks : " << effSec*100 << " % ("
980 << fNRecSec << "/" << fNAccSec <<")" << endl;
981 cout << "===== Ghost rate " << rateGhosts << " per event" << endl;
982 cout << "===== Clone rate " << rateClones << " per event" << endl;
983 cout << "============================================================"
984 << endl;
985
986 // Write histos to output
987 gDirectory->mkdir("STSReconstructionQA");
988 gDirectory->cd("STSReconstructionQA");
989 TIter next(fHistoList);
990 while ( TH1* histo = ((TH1*)next()) ) histo->Write();
991
992 // gDirectory->mkdir("STSOccupancy");
993 // gDirectory->cd("STSOccupancy");
994 // TIter nextO(fOccupHList);
995 // while ( TH1* histo = ((TH1*)nextO()) ) histo->Write();
996
997 // gDirectory->cd("..");
998 gDirectory->cd("..");
999}
1000// -------------------------------------------------------------------------
1001
1002
1003// ----- Private method GetGeometry ------------------------------------
1004InitStatus CbmStsReconstructionQa::GetGeometry() {
1005
1006 // Get target geometry
1007 if ( ! fPassGeo ) {
1008 cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!"
1009 <<endl;
1010 fTargetPos.SetXYZ(0., 0., 0.);
1011 return kERROR;
1012 }
1013 TObjArray* passNodes = fPassGeo->GetGeoPassiveNodes();
1014 if ( ! passNodes ) {
1015 cout << "-W- " << GetName() << "::GetGeometry: No passive node array"
1016 << endl;
1017 fTargetPos.SetXYZ(0., 0., 0.);
1018 return kERROR;
1019 }
1020 FairGeoNode* target = (FairGeoNode*) passNodes->FindObject("targ");
1021 if ( ! target ) {
1022 cout << "-E- " << GetName() << "::GetGeometry: No target node"
1023 << endl;
1024 fTargetPos.SetXYZ(0., 0., 0.);
1025 return kERROR;
1026 }
1027 FairGeoVector targetPos = target->getLabTransform()->getTranslation();
1028 FairGeoVector centerPos = target->getCenterPosition().getTranslation();
1029 Double_t targetX = targetPos.X() + centerPos.X();
1030 Double_t targetY = targetPos.Y() + centerPos.Y();
1031 Double_t targetZ = targetPos.Z() + centerPos.Z();
1032 fTargetPos.SetXYZ(targetX, targetY, targetZ);
1033
1034 // Get STS geometry
1035 if ( ! fStsGeo ) {
1036 cout << "-W- " << GetName() << "::GetGeometry: No passive geometry!"
1037 <<endl;
1038 fNStations = 0;
1039 return kERROR;
1040 }
1041 TObjArray* stsNodes = fStsGeo->GetGeoSensitiveNodes();
1042 if ( ! stsNodes ) {
1043 cout << "-E- " << GetName() << "::GetGeometry: No STS node array"
1044 << endl;
1045 fNStations = 0;
1046 return kERROR;
1047 }
1048 Int_t tempNofStations = stsNodes->GetEntries();
1049
1050 if ( fVerbose > 2 )
1051 cout << "There are " << tempNofStations << " nodes"
1052 << (tempNofStations > 10 ? "!!!" : "" ) << endl;
1053
1054 TString geoNodeName;
1055 fNStations = 0;
1056 TString stationNames[100];
1057 Bool_t stationKnown = kFALSE;
1058
1059 for ( Int_t istat = 0 ; istat < 20 ; istat++ )
1060 fNSectors[istat] = 0;
1061
1062 for ( Int_t ist = 0 ; ist < tempNofStations ; ist++ ) {
1063
1064 FairGeoNode* stsNode = (FairGeoNode*)stsNodes->At(ist);
1065 if ( ! stsNode ) {
1066 cout << "-W- CbmStsDigiScheme::Init: station#" << ist
1067 << " not found among sensitive nodes " << endl;
1068 continue;
1069 }
1070 geoNodeName = stsNode->GetName();
1071 // TArrayD* params = stsNode->getParameters();
1072
1073 stationKnown = kFALSE;
1074
1075 FairGeoVector* pt1 = stsNode->getPoint(0);
1076 FairGeoVector* pt2 = stsNode->getPoint(1);
1077 Double_t sectorWidth = TMath::Abs(pt1->getX()-pt2->getX());
1078 // cout << geoNodeName.Data() << " -> " << sectorWidth << " cm wide" << endl;
1079
1080 // check if the node belongs to some station, save the MCId
1081 for ( Int_t ikst = 0 ; ikst < fNStations ; ikst++ )
1082 if ( geoNodeName.Contains(stationNames[ikst]) ) {
1083 fStationNrFromMcId[stsNode->getMCid()] = ikst;
1084 stationKnown = kTRUE;
1085 fWidthSectors[ikst][fNSectors[ikst]] = sectorWidth;
1086 fNSectors[ikst]++;
1087 }
1088
1089 if ( stationKnown ) continue;
1090
1091 // if not known, register it and save MCId
1092 fStationNrFromMcId[stsNode->getMCid()] = fNStations;
1093 fWidthSectors[fNStations][fNSectors[fNStations]] = sectorWidth;
1094 fNSectors[fNStations]++;
1095
1096 // it will work only if the node name is organized as:
1097 // for station name is "stsstationXX", where XX is the station number (f.e. XX=07 for station number 7)
1098 // for sector name is "stsstationXXanythingHereToDistinguishDifferentSectors"
1099 geoNodeName.Remove(12,geoNodeName.Length()-12);
1100 stationNames[fNStations] = geoNodeName.Data();
1101 fNStations++;
1102
1103 if ( fVerbose > 3 )
1104 cout << "station #" << fNStations << " has MCID = "
1105 << stsNode->getMCid() << " and name " << stsNode->GetName() << endl;
1106
1107 // fStationsMCId[fNStations] = stsNode->getMCid(); // not used
1108 }
1109 cout << " Stations: " << fNStations << " contain" << flush;
1110 for ( Int_t istat = 0 ; istat < fNStations ; istat++ )
1111 cout << " " << fNSectors[istat] << flush;
1112 cout << " sensors" << endl;
1113
1114 return kSUCCESS;
1115
1116}
1117// -------------------------------------------------------------------------
1118
1119// ----- Private method CreateHistos -----------------------------------
1120void CbmStsReconstructionQa::CreateHistos() {
1121
1122 // Histogram list
1123 fHistoList = new TList();
1124 // fOccupHList = new TList();
1125
1126 // Momentum distributions
1127 Double_t minMom = 0.;
1128 Double_t maxMom = 10.;
1129 Int_t nBinsMom = 40;
1130 fhMomAccAll = new TH1F("hMomAccAll", "all reconstructable tracks",
1131 nBinsMom, minMom, maxMom);
1132 fhMomRecAll = new TH1F("hMomRecAll", "all reconstructed tracks",
1133 nBinsMom, minMom, maxMom);
1134 fhMomEffAll = new TH1F("hMomEffAll", "efficiency all tracks",
1135 nBinsMom, minMom, maxMom);
1136 fhMomEffAll->SetXTitle("p [GeV/c]");
1137 fhMomEffAll->SetYTitle("efficiency");
1138 fhMomAccPrim = new TH1F("hMomAccPrim", "reconstructable vertex tracks",
1139 nBinsMom, minMom, maxMom);
1140 fhMomRecPrim = new TH1F("hMomRecPrim", "reconstructed vertex tracks",
1141 nBinsMom, minMom, maxMom);
1142 fhMomEffPrim = new TH1F("hMomEffPrim", "efficiency vertex tracks",
1143 nBinsMom, minMom, maxMom);
1144 fhMomAccSec = new TH1F("hMomAccSec", "reconstructable non-vertex tracks",
1145 nBinsMom, minMom, maxMom);
1146 fhMomRecSec = new TH1F("hMomRecSec", "reconstructed non-vertex tracks",
1147 nBinsMom, minMom, maxMom);
1148 fhMomEffSec = new TH1F("hMomEffSec", "efficiency non-vertex tracks",
1149 nBinsMom, minMom, maxMom);
1150 fhMomGhosts = new TH1F("hMomGhosts", "momenta of ghosts",
1151 nBinsMom, minMom, maxMom);
1152 fhMomClones = new TH1F("hMomClones", "momenta of clones",
1153 nBinsMom, minMom, maxMom);
1154 fHistoList->Add(fhMomAccAll);
1155 fHistoList->Add(fhMomRecAll);
1156 fHistoList->Add(fhMomEffAll);
1157 fHistoList->Add(fhMomAccPrim);
1158 fHistoList->Add(fhMomRecPrim);
1159 fHistoList->Add(fhMomEffPrim);
1160 fHistoList->Add(fhMomAccSec);
1161 fHistoList->Add(fhMomRecSec);
1162 fHistoList->Add(fhMomEffSec);
1163 fHistoList->Add(fhMomClones);
1164 fHistoList->Add(fhMomGhosts);
1165
1166 // for ( Int_t itemp = 0 ; itemp < 25 ; itemp++ ) {
1167 // fhDirEmiPrimM[itemp] = new TH2F(Form("hDirEmiPrimM%d",itemp), "emitted vertex tracks",
1168 // 100,-1.,1.,100,-1.,1.);
1169 // fhDirAccPrimM[itemp] = new TH2F(Form("hDirAccPrimM%d",itemp), "reconstructable vertex tracks",
1170 // 100,-1.,1.,100,-1.,1.);
1171 // fhDirAcMPrimM[itemp] = new TH2F(Form("hDirAcMPrimM%d",itemp), "acceptance vertex tracks",
1172 // 100,-1.,1.,100,-1.,1.);
1173 // fhDirRecPrimM[itemp] = new TH2F(Form("hDirRecPrimM%d",itemp), "reconstructed vertex tracks",
1174 // 100,-1.,1.,100,-1.,1.);
1175 // fhDirEffPrimM[itemp] = new TH2F(Form("hDirEffPrimM%d",itemp), "efficiency vertex tracks",
1176 // 100,-1.,1.,100,-1.,1.);
1177 // fHistoList->Add(fhDirEmiPrimM[itemp]);
1178 // fHistoList->Add(fhDirAccPrimM[itemp]);
1179 // fHistoList->Add(fhDirAcMPrimM[itemp]);
1180 // fHistoList->Add(fhDirRecPrimM[itemp]);
1181 // fHistoList->Add(fhDirEffPrimM[itemp]);
1182 // fhDirEmiPrimP[itemp] = new TH2F(Form("hDirEmiPrimP%d",itemp), "emitted vertex tracks",
1183 // 100,-1.,1.,100,-1.,1.);
1184 // fhDirAccPrimP[itemp] = new TH2F(Form("hDirAccPrimP%d",itemp), "reconstructable vertex tracks",
1185 // 100,-1.,1.,100,-1.,1.);
1186 // fhDirAcMPrimP[itemp] = new TH2F(Form("hDirAcMPrimP%d",itemp), "acceptance vertex tracks",
1187 // 100,-1.,1.,100,-1.,1.);
1188 // fhDirRecPrimP[itemp] = new TH2F(Form("hDirRecPrimP%d",itemp), "reconstructed vertex tracks",
1189 // 100,-1.,1.,100,-1.,1.);
1190 // fhDirEffPrimP[itemp] = new TH2F(Form("hDirEffPrimP%d",itemp), "efficiency vertex tracks",
1191 // 100,-1.,1.,100,-1.,1.);
1192 // fHistoList->Add(fhDirEmiPrimP[itemp]);
1193 // fHistoList->Add(fhDirAccPrimP[itemp]);
1194 // fHistoList->Add(fhDirAcMPrimP[itemp]);
1195 // fHistoList->Add(fhDirRecPrimP[itemp]);
1196 // fHistoList->Add(fhDirEffPrimP[itemp]);
1197 // }
1198
1199
1200 for ( Int_t itemp = 0 ; itemp < 10 ; itemp++ ) {
1201 if ( fPartPdgTable[itemp] == -7777 ) break;
1202 if ( fVerbose > 3 )
1203 cout << "fpart pdg table content for itemp = " << itemp
1204 << " equals " << fPartPdgTable[itemp] << endl;
1205 fhMomAccPart[itemp] = new TH1F(Form("hMomAccPart%s%d",(fPartPdgTable[itemp]>0?"P":"M"),TMath::Abs(fPartPdgTable[itemp])),
1206 Form("reconstruable particle%d tracks",fPartPdgTable[itemp]),
1207 nBinsMom, minMom, maxMom);
1208 fhMomRecPart[itemp] = new TH1F(Form("hMomRecPart%s%d",(fPartPdgTable[itemp]>0?"P":"M"),TMath::Abs(fPartPdgTable[itemp])),
1209 Form("reconstructed particle%d tracks",fPartPdgTable[itemp]),
1210 nBinsMom, minMom, maxMom);
1211 fhMomEffPart[itemp] = new TH1F(Form("hMomEffPart%s%d",(fPartPdgTable[itemp]>0?"P":"M"),TMath::Abs(fPartPdgTable[itemp])),
1212 Form("efficiency particle%d tracks",fPartPdgTable[itemp]),
1213 nBinsMom, minMom, maxMom);
1214 fHistoList->Add(fhMomAccPart[itemp]);
1215 fHistoList->Add(fhMomRecPart[itemp]);
1216 fHistoList->Add(fhMomEffPart[itemp]);
1217 }
1218
1219
1220 // Number-of-points distributions
1221 Double_t minNp = -0.5;
1222 Double_t maxNp = 15.5;
1223 Int_t nBinsNp = 16;
1224 fhNpAccAll = new TH1F("hNpAccAll", "all reconstructable tracks",
1225 nBinsNp, minNp, maxNp);
1226 fhNpRecAll = new TH1F("hNpRecAll", "all reconstructed tracks",
1227 nBinsNp, minNp, maxNp);
1228 fhNpEffAll = new TH1F("hNpEffAll", "efficiency all tracks",
1229 nBinsNp, minNp, maxNp);
1230 fhNpAccPrim = new TH1F("hNpAccPrim", "reconstructable vertex tracks",
1231 nBinsNp, minNp, maxNp);
1232 fhNpRecPrim = new TH1F("hNpRecPrim", "reconstructed vertex tracks",
1233 nBinsNp, minNp, maxNp);
1234 fhNpEffPrim = new TH1F("hNpEffPrim", "efficiency vertex tracks",
1235 nBinsNp, minNp, maxNp);
1236 fhNpAccSec = new TH1F("hNpAccSec", "reconstructable non-vertex tracks",
1237 nBinsNp, minNp, maxNp);
1238 fhNpRecSec = new TH1F("hNpRecSec", "reconstructed non-vertex tracks",
1239 nBinsNp, minNp, maxNp);
1240 fhNpEffSec = new TH1F("hNpEffSec", "efficiency non-vertex tracks",
1241 nBinsNp, minNp, maxNp);
1242 fHistoList->Add(fhNpAccAll);
1243 fHistoList->Add(fhNpRecAll);
1244 fHistoList->Add(fhNpEffAll);
1245 fHistoList->Add(fhNpAccPrim);
1246 fHistoList->Add(fhNpRecPrim);
1247 fHistoList->Add(fhNpEffPrim);
1248 fHistoList->Add(fhNpAccSec);
1249 fHistoList->Add(fhNpRecSec);
1250 fHistoList->Add(fhNpEffSec);
1251
1252 // z(vertex) distributions
1253 Double_t minZ = 0.;
1254 Double_t maxZ = 50.;
1255 Int_t nBinsZ = 50;
1256 fhZAccSec = new TH1F("hZAccSec", "reconstructable non-vertex tracks",
1257 nBinsZ, minZ, maxZ);
1258 fhZRecSec = new TH1F("hZRecSecl", "reconstructed non-vertex tracks",
1259 nBinsZ, minZ, maxZ);
1260 fhZEffSec = new TH1F("hZEffRec", "efficiency non-vertex tracks",
1261 nBinsZ, minZ, maxZ);
1262 fHistoList->Add(fhZAccSec);
1263 fHistoList->Add(fhZRecSec);
1264 fHistoList->Add(fhZEffSec);
1265
1266 // Number-of-hit distributions
1267 fhNhClones = new TH1F("hNhClones", "number of hits for clones",
1268 nBinsNp, minNp, maxNp);
1269 fhNhGhosts = new TH1F("hNhGhosts", "number of hits for ghosts",
1270 nBinsNp, minNp, maxNp);
1271 fhNhGhosts->SetXTitle("# of hits");
1272 fhNhGhosts->SetYTitle("yield [a.u.]");
1273 fHistoList->Add(fhNhClones);
1274 fHistoList->Add(fhNhGhosts);
1275
1276 fhMomResAll = new TH2F("hMomResAll", "momentum resolution vs p for all tracks",
1277 nBinsMom,minMom,maxMom,
1278 20,-10.,10.);
1279 fhMomResPrim = new TH2F("hMomResPrim","momentum resolution vs p for vertex tracks",
1280 nBinsMom,minMom,maxMom,
1281 20,-10.,10.);
1282 fhMomResPrim->SetXTitle("p [GeV/c]");
1283 fhMomResPrim->SetYTitle("#delta p/p [%%]");
1284 fhMomResSec = new TH2F("hMomResSec", "momentum resolution vs p for non-vertex tracks",
1285 nBinsMom,minMom,maxMom,
1286 20,-10.,10.);
1287 fHistoList->Add(fhMomResAll);
1288 fHistoList->Add(fhMomResPrim);
1289 fHistoList->Add(fhMomResSec);
1290 fhLowBand = new TH1F("hLowBand","Lower Band",nBinsMom,minMom,maxMom);
1291 fhHigBand = new TH1F("hHigBand","Higher band",nBinsMom,minMom,maxMom);
1292 fHistoList->Add(fhLowBand);
1293 fHistoList->Add(fhHigBand);
1294
1295 /* for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
1296 fhHitPointCorrelation[ist] = new TH2F(Form("hHitPointCorrelation%i",ist+1),
1297 Form("Hit vs point correlation at station %i",ist+1),
1298 500,-.1, .1,500,-.1,.1);
1299 fhHitPointCorrelation[ist]->SetXTitle("#Delta x [cm]");
1300 fhHitPointCorrelation[ist]->SetYTitle("#Delta y [cm]");
1301 fHistoList->Add(fhHitPointCorrelation[ist]);
1302 }*/
1303
1304 fhPrimaryVertex = new TH3F("hPrimaryVertex","Primary vertex",200,-0.1,0.1,200,-0.1,0.1,200,-0.1,0.1);
1305 fHistoList->Add(fhPrimaryVertex);
1306
1307 fhRefTracks = new TH1F ("hRefTracks" ,"Nof reconstructed reference tracks",100,-0.5,999.5);
1308 fhRecRefTracks = new TH1F ("hRecRefTracks","Nof reconstruable reference tracks",100,-0.5,999.5);
1309 fHistoList->Add(fhRefTracks);
1310 fHistoList->Add(fhRecRefTracks);
1311
1312 Char_t lett[5] = {'X','Y','Z','x','y'};
1313
1314 for ( Int_t itemp = 0 ; itemp < 15 ; itemp++ ) {
1315 // if ( itemp < 6 || itemp == 7) {
1316 // fhStsTrackFCovEl[itemp] = new TH1F(Form("hStsTrackFCovEl%d",itemp),
1317 // Form("StsTrack ParamFirst cov. el. %d",itemp),
1318 // 200,-1.e-3,1.e-3);
1319 // fhStsTrackLCovEl[itemp] = new TH1F(Form("hStsTrackLCovEl%d",itemp),
1320 // Form("StsTrack ParamLast cov. el. %d",itemp),
1321 // 200,-1.e-3,1.e-3);
1322 // }
1323 // else {
1324 // if ( itemp == 6 ) {
1325 // fhStsTrackFCovEl[itemp] = new TH1F(Form("hStsTrackFCovEl%d",itemp),
1326 // Form("StsTrack ParamFirst cov. el. %d",itemp),
1327 // 200,-1.e-3,1.e-3);
1328 // fhStsTrackLCovEl[itemp] = new TH1F(Form("hStsTrackLCovEl%d",itemp),
1329 // Form("StsTrack ParamLast cov. el. %d",itemp),
1330 // 200,-1.e-3,1.e-3);
1331 // }
1332 // else {
1333 // if ( itemp == 8 ) {
1334 // fhStsTrackFCovEl[itemp] = new TH1F(Form("hStsTrackFCovEl%d",itemp),
1335 // Form("StsTrack ParamFirst cov. el. %d",itemp),
1336 // 200,-1.e-2,1.e-2);
1337 // fhStsTrackLCovEl[itemp] = new TH1F(Form("hStsTrackLCovEl%d",itemp),
1338 // Form("StsTrack ParamLast cov. el. %d",itemp),
1339 // 200,-1.e-2,1.e-2);
1340 // }
1341 // else {
1342 // fhStsTrackFCovEl[itemp] = new TH1F(Form("hStsTrackFCovEl%d",itemp),
1343 // Form("StsTrack ParamFirst cov. el. %d",itemp),
1344 // 200,-1.e-2,1.e-2);
1345 // fhStsTrackLCovEl[itemp] = new TH1F(Form("hStsTrackLCovEl%d",itemp),
1346 // Form("StsTrack ParamLast cov. el. %d",itemp),
1347 // 200,-1.e-2,1.e-2);
1348 // }
1349 // }
1350 // }
1351
1352 // fHistoList->Add(fhStsTrackFCovEl[itemp]);
1353 // fHistoList->Add(fhStsTrackLCovEl[itemp]);
1354
1355 if ( itemp >= 3 ) continue;
1356 Int_t nof = 100 ;
1357 Double_t beg = -50.;
1358 Double_t end = 50.;
1359 if ( itemp == 2 ) { nof = 120; beg = -10.; end = 110.; }
1360 fhStsTrackFPos[itemp] = new TH1F(Form("hStsTrackFPos%c",lett[itemp]),
1361 Form("StsTrack ParamFirst pos %c",lett[itemp]),
1362 nof,beg,end);
1363 fhStsTrackLPos[itemp] = new TH1F(Form("hStsTrackLPos%c",lett[itemp]),
1364 Form("StsTrack ParamLast pos %c",lett[itemp]),
1365 nof,beg,end);
1366 fHistoList->Add(fhStsTrackFPos[itemp]);
1367 fHistoList->Add(fhStsTrackLPos[itemp]);
1368
1369 if ( itemp >= 2 ) continue;
1370 fhStsTrackFDir[itemp] = new TH1F(Form("hStsTrackFDir%c",lett[itemp+3]),
1371 Form("StsTrack ParamFirst dir %c",lett[itemp+3]),
1372 10,-1.5,1.5);
1373 fhStsTrackLDir[itemp] = new TH1F(Form("hStsTrackLDir%c",lett[itemp+3]),
1374 Form("StsTrack ParamLast dir %c",lett[itemp+3]),
1375 10,-1.5,1.5);
1376 fHistoList->Add(fhStsTrackFDir[itemp]);
1377 fHistoList->Add(fhStsTrackLDir[itemp]);
1378 }
1379 fhStsTrackFMom = new TH1F("hStsTrackFMom", "Momentum of rec. tracks ParFirst", 100,-50.,50.);
1380 fhStsTrackLMom = new TH1F("hStsTrackLMom", "Momentum of rec. tracks ParLast" , 100,-50.,50.);
1381 fhStsTrackChiSq = new TH1F("hStsTrackChiSq","Chi square of rec. tracks",100,0.,1000.);
1382 fHistoList->Add(fhStsTrackFMom);
1383 fHistoList->Add(fhStsTrackLMom);
1384 fHistoList->Add(fhStsTrackChiSq);
1385
1386 /*Double_t binningNofHits[501];
1387 binningNofHits[0] = -0.5;
1388 for ( Int_t itemp = 0 ; itemp < 100 ; itemp++ ) {
1389 binningNofHits[itemp+ 1] = 0.5+(Double_t)itemp;
1390 binningNofHits[itemp+101] = 101.5+2.*(Double_t)itemp;
1391 binningNofHits[itemp+201] = 302.5+3.*(Double_t)itemp;
1392 binningNofHits[itemp+301] = 602.5+3.*(Double_t)itemp;
1393 binningNofHits[itemp+401] = 902.5+3.*(Double_t)itemp;
1394 }*/
1395 // for ( Int_t itemp = 0 ; itemp < 501 ; itemp++ )
1396 // cout << binningNofHits[itemp] << " " << flush;
1397 // cout << endl;
1398
1399 /* for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
1400 fhEnergyLoss[istat] = new TH2F(Form("hEnergyLossSt%d",istat+1),
1401 Form("Energy loss on station %d",istat+1),
1402 200,-100.,100.,200,-100.,100.);
1403 fOccupHList->Add(fhEnergyLoss[istat]);
1404 for ( Int_t isect = 0 ; isect < fNSectors[istat] ; isect++ ) {
1405 fhNofHits[istat][isect] = new TH1F(Form("hNofHitsSt%dSect%d",istat+1,isect+1),
1406 Form("Number of hits in sector %d of station %d",isect+1,istat+1),
1407 500,binningNofHits);
1408 // 500,-0.5,500.5);
1409 fOccupHList->Add(fhNofHits[istat][isect]);
1410
1411 Int_t nofChips = (Int_t)(TMath::Ceil(fWidthSectors[istat][isect]/7.5)); // fwidth in mm, 7.5mm = 125(channels)*60mum(pitch)
1412 Int_t lastChip = (Int_t)(TMath::Ceil(10.*fWidthSectors[istat][isect]));
1413 lastChip = lastChip%75;
1414 lastChip = (Int_t)(lastChip/.6);
1415
1416 TString addInfo = "";
1417 if ( nofChips != 8 ) {
1418 addInfo = Form(", only %d strips",lastChip);
1419 // cout << fWidthSectors[istat][isect] << " -> " << addInfo.Data() << endl;
1420 }
1421
1422 for ( Int_t iside = 0 ; iside < 2 ; iside++ ) {
1423 fhNofFiredDigis[istat][isect][iside] = new TH1F(Form("hNofFiredDigis%cSt%dSect%d",(iside==0?'F':'B'),istat+1,isect+1),
1424 Form("Number of digis on %s of sector %d of station %d",(iside==0?"front":"back"),isect+1,istat+1),
1425 501,-0.5,500.5);
1426 fOccupHList->Add(fhNofFiredDigis[istat][isect][iside]);
1427
1428 for ( Int_t ichip = 0 ; ichip < nofChips ; ichip++ ) {
1429 fhNofDigisPChip[istat][isect][iside][ichip] = new TH1F(Form("hNofFiredDigis%cSt%dSect%dChip%d",(iside==0?'F':'B'),istat+1,isect+1,ichip+1),
1430 Form("Number of digis on %s on chip %d of sector %d of station %d%s",(iside==0?"front":"back"),ichip+1,isect+1,istat+1,(ichip==nofChips-1?addInfo.Data():"")),
1431 101,-0.5,100.5);
1432 fOccupHList->Add(fhNofDigisPChip[istat][isect][iside][ichip]);
1433 }
1434 }
1435 }
1436 }*/
1437}
1438// -------------------------------------------------------------------------
1439
1440
1441
1442// ----- Private method Reset ------------------------------------------
1443void CbmStsReconstructionQa::Reset() {
1444
1445 TIter next(fHistoList);
1446 while ( TH1* histo = ((TH1*)next()) ) histo->Reset();
1447 TIter next0(fOccupHList);
1448 while ( TH1* histo = ((TH1*)next0()) ) histo->Reset();
1449
1450 fNAccAll = fNAccPrim = fNAccRef = fNAccSec = 0;
1451 fNRecAll = fNRecPrim = fNRecRef = fNRecSec = 0;
1452 fNGhosts = fNClones = fNEvents = 0;
1453 fNStsTracks = 0;
1454}
1455// -------------------------------------------------------------------------
1456
1457
1458// ----- Private method FillHitMap -------------------------------------
1459void CbmStsReconstructionQa::FillHitMap() {
1460 fHitMap.clear();
1461 Int_t nMC = fMCTracks->GetEntriesFast();
1462 for (Int_t i=0; i<nMC; i++) {
1463 for (Int_t j=0; j<10; j++){
1464 HitSt[i][j]=0;
1465 }
1466 }
1467 Int_t nHits = fStsHits->GetEntriesFast();
1468 for (Int_t iHit=0; iHit<nHits; iHit++) {
1469 CbmStsHit* hit = (CbmStsHit*) fStsHits->At(iHit);
1470 Int_t iMc = hit->GetRefIndex();
1471 if ( iMc < 0 ) continue;
1472 CbmStsPoint* stsPoint = (CbmStsPoint*) fStsPoints->At(iMc);
1473 Int_t iTrack = stsPoint->GetTrackID();
1474 HitSt [iTrack][hit->GetStationNr()]++;
1475 fHitMap[iTrack]++;
1476 }
1477}
1478// -------------------------------------------------------------------------
1479
1480
1481// ------ Private method FillMatchMap ----------------------------------
1482void CbmStsReconstructionQa::FillMatchMap(Int_t& nRec, Int_t& nGhosts,
1483 Int_t& nClones) {
1484
1485 // Clear matching maps
1486 fMatchMap.clear();
1487 fQualiMap.clear();
1488
1489 // Loop over StsTracks. Check matched MCtrack and fill maps.
1490 nGhosts = 0;
1491 nClones = 0;
1492 nRec = fStsTracks->GetEntriesFast();
1493 Int_t nMtc = fMatches->GetEntriesFast();
1494 if ( nMtc != nRec ) {
1495 cout << "-E- " << GetName() << "::Exec: Number of StsMatches ("
1496 << nMtc << ") does not equal number of StsTracks ("
1497 << nRec << ")" << endl;
1498 Fatal("Exec", "Inequal number of StsTrack and StsTrackMatch");
1499 }
1500 for (Int_t iRec=0; iRec<nRec; iRec++) {
1501
1502 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(iRec);
1503 if ( ! stsTrack ) {
1504 cout << "-E- " << GetName() << "::Exec: "
1505 << "No StsTrack at index " << iRec << endl;
1506 Fatal("Exec", "No StsTrack in array");
1507 }
1508 Int_t nHits = stsTrack->GetNStsHits();
1509
1510 FairTrackParam* trParF = (FairTrackParam*)stsTrack->GetParamFirst();
1511 FairTrackParam* trParL = (FairTrackParam*)stsTrack->GetParamLast();
1512 fhStsTrackFPos[0]->Fill(trParF->GetX());
1513 fhStsTrackFPos[1]->Fill(trParF->GetY());
1514 fhStsTrackFPos[2]->Fill(trParF->GetZ());
1515 fhStsTrackLPos[0]->Fill(trParL->GetX());
1516 fhStsTrackLPos[1]->Fill(trParL->GetY());
1517 fhStsTrackLPos[2]->Fill(trParL->GetZ());
1518 fhStsTrackFDir[0]->Fill(trParF->GetTx());fhStsTrackFDir[1]->Fill(trParF->GetTy());
1519 fhStsTrackLDir[0]->Fill(trParL->GetTx());fhStsTrackLDir[1]->Fill(trParL->GetTy());
1520 Int_t elToFill = 0;
1521 for ( Int_t ixx = 0 ; ixx < 5 ; ixx++ ) {
1522 for ( Int_t iyy = ixx ; iyy < 5 ; iyy++ ) {
1523 // fhStsTrackFCovEl[elToFill]->Fill(trParF->GetCovariance(ixx,iyy));
1524 // fhStsTrackLCovEl[elToFill]->Fill(trParL->GetCovariance(ixx,iyy));
1525 elToFill++;
1526 }
1527 }
1528 fhStsTrackFMom->Fill(trParF->GetQp());
1529 fhStsTrackLMom->Fill(trParL->GetQp());
1530 fhStsTrackChiSq->Fill(stsTrack->GetChi2());
1531
1532 CbmTrackMatch* match = (CbmTrackMatch*) fMatches->At(iRec);
1533 if ( ! match ) {
1534 cout << "-E- " << GetName() << "::Exec: "
1535 << "No StsTrackMatch at index " << iRec << endl;
1536 Fatal("Exec", "No StsTrackMatch in array");
1537 }
1538 Int_t nTrue = match->GetNofTrueHits();
1539
1540 Int_t iMC = match->GetMCTrackId();
1541 if (iMC == -1 ) { // no common point with MC, really ghastly!
1542 if ( fVerbose > 4 )
1543 cout << "-I- " << GetName() << ":"
1544 << "No MC match for StsTrack " << iRec << endl;
1545 fhNhGhosts->Fill(nHits);
1546 if ( stsTrack->GetParamFirst()->GetQp() )
1547 fhMomGhosts->Fill(1./TMath::Abs(stsTrack->GetParamFirst()->GetQp()));
1548 nGhosts++;
1549 continue;
1550 }
1551
1552 // Check matching criterion (quota)
1553 Double_t quali = 1.;
1554 if ( nHits )
1555 quali = Double_t(nTrue) / Double_t(nHits);
1556 if ( quali >= fQuota ) {
1557
1558 // No previous match for this MCTrack
1559 if ( fMatchMap.find(iMC) == fMatchMap.end() ) {
1560 fMatchMap[iMC] = iRec;
1561 fQualiMap[iMC] = quali;
1562 }
1563
1564 // Previous match; take the better one
1565 else {
1566 if ( fVerbose > 4 )
1567 cout << "-I- " << GetName() << ": "
1568 << "MCTrack " << iMC << " doubly matched."
1569 << "Current match " << iRec
1570 << ", previous match " << fMatchMap[iMC]
1571 << endl;
1572 if ( fQualiMap[iMC] < quali ) {
1573 CbmStsTrack* oldTrack
1574 = (CbmStsTrack*) fStsTracks->At(fMatchMap[iMC]);
1575 fhNhClones->Fill(Double_t(oldTrack->GetNStsHits()));
1576 if ( oldTrack->GetParamFirst()->GetQp() )
1577 fhMomClones->Fill(1./TMath::Abs(oldTrack->GetParamFirst()->GetQp()));
1578 fMatchMap[iMC] = iRec;
1579 fQualiMap[iMC] = quali;
1580 }
1581 else {
1582 fhNhClones->Fill(nHits);
1583 if ( stsTrack->GetParamFirst()->GetQp() )
1584 fhMomClones->Fill(1./TMath::Abs(stsTrack->GetParamFirst()->GetQp()));
1585 }
1586 nClones++;
1587 }
1588
1589 }
1590
1591 // If not matched, it's a ghost
1592 else {
1593 if ( fVerbose > 4 )
1594 cout << "-I- " << GetName() << ":"
1595 << "StsTrack " << iRec << " below matching criterion "
1596 << "(" << quali << ")" << endl;
1597 fhNhGhosts->Fill(nHits);
1598 if ( stsTrack->GetParamFirst()->GetQp() )
1599 fhMomGhosts->Fill(1./TMath::Abs(stsTrack->GetParamFirst()->GetQp()));
1600 nGhosts++;
1601 }
1602
1603 } // Loop over StsTracks
1604
1605}
1606// -------------------------------------------------------------------------
1607
1608// ----- Private method DivideHistos -----------------------------------
1609void CbmStsReconstructionQa::DivideHistos(TH1* histo1, TH1* histo2,
1610 TH1* histo3) {
1611
1612 if ( !histo1 || !histo2 || !histo3 ) {
1613 cout << "-E- " << GetName() << "::DivideHistos: "
1614 << "NULL histogram pointer" << endl;
1615 Fatal("DivideHistos", "Null histo pointer");
1616 }
1617
1618 Int_t nBins = histo1->GetNbinsX();
1619 if ( histo2->GetNbinsX() != nBins || histo3->GetNbinsX() != nBins ) {
1620 cout << "-E- " << GetName() << "::DivideHistos: "
1621 << "Different bin numbers in histos" << endl;
1622 cout << histo1->GetName() << " " << histo1->GetNbinsX() << endl;
1623 cout << histo2->GetName() << " " << histo2->GetNbinsX() << endl;
1624 cout << histo3->GetName() << " " << histo3->GetNbinsX() << endl;
1625 return;
1626 }
1627
1628 Double_t c1, c2, c3, ce;
1629 for (Int_t iBin=0; iBin<nBins; iBin++) {
1630 c1 = histo1->GetBinContent(iBin);
1631 c2 = histo2->GetBinContent(iBin);
1632 if ( c2 ) {
1633 c3 = c1 / c2;
1634 if ( c3 <= 1. )
1635 ce = TMath::Sqrt( c3 * ( 1. - c3 ) / c2 );
1636 else
1637 ce = TMath::Sqrt( c3 * ( 1. + c3 ) / c2 );
1638 }
1639 else {
1640 c3 = 0.;
1641 ce = 0.;
1642 }
1643 histo3->SetBinContent(iBin, c3);
1644 histo3->SetBinError(iBin, ce);
1645 }
1646
1647}
1648// -------------------------------------------------------------------------
i
int i
Definition P4_F32vec4.h:22
CbmGeoStsPar::GetGeoSensitiveNodes
TObjArray * GetGeoSensitiveNodes()
Definition CbmGeoStsPar.h:20
CbmMCTrack::GetMomentum
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:139
CbmMCTrack::GetStartVertex
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:149
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition CbmMCTrack.h:56
CbmStsDigi::GetSide
Int_t GetSide() const
Definition CbmStsDigi.h:79
CbmStsDigi::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsDigi.h:75
CbmStsDigi::GetChannelNr
Int_t GetChannelNr() const
Definition CbmStsDigi.h:83
CbmStsDigi::GetStationNr
Int_t GetStationNr() const
Definition CbmStsDigi.h:72
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:66
CbmStsHit::GetSectorNr
Int_t GetSectorNr() const
Definition CbmStsHit.h:68
CbmStsReconstructionQa::Exec
virtual void Exec(Option_t *opt)
Definition CbmStsReconstructionQa.cxx:415
CbmStsReconstructionQa::CbmStsReconstructionQa
CbmStsReconstructionQa(Int_t iVerbose=1)
Definition CbmStsReconstructionQa.cxx:44
CbmStsTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition CbmStsTrack.h:70
CbmStsTrack::GetNStsHits
Int_t GetNStsHits() const
Definition CbmStsTrack.h:60
CbmStsTrack::GetChi2
Double_t GetChi2() const
Definition CbmStsTrack.h:66
CbmStsTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition CbmStsTrack.h:69
FairGeoPassivePar::GetGeoPassiveNodes
TObjArray * GetGeoPassiveNodes()
Definition FairGeoPassivePar.h:32
run
-clang-format
std
STL namespace.