TPC/Base/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 #include "AliTPCPreprocessorOnline.h"
  49 #include "AliTPCCalibViewer.h"
  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(array->GetName(),array->GetName())
  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 Bool_t  AliTPCCalPad::MedianFilter(Int_t deltaRow, Int_t deltaPad, AliTPCCalPad*outlierPad,  Bool_t doEdge){
 163   //
 164   // replace constent with median in the neigborhood
 165   //
 166   Bool_t isOK=kTRUE;
 167   for (Int_t isec = 0; isec < kNsec; isec++) {
 168     AliTPCCalROC *outlierROC=(outlierPad==NULL)?NULL:outlierPad->GetCalROC(isec);
 169     if (fROC[isec]){
 170       isOK&=fROC[isec]->MedianFilter(deltaRow,deltaPad,outlierROC,doEdge);
 171     }
 172   }
 173   return isOK;
 174 }
 175
 176 Bool_t  AliTPCCalPad::LTMFilter(Int_t deltaRow, Int_t deltaPad, Float_t fraction, Int_t type, AliTPCCalPad*outlierPad,  Bool_t doEdge){
 177   //
 178   // replace constent with LTM statistic  in  neigborhood
 179   //
 180   Bool_t isOK=kTRUE;
 181   for (Int_t isec = 0; isec < kNsec; isec++) {
 182     AliTPCCalROC *outlierROC=(outlierPad==NULL)?NULL:outlierPad->GetCalROC(isec);
 183     if (fROC[isec]){
 184       isOK&=fROC[isec]->LTMFilter(deltaRow, deltaPad,fraction,type,outlierROC,doEdge);
 185     }
 186   }
 187   return isOK;
 188 }
 189
 190 Bool_t  AliTPCCalPad::Convolute(Double_t sigmaPad, Double_t sigmaRow,  AliTPCCalPad*outlierPad, TF1 *fpad, TF1 *frow){
 191   //
 192   // replace constent with median in the neigborhood
 193   //
 194   Bool_t isOK=kTRUE;
 195   for (Int_t isec = 0; isec < kNsec; isec++) {
 196     AliTPCCalROC *outlierROC=(outlierPad==NULL)?NULL:outlierPad->GetCalROC(isec);
 197     if (fROC[isec]){
 198       isOK&=fROC[isec]->Convolute(sigmaPad,sigmaRow,outlierROC,fpad,frow);
 199     }
 200   }
 201   return isOK;
 202 }
 203
 204
 205 //_____________________________________________________________________________
 206 void AliTPCCalPad::Add(Float_t c1)
 207 {
 208     //
 209     // add constant c1 to all channels of all ROCs
 210     //
 211
 212     for (Int_t isec = 0; isec < kNsec; isec++) {
 213         if (fROC[isec]){
 214             fROC[isec]->Add(c1);
 215         }
 216     }
 217 }
 218
 219 //_____________________________________________________________________________
 220 void AliTPCCalPad::Multiply(Float_t c1)
 221 {
 222   //
 223   // multiply each channel of all ROCs with c1
 224   //
 225     for (Int_t isec = 0; isec < kNsec; isec++) {
 226         if (fROC[isec]){
 227             fROC[isec]->Multiply(c1);
 228         }
 229     }
 230 }
 231
 232 //_____________________________________________________________________________
 233 void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
 234 {
 235   //
 236   // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
 237   //  - pad by pad -
 238   //
 239     for (Int_t isec = 0; isec < kNsec; isec++) {
 240         if (fROC[isec] && pad->GetCalROC(isec)){
 241             fROC[isec]->Add(pad->GetCalROC(isec),c1);
 242         }
 243     }
 244 }
 245
 246 //_____________________________________________________________________________
 247 void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
 248 {
 249   //
 250   // multiply each channel of all ROCs with the coresponding channel of 'pad'
 251   //     - pad by pad -
 252   //
 253    for (Int_t isec = 0; isec < kNsec; isec++) {
 254         if (fROC[isec]){
 255             fROC[isec]->Multiply(pad->GetCalROC(isec));
 256         }
 257     }
 258 }
 259
 260 //_____________________________________________________________________________
 261 void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
 262 {
 263   //
 264   // divide each channel of all ROCs by the coresponding channel of 'pad'
 265   //     - pad by pad -
 266   //
 267     for (Int_t isec = 0; isec < kNsec; isec++) {
 268         if (fROC[isec]){
 269             fROC[isec]->Divide(pad->GetCalROC(isec));
 270         }
 271     }
 272 }
 273
 274 //_____________________________________________________________________________
 275 void AliTPCCalPad::Reset()
 276 {
 277   //
 278   // Reset all cal Rocs
 279   //
 280   for (Int_t isec = 0; isec < kNsec; isec++) {
 281     if (fROC[isec]){
 282       fROC[isec]->Reset();
 283     }
 284   }
 285 }
 286
 287 //_____________________________________________________________________________
 288 TGraph  *  AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
 289   //
 290   //   type=1 - mean
 291   //        2 - median
 292   //        3 - LTM
 293   //
 294   Int_t npoints = 0;
 295   for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++;
 296   TGraph * graph = new TGraph(npoints);
 297   npoints=0;
 298   for (Int_t isec=0;isec<72;isec++){
 299     if (!fROC[isec]) continue;
 300     if (type==0)  graph->SetPoint(npoints,isec,fROC[isec]->GetMean());
 301     if (type==1)  graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
 302     if (type==2)  graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));
 303     npoints++;
 304   }
 305
 306   graph->GetXaxis()->SetTitle("Sector");
 307   if (type==0) {
 308     graph->GetYaxis()->SetTitle("Mean");
 309     graph->SetMarkerStyle(22);
 310   }
 311   if (type==1) {
 312     graph->GetYaxis()->SetTitle("Median");
 313     graph->SetMarkerStyle(22);
 314   }
 315   if (type==2) {
 316       graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));
 317       graph->SetMarkerStyle(24);
 318   }
 319
 320   return graph;
 321 }
 322
 323 //_____________________________________________________________________________
 324 Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
 325 {
 326     //
 327     // Calculates mean and RMS of all ROCs
 328     //
 329     Double_t sum = 0, sum2 = 0, n=0, val=0;
 330     for (Int_t isec = 0; isec < kNsec; isec++) {
 331         AliTPCCalROC *calRoc = fROC[isec];
 332         if ( calRoc ){
 333             for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
 334                 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
 335                     val = calRoc->GetValue(irow,ipad);
 336                     sum+=val;
 337                     sum2+=val*val;
 338                     n++;
 339                 }
 340             }
 341
 342         }
 343     }
 344     Double_t n1 = 1./n;
 345     Double_t mean = sum*n1;
 346     rms  = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
 347     return mean;
 348 }
 349
 350
 351 //_____________________________________________________________________________
 352 Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
 353 {
 354     //
 355     // return mean of the mean of all ROCs
 356     //
 357     Double_t arr[kNsec];
 358     Int_t n=0;
 359     for (Int_t isec = 0; isec < kNsec; isec++) {
 360        AliTPCCalROC *calRoc = fROC[isec];
 361        if ( calRoc ){
 362           AliTPCCalROC* outlierROC = 0;
 363           if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 364                arr[n] = calRoc->GetMean(outlierROC);
 365           n++;
 366        }
 367     }
 368     return TMath::Mean(n,arr);
 369 }
 370
 371 //_____________________________________________________________________________
 372 Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
 373 {
 374     //
 375     // return mean of the RMS of all ROCs
 376     //
 377     Double_t arr[kNsec];
 378     Int_t n=0;
 379     for (Int_t isec = 0; isec < kNsec; isec++) {
 380        AliTPCCalROC *calRoc = fROC[isec];
 381        if ( calRoc ){
 382           AliTPCCalROC* outlierROC = 0;
 383           if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 384           arr[n] = calRoc->GetRMS(outlierROC);
 385           n++;
 386        }
 387     }
 388     return TMath::Mean(n,arr);
 389 }
 390
 391 //_____________________________________________________________________________
 392 Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
 393 {
 394     //
 395     // return mean of the median of all ROCs
 396     //
 397     Double_t arr[kNsec];
 398     Int_t n=0;
 399     for (Int_t isec = 0; isec < kNsec; isec++) {
 400        AliTPCCalROC *calRoc = fROC[isec];
 401        if ( calRoc ){
 402           AliTPCCalROC* outlierROC = 0;
 403           if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 404           arr[n] = calRoc->GetMedian(outlierROC);
 405           n++;
 406        }
 407     }
 408     return TMath::Mean(n,arr);
 409 }
 410
 411 //_____________________________________________________________________________
 412 Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
 413 {
 414     //
 415     // return mean of the LTM and sigma of all ROCs
 416     //
 417     Double_t arrm[kNsec];
 418     Double_t arrs[kNsec];
 419     Double_t *sTemp=0x0;
 420     Int_t n=0;
 421
 422     for (Int_t isec = 0; isec < kNsec; isec++) {
 423         AliTPCCalROC *calRoc = fROC[isec];
 424         if ( calRoc ){
 425             if ( sigma ) sTemp=arrs+n;
 426        AliTPCCalROC* outlierROC = 0;
 427        if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 428             arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
 429             n++;
 430         }
 431     }
 432     if ( sigma ) *sigma = TMath::Mean(n,arrs);
 433     return TMath::Mean(n,arrm);
 434 }
 435
 436 //_____________________________________________________________________________
 437 TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type, Int_t side){
 438   //
 439   // make 1D histo
 440   // type -1 = user defined range
 441   //       0 = nsigma cut nsigma=min
 442   //
 443   if (type>=0){
 444     if (type==0){
 445       // nsigma range
 446       Float_t mean  = GetMean();
 447       Float_t sigma = GetRMS();
 448       Float_t nsigma = TMath::Abs(min);
 449       min = mean-nsigma*sigma;
 450       max = mean+nsigma*sigma;
 451     }
 452     if (type==1){
 453       // fixed range
 454       Float_t mean   = GetMedian();
 455       Float_t  delta = min;
 456       min = mean-delta;
 457       max = mean+delta;
 458     }
 459     if (type==2){
 460       //
 461       // LTM mean +- nsigma
 462       //
 463       Double_t sigma;
 464       Float_t mean  = GetLTM(&sigma,max);
 465       sigma*=min;
 466       min = mean-sigma;
 467       max = mean+sigma;
 468     }
 469   }
 470   TString name=Form("%s Pad 1D",GetTitle());
 471   TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
 472     for (Int_t isec = 0; isec < kNsec; isec++) {
 473       if (side==1 && isec%36>18) continue;
 474       if (side==-1 && isec%36<18) continue;
 475         if (fROC[isec]){
 476             for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
 477                 UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
 478                 for (UInt_t ipad=0; ipad<npads; ipad++){
 479                     his->Fill(fROC[isec]->GetValue(irow,ipad));
 480                 }
 481             }
 482         }
 483     }
 484   return his;
 485 }
 486
 487 //_____________________________________________________________________________
 488 TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
 489   //
 490   // Make 2D graph
 491   // side  -  specify the side A = 0 C = 1
 492   // type  -  used types of determination of boundaries in z
 493   //
 494   Float_t kEpsilon = 0.000000000001;
 495   TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
 496   AliTPCROC * roc  = AliTPCROC::Instance();
 497   for (Int_t isec=0; isec<72; isec++){
 498     if (side==0 && isec%36>=18) continue;
 499     if (side>0 && isec%36<18) continue;
 500     if (fROC[isec]){
 501       AliTPCCalROC * calRoc = fROC[isec];
 502       for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
 503         for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
 504           if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
 505             Float_t xyz[3];
 506             roc->GetPositionGlobal(isec,irow,ipad,xyz);
 507             Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
 508             Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
 509             Float_t value = calRoc->GetValue(irow,ipad);
 510             his->SetBinContent(binx,biny,value);
 511           }
 512     }
 513   }
 514   his->SetXTitle("x (cm)");
 515   his->SetYTitle("y (cm)");
 516   return his;
 517 }
 518
 519
 520 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 {
 521    //
 522    // Loops over all AliTPCCalROCs and performs a localFit in each ROC
 523    // AliTPCCalPad with fit-data is returned
 524    // rowRadius and padRadius specifies a window around a given pad in one sector.
 525    // The data of this window are fitted with a parabolic function.
 526    // This function is evaluated at the pad's position.
 527    // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
 528    // rowRadius  -  radius - rows to be used for smoothing
 529    // padradius  -  radius - pads to be used for smoothing
 530    // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 531    // robust     -  robust method of fitting  - (much slower)
 532    // chi2Threshold: Threshold for chi2 when EvalRobust is called
 533    // robustFraction: Fraction of data that will be used in EvalRobust
 534    //
 535    //
 536    AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
 537    for (Int_t isec = 0; isec < 72; isec++){
 538       if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
 539       if (PadOutliers)
 540          pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
 541       else
 542          pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
 543    }
 544    return pad;
 545 }
 546
 547
 548 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){
 549    //
 550    // Loops over all AliTPCCalROCs and performs a globalFit in each ROC
 551    // AliTPCCalPad with fit-data is returned
 552    // chi2Threshold: Threshold for chi2 when EvalRobust is called
 553    // robustFraction: Fraction of data that will be used in EvalRobust
 554    // chi2Threshold: Threshold for chi2 when EvalRobust is called
 555    // robustFraction: Fraction of data that will be used in EvalRobust
 556    // err: error of the data points
 557    // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array
 558    //
 559    AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
 560    TVectorD fitParam(0);
 561    TMatrixD covMatrix(0,0);
 562    Float_t chi2 = 0;
 563    for (Int_t isec = 0; isec < 72; isec++){
 564       if (PadOutliers)
 565          GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
 566       else
 567          GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err);
 568
 569       AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec);
 570       pad->SetCalROC(roc);
 571       delete roc;
 572       if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec);
 573       if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec);
 574    }
 575    return pad;
 576 }
 577 //_____________________________________________________________________________
 578 TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula)
 579 {
 580   //
 581   // create an array of TFormulas for the each parameter of the fit function
 582   //
 583
 584   // split fit string in single parameters
 585   // find dimension of the fit:
 586   TString fitString(fitFormula);
 587   fitString.ReplaceAll("++","#");
 588   fitString.ReplaceAll(" ","");
 589   TObjArray *arrFitParams = fitString.Tokenize("#");
 590   Int_t ndim = arrFitParams->GetEntries();
 591   //create array of TFormulas to evaluate the parameters
 592   TObjArray *arrFitFormulas = new TObjArray(ndim);
 593   arrFitFormulas->SetOwner(kTRUE);
 594   for (Int_t idim=0;idim<ndim;++idim){
 595     TString s=((TObjString*)arrFitParams->At(idim))->GetString();
 596     s.ReplaceAll("gx","[0]");
 597     s.ReplaceAll("gy","[1]");
 598     s.ReplaceAll("lx","[2]");
 599     s.ReplaceAll("ly","[3]");
 600     s.ReplaceAll("sector","[4]");
 601     arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim);
 602   }
 603   delete arrFitParams;
 604
 605   return arrFitFormulas;
 606 }
 607 //_____________________________________________________________________________
 608 void AliTPCCalPad::EvalFormulaArray(const TObjArray &arrFitFormulas, TVectorD &results,
 609                                     const Int_t sec, const Int_t row, const Int_t pad)
 610 {
 611   //
 612   // evaluate the fit formulas
 613   //
 614   Int_t ndim=arrFitFormulas.GetEntries();
 615   results.ResizeTo(ndim);
 616
 617   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();  // to calculate the pad's position
 618   Float_t localXYZ[3];
 619   Float_t globalXYZ[3];
 620   tpcROCinstance->GetPositionLocal(sec, row, pad, localXYZ);
 621   tpcROCinstance->GetPositionGlobal(sec, row, pad, globalXYZ);
 622   //calculate parameter values
 623   for (Int_t idim=0;idim<ndim;++idim){
 624     TFormula *f=(TFormula*)arrFitFormulas.At(idim);
 625     f->SetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],sec);
 626     results[idim]=f->Eval(0);
 627   }
 628 }
 629 //_____________________________________________________________________________
 630 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){
 631   //
 632   // Performs a fit on both sides.
 633   // Valid information for the fitFormula are the variables
 634   // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName
 635   // - sector:         the sector number.
 636   //  eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on
 637   //
 638   // PadOutliers - pads with value !=0 are not used in fitting procedure
 639   // chi2Threshold: Threshold for chi2 when EvalRobust is called
 640   // robustFraction: Fraction of data that will be used in EvalRobust
 641   //
 642
 643   TObjArray* arrFitFormulas=CreateFormulaArray(fitFormula);
 644   Int_t ndim = arrFitFormulas->GetEntries();
 645   //resize output data arrays
 646   fitParamSideA.ResizeTo(ndim+1);
 647   fitParamSideC.ResizeTo(ndim+1);
 648   covMatrixSideA.ResizeTo(ndim+1,ndim+1);
 649   covMatrixSideC.ResizeTo(ndim+1,ndim+1);
 650   // create linear fitter for A- and C- Side
 651   TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
 652   TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
 653   fitterGA->StoreData(kTRUE);
 654   fitterGC->StoreData(kTRUE);
 655   //parameter values
 656   TVectorD parValues(ndim);
 657
 658   AliTPCCalROC *rocErr=0x0;
 659
 660   for (UInt_t isec = 0; isec<kNsec; ++isec){
 661     AliTPCCalROC *rocOut=PadOutliers->GetCalROC(isec);
 662     AliTPCCalROC *rocData=GetCalROC(isec);
 663     if (pointError) rocErr=pointError->GetCalROC(isec);
 664     if (!rocData) continue;
 665     for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
 666       for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
 667         //check for outliers
 668         if (rocOut && rocOut->GetValue(irow,ipad)) continue;
 669         //calculate parameter values
 670         EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad);
 671         //get value
 672         Float_t value=rocData->GetValue(irow,ipad);
 673         //point error
 674         Int_t err=1;
 675         if (rocErr) {
 676           err=TMath::Nint(rocErr->GetValue(irow,ipad));
 677           if (err==0) err=1;
 678         }
 679         //add points to the fitters
 680         if (isec/18%2==0){
 681           fitterGA->AddPoint(parValues.GetMatrixArray(),value,err);
 682         }else{
 683           fitterGC->AddPoint(parValues.GetMatrixArray(),value,err);
 684         }
 685       }
 686     }
 687   }
 688   if (robust){
 689     fitterGA->EvalRobust(robustFraction);
 690     fitterGC->EvalRobust(robustFraction);
 691   } else {
 692     fitterGA->Eval();
 693     fitterGC->Eval();
 694   }
 695   chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1));
 696   chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1));
 697   fitterGA->GetParameters(fitParamSideA);
 698   fitterGC->GetParameters(fitParamSideC);
 699   fitterGA->GetCovarianceMatrix(covMatrixSideA);
 700   fitterGC->GetCovarianceMatrix(covMatrixSideC);
 701
 702   delete arrFitFormulas;
 703   delete fitterGA;
 704   delete fitterGC;
 705
 706 }
 707 //
 708 AliTPCCalPad *AliTPCCalPad::CreateCalPadFit(const char* fitFormula, const TVectorD &fitParamSideA, const TVectorD &fitParamSideC)
 709 {
 710   //
 711   //
 712   //
 713   TObjArray *arrFitFormulas=CreateFormulaArray(fitFormula);
 714   Int_t ndim = arrFitFormulas->GetEntries();
 715   //check if dimension of fit formula and fit parameters agree
 716   if (ndim!=fitParamSideA.GetNrows()||ndim!=fitParamSideC.GetNrows()){
 717     printf("AliTPCCalPad::CreateCalPadFit: Dimensions of fit formula and fit Parameters does not match!");
 718     return 0;
 719   }
 720   //create cal pad
 721   AliTPCCalPad *pad=new AliTPCCalPad("fitResultPad",Form("Fit result: %s",fitFormula));
 722   //fill cal pad with fit results if requested
 723   for (UInt_t isec = 0; isec<kNsec; ++isec){
 724     AliTPCCalROC *roc=pad->GetCalROC(isec);
 725     for (UInt_t irow = 0; irow < roc->GetNrows(); irow++) {
 726       for (UInt_t ipad = 0; ipad < roc->GetNPads(irow); ipad++) {
 727         const TVectorD *fitPar=0;
 728         TVectorD fitResArray;
 729         if (isec/18%2==0){
 730           fitPar=&fitParamSideA;
 731         }else{
 732           fitPar=&fitParamSideC;
 733         }
 734         EvalFormulaArray(*arrFitFormulas,fitResArray, isec, irow, ipad);
 735         for (Int_t idim=0;idim<ndim;++idim)
 736           fitResArray(idim)*=(*fitPar)(idim);
 737         roc->SetValue(irow,ipad,fitResArray.Sum());
 738       }
 739     }
 740   }
 741   delete arrFitFormulas;
 742   return pad;
 743 }
 744
 745
 746
 747 TCanvas * AliTPCCalPad::MakeReportPadSector(TTree *chain, const char* varName, const char*varTitle, const char *axisTitle, Float_t min, Float_t max, const char *cutUser){
 748   //
 749   // Make a report - cal pads per sector
 750   // mean valeus per sector and local X
 751   //
 752   TH1* his=0;
 753   TLegend *legend = 0;
 754   TCanvas *canvas = new TCanvas(Form("Sector: %s",varTitle),Form("Sector: %s",varTitle),1500,1100);
 755
 756   canvas->Divide(2);
 757   chain->SetAlias("lX","lx.fElements");
 758   //
 759   canvas->cd(1);
 760   TString strDraw=varName;
 761   strDraw+=":lX";
 762   legend = new TLegend(0.5,0.50,0.9,0.9, Form("%s TPC A side", varTitle));
 763   for (Int_t isec=-1; isec<18; isec+=1){
 764     TCut cutSec=Form("sector%%36==%d",isec);
 765     cutSec+=cutUser;
 766     if (isec==-1) cutSec="sector%36<18";
 767     chain->SetMarkerColor(1+(isec+2)%5);
 768     chain->SetLineColor(1+(isec+2)%5);
 769     chain->SetMarkerStyle(25+(isec+2)%4);
 770     //
 771     chain->Draw(strDraw.Data(),cutSec,"profgoff");
 772     his=(TH1*)chain->GetHistogram()->Clone();
 773     delete chain->GetHistogram();
 774     his->SetMaximum(max);
 775     his->SetMinimum(min);
 776     his->GetXaxis()->SetTitle("R (cm)");
 777     his->GetYaxis()->SetTitle(axisTitle);
 778     his->SetTitle(Form("%s- sector %d",varTitle, isec));
 779     his->SetName(Form("%s- sector %d",varTitle, isec));
 780     if (isec==-1) his->SetTitle(Form("%s A side",varTitle));
 781     if (isec==-1) his->Draw();
 782     his->Draw("same");
 783     legend->AddEntry(his);
 784   }
 785   legend->Draw();
 786   canvas->cd(2);
 787   //
 788   legend = new TLegend(0.5,0.50,0.9,0.9, Form("%s TPC C side", varTitle));
 789   for (Int_t isec=-1; isec<18; isec+=1){
 790     TCut cutSec=Form("(sector+18)%%36==%d",isec);
 791     cutSec+=cutUser;
 792     if (isec==-1) cutSec="sector%36>18";
 793     chain->SetMarkerColor(1+(isec+2)%5);
 794     chain->SetLineColor(1+(isec+2)%5);
 795     chain->SetMarkerStyle(25+isec%4);
 796     //
 797     chain->Draw(strDraw.Data(),cutSec,"profgoff");
 798     his=(TH1*)chain->GetHistogram()->Clone();
 799     delete chain->GetHistogram();
 800     his->SetMaximum(max);
 801     his->SetMinimum(min);
 802     his->GetXaxis()->SetTitle("R (cm)");
 803     his->GetYaxis()->SetTitle(axisTitle);
 804     his->SetTitle(Form("%s- sector %d",varTitle,isec));
 805     his->SetName(Form("%s- sector %d",varTitle,isec));
 806     if (isec==-1) his->SetTitle(Form("%s C side",varTitle));
 807     if (isec==-1) his->Draw();
 808     his->Draw("same");
 809     legend->AddEntry(his);
 810   }
 811   legend->Draw();
 812   //
 813   //
 814   return canvas;
 815 }
 816
 817
 818 TCanvas * AliTPCCalPad::MakeReportPadSector2D(TTree *chain, const char* varName, const char*varTitle, const char *axisTitle, Float_t min, Float_t max, const char *cutUser){
 819   //
 820   // Make a report - cal pads per sector
 821   // 2D view
 822   // Input tree should be created using AliPreprocesorOnline before
 823   //
 824   TH1* his=0;
 825   TCanvas *canvas = new TCanvas(Form("%s2D",varTitle),Form("%s2D",varTitle),1500,1100);
 826   canvas->Divide(2);
 827   //
 828   TString strDraw=varName;
 829   strDraw+=":gy.fElements:gx.fElements>>his(250,-250,250,250,-250,250)";
 830   //
 831   TVirtualPad * pad=0;
 832   pad=canvas->cd(1);
 833   pad->SetMargin(0.15,0.15,0.15,0.15);
 834   TCut cut=cutUser;
 835   chain->Draw(strDraw.Data(),"sector%36<18"+cut,"profgoffcolz2");
 836   his=(TH1*)chain->GetHistogram()->Clone();
 837   delete chain->GetHistogram();
 838   his->SetMaximum(max);
 839   his->SetMinimum(min);
 840   his->GetXaxis()->SetTitle("x (cm)");
 841   his->GetYaxis()->SetTitle("y (cm)");
 842   his->GetZaxis()->SetTitle(axisTitle);
 843   his->SetTitle(Form("%s A side",varTitle));
 844   his->SetName(Form("%s A side",varTitle));
 845   his->Draw("colz2");
 846   //
 847   pad=canvas->cd(2);
 848   pad->SetMargin(0.15,0.15,0.15,0.15);
 849
 850   chain->Draw(strDraw.Data(),"sector%36>=18"+cut,"profgoffcolz2");
 851   his=(TH1*)chain->GetHistogram()->Clone();
 852   delete chain->GetHistogram();
 853   his->SetMaximum(max);
 854   his->SetMinimum(min);
 855   his->GetXaxis()->SetTitle("x (cm)");
 856   his->GetYaxis()->SetTitle("y (cm)");
 857   his->GetZaxis()->SetTitle(axisTitle);
 858   his->SetTitle(Form("%s C side",varTitle));
 859   his->SetName(Form("%s C side",varTitle));
 860   his->Draw("colz2");
 861   //
 862   //
 863   return canvas;
 864 }
 865
 866 void  AliTPCCalPad::Draw(Option_t* option){
 867   //
 868   // Draw function - standard 2D view
 869   //
 870   TH1* his=0;
 871   TCanvas *canvas = new TCanvas(Form("%s2D",GetTitle()),Form("%s2D",GetTitle()),900,900);
 872   canvas->Divide(2,2);
 873   //
 874   //
 875   TVirtualPad * pad=0;
 876   pad=canvas->cd(1);
 877   pad->SetMargin(0.15,0.15,0.15,0.15);
 878   his=MakeHisto2D(0);
 879   his->GetXaxis()->SetTitle("x (cm)");
 880   his->GetYaxis()->SetTitle("y (cm)");
 881   his->GetZaxis()->SetTitle(GetTitle());
 882   his->SetTitle(Form("%s A side",GetTitle()));
 883   his->SetName(Form("%s A side",GetTitle()));
 884   his->Draw(option);
 885   //
 886   pad=canvas->cd(2);
 887   pad->SetMargin(0.15,0.15,0.15,0.15);
 888   his=MakeHisto2D(1);
 889   his->GetXaxis()->SetTitle("x (cm)");
 890   his->GetYaxis()->SetTitle("y (cm)");
 891   his->GetZaxis()->SetTitle(GetTitle());
 892   his->SetTitle(Form("%s C side",GetTitle()));
 893   his->SetName(Form("%s C side",GetTitle()));
 894   his->Draw(option);
 895   //
 896   pad=canvas->cd(3);
 897   pad->SetMargin(0.15,0.15,0.15,0.15);
 898   his=MakeHisto1D(-8,8,0,1);
 899   his->GetXaxis()->SetTitle(GetTitle());
 900   his->SetTitle(Form("%s A side",GetTitle()));
 901   his->SetName(Form("%s A side",GetTitle()));
 902   his->Draw("err");
 903   //
 904   pad=canvas->cd(4);
 905   pad->SetMargin(0.15,0.15,0.15,0.15);
 906   his=MakeHisto1D(-8,8,0,-1);
 907   his->GetXaxis()->SetTitle(GetTitle());
 908   his->SetTitle(Form("%s C side",GetTitle()));
 909   his->SetName(Form("%s C side",GetTitle()));
 910   his->Draw("err");
 911
 912
 913 }
 914
 915
 916 AliTPCCalPad * AliTPCCalPad::MakeCalPadFromHistoRPHI(TH2 * hisA, TH2* hisC){
 917   //
 918   // Make cal pad from r-phi histograms
 919   //
 920   AliTPCROC *proc= AliTPCROC::Instance();
 921   AliTPCCalPad *calPad = new AliTPCCalPad("his","his");
 922   Float_t globalPos[3];
 923   for (Int_t isec=0; isec<72; isec++){
 924     AliTPCCalROC* calRoc  = calPad->GetCalROC(isec);
 925     TH2 * his = ((isec%36<18) ? hisA:hisC);
 926     for (UInt_t irow=0; irow<calRoc->GetNrows(); irow+=1){
 927       Int_t jrow=irow;
 928       if (isec>=36) jrow+=63;
 929       for (UInt_t ipad=0;ipad<proc->GetNPads(isec,irow);ipad+=1){
 930         proc->GetPositionGlobal(isec,irow,ipad, globalPos);
 931         Double_t phi=TMath::ATan2(globalPos[1],globalPos[0]);
 932         //if (phi<0) phi+=TMath::Pi()*2;
 933         Int_t bin=his->FindBin(phi,jrow);
 934         Float_t value= his->GetBinContent(bin);
 935         calRoc->SetValue(irow,ipad,value);
 936       }
 937     }
 938   }
 939   return calPad;
 940 }
 941
 942 AliTPCCalPad *AliTPCCalPad::MakePadFromTree(TTree * treePad, const char *query, const char* name, Bool_t doFast){
 943   //
 944   // make cal pad from the tree
 945   //
 946   if (!treePad){
 947     ::Error("AliTPCCalPad::MakePadFromTree(TTree * treePad, const char *query, const char* name)","Input tree is missing");
 948     return 0;
 949   }
 950   if (treePad->GetEntries()!=kNsec) return 0;
 951   AliTPCCalPad * calPad= new AliTPCCalPad(name,name);
 952   if (name) calPad->SetName(name);
 953   if (!doFast){
 954     for (Int_t iSec=0; iSec<72; iSec++){
 955       AliTPCCalROC* calROC  = calPad->GetCalROC(iSec);
 956       UInt_t nchannels = (UInt_t)treePad->Draw(query,"1","goff",1,iSec);
 957       if (nchannels!=calROC->GetNchannels()) {
 958         ::Error("AliTPCCalPad::MakePad",TString::Format("%s\t:Wrong query sector\t%d\t%d",treePad->GetName(),iSec,nchannels).Data());
 959         break;
 960       }
 961       for (UInt_t index=0; index<nchannels; index++) calROC->SetValue(index,treePad->GetV1()[index]);
 962     }
 963   }else{
 964     UInt_t nchannelsTree = (UInt_t)treePad->Draw(query,"1","goff");
 965     UInt_t nchannelsAll=0;
 966     for (Int_t iSec=0; iSec<72; iSec++){
 967       AliTPCCalROC* calROC  = calPad->GetCalROC(iSec);
 968       UInt_t nchannels=calROC->GetNchannels();
 969       for (UInt_t index=0; index<nchannels; index++) {
 970         if (nchannelsAll<=nchannelsTree)calROC->SetValue(index,treePad->GetV1()[nchannelsAll]);
 971         nchannelsAll++;
 972       }
 973     }
 974     if (nchannelsAll>nchannelsTree){
 975       ::Error("AliTPCCalPad::MakePad",TString::Format("%s\t:Wrong query: cout mismatch\t%d\t%d",query, nchannelsAll,nchannelsTree).Data());
 976     }
 977   }
 978   return calPad;
 979 }
 980
 981 void AliTPCCalPad::AddFriend(TTree * treePad, const char *friendName, const char *fname){
 982   //
 983   //
 984   //
 985   TObjArray *fArray = new TObjArray(1);
 986   fArray->AddLast(this);
 987   this->SetName(friendName);
 988   AliTPCCalibViewer::MakeTree(fname, fArray,0);
 989   TFile * f = TFile::Open(fname);
 990   TTree * tree = (TTree*)f->Get("calPads");
 991   treePad->AddFriend(tree,friendName);
 992   //  tree->AddFriend(TString::Format("%s = calPads",friendName).Data(),fname);
 993 }