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