211 if ( ! fDigiScheme ) {
212 cerr <<
"-E- " << fName <<
"::Exec: No digi scheme!" << endl;
219 Int_t nofStsDigis = fStsDigis->GetEntriesFast();
220 for ( Int_t idigi = 0 ; idigi < nofStsDigis ; idigi++ ) {
226 Int_t nofStsHits = fStsHits->GetEntriesFast();
227 Int_t nofStsPoints = fStsPoints->GetEntriesFast();
228 Int_t hitStationLimits[2][100];
230 for ( Int_t ist = 0 ; ist < fNStations ; ist++ ) {
231 hitStationLimits[0][ist] = -1;
232 hitStationLimits[1][ist] = -1;
237 for ( Int_t ihit = 0 ; ihit < nofStsHits ; ihit++ ) {
240 if ( hitStationLimits[0][stsHit->
GetStationNr()-1] == -1 ) {
244 if ( hitStationLimits[1][stsHitBack->
GetStationNr()-1] == -1 ) {
245 hitStationLimits[1][stsHitBack->
GetStationNr()-1] = nofStsHits-ihit;
248 for ( Int_t istat = 0 ; istat < fDigiScheme->
GetNStations() ; istat++ ) {
250 for ( Int_t isect = 0 ; isect < station->
GetNSectors() ; isect++ ) {
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;
260 for ( Int_t isens = 0 ; isens < sector->
GetNSensors() ; isens++ ) {
262 fhNofHits[istat][isect][isens] -> Fill(fNofHits[istat][isect][isens]);
276 for ( Int_t ipnt = 0 ; ipnt < nofStsPoints ; ipnt++ ) {
284 TVector3 trackVertex;
285 TVector3 motherVertex;
288 Double_t xin = stsPoint->
GetXIn();
289 Double_t yin = stsPoint->
GetYIn();
290 Double_t zin = stsPoint->
GetZIn();
295 gGeoManager->FindNode(xin,yin,zin);
298 Int_t trackId = stsPoint->GetTrackID();
300 trackPdg = mcTrack -> GetPdgCode();
301 motherId = mcTrack -> GetMotherId();
303 Double_t mom = momentum.Mag();
318 Int_t startHit = hitStationLimits[0][stationNr-1];
319 Int_t finalHit = hitStationLimits[1][stationNr-1];
323 fhEnergyLoss[stationNr-1]->Fill(stsPoint->
GetXIn(),stsPoint->
GetYIn(),stsPoint->GetEnergyLoss());
325 fhIncAngle[stationNr-1]->Fill(stsPoint->
GetXIn(),stsPoint->
GetYIn(),TMath::Abs(TMath::ATan2((stsPoint->
GetZOut()-stsPoint->
GetZIn()),(stsPoint->
GetXOut()-stsPoint->
GetXIn()))*180./3.14));
330 Bool_t foundPnt = kFALSE;
332 if ( startHit == -1 && finalHit == -1 )
continue;
334 for ( Int_t ihit = startHit ; ihit < finalHit ; 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());
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()));
352 if ( ( TMath::Abs(stsHit->GetX()-stsPoint->
GetX(stsHit->GetZ())) < .01 ) &&
353 ( TMath::Abs(stsHit->GetY()-stsPoint->
GetY(stsHit->GetZ())) < .04 ) ) {
358 Int_t incAngle = (Int_t) (TMath::ATan2((stsPoint->
GetZOut()-stsPoint->
GetZIn()),(stsPoint->
GetXOut()-stsPoint->
GetXIn()))*180./3.14);
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;
365 fhPoints [stationNr-1]->Fill(stsPoint->
GetX(stsPoint->GetZ()),
366 stsPoint->
GetY(stsPoint->GetZ()));
368 if (trackPdg<10000&&motherId>=0) {
369 fNofPointsPdgSec[trackPdg] += 1;
373 if (trackPdg<10000&&motherId<0) {
374 fNofPointsPdgPrim[trackPdg] += 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()));
385 if (trackPdg<10000&&motherId>=0) {
386 fNofRecoPdgSec[trackPdg] += 1;
390 if (trackPdg<10000&&motherId<0) {
391 fNofRecoPdgPrim[trackPdg] += 1;
397 for (Int_t ang=0; ang<400; ang++) {
399 fhEffIncAng->Fill(ang,((Float_t)(fNofRecoPointsIncAng[ang]))/((Float_t)(fNofPointsIncAng[ang])));
401 for (Int_t momentum=0; momentum<50; momentum++) {
402 fhEffMom->Fill(momentum,((Float_t)(fNofRecoPointsMom[momentum]))/((Float_t)(fNofPointsMom[momentum])));
405 if ( fOnlineAnalysis ) {
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);
414 recoPad[0]->Update();
417 fhEffIncAng->SetLineWidth(2);
418 fhEffIncAng->SetLineColor(1);
419 fhEffIncAng->SetMarkerColor(kGreen);
420 fhEffIncAng->SetMarkerStyle(2);
421 fhEffIncAng->SetTitle(
"Efficiency");
423 TLine* oneLine =
new TLine(0.0,1.0,10.0,1.0);
424 oneLine->SetLineStyle(2);
426 recoPad[1]->Update();
432 Double_t resolution[2][10];
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);
451 Double_t resolutionPull[2];
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);
467 fhHitPointCorrelation[0]->Draw(
"col");
468 recoPad[3]->Update();
471 fhHitPointPull->Draw(
"col");
472 recoPad[4]->Update();
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]");
481 projY[0]->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");
492 recoPad[5]->Update();
495 projXPull->SetLineWidth(2);
496 projYPull->SetLineWidth(2);
497 projXPull->SetLineColor(2);
498 projYPull->SetLineColor(3);
499 projXPull->SetXTitle(
"Pull x, y ");
501 projYPull->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");
510 recoPad[6]->Update();
512 fhRecoPoints[0]->SetMarkerColor(kCyan);
513 fhRecoPoints[0]->Draw();
514 recoPad[2]->Update();
517 fhRecoPoints[1]->SetMarkerColor(kMagenta);
518 fhRecoPoints[1]->Draw();
519 recoPad[7]->Update();
522 fhRecoPoints[2]->SetMarkerColor(kOrange);
523 fhRecoPoints[2]->Draw();
524 recoPad[8]->Update();
527 fhRecoPoints[4]->SetMarkerColor(kGreen);
528 fhRecoPoints[4]->Draw();
529 recoPad[9]->Update();
532 fhRecoPoints[5]->SetMarkerColor(kBlue);
533 fhRecoPoints[5]->Draw();
534 recoPad[10]->Update();
537 fhRecoPoints[7]->SetMarkerColor(kRed);
538 fhRecoPoints[7]->Draw();
539 recoPad[11]->Update();
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)));
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)));
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.));
567 recoPad[12]->Clear();
568 printoutPave->Draw();
569 recoPad[12]->Update();
575 cout <<
"---------- StsFindHitsQa : Event " << fNEvents+1 <<
" summary ------------"
577 cout <<
"-----------------------------------------------------------"
582 cout <<
"\r+ " << setw(15) << left << fName <<
": event " << fNEvents+1 <<
" " << setprecision(4)
583 << setw(8) << fixed << right << fTimer.RealTime() <<
"s." << endl;
588 fTime1 += fTimer.RealTime();
598 fHistoList =
new TList();
599 for ( Int_t ist = 0 ; ist < fNStations ; ist++ )
600 fHistoListPS[ist] =
new TList();
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.);
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]");
639 fHistoListPS[ist]->Add(fhHitPointCorrelation[ist]);
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]);
651 fhPoints [ist] =
new TH2F(Form(
"hPoints%i",ist+1),Form(
"Points on station %i",ist+1),
654 fhRecoPoints[ist] =
new TH2F(Form(
"hRecoPoints%i",ist+1),Form(
"Reconstructed points on station %i",ist+1),
657 fHistoListPS[ist]->Add(fhPoints [ist]);
658 fHistoListPS[ist]->Add(fhRecoPoints[ist]);
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;
673 for ( Int_t istat = 0 ; istat < fNStations ; istat++ ) {
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++ ) {
686 for ( Int_t isens = 0 ; isens < sector->
GetNSensors() ; isens++ ) {
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),
693 fHistoListPS[istat]->Add(fhNofHits[istat][isect][isens]);