7de83cf7a85004256a687b329986b5858111f3a2
[u/mrichter/AliRoot.git] / TRD / AliTRDprobdist.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 "AliTRDprobdist.h"
25
26 ClassImp(AliTRDprobdist)
27
28 //_________________________________________________________________________
29 AliTRDprobdist::AliTRDprobdist(Char_t *responseFile)
30 {
31   //
32   //  The main constructor
33   //
34   Char_t *partName[5] = {"electron", "muon", "pion", "kaon", "proton"};
35   for(Int_t i=0; i<5; i++) {
36     fpartName[i] = partName[i]; 
37   }
38   // ADC Gain normalization
39   fADCNorm=1.0;
40   fNMom=kNMom;
41   // Set track momenta for which response functions are available
42   Double_t trackMomentum[kNMom]= {0.6, 0.8, 1, 1.5, 2, 3, 4, 5, 6, 8, 10};
43   for(Int_t imom=0; imom<fNMom; imom++) {
44     fTrackMomentum[imom] = trackMomentum[imom];
45   }
46   // Read the histogram file
47   ReadData(responseFile);
48   //
49 }
50
51 Bool_t AliTRDprobdist::ReadData(Char_t *responseFile)
52 {
53   //
54   // Read the TRD dEdx histograms.
55   //
56   // Read histogram Root file  
57   TFile *histFile = new TFile(responseFile);
58   if (!histFile || !histFile->IsOpen()) {
59     Error("AliTRDprobdist", "opening TRD histgram file %s failed", responseFile);
60     return kFALSE;
61   }
62
63   // Read histograms
64   Char_t text[200];
65   for (Int_t imom = 0; imom < kNMom; imom++) {
66     sprintf(text,"h1dEdxEL%01d",imom+1);
67     fh1dEdxEL[imom] = (TH1F*)histFile->Get(text);
68     fh1dEdxEL[imom]->Scale(1.0/fh1dEdxEL[imom]->Integral());
69
70     sprintf(text,"h1dEdxPI%01d",imom+1);
71     fh1dEdxPI[imom] = (TH1F*)histFile->Get(text);
72     fh1dEdxPI[imom]->Scale(1.0/fh1dEdxPI[imom]->Integral());
73
74     sprintf(text,"h1dEdxMU%01d",imom+1);
75     fh1dEdxMU[imom] = (TH1F*)histFile->Get(text);
76     fh1dEdxMU[imom]->Scale(1.0/fh1dEdxMU[imom]->Integral());
77
78     sprintf(text,"h1dEdxKA%01d",imom+1);
79     fh1dEdxKA[imom] = (TH1F*)histFile->Get(text);
80     fh1dEdxKA[imom]->Scale(1.0/fh1dEdxKA[imom]->Integral());
81
82     sprintf(text,"h1dEdxPR%01d",imom+1);
83     fh1dEdxPR[imom] = (TH1F*)histFile->Get(text);
84     fh1dEdxPR[imom]->Scale(1.0/fh1dEdxPR[imom]->Integral());
85
86     sprintf(text,"h1MaxTimBinEL%01d",imom+1);
87     fh1MaxTimBinEL[imom] = (TH1F*)histFile->Get(text);
88     fh1MaxTimBinEL[imom]->Scale(1.0/fh1MaxTimBinEL[imom]->Integral());
89
90     sprintf(text,"h1MaxTimBinPI%01d",imom+1);
91     fh1MaxTimBinPI[imom] = (TH1F*)histFile->Get(text);
92     fh1MaxTimBinPI[imom]->Scale(1.0/fh1MaxTimBinPI[imom]->Integral());
93   }
94   // Number of bins and bin size
95   fNbins = fh1dEdxPI[1]->GetNbinsX();
96   fBinSize = fh1dEdxPI[1]->GetBinWidth(1);
97   return kTRUE;
98 }
99
100 AliTRDprobdist::AliTRDprobdist(const AliTRDprobdist& pd):TNamed()
101 {
102   //
103   // Copy constructor.
104   //
105   for(Int_t i=0; i<5; i++) {
106     fpartName[i] = pd.fpartName[i]; 
107   }
108   // ADC Gain normalization
109   fADCNorm=pd.fADCNorm;
110   fNMom=pd.fNMom;
111   // Set track momenta for which response functions are available
112   for(Int_t imom=0; imom<fNMom; imom++) {
113     fTrackMomentum[imom] = pd.fTrackMomentum[imom];
114   }
115
116   fNbins=pd.fNbins;
117   fBinSize=pd.fBinSize;  
118
119   for (Int_t imom = 0; imom < kNMom; imom++) {
120     fh1dEdxEL[imom] = pd.fh1dEdxEL[imom];
121     fh1dEdxPI[imom] = pd.fh1dEdxPI[imom];
122     fh1dEdxMU[imom] = pd.fh1dEdxMU[imom];
123     fh1dEdxKA[imom] = pd.fh1dEdxKA[imom];
124     fh1dEdxPR[imom] = pd.fh1dEdxPR[imom];
125     fh1MaxTimBinEL[imom] = pd.fh1MaxTimBinEL[imom];
126     fh1MaxTimBinPI[imom] = pd.fh1MaxTimBinPI[imom];
127   }
128 }
129
130
131 AliTRDprobdist::~AliTRDprobdist()
132 {
133   // Destructor
134
135 }
136
137
138 //_________________________________________________________________________
139 Double_t  AliTRDprobdist::GetMean(Int_t k, Int_t ip) const
140 {
141   //
142   // Gets mean of de/dx dist. of e
143   printf("Mean for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
144   if(k==0) return fh1dEdxEL[ip]->GetMean();
145   if(k==1) return fh1dEdxMU[ip]->GetMean();
146   if(k==2) return fh1dEdxPI[ip]->GetMean();
147   if(k==3) return fh1dEdxKA[ip]->GetMean();
148   if(k==4) return fh1dEdxPR[ip]->GetMean();
149   return fh1dEdxPR[ip]->GetMean();
150 }
151
152 //_________________________________________________________________________
153 Double_t  AliTRDprobdist::GetNormalization(Int_t k, Int_t ip) const
154 {
155   //
156   // Gets Normalization of de/dx dist. of e
157
158   printf("Normalization for particle = %s and momentum = %.2f is:\n",fpartName[k], fTrackMomentum[ip]);
159   if(k==0) return fh1dEdxEL[ip]->Integral();
160   if(k==1) return fh1dEdxMU[ip]->Integral();
161   if(k==2) return fh1dEdxPI[ip]->Integral();
162   if(k==3) return fh1dEdxKA[ip]->Integral();
163   if(k==4) return fh1dEdxPR[ip]->Integral();
164   return fh1dEdxPR[ip]->Integral();
165 }
166
167 TH1F* AliTRDprobdist::GetHistogram(Int_t k, Int_t ip) const
168 {
169   //
170   //
171   printf("Histogram for particle = %s and momentum = %.2f is:\n", fpartName[k], fTrackMomentum[ip]);
172   if(k==0) return fh1dEdxEL[ip];
173   if(k==1) return fh1dEdxMU[ip];
174   if(k==2) return fh1dEdxPI[ip];
175   if(k==3) return fh1dEdxKA[ip];
176   if(k==4) return fh1dEdxPR[ip];
177   return fh1dEdxPR[ip];
178 }
179
180 //_________________________________________________________________________
181 Double_t AliTRDprobdist::GetProbability(Int_t k, Double_t mom, Double_t dedx1) const
182 {
183   //
184   // Gets the Probability of having dedx at a given momentum (mom)
185   // and particle type k (0 for e) and (2 for pi)
186   // from the precalculated de/dx distributions 
187   Double_t probability = 1.0;
188   Double_t dedx = dedx1/fADCNorm;
189   Int_t iEnBin= ((Int_t) (dedx/fBinSize+1));
190   if(iEnBin > fNbins) iEnBin = fNbins;
191
192   ////Electron//////////////////////////
193   if(k==0){    // electron 
194     Double_t slop;
195     // Lower limit
196     if(mom<=fTrackMomentum[0]) {
197       slop=(fh1dEdxEL[1]->GetBinContent(iEnBin)-fh1dEdxEL[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
198       probability= fh1dEdxEL[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
199       return probability;
200     }
201     // Upper Limit
202     if(mom>=fTrackMomentum[fNMom-1]) {
203       slop=(fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxEL[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
204       probability= fh1dEdxEL[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
205       return probability;
206     }
207     // In the range
208     for(Int_t ip=1; ip<fNMom; ip++) 
209       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
210         slop=(fh1dEdxEL[ip]->GetBinContent(iEnBin)-fh1dEdxEL[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
211         // Linear Interpolation
212         probability= fh1dEdxEL[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
213         return probability;
214       }
215   }
216
217   ////Pion//////////////////////////
218   if(k==2){    // Pion
219     Double_t slop;
220     // Lower limit
221     if(mom<=fTrackMomentum[0]) {
222       slop=(fh1dEdxPI[1]->GetBinContent(iEnBin)-fh1dEdxPI[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
223       probability= fh1dEdxPI[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
224       return probability;
225     }
226     // Upper Limit
227     if(mom>=fTrackMomentum[fNMom-1]) {
228       slop=(fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPI[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
229       probability= fh1dEdxPI[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
230       return probability;
231     }
232     // In the range
233     for(Int_t ip=1; ip<fNMom; ip++) 
234       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
235         slop=(fh1dEdxPI[ip]->GetBinContent(iEnBin)-fh1dEdxPI[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
236         // Linear Interpolation
237         probability= fh1dEdxPI[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
238         return probability;
239       }
240   }
241
242   ////Muon//////////////////////////
243   if(k==1){    // Muon
244     Double_t slop;
245     // Lower limit
246     if(mom<=fTrackMomentum[0]) {
247       slop=(fh1dEdxMU[1]->GetBinContent(iEnBin)-fh1dEdxMU[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
248       probability= fh1dEdxMU[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
249       return probability;
250     }
251     // Upper Limit
252     if(mom>=fTrackMomentum[fNMom-1]) {
253       slop=(fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxMU[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
254       probability= fh1dEdxMU[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
255       return probability;
256     }
257     // In the range
258     for(Int_t ip=1; ip<fNMom; ip++) 
259       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
260         slop=(fh1dEdxMU[ip]->GetBinContent(iEnBin)-fh1dEdxMU[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
261         // Linear Interpolation
262         probability= fh1dEdxMU[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
263         return probability;
264       }
265   }
266
267   ////Kaon//////////////////////////
268   if(k==3){    // Kaon
269     Double_t slop;
270     // Lower limit
271     if(mom<=fTrackMomentum[0]) {
272       slop=(fh1dEdxKA[1]->GetBinContent(iEnBin)-fh1dEdxKA[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
273       probability= fh1dEdxKA[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
274       return probability;
275     }
276     // Upper Limit
277     if(mom>=fTrackMomentum[fNMom-1]) {
278       slop=(fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxKA[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
279       probability= fh1dEdxKA[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
280       return probability;
281     }
282     // In the range
283     for(Int_t ip=1; ip<fNMom; ip++) 
284       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
285         slop=(fh1dEdxKA[ip]->GetBinContent(iEnBin)-fh1dEdxKA[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
286         // Linear Interpolation
287         probability= fh1dEdxKA[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
288         return probability;
289       }
290   }
291
292   ////Proton//////////////////////////
293   if(k==4){    // Proton
294     Double_t slop;
295     // Lower limit
296     if(mom<=fTrackMomentum[0]) {
297       slop=(fh1dEdxPR[1]->GetBinContent(iEnBin)-fh1dEdxPR[0]->GetBinContent(iEnBin))/(fTrackMomentum[1] - fTrackMomentum[0]);
298       probability= fh1dEdxPR[0]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[0]);
299       return probability;
300     }
301     // Upper Limit
302     if(mom>=fTrackMomentum[fNMom-1]) {
303       slop=(fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin)-fh1dEdxPR[fNMom-2]->GetBinContent(iEnBin))/(fTrackMomentum[fNMom-1] - fTrackMomentum[fNMom-2]);
304       probability= fh1dEdxPR[fNMom-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[fNMom-1]);
305       return probability;
306     }
307     // In the range
308     for(Int_t ip=1; ip<fNMom; ip++) 
309       if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
310         slop=(fh1dEdxPR[ip]->GetBinContent(iEnBin)-fh1dEdxPR[ip-1]->GetBinContent(iEnBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
311         // Linear Interpolation
312         probability= fh1dEdxPR[ip-1]->GetBinContent(iEnBin) + slop*(mom-fTrackMomentum[ip-1]);
313         return probability;
314       }
315   }
316
317   return probability;
318 }
319
320
321 //_________________________________________________________________________
322 Double_t AliTRDprobdist::GetProbabilityT(Int_t k, Double_t mom, Int_t timbin) const
323 {
324   //
325   // Gets the Probability of having timbin at a given momentum (mom)
326   // and particle type k (0 for e) and (2 for pi)
327   // from the precalculated timbin distributions 
328   Double_t probabilityT = 1.0;
329   if(timbin<=0) return 0.;
330   Int_t iTBin=timbin+1;
331
332   if(k==0){    // electron
333     if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinEL[0]->GetBinContent(iTBin);
334     if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinEL[fNMom-1]->GetBinContent(iTBin);
335   }
336   if(k==1||k==2||k==3||k==4){    // pion
337     if(mom<=fTrackMomentum[0]) probabilityT = fh1MaxTimBinPI[0]->GetBinContent(iTBin);
338     if(mom>=fTrackMomentum[fNMom-1]) probabilityT = fh1MaxTimBinPI[fNMom-1]->GetBinContent(iTBin);
339   }
340
341   if(k==0)    // electron
342   for(Int_t ip=1; ip<fNMom; ip++)
343   if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
344     Double_t slop=(fh1MaxTimBinEL[ip]->GetBinContent(iTBin)-fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
345     // Linear Interpolation
346     probabilityT= fh1MaxTimBinEL[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
347     return probabilityT;
348   }
349
350   if(k==1||k==2||k==3||k==4)   // pion and other particles
351   for(Int_t ip=1; ip<fNMom; ip++)
352   if((fTrackMomentum[ip-1]<= mom) && (mom<fTrackMomentum[ip])) {
353     Double_t slop=(fh1MaxTimBinPI[ip]->GetBinContent(iTBin)-fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin))/(fTrackMomentum[ip] - fTrackMomentum[ip-1]);
354     // Linear Interpolation
355     probabilityT= fh1MaxTimBinPI[ip-1]->GetBinContent(iTBin) + slop*(mom-fTrackMomentum[ip-1]);
356     return probabilityT;
357   }
358   return probabilityT;
359 }
360
361