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