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