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