Rename AliTRDCalPIDLQ
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPID.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 /* $Id$ */
17
18 //////////////////////////////////////////////////////////////////////////
19 //                                                                      //
20 // Container for the distributions of dE/dx and the time bin of the     //
21 // max. cluster for electrons and pions                                 //
22 //                                                                      //
23 // Authors:                                                             //
24 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
25 //   Alex Bercuci <a.bercuci@gsi.de>                                    //
26 //                                                                      //
27 //////////////////////////////////////////////////////////////////////////
28
29 #include <TH1F.h>
30 #include <TH2F.h>
31 #include <TFile.h>
32 #include <TROOT.h>
33
34 #include "AliLog.h"
35 #include "AliPID.h"
36 #include "AliESD.h"
37 #include "AliESDtrack.h"
38
39 #include "AliTRDCalPID.h"
40 #include "AliTRDcalibDB.h"
41
42 ClassImp(AliTRDCalPID)
43
44 Char_t* AliTRDCalPID::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
45
46 Char_t* AliTRDCalPID::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
47
48 Float_t AliTRDCalPID::fTrackMomentum[kNMom] = {0.6,  0.8,  1.0,  1.5,  2.0,  3.0,  4.0,  5.0,  6.0,  8.0,  10.0};
49   
50 Float_t AliTRDCalPID::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};
51
52     
53 //_________________________________________________________________________
54 AliTRDCalPID::AliTRDCalPID()
55   :TNamed("pid", "PID for TRD")
56   ,fMeanChargeRatio(0)
57   ,fHistdEdx(0x0)
58   ,fHistTimeBin(0x0)
59 {
60   //
61   //  The Default constructor
62   //
63
64   Init();
65
66 }
67
68 //_________________________________________________________________________
69 AliTRDCalPID::AliTRDCalPID(const Text_t *name, const Text_t *title) 
70   :TNamed(name,title)
71   ,fMeanChargeRatio(0)
72   ,fHistdEdx(0x0)
73   ,fHistTimeBin(0x0)
74 {
75   //
76   //  The main constructor
77   //
78   
79   Init();
80
81 }
82
83 //_____________________________________________________________________________
84 AliTRDCalPID::AliTRDCalPID(const AliTRDCalPID &c) 
85   :TNamed(c)
86   ,fMeanChargeRatio(c.fMeanChargeRatio)
87   ,fHistdEdx(0x0)
88   ,fHistTimeBin(0x0)
89 {
90   //
91   // Copy constructor
92   //
93
94   if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
95   
96 }
97
98 //_________________________________________________________________________
99 AliTRDCalPID::~AliTRDCalPID()
100 {
101   //
102   // Destructor
103   //
104   
105   CleanUp();
106
107 }
108
109 //_________________________________________________________________________
110 void AliTRDCalPID::CleanUp()
111 {
112   //
113   // Delets all newly created objects
114   //
115
116   if (fHistdEdx) {
117     delete fHistdEdx;
118     fHistdEdx = 0x0;
119   }
120   
121   if (fHistTimeBin) {
122     delete fHistTimeBin;
123     fHistTimeBin = 0x0;
124   }
125 }
126
127 //_____________________________________________________________________________
128 AliTRDCalPID &AliTRDCalPID::operator=(const AliTRDCalPID &c)
129 {
130   //
131   // Assignment operator
132   //
133
134   if (this != &c) ((AliTRDCalPID &) c).Copy(*this);
135   return *this;
136
137 }
138
139 //_____________________________________________________________________________
140 void AliTRDCalPID::Copy(TObject &c) const
141 {
142   //
143   // Copy function
144   //
145
146   AliTRDCalPID& target = (AliTRDCalPID &) c;
147   
148   target.CleanUp();
149   
150   target.fMeanChargeRatio = fMeanChargeRatio;
151
152   if (fHistdEdx) {
153     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
154   }
155   if (fHistTimeBin) {
156     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
157   }
158
159   TObject::Copy(c);
160
161 }
162
163 //_________________________________________________________________________
164 void AliTRDCalPID::Init()
165 {
166   //
167   // Initialization
168   //
169
170   fHistdEdx    = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
171   fHistdEdx->SetOwner();
172   fHistTimeBin = new TObjArray(2 * kNMom);
173   fHistTimeBin->SetOwner();  
174
175         // Initialization of estimator at object instantiation because late
176         // initialization in function GetProbability() is not working due to
177         // constantness of this function. 
178         // fEstimator = new AliTRDCalPIDRefMaker();
179         
180   // ADC Gain normalization
181   fMeanChargeRatio = 1.0;
182 }
183
184 //_________________________________________________________________________
185 Bool_t AliTRDCalPID::LoadLQReferences(Char_t *refFile)
186 {
187   //
188   // Read the TRD dEdx histograms.
189   //
190
191         Int_t nTimeBins = 22;
192         // Get number of time bins from CDB
193         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
194         if(!calibration){
195                 AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
196         }else{
197                 if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
198                 else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
199         }
200
201         
202   // Read histogram Root file  
203   TFile *histFile = TFile::Open(refFile, "READ");
204   if (!histFile || !histFile->IsOpen()) {
205     AliError(Form("Opening TRD histgram file %s failed", refFile));
206     return kFALSE;
207   }
208   gROOT->cd();
209
210   // Read histograms
211   for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
212     for (Int_t imom = 0; imom < kNMom; imom++){
213                         TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
214                         hist->Scale(1./hist->Integral());
215                         fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
216
217                         if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
218
219                         TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
220                         if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, nTimeBins));
221                         ht->Scale(1./ht->Integral());
222                         fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
223                 }
224   }
225   
226   histFile->Close();
227   delete histFile;
228   
229   // Number of bins and bin size
230   //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
231   //fNbins   = hist->GetNbinsX();
232   //fBinSize = hist->GetBinWidth(1);
233   
234   return kTRUE;
235
236 }
237
238 // //_________________________________________________________________________
239 // Double_t  AliTRDCalPID::GetMean(Int_t k, Int_t ip) const
240 // {
241 //   //
242 //   // Gets mean of de/dx dist. of e
243 //   //
244 // 
245 //   AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
246 //               ,fpartName[k]
247 //               ,fTrackMomentum[ip]));
248 //   if (k < 0 || k > AliPID::kSPECIES) {
249 //     return 0;
250 //   }
251 // 
252 //   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
253 // 
254 // }
255 // 
256 // //_________________________________________________________________________
257 // Double_t  AliTRDCalPID::GetNormalization(Int_t k, Int_t ip) const
258 // {
259 //   //
260 //   // Gets Normalization of de/dx dist. of e
261 //   //
262 // 
263 //   AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
264 //               ,fpartName[k]
265 //               ,fTrackMomentum[ip]));
266 //   if (k < 0 || k > AliPID::kSPECIES) {
267 //     return 0;
268 //   }
269 //   
270 //   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
271 // 
272 // }
273
274 //_________________________________________________________________________
275 TH1* AliTRDCalPID::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
276 {
277   //
278   // Returns one selected dEdx histogram
279   //
280
281   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
282         if(ip<0 || ip>= kNMom ) return 0x0;
283
284         AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
285   
286   return (TH1*)fHistdEdx->At(GetHistID(k, ip));
287
288 }
289
290 //_________________________________________________________________________
291 TH1* AliTRDCalPID::GetHistogramT(Int_t k, Int_t ip) const
292 {
293   //
294   // Returns one selected time bin max histogram
295   //
296
297   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
298         if(ip<0 || ip>= kNMom ) return 0x0;
299           
300         AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
301
302         return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*kNMom+ip);
303 }
304
305
306
307 //_________________________________________________________________________
308 Double_t AliTRDCalPID::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const
309 {
310   //
311         // Core function of AliTRDCalPID class for calculating the
312         // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
313         // given momentum "mom" and a given dE/dx (slice "dedx") yield per
314         // layer
315   //
316
317         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
318                 
319         //Double_t dedx   = dedx1/fMeanChargeRatio;
320         
321         // find the interval in momentum and track segment length which applies for this data
322         Int_t ilength = 1;
323   while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
324                 ilength++;
325         }
326         Int_t imom = 1;
327   while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
328         
329         Int_t nbinsx, nbinsy;
330         TAxis *ax = 0x0, *ay = 0x0;
331         Double_t LQ1, LQ2;
332         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
333         TH2 *hist = 0x0;
334         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom-1/*, ilength*/)))){
335                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
336                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom-1)));
337                 return 0.;
338         }
339         ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
340         ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
341         Float_t x = dedx[0]+dedx[1], y = dedx[2];
342   Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
343         Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
344         if(kX)
345                 if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y)); 
346     //fEstimator->Estimate2D2(hist, x, y);
347                 else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
348         else
349                 if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
350                 else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
351
352
353         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
354                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
355                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
356                 return LQ1;
357         }
358         if(kX)
359                 if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y)); 
360     //fEstimator->Estimate2D2(hist, x, y);
361                 else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
362         else
363                 if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
364                 else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
365
366         
367         // return interpolation over momentum binning
368   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
369
370 }
371
372 //_________________________________________________________________________
373 Double_t AliTRDCalPID::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
374 {
375   //
376   // Gets the Probability of having timbin at a given momentum (mom)
377   // and particle type (spec) (0 for e) and (2 for pi)
378   // from the precalculated timbin distributions 
379   //
380   
381         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
382
383   Int_t iTBin = timbin+1;
384   
385   // Everything which is not an electron counts as a pion for time bin max
386   //if(spec != AliPID::kElectron) spec = AliPID::kPion;
387   
388
389   
390         Int_t imom = 1;
391   while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
392
393         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
394         TH1F *hist = 0x0;
395         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom-1))){
396                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
397                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom-1));
398                 return 0.;
399         }
400         Double_t LQ1 = hist->GetBinContent(iTBin);
401
402         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom))){
403                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
404                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom));
405                 return LQ1;
406         }
407         Double_t LQ2 = hist->GetBinContent(iTBin);
408
409         // return interpolation over momentum binning
410   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
411 }
412