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