BmnRoot
Loading...
Searching...
No Matches
CbmStsFindHitsQa.cxx
Go to the documentation of this file.
1// -------------------------------------------------------------------------
2// ----- CbmStsFindHitsQa source file -----
3// ----- Created 01/07/2008 by R. Karabowicz -----
4// -------------------------------------------------------------------------
5
6#include "CbmStsFindHitsQa.h"
7
8#include "CbmGeoStsPar.h"
9#include "CbmStsDigi.h"
10#include "CbmStsCluster.h"
11#include "CbmStsDigiPar.h"
12#include "CbmStsDigiScheme.h"
13#include "CbmStsPoint.h"
14#include "CbmStsHit.h"
15#include "CbmStsSensor.h"
16#include "CbmStsSector.h"
17#include "CbmStsStation.h"
18#include "CbmMCTrack.h"
19
20#include "FairRootManager.h"
21#include "FairRunAna.h"
22#include "FairRuntimeDb.h"
23
24#include "TClonesArray.h"
25#include "TMath.h"
26#include "TGeoManager.h"
27#include "TCanvas.h"
28#include "TLine.h"
29#include "TPad.h"
30#include "TLegend.h"
31#include "TPaveText.h"
32#include "TFile.h"
33#include "TF1.h"
34#include "TH1F.h"
35#include "TH2F.h"
36#include "TH3F.h"
37#include "TList.h"
38#include "TVector3.h"
39
40#include <iomanip>
41using namespace std;
42
43// ----- Default constructor ------------------------------------------
45: FairTask("STS Hit Finder", 1),
46 fGeoPar(NULL),
47 fDigiPar(NULL),
48 fDigiScheme(NULL),
49 fMCTracks(NULL), // MCtrack
50 fStsPoints(NULL), // StsPoints
51 fStsHits(NULL), // StsHits
52 fMatches(NULL), // StsTrackMatch
53 fStsDigis(NULL), // StsDigi
54 fStsClusters(NULL), // StsCluster
55 fTimer(),
56 fhHitFindingEfficiency(),
57 fhEffIncAng(),
58 fhEffMom(),
59 fhEffPdgSec(),
60 fhEffPdgPrim(),
61 fhHitPointPull(),
62 fhHitPointCorr(),
63 fHistoList(),
64 fNStations(0),
65 fNEvents(0),
66 fTime1(0.),
67 fNofHits(),
68 fNofPointsPrim(),
69 fNofPointsSec(),
70 fNofRecoPrim(),
71 fNofRecoSec(),
72 fNofPointsMomSum(),
73 fNofRecoPointsMomSum(),
74 fOnlineAnalysis(kFALSE),
75 recoCanvas()
76{
77 //AZ fDigiScheme = new CbmStsDigiScheme();
78 fDigiScheme = CbmStsDigiScheme::Instance();
79}
80
81// ----- Standard constructor ------------------------------------------
82CbmStsFindHitsQa::CbmStsFindHitsQa(Bool_t visualizeBool, Int_t iVerbose)
83: FairTask("STSFindHitsQa", iVerbose),
84 fGeoPar(NULL),
85 fDigiPar(NULL),
86 fDigiScheme(NULL),
87 fMCTracks(NULL), // MCtrack
88 fStsPoints(NULL), // StsPoints
89 fStsHits(NULL), // StsHits
90 fMatches(NULL), // StsTrackMatch
91 fStsDigis(NULL), // StsDigi
92 fStsClusters(NULL), // StsCluster
93 fTimer(),
94 fhHitFindingEfficiency(),
95 fhEffIncAng(),
96 fhEffMom(),
97 fhEffPdgSec(),
98 fhEffPdgPrim(),
99 fhHitPointPull(),
100 fhHitPointCorr(),
101 fHistoList(),
102 fNStations(0),
103 fNEvents(0),
104 fTime1(0.),
105 fNofHits(),
106 fNofPointsPrim(),
107 fNofPointsSec(),
108 fNofRecoPrim(),
109 fNofRecoSec(),
110 fNofPointsMomSum(),
111 fNofRecoPointsMomSum(),
112 fOnlineAnalysis(visualizeBool),
113 recoCanvas()
114{
115 //AZ fDigiScheme = new CbmStsDigiScheme();
116 fDigiScheme = CbmStsDigiScheme::Instance();
117}
118
120: FairTask("STSFindHitsQa", iVerbose),
121 fGeoPar(NULL),
122 fDigiPar(NULL),
123 fDigiScheme(NULL),
124 fMCTracks(NULL), // MCtrack
125 fStsPoints(NULL), // StsPoints
126 fStsHits(NULL), // StsHits
127 fMatches(NULL), // StsTrackMatch
128 fStsDigis(NULL), // StsDigi
129 fStsClusters(NULL), // StsCluster
130 fTimer(),
131 fhHitFindingEfficiency(),
132 fhEffIncAng(),
133 fhEffMom(),
134 fhEffPdgSec(),
135 fhEffPdgPrim(),
136 fhHitPointPull(),
137 fhHitPointCorr(),
138 fHistoList(),
139 fNStations(0),
140 fNEvents(0),
141 fTime1(0.),
142 fNofHits(),
143 fNofPointsPrim(),
144 fNofPointsSec(),
145 fNofRecoPrim(),
146 fNofRecoSec(),
147 fNofPointsMomSum(),
148 fNofRecoPointsMomSum(),
149 recoCanvas()
150{
151 //AZ fDigiScheme = new CbmStsDigiScheme();
152 fDigiScheme = CbmStsDigiScheme::Instance();
153}
154
155// ----- Constructor with name -----------------------------------------
156CbmStsFindHitsQa::CbmStsFindHitsQa(const char* name, Int_t iVerbose)
157: FairTask(name, iVerbose),
158 fGeoPar(NULL),
159 fDigiPar(NULL),
160 fDigiScheme(NULL),
161 fMCTracks(NULL), // MCtrack
162 fStsPoints(NULL), // StsPoints
163 fStsHits(NULL), // StsHits
164 fMatches(NULL), // StsTrackMatch
165 fStsDigis(NULL), // StsDigi
166 fStsClusters(NULL), // StsCluster
167 fTimer(),
168 fhHitFindingEfficiency(),
169 fhEffIncAng(),
170 fhEffMom(),
171 fhEffPdgSec(),
172 fhEffPdgPrim(),
173 fhHitPointPull(),
174 fhHitPointCorr(),
175 fHistoList(),
176 fNStations(0),
177 fNEvents(0),
178 fTime1(0.),
179 fNofHits(),
180 fNofPointsPrim(),
181 fNofPointsSec(),
182 fNofRecoPrim(),
183 fNofRecoSec(),
184 fNofPointsMomSum(),
185 fNofRecoPointsMomSum(),
186 fOnlineAnalysis(kFALSE),
187 recoCanvas()
188{
189 //AZ fDigiScheme = new CbmStsDigiScheme();
190 fDigiScheme = CbmStsDigiScheme::Instance();
191}
192
193// ----- Destructor ----------------------------------------------------
195 fHistoList->Delete();
196 delete fHistoList;
197 for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
198 fHistoListPS[istat]->Delete();
199 }
200 delete[] fHistoList;
201}
202// -------------------------------------------------------------------------
203
204// ----- Public method Exec --------------------------------------------
205void CbmStsFindHitsQa::Exec(Option_t* opt) {
206
207 fTimer.Start();
208 //Bool_t warn = kFALSE;
209
210 // Check for digi scheme
211 if ( ! fDigiScheme ) {
212 cerr << "-E- " << fName << "::Exec: No digi scheme!" << endl;
213 return;
214 }
215
216 CbmStsStation* station = NULL;
217 CbmStsSector* sector = NULL;
218
219 Int_t nofStsDigis = fStsDigis->GetEntriesFast();
220 for ( Int_t idigi = 0 ; idigi < nofStsDigis ; idigi++ ) {
221 CbmStsDigi *stsDigi = (CbmStsDigi*)fStsDigis->At(idigi);
222 // cout << "digi side = " << stsDigi->GetSide() << endl;
223 // fNofFiredDigis[stsDigi->GetStationNr()-1][stsDigi->GetSectorNr()-1][stsDigi->GetSide()] += 1;
224 fNofDigisPChip[stsDigi->GetStationNr()-1][stsDigi->GetSectorNr()-1][stsDigi->GetSide()][(Int_t)(stsDigi->GetChannelNr()/128)] += 1;
225 }
226 Int_t nofStsHits = fStsHits->GetEntriesFast();
227 Int_t nofStsPoints = fStsPoints->GetEntriesFast();
228 Int_t hitStationLimits[2][100];
229
230 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
231 hitStationLimits[0][ist] = -1;
232 hitStationLimits[1][ist] = -1;
233 }
234
235 // check for limits of hit indices on different stations...
236
237 for ( Int_t ihit = 0 ; ihit < nofStsHits ; ihit++ ) {
238 CbmStsHit *stsHit = (CbmStsHit*)fStsHits->At(ihit);
239 fNofHits[stsHit->GetStationNr()-1][stsHit->GetSectorNr()-1][stsHit->GetSensorNr()-1] += 1; // count hits per sensor
240 if ( hitStationLimits[0][stsHit->GetStationNr()-1] == -1 ) {
241 hitStationLimits[0][stsHit->GetStationNr()-1] = ihit;
242 }
243 CbmStsHit *stsHitBack = (CbmStsHit*)fStsHits->At(nofStsHits-ihit-1);
244 if ( hitStationLimits[1][stsHitBack->GetStationNr()-1] == -1 ) {
245 hitStationLimits[1][stsHitBack->GetStationNr()-1] = nofStsHits-ihit;
246 }
247 }
248 for ( Int_t istat = 0 ; istat < fDigiScheme->GetNStations() ; istat++ ) {
249 station = fDigiScheme->GetStation(istat);
250 for ( Int_t isect = 0 ; isect < station->GetNSectors() ; isect++ ) {
251 sector = station->GetSector(isect);
252
253 CbmStsSensor* sensor = (CbmStsSensor*)sector->GetSensor(0);
254 for ( Int_t iside = 0 ; iside<2; iside++){
255 for (Int_t iChip = 0; iChip<8; iChip++){
256 if (fNEvents==999) cout<<" St "<<istat+1<<" sect "<<isect+1<<" nofsens "<<sector->GetNSensors()<<" sectX "<<sensor->GetX0()<<" sectY "<<sensor->GetY0()<<" side "<<iside<<" chip "<<iChip+1<<" fNofDigisPChip "<<fNofDigisPChip[istat][isect][iside][iChip]<<endl;
257 }
258 }
259
260 for ( Int_t isens = 0 ; isens < sector->GetNSensors() ; isens++ ) {
261 // cout << "filling " << istat << " " << isect << " " << isens << endl;
262 fhNofHits[istat][isect][isens] -> Fill(fNofHits[istat][isect][isens]);
263 // cout << "---> done " << isens << endl;
264
265 }
266 // for ( Int_t iside = 0 ; iside < 2 ; iside++ ) {
267 // fhNofFiredDigis[istat][isect][iside]->Fill(fNofFiredDigis[istat][isect][iside]);
268 // for ( Int_t ichip = 0 ; ichip < 8 ; ichip++ ) {
269 // if ( fhNofDigisPChip[istat][isect][iside][ichip] )
270 // fhNofDigisPChip[istat][isect][iside][ichip]->Fill(fNofDigisPChip[istat][isect][iside][ichip]);
271 // }
272 // }
273
274 }
275 }
276 for ( Int_t ipnt = 0 ; ipnt < nofStsPoints ; ipnt++ ) {
277 CbmStsPoint *stsPoint = (CbmStsPoint*)fStsPoints->At(ipnt);
278
279 // cout << "point " << ipnt << flush;
280
281 Int_t trackPdg = 0;
282 Int_t motherId = 0;
283 //Int_t motherPdg = 0;
284 TVector3 trackVertex;
285 TVector3 motherVertex;
286 TVector3 momentum;
287
288 Double_t xin = stsPoint->GetXIn();
289 Double_t yin = stsPoint->GetYIn();
290 Double_t zin = stsPoint->GetZIn();
291 //Double_t xout = stsPoint->GetXOut();
292 //Double_t yout = stsPoint->GetYOut();
293 //Double_t zout = stsPoint->GetZOut();
294
295 gGeoManager->FindNode(xin,yin,zin);
296 //TGeoNode* curNode = gGeoManager->GetCurrentNode();
297
298 Int_t trackId = stsPoint->GetTrackID();
299 CbmMCTrack *mcTrack = (CbmMCTrack*) fMCTracks->At(trackId);
300 trackPdg = mcTrack -> GetPdgCode();
301 motherId = mcTrack -> GetMotherId();
302 mcTrack->GetMomentum(momentum);
303 Double_t mom = momentum.Mag();
304
305 if (motherId>=0 ) {
306 CbmMCTrack *motherTrack = (CbmMCTrack*) fMCTracks->At(motherId);
307 //motherPdg = motherTrack->GetPdgCode();
308 mcTrack->GetStartVertex(trackVertex);
309 motherTrack->GetStartVertex(motherVertex);
310
311 }
312
313 // cout << " in node \"" << curNode->GetName() << "\"" << flush;
314 TString curPath = fDigiScheme->GetCurrentPath();
315 CbmStsSensor* sensor = fDigiScheme->GetSensorByName(curPath);
316 Int_t stationNr = sensor->GetStationNr();
317
318 Int_t startHit = hitStationLimits[0][stationNr-1];
319 Int_t finalHit = hitStationLimits[1][stationNr-1];
320
321 // cout << " at station " << stationNr << " ( " << startHit << " - " << finalHit << " )" << endl;
322
323 fhEnergyLoss[stationNr-1]->Fill(stsPoint->GetXIn(),stsPoint->GetYIn(),stsPoint->GetEnergyLoss());
324
325 fhIncAngle[stationNr-1]->Fill(stsPoint->GetXIn(),stsPoint->GetYIn(),TMath::Abs(TMath::ATan2((stsPoint->GetZOut()-stsPoint->GetZIn()),(stsPoint->GetXOut()-stsPoint->GetXIn()))*180./3.14));
326 // Float_t zP = stsPoint->GetZ();
327 // Float_t xP = stsPoint->GetX(zP);
328 // Float_t yP = stsPoint->GetY(zP);
329
330 Bool_t foundPnt = kFALSE;
331
332 if ( startHit == -1 && finalHit == -1 ) continue;
333
334 for ( Int_t ihit = startHit ; ihit < finalHit ; ihit++ ) {
335 CbmStsHit *stsHit= (CbmStsHit*)fStsHits->At(ihit);
336 if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ())) < .01 ) &&
337 ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsHit->GetZ())) < .04 ) ) {
338 fhHitPointCorrelation[stationNr-1]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
339 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
340 fhHitPointCorr->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
341 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
342 fhHitPointPull -> Fill( (stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()))/stsHit->GetDx() , (stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()))/stsHit->GetDy());
343
344 if ( stsPoint->GetY(stsHit->GetZ()) > 0. )
345 fhHitPointCorrelationU[stationNr-1]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
346 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
347 if ( stsPoint->GetY(stsHit->GetZ()) < 0. )
348 fhHitPointCorrelationB[stationNr-1]->Fill(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ()),
349 stsHit->GetY()-stsPoint->GetY(stsHit->GetZ()));
350 }
351
352 if ( ( TMath::Abs(stsHit->GetX()-stsPoint->GetX(stsHit->GetZ())) < .01 ) &&
353 ( TMath::Abs(stsHit->GetY()-stsPoint->GetY(stsHit->GetZ())) < .04 ) ) {
354 foundPnt = kTRUE;
355 }
356 // cout <<"St "<<stationNr<<" Sens "<<sensor->GetSensorNr()<< " chan "<<sensor->GetFrontChannel(stsHit->GetX(),stsHit->GetY(),stsHit->GetZ())<<" X "<<stsHit->GetX();
357 }
358 Int_t incAngle = (Int_t) (TMath::ATan2((stsPoint->GetZOut()-stsPoint->GetZIn()),(stsPoint->GetXOut()-stsPoint->GetXIn()))*180./3.14);
359 // Int_t incAngle = (Int_t) ( TMath::ATan((stsPoint->GetXOut()-stsPoint->GetXIn())/(stsPoint->GetZOut()-stsPoint->GetZIn()))*180./3.14);
360 fNofPoints [stationNr-1][sensor->GetSectorNr()-1][sensor->GetSensorNr()-1] += 1;
361 fNofPointsIncAng[TMath::Abs((Int_t)(incAngle))] +=1;
362 if (TMath::Abs((Int_t)(mom))<1) fNofPointsMom [TMath::Abs((Int_t)(mom))] +=1;
363 if (TMath::Abs((Int_t)(mom))<1) fNofPointsMomSum+=1;
364 // cout << stationNr-1 << " " << sensor->GetSectorNr()-1 << " " << sensor->GetSensorNr()-1 << "/" << flush;
365 fhPoints [stationNr-1]->Fill(stsPoint->GetX(stsPoint->GetZ()),
366 stsPoint->GetY(stsPoint->GetZ()));
367 fNofPointsSec = 1;
368 if (trackPdg<10000&&motherId>=0) {
369 fNofPointsPdgSec[trackPdg] += 1;
370 fNofPointsSec +=1;
371 }
372 fNofPointsPrim = 1;
373 if (trackPdg<10000&&motherId<0) {
374 fNofPointsPdgPrim[trackPdg] += 1;
375 fNofPointsPrim +=1;
376 }
377 if ( foundPnt ) {
378 fNofRecoPoints[stationNr-1][sensor->GetSectorNr()-1][sensor->GetSensorNr()-1] += 1;
379 fNofRecoPointsIncAng[TMath::Abs((Int_t)(incAngle))] +=1;
380 if (TMath::Abs((Int_t)(mom))<1) fNofRecoPointsMom [TMath::Abs((Int_t)(mom))] +=1;
381 if (TMath::Abs((Int_t)(mom))<1) fNofRecoPointsMomSum+=1;
382 fhRecoPoints [stationNr-1]->Fill(stsPoint->GetX(stsPoint->GetZ()),
383 stsPoint->GetY(stsPoint->GetZ()));
384 fNofRecoPrim = 1;
385 if (trackPdg<10000&&motherId>=0) {
386 fNofRecoPdgSec[trackPdg] += 1;
387 fNofRecoSec +=1;
388 }
389 fNofRecoPrim = 1;
390 if (trackPdg<10000&&motherId<0) {
391 fNofRecoPdgPrim[trackPdg] += 1;
392 fNofRecoPrim +=1;
393 }
394 }
395
396 }
397 for (Int_t ang=0; ang<400; ang++) {
398 //Double_t factor = (Float_t)(fNofRecoPointsIncAng[ang])/(Float_t)(fNofPointsIncAng[ang]);
399 fhEffIncAng->Fill(ang,((Float_t)(fNofRecoPointsIncAng[ang]))/((Float_t)(fNofPointsIncAng[ang])));
400 }
401 for (Int_t momentum=0; momentum<50; momentum++) {
402 fhEffMom->Fill(momentum,((Float_t)(fNofRecoPointsMom[momentum]))/((Float_t)(fNofPointsMom[momentum])));
403 }
404
405 if ( fOnlineAnalysis ) {
406 recoPad[0]->cd();
407 TLegend* brp = new TLegend(0.1,0.1,0.9,0.9,"STS hit finding");
408 brp->SetTextAlign(22);
409 brp->SetTextSize(0.6);
410 brp->SetTextColor(1);
411 brp->SetBorderSize(0);
412 brp->SetFillColor(0);
413 brp->Draw();
414 recoPad[0]->Update();
415
416 recoPad[1]->cd();
417 fhEffIncAng->SetLineWidth(2);
418 fhEffIncAng->SetLineColor(1);
419 fhEffIncAng->SetMarkerColor(kGreen);
420 fhEffIncAng->SetMarkerStyle(2);
421 fhEffIncAng->SetTitle("Efficiency");
422 fhEffIncAng->Draw();
423 TLine* oneLine = new TLine(0.0,1.0,10.0,1.0);
424 oneLine->SetLineStyle(2);
425 oneLine->Draw();
426 recoPad[1]->Update();
427
428 TH1F* projX[10];
429 TH1F* projY[10];
430 TF1* fitX[10];
431 TF1* fitY[10];
432 Double_t resolution[2][10];
433
434 for ( Int_t ist = 0 ; ist < 8 ; ist++ ) {
435 projX[ist] = (TH1F*)fhHitPointCorrelation[ist]->ProjectionX(Form("projX%i",ist+1));
436 projY[ist] = (TH1F*)fhHitPointCorrelation[ist]->ProjectionY(Form("projY%i",ist+1));
437 fitX[ist] = new TF1(Form("fitX%i",ist+1),"gaus",-0.02,0.02);
438 fitY[ist] = new TF1(Form("fitY%i",ist+1),"gaus",-0.02,0.02);
439 projX[ist]->SetAxisRange(-0.02,0.02,"X");
440 projY[ist]->SetAxisRange(-0.02,0.02,"X");
441 projX[ist]->Fit(fitX[ist],"QN","",-0.02,0.02);
442 projY[ist]->Fit(fitY[ist],"QN","",-0.02,0.02);
443 resolution[0][ist] = 10000.*fitX[ist]->GetParameter(2);
444 resolution[1][ist] = 10000.*fitY[ist]->GetParameter(2);
445 }
446
447 TH1F* projXPull;
448 TH1F* projYPull;
449 TF1* fitXPull;
450 TF1* fitYPull;
451 Double_t resolutionPull[2];
452
453
454 projXPull = (TH1F*)fhHitPointPull->ProjectionX("projXPull");
455 projYPull = (TH1F*)fhHitPointPull->ProjectionY("projYPull");
456 fitXPull = new TF1("fitXPull","gaus",-2,2);
457 fitYPull = new TF1("fitYPull","gaus",-2,2);
458 projXPull->SetAxisRange(-2,2,"X");
459 projYPull->SetAxisRange(-2,2,"X");
460 projXPull->Fit(fitXPull,"QN","",-2,2);
461 projYPull->Fit(fitYPull,"QN","",-2,2);
462 resolutionPull[0] = fitXPull->GetParameter(2);
463 resolutionPull[1] = fitYPull->GetParameter(2);
464
465
466 recoPad[3]->cd();
467 fhHitPointCorrelation[0]->Draw("col");
468 recoPad[3]->Update();
469
470 recoPad[4]->cd();
471 fhHitPointPull->Draw("col");
472 recoPad[4]->Update();
473
474 recoPad[5]->cd();
475 projX[0]->SetLineWidth(2);
476 projY[0]->SetLineWidth(2);
477 projX[0]->SetLineColor(2);
478 projY[0]->SetLineColor(3);
479 projY[0]->SetXTitle("#Delta x, y [cm]");
480 projX[0]->Draw();
481 projY[0]->Draw("same");
482 // fitY[fShowStation1]->Draw("same");
483 // fitX[fShowStation1]->Draw("same");
484 TLegend* legend1 = new TLegend(0.55,0.6,1.0,0.94);
485 legend1->SetBorderSize(0);
486 legend1->SetFillColor(0);
487 legend1->AddEntry(projX[0],
488 Form("X,#sigma=%3.2f#mum",resolution[0][0]),"l");
489 legend1->AddEntry(projY[0],
490 Form("Y,#sigma=%3.2f#mum",resolution[1][0]),"l");
491 legend1->Draw();
492 recoPad[5]->Update();
493
494 recoPad[6]->cd();
495 projXPull->SetLineWidth(2);
496 projYPull->SetLineWidth(2);
497 projXPull->SetLineColor(2);
498 projYPull->SetLineColor(3);
499 projXPull->SetXTitle("Pull x, y ");
500 projXPull->Draw();
501 projYPull->Draw("same");
502 // fitY[fShowStation2]->Draw("same");
503 // fitX[fShowStation2]->Draw("same");
504 TLegend* legend2a = new TLegend(0.55,0.6,1.0,0.94);
505 legend2a->SetBorderSize(0);
506 legend2a->SetFillColor(0);
507 legend2a->AddEntry(projXPull,Form("X,#sigma=%3.2f#mum",resolutionPull[0]),"l");
508 legend2a->AddEntry(projYPull,Form("Y,#sigma=%3.2f#mum",resolutionPull[1]),"l");
509 legend2a->Draw();
510 recoPad[6]->Update();
511 recoPad[2]->cd();
512 fhRecoPoints[0]->SetMarkerColor(kCyan);
513 fhRecoPoints[0]->Draw();
514 recoPad[2]->Update();
515
516 recoPad[7]->cd();
517 fhRecoPoints[1]->SetMarkerColor(kMagenta);
518 fhRecoPoints[1]->Draw();
519 recoPad[7]->Update();
520
521 recoPad[8]->cd();
522 fhRecoPoints[2]->SetMarkerColor(kOrange);
523 fhRecoPoints[2]->Draw();
524 recoPad[8]->Update();
525
526 recoPad[9]->cd();
527 fhRecoPoints[4]->SetMarkerColor(kGreen);
528 fhRecoPoints[4]->Draw();
529 recoPad[9]->Update();
530
531 recoPad[10]->cd();
532 fhRecoPoints[5]->SetMarkerColor(kBlue);
533 fhRecoPoints[5]->Draw();
534 recoPad[10]->Update();
535
536 recoPad[11]->cd();
537 fhRecoPoints[7]->SetMarkerColor(kRed);
538 fhRecoPoints[7]->Draw();
539 recoPad[11]->Update();
540
541 recoPad[12]->cd();
542 TPaveText* printoutPave = new TPaveText(0.0,0.0,1.0,1.0);
543 printoutPave->SetTextAlign(20);
544 printoutPave->SetTextSize(0.08);
545 printoutPave->SetTextColor(1);
546 printoutPave->SetBorderSize(0);
547 printoutPave->SetFillColor(0);
548 printoutPave->AddText(Form("%i events",fNEvents+1));
549 printoutPave->AddText(Form("Hits/Points %3.2f ",Double_t (nofStsHits)/Double_t (nofStsPoints)));
550
551 printoutPave->AddText(Form("%3.0f points prim, %3.0f reco prim",
552 Double_t (fNofPointsPrim)/Double_t (fNEvents+1),
553 Double_t (fNofRecoPrim)/Double_t (fNEvents+1)));
554 printoutPave->AddText(Form("%3.0f points sec, %3.0f reco sec",
555 Double_t (fNofPointsSec)/Double_t (fNEvents+1),
556 Double_t (fNofRecoSec)/Double_t (fNEvents+1)));
557
558 printoutPave->AddText(Form("Hit Finding Efficiency: %3.2f prim, %3.2f sec",
559 (Double_t (fNofRecoPrim)/Double_t (fNofPointsPrim))*100.,
560 (Double_t (fNofRecoSec) /Double_t (fNofPointsSec))*100.));
561 // printoutPave->AddText("Tracking efficiencies (p>1.0 GeV/c):");
562 // printoutPave->AddText(Form("all = %2.2f%%(%2.2f%%)",100.*effAll,100.*allEff));
563 // printoutPave->AddText(Form("vertex = %2.2f%%(%2.2f%%)",100.*effPrim,100.*primEff));
564 // printoutPave->AddText(Form("reference = %2.2f%%",100.*effRef));
565 // printoutPave->AddText(Form("non-vertex = %2.2f%%(%2.2f%%)",100.*effSec,100.*secEff));
566 // printoutPave->AddText(Form("Momentum resolution = %3.2f%%(%3.2f%%)",momentumResolutionAll,momentumResolutionPrim));
567 recoPad[12]->Clear();
568 printoutPave->Draw();
569 recoPad[12]->Update();
570
571 }
572
573 // Event summary
574 if ( fVerbose ) {
575 cout << "---------- StsFindHitsQa : Event " << fNEvents+1 << " summary ------------"
576 << endl;
577 cout << "-----------------------------------------------------------"
578 << endl;
579 cout << endl;
580 }
581 else
582 cout << "\r+ " << setw(15) << left << fName << ": event " << fNEvents+1 << " " << setprecision(4)
583 << setw(8) << fixed << right << fTimer.RealTime() << "s." << endl;
584
585 //cout << "-------- " << fNofRecoPoints[1][4][0] << " / " << fNofPoints[1][4][0] << endl;
586
587 fNEvents++;
588 fTime1 += fTimer.RealTime();
589}
590// -------------------------------------------------------------------------
591
592
593// ----- Private method CreateHistos -----------------------------------
595
596
597 // Histogram list
598 fHistoList = new TList();
599 for ( Int_t ist = 0 ; ist < fNStations ; ist++ )
600 fHistoListPS[ist] = new TList();
601
602 // Momentum distributions
603 //Double_t minMom = 0.;
604 //Double_t maxMom = 10.;
605 //Int_t nBinsMom = 40;
606
607 CbmStsStation* station = NULL;
608 CbmStsSector* sector = NULL;
609 //CbmStsSensor* sensor = NULL;
610
611 fhHitFindingEfficiency = new TH1F("hHitFindingEfficiency","Hit finding efficiency",1000000,
612 (Float_t)(2<<24)-0.5,(Float_t)(2<<24)-0.5+1000000.);
613 fhHitPointCorr = new TH2F(Form("hHitPointCorrelation"),
614 Form("Hit vs point correlation "),
615 500,-.1, .1,500,-.1,.1);
616 fhEffMom = new TH2F("hEffMom","hEffMom",500,0,25,100,0,1);
617 fhEffIncAng = new TH2F("hEffIncAng","hEffIncAng",100,0,100,100,0,1);
618 fhEffPdgSec = new TH2F("hEffPdgSec", "hEffPdgSec", 350,-1000,2500,10,0,1);
619 fhEffPdgPrim = new TH2F("hEffPdgPrim","hEffPdgPrim",350,-1000,2500,10,0,1);
620 fhHitPointPull = new TH2F("hHitPointPull","Hit vs sigma at station %i",500,-3., 3.,500,-3.,3.);
621
622 fHistoList->Add(fhEffPdgSec);
623 fHistoList->Add(fhEffPdgPrim);
624 fHistoList->Add(fhHitFindingEfficiency);
625 fHistoList->Add(fhEffIncAng);
626 fHistoList->Add(fhHitPointPull);
627 fHistoList->Add(fhEffMom);
628 fHistoList->Add(fhHitPointCorr);
629 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
630 fhHitPointCorrelation[ist] = new TH2F(Form("hHitPointCorrelation%i",ist+1),
631 Form("Hit vs point correlation at station %i",ist+1),
632 500,-.1, .1,500,-.1,.1);
633 fhHitPointCorrelation[ist]->SetXTitle("#Delta x [cm]");
634 fhHitPointCorrelation[ist]->SetYTitle("#Delta y [cm]");
635
636
637
638
639 fHistoListPS[ist]->Add(fhHitPointCorrelation[ist]);
640
641
642 fhHitPointCorrelationU[ist] = new TH2F(Form("hHitPointCorrelationU%i",ist+1),
643 Form("Hit vs point correlation at station %i",ist+1),
644 500,-.1, .1,500,-.1,.1);
645 fhHitPointCorrelationB[ist] = new TH2F(Form("hHitPointCorrelationB%i",ist+1),
646 Form("Hit vs point correlation at station %i",ist+1),
647 500,-.1, .1,500,-.1,.1);
648 fHistoListPS[ist]->Add(fhHitPointCorrelationU[ist]);
649 fHistoListPS[ist]->Add(fhHitPointCorrelationB[ist]);
650
651 fhPoints [ist] = new TH2F(Form("hPoints%i",ist+1),Form("Points on station %i",ist+1),
652 100,-100.,100.,
653 100,-100.,100.);
654 fhRecoPoints[ist] = new TH2F(Form("hRecoPoints%i",ist+1),Form("Reconstructed points on station %i",ist+1),
655 100,-100.,100.,
656 100,-100.,100.);
657 fHistoListPS[ist]->Add(fhPoints [ist]);
658 fHistoListPS[ist]->Add(fhRecoPoints[ist]);
659
660 }
661
662 Double_t binningNofHits[501];
663 binningNofHits[0] = -0.5;
664 for ( Int_t itemp = 0 ; itemp < 100 ; itemp++ ) {
665 binningNofHits[itemp+ 1] = 0.5+(Double_t)itemp;
666 binningNofHits[itemp+101] = 101.5+2.*(Double_t)itemp;
667 binningNofHits[itemp+201] = 302.5+3.*(Double_t)itemp;
668 binningNofHits[itemp+301] = 602.5+3.*(Double_t)itemp;
669 binningNofHits[itemp+401] = 902.5+3.*(Double_t)itemp;
670 }
671
672 //Double_t fSectorWidth = 0.;
673 for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
674 station = fDigiScheme->GetStation(istat);
675
676 fhEnergyLoss[istat] = new TH3F(Form("hEnergyLossSt%d",istat+1),
677 Form("Energy loss on station %d",istat+1),
678 20,-100.,100.,20,-100.,100.,100,0,0.01);
679 fhIncAngle [istat] = new TH3F(Form("hIncAngleSt%d",istat+1),
680 Form("Incident Angle on station %d;X;Y;Angle",istat+1),
681 100,-50,50,100,-40,40,20,0,100.);
682 fHistoListPS[istat]->Add(fhEnergyLoss[istat]);
683 fHistoListPS[istat]->Add(fhIncAngle[istat]);
684 for ( Int_t isect = 0 ; isect < station->GetNSectors() ; isect++ ) {
685 sector = station->GetSector(isect);
686 for ( Int_t isens = 0 ; isens < sector->GetNSensors() ; isens++ ) {
687 //sensor = sector->GetSensor(isens);
688 //fSectorWidth = 10.*sensor->GetLx();
689 fhNofHits[istat][isect][isens] = new TH1F(Form("hNofHitsSt%dSect%dSens%d",istat+1,isect+1,isens+1),
690 Form("Number of hits in sector %d, sensor %d of station %d",isect+1,isens+1,istat+1),
691 500,binningNofHits);
692 // 500,-0.5,500.5);
693 fHistoListPS[istat]->Add(fhNofHits[istat][isect][isens]);
694 }
695 }
696 }
697
698}
699// -------------------------------------------------------------------------
700
701// ----- Private method Reset ------------------------------------------
703
704
705 TIter next(fHistoList);
706 while ( TH1* histo = ((TH1*)next()) ) histo->Reset();
707 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
708 TIter next0(fHistoListPS[ist]);
709 while ( TH1* histo = ((TH1*)next0()) ) histo->Reset();
710 }
711
712}
713// -------------------------------------------------------------------------
714
715// ----- Private method SetParContainers -------------------------------
716void CbmStsFindHitsQa::SetParContainers() {
717
718
719 // Get run and runtime database
720 FairRunAna* run = FairRunAna::Instance();
721 if ( ! run ) Fatal("SetParContainers", "No analysis run");
722
723 FairRuntimeDb* db = run->GetRuntimeDb();
724 if ( ! db ) Fatal("SetParContainers", "No runtime database");
725
726 // Get STS geometry parameter container
727 fGeoPar = (CbmGeoStsPar*) db->getContainer("CbmGeoStsPar");
728
729 // Get STS digitisation parameter container
730 fDigiPar = (CbmStsDigiPar*) db->getContainer("CbmStsDigiPar");
731
732}
733// -------------------------------------------------------------------------
734
735
736// ----- Private method Init -------------------------------------------
737InitStatus CbmStsFindHitsQa::Init() {
738
739
740 // Get input array
741 FairRootManager* ioman = FairRootManager::Instance();
742 if ( ! ioman ) Fatal("Init", "No FairRootManager");
743
744 // Get MCTrack array
745 fMCTracks = (TClonesArray*) ioman->GetObject("MCTrack");
746 if ( ! fMCTracks ) {
747 cout << "-E- " << GetName() << "::Init: No MCTrack array!" << endl;
748 return kFATAL;
749 }
750 // Get StsPoint array
751 fStsPoints = (TClonesArray*) ioman->GetObject("StsPoint");
752 if ( ! fStsPoints ) {
753 //cout << "-E- " << GetName() << "::Init: No StsPoint array!" << endl;
754 return kFATAL;
755 }
756
757 // Get StsHit array
758 fStsHits = (TClonesArray*) ioman->GetObject("StsHit");
759 if ( ! fStsHits ) {
760 cout << "-E- " << GetName() << "::Init: No StsHit array!" << endl;
761 return kFATAL;
762 }
763
764 // Get StsDigis array
765 fStsDigis = (TClonesArray*) ioman->GetObject("StsDigi");
766 if ( ! fStsDigis ) {
767 cout << "-E- " << GetName() << "::Init: No StsDigi array!"
768 << endl;
769 return kERROR;
770 }
771
772 // Get StsClusters array
773 fStsClusters = (TClonesArray*) ioman->GetObject("StsCluster");
774 if ( ! fStsClusters ) {
775 cout << "-w- " << GetName() << "::Init: No StsCluster array!"
776 << endl;
777 // return kERROR;
778 }
779
780 // Build digitisation scheme
781 Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
782 if ( ! success ) return kERROR;
783
784 fNStations = fDigiScheme->GetNStations();
785
786 CbmStsStation* station = NULL;
787 CbmStsSector* sector = NULL;
788 for ( Int_t istat = fNStations ; istat > 0 ; ) {
789 station = fDigiScheme->GetStation(--istat);
790 for ( Int_t isect = station->GetNSectors() ; isect > 0 ; ) {
791 sector = station->GetSector(--isect);
792 for ( Int_t isens = sector->GetNSensors() ; isens > 0 ; ) {
793 fNofPoints [istat][isect][--isens] = 0;
794 fNofRecoPoints[istat][isect][ isens] = 0;
795 }
796 }
797 }
798 for (Int_t iang=0; iang<400; iang++) {
799 fNofPointsIncAng[iang] = 0;
800 fNofRecoPointsIncAng[iang] = 0;
801 }
802 for (Int_t ipdg=-10000; ipdg<10000; ipdg++) {
803 fNofRecoPdgSec [ipdg] = 0;
804 fNofPointsPdgSec [ipdg] = 0;
805 fNofRecoPdgPrim [ipdg] = 0;
806 fNofPointsPdgPrim [ipdg] = 0;
807 }
808
809 CreateHistos();
810 Reset();
811 if ( fOnlineAnalysis ) {
812 recoCanvas = new TCanvas("StsRecoCanvas","Sts reconstruction",10,10,600,900);
813 // TPad* recoPad[10];
814
815 recoPad[0] = new TPad("titlePad", "Title pad" ,0.00,0.90,1.00,1.00);
816 recoPad[1] = new TPad("efficiencyPad","Efficiency pad" ,0.00,0.60,0.40,0.90);
817 recoPad[3] = new TPad("hitVspoint1Pad","Hit vs point on st.1 pad" ,0.40,0.65,0.70,0.90);
818 recoPad[4] = new TPad("hitPullst1","Hit vs sigma on st.1 pad" ,0.70,0.65,1.00,0.90);
819 recoPad[5] = new TPad("hitVspntXY1Pad","Hit vs point XY on st.1 pad" ,0.40,0.40,0.70,0.65);
820 recoPad[6] = new TPad("hitVssigmaXY6Pad","Hit vs sigma XY on st.1 pad",0.70,0.40,1.00,0.65);
821 recoPad[2] = new TPad("Reco Points st1","Reco Points st1" ,0.00,0.20,0.33,0.40);
822 recoPad[7] = new TPad("Reco Points st2","Reco Points st2" ,0.33,0.20,0.66,0.40);
823 recoPad[8] = new TPad("Reco Points st3","Reco Points st3" ,0.66,0.20,1.00,0.40);
824 recoPad[9] = new TPad("Reco Points st5","Reco Points st5" ,0.00,0.00,0.33,0.20);
825 recoPad[10] = new TPad("Reco Points st6","Reco Points st6" ,0.33,0.00,0.66,0.20);
826 recoPad[11] = new TPad("Reco Points st8","Reco Points st8" ,0.66,0.00,1.00,0.20);
827 recoPad[12] = new TPad("printoutPad","Print information pad" ,0.00,0.40,0.40,0.60);
828 // recoPad[12]->SetLogz();
829 for ( Int_t ipad = 0 ; ipad < 13 ; ipad++ ) {
830 recoPad[ipad]->SetFillColor(0);
831 recoPad[ipad]->SetBorderMode(0);
832 recoPad[ipad]->Draw();
833 }
834 }
835 // Control output
836 if (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
837 else if (fVerbose > 2) fDigiScheme->Print(kTRUE);
838 cout << "-I- " << fName << "::Init: " << "STS digitisation scheme succesfully initialised" << endl;
839 cout << " Stations: " << fDigiScheme->GetNStations()
840 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: "
841 << fDigiScheme->GetNChannels() << endl;
842
843 return kSUCCESS;
844
845}
846// -------------------------------------------------------------------------
847
848// ----- Private method ReInit -----------------------------------------
849InitStatus CbmStsFindHitsQa::ReInit()
850{
851 // Clear digitisation scheme
852 fDigiScheme->Clear();
853
854 // Build new digitisation scheme
855 Bool_t success = fDigiScheme->Init(fGeoPar, fDigiPar);
856 if ( ! success ) return kERROR;
857
858 // Control output
859 if (fVerbose == 1 || fVerbose == 2) fDigiScheme->Print(kFALSE);
860 else if (fVerbose > 2) fDigiScheme->Print(kTRUE);
861 cout << "-I- " << fName << "::Init: "
862 << "STS digitisation scheme succesfully reinitialised" << endl;
863 cout << " Stations: " << fDigiScheme->GetNStations()
864 << ", Sectors: " << fDigiScheme->GetNSectors() << ", Channels: "
865 << fDigiScheme->GetNChannels() << endl;
866
867 return kSUCCESS;
868
869}
870// -------------------------------------------------------------------------
871
872// ----- Virtual method Finish -----------------------------------------
874
875 CbmStsStation* station = NULL;
876 CbmStsSector* sector = NULL;
877 //Int_t NofRecoPoints = 0;
878 Int_t HFE_reco = 0;
879 Int_t HFE_point = 0;
880 for ( Int_t istat = 0 ; istat < fDigiScheme->GetNStations() ; istat++ ) {
881 station = fDigiScheme->GetStation(istat);
882 for ( Int_t isect = 0 ; isect < station->GetNSectors() ; isect++ ) {
883 sector = station->GetSector(isect);
884 for ( Int_t isens = 0 ; isens < sector->GetNSensors() ; isens++ ) {
885 HFE_reco = 1;
886 HFE_point = 1;
887 if ( fNofPoints[istat][isect][isens] ) {
888 fhHitFindingEfficiency->Fill(2 << 24 | (istat+1) << 16 | (isect+1) << 4 | (isens+1) << 1,
889 ((Float_t)(fNofRecoPoints[istat][isect][isens]))/
890 ((Float_t)(fNofPoints[istat][isect][isens])));
891 HFE_reco += (fNofRecoPoints[istat][isect][isens]);
892 HFE_point += (fNofPoints[istat][isect][isens]);
893 }
894 }
895 }
896 }
897 for (Int_t iPdg=-10000; iPdg<10000; iPdg++) {
898 fhEffPdgSec -> Fill(iPdg, ((Float_t)(fNofRecoPdgSec[iPdg])) / ((Float_t)(fNofPointsPdgSec[iPdg])));
899 fhEffPdgPrim -> Fill(iPdg, ((Float_t)(fNofRecoPdgPrim[iPdg])) / ((Float_t)(fNofPointsPdgPrim[iPdg])));
900 }
901
902 cout << endl;
903 cout << "============================================================"
904 << endl;
905 cout << "===== " << fName << ": Run summary " << endl;
906 cout << "===== Hit finding efficiency = "<<(((Float_t)(HFE_reco))/((Float_t)(HFE_point)))*100.<<"%"<< endl;
907 cout << "============================================================"
908 << endl;
909
910 gDirectory->mkdir("STSFindHitsQA");
911 gDirectory->cd("STSFindHitsQA");
912 TIter next(fHistoList);
913 while ( TH1* histo = ((TH1*)next()) ) histo->Write();
914
915 for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
916 gDirectory->mkdir(Form("Station%d",istat+1));
917 gDirectory->cd(Form("Station%d",istat+1));
918 TIter nextO(fHistoListPS[istat]);
919 while ( TH1* histo = ((TH1*)nextO()) ) histo->Write();
920
921 gDirectory->cd("..");
922 }
923 gDirectory->cd("..");
924}
925// -------------------------------------------------------------------------
BmnSsdDigitizeParameters * fDigiPar
Digitisation parameters.
void GetMomentum(TVector3 &momentum) const
Definition CbmMCTrack.h:139
void GetStartVertex(TVector3 &vertex) const
Definition CbmMCTrack.h:149
void Print(Bool_t kLong=kFALSE)
static CbmStsDigiScheme * Instance(int version=1)
CbmStsSensor * GetSensorByName(TString sensorName)
CbmStsStation * GetStation(Int_t iStation)
Int_t GetSide() const
Definition CbmStsDigi.h:79
Int_t GetSectorNr() const
Definition CbmStsDigi.h:75
Int_t GetChannelNr() const
Definition CbmStsDigi.h:83
Int_t GetStationNr() const
Definition CbmStsDigi.h:72
virtual void Exec(Option_t *opt)
virtual void Finish()
Int_t GetSensorNr() const
Definition CbmStsHit.h:70
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:66
Int_t GetSectorNr() const
Definition CbmStsHit.h:68
Double_t GetY(Double_t z) const
Double_t GetX(Double_t z) const
Double_t GetXIn() const
Definition CbmStsPoint.h:69
Double_t GetZOut() const
Definition CbmStsPoint.h:74
Double_t GetZIn() const
Definition CbmStsPoint.h:71
Double_t GetYIn() const
Definition CbmStsPoint.h:70
Double_t GetXOut() const
Definition CbmStsPoint.h:72
CbmStsSensor * GetSensor(Int_t iSensor)
Int_t GetNSensors() const
Double_t GetX0() const
Int_t GetSectorNr() const
Int_t GetStationNr() const
Double_t GetY0() const
Int_t GetSensorNr() const
Int_t GetNSectors() const
CbmStsSector * GetSector(Int_t iSector)
-clang-format
STL namespace.