TRD/Cal/AliTRDCalDet.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 //  TRD calibration class for parameters which saved per detector            //
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TMath.h>
  25 #include <TH1F.h>
  26 #include <TH2F.h>
  27 #include <TStyle.h>
  28
  29 #include "AliTRDCalDet.h"
  30 #include "AliTRDgeometry.h"
  31 #include "AliMathBase.h"
  32 #include "AliTRDpadPlane.h"
  33
  34 ClassImp(AliTRDCalDet)
  35
  36 //_____________________________________________________________________________
  37 AliTRDCalDet::AliTRDCalDet():TNamed()
  38 {
  39   //
  40   // AliTRDCalDet default constructor
  41   //
  42
  43   for (Int_t idet = 0; idet < kNdet; idet++) {
  44     fData[idet] = 0;
  45   }
  46
  47 }
  48
  49 //_____________________________________________________________________________
  50 AliTRDCalDet::AliTRDCalDet(const Text_t *name, const Text_t *title)
  51                 :TNamed(name,title)
  52 {
  53   //
  54   // AliTRDCalDet constructor
  55   //
  56
  57   for (Int_t idet = 0; idet < kNdet; idet++) {
  58     fData[idet] = 0;
  59   }
  60
  61 }
  62
  63 //_____________________________________________________________________________
  64 AliTRDCalDet::AliTRDCalDet(const AliTRDCalDet &c):TNamed(c)
  65 {
  66   //
  67   // AliTRDCalDet copy constructor
  68   //
  69
  70   ((AliTRDCalDet &) c).Copy(*this);
  71
  72 }
  73
  74 ///_____________________________________________________________________________
  75 AliTRDCalDet::~AliTRDCalDet()
  76 {
  77   //
  78   // AliTRDCalDet destructor
  79   //
  80
  81 }
  82
  83 //_____________________________________________________________________________
  84 AliTRDCalDet &AliTRDCalDet::operator=(const AliTRDCalDet &c)
  85 {
  86   //
  87   // Assignment operator
  88   //
  89
  90   if (this != &c) ((AliTRDCalDet &) c).Copy(*this);
  91   return *this;
  92
  93 }
  94
  95 //_____________________________________________________________________________
  96 void AliTRDCalDet::Copy(TObject &c) const
  97 {
  98   //
  99   // Copy function
 100   //
 101
 102   for (Int_t idet = 0; idet < kNdet; idet++) {
 103     ((AliTRDCalDet &) c).fData[idet] = fData[idet];
 104   }
 105
 106   TObject::Copy(c);
 107
 108 }
 109
 110 //___________________________________________________________________________________
 111 Double_t AliTRDCalDet::GetMean(AliTRDCalDet * const outlierDet) const
 112 {
 113   //
 114   // Calculate the mean
 115   //
 116
 117   if (!outlierDet) return TMath::Mean(kNdet,fData);
 118   Double_t *ddata = new Double_t[kNdet];
 119   Int_t nPoints = 0;
 120   for (Int_t i=0;i<kNdet;i++) {
 121     if (!(outlierDet->GetValue(i))) {
 122       ddata[nPoints]= fData[nPoints];
 123       nPoints++;
 124     }
 125   }
 126   Double_t mean = TMath::Mean(nPoints,ddata);
 127   delete [] ddata;
 128   return mean;
 129 }
 130
 131 //_______________________________________________________________________________________
 132 Double_t AliTRDCalDet::GetMedian(AliTRDCalDet * const outlierDet) const
 133 {
 134   //
 135   // Calculate the median
 136   //
 137
 138   if (!outlierDet) return (Double_t) TMath::Median(kNdet,fData);
 139   Double_t *ddata = new Double_t[kNdet];
 140   Int_t nPoints = 0;
 141   for (Int_t i=0;i<kNdet;i++) {
 142     if (!(outlierDet->GetValue(i))) {
 143       ddata[nPoints]= fData[nPoints];
 144       nPoints++;
 145     }
 146   }
 147   Double_t mean = TMath::Median(nPoints,ddata);
 148   delete [] ddata;
 149   return mean;
 150
 151 }
 152
 153 //____________________________________________________________________________________________
 154 Double_t AliTRDCalDet::GetRMS(AliTRDCalDet * const outlierDet) const
 155 {
 156   //
 157   // Calculate the RMS
 158   //
 159
 160   if (!outlierDet) return TMath::RMS(kNdet,fData);
 161   Double_t *ddata = new Double_t[kNdet];
 162   Int_t nPoints = 0;
 163   for (Int_t i=0;i<kNdet;i++) {
 164     if (!(outlierDet->GetValue(i))) {
 165       ddata[nPoints]= fData[nPoints];
 166       nPoints++;
 167     }
 168   }
 169   Double_t mean = TMath::RMS(nPoints,ddata);
 170   delete [] ddata;
 171   return mean;
 172 }
 173
 174 //____________________________________________________________________________________________
 175 Double_t AliTRDCalDet::GetRMSRobust(Double_t robust) const
 176 {
 177   //
 178   // Calculate the robust RMS
 179   //
 180
 181   // sorted
 182   Int_t *index = new Int_t[kNdet];
 183   TMath::Sort((Int_t)kNdet,fData,index);
 184
 185   // reject
 186   Double_t reject = (Int_t) (kNdet*(1-robust)/2.0);
 187   if(reject <= 0.0) reject = 0.0;
 188   if(reject >= kNdet) reject = 0.0;
 189   //printf("Rejecter %f\n",reject);
 190
 191   Double_t *ddata = new Double_t[kNdet];
 192   Int_t nPoints = 0;
 193   for (Int_t i=0;i<kNdet;i++) {
 194     Bool_t rej = kFALSE;
 195     for(Int_t k = 0; k < reject; k++){
 196       if(i==index[k]) rej = kTRUE;
 197       if(i==index[kNdet-(k+1)]) rej  = kTRUE;
 198     }
 199     if(!rej){
 200       ddata[nPoints]= fData[i];
 201       nPoints++;
 202     }
 203   }
 204   //printf("Number of points %d\n",nPoints);
 205   Double_t mean = TMath::RMS(nPoints,ddata);
 206   delete [] ddata;
 207   delete [] index;
 208   return mean;
 209 }
 210
 211 //______________________________________________________________________________________________
 212 Double_t AliTRDCalDet::GetLTM(Double_t *sigma
 213                             , Double_t fraction
 214                             , AliTRDCalDet * const outlierDet)
 215 {
 216   //
 217   //  Calculate LTM mean and sigma
 218   //
 219
 220   Double_t *ddata = new Double_t[kNdet];
 221   Double_t mean=0, lsigma=0;
 222   UInt_t nPoints = 0;
 223   for (Int_t i=0;i<kNdet;i++) {
 224      if (!outlierDet || !(outlierDet->GetValue(i))) {
 225         ddata[nPoints]= fData[nPoints];
 226         nPoints++;
 227      }
 228   }
 229   Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
 230   AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
 231   if (sigma) *sigma=lsigma;
 232   delete [] ddata;
 233   return mean;
 234 }
 235
 236 //_________________________________________________________________________
 237 TH1F * AliTRDCalDet::MakeHisto1Distribution(Float_t min, Float_t max,Int_t type)
 238 {
 239   //
 240   // make 1D histo
 241   // type -1 = user defined range
 242   //       0 = nsigma cut nsigma=min
 243   //       1 = delta cut around median delta=min
 244   //
 245
 246   if (type>=0){
 247     if (type==0){
 248       // nsigma range
 249       Float_t mean  = GetMean();
 250       Float_t sigma = GetRMS();
 251       Float_t nsigma = TMath::Abs(min);
 252       min = mean-nsigma*sigma;
 253       max = mean+nsigma*sigma;
 254     }
 255     if (type==1){
 256       // fixed range
 257       Float_t mean   = GetMedian();
 258       Float_t  delta = min;
 259       min = mean-delta;
 260       max = mean+delta;
 261     }
 262     if (type==2){
 263       //
 264       // LTM mean +- nsigma
 265       //
 266       Double_t sigma;
 267       Float_t mean  = GetLTM(&sigma,max);
 268       sigma*=min;
 269       min = mean-sigma;
 270       max = mean+sigma;
 271     }
 272   }
 273   char  name[1000];
 274   snprintf(name,1000,"%s CalDet 1Distribution",GetTitle());
 275   TH1F * his = new TH1F(name,name,100, min,max);
 276   for (Int_t idet=0; idet<kNdet; idet++){
 277     his->Fill(GetValue(idet));
 278   }
 279   return his;
 280 }
 281
 282 //________________________________________________________________________________
 283 TH1F * AliTRDCalDet::MakeHisto1DAsFunctionOfDet(Float_t min, Float_t max,Int_t type)
 284 {
 285   //
 286   // make 1D histo
 287   // type -1 = user defined range
 288   //       0 = nsigma cut nsigma=min
 289   //       1 = delta cut around median delta=min
 290   //
 291
 292   if (type>=0){
 293     if (type==0){
 294       // nsigma range
 295       Float_t mean  = GetMean();
 296       Float_t sigma = GetRMS();
 297       Float_t nsigma = TMath::Abs(min);
 298       min = mean-nsigma*sigma;
 299       max = mean+nsigma*sigma;
 300     }
 301     if (type==1){
 302       // fixed range
 303       Float_t mean   = GetMedian();
 304       Float_t  delta = min;
 305       min = mean-delta;
 306       max = mean+delta;
 307     }
 308     if (type==2){
 309       Double_t sigma;
 310       Float_t mean  = GetLTM(&sigma,max);
 311       sigma*=min;
 312       min = mean-sigma;
 313       max = mean+sigma;
 314
 315     }
 316   }
 317
 318   char  name[1000];
 319   snprintf(name,1000,"%s CalDet as function of det",GetTitle());
 320   TH1F * his = new TH1F(name,name,kNdet, 0, kNdet);
 321   for(Int_t det = 0; det< kNdet; det++){
 322     his->Fill(det+0.5,GetValue(det));
 323   }
 324
 325   his->SetMaximum(max);
 326   his->SetMinimum(min);
 327   return his;
 328 }
 329
 330 //_____________________________________________________________________________
 331 TH2F *AliTRDCalDet::MakeHisto2DCh(Int_t ch, Float_t min, Float_t max, Int_t type)
 332 {
 333   //
 334   // Make 2D graph
 335   // ch    - chamber
 336   // type -1 = user defined range
 337   //       0 = nsigma cut nsigma=min
 338   //       1 = delta cut around median delta=min
 339   //
 340
 341   gStyle->SetPalette(1);
 342   if (type>=0){
 343     if (type==0){
 344       // nsigma range
 345       Float_t mean  = GetMean();
 346       Float_t sigma = GetRMS();
 347       Float_t nsigma = TMath::Abs(min);
 348       min = mean-nsigma*sigma;
 349       max = mean+nsigma*sigma;
 350     }
 351     if (type==1){
 352       // fixed range
 353       Float_t mean   = GetMedian();
 354       Float_t  delta = min;
 355       min = mean-delta;
 356       max = mean+delta;
 357     }
 358     if (type==2){
 359       Double_t sigma;
 360       Float_t mean  = GetLTM(&sigma,max);
 361       sigma*=min;
 362       min = mean-sigma;
 363       max = mean+sigma;
 364
 365     }
 366   }
 367
 368   AliTRDgeometry *trdGeo = new AliTRDgeometry();
 369
 370   Double_t poslocal[3]  = {0.0,0.0,0.0};
 371   Double_t posglobal[3] = {0.0,0.0,0.0};
 372
 373   char  name[1000];
 374   snprintf(name,1000,"%s CalDet 2D ch %d",GetTitle(),ch);
 375   TH2F * his = new TH2F(name, name, 400,-400.0,400.0,400,-400.0,400.0);
 376
 377   // Where we begin
 378   Int_t offsetch = 6*ch;
 379
 380
 381   for (Int_t isec = 0; isec < kNsect; isec++){
 382     for(Int_t ipl = 0; ipl < kNplan; ipl++){
 383       Int_t det   = offsetch+isec*30+ipl;
 384       AliTRDpadPlane *padPlane = trdGeo->GetPadPlane(ipl,ch);
 385       for (Int_t icol=0; icol<padPlane->GetNcols(); icol++){
 386         poslocal[0] = trdGeo->GetTime0(ipl);
 387         poslocal[2] = padPlane->GetRowPos(0);
 388         poslocal[1] = padPlane->GetColPos(icol);
 389         trdGeo->RotateBack(det,poslocal,posglobal);
 390         Int_t binx = 1+TMath::Nint((posglobal[0]+400.0)*0.5);
 391         Int_t biny = 1+TMath::Nint((posglobal[1]+400.0)*0.5);
 392         his->SetBinContent(binx,biny,fData[det]);
 393       }
 394     }
 395   }
 396   his->SetXTitle("x (cm)");
 397   his->SetYTitle("y (cm)");
 398   his->SetStats(0);
 399   his->SetMaximum(max);
 400   his->SetMinimum(min);
 401   delete trdGeo;
 402   return his;
 403 }
 404
 405 //_____________________________________________________________________________
 406 TH2F *AliTRDCalDet::MakeHisto2DSmPl(Int_t sm, Int_t pl, Float_t min, Float_t max, Int_t type)
 407 {
 408   //
 409   // Make 2D graph
 410   // sm    - supermodule number
 411   // pl    - plane number
 412   // type -1 = user defined range
 413   //       0 = nsigma cut nsigma=min
 414   //       1 = delta cut around median delta=min
 415   //
 416
 417   gStyle->SetPalette(1);
 418   if (type>=0){
 419     if (type==0){
 420       // nsigma range
 421       Float_t mean  = GetMean();
 422       Float_t sigma = GetRMS();
 423       Float_t nsigma = TMath::Abs(min);
 424       min = mean-nsigma*sigma;
 425       max = mean+nsigma*sigma;
 426     }
 427     if (type==1){
 428       // fixed range
 429       Float_t mean   = GetMedian();
 430       Float_t  delta = min;
 431       min = mean-delta;
 432       max = mean+delta;
 433     }
 434     if (type==2){
 435       Double_t sigma;
 436       Float_t mean  = GetLTM(&sigma,max);
 437       sigma*=min;
 438       min = mean-sigma;
 439       max = mean+sigma;
 440
 441     }
 442   }
 443
 444   AliTRDgeometry *trdGeo = new AliTRDgeometry();
 445   AliTRDpadPlane *padPlane0 = trdGeo->GetPadPlane(pl,0);
 446   Double_t row0    = padPlane0->GetRow0();
 447   Double_t col0    = padPlane0->GetCol0();
 448
 449   char  name[1000];
 450   snprintf(name,1000,"%s CalDet 2D sm %d and pl %d",GetTitle(),sm,pl);
 451   TH2F * his = new TH2F( name, name, 5,  -TMath::Abs(row0),  TMath::Abs(row0)
 452                                    , 4,-2*TMath::Abs(col0),2*TMath::Abs(col0));
 453
 454   // Where we begin
 455   Int_t offsetsmpl = 30*sm+pl;
 456
 457
 458   for (Int_t k = 0; k < kNcham; k++){
 459     Int_t det = offsetsmpl+k*6;
 460     Int_t kb  = kNcham-1-k;
 461     his->SetBinContent(kb+1,2,fData[det]);
 462     his->SetBinContent(kb+1,3,fData[det]);
 463   }
 464   his->SetXTitle("z (cm)");
 465   his->SetYTitle("y (cm)");
 466   his->SetStats(0);
 467   his->SetMaximum(max);
 468   his->SetMinimum(min);
 469   delete trdGeo;
 470   return his;
 471 }
 472
 473 //_____________________________________________________________________________
 474 void AliTRDCalDet::Add(Float_t c1)
 475 {
 476   //
 477   // Add constant for all detectors
 478   //
 479
 480   for (Int_t idet = 0; idet < kNdet; idet++) {
 481      fData[idet] += c1;
 482   }
 483 }
 484
 485 //_____________________________________________________________________________
 486 void AliTRDCalDet::Multiply(Float_t c1)
 487 {
 488     //
 489     // multiply constant for all detectors
 490     //
 491     for (Int_t idet = 0; idet < kNdet; idet++) {
 492       fData[idet] *= c1;
 493     }
 494 }
 495
 496 //_____________________________________________________________________________
 497 void AliTRDCalDet::Add(const AliTRDCalDet * calDet, Double_t c1)
 498 {
 499     //
 500     // add caldet channel by channel multiplied by c1
 501     //
 502     for (Int_t idet = 0; idet < kNdet; idet++) {
 503       fData[idet] += calDet->GetValue(idet)*c1;
 504     }
 505 }
 506
 507 //_____________________________________________________________________________
 508 void AliTRDCalDet::Multiply(const AliTRDCalDet * calDet)
 509 {
 510     //
 511     // multiply caldet channel by channel
 512     //
 513     for (Int_t idet = 0; idet < kNdet; idet++) {
 514       fData[idet] *= calDet->GetValue(idet);
 515     }
 516 }
 517
 518 //_____________________________________________________________________________
 519 void AliTRDCalDet::Divide(const AliTRDCalDet * calDet)
 520 {
 521     //
 522     // divide caldet channel by channel
 523     //
 524  Float_t kEpsilon=0.00000000000000001;
 525     for (Int_t idet = 0; idet < kNdet; idet++) {
 526       if (TMath::Abs(calDet->GetValue(idet))>kEpsilon){
 527         fData[idet] /= calDet->GetValue(idet);
 528         }
 529     }
 530 }
 531
 532