TPC/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 "TLinearFitter.h"
  34 #include "TArrayI.h"
  35 //
  36 //
  37 #include "AliTPCCalROC.h"
  38 #include "AliMathBase.h"
  39
  40 ClassImp(AliTPCCalROC)
  41
  42
  43 //_____________________________________________________________________________
  44 AliTPCCalROC::AliTPCCalROC()
  45              :TObject(),
  46               fSector(0),
  47               fNChannels(0),
  48               fNRows(0),
  49               fIndexes(0),
  50               fData(0)
  51 {
  52   //
  53   // Default constructor
  54   //
  55
  56 }
  57
  58 //_____________________________________________________________________________
  59 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
  60              :TObject(),
  61               fSector(0),
  62               fNChannels(0),
  63               fNRows(0),
  64               fIndexes(0),
  65               fData(0)
  66 {
  67   //
  68   // Constructor that initializes a given sector
  69   //
  70   fSector = sector;
  71   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
  72   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
  73   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
  74   fData = new Float_t[fNChannels];
  75   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
  80              :TObject(c),
  81               fSector(0),
  82               fNChannels(0),
  83               fNRows(0),
  84               fIndexes(0),
  85               fData(0)
  86 {
  87   //
  88   // AliTPCCalROC copy constructor
  89   //
  90   fSector = c.fSector;
  91   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
  92   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
  93   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
  94   //
  95   fData   = new Float_t[fNChannels];
  96   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
  97 }
  98 //____________________________________________________________________________
  99 AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
 100 {
 101   //
 102   // assignment operator - dummy
 103   //
 104   fData=param.fData;
 105   return (*this);
 106 }
 107
 108
 109 //_____________________________________________________________________________
 110 AliTPCCalROC::~AliTPCCalROC()
 111 {
 112   //
 113   // AliTPCCalROC destructor
 114   //
 115   if (fData) {
 116     delete [] fData;
 117     fData = 0;
 118   }
 119 }
 120
 121
 122
 123 void AliTPCCalROC::Streamer(TBuffer &R__b)
 124 {
 125    // Stream an object of class AliTPCCalROC.
 126    if (R__b.IsReading()) {
 127       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
 128       fIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
 129    } else {
 130       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
 131    }
 132 }
 133
 134 //  //
 135 //   // algebra
 136 //   void Add(Float_t c1);
 137 //   void Multiply(Float_t c1);
 138 //   void Add(const AliTPCCalROC * roc, Double_t c1 = 1);
 139 //   void Divide(const AliTPCCalROC * roc);
 140
 141 void AliTPCCalROC::Add(Float_t c1){
 142   //
 143   // add constant
 144   //
 145   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
 146 }
 147 void AliTPCCalROC::Multiply(Float_t c1){
 148   //
 149   // add constant
 150   //
 151   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
 152 }
 153
 154 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
 155   //
 156   // add values
 157   //
 158   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 159     fData[idata]+=roc->fData[idata]*c1;
 160   }
 161 }
 162
 163
 164 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
 165   //
 166   // multiply values - per by pad
 167   //
 168   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 169     fData[idata]*=roc->fData[idata];
 170   }
 171 }
 172
 173
 174 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
 175   //
 176   // divide values
 177   //
 178   Float_t kEpsilon=0.00000000000000001;
 179   for (UInt_t  idata = 0; idata< fNChannels; idata++){
 180     if (TMath::Abs(roc->fData[idata])>kEpsilon)
 181       fData[idata]/=roc->fData[idata];
 182   }
 183 }
 184
 185
 186
 187
 188 Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
 189   //
 190   //  Calculate LTM mean and sigma
 191   //
 192   Double_t *ddata = new Double_t[fNChannels];
 193   Double_t mean=0, lsigma=0;
 194   Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels));
 195   for (UInt_t i=0;i<fNChannels;i++) ddata[i]= fData[i];
 196   AliMathBase::EvaluateUni(UInt_t(fNChannels),ddata, mean, lsigma, hh);
 197   if (sigma) *sigma=lsigma;
 198   delete [] ddata;
 199   return mean;
 200 }
 201
 202 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
 203   //
 204   // make 1D histo
 205   // type -1 = user defined range
 206   //       0 = nsigma cut nsigma=min
 207   //       1 = delta cut around median delta=min
 208   if (type>=0){
 209     if (type==0){
 210       // nsigma range
 211       Float_t mean  = GetMean();
 212       Float_t sigma = GetRMS();
 213       Float_t nsigma = TMath::Abs(min);
 214       min = mean-nsigma*sigma;
 215       max = mean+nsigma*sigma;
 216     }
 217     if (type==1){
 218       // fixed range
 219       Float_t mean   = GetMedian();
 220       Float_t  delta = min;
 221       min = mean-delta;
 222       max = mean+delta;
 223     }
 224     if (type==2){
 225       //
 226       // LTM mean +- nsigma
 227       //
 228       Double_t sigma;
 229       Float_t mean  = GetLTM(&sigma,max);
 230       sigma*=min;
 231       min = mean-sigma;
 232       max = mean+sigma;
 233     }
 234   }
 235   char  name[1000];
 236   sprintf(name,"%s ROC 1D%d",GetTitle(),fSector);
 237   TH1F * his = new TH1F(name,name,100, min,max);
 238   for (UInt_t irow=0; irow<fNRows; irow++){
 239     UInt_t npads = (Int_t)GetNPads(irow);
 240     for (UInt_t ipad=0; ipad<npads; ipad++){
 241       his->Fill(GetValue(irow,ipad));
 242     }
 243   }
 244   return his;
 245 }
 246
 247
 248
 249 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
 250   //
 251   // make 2D histo
 252   // type -1 = user defined range
 253   //       0 = nsigma cut nsigma=min
 254   //       1 = delta cut around median delta=min
 255   if (type>=0){
 256     if (type==0){
 257       // nsigma range
 258       Float_t mean  = GetMean();
 259       Float_t sigma = GetRMS();
 260       Float_t nsigma = TMath::Abs(min);
 261       min = mean-nsigma*sigma;
 262       max = mean+nsigma*sigma;
 263     }
 264     if (type==1){
 265       // fixed range
 266       Float_t mean   = GetMedian();
 267       Float_t  delta = min;
 268       min = mean-delta;
 269       max = mean+delta;
 270     }
 271     if (type==2){
 272       Double_t sigma;
 273       Float_t mean  = GetLTM(&sigma,max);
 274       sigma*=min;
 275       min = mean-sigma;
 276       max = mean+sigma;
 277
 278     }
 279   }
 280   UInt_t maxPad = 0;
 281   for (UInt_t irow=0; irow<fNRows; irow++){
 282     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
 283   }
 284   char  name[1000];
 285   sprintf(name,"%s ROC%d",GetTitle(),fSector);
 286   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
 287   for (UInt_t irow=0; irow<fNRows; irow++){
 288     UInt_t npads = (Int_t)GetNPads(irow);
 289     for (UInt_t ipad=0; ipad<npads; ipad++){
 290       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
 291     }
 292   }
 293   his->SetMaximum(max);
 294   his->SetMinimum(min);
 295   return his;
 296 }
 297
 298 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
 299   //
 300   // Make Histogram with outliers
 301   // mode = 0 - sigma cut used
 302   // mode = 1 - absolute cut used
 303   // fraction - fraction of values used to define sigma
 304   // delta - in mode 0 - nsigma cut
 305   //            mode 1 - delta cut
 306   Double_t sigma;
 307   Float_t mean  = GetLTM(&sigma,fraction);
 308   if (type==0) delta*=sigma;
 309   UInt_t maxPad = 0;
 310   for (UInt_t irow=0; irow<fNRows; irow++){
 311     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
 312   }
 313
 314   char  name[1000];
 315   sprintf(name,"%s ROC Outliers%d",GetTitle(),fSector);
 316   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
 317   for (UInt_t irow=0; irow<fNRows; irow++){
 318     UInt_t npads = (Int_t)GetNPads(irow);
 319     for (UInt_t ipad=0; ipad<npads; ipad++){
 320       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
 321         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
 322     }
 323   }
 324   return his;
 325 }
 326
 327
 328
 329 void AliTPCCalROC::Draw(Option_t* opt){
 330   //
 331   // create histogram with values and draw it
 332   //
 333   TH1 * his=0;
 334   TString option=opt;
 335   option.ToUpper();
 336   if (option.Contains("1D")){
 337     his = MakeHisto1D();
 338   }
 339   else{
 340     his = MakeHisto2D();
 341   }
 342   his->Draw(option);
 343 }
 344
 345
 346
 347
 348
 349 void AliTPCCalROC::Test(){
 350   //
 351   // example function to show functionality and tes AliTPCCalROC
 352   //
 353   AliTPCCalROC  roc0(0);
 354   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
 355     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
 356       Float_t value  = irow+ipad/1000.;
 357       roc0.SetValue(irow,ipad,value);
 358     }
 359   }
 360   //
 361   AliTPCCalROC roc1(roc0);
 362   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
 363     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
 364       Float_t value  = irow+ipad/1000.;
 365       if (roc1.GetValue(irow,ipad)!=value){
 366         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
 367       }
 368     }
 369   }
 370   TFile f("calcTest.root","recreate");
 371   roc0.Write("Roc0");
 372   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
 373   f.Close();
 374   //
 375   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
 376     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
 377       printf("NPads - Read/Write error\trow=%d\n",irow);
 378     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
 379       Float_t value  = irow+ipad/1000.;
 380       if (roc2->GetValue(irow,ipad)!=value){
 381         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
 382       }
 383     }
 384   }
 385 }
 386
 387
 388
 389 AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
 390   //
 391   // MakeLocalFit - smoothing
 392   // rowRadius  -  radius - rows to be used for smoothing
 393   // padradius  -  radius - pads to be used for smoothing
 394   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 395   // robust     -  robust method of fitting  - (much slower)
 396
 397   AliTPCCalROC * ROCfitted = new AliTPCCalROC(fSector);
 398   TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
 399   fitterQ.StoreData(kTRUE);
 400   for (Int_t row=0; row < GetNrows(); row++) {
 401     //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
 402     for (Int_t pad=0; pad < GetNPads(row); pad++)
 403       ROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust));
 404   }
 405   return ROCfitted;
 406 }
 407
 408
 409 Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
 410   //
 411   //  AliTPCCalROC::GetNeighbourhoodValue - smoothing (PRIVATE)
 412   // rowRadius  -  radius - rows to be used for smoothing
 413   // padradius  -  radius - pads to be used for smoothing
 414   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
 415   // robust     -  robust method of fitting  - (much slower)
 416
 417
 418
 419   fitterQ->ClearPoints();
 420   TVectorD fitParam(6);
 421   Int_t    npoints=0;
 422   Double_t xx[6];
 423   Float_t dlx, dly;
 424   Float_t lPad[3] = {0};
 425   Float_t localXY[3] = {0};
 426
 427   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
 428   tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad);  // calculate position lPad by pad and row number
 429
 430   TArrayI *neighbourhoodRows = 0;
 431   TArrayI *neighbourhoodPads = 0;
 432   GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
 433
 434   Int_t r, p;
 435   for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
 436     r = neighbourhoodRows->At(i);
 437     p = neighbourhoodPads->At(i);
 438     if (r == -1 || p == -1) continue;
 439     tpcROCinstance->GetPositionLocal(fSector, r, p, localXY);   // calculate position localXY by pad and row number
 440     dlx = lPad[0] - localXY[0];
 441     dly = lPad[1] - localXY[1];
 442     xx[0] = 1;
 443     xx[1] = dlx;
 444     xx[2] = dly;
 445     xx[3] = dlx*dlx;
 446     xx[4] = dly*dly;
 447     xx[5] = dlx*dly;
 448     if (ROCoutliers && ROCoutliers->GetValue(r,p) != 1) {
 449       fitterQ->AddPoint(xx, GetValue(r, p), 1);
 450       npoints++;
 451     }
 452   }
 453   if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
 454     // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
 455     return 0.;  // for diagnostic
 456   }
 457   fitterQ->Eval();
 458   fitterQ->GetParameters(fitParam);
 459   Float_t chi2Q = 0;
 460   chi2Q = fitterQ->GetChisquare()/(npoints-6.);
 461   //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
 462   if (robust && chi2Q > 5) {
 463     //std::cout << "robust fitter called... " << std::endl;
 464     fitterQ->EvalRobust(0.7);
 465     fitterQ->GetParameters(fitParam);
 466   }
 467   Double_t value = fitParam[0];
 468
 469   delete neighbourhoodRows;
 470   delete neighbourhoodPads;
 471
 472   //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q <<  std::endl;
 473
 474   return value;
 475 }
 476
 477
 478
 479
 480 void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
 481   //
 482   //
 483   //
 484   rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
 485   padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
 486   Int_t* rowArrayTemp = rowArray->GetArray();
 487   Int_t* padArrayTemp = padArray->GetArray();
 488
 489   Int_t rmin = row - rRadius;
 490   Int_t rmax = row + rRadius;
 491
 492   // if window goes out of ROC
 493   if (rmin < 0) {
 494     rmax = rmax - rmin;
 495     rmin = 0;
 496   }
 497   if (rmax >= GetNrows()) {
 498     rmin = rmin - (rmax - GetNrows()+1);
 499     rmax = GetNrows() - 1;
 500       if (rmin  < 0 ) rmin = 0; // if the window is bigger than the ROC
 501   }
 502
 503   Int_t pmin, pmax;
 504   Int_t i = 0;
 505
 506   for (Int_t r = rmin; r <= rmax; r++) {
 507     pmin = pad - pRadius;
 508     pmax = pad + pRadius;
 509     if (pmin < 0) {
 510       pmax = pmax - pmin;
 511       pmin = 0;
 512     }
 513     if (pmax >= GetNPads(r)) {
 514       pmin = pmin - (pmax - GetNPads(r)+1);
 515       pmax = GetNPads(r) - 1;
 516       if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
 517     }
 518     for (Int_t p = pmin; p <= pmax; p++) {
 519       rowArrayTemp[i] = r;
 520       padArrayTemp[i] = p;
 521       i++;
 522     }
 523   }
 524   for (Int_t j = i; j < rowArray->GetSize(); j++){  // unused padArray-entries, in the case that the window is bigger than the ROC
 525     //std::cout << "trying to write -1" << std::endl;
 526     rowArrayTemp[j] = -1;
 527     padArrayTemp[j] = -1;
 528     //std::cout << "writing -1" << std::endl;
 529   }
 530 }
 531
 532