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