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 }