EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CbmRichRingFitterQa.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CbmRichRingFitterQa.cxx
1 
8 #include "CbmRichRingFitterQa.h"
9 #include "CbmRichRingFitterEllipseTau.h"
10 #include "CbmRichRingFitterEllipseMinuit.h"
11 #include "CbmRichRingFitterCOP.h"
12 
13 #include "TRandom.h"
14 #include "TMath.h"
15 #include "CbmRichRing.h"
16 #include "TCanvas.h"
17 #include "TMath.h"
18 #include "TMatrixD.h"
19 #include "TH1D.h"
20 
21 #include <iostream>
22 #include <vector>
23 
24 using std::vector;
25 using std::cout;
26 using std::endl;
27 
28 CbmRichRingFitterQa::CbmRichRingFitterQa():
29  fhErrorA(NULL),
30  fhErrorB(NULL),
31  fhErrorX(NULL),
32  fhErrorY(NULL),
33  fhErrorPhi(NULL),
34 
35  fhA(NULL),
36  fhB(NULL),
37  fhX(NULL),
38  fhY(NULL),
39  fhPhi(NULL),
40 
41  fhRadiusErr(NULL),
42  fhCircleXcErr(NULL),
43  fhCircleYcErr(NULL),
44 
45  fhRadius(NULL),
46  fhCircleXc(NULL),
47  fhCircleYc(NULL),
48 
49  fhRadiusPool(NULL),
50  fhCircleXcPool(NULL),
51  fhCircleYcPool(NULL)
52 {
53  fhErrorA = new TH1D("fhErrorA","fhErrorA;dA [cm];Counter", 100, -2., 2.);
54  fhErrorB = new TH1D("fhErrorB","fhErrorB;B [cm];Counter", 100, -2., 2.);
55  fhErrorX = new TH1D("fhErrorX","fhErrorX;X [cm];Counter", 100, -2., 2.);
56  fhErrorY = new TH1D("fhErrorY","fhErrorY;Y [cm];Counter", 100, -2., 2.);
57  fhErrorPhi = new TH1D("fhErrorPhi","fhErrorPhi;d#Phi [rad];Counter", 100, -2., 2.);
58 
59  fhA = new TH1D("fhA","fhA;A [cm];Counter", 100, 5., 7.);
60  fhB = new TH1D("fhB","fhB;B [cm];Counter", 100, 5., 7.);
61  fhX = new TH1D("fhX","fhX;X [cm];Counter", 100, -1., 1.);
62  fhY = new TH1D("fhY","fhY;Y [cm];Counter", 100, -1., 1.);
63  Double_t pi = TMath::Pi();
64  fhPhi = new TH1D("fhPhi","fhPhi;#Phi [rad];Counter", 100, - pi/2. - pi/6. , pi/2. + pi/6.);
65 
66  // circle fitting
67  fhRadiusErr = new TH1D("fhRadiusErr","fhRadiusErr;dR [cm];Counter", 100, -2., 2.);
68  fhCircleXcErr = new TH1D("fhCircleXcErr","fhCircleXcErr;dXc [cm];Counter", 100, -2., 2.);
69  fhCircleYcErr = new TH1D("fhCircleYcErr","fhCircleYcErr;dYc [cm];Counter", 100, -2., 2.);
70  fhRadius = new TH1D("fhRadius","fhRadius;Radius [cm];Counter", 100, 4., 8.);
71  fhCircleXc = new TH1D("fhCircleXc","fhCircleXc;Xc [cm];Counter", 100, -2., 2.);
72  fhCircleYc = new TH1D("fhCircleYc","fhCircleYc;Yc [cm];Counter", 100, -2., 2.);
73  fhRadiusPool = new TH1D("fhRadiusPool","fhRadiusPool;Pool R;Counter", 100, -5., 5.);
74  fhCircleXcPool = new TH1D("fhCircleXcPool","fhCircleXcPool;Pool Xc;Counter", 100, -5., 5.);
75  fhCircleYcPool = new TH1D("fhCircleYcPool","fhCircleYcPool;Pool Yc;Counter", 100, -5., 5.);
76 }
77 
78 CbmRichRingFitterQa::~CbmRichRingFitterQa()
79 {
80 
81 }
82 
83 void CbmRichRingFitterQa::GenerateEllipse()
84 {
85  Double_t maxX = 15;
86  Double_t maxY = 15;
87  Int_t nofHits = 50;
88  Double_t A = 6.;
89  Double_t B = 6.;
90  Double_t sigmaError = 0.2;
91  CbmRichRingLight ellipse;
92  Int_t nofBadFit = 0;
93  Double_t X0 = 0.;//gRandom->Rndm()*(maxX - A);
94  Double_t Y0 = 0.;//gRandom->Rndm()* (maxY - A);
95 
96  CbmRichRingFitterEllipseTau * fitEllipse = new CbmRichRingFitterEllipseTau();
97  CbmRichRingFitterCOP * fitCircle = new CbmRichRingFitterCOP();
98 
99  for (Int_t iR = 0; iR < 50000; iR++){
100  Double_t phi = 0.;//TMath::Pi()*(6./12.); //gRandom->Rndm()*TMath::Pi() - TMath::Pi()/2.;
101  ellipse.SetXYABP(X0, Y0, A, B, phi);
102  for (Int_t iH = 0; iH < nofHits; iH++){
103  Double_t alfa = gRandom->Rndm()*2.*TMath::Pi();
104 
105  Double_t errorX = gRandom->Gaus(0, sigmaError);
106  Double_t errorY = gRandom->Gaus(0, sigmaError);
107 
108  Double_t hx = A * cos(alfa);
109  Double_t hy = B * sin(alfa);
110 
111  Double_t hitXRot = hx * cos(phi) - hy * sin(phi);
112  Double_t hitYRot = hx * sin(phi) + hy * cos(phi);
113 
114  CbmRichHitLight hit(hitXRot + X0 + errorX, hitYRot + Y0 + errorY);
115  ellipse.AddHit(hit);
116  }
117  // ellipse fit
118  fitEllipse->DoFit(&ellipse);
119  fhErrorA->Fill(A - ellipse.GetAaxis());
120  fhErrorB->Fill(B - ellipse.GetBaxis());
121  fhErrorX->Fill(X0 - ellipse.GetCenterX());
122  fhErrorY->Fill(Y0 - ellipse.GetCenterY());
123  fhErrorPhi->Fill(phi - ellipse.GetPhi());
124  fhA->Fill(ellipse.GetAaxis());
125  fhB->Fill(ellipse.GetBaxis());
126  fhX->Fill(ellipse.GetCenterX());
127  fhY->Fill(ellipse.GetCenterY());
128  fhPhi->Fill(ellipse.GetPhi());
129 
130  // circle fit
131  fitCircle->DoFit(&ellipse);
132  TMatrixD cov(3,3);
133  CalculateFitErrors(&ellipse, sigmaError, cov);
134  Double_t mcR = (A + B) / 2.;
135  fhRadiusErr->Fill(mcR - ellipse.GetRadius());
136  fhCircleXcErr->Fill(X0 - ellipse.GetCenterX());
137  fhCircleYcErr->Fill(Y0 - ellipse.GetCenterY());
138  fhRadius->Fill(ellipse.GetRadius());
139  fhCircleXc->Fill(ellipse.GetCenterX());
140  fhCircleYc->Fill(ellipse.GetCenterY());
141  fhRadiusPool->Fill( (mcR - ellipse.GetRadius()) / sqrt(cov(2,2)) );
142  fhCircleXcPool->Fill( (X0 - ellipse.GetCenterX()) / sqrt(cov(0,0)) );
143  fhCircleYcPool->Fill( (Y0 - ellipse.GetCenterY()) / sqrt(cov(1,1)) );
144  }// iR
145 
146  Draw();
147  cout << nofBadFit << endl;
148 }
149 
150 void CbmRichRingFitterQa::Draw()
151 {
152  TCanvas * c = new TCanvas("rich_fitter_errors", "rich_fitter_errors", 900, 600);
153  c->Divide(3,2);
154  c->cd(1);
155  fhErrorA->Draw();
156  c->cd(2);
157  fhErrorB->Draw();
158  c->cd(3);
159  fhErrorX->Draw();
160  c->cd(4);
161  fhErrorY->Draw();
162  c->cd(5);
163  fhErrorPhi->Draw();
164  cout.precision(4);
165  cout << fhErrorA->GetMean() << " " << fhErrorA->GetRMS() << endl;
166  cout << fhErrorB->GetMean() << " " << fhErrorB->GetRMS() << endl;
167  cout << fhErrorX->GetMean() << " " << fhErrorX->GetRMS() << endl;
168  cout << fhErrorY->GetMean() << " " << fhErrorY->GetRMS() << endl;
169  cout << fhErrorPhi->GetMean() << " " << fhErrorPhi->GetRMS() << endl;
170 
171  TCanvas * c2 = new TCanvas("rich_fitter_params", "rich_fitter_params", 900, 600);
172  c2->Divide(3,2);
173  c2->cd(1);
174  fhA->Draw();
175  c2->cd(2);
176  fhB->Draw();
177  c2->cd(3);
178  fhX->Draw();
179  c2->cd(4);
180  fhY->Draw();
181  c2->cd(5);
182  fhPhi->Draw();
183 
184  TCanvas * c3 = new TCanvas("rich_fitter_circle", "rich_fitter_circle", 900, 900);
185  c3->Divide(3,3);
186  c3->cd(1);
187  fhRadiusErr->Draw();
188  c3->cd(2);
189  fhCircleXcErr->Draw();
190  c3->cd(3);
191  fhCircleYcErr->Draw();
192  c3->cd(4);
193  fhRadius->Draw();
194  c3->cd(5);
195  fhCircleXc->Draw();
196  c3->cd(6);
197  fhCircleYc->Draw();
198  c3->cd(7);
199  fhRadiusPool->Draw();
200  c3->cd(8);
201  fhCircleXcPool->Draw();
202  c3->cd(9);
203  fhCircleYcPool->Draw();
204 
205 }
206 
207 void CbmRichRingFitterQa::CalculateFitErrors(
208  CbmRichRingLight* ring,
209  Double_t sigma,
210  TMatrixD& cov)
211 {
212  //TMatrixD H3(3,3);
213  TMatrixD HY3(3,1);
214  //TMatrixD Cov3(3,3);
215  TMatrixD HC3(3,1);
216 
217  for (Int_t i = 0; i < 3; i++){
218  HY3(i,0) = 0;
219  HC3(i,0) = 0;
220  for (Int_t j = 0; j < 3; j++){
221  //H3(i,j) = 0;
222  cov(i,j) = 0;
223  }
224  }
225 
226  Double_t xc = ring->GetCenterX();
227  Double_t yc = ring->GetCenterY();
228  Double_t R = ring->GetRadius();
229  for (Int_t iHit = 0; iHit < ring->GetNofHits(); iHit++) {
230  Double_t xi = ring->GetHit(iHit).fX;
231  Double_t yi = ring->GetHit(iHit).fY;
232  Double_t ri = sqrt((xi - xc) * (xi - xc) + (yi - yc) * (yi - yc));
233  Double_t err = sigma;
234 
235  Double_t f1 = (-1.0 * (xi - xc)) / (ri * err);
236  Double_t f2 = (-1.0 * (yi - yc))/(ri * err);
237  Double_t f3 = (-1.) / err;
238  Double_t Y = (R - ri) / err;
239 
240  cov(0,0) = cov(0,0) + f1*f1;
241  cov(0,1) = cov(0,1) + f1*f2;
242  cov(0,2) = cov(0,2) + f1*f3;
243 
244  cov(1,0) = cov(0,1);
245  cov(1,1) = cov(1,1) + f2*f2;
246  cov(1,2) = cov(1,2) + f2*f3;
247 
248  cov(2,0) = cov(0,2);
249  cov(2,1) = cov(1,2);
250  cov(2,2) = cov(2,2) + f3*f3;
251 
252  HY3(0,0) = HY3(0,0) + Y*f1;
253  HY3(1,0) = HY3(1,0) + Y*f2;
254  HY3(2,0) = HY3(2,0) + Y*f3;
255  }// iHit
256 
257  //H3.Print();
258  Double_t det = 0.0;
259  cov.Invert(&det);
260  //Cov3 = H3;
261  //H3.Print();
262  //Cov3.Print();
263 
264  //HC3 = H3 * HY3;
265  //HC3.Print();
266 
267  //cout << "dX0= " << HC3(0,0) << " +- " << sqrt(Cov3(0,0)) << endl;
268  // cout << "dY0= " << HC3(1,0) << " +- " << sqrt(Cov3(1,1)) << endl;
269  //cout << "dR= " << HC3(2,0) << " +- " << sqrt(Cov3(2,2)) << endl;
270 
271  //cout << "dX0= " << HC3(0,0) << " +- " << sqrt(Cov3(0,0)) << endl;
272  // cout << "dY0= " << HC3(1,0) << " +- " << sqrt(Cov3(1,1)) << endl;
273  // cout << "dR= " << HC3(2,0) << " +- " << sqrt(Cov3(2,2)) << endl;
274 
275 
276 }
277 
278 ClassImp(CbmRichRingFitterQa);