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