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 AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
 388   //
 389   // MakeLocalFit - smoothing
 390   // rowRadius  -  radius - rows to be used for smoothing
 391   // padradius  -  radius - pads to be used for smoothing
 392   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 393   // robust     -  robust method of fitting  - (much slower)
 394
 395   AliTPCCalROC * ROCfitted = new AliTPCCalROC(fSector);
 396   TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
 397   fitterQ.StoreData(kTRUE);
 398   for (UInt_t row=0; row < GetNrows(); row++) {
 399     //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
 400     for (UInt_t pad=0; pad < GetNPads(row); pad++)
 401       ROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust));
 402   }
 403   return ROCfitted;
 404 }
 405
 406
 407 Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
 408   //
 409   //  AliTPCCalROC::GetNeighbourhoodValue - smoothing (PRIVATE)
 410   // rowRadius  -  radius - rows to be used for smoothing
 411   // padradius  -  radius - pads to be used for smoothing
 412   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 413   // robust     -  robust method of fitting  - (much slower)
 414
 415
 416
 417   fitterQ->ClearPoints();
 418   TVectorD fitParam(6);
 419   Int_t    npoints=0;
 420   Double_t xx[6];
 421   Float_t dlx, dly;
 422   Float_t lPad[3] = {0};
 423   Float_t localXY[3] = {0};
 424
 425   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 426   tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad);  // calculate position lPad by pad and row number
 427
 428   TArrayI *neighbourhoodRows = 0;
 429   TArrayI *neighbourhoodPads = 0;
 430   GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
 431
 432   Int_t r, p;
 433   for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
 434     r = neighbourhoodRows->At(i);
 435     p = neighbourhoodPads->At(i);
 436     if (r == -1 || p == -1) continue;
 437     tpcROCinstance->GetPositionLocal(fSector, r, p, localXY);   // calculate position localXY by pad and row number
 438     dlx = lPad[0] - localXY[0];
 439     dly = lPad[1] - localXY[1];
 440     xx[0] = 1;
 441     xx[1] = dlx;
 442     xx[2] = dly;
 443     xx[3] = dlx*dlx;
 444     xx[4] = dly*dly;
 445     xx[5] = dlx*dly;
 446     if (ROCoutliers && ROCoutliers->GetValue(r,p) != 1) {
 447       fitterQ->AddPoint(xx, GetValue(r, p), 1);
 448       npoints++;
 449     }
 450   }
 451   if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
 452     // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
 453     return 0.;  // for diagnostic
 454   }
 455   fitterQ->Eval();
 456   fitterQ->GetParameters(fitParam);
 457   Float_t chi2Q = 0;
 458   chi2Q = fitterQ->GetChisquare()/(npoints-6.);
 459   //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
 460   if (robust && chi2Q > 5) {
 461     //std::cout << "robust fitter called... " << std::endl;
 462     fitterQ->EvalRobust(0.7);
 463     fitterQ->GetParameters(fitParam);
 464   }
 465   Double_t value = fitParam[0];
 466
 467   delete neighbourhoodRows;
 468   delete neighbourhoodPads;
 469
 470   //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q <<  std::endl;
 471
 472   return value;
 473 }
 474
 475
 476
 477
 478 void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
 479   //
 480   //
 481   //
 482   rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
 483   padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
 484   Int_t* rowArrayTemp = rowArray->GetArray();
 485   Int_t* padArrayTemp = padArray->GetArray();
 486
 487   Int_t rmin = row - rRadius;
 488   UInt_t rmax = row + rRadius;
 489
 490   // if window goes out of ROC
 491   if (rmin < 0) {
 492     rmax = rmax - rmin;
 493     rmin = 0;
 494   }
 495   if (rmax >= GetNrows()) {
 496     rmin = rmin - (rmax - GetNrows()+1);
 497     rmax = GetNrows() - 1;
 498       if (rmin  < 0 ) rmin = 0; // if the window is bigger than the ROC
 499   }
 500
 501   Int_t pmin, pmax;
 502   Int_t i = 0;
 503
 504   for (UInt_t r = rmin; r <= rmax; r++) {
 505     pmin = pad - pRadius;
 506     pmax = pad + pRadius;
 507     if (pmin < 0) {
 508       pmax = pmax - pmin;
 509       pmin = 0;
 510     }
 511     if (pmax >= GetNPads(r)) {
 512       pmin = pmin - (pmax - GetNPads(r)+1);
 513       pmax = GetNPads(r) - 1;
 514       if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
 515     }
 516     for (Int_t p = pmin; p <= pmax; p++) {
 517       rowArrayTemp[i] = r;
 518       padArrayTemp[i] = p;
 519       i++;
 520     }
 521   }
 522   for (Int_t j = i; j < rowArray->GetSize(); j++){  // unused padArray-entries, in the case that the window is bigger than the ROC
 523     //std::cout << "trying to write -1" << std::endl;
 524     rowArrayTemp[j] = -1;
 525     padArrayTemp[j] = -1;
 526     //std::cout << "writing -1" << std::endl;
 527   }
 528 }
 529
 530
 531
 532 void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType){
 533   //
 534   // Makes global fit
 535   // do GlobalFit for given Secotr and return fit-parameters, covariance, and whatever
 536   // fitType == 0: fit plane
 537   // fitType == 1: fit parabolic
 538   // ROCoutliers - pads signed=1 - not used in fitting procedure
 539
 540   TLinearFitter* fitterG = 0;
 541   Double_t xx[6];
 542
 543   if (fitType  == 1)
 544     fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
 545   else
 546     fitterG = new TLinearFitter(3,"x0++x1++x2");
 547   fitterG->StoreData(kTRUE);
 548   fitterG->ClearPoints();
 549   Int_t    npoints=0;
 550
 551   Float_t dlx, dly;
 552   Float_t centerPad[3] = {0};
 553   Float_t localXY[3] = {0};
 554
 555   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 556   tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad);  // calculate center of ROC
 557
 558   // loop over all channels and read data into fitterG
 559   if (fitType == 1) {  // parabolic fit
 560     fitParam.ResizeTo(6);
 561     covMatrix.ResizeTo(6,6);
 562     for (UInt_t irow = 0; irow < GetNrows(); irow++) {
 563       for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
 564         // fill fitterG
 565         tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 566         dlx = centerPad[0] - localXY[0];
 567         dly = centerPad[1] - localXY[1];
 568         xx[0] = 1;
 569         xx[1] = dlx;
 570         xx[2] = dly;
 571         xx[3] = dlx*dlx;
 572         xx[4] = dly*dly;
 573         xx[5] = dlx*dly;
 574         if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) {
 575            npoints++;
 576            fitterG->AddPoint(xx, GetValue(irow, ipad), 1);
 577         }
 578       }
 579     }
 580   }
 581   else {   // linear fit
 582     fitParam.ResizeTo(3);
 583     covMatrix.ResizeTo(3,3);
 584     for (UInt_t irow = 0; irow < GetNrows(); irow++) {
 585       for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
 586         // fill fitterG
 587         tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 588         dlx = centerPad[0] - localXY[0];
 589         dly = centerPad[1] - localXY[1];
 590         xx[0] = 1;
 591         xx[1] = dlx;
 592         xx[2] = dly;
 593         if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) {
 594            npoints++;
 595            fitterG->AddPoint(xx, GetValue(irow, ipad), 1);
 596         }
 597       }
 598     }
 599   }
 600   fitterG->Eval();
 601   fitterG->GetParameters(fitParam);
 602   fitterG->GetCovarianceMatrix(covMatrix);
 603   if (fitType == 1)
 604     chi2 = fitterG->GetChisquare()/(npoints-6.);
 605   else chi2 = fitterG->GetChisquare()/(npoints-3.);
 606   if (robust && chi2 > 5) {
 607     //    std::cout << "robust fitter called... " << std::endl;
 608     fitterG->EvalRobust(0.7);
 609     fitterG->GetParameters(fitParam);
 610   }
 611   delete fitterG;
 612 }
 613
 614
 615 //
 616 AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector){
 617   //
 618   //
 619   // Create ROC with global fit parameters
 620   // fitType == 0: fit plane
 621   // fitType == 1: fit parabolic
 622   // loop over all channels and write fit values into ROCfitted
 623   //
 624   Float_t dlx, dly;
 625   Float_t centerPad[3] = {0};
 626   Float_t localXY[3] = {0};
 627   AliTPCCalROC * ROCfitted = new AliTPCCalROC(sector);
 628   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 629   tpcROCinstance->GetPositionLocal(sector, ROCfitted->GetNrows()/2, ROCfitted->GetNPads(ROCfitted->GetNrows()/2)/2, centerPad);  // calculate center of ROC
 630   Int_t fitType = 1;
 631   if (fitParam.GetNoElements() == 6) fitType = 1;
 632   else fitType = 0;
 633   Double_t value = 0;
 634   if (fitType == 1) { // parabolic fit
 635     for (UInt_t irow = 0; irow < ROCfitted->GetNrows(); irow++) {
 636       for (UInt_t ipad = 0; ipad < ROCfitted->GetNPads(irow); ipad++) {
 637         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 638         dlx = centerPad[0] - localXY[0];
 639         dly = centerPad[1] - localXY[1];
 640         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
 641         ROCfitted->SetValue(irow, ipad, value);
 642       }
 643     }
 644   }
 645   else {  // linear fit
 646     for (UInt_t irow = 0; irow < ROCfitted->GetNrows(); irow++) {
 647       for (UInt_t ipad = 0; ipad < ROCfitted->GetNPads(irow); ipad++) {
 648         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 649         dlx = centerPad[0] - localXY[0];
 650         dly = centerPad[1] - localXY[1];
 651         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
 652         ROCfitted->SetValue(irow, ipad, value);
 653       }
 654     }
 655   }
 656   return ROCfitted;
 657 }
 658
 659