TRD/Cal/AliTRDCalPIDLQ.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
  17 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 // Container for the distributions of dE/dx and the time bin of the          //
  20 // max. cluster for electrons and pions                                      //
  21 //                                                                           //
  22 // Author:                                                                   //
  23 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>                    //
  24 //                                                                           //
  25 ///////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TH1F.h>
  28 #include <TFile.h>
  29
  30 #include "AliLog.h"
  31 #include "AliPID.h"
  32
  33 #include "AliTRDCalPIDLQ.h"
  34
  35 ClassImp(AliTRDCalPIDLQ)
  36
  37 Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
  38
  39 //_________________________________________________________________________
  40 AliTRDCalPIDLQ::AliTRDCalPIDLQ()
  41   :TNamed()
  42   ,fNMom(0)
  43   ,fTrackMomentum(0)
  44   ,fMeanChargeRatio(0)
  45   ,fNbins(0)
  46   ,fBinSize(0)
  47   ,fHistdEdx(0)
  48   ,fHistTimeBin(0)
  49 {
  50   //
  51   //  The Default constructor
  52   //
  53
  54 }
  55
  56 //_________________________________________________________________________
  57 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
  58   :TNamed(name,title)
  59   ,fNMom(0)
  60   ,fTrackMomentum(0)
  61   ,fMeanChargeRatio(0)
  62   ,fNbins(0)
  63   ,fBinSize(0)
  64   ,fHistdEdx(0)
  65   ,fHistTimeBin(0)
  66 {
  67   //
  68   //  The main constructor
  69   //
  70
  71   Init();
  72
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c)
  77   :TNamed(c)
  78   ,fNMom(c.fNMom)
  79   ,fTrackMomentum(0)
  80   ,fMeanChargeRatio(c.fMeanChargeRatio)
  81   ,fNbins(c.fNbins)
  82   ,fBinSize(c.fBinSize)
  83   ,fHistdEdx(0)
  84   ,fHistTimeBin(0)
  85 {
  86   //
  87   // Copy constructor
  88   //
  89
  90   AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
  91
  92   target.fTrackMomentum = new Double_t[fNMom];
  93   for (Int_t i=0; i<fNMom; ++i) {
  94     target.fTrackMomentum[i] = fTrackMomentum[i];
  95   }
  96   if (fHistdEdx) {
  97     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
  98   }
  99
 100   if (fHistTimeBin) {
 101     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
 102   }
 103
 104 }
 105
 106 //_________________________________________________________________________
 107 AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
 108 {
 109   //
 110   // Destructor
 111   //
 112
 113   CleanUp();
 114
 115 }
 116
 117 //_________________________________________________________________________
 118 void AliTRDCalPIDLQ::CleanUp()
 119 {
 120   //
 121   // Delets all newly created objects
 122   //
 123
 124   if (fHistdEdx) {
 125     delete fHistdEdx;
 126     fHistdEdx = 0;
 127   }
 128
 129   if (fHistTimeBin) {
 130     delete fHistTimeBin;
 131     fHistTimeBin = 0;
 132   }
 133
 134   if (fTrackMomentum) {
 135     delete[] fTrackMomentum;
 136     fTrackMomentum = 0;
 137   }
 138
 139 }
 140
 141 //_____________________________________________________________________________
 142 AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
 143 {
 144   //
 145   // Assignment operator
 146   //
 147
 148   if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
 149   return *this;
 150
 151 }
 152
 153 //_____________________________________________________________________________
 154 void AliTRDCalPIDLQ::Copy(TObject &c) const
 155 {
 156   //
 157   // Copy function
 158   //
 159
 160   AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
 161
 162   target.CleanUp();
 163
 164   target.fNMom            = fNMom;
 165   target.fNbins           = fNbins;
 166   target.fBinSize         = fBinSize;
 167   target.fMeanChargeRatio = fMeanChargeRatio;
 168
 169   target.fTrackMomentum = new Double_t[fNMom];
 170   for (Int_t i=0; i<fNMom; ++i) {
 171     target.fTrackMomentum[i] = fTrackMomentum[i];
 172   }
 173
 174   if (fHistdEdx) {
 175     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
 176   }
 177   if (fHistTimeBin) {
 178     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
 179   }
 180
 181   TObject::Copy(c);
 182
 183 }
 184
 185 //_________________________________________________________________________
 186 void AliTRDCalPIDLQ::Init()
 187 {
 188   //
 189   // Initialization
 190   //
 191
 192   const Int_t kNMom = 11;
 193
 194   fNMom = kNMom;
 195   fTrackMomentum = new Double_t[fNMom];
 196   Double_t trackMomentum[kNMom] = {  0.6,  0.8,  1.0,  1.5,  2.0
 197                                   ,  3.0,  4.0,  5.0,  6.0,  8.0
 198                                   , 10.0 };
 199   for (Int_t imom = 0; imom < kNMom; imom++) {
 200     fTrackMomentum[imom] = trackMomentum[imom];
 201   }
 202
 203   fHistdEdx    = new TObjArray(AliPID::kSPECIES * fNMom);
 204   fHistdEdx->SetOwner();
 205   fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
 206   fHistTimeBin->SetOwner();
 207
 208   // ADC Gain normalization
 209   fMeanChargeRatio = 1.0;
 210
 211   // Number of bins and bin size
 212   fNbins   = 0;
 213   fBinSize = 0.0;
 214
 215 }
 216
 217 //_________________________________________________________________________
 218 Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
 219 {
 220   //
 221   // Read the TRD dEdx histograms.
 222   //
 223
 224   // Read histogram Root file
 225   TFile *histFile = new TFile(responseFile, "READ");
 226   if (!histFile || !histFile->IsOpen()) {
 227     AliError(Form("Opening TRD histgram file %s failed", responseFile));
 228     return kFALSE;
 229   }
 230   gROOT->cd();
 231
 232   // Read histograms
 233   Char_t text[200];
 234   for (Int_t particle = 0; particle < AliPID::kSPECIES; ++particle)
 235   {
 236     Char_t* particleKey = "";
 237     switch (particle)
 238     {
 239       case AliPID::kElectron: particleKey = "EL"; break;
 240       case AliPID::kPion: particleKey = "PI"; break;
 241       case AliPID::kMuon: particleKey = "MU"; break;
 242       case AliPID::kKaon: particleKey = "KA"; break;
 243       case AliPID::kProton: particleKey = "PR"; break;
 244     }
 245
 246     for (Int_t imom = 0; imom < fNMom; imom++)
 247     {
 248       sprintf(text, "h1dEdx%s%01d", particleKey, imom+1);
 249       TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
 250       hist->Scale(1.0/hist->Integral());
 251       fHistdEdx->AddAt(hist, GetHistID(particle, imom));
 252
 253       if (particle == AliPID::kElectron || particle == AliPID::kPion)
 254       {
 255         sprintf(text,"h1MaxTimBin%s%01d", particleKey, imom+1);
 256         TH1F* hist = (TH1F*)histFile->Get(text)->Clone();
 257         hist->Scale(1.0/hist->Integral());
 258         fHistTimeBin->AddAt(hist, GetHistID(particle,imom));
 259       }
 260     }
 261   }
 262
 263   histFile->Close();
 264   delete histFile;
 265
 266   // Number of bins and bin size
 267   TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
 268   fNbins   = hist->GetNbinsX();
 269   fBinSize = hist->GetBinWidth(1);
 270
 271   return kTRUE;
 272
 273 }
 274
 275 //_________________________________________________________________________
 276 Double_t  AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
 277 {
 278   //
 279   // Gets mean of de/dx dist. of e
 280   //
 281
 282   AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
 283               ,fpartName[k]
 284               ,fTrackMomentum[ip]));
 285   if (k < 0 || k > AliPID::kSPECIES) {
 286     return 0;
 287   }
 288
 289   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
 290
 291 }
 292
 293 //_________________________________________________________________________
 294 Double_t  AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
 295 {
 296   //
 297   // Gets Normalization of de/dx dist. of e
 298   //
 299
 300   AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
 301               ,fpartName[k]
 302               ,fTrackMomentum[ip]));
 303   if (k < 0 || k > AliPID::kSPECIES) {
 304     return 0;
 305   }
 306
 307   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
 308
 309 }
 310
 311 //_________________________________________________________________________
 312 TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
 313 {
 314   //
 315   // Returns one selected dEdx histogram
 316   //
 317
 318   AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
 319               ,fpartName[k]
 320               ,fTrackMomentum[ip]));
 321   if (k < 0 || k > AliPID::kSPECIES) {
 322     return 0;
 323   }
 324
 325   return (TH1F*) fHistdEdx->At(GetHistID(k,ip));
 326
 327 }
 328
 329 //_________________________________________________________________________
 330 TH1F* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
 331 {
 332   //
 333   // Returns one selected time bin max histogram
 334   //
 335
 336   AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
 337               ,fpartName[k]
 338               ,fTrackMomentum[ip]));
 339   if (k < 0 || k > AliPID::kSPECIES)
 340     return 0;
 341
 342   return (TH1F*) fHistTimeBin->At(GetHistID(k,ip));
 343
 344 }
 345
 346 //_________________________________________________________________________
 347 Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
 348 {
 349   //
 350   // Gets the Probability of having dedx at a given momentum (mom)
 351   // and particle type k (0 for e) and (2 for pi)
 352   // from the precalculated de/dx distributions
 353   //
 354
 355   Double_t dedx   = dedx1/fMeanChargeRatio;
 356   Int_t    iEnBin = ((Int_t) (dedx/fBinSize+1));
 357   if(iEnBin > fNbins) iEnBin = fNbins;
 358
 359   if (k < 0 || k > AliPID::kSPECIES) {
 360     return 1;
 361   }
 362
 363   TH1F* hist1 = 0;
 364   TH1F* hist2 = 0;
 365   Double_t mom1 = 0;
 366   Double_t mom2 = 0;
 367
 368   // Lower limit
 369   if (mom<=fTrackMomentum[0])  {
 370     hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,1));
 371     hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,0));
 372     mom1 = fTrackMomentum[1];
 373     mom2 = fTrackMomentum[0];
 374   }
 375
 376   // Upper Limit
 377   if(mom>=fTrackMomentum[fNMom-1]) {
 378     hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-1));
 379     hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,fNMom-2));
 380     mom2 = fTrackMomentum[fNMom-1];
 381     mom1 = fTrackMomentum[fNMom-2];
 382   }
 383
 384   // In the range
 385   for (Int_t ip=1; ip<fNMom; ip++) {
 386     if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 387       hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,ip));
 388       hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,ip-1));
 389       mom1 = fTrackMomentum[ip];
 390       mom2 = fTrackMomentum[ip-1];
 391     }
 392   }
 393
 394   Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin))
 395                 / (mom1 - mom2);
 396   return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);
 397
 398 }
 399
 400 //_________________________________________________________________________
 401 Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
 402 {
 403   //
 404   // Gets the Probability of having timbin at a given momentum (mom)
 405   // and particle type k (0 for e) and (2 for pi)
 406   // from the precalculated timbin distributions
 407   //
 408
 409   if (timbin<=0) {
 410     return 0.0;
 411   }
 412
 413   Int_t iTBin = timbin+1;
 414
 415   // Everything which is not an electron counts as a pion for time bin max
 416   if (k != AliPID::kElectron) {
 417     k = AliPID::kPion;
 418   }
 419
 420   if (mom<=fTrackMomentum[0]) {
 421     return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
 422   }
 423   if (mom>=fTrackMomentum[fNMom-1]) {
 424     return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
 425   }
 426
 427   for (Int_t ip=1; ip<fNMom; ip++) {
 428     if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
 429       Double_t slop = (((TH1F*) fHistTimeBin->At(GetHistID(k,ip)))->GetBinContent(iTBin)
 430                      - ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin))
 431                     / (fTrackMomentum[ip] - fTrackMomentum[ip-1]);
 432       // Linear interpolation
 433       return ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin)
 434               + slop * (mom - fTrackMomentum[ip-1]);
 435     }
 436   }
 437
 438   return -1;
 439
 440 }