First round of effc++ changes
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPIDLQ.cxx
CommitLineData
7754cd1f 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
2745a409 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///////////////////////////////////////////////////////////////////////////////
7754cd1f 26
7754cd1f 27#include <TH1F.h>
28#include <TFile.h>
29
b92cdef5 30#include "AliLog.h"
2745a409 31#include "AliPID.h"
b92cdef5 32
33#include "AliTRDCalPIDLQ.h"
34
7754cd1f 35ClassImp(AliTRDCalPIDLQ)
36
37Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
38
39//_________________________________________________________________________
2745a409 40AliTRDCalPIDLQ::AliTRDCalPIDLQ()
41 :TNamed()
42 ,fNMom(0)
43 ,fTrackMomentum(0)
44 ,fMeanChargeRatio(0)
45 ,fNbins(0)
46 ,fBinSize(0)
47 ,fHistdEdx(0)
48 ,fHistTimeBin(0)
7754cd1f 49{
50 //
51 // The Default constructor
52 //
53
54}
55
56//_________________________________________________________________________
2745a409 57AliTRDCalPIDLQ::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)
7754cd1f 66{
67 //
68 // The main constructor
69 //
70
71 Init();
b92cdef5 72
7754cd1f 73}
74
75//_____________________________________________________________________________
2745a409 76AliTRDCalPIDLQ::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)
7754cd1f 85{
86 //
b92cdef5 87 // Copy constructor
7754cd1f 88 //
89
2745a409 90 AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
7754cd1f 91
2745a409 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 }
b92cdef5 103
7754cd1f 104}
105
106//_________________________________________________________________________
107AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
108{
109 //
110 // Destructor
111 //
112
113 CleanUp();
b92cdef5 114
7754cd1f 115}
116
117//_________________________________________________________________________
118void AliTRDCalPIDLQ::CleanUp()
119{
b92cdef5 120 //
121 // Delets all newly created objects
122 //
123
2745a409 124 if (fHistdEdx) {
7754cd1f 125 delete fHistdEdx;
126 fHistdEdx = 0;
127 }
128
2745a409 129 if (fHistTimeBin) {
7754cd1f 130 delete fHistTimeBin;
131 fHistTimeBin = 0;
132 }
133
2745a409 134 if (fTrackMomentum) {
7754cd1f 135 delete[] fTrackMomentum;
136 fTrackMomentum = 0;
137 }
b92cdef5 138
7754cd1f 139}
140
141//_____________________________________________________________________________
142AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
143{
144 //
145 // Assignment operator
146 //
147
148 if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
149 return *this;
b92cdef5 150
7754cd1f 151}
152
153//_____________________________________________________________________________
154void AliTRDCalPIDLQ::Copy(TObject &c) const
155{
156 //
157 // Copy function
158 //
159
160 AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
161
162 target.CleanUp();
163
2745a409 164 target.fNMom = fNMom;
165 target.fNbins = fNbins;
166 target.fBinSize = fBinSize;
167 target.fMeanChargeRatio = fMeanChargeRatio;
7754cd1f 168
169 target.fTrackMomentum = new Double_t[fNMom];
2745a409 170 for (Int_t i=0; i<fNMom; ++i) {
7754cd1f 171 target.fTrackMomentum[i] = fTrackMomentum[i];
2745a409 172 }
7754cd1f 173
2745a409 174 if (fHistdEdx) {
7754cd1f 175 target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
2745a409 176 }
177 if (fHistTimeBin) {
7754cd1f 178 target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
2745a409 179 }
180
7754cd1f 181 TObject::Copy(c);
b92cdef5 182
7754cd1f 183}
184
185//_________________________________________________________________________
186void AliTRDCalPIDLQ::Init()
187{
b92cdef5 188 //
189 // Initialization
190 //
191
192 const Int_t kNMom = 11;
193
194 fNMom = kNMom;
7754cd1f 195 fTrackMomentum = new Double_t[fNMom];
b92cdef5 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++) {
7754cd1f 200 fTrackMomentum[imom] = trackMomentum[imom];
b92cdef5 201 }
202
203 fHistdEdx = new TObjArray(AliPID::kSPECIES * fNMom);
7754cd1f 204 fHistdEdx->SetOwner();
205 fHistTimeBin = new TObjArray(AliPID::kSPECIES * fNMom);
206 fHistTimeBin->SetOwner();
207
208 // ADC Gain normalization
b92cdef5 209 fMeanChargeRatio = 1.0;
7754cd1f 210
211 // Number of bins and bin size
b92cdef5 212 fNbins = 0;
213 fBinSize = 0.0;
214
7754cd1f 215}
216
217//_________________________________________________________________________
218Bool_t AliTRDCalPIDLQ::ReadData(Char_t *responseFile)
219{
220 //
221 // Read the TRD dEdx histograms.
222 //
b92cdef5 223
7754cd1f 224 // Read histogram Root file
225 TFile *histFile = new TFile(responseFile, "READ");
226 if (!histFile || !histFile->IsOpen()) {
2745a409 227 AliError(Form("Opening TRD histgram file %s failed", responseFile));
7754cd1f 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;
b92cdef5 272
7754cd1f 273}
274
275//_________________________________________________________________________
276Double_t AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
277{
278 //
279 // Gets mean of de/dx dist. of e
b92cdef5 280 //
281
2745a409 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) {
7754cd1f 286 return 0;
2745a409 287 }
288
7754cd1f 289 return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
b92cdef5 290
7754cd1f 291}
292
293//_________________________________________________________________________
294Double_t AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
295{
296 //
297 // Gets Normalization of de/dx dist. of e
b92cdef5 298 //
7754cd1f 299
2745a409 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) {
7754cd1f 304 return 0;
2745a409 305 }
7754cd1f 306
307 return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
b92cdef5 308
7754cd1f 309}
310
311//_________________________________________________________________________
312TH1F* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip) const
313{
314 //
b92cdef5 315 // Returns one selected dEdx histogram
7754cd1f 316 //
b92cdef5 317
2745a409 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) {
7754cd1f 322 return 0;
2745a409 323 }
7754cd1f 324
325 return (TH1F*) fHistdEdx->At(GetHistID(k,ip));
b92cdef5 326
327}
328
329//_________________________________________________________________________
330TH1F* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
331{
332 //
333 // Returns one selected time bin max histogram
334 //
335
2745a409 336 AliInfo(Form("Histogram for particle = %s and momentum = %.2f is:\n"
337 ,fpartName[k]
338 ,fTrackMomentum[ip]));
b92cdef5 339 if (k < 0 || k > AliPID::kSPECIES)
340 return 0;
341
342 return (TH1F*) fHistTimeBin->At(GetHistID(k,ip));
343
7754cd1f 344}
345
346//_________________________________________________________________________
347Double_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
b92cdef5 353 //
7754cd1f 354
2745a409 355 Double_t dedx = dedx1/fMeanChargeRatio;
356 Int_t iEnBin = ((Int_t) (dedx/fBinSize+1));
7754cd1f 357 if(iEnBin > fNbins) iEnBin = fNbins;
358
2745a409 359 if (k < 0 || k > AliPID::kSPECIES) {
7754cd1f 360 return 1;
2745a409 361 }
7754cd1f 362
363 TH1F* hist1 = 0;
364 TH1F* hist2 = 0;
365 Double_t mom1 = 0;
366 Double_t mom2 = 0;
367
368 // Lower limit
2745a409 369 if (mom<=fTrackMomentum[0]) {
7754cd1f 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
2745a409 377 if(mom>=fTrackMomentum[fNMom-1]) {
7754cd1f 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
2745a409 385 for (Int_t ip=1; ip<fNMom; ip++) {
386 if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
7754cd1f 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
b92cdef5 394 Double_t slop = (hist1->GetBinContent(iEnBin) - hist2->GetBinContent(iEnBin))
395 / (mom1 - mom2);
7754cd1f 396 return hist2->GetBinContent(iEnBin) + slop * (mom - mom2);
b92cdef5 397
7754cd1f 398}
399
400//_________________________________________________________________________
401Double_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
b92cdef5 407 //
7754cd1f 408
2745a409 409 if (timbin<=0) {
410 return 0.0;
411 }
412
413 Int_t iTBin = timbin+1;
7754cd1f 414
2745a409 415 // Everything which is not an electron counts as a pion for time bin max
416 if (k != AliPID::kElectron) {
7754cd1f 417 k = AliPID::kPion;
2745a409 418 }
7754cd1f 419
2745a409 420 if (mom<=fTrackMomentum[0]) {
7754cd1f 421 return ((TH1F*) fHistTimeBin->At(GetHistID(k,0)))->GetBinContent(iTBin);
2745a409 422 }
423 if (mom>=fTrackMomentum[fNMom-1]) {
7754cd1f 424 return ((TH1F*) fHistTimeBin->At(GetHistID(k,fNMom-1)))->GetBinContent(iTBin);
2745a409 425 }
426
427 for (Int_t ip=1; ip<fNMom; ip++) {
428 if ((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
b92cdef5 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]);
2745a409 432 // Linear interpolation
b92cdef5 433 return ((TH1F*) fHistTimeBin->At(GetHistID(k,ip-1)))->GetBinContent(iTBin)
434 + slop * (mom - fTrackMomentum[ip-1]);
7754cd1f 435 }
436 }
437
438 return -1;
b92cdef5 439
7754cd1f 440}