TRD/Cal/AliTRDCalPID.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 // Container for the distributions of dE/dx and the time bin of the     //
  21 // max. cluster for electrons and pions                                 //
  22 //                                                                      //
  23 // Authors:                                                             //
  24 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
  25 //   Alex Bercuci <a.bercuci@gsi.de>                                    //
  26 //                                                                      //
  27 //////////////////////////////////////////////////////////////////////////
  28
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TFile.h>
  32 #include <TROOT.h>
  33
  34 #include "AliLog.h"
  35 #include "AliPID.h"
  36 #include "AliESD.h"
  37 #include "AliESDtrack.h"
  38
  39 #include "AliTRDCalPID.h"
  40 #include "AliTRDcalibDB.h"
  41
  42 ClassImp(AliTRDCalPID)
  43
  44 Char_t* AliTRDCalPID::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
  45
  46 Char_t* AliTRDCalPID::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
  47
  48 Float_t AliTRDCalPID::fTrackMomentum[kNMom] = {0.6,  0.8,  1.0,  1.5,  2.0,  3.0,  4.0,  5.0,  6.0,  8.0,  10.0};
  49
  50 Float_t AliTRDCalPID::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};
  51
  52
  53 //_________________________________________________________________________
  54 AliTRDCalPID::AliTRDCalPID()
  55   :TNamed("pid", "PID for TRD")
  56   ,fMeanChargeRatio(0)
  57   ,fHistdEdx(0x0)
  58   ,fHistTimeBin(0x0)
  59 {
  60   //
  61   //  The Default constructor
  62   //
  63
  64   Init();
  65
  66 }
  67
  68 //_________________________________________________________________________
  69 AliTRDCalPID::AliTRDCalPID(const Text_t *name, const Text_t *title)
  70   :TNamed(name,title)
  71   ,fMeanChargeRatio(0)
  72   ,fHistdEdx(0x0)
  73   ,fHistTimeBin(0x0)
  74 {
  75   //
  76   //  The main constructor
  77   //
  78
  79   Init();
  80
  81 }
  82
  83 //_____________________________________________________________________________
  84 AliTRDCalPID::AliTRDCalPID(const AliTRDCalPID &c)
  85   :TNamed(c)
  86   ,fMeanChargeRatio(c.fMeanChargeRatio)
  87   ,fHistdEdx(0x0)
  88   ,fHistTimeBin(0x0)
  89 {
  90   //
  91   // Copy constructor
  92   //
  93
  94   if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
  95
  96 }
  97
  98 //_________________________________________________________________________
  99 AliTRDCalPID::~AliTRDCalPID()
 100 {
 101   //
 102   // Destructor
 103   //
 104
 105   CleanUp();
 106
 107 }
 108
 109 //_________________________________________________________________________
 110 void AliTRDCalPID::CleanUp()
 111 {
 112   //
 113   // Delets all newly created objects
 114   //
 115
 116   if (fHistdEdx) {
 117     delete fHistdEdx;
 118     fHistdEdx = 0x0;
 119   }
 120
 121   if (fHistTimeBin) {
 122     delete fHistTimeBin;
 123     fHistTimeBin = 0x0;
 124   }
 125 }
 126
 127 //_____________________________________________________________________________
 128 AliTRDCalPID &AliTRDCalPID::operator=(const AliTRDCalPID &c)
 129 {
 130   //
 131   // Assignment operator
 132   //
 133
 134   if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
 135   return *this;
 136
 137 }
 138
 139 //_____________________________________________________________________________
 140 void AliTRDCalPID::Copy(TObject &c) const
 141 {
 142   //
 143   // Copy function
 144   //
 145
 146   AliTRDCalPID& target = (AliTRDCalPID &) c;
 147
 148   target.CleanUp();
 149
 150   target.fMeanChargeRatio = fMeanChargeRatio;
 151
 152   if (fHistdEdx) {
 153     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
 154   }
 155   if (fHistTimeBin) {
 156     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
 157   }
 158
 159   TObject::Copy(c);
 160
 161 }
 162
 163 //_________________________________________________________________________
 164 void AliTRDCalPID::Init()
 165 {
 166   //
 167   // Initialization
 168   //
 169
 170   fHistdEdx    = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
 171   fHistdEdx->SetOwner();
 172   fHistTimeBin = new TObjArray(2 * kNMom);
 173   fHistTimeBin->SetOwner();
 174
 175         // Initialization of estimator at object instantiation because late
 176         // initialization in function GetProbability() is not working due to
 177         // constantness of this function.
 178         // fEstimator = new AliTRDCalPIDRefMaker();
 179
 180   // ADC Gain normalization
 181   fMeanChargeRatio = 1.0;
 182 }
 183
 184 //_________________________________________________________________________
 185 Bool_t AliTRDCalPID::LoadLQReferences(Char_t *refFile)
 186 {
 187   //
 188   // Read the TRD dEdx histograms.
 189   //
 190
 191         Int_t nTimeBins = 22;
 192         // Get number of time bins from CDB
 193         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 194         if(!calibration){
 195                 AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
 196         }else{
 197                 if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
 198                 else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
 199         }
 200
 201
 202   // Read histogram Root file
 203   TFile *histFile = TFile::Open(refFile, "READ");
 204   if (!histFile || !histFile->IsOpen()) {
 205     AliError(Form("Opening TRD histgram file %s failed", refFile));
 206     return kFALSE;
 207   }
 208   gROOT->cd();
 209
 210   // Read histograms
 211   for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
 212     for (Int_t imom = 0; imom < kNMom; imom++){
 213                         TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
 214                         hist->Scale(1./hist->Integral());
 215                         fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
 216
 217                         if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
 218
 219                         TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
 220                         if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, nTimeBins));
 221                         ht->Scale(1./ht->Integral());
 222                         fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
 223                 }
 224   }
 225
 226   histFile->Close();
 227   delete histFile;
 228
 229   // Number of bins and bin size
 230   //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
 231   //fNbins   = hist->GetNbinsX();
 232   //fBinSize = hist->GetBinWidth(1);
 233
 234   return kTRUE;
 235
 236 }
 237
 238 // //_________________________________________________________________________
 239 // Double_t  AliTRDCalPID::GetMean(Int_t k, Int_t ip) const
 240 // {
 241 //   //
 242 //   // Gets mean of de/dx dist. of e
 243 //   //
 244 //
 245 //   AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
 246 //               ,fpartName[k]
 247 //               ,fTrackMomentum[ip]));
 248 //   if (k < 0 || k > AliPID::kSPECIES) {
 249 //     return 0;
 250 //   }
 251 //
 252 //   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
 253 //
 254 // }
 255 //
 256 // //_________________________________________________________________________
 257 // Double_t  AliTRDCalPID::GetNormalization(Int_t k, Int_t ip) const
 258 // {
 259 //   //
 260 //   // Gets Normalization of de/dx dist. of e
 261 //   //
 262 //
 263 //   AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
 264 //               ,fpartName[k]
 265 //               ,fTrackMomentum[ip]));
 266 //   if (k < 0 || k > AliPID::kSPECIES) {
 267 //     return 0;
 268 //   }
 269 //
 270 //   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
 271 //
 272 // }
 273
 274 //_________________________________________________________________________
 275 TH1* AliTRDCalPID::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
 276 {
 277   //
 278   // Returns one selected dEdx histogram
 279   //
 280
 281   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
 282         if(ip<0 || ip>= kNMom ) return 0x0;
 283
 284         AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
 285
 286   return (TH1*)fHistdEdx->At(GetHistID(k, ip));
 287
 288 }
 289
 290 //_________________________________________________________________________
 291 TH1* AliTRDCalPID::GetHistogramT(Int_t k, Int_t ip) const
 292 {
 293   //
 294   // Returns one selected time bin max histogram
 295   //
 296
 297   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
 298         if(ip<0 || ip>= kNMom ) return 0x0;
 299
 300         AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
 301
 302         return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*kNMom+ip);
 303 }
 304
 305
 306
 307 //_________________________________________________________________________
 308 Double_t AliTRDCalPID::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const
 309 {
 310   //
 311         // Core function of AliTRDCalPID class for calculating the
 312         // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
 313         // given momentum "mom" and a given dE/dx (slice "dedx") yield per
 314         // layer
 315   //
 316
 317         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 318
 319         //Double_t dedx   = dedx1/fMeanChargeRatio;
 320
 321         // find the interval in momentum and track segment length which applies for this data
 322         Int_t ilength = 1;
 323   while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
 324                 ilength++;
 325         }
 326         Int_t imom = 1;
 327   while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
 328
 329         Int_t nbinsx, nbinsy;
 330         TAxis *ax = 0x0, *ay = 0x0;
 331         Double_t LQ1, LQ2;
 332         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 333         TH2 *hist = 0x0;
 334         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom-1/*, ilength*/)))){
 335                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 336                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom-1)));
 337                 return 0.;
 338         }
 339         ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
 340         ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
 341         Float_t x = dedx[0]+dedx[1], y = dedx[2];
 342   Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
 343         Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
 344         if(kX)
 345                 if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y));
 346     //fEstimator->Estimate2D2(hist, x, y);
 347                 else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 348         else
 349                 if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 350                 else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
 351
 352
 353         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
 354                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 355                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
 356                 return LQ1;
 357         }
 358         if(kX)
 359                 if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y));
 360     //fEstimator->Estimate2D2(hist, x, y);
 361                 else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 362         else
 363                 if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 364                 else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
 365
 366
 367         // return interpolation over momentum binning
 368   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 369
 370 }
 371
 372 //_________________________________________________________________________
 373 Double_t AliTRDCalPID::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
 374 {
 375   //
 376   // Gets the Probability of having timbin at a given momentum (mom)
 377   // and particle type (spec) (0 for e) and (2 for pi)
 378   // from the precalculated timbin distributions
 379   //
 380
 381         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 382
 383   Int_t iTBin = timbin+1;
 384
 385   // Everything which is not an electron counts as a pion for time bin max
 386   //if(spec != AliPID::kElectron) spec = AliPID::kPion;
 387
 388
 389
 390         Int_t imom = 1;
 391   while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
 392
 393         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 394         TH1F *hist = 0x0;
 395         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom-1))){
 396                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
 397                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom-1));
 398                 return 0.;
 399         }
 400         Double_t LQ1 = hist->GetBinContent(iTBin);
 401
 402         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom))){
 403                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
 404                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom));
 405                 return LQ1;
 406         }
 407         Double_t LQ2 = hist->GetBinContent(iTBin);
 408
 409         // return interpolation over momentum binning
 410   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 411 }
 412