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