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