TPC/AliTPCCalPad.cxx

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