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
  36 ClassImp(AliTPCCalPad)
  37
  38 //_____________________________________________________________________________
  39 AliTPCCalPad::AliTPCCalPad():TNamed()
  40 {
  41   //
  42   // AliTPCCalPad default constructor
  43   //
  44
  45   for (Int_t isec = 0; isec < kNsec; isec++) {
  46     fROC[isec] = 0;
  47   }
  48
  49 }
  50
  51 //_____________________________________________________________________________
  52 AliTPCCalPad::AliTPCCalPad(const Text_t *name, const Text_t *title)
  53                 :TNamed(name,title)
  54 {
  55   //
  56   // AliTPCCalPad constructor
  57   //
  58   for (Int_t isec = 0; isec < kNsec; isec++) {
  59     fROC[isec] = new AliTPCCalROC(isec);
  60   }
  61 }
  62
  63
  64 //_____________________________________________________________________________
  65 AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c)
  66 {
  67   //
  68   // AliTPCCalPad copy constructor
  69   //
  70
  71   for (Int_t isec = 0; isec < kNsec; isec++) {
  72          fROC[isec] = 0;
  73      if (c.fROC[isec])
  74        fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
  75   }
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed()
  80 {
  81   //
  82   // AliTPCCalPad default constructor
  83   //
  84
  85   for (Int_t isec = 0; isec < kNsec; isec++) {
  86     fROC[isec] = (AliTPCCalROC *)array->At(isec);
  87   }
  88
  89 }
  90
  91
  92 ///_____________________________________________________________________________
  93 AliTPCCalPad::~AliTPCCalPad()
  94 {
  95   //
  96   // AliTPCCalPad destructor
  97   //
  98
  99   for (Int_t isec = 0; isec < kNsec; isec++) {
 100     if (fROC[isec]) {
 101       delete fROC[isec];
 102       fROC[isec] = 0;
 103     }
 104   }
 105
 106 }
 107
 108 //_____________________________________________________________________________
 109 AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c)
 110 {
 111   //
 112   // Assignment operator
 113   //
 114
 115   if (this != &c) ((AliTPCCalPad &) c).Copy(*this);
 116   return *this;
 117
 118 }
 119
 120 //_____________________________________________________________________________
 121 void AliTPCCalPad::Copy(TObject &c) const
 122 {
 123   //
 124   // Copy function
 125   //
 126
 127   for (Int_t isec = 0; isec < kNsec; isec++) {
 128     if (fROC[isec]) {
 129       fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]);
 130     }
 131   }
 132   TObject::Copy(c);
 133 }
 134
 135
 136 void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){
 137    //
 138    // Set AliTPCCalROC copies values from 'roc'
 139    // if sector == -1 the sector specified in 'roc' is used
 140    // else sector specified in 'roc' is ignored and specified sector is filled
 141    //
 142    if (sector == -1) sector = roc->GetSector();
 143    for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++)
 144       fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel));
 145 }
 146
 147
 148
 149 //_____________________________________________________________________________
 150 void AliTPCCalPad::Add(Float_t c1)
 151 {
 152     //
 153     // add constant c1 to all channels of all ROCs
 154     //
 155
 156     for (Int_t isec = 0; isec < kNsec; isec++) {
 157         if (fROC[isec]){
 158             fROC[isec]->Add(c1);
 159         }
 160     }
 161 }
 162
 163 //_____________________________________________________________________________
 164 void AliTPCCalPad::Multiply(Float_t c1)
 165 {
 166   //
 167   // multiply each channel of all ROCs with c1
 168   //
 169     for (Int_t isec = 0; isec < kNsec; isec++) {
 170         if (fROC[isec]){
 171             fROC[isec]->Multiply(c1);
 172         }
 173     }
 174 }
 175
 176 //_____________________________________________________________________________
 177 void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
 178 {
 179   //
 180   // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
 181   //  - pad by pad -
 182   //
 183     for (Int_t isec = 0; isec < kNsec; isec++) {
 184         if (fROC[isec]){
 185             fROC[isec]->Add(pad->GetCalROC(isec),c1);
 186         }
 187     }
 188 }
 189
 190 //_____________________________________________________________________________
 191 void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
 192 {
 193   //
 194   // multiply each channel of all ROCs with the coresponding channel of 'pad'
 195   //     - pad by pad -
 196   //
 197    for (Int_t isec = 0; isec < kNsec; isec++) {
 198         if (fROC[isec]){
 199             fROC[isec]->Multiply(pad->GetCalROC(isec));
 200         }
 201     }
 202 }
 203
 204 //_____________________________________________________________________________
 205 void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
 206 {
 207   //
 208   // divide each channel of all ROCs by the coresponding channel of 'pad'
 209   //     - pad by pad -
 210   //
 211     for (Int_t isec = 0; isec < kNsec; isec++) {
 212         if (fROC[isec]){
 213             fROC[isec]->Divide(pad->GetCalROC(isec));
 214         }
 215     }
 216 }
 217
 218 //_____________________________________________________________________________
 219 TGraph  *  AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
 220   //
 221   //   type=1 - mean
 222   //        2 - median
 223   //        3 - LTM
 224   //
 225   Int_t npoints = 0;
 226   for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++;
 227   TGraph * graph = new TGraph(npoints);
 228   npoints=0;
 229   for (Int_t isec=0;isec<72;isec++){
 230     if (!fROC[isec]) continue;
 231     if (type==0)  graph->SetPoint(npoints,isec,fROC[isec]->GetMean());
 232     if (type==1)  graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
 233     if (type==2)  graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));
 234     npoints++;
 235   }
 236
 237   graph->GetXaxis()->SetTitle("Sector");
 238   if (type==0) {
 239     graph->GetYaxis()->SetTitle("Mean");
 240     graph->SetMarkerStyle(22);
 241   }
 242   if (type==1) {
 243     graph->GetYaxis()->SetTitle("Median");
 244     graph->SetMarkerStyle(22);
 245   }
 246   if (type==2) {
 247       graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));
 248       graph->SetMarkerStyle(24);
 249   }
 250
 251   return graph;
 252 }
 253
 254 //_____________________________________________________________________________
 255 Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
 256 {
 257     //
 258     // Calculates mean and RMS of all ROCs
 259     //
 260     Double_t sum = 0, sum2 = 0, n=0, val=0;
 261     for (Int_t isec = 0; isec < kNsec; isec++) {
 262         AliTPCCalROC *calRoc = fROC[isec];
 263         if ( calRoc ){
 264             for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
 265                 for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
 266                     val = calRoc->GetValue(irow,ipad);
 267                     sum+=val;
 268                     sum2+=val*val;
 269                     n++;
 270                 }
 271             }
 272
 273         }
 274     }
 275     Double_t n1 = 1./n;
 276     Double_t mean = sum*n1;
 277     rms  = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
 278     return mean;
 279 }
 280
 281
 282 //_____________________________________________________________________________
 283 Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
 284 {
 285     //
 286     // return mean of the mean of all ROCs
 287     //
 288     Double_t arr[kNsec];
 289     Int_t n=0;
 290     for (Int_t isec = 0; isec < kNsec; isec++) {
 291        AliTPCCalROC *calRoc = fROC[isec];
 292        if ( calRoc ){
 293           AliTPCCalROC* outlierROC = 0;
 294           if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 295                arr[n] = calRoc->GetMean(outlierROC);
 296           n++;
 297        }
 298     }
 299     return TMath::Mean(n,arr);
 300 }
 301
 302 //_____________________________________________________________________________
 303 Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
 304 {
 305     //
 306     // return mean of the RMS of all ROCs
 307     //
 308     Double_t arr[kNsec];
 309     Int_t n=0;
 310     for (Int_t isec = 0; isec < kNsec; isec++) {
 311        AliTPCCalROC *calRoc = fROC[isec];
 312        if ( calRoc ){
 313           AliTPCCalROC* outlierROC = 0;
 314           if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 315           arr[n] = calRoc->GetRMS(outlierROC);
 316           n++;
 317        }
 318     }
 319     return TMath::Mean(n,arr);
 320 }
 321
 322 //_____________________________________________________________________________
 323 Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
 324 {
 325     //
 326     // return mean of the median of all ROCs
 327     //
 328     Double_t arr[kNsec];
 329     Int_t n=0;
 330     for (Int_t isec = 0; isec < kNsec; isec++) {
 331        AliTPCCalROC *calRoc = fROC[isec];
 332        if ( calRoc ){
 333           AliTPCCalROC* outlierROC = 0;
 334           if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 335           arr[n] = calRoc->GetMedian(outlierROC);
 336           n++;
 337        }
 338     }
 339     return TMath::Mean(n,arr);
 340 }
 341
 342 //_____________________________________________________________________________
 343 Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
 344 {
 345     //
 346     // return mean of the LTM and sigma of all ROCs
 347     //
 348     Double_t arrm[kNsec];
 349     Double_t arrs[kNsec];
 350     Double_t *sTemp=0x0;
 351     Int_t n=0;
 352
 353     for (Int_t isec = 0; isec < kNsec; isec++) {
 354         AliTPCCalROC *calRoc = fROC[isec];
 355         if ( calRoc ){
 356             if ( sigma ) sTemp=arrs+n;
 357        AliTPCCalROC* outlierROC = 0;
 358        if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 359             arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
 360             n++;
 361         }
 362     }
 363     if ( sigma ) *sigma = TMath::Mean(n,arrs);
 364     return TMath::Mean(n,arrm);
 365 }
 366
 367 //_____________________________________________________________________________
 368 TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){
 369   //
 370   // make 1D histo
 371   // type -1 = user defined range
 372   //       0 = nsigma cut nsigma=min
 373   //
 374   if (type>=0){
 375     if (type==0){
 376       // nsigma range
 377       Float_t mean  = GetMean();
 378       Float_t sigma = GetRMS();
 379       Float_t nsigma = TMath::Abs(min);
 380       min = mean-nsigma*sigma;
 381       max = mean+nsigma*sigma;
 382     }
 383     if (type==1){
 384       // fixed range
 385       Float_t mean   = GetMedian();
 386       Float_t  delta = min;
 387       min = mean-delta;
 388       max = mean+delta;
 389     }
 390     if (type==2){
 391       //
 392       // LTM mean +- nsigma
 393       //
 394       Double_t sigma;
 395       Float_t mean  = GetLTM(&sigma,max);
 396       sigma*=min;
 397       min = mean-sigma;
 398       max = mean+sigma;
 399     }
 400   }
 401   char  name[1000];
 402   sprintf(name,"%s Pad 1D",GetTitle());
 403   TH1F * his = new TH1F(name,name,100, min,max);
 404     for (Int_t isec = 0; isec < kNsec; isec++) {
 405         if (fROC[isec]){
 406             for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
 407                 UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
 408                 for (UInt_t ipad=0; ipad<npads; ipad++){
 409                     his->Fill(fROC[isec]->GetValue(irow,ipad));
 410                 }
 411             }
 412         }
 413     }
 414   return his;
 415 }
 416
 417 //_____________________________________________________________________________
 418 TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
 419   //
 420   // Make 2D graph
 421   // side  -  specify the side A = 0 C = 1
 422   // type  -  used types of determination of boundaries in z
 423   //
 424   Float_t kEpsilon = 0.000000000001;
 425   TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
 426   AliTPCROC * roc  = AliTPCROC::Instance();
 427   for (Int_t isec=0; isec<72; isec++){
 428     if (side==0 && isec%36>=18) continue;
 429     if (side>0 && isec%36<18) continue;
 430     if (fROC[isec]){
 431       AliTPCCalROC * calRoc = fROC[isec];
 432       for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
 433         for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
 434           if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
 435             Float_t xyz[3];
 436             roc->GetPositionGlobal(isec,irow,ipad,xyz);
 437             Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
 438             Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
 439             Float_t value = calRoc->GetValue(irow,ipad);
 440             his->SetBinContent(binx,biny,value);
 441           }
 442     }
 443   }
 444   his->SetXTitle("x (cm)");
 445   his->SetYTitle("y (cm)");
 446   return his;
 447 }
 448
 449
 450 AliTPCCalPad* AliTPCCalPad::LocalFit(const char* padName, Int_t rowRadius, Int_t padRadius, AliTPCCalPad* PadOutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction){
 451    //
 452    // Loops over all AliTPCCalROCs and performs a localFit in each ROC
 453    // AliTPCCalPad with fit-data is returned
 454    // rowRadius and padRadius specifies a window around a given pad in one sector.
 455    // The data of this window are fitted with a parabolic function.
 456    // This function is evaluated at the pad's position.
 457    // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
 458    // rowRadius  -  radius - rows to be used for smoothing
 459    // padradius  -  radius - pads to be used for smoothing
 460    // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 461    // robust     -  robust method of fitting  - (much slower)
 462    // chi2Threshold: Threshold for chi2 when EvalRobust is called
 463    // robustFraction: Fraction of data that will be used in EvalRobust
 464    //
 465    //
 466    AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
 467    for (Int_t isec = 0; isec < 72; isec++){
 468       if (PadOutliers)
 469          SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
 470       else
 471          SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
 472    }
 473    return pad;
 474 }
 475
 476
 477 AliTPCCalPad* AliTPCCalPad::GlobalFit(const char* padName, AliTPCCalPad* PadOutliers, Bool_t robust, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction){
 478    //
 479    // Loops over all AliTPCCalROCs and performs a globalFit in each ROC
 480    // AliTPCCalPad with fit-data is returned
 481    // chi2Threshold: Threshold for chi2 when EvalRobust is called
 482    // robustFraction: Fraction of data that will be used in EvalRobust
 483    // chi2Threshold: Threshold for chi2 when EvalRobust is called
 484    // robustFraction: Fraction of data that will be used in EvalRobust
 485    //
 486    AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
 487    TVectorD fitParam(0);
 488    TMatrixD covMatrix(0,0);
 489    Float_t chi2 = 0;
 490    for (Int_t isec = 0; isec < 72; isec++){
 491       if (PadOutliers)
 492          GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction);
 493       else
 494          GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction);
 495       pad->SetCalROC(AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec));
 496    }
 497    return pad;
 498 }
 499
 500
 501 void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, Int_t fitType, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction){
 502   //
 503   // Makes a  GlobalFit over each side and return fit-parameters, covariance and chi2 for each side
 504   // fitType == 0: fit plane function
 505   // fitType == 1: fit parabolic function
 506   // PadOutliers - pads with value !=0 are not used in fitting procedure
 507   // chi2Threshold: Threshold for chi2 when EvalRobust is called
 508   // robustFraction: Fraction of data that will be used in EvalRobust
 509   //
 510   TLinearFitter* fitterGA = 0;
 511   TLinearFitter* fitterGC = 0;
 512
 513   if (fitType  == 1) {
 514     fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
 515     fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
 516   }
 517   else {
 518     fitterGA = new TLinearFitter(3,"x0++x1++x2");
 519     fitterGC = new TLinearFitter(3,"x0++x1++x2");
 520   }
 521   fitterGA->StoreData(kTRUE);
 522   fitterGC->StoreData(kTRUE);
 523   fitterGA->ClearPoints();
 524   fitterGC->ClearPoints();
 525   Double_t xx[6];
 526   Int_t    npointsA=0;
 527   Int_t    npointsC=0;
 528
 529   Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position
 530   Float_t lx, ly;  // pads position
 531
 532   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();  // to calculate the pad's position
 533
 534   // loop over all sectors and pads and read data into fitterGA and fitterGC
 535   if (fitType == 1) {
 536   // parabolic fit
 537     fitParamSideA.ResizeTo(6);
 538     fitParamSideC.ResizeTo(6);
 539     covMatrixSideA.ResizeTo(6,6);
 540     covMatrixSideC.ResizeTo(6,6);
 541     for (UInt_t isec = 0; isec<72; isec++){
 542       for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
 543          for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
 544             // fill fitterG
 545             tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY);   // calculate position localXY by sector, pad and row number
 546             lx = localXY[0];
 547             ly = localXY[1];
 548             xx[0] = 1;
 549             xx[1] = lx;
 550             xx[2] = ly;
 551             xx[3] = lx*lx;
 552             xx[4] = ly*ly;
 553             xx[5] = lx*ly;
 554             if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
 555             // if given pad is no outlier, add it to TLinearFitter, decide to which of both
 556 //                sector  0 - 17: IROC, A
 557 //                sector 18 - 35: IROC, C
 558 //                sector 36 - 53: OROC, A
 559 //                sector 54 - 71: CROC, C
 560                if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A
 561                   npointsA++;
 562                   fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
 563                }
 564                else { // side C
 565                   npointsC++;
 566                   fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
 567                }
 568             }
 569          }
 570       }
 571     }
 572   }
 573   else {
 574   // linear fit
 575     fitParamSideA.ResizeTo(3);
 576     fitParamSideC.ResizeTo(3);
 577     covMatrixSideA.ResizeTo(3,3);
 578     covMatrixSideC.ResizeTo(3,3);
 579
 580     for (UInt_t isec = 0; isec<72; isec++){
 581       for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
 582          for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
 583             // fill fitterG
 584             tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY);   // calculate position localXY by sector, pad and row number
 585             lx = localXY[0];
 586             ly = localXY[1];
 587             xx[0] = 1;
 588             xx[1] = lx;
 589             xx[2] = ly;
 590             if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
 591             // if given pad is no outlier, add it to TLinearFitter, decide to which of both
 592 //                sector  0 - 17: IROC, A
 593 //                sector 18 - 35: IROC, C
 594 //                sector 36 - 53: OROC, A
 595 //                sector 54 - 71: CROC, C
 596                if (isec <= 17 || (isec >= 36 && isec <= 53)) {
 597                // Side A
 598                   npointsA++;
 599                   fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
 600                }
 601                else {
 602                // side C
 603                   npointsC++;
 604                   fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
 605                }
 606             }
 607          }
 608       }
 609     }
 610   }
 611
 612   fitterGA->Eval();
 613   fitterGC->Eval();
 614   fitterGA->GetParameters(fitParamSideA);
 615   fitterGC->GetParameters(fitParamSideC);
 616   fitterGA->GetCovarianceMatrix(covMatrixSideA);
 617   fitterGC->GetCovarianceMatrix(covMatrixSideC);
 618   if (fitType == 1){
 619     chi2SideA = fitterGA->GetChisquare()/(npointsA-6.);
 620     chi2SideC = fitterGC->GetChisquare()/(npointsC-6.);
 621   }
 622   else {
 623    chi2SideA = fitterGA->GetChisquare()/(npointsA-3.);
 624    chi2SideC = fitterGC->GetChisquare()/(npointsC-3.);
 625   }
 626   if (robust && chi2SideA > chi2Threshold) {
 627     //    std::cout << "robust fitter called... " << std::endl;
 628     fitterGA->EvalRobust(robustFraction);
 629     fitterGA->GetParameters(fitParamSideA);
 630   }
 631   if (robust && chi2SideC > chi2Threshold) {
 632     //    std::cout << "robust fitter called... " << std::endl;
 633     fitterGC->EvalRobust(robustFraction);
 634     fitterGC->GetParameters(fitParamSideC);
 635   }
 636   delete fitterGA;
 637   delete fitterGC;
 638 }
 639
 640
 641 void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array, const char * mapFileName, AliTPCCalPad* outlierPad, Float_t ltmFraction) {
 642   //
 643   // Write a tree with all available information
 644   // im mapFileName is speciefied, the Map information are also written to the tree
 645   // pads specified in outlierPad are not used for calculating statistics
 646   //  - this function will be moved to AliTPCCalibViewer -
 647   //  - DO NOT USE THIS FUNCTION ANY MORE -
 648   //
 649    AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 650
 651    TObjArray* mapIROCs = 0;
 652    TObjArray* mapOROCs = 0;
 653    TVectorF *mapIROCArray = 0;
 654    TVectorF *mapOROCArray = 0;
 655    Int_t mapEntries = 0;
 656    TString* mapNames = 0;
 657
 658    if (mapFileName) {
 659       TFile mapFile(mapFileName, "read");
 660
 661       TList* listOfROCs = mapFile.GetListOfKeys();
 662       mapEntries = listOfROCs->GetEntries()/2;
 663       mapIROCs = new TObjArray(mapEntries*2);
 664       mapOROCs = new TObjArray(mapEntries*2);
 665       mapIROCArray = new TVectorF[mapEntries];
 666       mapOROCArray = new TVectorF[mapEntries];
 667
 668       mapNames = new TString[mapEntries];
 669       for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
 670          TString ROCname(((TKey*)(listOfROCs->At(ivalue*2)))->GetName());
 671          ROCname.Remove(ROCname.Length()-4, 4);
 672          mapIROCs->AddAt((AliTPCCalROC*)mapFile.Get((ROCname + "IROC").Data()), ivalue);
 673          mapOROCs->AddAt((AliTPCCalROC*)mapFile.Get((ROCname + "OROC").Data()), ivalue);
 674          mapNames[ivalue].Append(ROCname);
 675       }
 676
 677       for (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
 678          mapIROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(0));
 679          mapOROCArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(36));
 680
 681          for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(0); ichannel++)
 682             (mapIROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapIROCs->At(ivalue)))->GetValue(ichannel);
 683          for (UInt_t ichannel = 0; ichannel < tpcROCinstance->GetNChannels(36); ichannel++)
 684             (mapOROCArray[ivalue])[ichannel] = ((AliTPCCalROC*)(mapOROCs->At(ivalue)))->GetValue(ichannel);
 685       }
 686
 687    } //  if (mapFileName)
 688
 689    TTreeSRedirector cstream(fileName);
 690    Int_t arrayEntries = array->GetEntries();
 691
 692    TString* names = new TString[arrayEntries];
 693    for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
 694       names[ivalue].Append(((AliTPCCalPad*)array->At(ivalue))->GetName());
 695
 696    for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++) {
 697       //
 698       // get statistic for given sector
 699       //
 700       TVectorF median(arrayEntries);
 701       TVectorF mean(arrayEntries);
 702       TVectorF rms(arrayEntries);
 703       TVectorF ltm(arrayEntries);
 704       TVectorF ltmrms(arrayEntries);
 705       TVectorF medianWithOut(arrayEntries);
 706       TVectorF meanWithOut(arrayEntries);
 707       TVectorF rmsWithOut(arrayEntries);
 708       TVectorF ltmWithOut(arrayEntries);
 709       TVectorF ltmrmsWithOut(arrayEntries);
 710
 711       TVectorF *vectorArray = new TVectorF[arrayEntries];
 712       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++)
 713          vectorArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
 714
 715       for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
 716          AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
 717          AliTPCCalROC* calROC = calPad->GetCalROC(isector);
 718          AliTPCCalROC* outlierROC = 0;
 719          if (outlierPad) outlierROC = outlierPad->GetCalROC(isector);
 720          if (calROC) {
 721             median[ivalue] = calROC->GetMedian();
 722             mean[ivalue] = calROC->GetMean();
 723             rms[ivalue] = calROC->GetRMS();
 724             Double_t ltmrmsValue = 0;
 725             ltm[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction);
 726             ltmrms[ivalue] = ltmrmsValue;
 727             if (outlierROC) {
 728                medianWithOut[ivalue] = calROC->GetMedian(outlierROC);
 729                meanWithOut[ivalue] = calROC->GetMean(outlierROC);
 730                rmsWithOut[ivalue] = calROC->GetRMS(outlierROC);
 731                ltmrmsValue = 0;
 732                ltmWithOut[ivalue] = calROC->GetLTM(&ltmrmsValue, ltmFraction, outlierROC);
 733                ltmrmsWithOut[ivalue] = ltmrmsValue;
 734             }
 735          }
 736          else {
 737             median[ivalue] = 0.;
 738             mean[ivalue] = 0.;
 739             rms[ivalue] = 0.;
 740             ltm[ivalue] = 0.;
 741             ltmrms[ivalue] = 0.;
 742             medianWithOut[ivalue] = 0.;
 743             meanWithOut[ivalue] = 0.;
 744             rmsWithOut[ivalue] = 0.;
 745             ltmWithOut[ivalue] = 0.;
 746             ltmrmsWithOut[ivalue] = 0.;
 747          }
 748       }
 749
 750       //
 751       // fill vectors of variable per pad
 752       //
 753       TVectorF *posArray = new TVectorF[8];
 754       for (Int_t ivalue = 0; ivalue < 8; ivalue++)
 755          posArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector));
 756
 757       Float_t posG[3] = {0};
 758       Float_t posL[3] = {0};
 759       Int_t ichannel = 0;
 760       for (UInt_t irow = 0; irow < tpcROCinstance->GetNRows(isector); irow++) {
 761          for (UInt_t ipad = 0; ipad < tpcROCinstance->GetNPads(isector, irow); ipad++) {
 762             tpcROCinstance->GetPositionLocal(isector, irow, ipad, posL);
 763             tpcROCinstance->GetPositionGlobal(isector, irow, ipad, posG);
 764             posArray[0][ichannel] = irow;
 765             posArray[1][ichannel] = ipad;
 766             posArray[2][ichannel] = posL[0];
 767             posArray[3][ichannel] = posL[1];
 768             posArray[4][ichannel] = posG[0];
 769             posArray[5][ichannel] = posG[1];
 770             posArray[6][ichannel] = (Int_t)(ipad - (Double_t)(tpcROCinstance->GetNPads(isector, irow))/2);
 771             posArray[7][ichannel] = ichannel;
 772
 773             // loop over array containing AliTPCCalPads
 774             for (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
 775                AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue);
 776                AliTPCCalROC* calROC = calPad->GetCalROC(isector);
 777                if (calROC)
 778                   (vectorArray[ivalue])[ichannel] = calROC->GetValue(irow, ipad);
 779                else
 780                   (vectorArray[ivalue])[ichannel] = 0;
 781             }
 782             ichannel++;
 783          }
 784       }
 785
 786       cstream << "calPads" <<
 787          "sector=" << isector;
 788
 789       for  (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
 790          cstream << "calPads" <<
 791             (Char_t*)((names[ivalue] + "_Median=").Data()) << median[ivalue] <<
 792             (Char_t*)((names[ivalue] + "_Mean=").Data()) << mean[ivalue] <<
 793             (Char_t*)((names[ivalue] + "_RMS=").Data()) << rms[ivalue] <<
 794             (Char_t*)((names[ivalue] + "_LTM=").Data()) << ltm[ivalue] <<
 795             (Char_t*)((names[ivalue] + "_RMS_LTM=").Data()) << ltmrms[ivalue];
 796          if (outlierPad) {
 797             cstream << "calPads" <<
 798                (Char_t*)((names[ivalue] + "_Median_OutlierCutted=").Data()) << medianWithOut[ivalue] <<
 799                (Char_t*)((names[ivalue] + "_Mean_OutlierCutted=").Data()) << meanWithOut[ivalue] <<
 800                (Char_t*)((names[ivalue] + "_RMS_OutlierCutted=").Data()) << rmsWithOut[ivalue] <<
 801                (Char_t*)((names[ivalue] + "_LTM_OutlierCutted=").Data()) << ltmWithOut[ivalue] <<
 802                (Char_t*)((names[ivalue] + "_RMS_LTM_OutlierCutted=").Data()) << ltmrmsWithOut[ivalue];
 803          }
 804       }
 805
 806       for  (Int_t ivalue = 0; ivalue < arrayEntries; ivalue++) {
 807          cstream << "calPads" <<
 808             (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue];
 809       }
 810
 811       if (mapFileName) {
 812          for  (Int_t ivalue = 0; ivalue < mapEntries; ivalue++) {
 813             if (isector < 36)
 814                cstream << "calPads" <<
 815                   (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapIROCArray[ivalue];
 816             else
 817                cstream << "calPads" <<
 818                   (Char_t*)((mapNames[ivalue] + ".=").Data()) << &mapOROCArray[ivalue];
 819          }
 820       }
 821
 822       cstream << "calPads" <<
 823          "row.=" << &posArray[0] <<
 824          "pad.=" << &posArray[1] <<
 825          "lx.=" << &posArray[2] <<
 826          "ly.=" << &posArray[3] <<
 827          "gx.=" << &posArray[4] <<
 828          "gy.=" << &posArray[5] <<
 829          "rpad.=" << &posArray[6] <<
 830          "channel.=" << &posArray[7];
 831
 832       cstream << "calPads" <<
 833          "\n";
 834
 835       delete[] posArray;
 836       delete[] vectorArray;
 837    }
 838    delete[] names;
 839    if (mapFileName) {
 840       delete mapIROCs;
 841       delete mapOROCs;
 842       delete[] mapIROCArray;
 843       delete[] mapOROCArray;
 844       delete[] mapNames;
 845    }
 846 }
 847
 848