CbmRoot
CbmRichAlignment.cxx
Go to the documentation of this file.
1 // ---------- Original Headers ---------- //
2 #include "CbmRichAlignment.h"
3 #include "FairLogger.h"
4 #include "FairRootManager.h"
5 
6 #include "TClonesArray.h"
7 
8 #include "CbmDrawHist.h"
9 #include "CbmRichHit.h"
10 #include "TCanvas.h"
11 #include "TF1.h"
12 #include "TH1D.h"
13 #include "TH2D.h"
14 
15 #include <iostream>
16 #include <vector>
17 
18 // ---------- Included Headers ---------- //
19 #include "CbmMCTrack.h"
20 #include "CbmRichConverter.h"
21 #include "CbmRichPoint.h"
22 #include "CbmRichRingFitterCOP.h"
24 #include "CbmRichRingLight.h"
25 #include "CbmTrackMatchNew.h"
26 #include "CbmUtils.h"
27 #include "FairTrackParam.h"
28 #include "FairVolume.h"
29 #include "TEllipse.h"
30 #include "TGeoManager.h"
31 #include <algorithm>
32 #include <fstream>
33 #include <iomanip>
34 #include <stdlib.h>
35 //#include <stdio.h>
36 #include "CbmGlobalTrack.h"
37 #include "CbmRichGeoManager.h"
38 #include "CbmRichHitProducer.h"
39 
40 //#include "TLorentzVector.h"
41 #include "TGeoSphere.h"
42 #include "TVirtualMC.h"
43 class TGeoNode;
44 class TGeoVolume;
45 class TGeoShape;
46 class TGeoMatrix;
47 #include <sstream>
48 
50  : FairTask()
51  , fRichHits(NULL)
52  , fRichRings(NULL)
53  , fRichProjections(NULL)
54  , fRichMirrorPoints(NULL)
55  , fMCTracks(NULL)
56  , fRichRingMatches(NULL)
57  , fRichPoints(NULL)
58  , fHM(NULL)
59  ,
60  //fGP(),
61  fEventNum(0)
62  , fOutputDir("")
63  , fRunTitle("")
64  , fAxisRotTitle("")
65  , fNumbAxis(0)
66  , fTile(0)
67  , fDrawAlignment(kTRUE)
68  , fCopFit(NULL)
69  , fTauFit(NULL)
70  , fPhi() {}
71 
73 
74 InitStatus CbmRichAlignment::Init() {
75  FairRootManager* manager = FairRootManager::Instance();
76 
77  fRichHits = (TClonesArray*) manager->GetObject("RichHit");
78  if (NULL == fRichHits) {
79  Fatal("CbmRichAlignment::Init", "No RichHit array !");
80  }
81 
82  fRichRings = (TClonesArray*) manager->GetObject("RichRing");
83  if (NULL == fRichRings) {
84  Fatal("CbmRichAlignment::Init", "No RichRing array !");
85  }
86 
87  fRichProjections = (TClonesArray*) manager->GetObject("RichProjection");
88  if (NULL == fRichProjections) {
89  Fatal("CbmRichAlignment::Init", "No RichProjection array !");
90  }
91 
92  fRichPoints = (TClonesArray*) manager->GetObject("RichPoint");
93  if (NULL == fRichPoints) {
94  Fatal("CbmRichAlignment::Init", "No RichPoint array !");
95  }
96 
97  fMCTracks = (TClonesArray*) manager->GetObject("MCTrack");
98  if (NULL == fMCTracks) {
99  Fatal("CbmRichAlignment::Init", "No MCTracks array !");
100  }
101 
102  fRichRingMatches = (TClonesArray*) manager->GetObject("RichRingMatch");
103  if (NULL == fRichRingMatches) {
104  Fatal("CbmRichAlignment::Init", "No RichRingMatches array !");
105  }
106 
107  fRichMirrorPoints = (TClonesArray*) manager->GetObject("RichMirrorPoint");
108  if (NULL == fRichMirrorPoints) {
109  Fatal("CbmRichAlignment::Init", "No RichMirrorPoints array !");
110  }
111 
112  /* fRichRefPlanePoints = (TClonesArray*) manager->GetObject("RefPlanePoint");
113  if (NULL == fRichRefPlanePoints) { Fatal("CbmRichAlignment::Init", "No RichRefPlanePoint array !"); }
114 
115  fRichMCPoints = (TClonesArray*) manager->GetObject("RichPoint");
116  if (NULL == fRichMCPoints) { Fatal("CbmRichAlignment::Init", "No RichMCPoints array !"); }
117 
118  fGlobalTracks = (TClonesArray*) manager->GetObject("GlobalTrack");
119  if (NULL == fGlobalTracks) { Fatal("CbmAnaDielectronTask::Init","No GlobalTrack array!"); }
120 */
121 
125 
127 
128  return kSUCCESS;
129 }
130 
132  fHM = new CbmHistManager();
133 
134  fHM->Create1<TH1D>("fHCenterDistance",
135  "fHCenterDistance;Distance C-C';Nb of entries",
136  100,
137  -0.1,
138  5.);
139  fHM->Create1<TH1D>(
140  "fHPhi", "fHPhi;Phi_Ch [rad];Nb of entries", 200, -3.4, 3.4);
141  fHM->Create1<TH1D>(
142  "fHThetaDiff", "fHThetaDiff;Th_Ch-Th_0 [cm];Nb of entries", 252, -5., 5.);
143  fHM->Create2<TH2D>(
144  "fHCherenkovHitsDistribTheta0",
145  "fHCherenkovHitsDistribTheta0;Phi_0 [rad];Theta_0 [cm];Entries",
146  200,
147  -2.,
148  2.,
149  600,
150  2.,
151  8.);
152  fHM->Create2<TH2D>(
153  "fHCherenkovHitsDistribThetaCh",
154  "fHCherenkovHitsDistribThetaCh;Phi_Ch [rad];Theta_Ch [cm];Entries",
155  200,
156  -3.4,
157  3.4,
158  600,
159  0.,
160  20);
161  fHM->Create2<TH2D>(
162  "fHCherenkovHitsDistribReduced",
163  "fHCherenkovHitsDistribReduced;Phi_Ch [rad];Th_Ch-Th_0 [cm];Entries",
164  200,
165  -3.4,
166  3.4,
167  500,
168  -5.,
169  5.);
170 }
171 
172 void CbmRichAlignment::Exec(Option_t* option) {
173  cout << endl
174  << "//"
175  "--------------------------------------------------------------------"
176  "-----------------------------------------------//"
177  << endl;
178  cout << "//------------------------------------------ CbmRichAlignment: EXEC "
179  "Function ------------------------------------------//"
180  << endl;
181  cout << "//"
182  "--------------------------------------------------------------------"
183  "--------------------------------------------------//"
184  << endl;
185  fEventNum++;
186  //LOG(debug2) << "CbmRichAlignment : Event #" << fEventNum;
187  cout << "CbmRichAlignment : Event #" << fEventNum << endl;
188 
189  Int_t nofRingsInEvent = fRichRings->GetEntries();
190  Int_t nofMirrorPoints = fRichMirrorPoints->GetEntries();
191  Int_t nofHitsInEvent = fRichHits->GetEntries();
192  Int_t NofMCTracks = fMCTracks->GetEntriesFast();
193  // Int_t NofMCPoints = fRichMCPoints->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
198  << endl; //", nb of Monte-Carlo points = " << NofMCPoints <<
199 
200  if (nofRingsInEvent == 0) {
201  cout << "Error no rings registered in event." << endl << endl;
202  } else {
204  }
205 }
206 
208  cout << "//------------------------------ CbmRichAlignment: Calculate Angles "
209  "& Draw Distrib ------------------------------//"
210  << endl
211  << endl;
212 
213  Double_t trackX = 0., trackY = 0.;
214  GetTrackPosition(trackX, trackY);
215 
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;
220  fPhi.resize(kMAX_NOF_HITS);
221 
222  // ------------------------- Loop to get ring center coordinates and photon hit coordinates per ring and per event ------------------------- //
223  if (nofRingsInEvent >= 1) {
224  cout << "Number of Rings in event: " << nofRingsInEvent << endl;
225  //sleep(2);
226  for (Int_t iR = 0; iR < nofRingsInEvent; iR++) {
227  // ----- Convert Ring to Ring Light ----- //
228  CbmRichRing* ring = static_cast<CbmRichRing*>(fRichRings->At(iR));
229  CbmRichRingLight ringL;
231  fCopFit->DoFit(&ringL);
232  // ----- Distance between mean center and fitted center calculation ----- //
233  DistCenters = sqrt(TMath::Power(ringL.GetCenterX() - trackX, 2)
234  + TMath::Power(ringL.GetCenterY() - trackY, 2));
235  fHM->H1("fHCenterDistance")->Fill(DistCenters);
236  // ----- Declaration of new variables ----- //
237  Int_t NofHits = ringL.GetNofHits();
238  Float_t xRing = ringL.GetCenterX();
239  Float_t yRing = ringL.GetCenterY();
240 
241  for (Int_t iH = 0; iH < NofHits; iH++) {
242  // ----- Phi angle calculation ----- //
243  Int_t HitIndex = ring->GetHit(iH);
244  CbmRichHit* hit = static_cast<CbmRichHit*>(fRichHits->At(HitIndex));
245  Float_t xHit = hit->GetX();
246  Float_t yHit = hit->GetY();
247  Angles_0 = TMath::ATan2(
248  (hit->GetX() - ringL.GetCenterX()),
249  (ringL.GetCenterY() - hit->GetY())); //* TMath::RadToDeg();
250  //cout << "Angles_0 = " << Angles_0[iH] << endl;
251 
252  if (xRing - xHit == 0 || yRing - yHit == 0) continue;
253  fPhi[iH] = TMath::ATan2(yHit - yRing, xHit - xRing);
254  fHM->H1("fHPhi")->Fill(fPhi[iH]);
255 
256  // ----- Theta_Ch and Theta_0 calculations ----- //
257  Theta_Ch = sqrt(TMath::Power(trackX - hit->GetX(), 2)
258  + TMath::Power(trackY - hit->GetY(), 2));
259  Theta_0 = sqrt(TMath::Power(ringL.GetCenterX() - hit->GetX(), 2)
260  + TMath::Power(ringL.GetCenterY() - hit->GetY(), 2));
261  //cout << "Theta_0 = " << Theta_0 << endl;
262  fHM->H1("fHThetaDiff")->Fill(Theta_Ch - Theta_0);
263 
264  // ----- Filling of final histograms ----- //
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));
269  }
270  //cout << endl;
271  }
272  }
273 }
274 
275 void CbmRichAlignment::GetTrackPosition(Double_t& x, Double_t& y) {
276  Int_t NofProjections = fRichProjections->GetEntries();
277  //cout << "!!! NB PROJECTIONS !!! " << NofProjections << endl;
278  for (Int_t iP = 0; iP < NofProjections; iP++) {
279  FairTrackParam* pr = (FairTrackParam*) fRichProjections->At(iP);
280  if (NULL == pr) {
281  x = 0.;
282  y = 0.;
283  cout << "Error: CbmRichAlignment::GetTrackPosition. No fair track param "
284  "found."
285  << endl;
286  }
287  x = pr->GetX();
288  y = pr->GetY();
289  //cout << "Center X: " << *x << " and Center y: " << *y << endl;
290  }
291 }
292 
294  TCanvas* c1 = new TCanvas(fRunTitle + "_Data_Histograms_" + fAxisRotTitle,
295  fRunTitle + "_Data_Histograms_" + fAxisRotTitle,
296  800,
297  400);
298  c1->Divide(2, 1);
299  c1->cd(1);
300  /*c1->SetGridx(1);
301  c1->SetGridy(1);
302  fHM2->H1("fHCenterDistance")->GetXaxis()->SetLabelSize(0.03);
303  fHM2->H1("fHCenterDistance")->GetXaxis()->SetTitleSize(0.03);
304  fHM2->H1("fHCenterDistance")->GetXaxis()->CenterTitle();
305  fHM2->H1("fHCenterDistance")->GetXaxis()->SetNdivisions(612,kTRUE);
306  fHM2->H1("fHCenterDistance")->GetYaxis()->CenterTitle();
307  fHM2->H1("fHCenterDistance")->GetYaxis()->SetTitleSize(0.03);
308  fHM2->H1("fHCenterDistance")->GetYaxis()->SetTitleOffset(1.);
309  fHM2->H1("fHCenterDistance")->Draw();*/
310  DrawH1(fHM->H1("fHCenterDistance"));
311  c1->cd(2);
312  DrawH2(fHM->H2("fHCherenkovHitsDistribReduced"));
313  /* c1->cd(3);
314  DrawH1(fHM->H1("fHThetaDiff"));
315  c1->cd(4);
316  DrawH2(fHM->H2("fHCherenkovHitsDistribTheta0"));
317  c1->cd(5);
318  DrawH1(fHM->H1("fHPhi"));
319  c1->cd(6);
320  DrawH2(fHM->H2("fHCherenkovHitsDistribThetaCh"));
321 */
322  Cbm::SaveCanvasAsImage(c1, string(fOutputDir.Data()), "png");
323 }
324 
325 void CbmRichAlignment::DrawFit(vector<Double_t>& outputFit, Int_t thresh) {
326  //vector<Double_t> paramVect;
327  //paramVect.reserve(5);
328 
329  TCanvas* c3 = new TCanvas(fRunTitle + "_Fit_Histograms_" + fAxisRotTitle,
330  fRunTitle + "_Fit_Histograms_" + fAxisRotTitle,
331  1100,
332  600);
333  c3->SetFillColor(42);
334  c3->Divide(4, 2);
335  gPad->SetTopMargin(0.1);
336  gPad->SetFillColor(33);
337  c3->cd(1);
338  // TH2D* CloneArr = (TH2D*)fHM2->H2("fHCherenkovHitsDistribThetaCh")->Clone();
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();
348  //CloneArr->GetYaxis()->SetRangeUser(-2.5,2.5);
349  CloneArr->GetZaxis()->SetLabelSize(0.03);
350  CloneArr->GetZaxis()->SetTitleSize(0.03);
351  CloneArr->GetYaxis()->SetTitleOffset(1.0);
352  CloneArr->Draw("colz");
353  //Double_t ymax = CloneArr->GetYaxis()->GetXmax();
354  //Double_t ymin = CloneArr->GetYaxis()->GetXmin();
355  //TF1 *fgauss = TF1 fgauss("gauss", "[0]*exp(-0.5*((x-[1])/[2])**2)", 0, 100);
356 
357  // ------------------------------ APPLY THRESHOLD TO 2D-HISTO ------------------------------ //
358  // TH2D* CloneArr_2 = (TH2D*)fHM2->H2("fHCherenkovHitsDistribThetaCh")->Clone();
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++) {
362  /*if (CloneArr_2->GetBinContent(x_bin, y_bin)!=0) {
363  cout << "Bin Content: " << CloneArr_2->GetBinContent(x_bin, y_bin) << endl;
364  sleep(1);
365  }
366  else;*/
367  if (CloneArr_2->GetBinContent(x_bin, y_bin) < thresh) {
368  CloneArr_2->SetBinContent(x_bin, y_bin, 0);
369  }
370  }
371  }
372  c3->cd(2);
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();
381  //CloneArr_2->GetYaxis()->SetRangeUser(-2.5,2.5);
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");
386  CloneArr_2->Write();
387 
388  // -------------------- FIT SLICES AND FIT THE MEAN OF THE RESULT TO A SIN FUNCTION -------------------- //
389  CloneArr_2->FitSlicesY(0, 0, -1, 1);
390  c3->cd(3);
391  TH1D* histo_0 = (TH1D*) gDirectory->Get("fHCherenkovHitsDistribReduced_0");
392  histo_0->Draw();
393  c3->cd(4);
394  TH1D* histo_1 = (TH1D*) gDirectory->Get("fHCherenkovHitsDistribReduced_1");
395  //histo_1->GetYaxis()->SetRangeUser(-2.5, 2.5);
396  histo_1->Draw();
397  c3->cd(5);
398  TH1D* histo_2 = (TH1D*) gDirectory->Get("fHCherenkovHitsDistribReduced_2");
399  histo_2->Draw();
400  c3->cd(6);
401  TH1D* histo_chi2 =
402  (TH1D*) gDirectory->Get("fHCherenkovHitsDistribReduced_chi2");
403  histo_chi2->Draw();
404 
405  c3->cd(7);
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();
415  //cout << setprecision(6) << "Delta_phi = " << fit->GetParameter(0) << " and delta_lambda = " << fit->GetParameter(1) << endl;
416  //cout << "Delta_phi error = " << fit->GetParError(0) << " and delta_lambda error = " << fit->GetParError(1) << endl;
417  //cout << endl << "Chi2: " << chi2 << endl;
418 
419  /* paramVect.push_back(fit->GetParameter("Delta_phi"));
420  paramVect.push_back(fit->GetParameter("Delta_lambda"));
421  paramVect.push_back(fit->GetParameter("Offset"));
422  paramVect.push_back(fit->GetChisquare());
423  //cout << "Vectors: Delta_phi = " << paramVect[0] << ", Delta_lambda = " << paramVect[1] << ", Offset = " << paramVect[2] << endl;
424 */
425  f1->SetParameters(fit->GetParameter(0), fit->GetParameter(1));
426  char leg[128];
427  f1->SetLineColor(2);
428  f1->Draw();
429  f1->Write();
430 
431  // ------------------------------ CALCULATION OF MISALIGNMENT ANGLE ------------------------------ //
432  Double_t Focal_length = 150., q = 0., A = 0., Alpha = 0., mis_x = 0.,
433  mis_y = 0.;
434  // mis_x && mis_y corresponds respect. to rotation angles around the Y and X axes.
435  // !!! BEWARE: AXES INDEXES ARE SWITCHED !!!
436  q = TMath::ATan(fit->GetParameter(0) / fit->GetParameter(1));
437  cout << endl
438  << "fit_1 = " << fit->GetParameter(0)
439  << " and fit_2 = " << fit->GetParameter(1) << endl;
440  //cout << "q = " << q << endl;
441  A = fit->GetParameter(1) / TMath::Cos(q);
442  //cout << "Parameter a = " << A << endl;
443  Alpha =
444  TMath::ATan(A / 1.5) * 0.5
445  * TMath::Power(
446  10,
447  3); // *0.5, because a mirror rotation of alpha implies a rotation in the particle trajectory of 2*alpha ; 1.5 meters = Focal length = Radius_of_curvature/2
448  //cout << setprecision(6) << "Total angle of misalignment alpha = " << Alpha << endl; // setprecision(#) gives the number of digits in the cout.
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);
453  //cout << "Horizontal displacement = " << outputFit[0] << " [mrad] and vertical displacement = " << outputFit[1] << " [mrad]." << endl;
454 
455  TLegend* LEG = new TLegend(0.27, 0.7, 0.85, 0.87); // Set legend position
456  LEG->SetBorderSize(1);
457  LEG->SetFillColor(0);
458  LEG->SetMargin(0.2);
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");
468  LEG->Draw();
469  //Cbm::SaveCanvasAsImage(c3, string(fOutputDir.Data()), "png");
470 
471  TCanvas* c4 = new TCanvas(fRunTitle + "_Sinus_Fit_" + fAxisRotTitle,
472  fRunTitle + "_Sinus_Fit_" + fAxisRotTitle,
473  400,
474  400);
475  c4->SetGrid(1, 1);
476  CloneArr_2->Draw("colz");
477  f1->Draw("same");
478  TLegend* LEG1 = new TLegend(0.35, 0.7, 0.72, 0.85); // Set legend position
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");
491  LEG1->Draw();
492  Cbm::SaveCanvasAsImage(c4, string(fOutputDir.Data()), "png");
493 
494  // ------------------------------ APPLY SECOND FIT USING LOG-LIKELIHOOD METHOD ------------------------------ //
495  /* TCanvas* c4 = new TCanvas(fRunTitle + "_Second_Fit_" + fAxisRotTitle, fRunTitle + "_Second_Fit_" + fAxisRotTitle, 600, 600);
496  c4->SetFillColor(42);
497  gPad->SetTopMargin(0.1);
498  gPad->SetFillColor(33);
499  f1->SetParameters(fit->GetParameter(0), fit->GetParameter(1), fit->GetParameter(2));
500  histo_1->Fit("f1","L","");
501  TF1 *fit2 = histo_1->GetFunction("f1");
502  f1->SetParameters(fit2->GetParameter(0), fit2->GetParameter(1), fit2->GetParameter(2));
503  f1->Draw();
504 
505  Double_t q_2 = TMath::ATan(fit2->GetParameter(0)/fit2->GetParameter(1));
506  //cout << endl << "fit2_1 = " << fit2->GetParameter(0) << " and fit2_2 = " << fit2->GetParameter(1) << endl;
507  //cout << "q_2 = " << q_2 << endl;
508  Double_t A_2 = fit2->GetParameter(1)/TMath::Cos(q_2);
509  //cout << "Parameter a_2 = " << A_2 << endl;
510  Double_t Alpha_2 = TMath::ATan(A_2/1.5)*0.5*TMath::Power(10,3);
511  //cout << setprecision(6) << "Total angle of misalignment alpha_2 = " << Alpha_2 << endl;
512  Double_t mis_x_2 = TMath::ATan(fit2->GetParameter(0)/Focal_length)*0.5*TMath::Power(10,3);
513  Double_t mis_y_2 = TMath::ATan(fit2->GetParameter(1)/Focal_length)*0.5*TMath::Power(10,3);
514 
515  TLegend* LEG2= new TLegend(0.31,0.7,0.72,0.85); // Set legend position
516  LEG2->SetBorderSize(1);
517  LEG2->SetFillColor(0);
518  LEG2->SetMargin(0.2);
519  LEG2->SetTextSize(0.03);
520  sprintf(leg, "Fitted sinusoid");
521  LEG2->AddEntry(f1, leg, "l");
522  sprintf(leg, "Misalign in X = %f", mis_x_2);
523  LEG2->AddEntry("", leg, "l");
524  sprintf(leg, "Misalign in Y = %f", mis_y_2);
525  LEG2->AddEntry("", leg, "l");
526  sprintf(leg, "Offset = %f", fit2->GetParameter(2));
527  LEG2->AddEntry("", leg, "l");
528  LEG2->Draw();
529  Cbm::SaveCanvasAsImage(c4, string(fOutputDir.Data()), "png");*/
530 
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);
535 }
536 
537 void CbmRichAlignment::DrawHistFromFile(TString fileName) {
538  fHM = new CbmHistManager();
539  TFile* file = new TFile(fileName, "READ");
540  fHM->ReadFromFile(file);
542 }
543 
545  cout
546  << endl
547  << "// "
548  "-----------------------------------------------------------------------"
549  "----------------------------------------------------------------- //"
550  << endl;
551  cout
552  << "// -------------------------------------------------- CbmRichAlignment "
553  "- Finish Function -------------------------------------------------- //"
554  << endl;
555  cout
556  << "// "
557  "-----------------------------------------------------------------------"
558  "----------------------------------------------------------------- //"
559  << endl
560  << endl;
561 
562  if (fDrawAlignment) {
564  Int_t thresh = 5;
565  vector<Double_t> outputFit(4);
566  DrawFit(outputFit, thresh);
567  cout << setprecision(6) << endl;
568  cout << "Horizontal displacement = " << outputFit[0]
569  << " [mrad] and vertical displacement = " << outputFit[1] << " [mrad]."
570  << endl;
571 
572  fHM->Create2<TH2D>(
573  "fHCherenkovHitsDistribReduced",
574  "fHCherenkovHitsDistribReduced;Phi_Ch [rad];Th_Ch-Th_0 [cm];Entries",
575  200,
576  -3.4,
577  3.4,
578  500,
579  -5.,
580  5.);
581 
582  ofstream corr_file;
583  /* // Converting double to string.
584  TString s;
585  std::ostringstream strs;
586  strs << fNumb;
587  std::string str = strs.str();
588  s = "/data/misalignment_correction/Sim_Outputs/Alignment_Correction/correction_param" + str + ".txt";
589 */
590  TString s = fOutputDir + "correction_param_" + fNumbAxis + fTile + ".txt";
591  corr_file.open(s);
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";
597  corr_file.close();
598  cout << "Wrote correction parameters to: " << s << endl;
599  } else {
600  cout
601  << "Error in CbmRichAlignment::Finish ; unable to open parameter file!"
602  << endl;
603  }
604  }
605 
606  //cout << setprecision(6) << endl;
607 }
CbmRichAlignment::fRichProjections
TClonesArray * fRichProjections
Definition: CbmRichAlignment.h:113
CbmRichPoint.h
CbmRichAlignment::DrawHistAlignment
void DrawHistAlignment()
Definition: CbmRichAlignment.cxx:293
CbmRichAlignment
Definition: CbmRichAlignment.h:18
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmRichRingFitterEllipseTau
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
Definition: CbmRichRingFitterEllipseTau.h:35
CbmRichAlignment::CalculateAnglesAndDrawDistrib
void CalculateAnglesAndDrawDistrib()
Definition: CbmRichAlignment.cxx:207
CbmRichAlignment::fEventNum
UInt_t fEventNum
Definition: CbmRichAlignment.h:123
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmRichRingFitterEllipseTau.h
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmRichAlignment::fNumbAxis
TString fNumbAxis
Definition: CbmRichAlignment.h:124
CbmHistManager::ReadFromFile
void ReadFromFile(TFile *file)
Read histograms from file.
Definition: core/base/CbmHistManager.cxx:110
CbmHistManager::Create2
void Create2(const std::string &name, const std::string &title, Int_t nofBinsX, Double_t minBinX, Double_t maxBinX, Int_t nofBinsY, Double_t minBinY, Double_t maxBinY)
Helper function for creation of 2-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:104
CbmHistManager::H2
TH2 * H2(const std::string &name) const
Return pointer to TH2 histogram.
Definition: CbmHistManager.h:190
CbmGlobalTrack.h
CbmRichAlignment::~CbmRichAlignment
virtual ~CbmRichAlignment()
Definition: CbmRichAlignment.cxx:72
CbmRichAlignment::fRichPoints
TClonesArray * fRichPoints
Definition: CbmRichAlignment.h:114
CbmRichRing
Definition: CbmRichRing.h:17
CbmRichAlignment::CbmRichAlignment
CbmRichAlignment()
Definition: CbmRichAlignment.cxx:49
CbmDrawHist.h
Helper functions for drawing 1D and 2D histograms and graphs.
CbmRichAlignment::GetTrackPosition
void GetTrackPosition(Double_t &x, Double_t &y)
Definition: CbmRichAlignment.cxx:275
CbmRichAlignment::fPhi
vector< Float_t > fPhi
Definition: CbmRichAlignment.h:127
CbmRichAlignment::fRichRingMatches
TClonesArray * fRichRingMatches
Definition: CbmRichAlignment.h:116
CbmRichRingFitterCOP
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
Definition: CbmRichRingFitterCOP.h:28
CbmRichGeoManager.h
DrawH1
void DrawH1(TH1 *hist, HistScale logx, HistScale logy, const string &drawOpt, Int_t color, Int_t lineWidth, Int_t lineStyle, Int_t markerSize, Int_t markerStyle)
Definition: CbmDrawHist.cxx:49
CbmRichRingLight::GetNofHits
int GetNofHits() const
Return number of hits in ring.
Definition: CbmRichRingLight.h:108
CbmRichRing::GetHit
UInt_t GetHit(Int_t i) const
Definition: CbmRichRing.h:42
CbmRichAlignment::fOutputDir
TString fOutputDir
Definition: CbmRichAlignment.h:129
CbmRichAlignment::fHM
CbmHistManager * fHM
Definition: CbmRichAlignment.h:121
CbmHistManager
Histogram manager.
Definition: CbmHistManager.h:41
CbmRichAlignment::DrawHistFromFile
void DrawHistFromFile(TString fileName)
Definition: CbmRichAlignment.cxx:537
CbmRichConverter.h
Convert internal data classes to cbmroot common data classes.
CbmRichRingLight::GetCenterY
float GetCenterY() const
Definition: CbmRichRingLight.h:160
CbmRichRingLight.h
CbmRichAlignment.h
CbmRichAlignment::fRichHits
TClonesArray * fRichHits
Definition: CbmRichAlignment.h:111
CbmTrackMatchNew.h
CbmRichAlignment::fMCTracks
TClonesArray * fMCTracks
Definition: CbmRichAlignment.h:115
CbmRichAlignment::InitHistAlignment
void InitHistAlignment()
Definition: CbmRichAlignment.cxx:131
CbmHistManager::Create1
void Create1(const std::string &name, const std::string &title, Int_t nofBins, Double_t minBin, Double_t maxBin)
Helper function for creation of 1-dimensional histograms and profiles. Template argument is a real ob...
Definition: CbmHistManager.h:81
CbmRichAlignment::fTile
TString fTile
Definition: CbmRichAlignment.h:125
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmHistManager::H1
TH1 * H1(const std::string &name) const
Return pointer to TH1 histogram.
Definition: CbmHistManager.h:170
CbmRichAlignment::fTauFit
CbmRichRingFitterEllipseTau * fTauFit
Definition: CbmRichAlignment.h:134
CbmUtils.h
CbmRichRingFitterCOP.h
Here the ring is fitted with the COP algorithm from A. Ayriyan/G. Ososkov.
CbmRichConverter::Init
static void Init()
Initialize array of RICH hits.
Definition: CbmRichConverter.h:93
CbmRichAlignment::kMAX_NOF_HITS
static const int kMAX_NOF_HITS
Definition: CbmRichAlignment.h:20
CbmRichAlignment::DrawFit
void DrawFit(vector< Double_t > &outputFit, Int_t thresh)
Definition: CbmRichAlignment.cxx:325
CbmRichRingFitterCOP::DoFit
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition: CbmRichRingFitterCOP.cxx:19
CbmMCTrack.h
CbmRichAlignment::fDrawAlignment
Bool_t fDrawAlignment
Definition: CbmRichAlignment.h:126
CbmRichAlignment::Init
virtual InitStatus Init()
Inherited from FairTask.
Definition: CbmRichAlignment.cxx:74
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichHitProducer.h
Class for producing RICH hits directly from MCPoints.
CbmRichConverter::CopyHitsToRingLight
static void CopyHitsToRingLight(const CbmRichRing *ring1, CbmRichRingLight *ring2)
Copy hits from CbmRichRing to CbmRichRingLight.
Definition: CbmRichConverter.h:41
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichAlignment::Exec
virtual void Exec(Option_t *option)
Inherited from FairTask.
Definition: CbmRichAlignment.cxx:172
DrawH2
void DrawH2(TH2 *hist, HistScale logx, HistScale logy, HistScale logz, const string &drawOpt)
Definition: CbmDrawHist.cxx:84
Cbm::SaveCanvasAsImage
void SaveCanvasAsImage(TCanvas *c, const std::string &dir, const std::string &option)
Definition: CbmUtils.cxx:20
CbmRichAlignment::fRunTitle
TString fRunTitle
Definition: CbmRichAlignment.h:130
CbmRichAlignment::fAxisRotTitle
TString fAxisRotTitle
Definition: CbmRichAlignment.h:131
CbmRichHit.h
CbmRichRingLight::GetCenterX
float GetCenterX() const
Definition: CbmRichRingLight.h:159
CbmRichRingLight
Definition: CbmRichRingLight.h:39
CbmRichAlignment::Finish
virtual void Finish()
Inherited from FairTask.
Definition: CbmRichAlignment.cxx:544
CbmRichHit
Definition: CbmRichHit.h:19
CbmRichAlignment::fRichMirrorPoints
TClonesArray * fRichMirrorPoints
Definition: CbmRichAlignment.h:117
CbmRichAlignment::fCopFit
CbmRichRingFitterCOP * fCopFit
Definition: CbmRichAlignment.h:133
CbmRichAlignment::fRichRings
TClonesArray * fRichRings
Definition: CbmRichAlignment.h:112