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