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