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