TPC/AliTPCCalROC.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
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //     Calibration base class for a single ROC                               //
  20 //     Contains one float value per pad                                      //
  21 //     mapping of the pads taken form AliTPCROC                              //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 //
  26 // ROOT includes
  27 //
  28 #include "TMath.h"
  29 #include "TClass.h"
  30 #include "TFile.h"
  31 #include "TH1F.h"
  32 #include "TH2F.h"
  33 #include "TArrayI.h"
  34
  35 //
  36 //
  37 #include "AliTPCCalROC.h"
  38 #include "AliMathBase.h"
  39
  40 #include "TRandom3.h"      // only needed by the AliTPCCalROCTest() method
  41
  42 ClassImp(AliTPCCalROC)
  43
  44
  45 //_____________________________________________________________________________
  46 AliTPCCalROC::AliTPCCalROC()
  47              :TNamed(),
  48               fSector(0),
  49               fNChannels(0),
  50               fNRows(0),
  51               fkIndexes(0),
  52               fData(0)
  53 {
  54   //
  55   // Default constructor
  56   //
  57
  58 }
  59
  60 //_____________________________________________________________________________
  61 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
  62              :TNamed(),
  63               fSector(0),
  64               fNChannels(0),
  65               fNRows(0),
  66               fkIndexes(0),
  67               fData(0)
  68 {
  69   //
  70   // Constructor that initializes a given sector
  71   //
  72   fSector = sector;
  73   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
  74   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
  75   fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
  76   fData = new Float_t[fNChannels];
  77   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
  78 }
  79
  80 //_____________________________________________________________________________
  81 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
  82              :TNamed(c),
  83               fSector(0),
  84               fNChannels(0),
  85               fNRows(0),
  86               fkIndexes(0),
  87               fData(0)
  88 {
  89   //
  90   // AliTPCCalROC copy constructor
  91   //
  92   fSector = c.fSector;
  93   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
  94   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
  95   fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
  96   //
  97   fData   = new Float_t[fNChannels];
  98   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
  99 }
 100 //____________________________________________________________________________
 101 AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
 102 {
 103   //
 104   // assignment operator - dummy
 105   //
 106   if (this == &param) return (*this);
 107   fData=param.fData;
 108   return (*this);
 109 }
 110
 111
 112 //_____________________________________________________________________________
 113 AliTPCCalROC::~AliTPCCalROC()
 114 {
 115   //
 116   // AliTPCCalROC destructor
 117   //
 118   if (fData) {
 119     delete [] fData;
 120     fData = 0;
 121   }
 122 }
 123
 124
 125
 126 void AliTPCCalROC::Streamer(TBuffer &R__b)
 127 {
 128    //
 129    // Stream an object of class AliTPCCalROC.
 130    //
 131    if (R__b.IsReading()) {
 132       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
 133       fkIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
 134    } else {
 135       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
 136    }
 137 }
 138
 139
 140 // algebra fuctions:
 141
 142 void AliTPCCalROC::Add(Float_t c1){
 143   //
 144   // add c1 to each channel of the ROC
 145   //
 146   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
 147 }
 148
 149
 150 void AliTPCCalROC::Multiply(Float_t c1){
 151   //
 152   // multiply each channel of the ROC with c1
 153   //
 154   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
 155 }
 156
 157
 158 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
 159   //
 160   // multiply AliTPCCalROC roc by c1 and add each channel to the coresponing channel in the ROC
 161   //  - pad by pad
 162   //
 163   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 164     fData[idata]+=roc->fData[idata]*c1;
 165   }
 166 }
 167
 168
 169 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
 170   //
 171   // multiply each channel of the ROC with the coresponding channel of 'roc'
 172   //     - pad by pad -
 173   //
 174   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 175     fData[idata]*=roc->fData[idata];
 176   }
 177 }
 178
 179
 180 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
 181   //
 182   // divide each channel of the ROC by the coresponding channel of 'roc'
 183   //     - pad by pad -
 184   //
 185   Float_t kEpsilon=0.00000000000000001;
 186   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 187     if (TMath::Abs(roc->fData[idata])>kEpsilon)
 188       fData[idata]/=roc->fData[idata];
 189   }
 190 }
 191
 192 Double_t AliTPCCalROC::GetMean(AliTPCCalROC *const outlierROC) const {
 193    //
 194    //  returns the mean value of the ROC
 195    //  pads with value != 0 in outlierROC are not used for the calculation
 196    //  return 0 if no data is accepted by the outlier cuts
 197    //
 198    if (!outlierROC) return TMath::Mean(fNChannels, fData);
 199    Double_t *ddata = new Double_t[fNChannels];
 200    Int_t nPoints = 0;
 201    for (UInt_t i=0;i<fNChannels;i++) {
 202       if (!(outlierROC->GetValue(i))) {
 203          ddata[nPoints]= fData[i];
 204          nPoints++;
 205       }
 206    }
 207    Double_t mean = 0;
 208    if(nPoints>0)
 209      mean = TMath::Mean(nPoints, ddata);
 210    delete [] ddata;
 211    return mean;
 212 }
 213
 214 Double_t AliTPCCalROC::GetMedian(AliTPCCalROC *const outlierROC) const {
 215    //
 216    //  returns the median value of the ROC
 217    //  pads with value != 0 in outlierROC are not used for the calculation
 218    //  return 0 if no data is accepted by the outlier cuts
 219    //
 220    if (!outlierROC) return TMath::Median(fNChannels, fData);
 221    Double_t *ddata = new Double_t[fNChannels];
 222    Int_t nPoints = 0;
 223    for (UInt_t i=0;i<fNChannels;i++) {
 224       if (!(outlierROC->GetValue(i))) {
 225          ddata[nPoints]= fData[i];
 226          nPoints++;
 227       }
 228    }
 229    Double_t median = 0;
 230    if(nPoints>0)
 231      median = TMath::Median(nPoints, ddata);
 232    delete [] ddata;
 233    return median;
 234 }
 235
 236 Double_t AliTPCCalROC::GetRMS(AliTPCCalROC *const outlierROC) const {
 237    //
 238    //  returns the RMS value of the ROC
 239    //  pads with value != 0 in outlierROC are not used for the calculation
 240    //  return 0 if no data is accepted by the outlier cuts
 241    //
 242    if (!outlierROC) return TMath::RMS(fNChannels, fData);
 243    Double_t *ddata = new Double_t[fNChannels];
 244    Int_t nPoints = 0;
 245    for (UInt_t i=0;i<fNChannels;i++) {
 246       if (!(outlierROC->GetValue(i))) {
 247          ddata[nPoints]= fData[i];
 248          nPoints++;
 249       }
 250    }
 251    Double_t rms = 0;
 252    if(nPoints>0)
 253      rms = TMath::RMS(nPoints, ddata);
 254    delete [] ddata;
 255    return rms;
 256 }
 257
 258 Double_t AliTPCCalROC::GetLTM(Double_t *const sigma, Double_t fraction, AliTPCCalROC *const outlierROC){
 259   //
 260   //  returns the LTM and sigma
 261   //  pads with value != 0 in outlierROC are not used for the calculation
 262    //  return 0 if no data is accepted by the outlier cuts
 263   //
 264   Double_t *ddata = new Double_t[fNChannels];
 265   UInt_t nPoints = 0;
 266   for (UInt_t i=0;i<fNChannels;i++) {
 267      if (!outlierROC || !(outlierROC->GetValue(i))) {
 268         ddata[nPoints]= fData[i];
 269         nPoints++;
 270      }
 271   }
 272
 273   Double_t ltm =0, lsigma=0;
 274   if(nPoints>0) {
 275     Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
 276     AliMathBase::EvaluateUni(nPoints,ddata, ltm, lsigma, hh);
 277     if (sigma) *sigma=lsigma;
 278   }
 279
 280   delete [] ddata;
 281   return ltm;
 282 }
 283
 284 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
 285   //
 286   // make 1D histo
 287   // type -1 = user defined range
 288   //       0 = nsigma cut nsigma=min
 289   //       1 = delta cut around median delta=min
 290   //
 291   if (type>=0){
 292     if (type==0){
 293       // nsigma range
 294       Float_t mean  = GetMean();
 295       Float_t sigma = GetRMS();
 296       Float_t nsigma = TMath::Abs(min);
 297       min = mean-nsigma*sigma;
 298       max = mean+nsigma*sigma;
 299     }
 300     if (type==1){
 301       // fixed range
 302       Float_t mean   = GetMedian();
 303       Float_t  delta = min;
 304       min = mean-delta;
 305       max = mean+delta;
 306     }
 307     if (type==2){
 308       //
 309       // LTM mean +- nsigma
 310       //
 311       Double_t sigma;
 312       Float_t mean  = GetLTM(&sigma,max);
 313       sigma*=min;
 314       min = mean-sigma;
 315       max = mean+sigma;
 316     }
 317   }
 318   TString name=Form("%s ROC 1D%d",GetTitle(),fSector);
 319   TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
 320   for (UInt_t irow=0; irow<fNRows; irow++){
 321     UInt_t npads = (Int_t)GetNPads(irow);
 322     for (UInt_t ipad=0; ipad<npads; ipad++){
 323       his->Fill(GetValue(irow,ipad));
 324     }
 325   }
 326   return his;
 327 }
 328
 329
 330
 331 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
 332   //
 333   // make 2D histo
 334   // type -1 = user defined range
 335   //       0 = nsigma cut nsigma=min
 336   //       1 = delta cut around median delta=min
 337   //
 338   if (type>=0){
 339     if (type==0){
 340       // nsigma range
 341       Float_t mean  = GetMean();
 342       Float_t sigma = GetRMS();
 343       Float_t nsigma = TMath::Abs(min);
 344       min = mean-nsigma*sigma;
 345       max = mean+nsigma*sigma;
 346     }
 347     if (type==1){
 348       // fixed range
 349       Float_t mean   = GetMedian();
 350       Float_t  delta = min;
 351       min = mean-delta;
 352       max = mean+delta;
 353     }
 354     if (type==2){
 355       Double_t sigma;
 356       Float_t mean  = GetLTM(&sigma,max);
 357       sigma*=min;
 358       min = mean-sigma;
 359       max = mean+sigma;
 360
 361     }
 362   }
 363   UInt_t maxPad = 0;
 364   for (UInt_t irow=0; irow<fNRows; irow++){
 365     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
 366   }
 367
 368   TString name=Form("%s ROC%d",GetTitle(),fSector);
 369   TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
 370   for (UInt_t irow=0; irow<fNRows; irow++){
 371     UInt_t npads = (Int_t)GetNPads(irow);
 372     for (UInt_t ipad=0; ipad<npads; ipad++){
 373       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
 374     }
 375   }
 376   his->SetMaximum(max);
 377   his->SetMinimum(min);
 378   return his;
 379 }
 380
 381 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
 382   //
 383   // Make Histogram with outliers
 384   // mode = 0 - sigma cut used
 385   // mode = 1 - absolute cut used
 386   // fraction - fraction of values used to define sigma
 387   // delta - in mode 0 - nsigma cut
 388   //            mode 1 - delta cut
 389   //
 390   Double_t sigma;
 391   Float_t mean  = GetLTM(&sigma,fraction);
 392   if (type==0) delta*=sigma;
 393   UInt_t maxPad = 0;
 394   for (UInt_t irow=0; irow<fNRows; irow++){
 395     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
 396   }
 397
 398   TString name=Form("%s ROC Outliers%d",GetTitle(),fSector);
 399   TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
 400   for (UInt_t irow=0; irow<fNRows; irow++){
 401     UInt_t npads = (Int_t)GetNPads(irow);
 402     for (UInt_t ipad=0; ipad<npads; ipad++){
 403       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
 404         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
 405     }
 406   }
 407   return his;
 408 }
 409
 410
 411
 412 void AliTPCCalROC::Draw(Option_t* opt){
 413   //
 414   // create histogram with values and draw it
 415   //
 416   TH1 * his=0;
 417   TString option=opt;
 418   option.ToUpper();
 419   if (option.Contains("1D")){
 420     his = MakeHisto1D();
 421   }
 422   else{
 423     his = MakeHisto2D();
 424   }
 425   his->Draw(option);
 426 }
 427
 428
 429
 430
 431
 432 void AliTPCCalROC::Test() {
 433   //
 434   // example function to show functionality and test AliTPCCalROC
 435   //
 436
 437   Float_t kEpsilon=0.00001;
 438
 439   // create CalROC with defined values
 440   AliTPCCalROC  roc0(0);
 441   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
 442     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
 443       Float_t value  = irow+ipad/1000.;
 444       roc0.SetValue(irow,ipad,value);
 445     }
 446   }
 447
 448   // copy CalROC, readout values and compare to original
 449   AliTPCCalROC roc1(roc0);
 450   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
 451     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
 452       Float_t value  = irow+ipad/1000.;
 453       if (roc1.GetValue(irow,ipad)!=value){
 454         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
 455       }
 456     }
 457   }
 458
 459   // write original CalROC to file
 460   TFile f("calcTest.root","recreate");
 461   roc0.Write("Roc0");
 462   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
 463   f.Close();
 464
 465   // read CalROC from file and compare to original CalROC
 466   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
 467     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
 468       printf("NPads - Read/Write error\trow=%d\n",irow);
 469     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
 470       Float_t value  = irow+ipad/1000.;
 471       if (roc2->GetValue(irow,ipad)!=value){
 472         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
 473       }
 474     }
 475   }
 476
 477   //
 478   // Algebra Tests
 479   //
 480
 481   // Add constant
 482   AliTPCCalROC roc3(roc0);
 483   roc3.Add(1.5);
 484   for (UInt_t irow = 0; irow <roc3.GetNrows(); irow++){
 485     for (UInt_t ipad = 0; ipad <roc3.GetNPads(irow); ipad++){
 486       Float_t value  = irow+ipad/1000. + 1.5;
 487       if (TMath::Abs(roc3.GetValue(irow,ipad)-value) > kEpsilon){
 488         printf("Add constant - error\trow=%d\tpad=%d\n",irow,ipad);
 489       }
 490     }
 491   }
 492
 493   // Add another CalROC
 494   AliTPCCalROC roc4(roc0);
 495   roc4.Add(&roc0, -1.5);
 496   for (UInt_t irow = 0; irow <roc4.GetNrows(); irow++){
 497     for (UInt_t ipad = 0; ipad <roc4.GetNPads(irow); ipad++){
 498       Float_t value  = irow+ipad/1000. - 1.5 * (irow+ipad/1000.);
 499       if (TMath::Abs(roc4.GetValue(irow,ipad)-value) > kEpsilon){
 500         printf("Add CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
 501       }
 502     }
 503   }
 504
 505   // Multiply with constant
 506   AliTPCCalROC roc5(roc0);
 507   roc5.Multiply(-1.4);
 508   for (UInt_t irow = 0; irow <roc5.GetNrows(); irow++){
 509     for (UInt_t ipad = 0; ipad <roc5.GetNPads(irow); ipad++){
 510       Float_t value  = (irow+ipad/1000.) * (-1.4);
 511       if (TMath::Abs(roc5.GetValue(irow,ipad)-value) > kEpsilon){
 512         printf("Multiply with constant - error\trow=%d\tpad=%d\n",irow,ipad);
 513       }
 514     }
 515   }
 516
 517   // Multiply another CalROC
 518   AliTPCCalROC roc6(roc0);
 519   roc6.Multiply(&roc0);
 520   for (UInt_t irow = 0; irow <roc6.GetNrows(); irow++){
 521     for (UInt_t ipad = 0; ipad <roc6.GetNPads(irow); ipad++){
 522       Float_t value  = (irow+ipad/1000.) * (irow+ipad/1000.);
 523       if (TMath::Abs(roc6.GetValue(irow,ipad)-value) > kEpsilon*100){
 524         printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
 525       }
 526     }
 527   }
 528
 529
 530   // Divide by CalROC
 531   AliTPCCalROC roc7(roc0);
 532   roc7.Divide(&roc0);
 533   for (UInt_t irow = 0; irow <roc7.GetNrows(); irow++){
 534     for (UInt_t ipad = 0; ipad <roc7.GetNPads(irow); ipad++){
 535       Float_t value  = 1.;
 536       if (irow+ipad == 0) value = roc0.GetValue(irow,ipad);
 537       if (TMath::Abs(roc7.GetValue(irow,ipad)-value) > kEpsilon){
 538         printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
 539       }
 540     }
 541   }
 542
 543   //
 544   // Statistics Test
 545   //
 546
 547   // create CalROC with defined values
 548   TRandom3 rnd(0);
 549   AliTPCCalROC sroc0(0);
 550   for (UInt_t ichannel = 0; ichannel < sroc0.GetNchannels(); ichannel++){
 551     Float_t value  = rnd.Gaus(10., 2.);
 552     sroc0.SetValue(ichannel,value);
 553   }
 554
 555   printf("Mean   (should be close to 10): %f\n", sroc0.GetMean());
 556   printf("RMS    (should be close to  2): %f\n", sroc0.GetRMS());
 557   printf("Median (should be close to 10): %f\n", sroc0.GetMedian());
 558   printf("LTM    (should be close to 10): %f\n", sroc0.GetLTM());
 559
 560   //AliTPCCalROC* sroc1 = sroc0.LocalFit(4, 8);
 561
 562   //delete sroc1;
 563
 564 //        std::cout << TMath::Abs(roc5.GetValue(irow,ipad)-value) << std::endl;
 565 }
 566
 567
 568 AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
 569   //
 570   // MakeLocalFit - smoothing
 571   // returns a AliTPCCalROC with smoothed data
 572   // rowRadius and padRadius specifies a window around a given pad.
 573   // The data of this window are fitted with a parabolic function.
 574   // This function is evaluated at the pad's position.
 575   // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
 576   // rowRadius  -  radius - rows to be used for smoothing
 577   // padradius  -  radius - pads to be used for smoothing
 578   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 579   // robust     -  robust method of fitting  - (much slower)
 580   // chi2Threshold: Threshold for chi2 when EvalRobust is called
 581   // robustFraction: Fraction of data that will be used in EvalRobust
 582   //
 583   AliTPCCalROC * xROCfitted = new AliTPCCalROC(fSector);
 584   TLinearFitter fitterQ(6,"hyp5");
 585    // TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
 586   fitterQ.StoreData(kTRUE);
 587   for (UInt_t row=0; row < GetNrows(); row++) {
 588     //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
 589     for (UInt_t pad=0; pad < GetNPads(row); pad++)
 590       xROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
 591   }
 592   return xROCfitted;
 593 }
 594
 595
 596 Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC *const ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
 597   //
 598   // AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
 599   // in this function the fit for LocalFit is done
 600   //
 601
 602   fitterQ->ClearPoints();
 603   TVectorD fitParam(6);
 604   Int_t    npoints=0;
 605   Double_t xx[6];
 606   Float_t dlx, dly;
 607   Float_t lPad[3] = {0};
 608   Float_t localXY[3] = {0};
 609
 610   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 611   tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad);  // calculate position lPad by pad and row number
 612
 613   TArrayI *neighbourhoodRows = 0;
 614   TArrayI *neighbourhoodPads = 0;
 615
 616   //std::cerr << "Trying to get neighbourhood for row " << row << ", pad " << pad << std::endl;
 617   GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
 618   //std::cerr << "Got neighbourhood for row " << row << ", pad " << pad << std::endl;
 619
 620   Int_t r, p;
 621   for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
 622     r = neighbourhoodRows->At(i);
 623     p = neighbourhoodPads->At(i);
 624     if (r == -1 || p == -1) continue;   // window is bigger than ROC
 625     tpcROCinstance->GetPositionLocal(fSector, r, p, localXY);   // calculate position localXY by pad and row number
 626     dlx = lPad[0] - localXY[0];
 627     dly = lPad[1] - localXY[1];
 628     //xx[0] = 1;
 629     xx[1] = dlx;
 630     xx[2] = dly;
 631     xx[3] = dlx*dlx;
 632     xx[4] = dly*dly;
 633     xx[5] = dlx*dly;
 634     if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
 635       fitterQ->AddPoint(&xx[1], GetValue(r, p), 1);
 636       npoints++;
 637     }
 638   }
 639
 640   delete neighbourhoodRows;
 641   delete neighbourhoodPads;
 642
 643   if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
 644     // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
 645     return 0.;  // for diagnostic
 646   }
 647   fitterQ->Eval();
 648   fitterQ->GetParameters(fitParam);
 649   Float_t chi2Q = 0;
 650   if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
 651   //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
 652   if (robust && chi2Q > chi2Threshold) {
 653     //std::cout << "robust fitter called... " << std::endl;
 654     fitterQ->EvalRobust(robustFraction);
 655     fitterQ->GetParameters(fitParam);
 656   }
 657   Double_t value = fitParam[0];
 658
 659   //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q <<  std::endl;
 660   return value;
 661 }
 662
 663
 664
 665
 666 void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
 667   //
 668   // AliTPCCalROC::GetNeighbourhood - PRIVATE
 669   // in this function the window for LocalFit is determined
 670   //
 671   rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
 672   padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
 673
 674   Int_t rmin = row - rRadius;
 675   UInt_t rmax = row + rRadius;
 676
 677   // if window goes out of ROC
 678   if (rmin < 0) {
 679     rmax = rmax - rmin;
 680     rmin = 0;
 681   }
 682   if (rmax >= GetNrows()) {
 683     rmin = rmin - (rmax - GetNrows()+1);
 684     rmax = GetNrows() - 1;
 685       if (rmin  < 0 ) rmin = 0; // if the window is bigger than the ROC
 686   }
 687
 688   Int_t pmin, pmax;
 689   Int_t i = 0;
 690
 691   for (UInt_t r = rmin; r <= rmax; r++) {
 692     pmin = pad - pRadius;
 693     pmax = pad + pRadius;
 694     if (pmin < 0) {
 695       pmax = pmax - pmin;
 696       pmin = 0;
 697     }
 698     if (pmax >= (Int_t)GetNPads(r)) {
 699       pmin = pmin - (pmax - GetNPads(r)+1);
 700       pmax = GetNPads(r) - 1;
 701       if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
 702     }
 703     for (Int_t p = pmin; p <= pmax; p++) {
 704       (*rowArray)[i] = r;
 705       (*padArray)[i] = p;
 706       i++;
 707     }
 708   }
 709   for (Int_t j = i; j < rowArray->GetSize(); j++){  // unused padArray-entries, in the case that the window is bigger than the ROC
 710     //std::cout << "trying to write -1" << std::endl;
 711     (*rowArray)[j] = -1;
 712     (*padArray)[j] = -1;
 713     //std::cout << "writing -1" << std::endl;
 714   }
 715 }
 716
 717
 718
 719 void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err){
 720   //
 721   // Makes a  GlobalFit for the given secotr and return fit-parameters, covariance and chi2
 722   // The origin of the fit function is the center of the ROC!
 723   // fitType == 0: fit plane function
 724   // fitType == 1: fit parabolic function
 725   // ROCoutliers - pads with value !=0 are not used in fitting procedure
 726   // chi2Threshold: Threshold for chi2 when EvalRobust is called
 727   // robustFraction: Fraction of data that will be used in EvalRobust
 728   // err: error of the data points
 729   //
 730   TLinearFitter* fitterG = 0;
 731   Double_t xx[6];
 732
 733   if (fitType  == 1) {
 734     fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
 735     fitParam.ResizeTo(6);
 736     covMatrix.ResizeTo(6,6);
 737   } else {
 738     fitterG = new TLinearFitter(3,"x0++x1++x2");
 739     fitParam.ResizeTo(3);
 740     covMatrix.ResizeTo(3,3);
 741   }
 742   fitterG->StoreData(kTRUE);
 743   fitterG->ClearPoints();
 744   Int_t    npoints=0;
 745
 746   Float_t dlx, dly;
 747   Float_t centerPad[3] = {0};
 748   Float_t localXY[3] = {0};
 749
 750   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 751   tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad);  // calculate center of ROC
 752
 753   // loop over all channels and read data into fitterG
 754   for (UInt_t irow = 0; irow < GetNrows(); irow++) {
 755     for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
 756       // fill fitterG
 757       if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 0) continue;
 758       tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 759       dlx = localXY[0] - centerPad[0];
 760       dly = localXY[1] - centerPad[1];
 761       xx[0] = 1;
 762       xx[1] = dlx;
 763       xx[2] = dly;
 764       xx[3] = dlx*dlx;
 765       xx[4] = dly*dly;
 766       xx[5] = dlx*dly;
 767       npoints++;
 768       fitterG->AddPoint(xx, GetValue(irow, ipad), err);
 769     }
 770   }
 771   if(npoints>10) { // make sure there is something to fit
 772     fitterG->Eval();
 773     fitterG->GetParameters(fitParam);
 774     fitterG->GetCovarianceMatrix(covMatrix);
 775     if (fitType == 1)
 776       chi2 = fitterG->GetChisquare()/(npoints-6.);
 777     else chi2 = fitterG->GetChisquare()/(npoints-3.);
 778     if (robust && chi2 > chi2Threshold) {
 779       //    std::cout << "robust fitter called... " << std::endl;
 780       fitterG->EvalRobust(robustFraction);
 781       fitterG->GetParameters(fitParam);
 782     }
 783   } else {
 784     // set parameteres to 0
 785     Int_t nParameters = 3;
 786     if (fitType  == 1)
 787       nParameters = 6;
 788
 789     for(Int_t i = 0; i < nParameters; i++)
 790       fitParam[i] = 0;
 791   }
 792
 793   delete fitterG;
 794 }
 795
 796
 797 AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector){
 798   //
 799   // Create ROC with global fit parameters
 800   // The origin of the fit function is the center of the ROC!
 801   // loop over all channels, write fit values into new ROC and return it
 802   //
 803   Float_t dlx, dly;
 804   Float_t centerPad[3] = {0};
 805   Float_t localXY[3] = {0};
 806   AliTPCCalROC * xROCfitted = new AliTPCCalROC(sector);
 807   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 808   tpcROCinstance->GetPositionLocal(sector, xROCfitted->GetNrows()/2, xROCfitted->GetNPads(xROCfitted->GetNrows()/2)/2, centerPad);  // calculate center of ROC
 809   Int_t fitType = 1;
 810   if (fitParam.GetNoElements() == 6) fitType = 1;
 811   else fitType = 0;
 812   Double_t value = 0;
 813   if (fitType == 1) { // parabolic fit
 814     for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
 815       for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
 816         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 817         dlx = localXY[0] - centerPad[0];
 818         dly = localXY[1] - centerPad[1];
 819         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
 820         xROCfitted->SetValue(irow, ipad, value);
 821       }
 822     }
 823   }
 824   else {  // linear fit
 825     for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
 826       for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
 827         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
 828         dlx = localXY[0] - centerPad[0];
 829         dly = localXY[1] - centerPad[1];
 830         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
 831         xROCfitted->SetValue(irow, ipad, value);
 832       }
 833     }
 834   }
 835   return xROCfitted;
 836 }
 837