Add global fit of values (Marian, S.Gaertner, L.Bozyk)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalROC.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
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  Calibration base class for a single ROC                                  //
20 //  Contains one float value per pad                                         //
21 //     mapping of the pads taken form AliTPCROC                              //
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 //
26 // ROOT includes 
27 //
28 #include "TMath.h"
29 #include "TClass.h"
30 #include "TFile.h"
31 #include "TH1F.h"
32 #include "TH2F.h"
33 #include "TArrayI.h"
34 //
35 //
36 #include "AliTPCCalROC.h"
37 #include "AliMathBase.h"
38
39 ClassImp(AliTPCCalROC)
40
41
42 //_____________________________________________________________________________
43 AliTPCCalROC::AliTPCCalROC()
44              :TObject(),
45               fSector(0),
46               fNChannels(0),
47               fNRows(0),
48               fIndexes(0),
49               fData(0)
50 {
51   //
52   // Default constructor
53   //
54
55 }
56
57 //_____________________________________________________________________________
58 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
59              :TObject(),
60               fSector(0),
61               fNChannels(0),
62               fNRows(0),
63               fIndexes(0),
64               fData(0)
65 {
66   //
67   // Constructor that initializes a given sector
68   //
69   fSector = sector;
70   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
71   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
72   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
73   fData = new Float_t[fNChannels];
74   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
75 }
76
77 //_____________________________________________________________________________
78 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
79              :TObject(c),
80               fSector(0),
81               fNChannels(0),
82               fNRows(0),
83               fIndexes(0),
84               fData(0)
85 {
86   //
87   // AliTPCCalROC copy constructor
88   //
89   fSector = c.fSector;
90   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
91   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
92   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
93   //
94   fData   = new Float_t[fNChannels];
95   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
96 }
97 //____________________________________________________________________________
98 AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
99 {
100   //
101   // assignment operator - dummy
102   //
103   fData=param.fData;
104   return (*this);
105 }
106
107
108 //_____________________________________________________________________________
109 AliTPCCalROC::~AliTPCCalROC()
110 {
111   //
112   // AliTPCCalROC destructor
113   //
114   if (fData) {
115     delete [] fData;
116     fData = 0;
117   }
118 }
119
120
121
122 void AliTPCCalROC::Streamer(TBuffer &R__b)
123 {
124    // Stream an object of class AliTPCCalROC.
125    if (R__b.IsReading()) {
126       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
127       fIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
128    } else {
129       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
130    }
131 }
132
133 //  //
134 //   // algebra
135 //   void Add(Float_t c1);
136 //   void Multiply(Float_t c1);
137 //   void Add(const AliTPCCalROC * roc, Double_t c1 = 1);
138 //   void Divide(const AliTPCCalROC * roc);   
139
140 void AliTPCCalROC::Add(Float_t c1){
141   //
142   // add constant
143   //
144   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
145 }
146 void AliTPCCalROC::Multiply(Float_t c1){
147   //
148   // add constant
149   //
150   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
151 }
152
153 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
154   //
155   // add values 
156   //
157   for (UInt_t  idata = 0; idata< fNChannels; idata++){
158     fData[idata]+=roc->fData[idata]*c1;
159   }
160 }
161
162
163 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
164   //
165   // multiply values - per by pad
166   //
167   for (UInt_t  idata = 0; idata< fNChannels; idata++){
168     fData[idata]*=roc->fData[idata];
169   }
170 }
171
172
173 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
174   //
175   // divide values 
176   //
177   Float_t kEpsilon=0.00000000000000001;
178   for (UInt_t  idata = 0; idata< fNChannels; idata++){
179     if (TMath::Abs(roc->fData[idata])>kEpsilon)
180       fData[idata]/=roc->fData[idata];
181   }
182 }
183
184
185
186
187 Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
188   //
189   //  Calculate LTM mean and sigma
190   //
191   Double_t *ddata = new Double_t[fNChannels];
192   Double_t mean=0, lsigma=0;
193   Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels));
194   for (UInt_t i=0;i<fNChannels;i++) ddata[i]= fData[i];
195   AliMathBase::EvaluateUni(UInt_t(fNChannels),ddata, mean, lsigma, hh);
196   if (sigma) *sigma=lsigma;
197   delete [] ddata;
198   return mean;
199 }
200
201 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
202   //
203   // make 1D histo
204   // type -1 = user defined range
205   //       0 = nsigma cut nsigma=min
206   //       1 = delta cut around median delta=min
207   if (type>=0){
208     if (type==0){
209       // nsigma range
210       Float_t mean  = GetMean();
211       Float_t sigma = GetRMS();
212       Float_t nsigma = TMath::Abs(min);
213       min = mean-nsigma*sigma;
214       max = mean+nsigma*sigma;
215     }
216     if (type==1){
217       // fixed range
218       Float_t mean   = GetMedian();
219       Float_t  delta = min;
220       min = mean-delta;
221       max = mean+delta;
222     }
223     if (type==2){
224       //
225       // LTM mean +- nsigma
226       //
227       Double_t sigma;
228       Float_t mean  = GetLTM(&sigma,max);
229       sigma*=min;
230       min = mean-sigma;
231       max = mean+sigma;
232     }
233   }
234   char  name[1000];
235   sprintf(name,"%s ROC 1D%d",GetTitle(),fSector);
236   TH1F * his = new TH1F(name,name,100, min,max);
237   for (UInt_t irow=0; irow<fNRows; irow++){
238     UInt_t npads = (Int_t)GetNPads(irow);
239     for (UInt_t ipad=0; ipad<npads; ipad++){
240       his->Fill(GetValue(irow,ipad));
241     }
242   }
243   return his;
244 }
245
246
247
248 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
249   //
250   // make 2D histo
251   // type -1 = user defined range
252   //       0 = nsigma cut nsigma=min
253   //       1 = delta cut around median delta=min
254   if (type>=0){
255     if (type==0){
256       // nsigma range
257       Float_t mean  = GetMean();
258       Float_t sigma = GetRMS();
259       Float_t nsigma = TMath::Abs(min);
260       min = mean-nsigma*sigma;
261       max = mean+nsigma*sigma;
262     }
263     if (type==1){
264       // fixed range
265       Float_t mean   = GetMedian();
266       Float_t  delta = min;
267       min = mean-delta;
268       max = mean+delta;
269     }
270     if (type==2){
271       Double_t sigma;
272       Float_t mean  = GetLTM(&sigma,max);
273       sigma*=min;
274       min = mean-sigma;
275       max = mean+sigma;
276
277     }
278   }
279   UInt_t maxPad = 0;
280   for (UInt_t irow=0; irow<fNRows; irow++){
281     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
282   }
283   char  name[1000];
284   sprintf(name,"%s ROC%d",GetTitle(),fSector);
285   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
286   for (UInt_t irow=0; irow<fNRows; irow++){
287     UInt_t npads = (Int_t)GetNPads(irow);
288     for (UInt_t ipad=0; ipad<npads; ipad++){
289       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
290     }
291   }
292   his->SetMaximum(max);
293   his->SetMinimum(min);
294   return his;
295 }
296
297 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
298   //
299   // Make Histogram with outliers
300   // mode = 0 - sigma cut used
301   // mode = 1 - absolute cut used
302   // fraction - fraction of values used to define sigma
303   // delta - in mode 0 - nsigma cut
304   //            mode 1 - delta cut
305   Double_t sigma;
306   Float_t mean  = GetLTM(&sigma,fraction);  
307   if (type==0) delta*=sigma; 
308   UInt_t maxPad = 0;
309   for (UInt_t irow=0; irow<fNRows; irow++){
310     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
311   }
312
313   char  name[1000];
314   sprintf(name,"%s ROC Outliers%d",GetTitle(),fSector);
315   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
316   for (UInt_t irow=0; irow<fNRows; irow++){
317     UInt_t npads = (Int_t)GetNPads(irow);
318     for (UInt_t ipad=0; ipad<npads; ipad++){
319       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
320         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
321     }
322   }
323   return his;
324 }
325
326
327
328 void AliTPCCalROC::Draw(Option_t* opt){
329   //
330   // create histogram with values and draw it
331   //
332   TH1 * his=0; 
333   TString option=opt;
334   option.ToUpper();
335   if (option.Contains("1D")){
336     his = MakeHisto1D();
337   }
338   else{
339     his = MakeHisto2D();
340   }
341   his->Draw(option);
342 }
343
344
345
346
347
348 void AliTPCCalROC::Test(){
349   //
350   // example function to show functionality and tes AliTPCCalROC
351   //
352   AliTPCCalROC  roc0(0);  
353   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
354     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
355       Float_t value  = irow+ipad/1000.;
356       roc0.SetValue(irow,ipad,value);
357     }
358   }
359   //
360   AliTPCCalROC roc1(roc0);
361   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
362     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
363       Float_t value  = irow+ipad/1000.;
364       if (roc1.GetValue(irow,ipad)!=value){
365         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
366       }
367     }
368   }  
369   TFile f("calcTest.root","recreate");
370   roc0.Write("Roc0");
371   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
372   f.Close();
373   //
374   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
375     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
376       printf("NPads - Read/Write error\trow=%d\n",irow);
377     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
378       Float_t value  = irow+ipad/1000.;
379       if (roc2->GetValue(irow,ipad)!=value){
380         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
381       }
382     }
383   }   
384 }
385
386
387
388 AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
389   //
390   // MakeLocalFit - smoothing
391   // rowRadius  -  radius - rows to be used for smoothing
392   // padradius  -  radius - pads to be used for smoothing
393   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
394   // robust     -  robust method of fitting  - (much slower)
395   
396   AliTPCCalROC * ROCfitted = new AliTPCCalROC(fSector);
397   TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
398   fitterQ.StoreData(kTRUE);
399   for (UInt_t row=0; row < GetNrows(); row++) {
400     //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
401     for (UInt_t pad=0; pad < GetNPads(row); pad++)
402       ROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust));
403   }
404   return ROCfitted;
405 }
406
407
408 Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
409   //
410   //  AliTPCCalROC::GetNeighbourhoodValue - smoothing (PRIVATE)
411   // rowRadius  -  radius - rows to be used for smoothing
412   // padradius  -  radius - pads to be used for smoothing
413   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
414   // robust     -  robust method of fitting  - (much slower)
415   
416
417
418   fitterQ->ClearPoints();
419   TVectorD fitParam(6);
420   Int_t    npoints=0;
421   Double_t xx[6];
422   Float_t dlx, dly;
423   Float_t lPad[3] = {0};
424   Float_t localXY[3] = {0};
425   
426   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
427   tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad);  // calculate position lPad by pad and row number
428   
429   TArrayI *neighbourhoodRows = 0;
430   TArrayI *neighbourhoodPads = 0;
431   GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
432   
433   Int_t r, p;
434   for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
435     r = neighbourhoodRows->At(i);
436     p = neighbourhoodPads->At(i);
437     if (r == -1 || p == -1) continue;
438     tpcROCinstance->GetPositionLocal(fSector, r, p, localXY);   // calculate position localXY by pad and row number
439     dlx = lPad[0] - localXY[0];
440     dly = lPad[1] - localXY[1];
441     xx[0] = 1;
442     xx[1] = dlx;
443     xx[2] = dly;
444     xx[3] = dlx*dlx;
445     xx[4] = dly*dly;
446     xx[5] = dlx*dly;
447     if (ROCoutliers && ROCoutliers->GetValue(r,p) != 1) {
448       fitterQ->AddPoint(xx, GetValue(r, p), 1);
449       npoints++;
450     }
451   }
452   if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
453     // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
454     return 0.;  // for diagnostic
455   }
456   fitterQ->Eval();
457   fitterQ->GetParameters(fitParam);
458   Float_t chi2Q = 0;
459   chi2Q = fitterQ->GetChisquare()/(npoints-6.);
460   //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
461   if (robust && chi2Q > 5) {
462     //std::cout << "robust fitter called... " << std::endl;
463     fitterQ->EvalRobust(0.7);
464     fitterQ->GetParameters(fitParam);
465   }
466   Double_t value = fitParam[0];
467   
468   delete neighbourhoodRows;
469   delete neighbourhoodPads;
470   
471   //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q <<  std::endl;
472   
473   return value;
474 }
475
476
477
478
479 void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
480   //
481   //
482   //
483   rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
484   padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
485   Int_t* rowArrayTemp = rowArray->GetArray();
486   Int_t* padArrayTemp = padArray->GetArray();
487   
488   Int_t rmin = row - rRadius;
489   UInt_t rmax = row + rRadius;
490   
491   // if window goes out of ROC
492   if (rmin < 0) {
493     rmax = rmax - rmin;
494     rmin = 0;
495   }
496   if (rmax >= GetNrows()) {
497     rmin = rmin - (rmax - GetNrows()+1);
498     rmax = GetNrows() - 1;
499       if (rmin  < 0 ) rmin = 0; // if the window is bigger than the ROC
500   }
501   
502   Int_t pmin, pmax;
503   Int_t i = 0;
504   
505   for (UInt_t r = rmin; r <= rmax; r++) {
506     pmin = pad - pRadius;
507     pmax = pad + pRadius;
508     if (pmin < 0) {
509       pmax = pmax - pmin;
510       pmin = 0;
511     }
512     if (pmax >= GetNPads(r)) {
513       pmin = pmin - (pmax - GetNPads(r)+1);
514       pmax = GetNPads(r) - 1;
515       if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
516     }
517     for (Int_t p = pmin; p <= pmax; p++) {
518       rowArrayTemp[i] = r;
519       padArrayTemp[i] = p;
520       i++;
521     }
522   }
523   for (Int_t j = i; j < rowArray->GetSize(); j++){  // unused padArray-entries, in the case that the window is bigger than the ROC
524     //std::cout << "trying to write -1" << std::endl;
525     rowArrayTemp[j] = -1;
526     padArrayTemp[j] = -1;
527     //std::cout << "writing -1" << std::endl;
528   } 
529 }
530
531
532
533 void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType){
534   //
535   // Makes global fit  
536   // do GlobalFit for given Secotr and return fit-parameters, covariance, and whatever
537   // fitType == 0: fit plane
538   // fitType == 1: fit parabolic
539   // ROCoutliers - pads signed=1 - not used in fitting procedure
540    
541   TLinearFitter* fitterG = 0;
542   Double_t xx[6];
543   
544   if (fitType  == 1) 
545     fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
546   else 
547     fitterG = new TLinearFitter(3,"x0++x1++x2");
548   //fitterG->StoreData(kTRUE);   
549   fitterG->ClearPoints();
550   Int_t    npoints=0;
551   
552   Float_t dlx, dly;
553   Float_t centerPad[3] = {0};
554   Float_t localXY[3] = {0};
555   
556   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
557   tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad);  // calculate center of ROC 
558   
559   // loop over all channels and read data into fitterG
560   if (fitType == 1) {  // parabolic fit
561     fitParam.ResizeTo(6);
562     covMatrix.ResizeTo(6,6);
563     for (UInt_t irow = 0; irow < GetNrows(); irow++) {
564       for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
565         // fill fitterG
566         tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
567         dlx = centerPad[0] - localXY[0];
568         dly = centerPad[1] - localXY[1];
569         xx[0] = 1;
570         xx[1] = dlx;
571         xx[2] = dly;
572         xx[3] = dlx*dlx;
573         xx[4] = dly*dly;
574         xx[5] = dlx*dly;
575         if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) 
576           fitterG->AddPoint(xx, GetValue(irow, ipad), 1);  
577       }
578     }
579   }
580   else {   // linear fit
581     fitParam.ResizeTo(3);
582     covMatrix.ResizeTo(3,3);
583     for (UInt_t irow = 0; irow < GetNrows(); irow++) {
584       for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
585         // fill fitterG
586         tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
587         dlx = centerPad[0] - localXY[0];
588         dly = centerPad[1] - localXY[1];
589         xx[0] = 1;
590         xx[1] = dlx;
591         xx[2] = dly;
592         if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) 
593           fitterG->AddPoint(xx, GetValue(irow, ipad), 1);  
594       }
595     }
596   }
597   fitterG->Eval();
598   fitterG->GetParameters(fitParam);
599   fitterG->GetCovarianceMatrix(covMatrix);
600   if (fitType == 1)
601     chi2 = fitterG->GetChisquare()/(npoints-6.);
602   else chi2 = fitterG->GetChisquare()/(npoints-3.);
603   if (robust && chi2 > 5) {
604     //    std::cout << "robust fitter called... " << std::endl;
605     fitterG->EvalRobust(0.7);
606     fitterG->GetParameters(fitParam);
607   }
608   delete fitterG;
609 }
610
611
612 //
613 AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector)
614 {
615   //
616   //
617   // Create ROC with global fit parameters
618   // fitType == 0: fit plane
619   // fitType == 1: fit parabolic
620   // loop over all channels and write fit values into ROCfitted
621   //
622   Float_t dlx, dly;
623   Float_t centerPad[3] = {0};
624   Float_t localXY[3] = {0};
625   AliTPCCalROC * ROCfitted = new AliTPCCalROC(sector);
626   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
627   Int_t fitType = 1;
628   if (fitParam.GetNoElements() == 6) fitType = 1;
629   else fitType = 0;
630   Double_t value = 0;
631   if (fitType == 1) { // parabolic fit
632     for (UInt_t irow = 0; irow < ROCfitted->GetNrows(); irow++) {
633       for (UInt_t ipad = 0; ipad < ROCfitted->GetNPads(irow); ipad++) {
634         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
635         dlx = centerPad[0] - localXY[0];
636         dly = centerPad[1] - localXY[1];
637         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
638         ROCfitted->SetValue(irow, ipad, value);
639       }
640     }   
641   }
642   else {  // linear fit
643     for (UInt_t irow = 0; irow < ROCfitted->GetNrows(); irow++) {
644       for (UInt_t ipad = 0; ipad < ROCfitted->GetNPads(irow); ipad++) {
645         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
646         dlx = centerPad[0] - localXY[0];
647         dly = centerPad[1] - localXY[1];
648         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
649         ROCfitted->SetValue(irow, ipad, value);
650       }
651     }   
652   }
653   return ROCfitted;
654 }
655
656
657