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::GetLTM(Double_t *sigma
 176                             , Double_t fraction
 177                             , AliTRDCalDet * const outlierDet)
 178 {
 179   //
 180   //  Calculate LTM mean and sigma
 181   //
 182
 183   Double_t *ddata = new Double_t[kNdet];
 184   Double_t mean=0, lsigma=0;
 185   UInt_t nPoints = 0;
 186   for (Int_t i=0;i<kNdet;i++) {
 187      if (!outlierDet || !(outlierDet->GetValue(i))) {
 188         ddata[nPoints]= fData[nPoints];
 189         nPoints++;
 190      }
 191   }
 192   Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
 193   AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
 194   if (sigma) *sigma=lsigma;
 195   delete [] ddata;
 196   return mean;
 197 }
 198
 199 //_________________________________________________________________________
 200 TH1F * AliTRDCalDet::MakeHisto1Distribution(Float_t min, Float_t max,Int_t type)
 201 {
 202   //
 203   // make 1D histo
 204   // type -1 = user defined range
 205   //       0 = nsigma cut nsigma=min
 206   //       1 = delta cut around median delta=min
 207   //
 208
 209   if (type>=0){
 210     if (type==0){
 211       // nsigma range
 212       Float_t mean  = GetMean();
 213       Float_t sigma = GetRMS();
 214       Float_t nsigma = TMath::Abs(min);
 215       min = mean-nsigma*sigma;
 216       max = mean+nsigma*sigma;
 217     }
 218     if (type==1){
 219       // fixed range
 220       Float_t mean   = GetMedian();
 221       Float_t  delta = min;
 222       min = mean-delta;
 223       max = mean+delta;
 224     }
 225     if (type==2){
 226       //
 227       // LTM mean +- nsigma
 228       //
 229       Double_t sigma;
 230       Float_t mean  = GetLTM(&sigma,max);
 231       sigma*=min;
 232       min = mean-sigma;
 233       max = mean+sigma;
 234     }
 235   }
 236   char  name[1000];
 237   sprintf(name,"%s CalDet 1Distribution",GetTitle());
 238   TH1F * his = new TH1F(name,name,100, min,max);
 239   for (Int_t idet=0; idet<kNdet; idet++){
 240     his->Fill(GetValue(idet));
 241   }
 242   return his;
 243 }
 244
 245 //________________________________________________________________________________
 246 TH1F * AliTRDCalDet::MakeHisto1DAsFunctionOfDet(Float_t min, Float_t max,Int_t type)
 247 {
 248   //
 249   // make 1D histo
 250   // type -1 = user defined range
 251   //       0 = nsigma cut nsigma=min
 252   //       1 = delta cut around median delta=min
 253   //
 254
 255   if (type>=0){
 256     if (type==0){
 257       // nsigma range
 258       Float_t mean  = GetMean();
 259       Float_t sigma = GetRMS();
 260       Float_t nsigma = TMath::Abs(min);
 261       min = mean-nsigma*sigma;
 262       max = mean+nsigma*sigma;
 263     }
 264     if (type==1){
 265       // fixed range
 266       Float_t mean   = GetMedian();
 267       Float_t  delta = min;
 268       min = mean-delta;
 269       max = mean+delta;
 270     }
 271     if (type==2){
 272       Double_t sigma;
 273       Float_t mean  = GetLTM(&sigma,max);
 274       sigma*=min;
 275       min = mean-sigma;
 276       max = mean+sigma;
 277
 278     }
 279   }
 280
 281   char  name[1000];
 282   sprintf(name,"%s CalDet as function of det",GetTitle());
 283   TH1F * his = new TH1F(name,name,kNdet, 0, kNdet);
 284   for(Int_t det = 0; det< kNdet; det++){
 285     his->Fill(det+0.5,GetValue(det));
 286   }
 287
 288   his->SetMaximum(max);
 289   his->SetMinimum(min);
 290   return his;
 291 }
 292
 293 //_____________________________________________________________________________
 294 TH2F *AliTRDCalDet::MakeHisto2DCh(Int_t ch, Float_t min, Float_t max, Int_t type)
 295 {
 296   //
 297   // Make 2D graph
 298   // ch    - chamber
 299   // type -1 = user defined range
 300   //       0 = nsigma cut nsigma=min
 301   //       1 = delta cut around median delta=min
 302   //
 303
 304   gStyle->SetPalette(1);
 305   if (type>=0){
 306     if (type==0){
 307       // nsigma range
 308       Float_t mean  = GetMean();
 309       Float_t sigma = GetRMS();
 310       Float_t nsigma = TMath::Abs(min);
 311       min = mean-nsigma*sigma;
 312       max = mean+nsigma*sigma;
 313     }
 314     if (type==1){
 315       // fixed range
 316       Float_t mean   = GetMedian();
 317       Float_t  delta = min;
 318       min = mean-delta;
 319       max = mean+delta;
 320     }
 321     if (type==2){
 322       Double_t sigma;
 323       Float_t mean  = GetLTM(&sigma,max);
 324       sigma*=min;
 325       min = mean-sigma;
 326       max = mean+sigma;
 327
 328     }
 329   }
 330
 331   AliTRDgeometry *trdGeo = new AliTRDgeometry();
 332
 333   Double_t poslocal[3]  = {0.0,0.0,0.0};
 334   Double_t posglobal[3] = {0.0,0.0,0.0};
 335
 336   char  name[1000];
 337   sprintf(name,"%s CalDet 2D ch %d",GetTitle(),ch);
 338   TH2F * his = new TH2F(name, name, 400,-400.0,400.0,400,-400.0,400.0);
 339
 340   // Where we begin
 341   Int_t offsetch = 6*ch;
 342
 343
 344   for (Int_t isec = 0; isec < kNsect; isec++){
 345     for(Int_t ipl = 0; ipl < kNplan; ipl++){
 346       Int_t det   = offsetch+isec*30+ipl;
 347       AliTRDpadPlane *padPlane = trdGeo->GetPadPlane(ipl,ch);
 348       for (Int_t icol=0; icol<padPlane->GetNcols(); icol++){
 349         poslocal[0] = trdGeo->GetTime0(ipl);
 350         poslocal[2] = padPlane->GetRowPos(0);
 351         poslocal[1] = padPlane->GetColPos(icol);
 352         trdGeo->RotateBack(det,poslocal,posglobal);
 353         Int_t binx = 1+TMath::Nint((posglobal[0]+400.0)*0.5);
 354         Int_t biny = 1+TMath::Nint((posglobal[1]+400.0)*0.5);
 355         his->SetBinContent(binx,biny,fData[det]);
 356       }
 357     }
 358   }
 359   his->SetXTitle("x (cm)");
 360   his->SetYTitle("y (cm)");
 361   his->SetStats(0);
 362   his->SetMaximum(max);
 363   his->SetMinimum(min);
 364   delete trdGeo;
 365   return his;
 366 }
 367
 368 //_____________________________________________________________________________
 369 TH2F *AliTRDCalDet::MakeHisto2DSmPl(Int_t sm, Int_t pl, Float_t min, Float_t max, Int_t type)
 370 {
 371   //
 372   // Make 2D graph
 373   // sm    - supermodule number
 374   // pl    - plane number
 375   // type -1 = user defined range
 376   //       0 = nsigma cut nsigma=min
 377   //       1 = delta cut around median delta=min
 378   //
 379
 380   gStyle->SetPalette(1);
 381   if (type>=0){
 382     if (type==0){
 383       // nsigma range
 384       Float_t mean  = GetMean();
 385       Float_t sigma = GetRMS();
 386       Float_t nsigma = TMath::Abs(min);
 387       min = mean-nsigma*sigma;
 388       max = mean+nsigma*sigma;
 389     }
 390     if (type==1){
 391       // fixed range
 392       Float_t mean   = GetMedian();
 393       Float_t  delta = min;
 394       min = mean-delta;
 395       max = mean+delta;
 396     }
 397     if (type==2){
 398       Double_t sigma;
 399       Float_t mean  = GetLTM(&sigma,max);
 400       sigma*=min;
 401       min = mean-sigma;
 402       max = mean+sigma;
 403
 404     }
 405   }
 406
 407   AliTRDgeometry *trdGeo = new AliTRDgeometry();
 408   AliTRDpadPlane *padPlane0 = trdGeo->GetPadPlane(pl,0);
 409   Double_t row0    = padPlane0->GetRow0();
 410   Double_t col0    = padPlane0->GetCol0();
 411
 412   char  name[1000];
 413   sprintf(name,"%s CalDet 2D sm %d and pl %d",GetTitle(),sm,pl);
 414   TH2F * his = new TH2F( name, name, 5,  -TMath::Abs(row0),  TMath::Abs(row0)
 415                                    , 4,-2*TMath::Abs(col0),2*TMath::Abs(col0));
 416
 417   // Where we begin
 418   Int_t offsetsmpl = 30*sm+pl;
 419
 420
 421   for (Int_t k = 0; k < kNcham; k++){
 422     Int_t det = offsetsmpl+k*6;
 423     Int_t kb  = kNcham-1-k;
 424     his->SetBinContent(kb+1,2,fData[det]);
 425     his->SetBinContent(kb+1,3,fData[det]);
 426   }
 427   his->SetXTitle("z (cm)");
 428   his->SetYTitle("y (cm)");
 429   his->SetStats(0);
 430   his->SetMaximum(max);
 431   his->SetMinimum(min);
 432   delete trdGeo;
 433   return his;
 434 }
 435
 436 //_____________________________________________________________________________
 437 void AliTRDCalDet::Add(Float_t c1)
 438 {
 439   //
 440   // Add constant for all detectors
 441   //
 442
 443   for (Int_t idet = 0; idet < kNdet; idet++) {
 444      fData[idet] += c1;
 445   }
 446 }
 447
 448 //_____________________________________________________________________________
 449 void AliTRDCalDet::Multiply(Float_t c1)
 450 {
 451     //
 452     // multiply constant for all detectors
 453     //
 454     for (Int_t idet = 0; idet < kNdet; idet++) {
 455       fData[idet] *= c1;
 456     }
 457 }
 458
 459 //_____________________________________________________________________________
 460 void AliTRDCalDet::Add(const AliTRDCalDet * calDet, Double_t c1)
 461 {
 462     //
 463     // add caldet channel by channel multiplied by c1
 464     //
 465     for (Int_t idet = 0; idet < kNdet; idet++) {
 466       fData[idet] += calDet->GetValue(idet)*c1;
 467     }
 468 }
 469
 470 //_____________________________________________________________________________
 471 void AliTRDCalDet::Multiply(const AliTRDCalDet * calDet)
 472 {
 473     //
 474     // multiply caldet channel by channel
 475     //
 476     for (Int_t idet = 0; idet < kNdet; idet++) {
 477       fData[idet] *= calDet->GetValue(idet);
 478     }
 479 }
 480
 481 //_____________________________________________________________________________
 482 void AliTRDCalDet::Divide(const AliTRDCalDet * calDet)
 483 {
 484     //
 485     // divide caldet channel by channel
 486     //
 487  Float_t kEpsilon=0.00000000000000001;
 488     for (Int_t idet = 0; idet < kNdet; idet++) {
 489       if (TMath::Abs(calDet->GetValue(idet))>kEpsilon){
 490         fData[idet] /= calDet->GetValue(idet);
 491         }
 492     }
 493 }