AliTPCcalibDButil.cxx.diff Add Time0 creation using a combination of CE and...
[u/mrichter/AliRoot.git] / TPC / AliTPCCalPad.cxx
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 //                                                                           //
20 //  TPC calibration class for parameters which are saved per pad             //
21 //  Each AliTPCCalPad consists of 72 AliTPCCalROC-objects                    //
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 #include "AliTPCCalPad.h"
26 #include "AliTPCCalROC.h"
27 #include <TObjArray.h>
28 #include <TAxis.h>
29 #include <TGraph.h>
30 #include <TGraph2D.h>
31 #include <TH2F.h>
32 #include "TTreeStream.h"
33 #include "TFile.h"
34 #include "TKey.h"
35 #include <TFormula.h>
36 #include <TString.h>
37 #include <TObjString.h>
38 #include <iostream>
39 #include <AliLog.h>
40
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
76   for (Int_t isec = 0; isec < kNsec; isec++) {
77          fROC[isec] = 0;
78      if (c.fROC[isec])
79        fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
80   }
81 }
82
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
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 }
139
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();
148    if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector);
149    for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++) 
150       fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel));
151 }
152
153
154
155 //_____________________________________________________________________________
156 void AliTPCCalPad::Add(Float_t c1)
157 {
158     //
159     // add constant c1 to all channels of all ROCs
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 {
172   //
173   // multiply each channel of all ROCs with c1
174   //    
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 {
185   //
186   // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
187   //  - pad by pad -
188   //
189     for (Int_t isec = 0; isec < kNsec; isec++) {
190         if (fROC[isec] && pad->GetCalROC(isec)){
191             fROC[isec]->Add(pad->GetCalROC(isec),c1);
192         }
193     }
194 }
195
196 //_____________________________________________________________________________
197 void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
198 {
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++) {
204         if (fROC[isec]){
205             fROC[isec]->Multiply(pad->GetCalROC(isec));
206         }
207     }
208 }
209
210 //_____________________________________________________________________________
211 void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
212 {
213   //
214   // divide each channel of all ROCs by the coresponding channel of 'pad'
215   //     - pad by pad -
216   //    
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 //_____________________________________________________________________________
225 TGraph  *  AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
226   //
227   //   type=1 - mean
228   //        2 - median
229   //        3 - LTM
230   //
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 }
259
260 //_____________________________________________________________________________
261 Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
262 {
263     //
264     // Calculates mean and RMS of all ROCs
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 //_____________________________________________________________________________
289 Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
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++) {
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        }
304     }
305     return TMath::Mean(n,arr);
306 }
307
308 //_____________________________________________________________________________
309 Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
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++) {
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        }
324     }
325     return TMath::Mean(n,arr);
326 }
327
328 //_____________________________________________________________________________
329 Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
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++) {
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        }
344     }
345     return TMath::Mean(n,arr);
346 }
347
348 //_____________________________________________________________________________
349 Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
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;
363        AliTPCCalROC* outlierROC = 0;
364        if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
365             arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
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
379   //
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 //_____________________________________________________________________________
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
429   //
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
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 {
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++){
474       if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
475       if (PadOutliers)
476          pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
477       else 
478          pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
479    }
480    return pad;
481 }
482
483
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){
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
492    // err: error of the data points
493    // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array
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)
501          GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
502       else 
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);
510    }
511    return pad;
512 }
513 //_____________________________________________________________________________
514 TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula)
515 {
516   //
517   // create an array of TFormulas for the each parameter of the fit function
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();
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]");
536     s.ReplaceAll("sector","[4]");
537     arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim);
538   }
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   
553   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();  // to calculate the pad's position
554   Float_t localXYZ[3];
555   Float_t globalXYZ[3];
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
592   TVectorD parValues(ndim);
593
594   AliTPCCalROC *rocErr=0x0;
595   
596   for (UInt_t isec = 0; isec<kNsec; ++isec){
597     AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec);
598     AliTPCCalROC *rocData=GetCalROC(isec);
599     if (pointError) rocErr=pointError->GetCalROC(isec);
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;
605         //calculate parameter values
606         EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad);
607         //get value
608         Float_t value=rocData->GetValue(irow,ipad);
609         //point error
610         Int_t err=1;
611         if (rocErr) {
612           err=rocErr->GetValue(irow,ipad);
613           if (err==0) err=1;
614         }
615         //add points to the fitters
616         if (isec/18%2==0){
617           fitterGA->AddPoint(parValues.GetMatrixArray(),value,err);
618         }else{
619           fitterGC->AddPoint(parValues.GetMatrixArray(),value,err);
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   
638   delete arrFitFormulas;
639   delete fitterGA;
640   delete fitterGC;
641   
642 }
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 }
680 /*
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 }
819 */
820