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