]>
Commit | Line | Data |
---|---|---|
07627591 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id$ */ | |
17 | ||
18 | /////////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
a6d2bd0c | 20 | // TPC calibration class for parameters which are saved per pad // |
21 | // Each AliTPCCalPad consists of 72 AliTPCCalROC-objects // | |
07627591 | 22 | // // |
23 | /////////////////////////////////////////////////////////////////////////////// | |
24 | ||
25 | #include "AliTPCCalPad.h" | |
26 | #include "AliTPCCalROC.h" | |
184bcc16 | 27 | #include <TObjArray.h> |
28 | #include <TAxis.h> | |
29 | #include <TGraph.h> | |
200be8a6 | 30 | #include <TGraph2D.h> |
31 | #include <TH2F.h> | |
586007f3 | 32 | #include "TTreeStream.h" |
ca5dca67 | 33 | #include "TFile.h" |
34 | #include "TKey.h" | |
72d0ab7e | 35 | #include <iostream> |
586007f3 | 36 | |
07627591 | 37 | ClassImp(AliTPCCalPad) |
38 | ||
39 | //_____________________________________________________________________________ | |
40 | AliTPCCalPad::AliTPCCalPad():TNamed() | |
41 | { | |
42 | // | |
43 | // AliTPCCalPad default constructor | |
44 | // | |
45 | ||
46 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
47 | fROC[isec] = 0; | |
48 | } | |
49 | ||
50 | } | |
51 | ||
52 | //_____________________________________________________________________________ | |
53 | AliTPCCalPad::AliTPCCalPad(const Text_t *name, const Text_t *title) | |
54 | :TNamed(name,title) | |
55 | { | |
56 | // | |
57 | // AliTPCCalPad constructor | |
58 | // | |
59 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
60 | fROC[isec] = new AliTPCCalROC(isec); | |
61 | } | |
62 | } | |
63 | ||
64 | ||
65 | //_____________________________________________________________________________ | |
66 | AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c) | |
67 | { | |
68 | // | |
69 | // AliTPCCalPad copy constructor | |
70 | // | |
71 | ||
e6c51e6e | 72 | for (Int_t isec = 0; isec < kNsec; isec++) { |
7fb8d357 | 73 | fROC[isec] = 0; |
74 | if (c.fROC[isec]) | |
75 | fROC[isec] = new AliTPCCalROC(*(c.fROC[isec])); | |
e6c51e6e | 76 | } |
07627591 | 77 | } |
78 | ||
184bcc16 | 79 | //_____________________________________________________________________________ |
80 | AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed() | |
81 | { | |
82 | // | |
83 | // AliTPCCalPad default constructor | |
84 | // | |
85 | ||
86 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
87 | fROC[isec] = (AliTPCCalROC *)array->At(isec); | |
88 | } | |
89 | ||
90 | } | |
91 | ||
92 | ||
07627591 | 93 | ///_____________________________________________________________________________ |
94 | AliTPCCalPad::~AliTPCCalPad() | |
95 | { | |
96 | // | |
97 | // AliTPCCalPad destructor | |
98 | // | |
99 | ||
100 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
101 | if (fROC[isec]) { | |
102 | delete fROC[isec]; | |
103 | fROC[isec] = 0; | |
104 | } | |
105 | } | |
106 | ||
107 | } | |
108 | ||
109 | //_____________________________________________________________________________ | |
110 | AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c) | |
111 | { | |
112 | // | |
113 | // Assignment operator | |
114 | // | |
115 | ||
116 | if (this != &c) ((AliTPCCalPad &) c).Copy(*this); | |
117 | return *this; | |
118 | ||
119 | } | |
120 | ||
121 | //_____________________________________________________________________________ | |
122 | void AliTPCCalPad::Copy(TObject &c) const | |
123 | { | |
124 | // | |
125 | // Copy function | |
126 | // | |
127 | ||
128 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
129 | if (fROC[isec]) { | |
130 | fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]); | |
131 | } | |
132 | } | |
133 | TObject::Copy(c); | |
134 | } | |
184bcc16 | 135 | |
a6d2bd0c | 136 | |
137 | void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){ | |
138 | // | |
139 | // Set AliTPCCalROC copies values from 'roc' | |
140 | // if sector == -1 the sector specified in 'roc' is used | |
141 | // else sector specified in 'roc' is ignored and specified sector is filled | |
142 | // | |
143 | if (sector == -1) sector = roc->GetSector(); | |
72d0ab7e | 144 | if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector); |
a6d2bd0c | 145 | for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++) |
146 | fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel)); | |
147 | } | |
148 | ||
149 | ||
150 | ||
90127643 | 151 | //_____________________________________________________________________________ |
152 | void AliTPCCalPad::Add(Float_t c1) | |
153 | { | |
154 | // | |
a6d2bd0c | 155 | // add constant c1 to all channels of all ROCs |
90127643 | 156 | // |
157 | ||
158 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
159 | if (fROC[isec]){ | |
160 | fROC[isec]->Add(c1); | |
161 | } | |
162 | } | |
163 | } | |
164 | ||
165 | //_____________________________________________________________________________ | |
166 | void AliTPCCalPad::Multiply(Float_t c1) | |
167 | { | |
a6d2bd0c | 168 | // |
169 | // multiply each channel of all ROCs with c1 | |
170 | // | |
90127643 | 171 | for (Int_t isec = 0; isec < kNsec; isec++) { |
172 | if (fROC[isec]){ | |
173 | fROC[isec]->Multiply(c1); | |
174 | } | |
175 | } | |
176 | } | |
177 | ||
178 | //_____________________________________________________________________________ | |
179 | void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1) | |
180 | { | |
a6d2bd0c | 181 | // |
182 | // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs | |
183 | // - pad by pad - | |
184 | // | |
90127643 | 185 | for (Int_t isec = 0; isec < kNsec; isec++) { |
79c0f359 | 186 | if (fROC[isec] && pad->GetCalROC(isec)){ |
90127643 | 187 | fROC[isec]->Add(pad->GetCalROC(isec),c1); |
188 | } | |
189 | } | |
190 | } | |
191 | ||
192 | //_____________________________________________________________________________ | |
193 | void AliTPCCalPad::Multiply(const AliTPCCalPad * pad) | |
194 | { | |
a6d2bd0c | 195 | // |
196 | // multiply each channel of all ROCs with the coresponding channel of 'pad' | |
197 | // - pad by pad - | |
198 | // | |
199 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
90127643 | 200 | if (fROC[isec]){ |
201 | fROC[isec]->Multiply(pad->GetCalROC(isec)); | |
202 | } | |
203 | } | |
204 | } | |
205 | ||
206 | //_____________________________________________________________________________ | |
207 | void AliTPCCalPad::Divide(const AliTPCCalPad * pad) | |
208 | { | |
a6d2bd0c | 209 | // |
210 | // divide each channel of all ROCs by the coresponding channel of 'pad' | |
211 | // - pad by pad - | |
212 | // | |
90127643 | 213 | for (Int_t isec = 0; isec < kNsec; isec++) { |
214 | if (fROC[isec]){ | |
215 | fROC[isec]->Divide(pad->GetCalROC(isec)); | |
216 | } | |
217 | } | |
218 | } | |
219 | ||
220 | //_____________________________________________________________________________ | |
184bcc16 | 221 | TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){ |
222 | // | |
223 | // type=1 - mean | |
224 | // 2 - median | |
225 | // 3 - LTM | |
a6d2bd0c | 226 | // |
184bcc16 | 227 | Int_t npoints = 0; |
228 | for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++; | |
229 | TGraph * graph = new TGraph(npoints); | |
230 | npoints=0; | |
231 | for (Int_t isec=0;isec<72;isec++){ | |
232 | if (!fROC[isec]) continue; | |
233 | if (type==0) graph->SetPoint(npoints,isec,fROC[isec]->GetMean()); | |
234 | if (type==1) graph->SetPoint(npoints,isec,fROC[isec]->GetMedian()); | |
235 | if (type==2) graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio)); | |
236 | npoints++; | |
237 | } | |
238 | ||
239 | graph->GetXaxis()->SetTitle("Sector"); | |
240 | if (type==0) { | |
241 | graph->GetYaxis()->SetTitle("Mean"); | |
242 | graph->SetMarkerStyle(22); | |
243 | } | |
244 | if (type==1) { | |
245 | graph->GetYaxis()->SetTitle("Median"); | |
246 | graph->SetMarkerStyle(22); | |
247 | } | |
248 | if (type==2) { | |
249 | graph->GetYaxis()->SetTitle(Form("Mean%f",ratio)); | |
250 | graph->SetMarkerStyle(24); | |
251 | } | |
252 | ||
253 | return graph; | |
254 | } | |
200be8a6 | 255 | |
90127643 | 256 | //_____________________________________________________________________________ |
257 | Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms) | |
258 | { | |
259 | // | |
a6d2bd0c | 260 | // Calculates mean and RMS of all ROCs |
90127643 | 261 | // |
262 | Double_t sum = 0, sum2 = 0, n=0, val=0; | |
263 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
264 | AliTPCCalROC *calRoc = fROC[isec]; | |
265 | if ( calRoc ){ | |
266 | for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){ | |
267 | for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){ | |
268 | val = calRoc->GetValue(irow,ipad); | |
269 | sum+=val; | |
270 | sum2+=val*val; | |
271 | n++; | |
272 | } | |
273 | } | |
274 | ||
275 | } | |
276 | } | |
277 | Double_t n1 = 1./n; | |
278 | Double_t mean = sum*n1; | |
279 | rms = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean)); | |
280 | return mean; | |
281 | } | |
282 | ||
283 | ||
284 | //_____________________________________________________________________________ | |
ca5dca67 | 285 | Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad) |
90127643 | 286 | { |
287 | // | |
288 | // return mean of the mean of all ROCs | |
289 | // | |
290 | Double_t arr[kNsec]; | |
291 | Int_t n=0; | |
292 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
ca5dca67 | 293 | AliTPCCalROC *calRoc = fROC[isec]; |
294 | if ( calRoc ){ | |
295 | AliTPCCalROC* outlierROC = 0; | |
296 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
297 | arr[n] = calRoc->GetMean(outlierROC); | |
298 | n++; | |
299 | } | |
90127643 | 300 | } |
301 | return TMath::Mean(n,arr); | |
302 | } | |
303 | ||
304 | //_____________________________________________________________________________ | |
ca5dca67 | 305 | Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad) |
90127643 | 306 | { |
307 | // | |
308 | // return mean of the RMS of all ROCs | |
309 | // | |
310 | Double_t arr[kNsec]; | |
311 | Int_t n=0; | |
312 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
ca5dca67 | 313 | AliTPCCalROC *calRoc = fROC[isec]; |
314 | if ( calRoc ){ | |
315 | AliTPCCalROC* outlierROC = 0; | |
316 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
317 | arr[n] = calRoc->GetRMS(outlierROC); | |
318 | n++; | |
319 | } | |
90127643 | 320 | } |
321 | return TMath::Mean(n,arr); | |
322 | } | |
323 | ||
324 | //_____________________________________________________________________________ | |
ca5dca67 | 325 | Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad) |
90127643 | 326 | { |
327 | // | |
328 | // return mean of the median of all ROCs | |
329 | // | |
330 | Double_t arr[kNsec]; | |
331 | Int_t n=0; | |
332 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
ca5dca67 | 333 | AliTPCCalROC *calRoc = fROC[isec]; |
334 | if ( calRoc ){ | |
335 | AliTPCCalROC* outlierROC = 0; | |
336 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
337 | arr[n] = calRoc->GetMedian(outlierROC); | |
338 | n++; | |
339 | } | |
90127643 | 340 | } |
341 | return TMath::Mean(n,arr); | |
342 | } | |
343 | ||
344 | //_____________________________________________________________________________ | |
ca5dca67 | 345 | Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad) |
90127643 | 346 | { |
347 | // | |
348 | // return mean of the LTM and sigma of all ROCs | |
349 | // | |
350 | Double_t arrm[kNsec]; | |
351 | Double_t arrs[kNsec]; | |
352 | Double_t *sTemp=0x0; | |
353 | Int_t n=0; | |
354 | ||
355 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
356 | AliTPCCalROC *calRoc = fROC[isec]; | |
357 | if ( calRoc ){ | |
358 | if ( sigma ) sTemp=arrs+n; | |
ca5dca67 | 359 | AliTPCCalROC* outlierROC = 0; |
360 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
361 | arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC); | |
90127643 | 362 | n++; |
363 | } | |
364 | } | |
365 | if ( sigma ) *sigma = TMath::Mean(n,arrs); | |
366 | return TMath::Mean(n,arrm); | |
367 | } | |
368 | ||
369 | //_____________________________________________________________________________ | |
370 | TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){ | |
371 | // | |
372 | // make 1D histo | |
373 | // type -1 = user defined range | |
374 | // 0 = nsigma cut nsigma=min | |
a6d2bd0c | 375 | // |
90127643 | 376 | if (type>=0){ |
377 | if (type==0){ | |
378 | // nsigma range | |
379 | Float_t mean = GetMean(); | |
380 | Float_t sigma = GetRMS(); | |
381 | Float_t nsigma = TMath::Abs(min); | |
382 | min = mean-nsigma*sigma; | |
383 | max = mean+nsigma*sigma; | |
384 | } | |
385 | if (type==1){ | |
386 | // fixed range | |
387 | Float_t mean = GetMedian(); | |
388 | Float_t delta = min; | |
389 | min = mean-delta; | |
390 | max = mean+delta; | |
391 | } | |
392 | if (type==2){ | |
393 | // | |
394 | // LTM mean +- nsigma | |
395 | // | |
396 | Double_t sigma; | |
397 | Float_t mean = GetLTM(&sigma,max); | |
398 | sigma*=min; | |
399 | min = mean-sigma; | |
400 | max = mean+sigma; | |
401 | } | |
402 | } | |
403 | char name[1000]; | |
404 | sprintf(name,"%s Pad 1D",GetTitle()); | |
405 | TH1F * his = new TH1F(name,name,100, min,max); | |
406 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
407 | if (fROC[isec]){ | |
408 | for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){ | |
409 | UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow); | |
410 | for (UInt_t ipad=0; ipad<npads; ipad++){ | |
411 | his->Fill(fROC[isec]->GetValue(irow,ipad)); | |
412 | } | |
413 | } | |
414 | } | |
415 | } | |
416 | return his; | |
417 | } | |
418 | ||
419 | //_____________________________________________________________________________ | |
200be8a6 | 420 | TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){ |
421 | // | |
422 | // Make 2D graph | |
423 | // side - specify the side A = 0 C = 1 | |
424 | // type - used types of determination of boundaries in z | |
a6d2bd0c | 425 | // |
200be8a6 | 426 | Float_t kEpsilon = 0.000000000001; |
427 | TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250); | |
428 | AliTPCROC * roc = AliTPCROC::Instance(); | |
429 | for (Int_t isec=0; isec<72; isec++){ | |
430 | if (side==0 && isec%36>=18) continue; | |
431 | if (side>0 && isec%36<18) continue; | |
432 | if (fROC[isec]){ | |
433 | AliTPCCalROC * calRoc = fROC[isec]; | |
434 | for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++) | |
435 | for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++) | |
436 | if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){ | |
437 | Float_t xyz[3]; | |
438 | roc->GetPositionGlobal(isec,irow,ipad,xyz); | |
439 | Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5); | |
440 | Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5); | |
441 | Float_t value = calRoc->GetValue(irow,ipad); | |
442 | his->SetBinContent(binx,biny,value); | |
443 | } | |
444 | } | |
445 | } | |
446 | his->SetXTitle("x (cm)"); | |
447 | his->SetYTitle("y (cm)"); | |
448 | return his; | |
449 | } | |
450 | ||
451 | ||
72d0ab7e | 452 | AliTPCCalPad* AliTPCCalPad::LocalFit(const char* padName, Int_t rowRadius, Int_t padRadius, AliTPCCalPad* PadOutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction, Bool_t printCurrentSector) const { |
a6d2bd0c | 453 | // |
454 | // Loops over all AliTPCCalROCs and performs a localFit in each ROC | |
455 | // AliTPCCalPad with fit-data is returned | |
456 | // rowRadius and padRadius specifies a window around a given pad in one sector. | |
457 | // The data of this window are fitted with a parabolic function. | |
458 | // This function is evaluated at the pad's position. | |
459 | // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window. | |
460 | // rowRadius - radius - rows to be used for smoothing | |
461 | // padradius - radius - pads to be used for smoothing | |
462 | // ROCoutlier - map of outliers - pads not to be used for local smoothing | |
463 | // robust - robust method of fitting - (much slower) | |
464 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
465 | // robustFraction: Fraction of data that will be used in EvalRobust | |
466 | // | |
467 | // | |
468 | AliTPCCalPad* pad = new AliTPCCalPad(padName, padName); | |
469 | for (Int_t isec = 0; isec < 72; isec++){ | |
72d0ab7e | 470 | if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush; |
a6d2bd0c | 471 | if (PadOutliers) |
72d0ab7e | 472 | pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction)); |
a6d2bd0c | 473 | else |
72d0ab7e | 474 | pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction)); |
a6d2bd0c | 475 | } |
476 | return pad; | |
477 | } | |
478 | ||
479 | ||
6d8681ef | 480 | AliTPCCalPad* AliTPCCalPad::GlobalFit(const char* padName, AliTPCCalPad* PadOutliers, Bool_t robust, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err, TObjArray *fitParArr, TObjArray *fitCovArr){ |
a6d2bd0c | 481 | // |
482 | // Loops over all AliTPCCalROCs and performs a globalFit in each ROC | |
483 | // AliTPCCalPad with fit-data is returned | |
484 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
485 | // robustFraction: Fraction of data that will be used in EvalRobust | |
486 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
487 | // robustFraction: Fraction of data that will be used in EvalRobust | |
6d8681ef | 488 | // err: error of the data points |
489 | // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array | |
a6d2bd0c | 490 | // |
491 | AliTPCCalPad* pad = new AliTPCCalPad(padName, padName); | |
492 | TVectorD fitParam(0); | |
493 | TMatrixD covMatrix(0,0); | |
494 | Float_t chi2 = 0; | |
495 | for (Int_t isec = 0; isec < 72; isec++){ | |
496 | if (PadOutliers) | |
6d8681ef | 497 | GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err); |
a6d2bd0c | 498 | else |
6d8681ef | 499 | GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err); |
500 | ||
501 | AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec); | |
502 | pad->SetCalROC(roc); | |
503 | delete roc; | |
504 | if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec); | |
505 | if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec); | |
a6d2bd0c | 506 | } |
507 | return pad; | |
508 | } | |
509 | ||
510 | ||
511 | void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, Int_t fitType, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction){ | |
512 | // | |
513 | // Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side | |
514 | // fitType == 0: fit plane function | |
515 | // fitType == 1: fit parabolic function | |
516 | // PadOutliers - pads with value !=0 are not used in fitting procedure | |
517 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
518 | // robustFraction: Fraction of data that will be used in EvalRobust | |
519 | // | |
520 | TLinearFitter* fitterGA = 0; | |
521 | TLinearFitter* fitterGC = 0; | |
522 | ||
523 | if (fitType == 1) { | |
524 | fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5"); | |
525 | fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5"); | |
526 | } | |
527 | else { | |
528 | fitterGA = new TLinearFitter(3,"x0++x1++x2"); | |
529 | fitterGC = new TLinearFitter(3,"x0++x1++x2"); | |
530 | } | |
531 | fitterGA->StoreData(kTRUE); | |
532 | fitterGC->StoreData(kTRUE); | |
533 | fitterGA->ClearPoints(); | |
534 | fitterGC->ClearPoints(); | |
535 | Double_t xx[6]; | |
536 | Int_t npointsA=0; | |
537 | Int_t npointsC=0; | |
538 | ||
539 | Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position | |
540 | Float_t lx, ly; // pads position | |
541 | ||
542 | AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position | |
543 | ||
544 | // loop over all sectors and pads and read data into fitterGA and fitterGC | |
545 | if (fitType == 1) { | |
546 | // parabolic fit | |
547 | fitParamSideA.ResizeTo(6); | |
548 | fitParamSideC.ResizeTo(6); | |
549 | covMatrixSideA.ResizeTo(6,6); | |
550 | covMatrixSideC.ResizeTo(6,6); | |
551 | for (UInt_t isec = 0; isec<72; isec++){ | |
552 | for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { | |
553 | for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { | |
554 | // fill fitterG | |
555 | tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number | |
556 | lx = localXY[0]; | |
557 | ly = localXY[1]; | |
558 | xx[0] = 1; | |
559 | xx[1] = lx; | |
560 | xx[2] = ly; | |
561 | xx[3] = lx*lx; | |
562 | xx[4] = ly*ly; | |
563 | xx[5] = lx*ly; | |
564 | if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) { | |
565 | // if given pad is no outlier, add it to TLinearFitter, decide to which of both | |
566 | // sector 0 - 17: IROC, A | |
567 | // sector 18 - 35: IROC, C | |
568 | // sector 36 - 53: OROC, A | |
569 | // sector 54 - 71: CROC, C | |
570 | if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A | |
571 | npointsA++; | |
572 | fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
573 | } | |
574 | else { // side C | |
575 | npointsC++; | |
576 | fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
577 | } | |
578 | } | |
579 | } | |
580 | } | |
581 | } | |
582 | } | |
583 | else { | |
584 | // linear fit | |
585 | fitParamSideA.ResizeTo(3); | |
586 | fitParamSideC.ResizeTo(3); | |
587 | covMatrixSideA.ResizeTo(3,3); | |
588 | covMatrixSideC.ResizeTo(3,3); | |
589 | ||
590 | for (UInt_t isec = 0; isec<72; isec++){ | |
591 | for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { | |
592 | for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { | |
593 | // fill fitterG | |
594 | tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number | |
595 | lx = localXY[0]; | |
596 | ly = localXY[1]; | |
597 | xx[0] = 1; | |
598 | xx[1] = lx; | |
599 | xx[2] = ly; | |
600 | if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) { | |
601 | // if given pad is no outlier, add it to TLinearFitter, decide to which of both | |
602 | // sector 0 - 17: IROC, A | |
603 | // sector 18 - 35: IROC, C | |
604 | // sector 36 - 53: OROC, A | |
605 | // sector 54 - 71: CROC, C | |
606 | if (isec <= 17 || (isec >= 36 && isec <= 53)) { | |
607 | // Side A | |
608 | npointsA++; | |
609 | fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
610 | } | |
611 | else { | |
612 | // side C | |
613 | npointsC++; | |
614 | fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
615 | } | |
616 | } | |
617 | } | |
618 | } | |
619 | } | |
620 | } | |
621 | ||
622 | fitterGA->Eval(); | |
623 | fitterGC->Eval(); | |
624 | fitterGA->GetParameters(fitParamSideA); | |
625 | fitterGC->GetParameters(fitParamSideC); | |
626 | fitterGA->GetCovarianceMatrix(covMatrixSideA); | |
627 | fitterGC->GetCovarianceMatrix(covMatrixSideC); | |
628 | if (fitType == 1){ | |
629 | chi2SideA = fitterGA->GetChisquare()/(npointsA-6.); | |
630 | chi2SideC = fitterGC->GetChisquare()/(npointsC-6.); | |
631 | } | |
632 | else { | |
633 | chi2SideA = fitterGA->GetChisquare()/(npointsA-3.); | |
634 | chi2SideC = fitterGC->GetChisquare()/(npointsC-3.); | |
635 | } | |
636 | if (robust && chi2SideA > chi2Threshold) { | |
637 | // std::cout << "robust fitter called... " << std::endl; | |
638 | fitterGA->EvalRobust(robustFraction); | |
639 | fitterGA->GetParameters(fitParamSideA); | |
640 | } | |
641 | if (robust && chi2SideC > chi2Threshold) { | |
642 | // std::cout << "robust fitter called... " << std::endl; | |
643 | fitterGC->EvalRobust(robustFraction); | |
644 | fitterGC->GetParameters(fitParamSideC); | |
645 | } | |
646 | delete fitterGA; | |
647 | delete fitterGC; | |
648 | } | |
200be8a6 | 649 | |
586007f3 | 650 |