]>
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" | |
5312f439 | 35 | #include <TFormula.h> |
36 | #include <TString.h> | |
37 | #include <TObjString.h> | |
72d0ab7e | 38 | #include <iostream> |
732e90a8 | 39 | #include <AliLog.h> |
586007f3 | 40 | |
07627591 | 41 | ClassImp(AliTPCCalPad) |
42 | ||
43 | //_____________________________________________________________________________ | |
44 | AliTPCCalPad::AliTPCCalPad():TNamed() | |
45 | { | |
46 | // | |
47 | // AliTPCCalPad default constructor | |
48 | // | |
49 | ||
50 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
51 | fROC[isec] = 0; | |
52 | } | |
53 | ||
54 | } | |
55 | ||
56 | //_____________________________________________________________________________ | |
57 | AliTPCCalPad::AliTPCCalPad(const Text_t *name, const Text_t *title) | |
58 | :TNamed(name,title) | |
59 | { | |
60 | // | |
61 | // AliTPCCalPad constructor | |
62 | // | |
63 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
64 | fROC[isec] = new AliTPCCalROC(isec); | |
65 | } | |
66 | } | |
67 | ||
68 | ||
69 | //_____________________________________________________________________________ | |
70 | AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c) | |
71 | { | |
72 | // | |
73 | // AliTPCCalPad copy constructor | |
74 | // | |
75 | ||
e6c51e6e | 76 | for (Int_t isec = 0; isec < kNsec; isec++) { |
7fb8d357 | 77 | fROC[isec] = 0; |
78 | if (c.fROC[isec]) | |
79 | fROC[isec] = new AliTPCCalROC(*(c.fROC[isec])); | |
e6c51e6e | 80 | } |
07627591 | 81 | } |
82 | ||
184bcc16 | 83 | //_____________________________________________________________________________ |
84 | AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed() | |
85 | { | |
86 | // | |
87 | // AliTPCCalPad default constructor | |
88 | // | |
89 | ||
90 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
91 | fROC[isec] = (AliTPCCalROC *)array->At(isec); | |
92 | } | |
93 | ||
94 | } | |
95 | ||
96 | ||
07627591 | 97 | ///_____________________________________________________________________________ |
98 | AliTPCCalPad::~AliTPCCalPad() | |
99 | { | |
100 | // | |
101 | // AliTPCCalPad destructor | |
102 | // | |
103 | ||
104 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
105 | if (fROC[isec]) { | |
106 | delete fROC[isec]; | |
107 | fROC[isec] = 0; | |
108 | } | |
109 | } | |
110 | ||
111 | } | |
112 | ||
113 | //_____________________________________________________________________________ | |
114 | AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c) | |
115 | { | |
116 | // | |
117 | // Assignment operator | |
118 | // | |
119 | ||
120 | if (this != &c) ((AliTPCCalPad &) c).Copy(*this); | |
121 | return *this; | |
122 | ||
123 | } | |
124 | ||
125 | //_____________________________________________________________________________ | |
126 | void AliTPCCalPad::Copy(TObject &c) const | |
127 | { | |
128 | // | |
129 | // Copy function | |
130 | // | |
131 | ||
132 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
133 | if (fROC[isec]) { | |
134 | fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]); | |
135 | } | |
136 | } | |
137 | TObject::Copy(c); | |
138 | } | |
184bcc16 | 139 | |
a6d2bd0c | 140 | |
141 | void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){ | |
142 | // | |
143 | // Set AliTPCCalROC copies values from 'roc' | |
144 | // if sector == -1 the sector specified in 'roc' is used | |
145 | // else sector specified in 'roc' is ignored and specified sector is filled | |
146 | // | |
147 | if (sector == -1) sector = roc->GetSector(); | |
72d0ab7e | 148 | if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector); |
a6d2bd0c | 149 | for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++) |
150 | fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel)); | |
151 | } | |
152 | ||
153 | ||
154 | ||
90127643 | 155 | //_____________________________________________________________________________ |
156 | void AliTPCCalPad::Add(Float_t c1) | |
157 | { | |
158 | // | |
a6d2bd0c | 159 | // add constant c1 to all channels of all ROCs |
90127643 | 160 | // |
161 | ||
162 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
163 | if (fROC[isec]){ | |
164 | fROC[isec]->Add(c1); | |
165 | } | |
166 | } | |
167 | } | |
168 | ||
169 | //_____________________________________________________________________________ | |
170 | void AliTPCCalPad::Multiply(Float_t c1) | |
171 | { | |
a6d2bd0c | 172 | // |
173 | // multiply each channel of all ROCs with c1 | |
174 | // | |
90127643 | 175 | for (Int_t isec = 0; isec < kNsec; isec++) { |
176 | if (fROC[isec]){ | |
177 | fROC[isec]->Multiply(c1); | |
178 | } | |
179 | } | |
180 | } | |
181 | ||
182 | //_____________________________________________________________________________ | |
183 | void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1) | |
184 | { | |
a6d2bd0c | 185 | // |
186 | // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs | |
187 | // - pad by pad - | |
188 | // | |
90127643 | 189 | for (Int_t isec = 0; isec < kNsec; isec++) { |
79c0f359 | 190 | if (fROC[isec] && pad->GetCalROC(isec)){ |
90127643 | 191 | fROC[isec]->Add(pad->GetCalROC(isec),c1); |
192 | } | |
193 | } | |
194 | } | |
195 | ||
196 | //_____________________________________________________________________________ | |
197 | void AliTPCCalPad::Multiply(const AliTPCCalPad * pad) | |
198 | { | |
a6d2bd0c | 199 | // |
200 | // multiply each channel of all ROCs with the coresponding channel of 'pad' | |
201 | // - pad by pad - | |
202 | // | |
203 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
90127643 | 204 | if (fROC[isec]){ |
205 | fROC[isec]->Multiply(pad->GetCalROC(isec)); | |
206 | } | |
207 | } | |
208 | } | |
209 | ||
210 | //_____________________________________________________________________________ | |
211 | void AliTPCCalPad::Divide(const AliTPCCalPad * pad) | |
212 | { | |
a6d2bd0c | 213 | // |
214 | // divide each channel of all ROCs by the coresponding channel of 'pad' | |
215 | // - pad by pad - | |
216 | // | |
90127643 | 217 | for (Int_t isec = 0; isec < kNsec; isec++) { |
218 | if (fROC[isec]){ | |
219 | fROC[isec]->Divide(pad->GetCalROC(isec)); | |
220 | } | |
221 | } | |
222 | } | |
223 | ||
224 | //_____________________________________________________________________________ | |
184bcc16 | 225 | TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){ |
226 | // | |
227 | // type=1 - mean | |
228 | // 2 - median | |
229 | // 3 - LTM | |
a6d2bd0c | 230 | // |
184bcc16 | 231 | Int_t npoints = 0; |
232 | for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++; | |
233 | TGraph * graph = new TGraph(npoints); | |
234 | npoints=0; | |
235 | for (Int_t isec=0;isec<72;isec++){ | |
236 | if (!fROC[isec]) continue; | |
237 | if (type==0) graph->SetPoint(npoints,isec,fROC[isec]->GetMean()); | |
238 | if (type==1) graph->SetPoint(npoints,isec,fROC[isec]->GetMedian()); | |
239 | if (type==2) graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio)); | |
240 | npoints++; | |
241 | } | |
242 | ||
243 | graph->GetXaxis()->SetTitle("Sector"); | |
244 | if (type==0) { | |
245 | graph->GetYaxis()->SetTitle("Mean"); | |
246 | graph->SetMarkerStyle(22); | |
247 | } | |
248 | if (type==1) { | |
249 | graph->GetYaxis()->SetTitle("Median"); | |
250 | graph->SetMarkerStyle(22); | |
251 | } | |
252 | if (type==2) { | |
253 | graph->GetYaxis()->SetTitle(Form("Mean%f",ratio)); | |
254 | graph->SetMarkerStyle(24); | |
255 | } | |
256 | ||
257 | return graph; | |
258 | } | |
200be8a6 | 259 | |
90127643 | 260 | //_____________________________________________________________________________ |
261 | Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms) | |
262 | { | |
263 | // | |
a6d2bd0c | 264 | // Calculates mean and RMS of all ROCs |
90127643 | 265 | // |
266 | Double_t sum = 0, sum2 = 0, n=0, val=0; | |
267 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
268 | AliTPCCalROC *calRoc = fROC[isec]; | |
269 | if ( calRoc ){ | |
270 | for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){ | |
271 | for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){ | |
272 | val = calRoc->GetValue(irow,ipad); | |
273 | sum+=val; | |
274 | sum2+=val*val; | |
275 | n++; | |
276 | } | |
277 | } | |
278 | ||
279 | } | |
280 | } | |
281 | Double_t n1 = 1./n; | |
282 | Double_t mean = sum*n1; | |
283 | rms = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean)); | |
284 | return mean; | |
285 | } | |
286 | ||
287 | ||
288 | //_____________________________________________________________________________ | |
ca5dca67 | 289 | Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad) |
90127643 | 290 | { |
291 | // | |
292 | // return mean of the mean of all ROCs | |
293 | // | |
294 | Double_t arr[kNsec]; | |
295 | Int_t n=0; | |
296 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
ca5dca67 | 297 | AliTPCCalROC *calRoc = fROC[isec]; |
298 | if ( calRoc ){ | |
299 | AliTPCCalROC* outlierROC = 0; | |
300 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
301 | arr[n] = calRoc->GetMean(outlierROC); | |
302 | n++; | |
303 | } | |
90127643 | 304 | } |
305 | return TMath::Mean(n,arr); | |
306 | } | |
307 | ||
308 | //_____________________________________________________________________________ | |
ca5dca67 | 309 | Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad) |
90127643 | 310 | { |
311 | // | |
312 | // return mean of the RMS of all ROCs | |
313 | // | |
314 | Double_t arr[kNsec]; | |
315 | Int_t n=0; | |
316 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
ca5dca67 | 317 | AliTPCCalROC *calRoc = fROC[isec]; |
318 | if ( calRoc ){ | |
319 | AliTPCCalROC* outlierROC = 0; | |
320 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
321 | arr[n] = calRoc->GetRMS(outlierROC); | |
322 | n++; | |
323 | } | |
90127643 | 324 | } |
325 | return TMath::Mean(n,arr); | |
326 | } | |
327 | ||
328 | //_____________________________________________________________________________ | |
ca5dca67 | 329 | Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad) |
90127643 | 330 | { |
331 | // | |
332 | // return mean of the median of all ROCs | |
333 | // | |
334 | Double_t arr[kNsec]; | |
335 | Int_t n=0; | |
336 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
ca5dca67 | 337 | AliTPCCalROC *calRoc = fROC[isec]; |
338 | if ( calRoc ){ | |
339 | AliTPCCalROC* outlierROC = 0; | |
340 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
341 | arr[n] = calRoc->GetMedian(outlierROC); | |
342 | n++; | |
343 | } | |
90127643 | 344 | } |
345 | return TMath::Mean(n,arr); | |
346 | } | |
347 | ||
348 | //_____________________________________________________________________________ | |
ca5dca67 | 349 | Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad) |
90127643 | 350 | { |
351 | // | |
352 | // return mean of the LTM and sigma of all ROCs | |
353 | // | |
354 | Double_t arrm[kNsec]; | |
355 | Double_t arrs[kNsec]; | |
356 | Double_t *sTemp=0x0; | |
357 | Int_t n=0; | |
358 | ||
359 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
360 | AliTPCCalROC *calRoc = fROC[isec]; | |
361 | if ( calRoc ){ | |
362 | if ( sigma ) sTemp=arrs+n; | |
ca5dca67 | 363 | AliTPCCalROC* outlierROC = 0; |
364 | if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); | |
365 | arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC); | |
90127643 | 366 | n++; |
367 | } | |
368 | } | |
369 | if ( sigma ) *sigma = TMath::Mean(n,arrs); | |
370 | return TMath::Mean(n,arrm); | |
371 | } | |
372 | ||
373 | //_____________________________________________________________________________ | |
374 | TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){ | |
375 | // | |
376 | // make 1D histo | |
377 | // type -1 = user defined range | |
378 | // 0 = nsigma cut nsigma=min | |
a6d2bd0c | 379 | // |
90127643 | 380 | if (type>=0){ |
381 | if (type==0){ | |
382 | // nsigma range | |
383 | Float_t mean = GetMean(); | |
384 | Float_t sigma = GetRMS(); | |
385 | Float_t nsigma = TMath::Abs(min); | |
386 | min = mean-nsigma*sigma; | |
387 | max = mean+nsigma*sigma; | |
388 | } | |
389 | if (type==1){ | |
390 | // fixed range | |
391 | Float_t mean = GetMedian(); | |
392 | Float_t delta = min; | |
393 | min = mean-delta; | |
394 | max = mean+delta; | |
395 | } | |
396 | if (type==2){ | |
397 | // | |
398 | // LTM mean +- nsigma | |
399 | // | |
400 | Double_t sigma; | |
401 | Float_t mean = GetLTM(&sigma,max); | |
402 | sigma*=min; | |
403 | min = mean-sigma; | |
404 | max = mean+sigma; | |
405 | } | |
406 | } | |
407 | char name[1000]; | |
408 | sprintf(name,"%s Pad 1D",GetTitle()); | |
409 | TH1F * his = new TH1F(name,name,100, min,max); | |
410 | for (Int_t isec = 0; isec < kNsec; isec++) { | |
411 | if (fROC[isec]){ | |
412 | for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){ | |
413 | UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow); | |
414 | for (UInt_t ipad=0; ipad<npads; ipad++){ | |
415 | his->Fill(fROC[isec]->GetValue(irow,ipad)); | |
416 | } | |
417 | } | |
418 | } | |
419 | } | |
420 | return his; | |
421 | } | |
422 | ||
423 | //_____________________________________________________________________________ | |
200be8a6 | 424 | TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){ |
425 | // | |
426 | // Make 2D graph | |
427 | // side - specify the side A = 0 C = 1 | |
428 | // type - used types of determination of boundaries in z | |
a6d2bd0c | 429 | // |
200be8a6 | 430 | Float_t kEpsilon = 0.000000000001; |
431 | TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250); | |
432 | AliTPCROC * roc = AliTPCROC::Instance(); | |
433 | for (Int_t isec=0; isec<72; isec++){ | |
434 | if (side==0 && isec%36>=18) continue; | |
435 | if (side>0 && isec%36<18) continue; | |
436 | if (fROC[isec]){ | |
437 | AliTPCCalROC * calRoc = fROC[isec]; | |
438 | for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++) | |
439 | for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++) | |
440 | if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){ | |
441 | Float_t xyz[3]; | |
442 | roc->GetPositionGlobal(isec,irow,ipad,xyz); | |
443 | Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5); | |
444 | Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5); | |
445 | Float_t value = calRoc->GetValue(irow,ipad); | |
446 | his->SetBinContent(binx,biny,value); | |
447 | } | |
448 | } | |
449 | } | |
450 | his->SetXTitle("x (cm)"); | |
451 | his->SetYTitle("y (cm)"); | |
452 | return his; | |
453 | } | |
454 | ||
455 | ||
72d0ab7e | 456 | 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 | 457 | // |
458 | // Loops over all AliTPCCalROCs and performs a localFit in each ROC | |
459 | // AliTPCCalPad with fit-data is returned | |
460 | // rowRadius and padRadius specifies a window around a given pad in one sector. | |
461 | // The data of this window are fitted with a parabolic function. | |
462 | // This function is evaluated at the pad's position. | |
463 | // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window. | |
464 | // rowRadius - radius - rows to be used for smoothing | |
465 | // padradius - radius - pads to be used for smoothing | |
466 | // ROCoutlier - map of outliers - pads not to be used for local smoothing | |
467 | // robust - robust method of fitting - (much slower) | |
468 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
469 | // robustFraction: Fraction of data that will be used in EvalRobust | |
470 | // | |
471 | // | |
472 | AliTPCCalPad* pad = new AliTPCCalPad(padName, padName); | |
473 | for (Int_t isec = 0; isec < 72; isec++){ | |
72d0ab7e | 474 | if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush; |
a6d2bd0c | 475 | if (PadOutliers) |
72d0ab7e | 476 | pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction)); |
a6d2bd0c | 477 | else |
72d0ab7e | 478 | pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction)); |
a6d2bd0c | 479 | } |
480 | return pad; | |
481 | } | |
482 | ||
483 | ||
6d8681ef | 484 | 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 | 485 | // |
486 | // Loops over all AliTPCCalROCs and performs a globalFit in each ROC | |
487 | // AliTPCCalPad with fit-data is returned | |
488 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
489 | // robustFraction: Fraction of data that will be used in EvalRobust | |
490 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
491 | // robustFraction: Fraction of data that will be used in EvalRobust | |
6d8681ef | 492 | // err: error of the data points |
493 | // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array | |
a6d2bd0c | 494 | // |
495 | AliTPCCalPad* pad = new AliTPCCalPad(padName, padName); | |
496 | TVectorD fitParam(0); | |
497 | TMatrixD covMatrix(0,0); | |
498 | Float_t chi2 = 0; | |
499 | for (Int_t isec = 0; isec < 72; isec++){ | |
500 | if (PadOutliers) | |
6d8681ef | 501 | GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err); |
a6d2bd0c | 502 | else |
6d8681ef | 503 | GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err); |
504 | ||
505 | AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec); | |
506 | pad->SetCalROC(roc); | |
507 | delete roc; | |
508 | if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec); | |
509 | if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec); | |
a6d2bd0c | 510 | } |
511 | return pad; | |
512 | } | |
732e90a8 | 513 | //_____________________________________________________________________________ |
514 | TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula) | |
515 | { | |
5312f439 | 516 | // |
732e90a8 | 517 | // create an array of TFormulas for the each parameter of the fit function |
5312f439 | 518 | // |
519 | ||
520 | // split fit string in single parameters | |
521 | // find dimension of the fit: | |
522 | TString fitString(fitFormula); | |
523 | fitString.ReplaceAll("++","#"); | |
524 | fitString.ReplaceAll(" ",""); | |
525 | TObjArray *arrFitParams = fitString.Tokenize("#"); | |
526 | Int_t ndim = arrFitParams->GetEntries(); | |
5312f439 | 527 | //create array of TFormulas to evaluate the parameters |
528 | TObjArray *arrFitFormulas = new TObjArray(ndim); | |
529 | arrFitFormulas->SetOwner(kTRUE); | |
530 | for (Int_t idim=0;idim<ndim;++idim){ | |
531 | TString s=((TObjString*)arrFitParams->At(idim))->GetString(); | |
532 | s.ReplaceAll("gx","[0]"); | |
533 | s.ReplaceAll("gy","[1]"); | |
534 | s.ReplaceAll("lx","[2]"); | |
535 | s.ReplaceAll("ly","[3]"); | |
7390f655 | 536 | s.ReplaceAll("sector","[4]"); |
5312f439 | 537 | arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim); |
538 | } | |
732e90a8 | 539 | delete arrFitParams; |
540 | ||
541 | return arrFitFormulas; | |
542 | } | |
543 | //_____________________________________________________________________________ | |
544 | void AliTPCCalPad::EvalFormulaArray(const TObjArray &arrFitFormulas, TVectorD &results, | |
545 | const Int_t sec, const Int_t row, const Int_t pad) | |
546 | { | |
547 | // | |
548 | // evaluate the fit formulas | |
549 | // | |
550 | Int_t ndim=arrFitFormulas.GetEntries(); | |
551 | results.ResizeTo(ndim); | |
552 | ||
5312f439 | 553 | AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position |
554 | Float_t localXYZ[3]; | |
555 | Float_t globalXYZ[3]; | |
732e90a8 | 556 | tpcROCinstance->GetPositionLocal(sec, row, pad, localXYZ); |
557 | tpcROCinstance->GetPositionGlobal(sec, row, pad, globalXYZ); | |
558 | //calculate parameter values | |
559 | for (Int_t idim=0;idim<ndim;++idim){ | |
560 | TFormula *f=(TFormula*)arrFitFormulas.At(idim); | |
561 | f->SetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],sec); | |
562 | results[idim]=f->Eval(0); | |
563 | } | |
564 | } | |
565 | //_____________________________________________________________________________ | |
566 | void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, const char* fitFormula, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, AliTPCCalPad *pointError, Bool_t robust, Double_t robustFraction){ | |
567 | // | |
568 | // Performs a fit on both sides. | |
569 | // Valid information for the fitFormula are the variables | |
570 | // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName | |
571 | // - sector: the sector number. | |
572 | // eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on | |
573 | // | |
574 | // PadOutliers - pads with value !=0 are not used in fitting procedure | |
575 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
576 | // robustFraction: Fraction of data that will be used in EvalRobust | |
577 | // | |
578 | ||
579 | TObjArray* arrFitFormulas=CreateFormulaArray(fitFormula); | |
580 | Int_t ndim = arrFitFormulas->GetEntries(); | |
581 | //resize output data arrays | |
582 | fitParamSideA.ResizeTo(ndim+1); | |
583 | fitParamSideC.ResizeTo(ndim+1); | |
584 | covMatrixSideA.ResizeTo(ndim+1,ndim+1); | |
585 | covMatrixSideC.ResizeTo(ndim+1,ndim+1); | |
586 | // create linear fitter for A- and C- Side | |
587 | TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim)); | |
588 | TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim)); | |
589 | fitterGA->StoreData(kTRUE); | |
590 | fitterGC->StoreData(kTRUE); | |
591 | //parameter values | |
5312f439 | 592 | TVectorD parValues(ndim); |
732e90a8 | 593 | |
594 | AliTPCCalROC *rocErr=0x0; | |
5312f439 | 595 | |
596 | for (UInt_t isec = 0; isec<kNsec; ++isec){ | |
597 | AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec); | |
598 | AliTPCCalROC *rocData=GetCalROC(isec); | |
732e90a8 | 599 | if (pointError) rocErr=pointError->GetCalROC(isec); |
5312f439 | 600 | if (!rocData) continue; |
601 | for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { | |
602 | for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { | |
603 | //check for outliers | |
604 | if (rocOut && rocOut->GetValue(irow,ipad)) continue; | |
5312f439 | 605 | //calculate parameter values |
732e90a8 | 606 | EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad); |
5312f439 | 607 | //get value |
608 | Float_t value=rocData->GetValue(irow,ipad); | |
732e90a8 | 609 | //point error |
610 | Int_t err=1; | |
611 | if (rocErr) { | |
c018bcd4 | 612 | err=TMath::Nint(rocErr->GetValue(irow,ipad)); |
732e90a8 | 613 | if (err==0) err=1; |
614 | } | |
5312f439 | 615 | //add points to the fitters |
616 | if (isec/18%2==0){ | |
732e90a8 | 617 | fitterGA->AddPoint(parValues.GetMatrixArray(),value,err); |
5312f439 | 618 | }else{ |
732e90a8 | 619 | fitterGC->AddPoint(parValues.GetMatrixArray(),value,err); |
5312f439 | 620 | } |
621 | } | |
622 | } | |
623 | } | |
624 | if (robust){ | |
625 | fitterGA->EvalRobust(robustFraction); | |
626 | fitterGC->EvalRobust(robustFraction); | |
627 | } else { | |
628 | fitterGA->Eval(); | |
629 | fitterGC->Eval(); | |
630 | } | |
631 | chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1)); | |
632 | chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1)); | |
633 | fitterGA->GetParameters(fitParamSideA); | |
634 | fitterGC->GetParameters(fitParamSideC); | |
635 | fitterGA->GetCovarianceMatrix(covMatrixSideA); | |
636 | fitterGC->GetCovarianceMatrix(covMatrixSideC); | |
637 | ||
5312f439 | 638 | delete arrFitFormulas; |
639 | delete fitterGA; | |
640 | delete fitterGC; | |
641 | ||
642 | } | |
732e90a8 | 643 | // |
644 | AliTPCCalPad *AliTPCCalPad::CreateCalPadFit(const char* fitFormula, const TVectorD &fitParamSideA, const TVectorD &fitParamSideC) | |
645 | { | |
646 | // | |
647 | // | |
648 | // | |
649 | TObjArray *arrFitFormulas=CreateFormulaArray(fitFormula); | |
650 | Int_t ndim = arrFitFormulas->GetEntries(); | |
651 | //check if dimension of fit formula and fit parameters agree | |
652 | if (ndim!=fitParamSideA.GetNrows()||ndim!=fitParamSideC.GetNrows()){ | |
653 | printf("AliTPCCalPad::CreateCalPadFit: Dimensions of fit formula and fit Parameters does not match!"); | |
654 | return 0; | |
655 | } | |
656 | //create cal pad | |
657 | AliTPCCalPad *pad=new AliTPCCalPad("fitResultPad",Form("Fit result: %s",fitFormula)); | |
658 | //fill cal pad with fit results if requested | |
659 | for (UInt_t isec = 0; isec<kNsec; ++isec){ | |
660 | AliTPCCalROC *roc=pad->GetCalROC(isec); | |
661 | for (UInt_t irow = 0; irow < roc->GetNrows(); irow++) { | |
662 | for (UInt_t ipad = 0; ipad < roc->GetNPads(irow); ipad++) { | |
663 | const TVectorD *fitPar=0; | |
664 | TVectorD fitResArray; | |
665 | if (isec/18%2==0){ | |
666 | fitPar=&fitParamSideA; | |
667 | }else{ | |
668 | fitPar=&fitParamSideC; | |
669 | } | |
670 | EvalFormulaArray(*arrFitFormulas,fitResArray, isec, irow, ipad); | |
671 | for (Int_t idim=0;idim<ndim;++idim) | |
672 | fitResArray(idim)*=(*fitPar)(idim); | |
673 | roc->SetValue(irow,ipad,fitResArray.Sum()); | |
674 | } | |
675 | } | |
676 | } | |
677 | delete arrFitFormulas; | |
678 | return pad; | |
679 | } | |
5312f439 | 680 | /* |
a6d2bd0c | 681 | 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){ |
682 | // | |
683 | // Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side | |
684 | // fitType == 0: fit plane function | |
685 | // fitType == 1: fit parabolic function | |
686 | // PadOutliers - pads with value !=0 are not used in fitting procedure | |
687 | // chi2Threshold: Threshold for chi2 when EvalRobust is called | |
688 | // robustFraction: Fraction of data that will be used in EvalRobust | |
689 | // | |
690 | TLinearFitter* fitterGA = 0; | |
691 | TLinearFitter* fitterGC = 0; | |
692 | ||
693 | if (fitType == 1) { | |
694 | fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5"); | |
695 | fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5"); | |
696 | } | |
697 | else { | |
698 | fitterGA = new TLinearFitter(3,"x0++x1++x2"); | |
699 | fitterGC = new TLinearFitter(3,"x0++x1++x2"); | |
700 | } | |
701 | fitterGA->StoreData(kTRUE); | |
702 | fitterGC->StoreData(kTRUE); | |
703 | fitterGA->ClearPoints(); | |
704 | fitterGC->ClearPoints(); | |
705 | Double_t xx[6]; | |
706 | Int_t npointsA=0; | |
707 | Int_t npointsC=0; | |
708 | ||
709 | Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position | |
710 | Float_t lx, ly; // pads position | |
711 | ||
712 | AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position | |
713 | ||
714 | // loop over all sectors and pads and read data into fitterGA and fitterGC | |
715 | if (fitType == 1) { | |
716 | // parabolic fit | |
717 | fitParamSideA.ResizeTo(6); | |
718 | fitParamSideC.ResizeTo(6); | |
719 | covMatrixSideA.ResizeTo(6,6); | |
720 | covMatrixSideC.ResizeTo(6,6); | |
721 | for (UInt_t isec = 0; isec<72; isec++){ | |
722 | for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { | |
723 | for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { | |
724 | // fill fitterG | |
725 | tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number | |
726 | lx = localXY[0]; | |
727 | ly = localXY[1]; | |
728 | xx[0] = 1; | |
729 | xx[1] = lx; | |
730 | xx[2] = ly; | |
731 | xx[3] = lx*lx; | |
732 | xx[4] = ly*ly; | |
733 | xx[5] = lx*ly; | |
734 | if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) { | |
735 | // if given pad is no outlier, add it to TLinearFitter, decide to which of both | |
736 | // sector 0 - 17: IROC, A | |
737 | // sector 18 - 35: IROC, C | |
738 | // sector 36 - 53: OROC, A | |
739 | // sector 54 - 71: CROC, C | |
740 | if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A | |
741 | npointsA++; | |
742 | fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
743 | } | |
744 | else { // side C | |
745 | npointsC++; | |
746 | fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
747 | } | |
748 | } | |
749 | } | |
750 | } | |
751 | } | |
752 | } | |
753 | else { | |
754 | // linear fit | |
755 | fitParamSideA.ResizeTo(3); | |
756 | fitParamSideC.ResizeTo(3); | |
757 | covMatrixSideA.ResizeTo(3,3); | |
758 | covMatrixSideC.ResizeTo(3,3); | |
759 | ||
760 | for (UInt_t isec = 0; isec<72; isec++){ | |
761 | for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { | |
762 | for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { | |
763 | // fill fitterG | |
764 | tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number | |
765 | lx = localXY[0]; | |
766 | ly = localXY[1]; | |
767 | xx[0] = 1; | |
768 | xx[1] = lx; | |
769 | xx[2] = ly; | |
770 | if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) { | |
771 | // if given pad is no outlier, add it to TLinearFitter, decide to which of both | |
772 | // sector 0 - 17: IROC, A | |
773 | // sector 18 - 35: IROC, C | |
774 | // sector 36 - 53: OROC, A | |
775 | // sector 54 - 71: CROC, C | |
776 | if (isec <= 17 || (isec >= 36 && isec <= 53)) { | |
777 | // Side A | |
778 | npointsA++; | |
779 | fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
780 | } | |
781 | else { | |
782 | // side C | |
783 | npointsC++; | |
784 | fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); | |
785 | } | |
786 | } | |
787 | } | |
788 | } | |
789 | } | |
790 | } | |
791 | ||
792 | fitterGA->Eval(); | |
793 | fitterGC->Eval(); | |
794 | fitterGA->GetParameters(fitParamSideA); | |
795 | fitterGC->GetParameters(fitParamSideC); | |
796 | fitterGA->GetCovarianceMatrix(covMatrixSideA); | |
797 | fitterGC->GetCovarianceMatrix(covMatrixSideC); | |
798 | if (fitType == 1){ | |
799 | chi2SideA = fitterGA->GetChisquare()/(npointsA-6.); | |
800 | chi2SideC = fitterGC->GetChisquare()/(npointsC-6.); | |
801 | } | |
802 | else { | |
803 | chi2SideA = fitterGA->GetChisquare()/(npointsA-3.); | |
804 | chi2SideC = fitterGC->GetChisquare()/(npointsC-3.); | |
805 | } | |
806 | if (robust && chi2SideA > chi2Threshold) { | |
807 | // std::cout << "robust fitter called... " << std::endl; | |
808 | fitterGA->EvalRobust(robustFraction); | |
809 | fitterGA->GetParameters(fitParamSideA); | |
810 | } | |
811 | if (robust && chi2SideC > chi2Threshold) { | |
812 | // std::cout << "robust fitter called... " << std::endl; | |
813 | fitterGC->EvalRobust(robustFraction); | |
814 | fitterGC->GetParameters(fitParamSideC); | |
815 | } | |
816 | delete fitterGA; | |
817 | delete fitterGC; | |
818 | } | |
5312f439 | 819 | */ |
586007f3 | 820 |