c82b17f2be670b602908b427fd12c68754e04957
[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 // 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 }