Using the BxByBz tracking functions in AliReconstruction, and the corresponding chang...
[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.
7390f655 516 // Valid information for the fitFormula are the variables
517 // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName
518 // - sector: the sector number.
519 // eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on
5312f439 520 //
521 // PadOutliers - pads with value !=0 are not used in fitting procedure
522 // chi2Threshold: Threshold for chi2 when EvalRobust is called
523 // robustFraction: Fraction of data that will be used in EvalRobust
524 //
525
526 // split fit string in single parameters
527 // find dimension of the fit:
528 TString fitString(fitFormula);
529 fitString.ReplaceAll("++","#");
530 fitString.ReplaceAll(" ","");
531 TObjArray *arrFitParams = fitString.Tokenize("#");
532 Int_t ndim = arrFitParams->GetEntries();
533 //resize output data arrays
534 fitParamSideA.ResizeTo(ndim+1);
535 fitParamSideC.ResizeTo(ndim+1);
536 covMatrixSideA.ResizeTo(ndim+1,ndim+1);
537 covMatrixSideC.ResizeTo(ndim+1,ndim+1);
538 // create linear fitter for A- and C- Side
539 TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
540 TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
541 fitterGA->StoreData(kTRUE);
542 fitterGC->StoreData(kTRUE);
543 //create array of TFormulas to evaluate the parameters
544 TObjArray *arrFitFormulas = new TObjArray(ndim);
545 arrFitFormulas->SetOwner(kTRUE);
546 for (Int_t idim=0;idim<ndim;++idim){
547 TString s=((TObjString*)arrFitParams->At(idim))->GetString();
548 s.ReplaceAll("gx","[0]");
549 s.ReplaceAll("gy","[1]");
550 s.ReplaceAll("lx","[2]");
551 s.ReplaceAll("ly","[3]");
7390f655 552 s.ReplaceAll("sector","[4]");
5312f439 553 arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim);
554 }
555 //loop over data and add points to the fitter
556 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
557 Float_t localXYZ[3];
558 Float_t globalXYZ[3];
559 TVectorD parValues(ndim);
560
561 for (UInt_t isec = 0; isec<kNsec; ++isec){
562 AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec);
563 AliTPCCalROC *rocData=GetCalROC(isec);
564 if (!rocData) continue;
565 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
566 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
567 //check for outliers
568 if (rocOut && rocOut->GetValue(irow,ipad)) continue;
569 //calculate local and global pad positions
570 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXYZ);
571 tpcROCinstance->GetPositionGlobal(isec, irow, ipad, globalXYZ);
572 //calculate parameter values
573 for (Int_t idim=0;idim<ndim;++idim){
574 TFormula *f=(TFormula*)arrFitFormulas->At(idim);
7390f655 575 f->SetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],isec);
5312f439 576 parValues[idim]=f->Eval(0);
577 }
578 //get value
579 Float_t value=rocData->GetValue(irow,ipad);
580 //add points to the fitters
581 if (isec/18%2==0){
582 fitterGA->AddPoint(parValues.GetMatrixArray(),value,pointError);
583 }else{
584 fitterGC->AddPoint(parValues.GetMatrixArray(),value,pointError);
585 }
586 }
587 }
588 }
589 if (robust){
590 fitterGA->EvalRobust(robustFraction);
591 fitterGC->EvalRobust(robustFraction);
592 } else {
593 fitterGA->Eval();
594 fitterGC->Eval();
595 }
596 chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1));
597 chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1));
598 fitterGA->GetParameters(fitParamSideA);
599 fitterGC->GetParameters(fitParamSideC);
600 fitterGA->GetCovarianceMatrix(covMatrixSideA);
601 fitterGC->GetCovarianceMatrix(covMatrixSideC);
602
603 delete arrFitParams;
604 delete arrFitFormulas;
605 delete fitterGA;
606 delete fitterGC;
607
608}
a6d2bd0c 609
5312f439 610/*
a6d2bd0c 611void 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){
612 //
613 // Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side
614 // fitType == 0: fit plane function
615 // fitType == 1: fit parabolic function
616 // PadOutliers - pads with value !=0 are not used in fitting procedure
617 // chi2Threshold: Threshold for chi2 when EvalRobust is called
618 // robustFraction: Fraction of data that will be used in EvalRobust
619 //
620 TLinearFitter* fitterGA = 0;
621 TLinearFitter* fitterGC = 0;
622
623 if (fitType == 1) {
624 fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
625 fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
626 }
627 else {
628 fitterGA = new TLinearFitter(3,"x0++x1++x2");
629 fitterGC = new TLinearFitter(3,"x0++x1++x2");
630 }
631 fitterGA->StoreData(kTRUE);
632 fitterGC->StoreData(kTRUE);
633 fitterGA->ClearPoints();
634 fitterGC->ClearPoints();
635 Double_t xx[6];
636 Int_t npointsA=0;
637 Int_t npointsC=0;
638
639 Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position
640 Float_t lx, ly; // pads position
641
642 AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
643
644 // loop over all sectors and pads and read data into fitterGA and fitterGC
645 if (fitType == 1) {
646 // parabolic fit
647 fitParamSideA.ResizeTo(6);
648 fitParamSideC.ResizeTo(6);
649 covMatrixSideA.ResizeTo(6,6);
650 covMatrixSideC.ResizeTo(6,6);
651 for (UInt_t isec = 0; isec<72; isec++){
652 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
653 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
654 // fill fitterG
655 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
656 lx = localXY[0];
657 ly = localXY[1];
658 xx[0] = 1;
659 xx[1] = lx;
660 xx[2] = ly;
661 xx[3] = lx*lx;
662 xx[4] = ly*ly;
663 xx[5] = lx*ly;
664 if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
665 // if given pad is no outlier, add it to TLinearFitter, decide to which of both
666// sector 0 - 17: IROC, A
667// sector 18 - 35: IROC, C
668// sector 36 - 53: OROC, A
669// sector 54 - 71: CROC, C
670 if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A
671 npointsA++;
672 fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
673 }
674 else { // side C
675 npointsC++;
676 fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
677 }
678 }
679 }
680 }
681 }
682 }
683 else {
684 // linear fit
685 fitParamSideA.ResizeTo(3);
686 fitParamSideC.ResizeTo(3);
687 covMatrixSideA.ResizeTo(3,3);
688 covMatrixSideC.ResizeTo(3,3);
689
690 for (UInt_t isec = 0; isec<72; isec++){
691 for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
692 for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
693 // fill fitterG
694 tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
695 lx = localXY[0];
696 ly = localXY[1];
697 xx[0] = 1;
698 xx[1] = lx;
699 xx[2] = ly;
700 if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
701 // if given pad is no outlier, add it to TLinearFitter, decide to which of both
702// sector 0 - 17: IROC, A
703// sector 18 - 35: IROC, C
704// sector 36 - 53: OROC, A
705// sector 54 - 71: CROC, C
706 if (isec <= 17 || (isec >= 36 && isec <= 53)) {
707 // Side A
708 npointsA++;
709 fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
710 }
711 else {
712 // side C
713 npointsC++;
714 fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
715 }
716 }
717 }
718 }
719 }
720 }
721
722 fitterGA->Eval();
723 fitterGC->Eval();
724 fitterGA->GetParameters(fitParamSideA);
725 fitterGC->GetParameters(fitParamSideC);
726 fitterGA->GetCovarianceMatrix(covMatrixSideA);
727 fitterGC->GetCovarianceMatrix(covMatrixSideC);
728 if (fitType == 1){
729 chi2SideA = fitterGA->GetChisquare()/(npointsA-6.);
730 chi2SideC = fitterGC->GetChisquare()/(npointsC-6.);
731 }
732 else {
733 chi2SideA = fitterGA->GetChisquare()/(npointsA-3.);
734 chi2SideC = fitterGC->GetChisquare()/(npointsC-3.);
735 }
736 if (robust && chi2SideA > chi2Threshold) {
737 // std::cout << "robust fitter called... " << std::endl;
738 fitterGA->EvalRobust(robustFraction);
739 fitterGA->GetParameters(fitParamSideA);
740 }
741 if (robust && chi2SideC > chi2Threshold) {
742 // std::cout << "robust fitter called... " << std::endl;
743 fitterGC->EvalRobust(robustFraction);
744 fitterGC->GetParameters(fitParamSideC);
745 }
746 delete fitterGA;
747 delete fitterGC;
748}
5312f439 749*/
586007f3 750