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