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