Corrected parameterization of the error. Overlooked in the previous commit.
[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>
29#include <TObjArray.h>
30#include <TString.h>
31#include <TSystem.h>
32
33#include "AliLog.h"
34
35#include "AliTRDPIDResponse.h"
36
37ClassImp(AliTRDPIDResponse)
38
39const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
40
41//____________________________________________________________
42AliTRDPIDResponse::AliTRDPIDResponse(const Char_t * filename):
43 TObject()
44 ,fReferences(NULL)
45 ,fPIDmethod(2)
46{
47 //
48 // Default constructor
49 //
338c9979 50 Int_t size = (AliPID::kSPECIES+1)*(kNPBins+1);
51 memset(fMapRefHists, 0, sizeof(Int_t) * size);
ffb1ee30 52 Load(filename);
53}
54
55//____________________________________________________________
56AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
57 TObject(ref)
58 ,fReferences(ref.fReferences)
59 ,fPIDmethod(ref.fPIDmethod)
60{
61 //
62 // Copy constructor
63 // Flat copy of the reference histos
64 // For creating a deep copy call SetOwner
65 //
338c9979 66 Int_t size = (AliPID::kSPECIES+1)*(kNPBins+1);
67 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
ffb1ee30 68}
69
70//____________________________________________________________
71AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
72 //
73 // Assignment operator
74 // Performs a flat copy of the reference histos
75 //
76 if(this == &ref) return *this;
77
78 // Make copy
79 TObject::operator=(ref);
80 if(IsOwner()){
81 if(fReferences){
82 fReferences->Delete();
83 delete fReferences;
84 }
85 SetBit(kIsOwner, kFALSE);
86 }
87 fReferences = ref.fReferences;
338c9979 88 Int_t size = (AliPID::kSPECIES+1)*(kNPBins+1);
89 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
ffb1ee30 90 fPIDmethod = ref.fPIDmethod;
91
92 return *this;
93}
94
95//____________________________________________________________
96AliTRDPIDResponse::~AliTRDPIDResponse(){
97 //
98 // Destructor
99 //
100 if(fReferences && IsOwner()){
101 fReferences->Delete();
102 delete fReferences;
103 }
104}
105
106//____________________________________________________________
107Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
108 //
109 // Load References into the toolkit
110 //
111 if(!filename)
112 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
113 TFile *in = TFile::Open(filename);
114 if(!in){
115 AliError("Ref file not available.");
116 return kFALSE;
117 }
118
119 fReferences = new TObjArray();
120 fReferences->SetOwner();
121 TIter iter(in->GetListOfKeys());
122 TObject *tmp = NULL;
123 Int_t arrayPos = 0, pbin = 0;
124 AliPID::EParticleType species;
125 TString histname;
126 while((tmp = iter())){
127 histname = tmp->GetName();
128 if(histname.BeginsWith("fHQel")){ // Electron histogram
129 histname.ReplaceAll("fHQel_p","");
130 species = AliPID::kElectron;
131 } else { // Pion histogram
132 histname.ReplaceAll("fHQpi_p","");
133 species = AliPID::kPion;
134 }
135 pbin = histname.Atoi();
136 fMapRefHists[species][pbin] = arrayPos;
137 fReferences->AddAt(tmp->Clone(), arrayPos);
138 arrayPos++;
139 }
140 in->Close();
141
142 AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
143 delete in;
144 return kTRUE;
145}
146
147//____________________________________________________________
b439f460 148Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
ffb1ee30 149{
150//
151// Calculate TRD likelihood values for the track based on dedx and
152// momentum values. The likelihoods are calculated by query the
153// reference data depending on the PID method selected
154//
155// Input parameter :
156// n - number of dedx slices/chamber
157// dedx - array of dedx slices organized layer wise
158// p - array of momentum measurements organized layer wise
159//
160// Return parameters
161// prob - probabilities allocated by TRD for particle specis
162// kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
163//
164// Return value
165// true if calculation success
166//
167
168 if(!fReferences){
169 AliWarning("Missing reference data. PID calculation not possible.");
170 return kFALSE;
171 }
172
173 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
174 Double_t prLayer[AliPID::kSPECIES];
175 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
176 Double_t DE[10], s(0.);
177 for(Int_t il(kNlayer); il--;){
178 if(!CookdEdx(&dedx[il*n], &DE[0])) continue;
179
180 s=0.;
181 for(Int_t is(AliPID::kSPECIES); is--;){
182 prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
183 s+=prLayer[is];
184 }
185 if(s<1.e-30){
186 AliDebug(2, Form("Null all species prob layer[%d].", il));
187 continue;
188 }
189 for(Int_t is(AliPID::kSPECIES); is--;){
190 if(kNorm) prLayer[is] /= s;
191 prob[is] *= prLayer[is];
192 }
193 }
194 if(!kNorm) return kTRUE;
195
196 s=0.;
197 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
198 if(s<1.e-30){
199 AliDebug(2, "Null total prob.");
200 return kFALSE;
201 }
202 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
203 return kTRUE;
204}
205
206
207//____________________________________________________________
b439f460 208Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
ffb1ee30 209 //
210 // Get the non-normalized probability for a certain particle species as coming
211 // from the reference histogram
212 // Interpolation between momentum bins
213 //
214 Int_t pbin = GetLowerMomentumBin(plocal);
215 Double_t pLayer = 0.;
216 // Do Interpolation exept for underflow and overflow
217 if(pbin >= 0 && pbin < kNPBins-1){
218 TH1 *refHistUpper = NULL, *refHistLower = NULL;
219 refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
220 refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
221
338c9979 222 if (refHistLower && refHistUpper ) {
223 Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
224 Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
ffb1ee30 225
338c9979 226 pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * plocal;
227 }
ffb1ee30 228 }
229 else{
230 TH1 *refHist = NULL;
231 if(pbin < 0){
232 // underflow
233 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[0][0]));
234 } else {
235 // overflow
236 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins]));
237 }
338c9979 238 if (refHist)
239 pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
ffb1ee30 240 }
241 return pLayer;
242}
243
244//____________________________________________________________
b439f460 245Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
ffb1ee30 246 //
247 // Get the momentum bin for a given momentum value
248 //
249 Int_t bin = -1;
250 if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
251 for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
252 if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
253 bin = ibin;
254 break;
255 }
256 }
257 return bin;
258}
259
260//____________________________________________________________
261void AliTRDPIDResponse::SetOwner(){
262 //
263 // Make Deep Copy of the Reference Histograms
264 //
265 if(!fReferences || IsOwner()) return;
266 TObjArray *ctmp = new TObjArray();
267 for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
268 ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
269 fReferences = ctmp;
270 SetBit(kIsOwner, kTRUE);
271}
272
273//____________________________________________________________
b439f460 274Bool_t AliTRDPIDResponse::CookdEdx(Double_t *in, Double_t *out) const
ffb1ee30 275{
276 switch(fPIDmethod){
277 case 0: // NN
278 break;
279 case 1: // 2D LQ
280 break;
281 case 2: // 1D LQ
282 out[0]= 0.;
283 for(Int_t islice = 0; islice < 8; islice++) out[0] += in[islice];
284 if(out[0] < 1e-6) return kFALSE;
285 break;
286 default:
287 return kFALSE;
288 }
289 return kTRUE;
290}
291