First round of effc++ changes
[u/mrichter/AliRoot.git] / 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 }