1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------
17 // Class for dE/dx and Time Bin of Max. Cluster for Electrons and
19 // It is instantiated in class AliTRDpidESD for particle identification
21 // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>
22 //-----------------------------------------------------------------
29 #include "AliTRDCalPIDLQ.h"
31 ClassImp(AliTRDCalPIDLQ)
33 Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
35 //_________________________________________________________________________
36 AliTRDCalPIDLQ::AliTRDCalPIDLQ():
47 // The Default constructor
52 //_________________________________________________________________________
53 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) : TNamed(name, title)
56 // The main constructor
63 //_____________________________________________________________________________
64 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) : TNamed(c)
72 ((AliTRDCalPIDLQ &) c).Copy(*this);
76 //_________________________________________________________________________
77 AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
87 //_________________________________________________________________________
88 void AliTRDCalPIDLQ::CleanUp()
91 // Delets all newly created objects
108 delete[] fTrackMomentum;
114 //_____________________________________________________________________________
115 AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
118 // Assignment operator
121 if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
126 //_____________________________________________________________________________
127 void AliTRDCalPIDLQ::Copy(TObject &c) const
133 AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
137 target.fNMom = fNMom;
139 target.fTrackMomentum = new Double_t[fNMom];
140 for (Int_t i=0; i<fNMom; ++i)
141 target.fTrackMomentum[i] = fTrackMomentum[i];
143 target.fMeanChargeRatio = fMeanChargeRatio;
145 target.fNbins = fNbins;
146 target.fBinSize = fBinSize;
149 target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
152 target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
158 //_________________________________________________________________________
159 void AliTRDCalPIDLQ::Init()
165 const Int_t kNMom = 11;
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
172 for (Int_t imom = 0; imom < kNMom; imom++) {
173 fTrackMomentum[imom] = trackMomentum[imom];
176 fHistdEdx = new TObjArray(AliPID::kSPECIES * fNMom);
177 fHistdEdx->SetOwner();
178 fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
179 fHistTimeBin->SetOwner();
181 // ADC Gain normalization
182 fMeanChargeRatio = 1.0;
184 // Number of bins and bin size
190 //_________________________________________________________________________
191 Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
194 // Read the TRD dEdx histograms.
197 // Read histogram Root file
198 TFile *histFile = new TFile(responseFile, "READ");
199 if (!histFile || !histFile->IsOpen()) {
201 error.Form("Opening TRD histgram file %s failed", responseFile);
209 for (Int_t particle = 0; particle < AliPID::kSPECIES; ++particle)
211 Char_t* particleKey = "";
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;
221 for (Int_t imom = 0; imom < fNMom; imom++)
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));
228 if (particle == AliPID::kElectron || particle == AliPID::kPion)
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));
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);
250 //_________________________________________________________________________
251 Double_t AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
254 // Gets mean of de/dx dist. of e
257 printf("Mean for particle = %s and momentum = %.2f is:\n"
259 ,fTrackMomentum[ip]);
260 if (k < 0 || k > AliPID::kSPECIES)
263 return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
267 //_________________________________________________________________________
268 Double_t AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
271 // Gets Normalization of de/dx dist. of e
274 printf("Normalization for particle = %s and momentum = %.2f is:\n"
276 ,fTrackMomentum[ip]);
277 if (k < 0 || k > AliPID::kSPECIES)
280 return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
284 //_________________________________________________________________________
285 TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
288 // Returns one selected dEdx histogram
291 printf("Histogram for particle = %s and momentum = %.2f is:\n"
293 ,fTrackMomentum[ip]);
294 if (k < 0 || k > AliPID::kSPECIES)
297 return (TH1F*) fHistdEdx->At(GetHistID(k,ip));
301 //_________________________________________________________________________
302 TH1F* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
305 // Returns one selected time bin max histogram
308 printf("Histogram for particle = %s and momentum = %.2f is:\n"
310 ,fTrackMomentum[ip]);
311 if (k < 0 || k > AliPID::kSPECIES)
314 return (TH1F*) fHistTimeBin->At(GetHistID(k,ip));
318 //_________________________________________________________________________
319 Double_t AliTRDCalPIDLQ::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
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
327 Double_t dedx = dedx1/fMeanChargeRatio;
328 Int_t iEnBin= ((Int_t) (dedx/fBinSize+1));
329 if(iEnBin > fNbins) iEnBin = fNbins;
331 if (k < 0 || k > AliPID::kSPECIES)
340 if (mom<=fTrackMomentum[0])
342 hist1 = (TH1F*) fHistdEdx->At(GetHistID(k,1));
343 hist2 = (TH1F*) fHistdEdx->At(GetHistID(k,0));
344 mom1 = fTrackMomentum[1];
345 mom2 = fTrackMomentum[0];
349 if(mom>=fTrackMomentum[fNMom-1])
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];
358 for (Int_t ip=1; ip<fNMom; ip++)
360 if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip]))
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];
369 Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin))
371 return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);
375 //_________________________________________________________________________
376 Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
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
386 Int_t iTBin=timbin+1;
388 // everything which is not electron counts as pion for time bin
389 if (k != AliPID::kElectron)
392 if (mom<=fTrackMomentum[0])
393 return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
395 if (mom>=fTrackMomentum[fNMom-1])
396 return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
398 for (Int_t ip=1; ip<fNMom; ip++)
400 if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip]))
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]);