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