TRD/TRDbase/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) {
 190     return *this;
 191   }
 192
 193   fPla       = c.fPla;
 194   fCha       = c.fCha;
 195   fNrows     = c.fNrows;
 196   fNcols     = c.fNcols;
 197   fNchannels = c.fNchannels;
 198
 199   if (fData) {
 200     delete [] fData;
 201   }
 202   fData = new UShort_t[fNchannels];
 203   for (Int_t iBin = 0; iBin < fNchannels; iBin++) {
 204     fData[iBin] = ((AliTRDCalROC &) c).fData[iBin];
 205   }
 206
 207   return *this;
 208
 209 }
 210
 211 //_____________________________________________________________________________
 212 void AliTRDCalROC::Copy(TObject &c) const
 213 {
 214   //
 215   // Copy function
 216   //
 217
 218   ((AliTRDCalROC &) c).fPla          = fPla;
 219   ((AliTRDCalROC &) c).fCha          = fCha;
 220
 221   ((AliTRDCalROC &) c).fNrows        = fNrows;
 222   ((AliTRDCalROC &) c).fNcols        = fNcols;
 223
 224   Int_t iBin = 0;
 225
 226   ((AliTRDCalROC &) c).fNchannels = fNchannels;
 227
 228   if (((AliTRDCalROC &) c).fData) delete [] ((AliTRDCalROC &) c).fData;
 229   ((AliTRDCalROC &) c).fData = new UShort_t[fNchannels];
 230   for (iBin = 0; iBin < fNchannels; iBin++) {
 231     ((AliTRDCalROC &) c).fData[iBin] = fData[iBin];
 232   }
 233
 234   TObject::Copy(c);
 235
 236 }
 237
 238 //___________________________________________________________________________________
 239 Double_t AliTRDCalROC::GetMean(AliTRDCalROC* const outlierROC) const
 240 {
 241   //
 242   // Calculate the mean
 243   //
 244
 245    Double_t *ddata = new Double_t[fNchannels];
 246    Int_t nPoints = 0;
 247    for (Int_t i=0;i<fNchannels;i++) {
 248       if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 249         //if(fData[i] > 0.000000000000001){
 250          ddata[nPoints]= (Double_t) fData[i]/10000;
 251          nPoints++;
 252          //}
 253       }
 254    }
 255    Double_t mean = TMath::Mean(nPoints,ddata);
 256    delete [] ddata;
 257    return mean;
 258 }
 259
 260 //_______________________________________________________________________________________
 261 Double_t AliTRDCalROC::GetMedian(AliTRDCalROC* const outlierROC) const
 262 {
 263   //
 264   // Calculate the median
 265   //
 266
 267   Double_t *ddata = new Double_t[fNchannels];
 268    Int_t nPoints = 0;
 269    for (Int_t i=0;i<fNchannels;i++) {
 270        if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 271          if(fData[i] > 0.000000000000001){
 272            ddata[nPoints]= (Double_t) fData[i]/10000;
 273            nPoints++;
 274          }
 275        }
 276    }
 277    Double_t mean = TMath::Median(nPoints,ddata);
 278    delete [] ddata;
 279    return mean;
 280 }
 281
 282 //____________________________________________________________________________________________
 283 Double_t AliTRDCalROC::GetRMS(AliTRDCalROC* const outlierROC) const
 284 {
 285   //
 286   // Calculate the RMS
 287   //
 288
 289   Double_t *ddata = new Double_t[fNchannels];
 290   Int_t nPoints = 0;
 291   for (Int_t i=0;i<fNchannels;i++) {
 292     if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 293       //if(fData[i] > 0.000000000000001){
 294          ddata[nPoints]= (Double_t)fData[i]/10000;
 295          nPoints++;
 296          //}
 297     }
 298   }
 299   Double_t mean = TMath::RMS(nPoints,ddata);
 300   delete [] ddata;
 301   return mean;
 302 }
 303
 304 //______________________________________________________________________________________________
 305 Double_t AliTRDCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTRDCalROC* const outlierROC)
 306 {
 307   //
 308   //  Calculate LTM mean and sigma
 309   //
 310
 311   Double_t *ddata = new Double_t[fNchannels];
 312   Double_t mean=0, lsigma=0;
 313   UInt_t nPoints = 0;
 314   for (Int_t i=0;i<fNchannels;i++) {
 315      if (!outlierROC || !(outlierROC->GetValue(i))) {
 316        if(fData[i] > 0.000000000000001){
 317          ddata[nPoints]= (Double_t) fData[i]/10000;
 318          nPoints++;
 319        }
 320      }
 321   }
 322   Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
 323   AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
 324   if (sigma) *sigma=lsigma;
 325   delete [] ddata;
 326   return mean;
 327 }
 328
 329 //___________________________________________________________________________________
 330 Bool_t AliTRDCalROC::Add(Float_t c1)
 331 {
 332   //
 333   // add constant
 334   //
 335
 336   Bool_t result = kTRUE;
 337   for (Int_t  idata = 0; idata< fNchannels; idata++) {
 338     if(((GetValue(idata)+c1) <= 6.5535) && ((GetValue(idata)+c1) >= 0.0)) SetValue(idata,GetValue(idata)+c1);
 339     else {
 340       SetValue(idata,0.0);
 341       result = kFALSE;
 342     }
 343   }
 344   return result;
 345 }
 346
 347 //_______________________________________________________________________________________
 348 Bool_t AliTRDCalROC::Multiply(Float_t c1)
 349 {
 350   //
 351   // multiply constant
 352   //
 353
 354   Bool_t result = kTRUE;
 355   if(c1 < 0) return kFALSE;
 356   for (Int_t  idata = 0; idata< fNchannels; idata++) {
 357     if((GetValue(idata)*c1) <= 6.5535) SetValue(idata,GetValue(idata)*c1);
 358     else {
 359       SetValue(idata,0.0);
 360       result = kFALSE;
 361     }
 362   }
 363   return result;
 364 }
 365
 366 //____________________________________________________________________________________________
 367 Bool_t AliTRDCalROC::Add(const AliTRDCalROC * roc, Double_t c1)
 368 {
 369   //
 370   // add values
 371   //
 372
 373   Bool_t result = kTRUE;
 374   for (Int_t  idata = 0; idata< fNchannels; idata++){
 375     if(((GetValue(idata)+roc->GetValue(idata)*c1) <= 6.5535) &&
 376        ((GetValue(idata)+roc->GetValue(idata)*c1) >= 0.0))
 377           SetValue(idata,GetValue(idata)+roc->GetValue(idata)*c1);
 378     else {
 379       SetValue(idata,0.0);
 380       result = kFALSE;
 381     }
 382   }
 383   return result;
 384 }
 385
 386 //____________________________________________________________________________________________
 387 Bool_t AliTRDCalROC::Multiply(const AliTRDCalROC*  roc)
 388 {
 389   //
 390   // multiply values - per by pad
 391   //
 392
 393   Bool_t result = kTRUE;
 394   for (Int_t  idata = 0; idata< fNchannels; idata++){
 395     if((GetValue(idata)*roc->GetValue(idata)) <= 6.5535)
 396       SetValue(idata,GetValue(idata)*roc->GetValue(idata));
 397     else {
 398       SetValue(idata,0.0);
 399       result = kFALSE;
 400     }
 401   }
 402   return result;
 403 }
 404
 405 //______________________________________________________________________________________________
 406 Bool_t AliTRDCalROC::Divide(const AliTRDCalROC*  roc)
 407 {
 408   //
 409   // divide values
 410   //
 411
 412   Bool_t result = kTRUE;
 413   Float_t kEpsilon=0.00000000000000001;
 414   for (Int_t  idata = 0; idata< fNchannels; idata++){
 415     if (TMath::Abs(roc->GetValue(idata))>kEpsilon){
 416       if((GetValue(idata)/roc->GetValue(idata)) <= 6.5535)
 417         SetValue(idata,GetValue(idata)/roc->GetValue(idata));
 418       else {
 419         SetValue(idata,0.0);
 420         result = kFALSE;
 421       }
 422     }
 423     else result = kFALSE;
 424   }
 425   return result;
 426 }
 427 //______________________________________________________________________________________________
 428 Bool_t AliTRDCalROC::Unfold()
 429 {
 430   //
 431   // Compute the mean value per pad col
 432   // Divide with this value each pad col
 433   // This is for the noise study
 434   // Return kFALSE if one or more of the pad col was not normalised
 435   //
 436
 437   Bool_t result = kTRUE;
 438   Float_t kEpsilon=0.00000000000000001;
 439
 440   // calcul the mean value per col
 441   for(Int_t icol = 0; icol < fNcols; icol++){
 442
 443     Float_t mean = 0.0;
 444     Float_t nb   = 0.0;
 445
 446     for(Int_t irow = 0; irow < fNrows; irow++){
 447       if((GetValue(icol,irow) > 0.06) && (GetValue(icol,irow) < 0.15)){
 448         mean += GetValue(icol,irow);
 449         nb += 1.0;
 450       }
 451     }
 452
 453     if(nb > kEpsilon) {
 454
 455       mean = mean/nb;
 456
 457       if(mean > kEpsilon){
 458         for(Int_t irow = 0; irow < fNrows; irow++){
 459           Float_t value = GetValue(icol,irow);
 460           SetValue(icol,irow,(Float_t)(value/mean));
 461         }
 462       }
 463       else result = kFALSE;
 464
 465     }
 466
 467     else result = kFALSE;
 468
 469
 470   }
 471
 472
 473   return result;
 474 }
 475
 476 //__________________________________________________________________________________
 477 TH2F * AliTRDCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type,  Float_t mu)
 478 {
 479   //
 480   // make 2D histo
 481   // type -1 = user defined range
 482   //       0 = nsigma cut nsigma=min
 483   //       1 = delta cut around median delta=min
 484
 485   Float_t kEpsilonr = 0.005;
 486   gStyle->SetPalette(1);
 487
 488   if (type>=0){
 489     if (type==0){
 490       // nsigma range
 491       Float_t mean  = GetMean();
 492       Float_t sigma = GetRMS();
 493       Float_t nsigma = TMath::Abs(min);
 494       sigma *= nsigma;
 495       if(sigma < kEpsilonr) sigma = kEpsilonr;
 496       min = mean-sigma;
 497       max = mean+sigma;
 498     }
 499     if (type==1){
 500       // fixed range
 501       Float_t mean   = GetMedian();
 502       Float_t  delta = min;
 503       min = mean-delta;
 504       max = mean+delta;
 505     }
 506     if (type==2){
 507       Double_t sigma;
 508       Float_t mean  = GetLTM(&sigma,max);
 509       sigma*=min;
 510       if(sigma < kEpsilonr) sigma = kEpsilonr;
 511       min = mean-sigma;
 512       max = mean+sigma;
 513
 514     }
 515   }
 516
 517   char  name[1000];
 518   snprintf(name,1000,"%s 2D Plane %d Chamber %d",GetTitle(),fPla, fCha);
 519   TH2F * his = new TH2F(name,name,fNrows,0, fNrows, fNcols, 0, fNcols);
 520   for (Int_t irow=0; irow<fNrows; irow++){
 521     for (Int_t icol=0; icol<fNcols; icol++){
 522       his->Fill(irow+0.5,icol+0.5,GetValue(icol,irow)*mu);
 523     }
 524   }
 525   his->SetStats(0);
 526   his->SetMaximum(max);
 527   his->SetMinimum(min);
 528   return his;
 529 }
 530
 531 //_______________________________________________________________________________________
 532 TH1F * AliTRDCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type,  Float_t mu)
 533 {
 534   //
 535   // make 1D histo
 536   // type -1 = user defined range
 537   //       0 = nsigma cut nsigma=min
 538   //       1 = delta cut around median delta=min
 539
 540   Float_t kEpsilonr = 0.005;
 541
 542   if (type>=0){
 543     if (type==0){
 544       // nsigma range
 545       Float_t mean  = GetMean();
 546       Float_t sigma = GetRMS();
 547       Float_t nsigma = TMath::Abs(min);
 548       sigma *= nsigma;
 549       if(sigma < kEpsilonr) sigma = kEpsilonr;
 550       min = mean-sigma;
 551       max = mean+sigma;
 552     }
 553     if (type==1){
 554       // fixed range
 555       Float_t mean   = GetMedian();
 556       Float_t  delta = min;
 557       min = mean-delta;
 558       max = mean+delta;
 559     }
 560     if (type==2){
 561       //
 562       // LTM mean +- nsigma
 563       //
 564       Double_t sigma;
 565       Float_t mean  = GetLTM(&sigma,max);
 566       sigma*=min;
 567       if(sigma < kEpsilonr) sigma = kEpsilonr;
 568       min = mean-sigma;
 569       max = mean+sigma;
 570     }
 571   }
 572   char  name[1000];
 573   snprintf(name,1000,"%s 1D Plane %d Chamber %d",GetTitle(),fPla, fCha);
 574   TH1F * his = new TH1F(name,name,100, min,max);
 575   for (Int_t irow=0; irow<fNrows; irow++){
 576     for (Int_t icol=0; icol<fNcols; icol++){
 577       his->Fill(GetValue(icol,irow)*mu);
 578     }
 579   }
 580   return his;
 581 }