3 #include "FairLogger.h"
4 #include "FairRootManager.h"
6 #include "TClonesArray.h"
27 #include "FairTrackParam.h"
28 #include "FairVolume.h"
30 #include "TGeoManager.h"
41 #include "TGeoSphere.h"
42 #include "TVirtualMC.h"
53 , fRichProjections(NULL)
54 , fRichMirrorPoints(NULL)
56 , fRichRingMatches(NULL)
67 , fDrawAlignment(kTRUE)
75 FairRootManager* manager = FairRootManager::Instance();
77 fRichHits = (TClonesArray*) manager->GetObject(
"RichHit");
79 Fatal(
"CbmRichAlignment::Init",
"No RichHit array !");
82 fRichRings = (TClonesArray*) manager->GetObject(
"RichRing");
84 Fatal(
"CbmRichAlignment::Init",
"No RichRing array !");
89 Fatal(
"CbmRichAlignment::Init",
"No RichProjection array !");
92 fRichPoints = (TClonesArray*) manager->GetObject(
"RichPoint");
94 Fatal(
"CbmRichAlignment::Init",
"No RichPoint array !");
97 fMCTracks = (TClonesArray*) manager->GetObject(
"MCTrack");
99 Fatal(
"CbmRichAlignment::Init",
"No MCTracks array !");
104 Fatal(
"CbmRichAlignment::Init",
"No RichRingMatches array !");
109 Fatal(
"CbmRichAlignment::Init",
"No RichMirrorPoints array !");
135 "fHCenterDistance;Distance C-C';Nb of entries",
140 "fHPhi",
"fHPhi;Phi_Ch [rad];Nb of entries", 200, -3.4, 3.4);
142 "fHThetaDiff",
"fHThetaDiff;Th_Ch-Th_0 [cm];Nb of entries", 252, -5., 5.);
144 "fHCherenkovHitsDistribTheta0",
145 "fHCherenkovHitsDistribTheta0;Phi_0 [rad];Theta_0 [cm];Entries",
153 "fHCherenkovHitsDistribThetaCh",
154 "fHCherenkovHitsDistribThetaCh;Phi_Ch [rad];Theta_Ch [cm];Entries",
162 "fHCherenkovHitsDistribReduced",
163 "fHCherenkovHitsDistribReduced;Phi_Ch [rad];Th_Ch-Th_0 [cm];Entries",
175 "--------------------------------------------------------------------"
176 "-----------------------------------------------//"
178 cout <<
"//------------------------------------------ CbmRichAlignment: EXEC "
179 "Function ------------------------------------------//"
182 "--------------------------------------------------------------------"
183 "--------------------------------------------------//"
187 cout <<
"CbmRichAlignment : Event #" <<
fEventNum << endl;
189 Int_t nofRingsInEvent =
fRichRings->GetEntries();
191 Int_t nofHitsInEvent =
fRichHits->GetEntries();
192 Int_t NofMCTracks =
fMCTracks->GetEntriesFast();
194 cout <<
"Nb of rings in evt = " << nofRingsInEvent
195 <<
", nb of mirror points = " << nofMirrorPoints
196 <<
", nb of hits in evt = " << nofHitsInEvent
197 <<
" and nb of Monte-Carlo tracks = " << NofMCTracks << endl
200 if (nofRingsInEvent == 0) {
201 cout <<
"Error no rings registered in event." << endl << endl;
208 cout <<
"//------------------------------ CbmRichAlignment: Calculate Angles "
209 "& Draw Distrib ------------------------------//"
213 Double_t trackX = 0., trackY = 0.;
216 Int_t nofRingsInEvent =
fRichRings->GetEntries();
217 Float_t DistCenters, Theta_Ch, Theta_0, Angles_0;
218 Float_t
Pi = 3.14159265;
219 Float_t TwoPi = 2. * 3.14159265;
223 if (nofRingsInEvent >= 1) {
224 cout <<
"Number of Rings in event: " << nofRingsInEvent << endl;
226 for (Int_t iR = 0; iR < nofRingsInEvent; iR++) {
234 + TMath::Power(ringL.
GetCenterY() - trackY, 2));
235 fHM->
H1(
"fHCenterDistance")->Fill(DistCenters);
241 for (Int_t iH = 0; iH < NofHits; iH++) {
243 Int_t HitIndex = ring->
GetHit(iH);
245 Float_t xHit = hit->
GetX();
246 Float_t yHit = hit->
GetY();
247 Angles_0 = TMath::ATan2(
252 if (xRing - xHit == 0 || yRing - yHit == 0)
continue;
253 fPhi[iH] = TMath::ATan2(yHit - yRing, xHit - xRing);
257 Theta_Ch =
sqrt(TMath::Power(trackX - hit->
GetX(), 2)
258 + TMath::Power(trackY - hit->
GetY(), 2));
262 fHM->
H1(
"fHThetaDiff")->Fill(Theta_Ch - Theta_0);
265 fHM->
H2(
"fHCherenkovHitsDistribTheta0")->Fill(Angles_0, Theta_0);
266 fHM->
H2(
"fHCherenkovHitsDistribThetaCh")->Fill(
fPhi[iH], Theta_Ch);
267 fHM->
H2(
"fHCherenkovHitsDistribReduced")
268 ->Fill(
fPhi[iH], (Theta_Ch - Theta_0));
278 for (Int_t iP = 0; iP < NofProjections; iP++) {
283 cout <<
"Error: CbmRichAlignment::GetTrackPosition. No fair track param "
333 c3->SetFillColor(42);
335 gPad->SetTopMargin(0.1);
336 gPad->SetFillColor(33);
339 TH2D* CloneArr = (TH2D*)
fHM->
H2(
"fHCherenkovHitsDistribReduced")->Clone();
340 CloneArr->GetXaxis()->SetLabelSize(0.03);
341 CloneArr->GetXaxis()->SetTitleSize(0.03);
342 CloneArr->GetXaxis()->CenterTitle();
343 CloneArr->GetXaxis()->SetNdivisions(612, kTRUE);
344 CloneArr->GetYaxis()->SetLabelSize(0.03);
345 CloneArr->GetYaxis()->SetTitleSize(0.03);
346 CloneArr->GetYaxis()->SetNdivisions(612, kTRUE);
347 CloneArr->GetYaxis()->CenterTitle();
349 CloneArr->GetZaxis()->SetLabelSize(0.03);
350 CloneArr->GetZaxis()->SetTitleSize(0.03);
351 CloneArr->GetYaxis()->SetTitleOffset(1.0);
352 CloneArr->Draw(
"colz");
359 TH2D* CloneArr_2 = (TH2D*)
fHM->
H2(
"fHCherenkovHitsDistribReduced")->Clone();
360 for (Int_t y_bin = 1; y_bin <= 500; y_bin++) {
361 for (Int_t x_bin = 1; x_bin <= 200; x_bin++) {
367 if (CloneArr_2->GetBinContent(x_bin, y_bin) < thresh) {
368 CloneArr_2->SetBinContent(x_bin, y_bin, 0);
373 CloneArr_2->GetXaxis()->SetLabelSize(0.03);
374 CloneArr_2->GetXaxis()->SetTitleSize(0.03);
375 CloneArr_2->GetXaxis()->CenterTitle();
376 CloneArr_2->GetXaxis()->SetNdivisions(612, kTRUE);
377 CloneArr_2->GetYaxis()->SetLabelSize(0.03);
378 CloneArr_2->GetYaxis()->SetTitleSize(0.03);
379 CloneArr_2->GetYaxis()->SetNdivisions(612, kTRUE);
380 CloneArr_2->GetYaxis()->CenterTitle();
382 CloneArr_2->GetZaxis()->SetLabelSize(0.03);
383 CloneArr_2->GetZaxis()->SetTitleSize(0.03);
384 CloneArr_2->GetYaxis()->SetTitleOffset(1.0);
385 CloneArr_2->Draw(
"colz");
389 CloneArr_2->FitSlicesY(0, 0, -1, 1);
391 TH1D* histo_0 = (TH1D*) gDirectory->Get(
"fHCherenkovHitsDistribReduced_0");
394 TH1D* histo_1 = (TH1D*) gDirectory->Get(
"fHCherenkovHitsDistribReduced_1");
398 TH1D* histo_2 = (TH1D*) gDirectory->Get(
"fHCherenkovHitsDistribReduced_2");
402 (TH1D*) gDirectory->Get(
"fHCherenkovHitsDistribReduced_chi2");
406 TF1* f1 =
new TF1(
"f1",
"[2]+[0]*cos(x)+[1]*sin(x)", -3.5, 3.5);
407 f1->SetParameters(0, 0, 0);
408 f1->SetParNames(
"Delta_phi",
"Delta_lambda",
"Offset");
409 histo_1->Fit(
"f1",
"",
"");
410 TF1* fit = histo_1->GetFunction(
"f1");
411 Double_t p1 = fit->GetParameter(
"Delta_phi");
412 Double_t p2 = fit->GetParameter(
"Delta_lambda");
413 Double_t p3 = fit->GetParameter(
"Offset");
414 Double_t chi2 = fit->GetChisquare();
425 f1->SetParameters(fit->GetParameter(0), fit->GetParameter(1));
432 Double_t Focal_length = 150., q = 0., A = 0., Alpha = 0., mis_x = 0.,
436 q = TMath::ATan(fit->GetParameter(0) / fit->GetParameter(1));
438 <<
"fit_1 = " << fit->GetParameter(0)
439 <<
" and fit_2 = " << fit->GetParameter(1) << endl;
441 A = fit->GetParameter(1) / TMath::Cos(q);
444 TMath::ATan(A / 1.5) * 0.5
449 mis_x = TMath::ATan(fit->GetParameter(0) / Focal_length) * 0.5
450 * TMath::Power(10, 3);
451 mis_y = TMath::ATan(fit->GetParameter(1) / Focal_length) * 0.5
452 * TMath::Power(10, 3);
455 TLegend* LEG =
new TLegend(0.27, 0.7, 0.85, 0.87);
456 LEG->SetBorderSize(1);
457 LEG->SetFillColor(0);
459 LEG->SetTextSize(0.03);
460 sprintf(leg,
"Fitted sinusoid");
461 LEG->AddEntry(f1, leg,
"l");
462 sprintf(leg,
"Rotation angle around X = %f", mis_y);
463 LEG->AddEntry(
"", leg,
"l");
464 sprintf(leg,
"Rotation angle around Y = %f", mis_x);
465 LEG->AddEntry(
"", leg,
"l");
466 sprintf(leg,
"Offset = %f", fit->GetParameter(2));
467 LEG->AddEntry(
"", leg,
"l");
476 CloneArr_2->Draw(
"colz");
478 TLegend* LEG1 =
new TLegend(0.35, 0.7, 0.72, 0.85);
479 LEG1->SetBorderSize(1);
480 LEG1->SetFillColor(0);
481 LEG1->SetMargin(0.2);
482 LEG1->SetTextSize(0.03);
483 sprintf(leg,
"Fitted sinusoid");
484 LEG1->AddEntry(f1, leg,
"l");
485 sprintf(leg,
"Misalign in X = %f", mis_x);
486 LEG1->AddEntry(
"", leg,
"l");
487 sprintf(leg,
"Misalign in Y = %f", mis_y);
488 LEG1->AddEntry(
"", leg,
"l");
489 sprintf(leg,
"Offset = %f", fit->GetParameter(2));
490 LEG1->AddEntry(
"", leg,
"l");
531 outputFit.at(0) = mis_y;
532 outputFit.at(1) = mis_x;
533 outputFit.at(2) = fit->GetParameter(1);
534 outputFit.at(3) = fit->GetParameter(0);
539 TFile* file =
new TFile(fileName,
"READ");
548 "-----------------------------------------------------------------------"
549 "----------------------------------------------------------------- //"
552 <<
"// -------------------------------------------------- CbmRichAlignment "
553 "- Finish Function -------------------------------------------------- //"
557 "-----------------------------------------------------------------------"
558 "----------------------------------------------------------------- //"
565 vector<Double_t> outputFit(4);
567 cout << setprecision(6) << endl;
568 cout <<
"Horizontal displacement = " << outputFit[0]
569 <<
" [mrad] and vertical displacement = " << outputFit[1] <<
" [mrad]."
573 "fHCherenkovHitsDistribReduced",
574 "fHCherenkovHitsDistribReduced;Phi_Ch [rad];Th_Ch-Th_0 [cm];Entries",
592 if (corr_file.is_open()) {
593 corr_file << setprecision(7) << outputFit[0] <<
"\t";
594 corr_file << setprecision(7) << outputFit[1] <<
"\t";
595 corr_file << setprecision(7) << outputFit[2] <<
"\t";
596 corr_file << setprecision(7) << outputFit[3] <<
"\t";
598 cout <<
"Wrote correction parameters to: " << s << endl;
601 <<
"Error in CbmRichAlignment::Finish ; unable to open parameter file!"