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