Updating AOD filter
[u/mrichter/AliRoot.git] / STEER / AliTRDPIDResponse.cxx
CommitLineData
ffb1ee30 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>
ce5d6908 29#include <TKey.h>
ffb1ee30 30#include <TObjArray.h>
31#include <TString.h>
ce5d6908 32#include <TROOT.h>
ffb1ee30 33#include <TSystem.h>
2ad33025 34#include <TDirectory.h>
ffb1ee30 35
36#include "AliLog.h"
37
38#include "AliTRDPIDResponse.h"
39
40ClassImp(AliTRDPIDResponse)
41
42const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
43
44//____________________________________________________________
e0de37e9 45AliTRDPIDResponse::AliTRDPIDResponse():
46 TObject()
47 ,fReferences(NULL)
48 ,fGainNormalisationFactor(1.)
49 ,fPIDmethod(kLQ1D)
ffb1ee30 50{
e0de37e9 51 //
52 // Default constructor
53 //
ce5d6908 54 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
55 for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
56 fMapRefHists[ispec][ipbin] = -1;
ffb1ee30 57}
58
59//____________________________________________________________
60AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9 61 TObject(ref)
62 ,fReferences(ref.fReferences)
63 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
64 ,fPIDmethod(ref.fPIDmethod)
ffb1ee30 65{
e0de37e9 66 //
67 // Copy constructor
68 // Flat copy of the reference histos
69 // For creating a deep copy call SetOwner
70 //
ce5d6908 71 Int_t size = (AliPID::kSPECIES)*(kNPBins);
43b986eb 72 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Int_t) * size);
ffb1ee30 73}
74
75//____________________________________________________________
76AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9 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);
43b986eb 94 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Int_t) * size);
e0de37e9 95 fPIDmethod = ref.fPIDmethod;
96
97 return *this;
ffb1ee30 98}
99
100//____________________________________________________________
101AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9 102 //
103 // Destructor
104 //
105 if(IsOwner()){
106 // Destroy histos
107 if(fReferences) fReferences->Delete();
108 }
109 if(fReferences) delete fReferences;
ffb1ee30 110}
111
112//____________________________________________________________
113Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
e0de37e9 114 //
115 // Load References into the toolkit
116 //
117 AliDebug(1, "Loading reference histos from root file");
2ad33025 118 TDirectory *owd = gDirectory;// store old working directory
e0de37e9 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}
2ad33025 171
e0de37e9 172//____________________________________________________________
173Bool_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;
ffb1ee30 220}
221
222//____________________________________________________________
b439f460 223Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
ffb1ee30 224{
e0de37e9 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 //
ffb1ee30 242
e0de37e9 243 if(!fReferences){
244 AliWarning("Missing reference data. PID calculation not possible.");
245 return kFALSE;
246 }
ffb1ee30 247
e0de37e9 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;
ffb1ee30 254
e0de37e9 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;
ffb1ee30 271
e0de37e9 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;
ffb1ee30 280}
281
ffb1ee30 282//____________________________________________________________
b439f460 283Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
e0de37e9 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;
ffb1ee30 326}
327
328//____________________________________________________________
b439f460 329Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
e0de37e9 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;
ffb1ee30 342}
343
344//____________________________________________________________
345void AliTRDPIDResponse::SetOwner(){
e0de37e9 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);
ffb1ee30 356}
357
358//____________________________________________________________
e0de37e9 359Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, Double_t *in, Double_t *out) const
ffb1ee30 360{
e0de37e9 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;
ffb1ee30 375}
376