Coverity 16571
[u/mrichter/AliRoot.git] / STEER / AliTRDPIDResponse.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 // PID Response class for the TRD detector
17 // Based on 1D Likelihood approach
18 // Calculation of probabilities using Bayesian approach
19 // Attention: This method is only used to separate electrons from pions
20 //
21 // Authors:
22 //  Markus Fasel <M.Fasel@gsi.de>
23 //  Anton Andronic <A.Andronic@gsi.de>
24 //
25 #include <TClass.h>
26 #include <TAxis.h>
27 #include <TFile.h>
28 #include <TH1.h>
29 #include <TKey.h>
30 #include <TObjArray.h>
31 #include <TString.h>
32 #include <TROOT.h> 
33 #include <TSystem.h>
34 #include <TDirectory.h>
35
36 #include "AliLog.h"
37
38 #include "AliTRDPIDResponse.h"
39
40 ClassImp(AliTRDPIDResponse)
41
42 const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
43
44 //____________________________________________________________
45 AliTRDPIDResponse::AliTRDPIDResponse():
46   TObject()
47   ,fReferences(NULL)
48   ,fGainNormalisationFactor(1.)
49   ,fPIDmethod(kLQ1D)
50 {
51   //
52   // Default constructor
53   //
54   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
55     for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
56       fMapRefHists[ispec][ipbin] = -1;
57 }
58
59 //____________________________________________________________
60 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
61   TObject(ref)
62   ,fReferences(ref.fReferences)
63   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
64   ,fPIDmethod(ref.fPIDmethod)
65 {
66   //
67   // Copy constructor
68   // Flat copy of the reference histos
69   // For creating a deep copy call SetOwner
70   //
71   Int_t size  = (AliPID::kSPECIES)*(kNPBins);
72   memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
73 }
74
75 //____________________________________________________________
76 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
77   //
78   // Assignment operator
79   // Performs a flat copy of the reference histos
80   //
81   if(this == &ref) return *this;
82   
83   // Make copy
84   TObject::operator=(ref);
85   if(IsOwner()){
86     if(fReferences){
87       fReferences->Delete();
88       delete fReferences;
89     }
90     SetBit(kIsOwner, kFALSE);
91   } else if(fReferences) delete fReferences;
92   fReferences = ref.fReferences;
93   Int_t size  = (AliPID::kSPECIES)*(kNPBins);
94   memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
95   fPIDmethod = ref.fPIDmethod;
96   
97   return *this;
98 }
99
100 //____________________________________________________________
101 AliTRDPIDResponse::~AliTRDPIDResponse(){
102   //
103   // Destructor
104   //
105   if(IsOwner()){
106     // Destroy histos
107     if(fReferences) fReferences->Delete();
108   }
109   if(fReferences) delete fReferences;
110 }
111
112 //____________________________________________________________
113 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
114   //
115   // Load References into the toolkit
116   //
117   AliDebug(1, "Loading reference histos from root file");
118   TDirectory *owd = gDirectory;// store old working directory
119   
120   if(!filename)
121     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
122   TFile *in = TFile::Open(filename);
123   if(!in){
124     AliError("Ref file not available.");
125     return kFALSE;
126   }
127   
128   gROOT->cd();
129   fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
130   fReferences->SetOwner();
131   TIter iter(in->GetListOfKeys());
132   TKey *histkey = NULL;
133   TObject *tmp = NULL;
134   Int_t arrayPos = 0, pbin = 0; 
135   AliPID::EParticleType species;
136   TString histname;
137   TH1 *hnew = NULL;
138   while((histkey = dynamic_cast<TKey *>(iter()))){
139     tmp = histkey->ReadObj();
140     histname = tmp->GetName();
141     if(histname.BeginsWith("fHQel")){ // Electron histogram
142       histname.ReplaceAll("fHQel_p","");
143       species = AliPID::kElectron;
144     } else if(histname.BeginsWith("fHQmu")){ // Muon histogram
145       histname.ReplaceAll("fHQmu_p","");
146       species = AliPID::kMuon;
147     } else if(histname.BeginsWith("fHQpi")){ // Pion histogram
148       histname.ReplaceAll("fHQpi_p","");
149       species = AliPID::kPion;
150     }else if(histname.BeginsWith("fHQka")){ // Kaon histogram
151       histname.ReplaceAll("fHQka_p","");
152       species = AliPID::kKaon;
153     }else if(histname.BeginsWith("fHQpr")){ // Proton histogram
154       histname.ReplaceAll("fHQpr_p","");
155       species = AliPID::kProton;
156     } else continue;
157     pbin = histname.Atoi() - 1;
158     AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
159     fMapRefHists[species][pbin] = arrayPos;
160     fReferences->AddAt((hnew =new TH1F(*dynamic_cast<TH1F *>(tmp))), arrayPos);
161     hnew->SetDirectory(gROOT);
162     arrayPos++;
163   } 
164   in->Close();
165   owd->cd();
166   AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
167   SetBit(kIsOwner, kTRUE);
168   delete in;
169   return kTRUE;
170 }
171
172 //____________________________________________________________
173 Bool_t AliTRDPIDResponse::Load(const TObjArray *histos){
174   //
175   // Load Reference into the PID Response (from OCDB/OADB)
176   //
177   AliDebug(1, "Setting references (from the database)");
178   if(fReferences) fReferences->Clear();
179   else
180     fReferences = new TObjArray(AliPID::kSPECIES*kNPBins);
181   fReferences->SetOwner();
182   TIter iter(histos);
183   Int_t arrayPos = 0, pbin = 0; 
184   AliPID::EParticleType species;
185   TObject *tmp = NULL;
186   TString histname;
187   TH1 *hnew = NULL;
188   while((tmp = iter())){
189     histname = tmp->GetName();
190     if(histname.BeginsWith("fHQel")){ // Electron histogram
191       histname.ReplaceAll("fHQel_p","");
192       species = AliPID::kElectron;
193     } else if(histname.BeginsWith("fHQmu")){ // Muon histogram
194       histname.ReplaceAll("fHQmu_p","");
195       species = AliPID::kMuon;
196     } else if(histname.BeginsWith("fHQpi")){ // Pion histogram
197       histname.ReplaceAll("fHQpi_p","");
198       species = AliPID::kPion;
199     }else if(histname.BeginsWith("fHQka")){ // Kaon histogram
200       histname.ReplaceAll("fHQka_p","");
201       species = AliPID::kKaon;
202     }else if(histname.BeginsWith("fHQpr")){ // Proton histogram
203       histname.ReplaceAll("fHQpr_p","");
204       species = AliPID::kProton;
205     } else continue;
206     pbin = histname.Atoi() - 1;
207     AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
208     fMapRefHists[species][pbin] = arrayPos;
209     fReferences->AddAt((hnew = new TH1F(*dynamic_cast<TH1F *>(tmp))), arrayPos);
210     hnew->SetDirectory(gROOT);
211     arrayPos++;
212   }
213   if(!arrayPos){
214     AliError("Failed loading References");
215     return kFALSE;
216   }
217   AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
218   SetBit(kIsOwner, kTRUE);
219   return kTRUE;
220 }
221
222 //____________________________________________________________
223 Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
224 {
225   //
226   // Calculate TRD likelihood values for the track based on dedx and 
227   // momentum values. The likelihoods are calculated by query the 
228   // reference data depending on the PID method selected
229   //
230   // Input parameter :
231   //   n - number of dedx slices/chamber
232   //   dedx - array of dedx slices organized layer wise
233   //   p - array of momentum measurements organized layer wise
234   // 
235   // Return parameters
236   //   prob - probabilities allocated by TRD for particle specis
237   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
238   // 
239   // Return value
240   //   true if calculation success
241   // 
242
243   if(!fReferences){
244     AliWarning("Missing reference data. PID calculation not possible.");
245     return kFALSE;
246   }
247
248   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
249   Double_t prLayer[AliPID::kSPECIES];
250   Double_t DE[10], s(0.);
251   for(Int_t il(kNlayer); il--;){
252     memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
253     if(!CookdEdx(n, &dedx[il*n], &DE[0])) continue;
254
255     s=0.;
256     for(Int_t is(AliPID::kSPECIES); is--;){
257       if((DE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
258       AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
259       s+=prLayer[is];
260     }
261     if(s<1.e-30){
262       AliDebug(2, Form("Null all species prob layer[%d].", il));
263       continue;
264     }
265     for(Int_t is(AliPID::kSPECIES); is--;){
266       if(kNorm) prLayer[is] /= s;
267       prob[is] *= prLayer[is];
268     }
269   }
270   if(!kNorm) return kTRUE;
271
272   s=0.;
273   for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
274   if(s<1.e-30){
275     AliDebug(2, "Null total prob.");
276     return kFALSE;
277   }
278   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
279   return kTRUE;
280 }
281
282 //____________________________________________________________
283 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
284   //
285   // Get the non-normalized probability for a certain particle species as coming
286   // from the reference histogram
287   // Interpolation between momentum bins
288   //
289   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
290   Int_t pbin = GetLowerMomentumBin(plocal);  
291   AliDebug(1, Form("Bin %d", pbin));
292   Double_t pLayer = 0.;
293   // Do Interpolation exept for underflow and overflow
294   if(pbin >= 0 && pbin < kNPBins-1){
295     TH1 *refHistUpper = NULL, *refHistLower = NULL;
296     if(fMapRefHists[species][pbin] >= 0)
297       refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
298     if(fMapRefHists[species][pbin+1] >= 0)
299       refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
300     AliDebug(1, Form("Reference Histos (Upper/Lower): [%p|%p]", refHistUpper, refHistLower));
301   
302     if (refHistLower && refHistUpper ) {
303       Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
304       Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
305   
306       pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]); 
307     }
308   }
309   else{
310     TH1 *refHist = NULL;
311     if(pbin < 0){
312       // underflow
313       if(fMapRefHists[species][0] >= 0){
314         refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0])); 
315       }
316     } else {
317       // overflow
318       if(fMapRefHists[species][kNPBins-1] > 0)
319         refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1])); 
320     }
321     if (refHist)
322       pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
323   }
324   AliDebug(1, Form("Probability %f", pLayer));
325   return pLayer;
326 }
327
328 //____________________________________________________________
329 Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
330   //
331   // Get the momentum bin for a given momentum value
332   //
333   Int_t bin = -1;
334   if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
335   for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
336     if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
337       bin = ibin;
338       break;
339     }
340   }
341   return bin;
342 }
343
344 //____________________________________________________________
345 void AliTRDPIDResponse::SetOwner(){
346   //
347   // Make Deep Copy of the Reference Histograms
348   //
349   if(!fReferences || IsOwner()) return;
350   TObjArray *ctmp = new TObjArray();
351   for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
352     ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
353   delete fReferences;
354   fReferences = ctmp;
355   SetBit(kIsOwner, kTRUE);
356 }
357
358 //____________________________________________________________
359 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, Double_t *in, Double_t *out) const
360 {
361   switch(fPIDmethod){
362   case kNN: // NN 
363     break;
364   case kLQ2D: // 2D LQ 
365     break;
366   case kLQ1D: // 1D LQ 
367     out[0]= 0.;
368     for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
369     if(out[0] < 1e-6) return kFALSE;
370     break;
371   default:
372     return kFALSE;
373   }
374   return kTRUE;
375 }
376