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          ddata[NPoints]= (Double_t) fData[NPoints]/10000;
 245          NPoints++;
 246       }
 247    }
 248    Double_t mean = TMath::Mean(NPoints,ddata);
 249    delete [] ddata;
 250    return mean;
 251 }
 252
 253 //_______________________________________________________________________________________
 254 Double_t AliTRDCalROC::GetMedian(AliTRDCalROC* outlierROC)
 255 {
 256   //
 257   // Calculate the median
 258   //
 259
 260   Double_t *ddata = new Double_t[fNchannels];
 261    Int_t NPoints = 0;
 262    for (Int_t i=0;i<fNchannels;i++) {
 263        if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 264          ddata[NPoints]= (Double_t) fData[NPoints]/10000;
 265          NPoints++;
 266       }
 267    }
 268    Double_t mean = TMath::Median(NPoints,ddata);
 269    delete [] ddata;
 270    return mean;
 271 }
 272
 273 //____________________________________________________________________________________________
 274 Double_t AliTRDCalROC::GetRMS(AliTRDCalROC* outlierROC)
 275 {
 276   //
 277   // Calculate the RMS
 278   //
 279
 280   Double_t *ddata = new Double_t[fNchannels];
 281    Int_t NPoints = 0;
 282    for (Int_t i=0;i<fNchannels;i++) {
 283      if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
 284          ddata[NPoints]= (Double_t)fData[NPoints]/10000;
 285          NPoints++;
 286       }
 287    }
 288    Double_t mean = TMath::RMS(NPoints,ddata);
 289    delete [] ddata;
 290    return mean;
 291 }
 292
 293 //______________________________________________________________________________________________
 294 Double_t AliTRDCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTRDCalROC* outlierROC)
 295 {
 296   //
 297   //  Calculate LTM mean and sigma
 298   //
 299
 300   Double_t *ddata = new Double_t[fNchannels];
 301   Double_t mean=0, lsigma=0;
 302   UInt_t NPoints = 0;
 303   for (Int_t i=0;i<fNchannels;i++) {
 304      if (!outlierROC || !(outlierROC->GetValue(i))) {
 305         ddata[NPoints]= (Double_t) fData[NPoints]/10000;
 306         NPoints++;
 307      }
 308   }
 309   Int_t hh = TMath::Min(TMath::Nint(fraction *NPoints), Int_t(NPoints));
 310   AliMathBase::EvaluateUni(NPoints,ddata, mean, lsigma, hh);
 311   if (sigma) *sigma=lsigma;
 312   delete [] ddata;
 313   return mean;
 314 }
 315
 316 //___________________________________________________________________________________
 317 Bool_t AliTRDCalROC::Add(Float_t c1)
 318 {
 319   //
 320   // add constant
 321   //
 322   Bool_t result = kTRUE;
 323   for (Int_t  idata = 0; idata< fNchannels; idata++) {
 324     if(((GetValue(idata)+c1) <= 6.5535) && ((GetValue(idata)+c1) >= 0.0)) SetValue(idata,GetValue(idata)+c1);
 325     else {
 326       SetValue(idata,0.0);
 327       result = kFALSE;
 328     }
 329   }
 330   return result;
 331 }
 332
 333 //_______________________________________________________________________________________
 334 Bool_t AliTRDCalROC::Multiply(Float_t c1)
 335 {
 336   //
 337   // multiply constant
 338   //
 339   Bool_t result = kTRUE;
 340   if(c1 < 0) return kFALSE;
 341   for (Int_t  idata = 0; idata< fNchannels; idata++) {
 342     if((GetValue(idata)*c1) <= 6.5535) SetValue(idata,GetValue(idata)*c1);
 343     else {
 344       SetValue(idata,0.0);
 345       result = kFALSE;
 346     }
 347   }
 348   return result;
 349 }
 350
 351 //____________________________________________________________________________________________
 352 Bool_t AliTRDCalROC::Add(const AliTRDCalROC * roc, Double_t c1)
 353 {
 354   //
 355   // add values
 356   //
 357   Bool_t result = kTRUE;
 358   for (Int_t  idata = 0; idata< fNchannels; idata++){
 359     if(((GetValue(idata)+roc->GetValue(idata)*c1) <= 6.5535) &&
 360        ((GetValue(idata)+roc->GetValue(idata)*c1) >= 0.0))
 361           SetValue(idata,GetValue(idata)+roc->GetValue(idata)*c1);
 362     else {
 363       SetValue(idata,0.0);
 364       result = kFALSE;
 365     }
 366   }
 367   return result;
 368 }
 369
 370 //____________________________________________________________________________________________
 371 Bool_t AliTRDCalROC::Multiply(const AliTRDCalROC*  roc)
 372 {
 373   //
 374   // multiply values - per by pad
 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   Bool_t result = kTRUE;
 395   Float_t kEpsilon=0.00000000000000001;
 396   for (Int_t  idata = 0; idata< fNchannels; idata++){
 397     if (TMath::Abs(roc->GetValue(idata))>kEpsilon){
 398       if((GetValue(idata)/roc->GetValue(idata)) <= 6.5535)
 399         SetValue(idata,GetValue(idata)/roc->GetValue(idata));
 400       else {
 401         SetValue(idata,0.0);
 402         result = kFALSE;
 403       }
 404     }
 405     else result = kFALSE;
 406   }
 407   return result;
 408 }
 409
 410 //__________________________________________________________________________________
 411 TH2F * AliTRDCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type,  Float_t mu)
 412 {
 413   //
 414   // make 2D histo
 415   // type -1 = user defined range
 416   //       0 = nsigma cut nsigma=min
 417   //       1 = delta cut around median delta=min
 418   Float_t kEpsilonr = 0.005;
 419   gStyle->SetPalette(1);
 420
 421   if (type>=0){
 422     if (type==0){
 423       // nsigma range
 424       Float_t mean  = GetMean();
 425       Float_t sigma = GetRMS();
 426       Float_t nsigma = TMath::Abs(min);
 427       sigma *= nsigma;
 428       if(sigma < kEpsilonr) sigma = kEpsilonr;
 429       min = mean-sigma;
 430       max = mean+sigma;
 431     }
 432     if (type==1){
 433       // fixed range
 434       Float_t mean   = GetMedian();
 435       Float_t  delta = min;
 436       min = mean-delta;
 437       max = mean+delta;
 438     }
 439     if (type==2){
 440       Double_t sigma;
 441       Float_t mean  = GetLTM(&sigma,max);
 442       sigma*=min;
 443       if(sigma < kEpsilonr) sigma = kEpsilonr;
 444       min = mean-sigma;
 445       max = mean+sigma;
 446
 447     }
 448   }
 449
 450   char  name[1000];
 451   sprintf(name,"%s 2D Plane %d Chamber %d",GetTitle(),fPla, fCha);
 452   TH2F * his = new TH2F(name,name,fNrows,0, fNrows, fNcols, 0, fNcols);
 453   for (Int_t irow=0; irow<fNrows; irow++){
 454     for (Int_t icol=0; icol<fNcols; icol++){
 455       his->Fill(irow+0.5,icol+0.5,GetValue(icol,irow)*mu);
 456     }
 457   }
 458   his->SetStats(0);
 459   his->SetMaximum(max);
 460   his->SetMinimum(min);
 461   return his;
 462 }
 463
 464 //_______________________________________________________________________________________
 465 TH1F * AliTRDCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type,  Float_t mu)
 466 {
 467   //
 468   // make 1D histo
 469   // type -1 = user defined range
 470   //       0 = nsigma cut nsigma=min
 471   //       1 = delta cut around median delta=min
 472   Float_t kEpsilonr = 0.005;
 473
 474
 475   if (type>=0){
 476     if (type==0){
 477       // nsigma range
 478       Float_t mean  = GetMean();
 479       Float_t sigma = GetRMS();
 480       Float_t nsigma = TMath::Abs(min);
 481       sigma *= nsigma;
 482       if(sigma < kEpsilonr) sigma = kEpsilonr;
 483       min = mean-sigma;
 484       max = mean+sigma;
 485     }
 486     if (type==1){
 487       // fixed range
 488       Float_t mean   = GetMedian();
 489       Float_t  delta = min;
 490       min = mean-delta;
 491       max = mean+delta;
 492     }
 493     if (type==2){
 494       //
 495       // LTM mean +- nsigma
 496       //
 497       Double_t sigma;
 498       Float_t mean  = GetLTM(&sigma,max);
 499       sigma*=min;
 500       if(sigma < kEpsilonr) sigma = kEpsilonr;
 501       min = mean-sigma;
 502       max = mean+sigma;
 503     }
 504   }
 505   char  name[1000];
 506   sprintf(name,"%s 1D Plane %d Chamber %d",GetTitle(),fPla, fCha);
 507   TH1F * his = new TH1F(name,name,100, min,max);
 508   for (Int_t irow=0; irow<fNrows; irow++){
 509     for (Int_t icol=0; icol<fNcols; icol++){
 510       his->Fill(GetValue(icol,irow)*mu);
 511     }
 512   }
 513   return his;
 514 }