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