TRD/Cal/AliTRDCalROC.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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Calibration base class for a single ROC                                  //
  21 //  Contains one UShort_t value per pad                                      //
  22 //  However, values are set and get as float, there are stored internally as //
  23 //  (UShort_t) value * 10000                                                 //
  24 //                                                                           //
  25 ///////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TStyle.h>
  28 #include <TMath.h>
  29 #include <TH2F.h>
  30 #include <TH1F.h>
  31
  32 #include "AliMathBase.h"
  33
  34 #include "AliTRDCalROC.h"
  35
  36 ClassImp(AliTRDCalROC)
  37
  38 //_____________________________________________________________________________
  39 AliTRDCalROC::AliTRDCalROC()
  40   :TObject()
  41   ,fPla(0)
  42   ,fCha(0)
  43   ,fNrows(0)
  44   ,fNcols(0)
  45   ,fNchannels(0)
  46   ,fData(0)
  47 {
  48   //
  49   // Default constructor
  50   //
  51
  52 }
  53
  54 //_____________________________________________________________________________
  55 AliTRDCalROC::AliTRDCalROC(Int_t p, Int_t c)
  56   :TObject()
  57   ,fPla(p)
  58   ,fCha(c)
  59   ,fNrows(0)
  60   ,fNcols(144)
  61   ,fNchannels(0)
  62   ,fData(0)
  63 {
  64   //
  65   // Constructor that initializes a given pad plane type
  66   //
  67
  68   //
  69   // The pad plane parameter
  70   //
  71   switch (p) {
  72   case 0:
  73     if (c == 2) {
  74       // L0C0 type
  75       fNrows        =  12;
  76     }
  77     else {
  78       // L0C1 type
  79       fNrows        =  16;
  80     }
  81     break;
  82   case 1:
  83     if (c == 2) {
  84       // L1C0 type
  85       fNrows        =  12;
  86     }
  87     else {
  88       // L1C1 type
  89       fNrows        =  16;
  90     }
  91     break;
  92   case 2:
  93     if (c == 2) {
  94       // L2C0 type
  95       fNrows        =  12;
  96     }
  97     else {
  98       // L2C1 type
  99       fNrows        =  16;
 100     }
 101     break;
 102   case 3:
 103     if (c == 2) {
 104       // L3C0 type
 105       fNrows        =  12;
 106     }
 107     else {
 108       // L3C1 type
 109       fNrows        =  16;
 110     }
 111     break;
 112   case 4:
 113     if (c == 2) {
 114       // L4C0 type
 115       fNrows        =  12;
 116     }
 117     else {
 118       // L4C1 type
 119       fNrows        =  16;
 120     }
 121     break;
 122   case 5:
 123     if (c == 2) {
 124       // L5C0 type
 125       fNrows        =  12;
 126     }
 127     else {
 128       // L5C1 type
 129       fNrows        =  16;
 130     }
 131     break;
 132   };
 133
 134   fNchannels = fNrows * fNcols;
 135   if (fNchannels != 0) {
 136     fData = new UShort_t[fNchannels];
 137   }
 138
 139   for (Int_t i=0; i<fNchannels; ++i) {
 140     fData[i] = 0;
 141   }
 142
 143 }
 144
 145 //_____________________________________________________________________________
 146 AliTRDCalROC::AliTRDCalROC(const AliTRDCalROC &c)
 147   :TObject(c)
 148   ,fPla(c.fPla)
 149   ,fCha(c.fCha)
 150   ,fNrows(c.fNrows)
 151   ,fNcols(c.fNcols)
 152   ,fNchannels(c.fNchannels)
 153   ,fData(0)
 154 {
 155   //
 156   // AliTRDCalROC copy constructor
 157   //
 158
 159   Int_t iBin = 0;
 160
 161   fData = new UShort_t[fNchannels];
 162   for (iBin = 0; iBin < fNchannels; iBin++) {
 163     fData[iBin] = ((AliTRDCalROC &) c).fData[iBin];
 164   }
 165
 166 }
 167
 168 //_____________________________________________________________________________
 169 AliTRDCalROC::~AliTRDCalROC()
 170 {
 171   //
 172   // AliTRDCalROC destructor
 173   //
 174
 175   if (fData) {
 176     delete [] fData;
 177     fData = 0;
 178   }
 179
 180 }
 181
 182 //_____________________________________________________________________________
 183 AliTRDCalROC &AliTRDCalROC::operator=(const AliTRDCalROC &c)
 184 {
 185   //
 186   // Assignment operator
 187   //
 188
 189   if (this != &c) ((AliTRDCalROC &) c).Copy(*this);
 190   return *this;
 191
 192 }
 193
 194 //_____________________________________________________________________________
 195 void AliTRDCalROC::Copy(TObject &c) const
 196 {
 197   //
 198   // Copy function
 199   //
 200
 201   ((AliTRDCalROC &) c).fPla          = fPla;
 202   ((AliTRDCalROC &) c).fCha          = fCha;
 203
 204   ((AliTRDCalROC &) c).fNrows        = fNrows;
 205   ((AliTRDCalROC &) c).fNcols        = fNcols;
 206
 207   Int_t iBin = 0;
 208
 209   ((AliTRDCalROC &) c).fNchannels = fNchannels;
 210
 211   if (((AliTRDCalROC &) c).fData) delete [] ((AliTRDCalROC &) c).fData;
 212   ((AliTRDCalROC &) c).fData = new UShort_t[fNchannels];
 213   for (iBin = 0; iBin < fNchannels; iBin++) {
 214     ((AliTRDCalROC &) c).fData[iBin] = fData[iBin];
 215   }
 216
 217   TObject::Copy(c);
 218
 219 }
 220
 221 //___________________________________________________________________________________
 222 Double_t AliTRDCalROC::GetMean(AliTRDCalROC* const outlierROC) const
 223 {
 224   //
 225   // Calculate the mean
 226   //
 227
 228    Double_t *ddata = new Double_t[fNchannels];
 229    Int_t nPoints = 0;
 230    for (Int_t i=0;i<fNchannels;i++) {
 231       if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 232         if(fData[i] > 0.000000000000001){
 233          ddata[nPoints]= (Double_t) fData[i]/10000;
 234          nPoints++;
 235         }
 236       }
 237    }
 238    Double_t mean = TMath::Mean(nPoints,ddata);
 239    delete [] ddata;
 240    return mean;
 241 }
 242
 243 //_______________________________________________________________________________________
 244 Double_t AliTRDCalROC::GetMedian(AliTRDCalROC* const outlierROC) const
 245 {
 246   //
 247   // Calculate the median
 248   //
 249
 250   Double_t *ddata = new Double_t[fNchannels];
 251    Int_t nPoints = 0;
 252    for (Int_t i=0;i<fNchannels;i++) {
 253        if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 254          if(fData[i] > 0.000000000000001){
 255            ddata[nPoints]= (Double_t) fData[i]/10000;
 256            nPoints++;
 257          }
 258        }
 259    }
 260    Double_t mean = TMath::Median(nPoints,ddata);
 261    delete [] ddata;
 262    return mean;
 263 }
 264
 265 //____________________________________________________________________________________________
 266 Double_t AliTRDCalROC::GetRMS(AliTRDCalROC* const outlierROC) const
 267 {
 268   //
 269   // Calculate the RMS
 270   //
 271
 272   Double_t *ddata = new Double_t[fNchannels];
 273   Int_t nPoints = 0;
 274   for (Int_t i=0;i<fNchannels;i++) {
 275     if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 276        if(fData[i] > 0.000000000000001){
 277          ddata[nPoints]= (Double_t)fData[i]/10000;
 278          nPoints++;
 279        }
 280     }
 281   }
 282   Double_t mean = TMath::RMS(nPoints,ddata);
 283   delete [] ddata;
 284   return mean;
 285 }
 286
 287 //______________________________________________________________________________________________
 288 Double_t AliTRDCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTRDCalROC* const outlierROC)
 289 {
 290   //
 291   //  Calculate LTM mean and sigma
 292   //
 293
 294   Double_t *ddata = new Double_t[fNchannels];
 295   Double_t mean=0, lsigma=0;
 296   UInt_t nPoints = 0;
 297   for (Int_t i=0;i<fNchannels;i++) {
 298      if (!outlierROC || !(outlierROC->GetValue(i))) {
 299        if(fData[i] > 0.000000000000001){
 300          ddata[nPoints]= (Double_t) fData[i]/10000;
 301          nPoints++;
 302        }
 303      }
 304   }
 305   Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
 306   AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
 307   if (sigma) *sigma=lsigma;
 308   delete [] ddata;
 309   return mean;
 310 }
 311
 312 //___________________________________________________________________________________
 313 Bool_t AliTRDCalROC::Add(Float_t c1)
 314 {
 315   //
 316   // add constant
 317   //
 318
 319   Bool_t result = kTRUE;
 320   for (Int_t  idata = 0; idata< fNchannels; idata++) {
 321     if(((GetValue(idata)+c1) <= 6.5535) && ((GetValue(idata)+c1) >= 0.0)) SetValue(idata,GetValue(idata)+c1);
 322     else {
 323       SetValue(idata,0.0);
 324       result = kFALSE;
 325     }
 326   }
 327   return result;
 328 }
 329
 330 //_______________________________________________________________________________________
 331 Bool_t AliTRDCalROC::Multiply(Float_t c1)
 332 {
 333   //
 334   // multiply constant
 335   //
 336
 337   Bool_t result = kTRUE;
 338   if(c1 < 0) return kFALSE;
 339   for (Int_t  idata = 0; idata< fNchannels; idata++) {
 340     if((GetValue(idata)*c1) <= 6.5535) SetValue(idata,GetValue(idata)*c1);
 341     else {
 342       SetValue(idata,0.0);
 343       result = kFALSE;
 344     }
 345   }
 346   return result;
 347 }
 348
 349 //____________________________________________________________________________________________
 350 Bool_t AliTRDCalROC::Add(const AliTRDCalROC * roc, Double_t c1)
 351 {
 352   //
 353   // add values
 354   //
 355
 356   Bool_t result = kTRUE;
 357   for (Int_t  idata = 0; idata< fNchannels; idata++){
 358     if(((GetValue(idata)+roc->GetValue(idata)*c1) <= 6.5535) &&
 359        ((GetValue(idata)+roc->GetValue(idata)*c1) >= 0.0))
 360           SetValue(idata,GetValue(idata)+roc->GetValue(idata)*c1);
 361     else {
 362       SetValue(idata,0.0);
 363       result = kFALSE;
 364     }
 365   }
 366   return result;
 367 }
 368
 369 //____________________________________________________________________________________________
 370 Bool_t AliTRDCalROC::Multiply(const AliTRDCalROC*  roc)
 371 {
 372   //
 373   // multiply values - per by pad
 374   //
 375
 376   Bool_t result = kTRUE;
 377   for (Int_t  idata = 0; idata< fNchannels; idata++){
 378     if((GetValue(idata)*roc->GetValue(idata)) <= 6.5535)
 379       SetValue(idata,GetValue(idata)*roc->GetValue(idata));
 380     else {
 381       SetValue(idata,0.0);
 382       result = kFALSE;
 383     }
 384   }
 385   return result;
 386 }
 387
 388 //______________________________________________________________________________________________
 389 Bool_t AliTRDCalROC::Divide(const AliTRDCalROC*  roc)
 390 {
 391   //
 392   // divide values
 393   //
 394
 395   Bool_t result = kTRUE;
 396   Float_t kEpsilon=0.00000000000000001;
 397   for (Int_t  idata = 0; idata< fNchannels; idata++){
 398     if (TMath::Abs(roc->GetValue(idata))>kEpsilon){
 399       if((GetValue(idata)/roc->GetValue(idata)) <= 6.5535)
 400         SetValue(idata,GetValue(idata)/roc->GetValue(idata));
 401       else {
 402         SetValue(idata,0.0);
 403         result = kFALSE;
 404       }
 405     }
 406     else result = kFALSE;
 407   }
 408   return result;
 409 }
 410 //______________________________________________________________________________________________
 411 Bool_t AliTRDCalROC::Unfold()
 412 {
 413   //
 414   // Compute the mean value per pad col
 415   // Divide with this value each pad col
 416   // This is for the noise study
 417   // Return kFALSE if one or more of the pad col was not normalised
 418   //
 419
 420   Bool_t result = kTRUE;
 421   Float_t kEpsilon=0.00000000000000001;
 422
 423   // calcul the mean value per col
 424   for(Int_t icol = 0; icol < fNcols; icol++){
 425
 426     Float_t mean = 0.0;
 427     Float_t nb   = 0.0;
 428
 429     for(Int_t irow = 0; irow < fNrows; irow++){
 430       if((GetValue(icol,irow) > 0.06) && (GetValue(icol,irow) < 0.15)){
 431         mean += GetValue(icol,irow);
 432         nb += 1.0;
 433       }
 434     }
 435
 436     if(nb > kEpsilon) {
 437
 438       mean = mean/nb;
 439
 440       if(mean > kEpsilon){
 441         for(Int_t irow = 0; irow < fNrows; irow++){
 442           Float_t value = GetValue(icol,irow);
 443           SetValue(icol,irow,(Float_t)(value/mean));
 444         }
 445       }
 446       else result = kFALSE;
 447
 448     }
 449
 450     else result = kFALSE;
 451
 452
 453   }
 454
 455
 456   return result;
 457 }
 458
 459 //__________________________________________________________________________________
 460 TH2F * AliTRDCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type,  Float_t mu)
 461 {
 462   //
 463   // make 2D histo
 464   // type -1 = user defined range
 465   //       0 = nsigma cut nsigma=min
 466   //       1 = delta cut around median delta=min
 467
 468   Float_t kEpsilonr = 0.005;
 469   gStyle->SetPalette(1);
 470
 471   if (type>=0){
 472     if (type==0){
 473       // nsigma range
 474       Float_t mean  = GetMean();
 475       Float_t sigma = GetRMS();
 476       Float_t nsigma = TMath::Abs(min);
 477       sigma *= nsigma;
 478       if(sigma < kEpsilonr) sigma = kEpsilonr;
 479       min = mean-sigma;
 480       max = mean+sigma;
 481     }
 482     if (type==1){
 483       // fixed range
 484       Float_t mean   = GetMedian();
 485       Float_t  delta = min;
 486       min = mean-delta;
 487       max = mean+delta;
 488     }
 489     if (type==2){
 490       Double_t sigma;
 491       Float_t mean  = GetLTM(&sigma,max);
 492       sigma*=min;
 493       if(sigma < kEpsilonr) sigma = kEpsilonr;
 494       min = mean-sigma;
 495       max = mean+sigma;
 496
 497     }
 498   }
 499
 500   char  name[1000];
 501   sprintf(name,"%s 2D Plane %d Chamber %d",GetTitle(),fPla, fCha);
 502   TH2F * his = new TH2F(name,name,fNrows,0, fNrows, fNcols, 0, fNcols);
 503   for (Int_t irow=0; irow<fNrows; irow++){
 504     for (Int_t icol=0; icol<fNcols; icol++){
 505       his->Fill(irow+0.5,icol+0.5,GetValue(icol,irow)*mu);
 506     }
 507   }
 508   his->SetStats(0);
 509   his->SetMaximum(max);
 510   his->SetMinimum(min);
 511   return his;
 512 }
 513
 514 //_______________________________________________________________________________________
 515 TH1F * AliTRDCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type,  Float_t mu)
 516 {
 517   //
 518   // make 1D histo
 519   // type -1 = user defined range
 520   //       0 = nsigma cut nsigma=min
 521   //       1 = delta cut around median delta=min
 522
 523   Float_t kEpsilonr = 0.005;
 524
 525   if (type>=0){
 526     if (type==0){
 527       // nsigma range
 528       Float_t mean  = GetMean();
 529       Float_t sigma = GetRMS();
 530       Float_t nsigma = TMath::Abs(min);
 531       sigma *= nsigma;
 532       if(sigma < kEpsilonr) sigma = kEpsilonr;
 533       min = mean-sigma;
 534       max = mean+sigma;
 535     }
 536     if (type==1){
 537       // fixed range
 538       Float_t mean   = GetMedian();
 539       Float_t  delta = min;
 540       min = mean-delta;
 541       max = mean+delta;
 542     }
 543     if (type==2){
 544       //
 545       // LTM mean +- nsigma
 546       //
 547       Double_t sigma;
 548       Float_t mean  = GetLTM(&sigma,max);
 549       sigma*=min;
 550       if(sigma < kEpsilonr) sigma = kEpsilonr;
 551       min = mean-sigma;
 552       max = mean+sigma;
 553     }
 554   }
 555   char  name[1000];
 556   sprintf(name,"%s 1D Plane %d Chamber %d",GetTitle(),fPla, fCha);
 557   TH1F * his = new TH1F(name,name,100, min,max);
 558   for (Int_t irow=0; irow<fNrows; irow++){
 559     for (Int_t icol=0; icol<fNcols; icol++){
 560       his->Fill(GetValue(icol,irow)*mu);
 561     }
 562   }
 563   return his;
 564 }