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