Added TPC crossed rows customizable cuts and values; added TOF deltas values
[u/mrichter/AliRoot.git] / STEER / STEERBase / 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//
ffb1ee30 25#include <TAxis.h>
ea235c90 26#include <TClass.h>
27#include <TDirectory.h>
ffb1ee30 28#include <TFile.h>
29#include <TH1.h>
ce5d6908 30#include <TKey.h>
ea235c90 31#include <TMath.h>
ffb1ee30 32#include <TObjArray.h>
ce5d6908 33#include <TROOT.h>
ea235c90 34#include <TString.h>
ffb1ee30 35#include <TSystem.h>
36
37#include "AliLog.h"
38
db0e2c5f 39#include "AliTRDPIDResponseObject.h"
ffb1ee30 40#include "AliTRDPIDResponse.h"
db0e2c5f 41#include "AliTRDTKDInterpolator.h"
ffb1ee30 42
43ClassImp(AliTRDPIDResponse)
44
ffb1ee30 45//____________________________________________________________
e0de37e9 46AliTRDPIDResponse::AliTRDPIDResponse():
47 TObject()
db0e2c5f 48 ,fkPIDResponseObject(NULL)
e0de37e9 49 ,fGainNormalisationFactor(1.)
ffb1ee30 50{
e0de37e9 51 //
52 // Default constructor
53 //
ffb1ee30 54}
55
56//____________________________________________________________
57AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9 58 TObject(ref)
db0e2c5f 59 ,fkPIDResponseObject(NULL)
e0de37e9 60 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
ffb1ee30 61{
e0de37e9 62 //
63 // Copy constructor
e0de37e9 64 //
ffb1ee30 65}
66
67//____________________________________________________________
68AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9 69 //
70 // Assignment operator
e0de37e9 71 //
72 if(this == &ref) return *this;
73
74 // Make copy
75 TObject::operator=(ref);
51a0ce25 76 fGainNormalisationFactor = ref.fGainNormalisationFactor;
db0e2c5f 77 fkPIDResponseObject = ref.fkPIDResponseObject;
e0de37e9 78
79 return *this;
ffb1ee30 80}
81
82//____________________________________________________________
83AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9 84 //
85 // Destructor
86 //
db0e2c5f 87 if(IsOwner()) delete fkPIDResponseObject;
ffb1ee30 88}
89
90//____________________________________________________________
db0e2c5f 91Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
e0de37e9 92 //
93 // Load References into the toolkit
94 //
95 AliDebug(1, "Loading reference histos from root file");
2ad33025 96 TDirectory *owd = gDirectory;// store old working directory
e0de37e9 97
98 if(!filename)
99 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
100 TFile *in = TFile::Open(filename);
101 if(!in){
102 AliError("Ref file not available.");
103 return kFALSE;
104 }
105
106 gROOT->cd();
db0e2c5f 107 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
51a0ce25 108 in->Close(); delete in;
e0de37e9 109 owd->cd();
e0de37e9 110 SetBit(kIsOwner, kTRUE);
db0e2c5f 111 AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
e0de37e9 112 return kTRUE;
ffb1ee30 113}
114
115//____________________________________________________________
bd58d4b9 116Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const
ffb1ee30 117{
e0de37e9 118 //
119 // Calculate TRD likelihood values for the track based on dedx and
120 // momentum values. The likelihoods are calculated by query the
121 // reference data depending on the PID method selected
122 //
123 // Input parameter :
124 // n - number of dedx slices/chamber
125 // dedx - array of dedx slices organized layer wise
126 // p - array of momentum measurements organized layer wise
127 //
128 // Return parameters
129 // prob - probabilities allocated by TRD for particle specis
130 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
131 //
132 // Return value
133 // true if calculation success
bd58d4b9 134 //
135 AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
136
ffb1ee30 137
db0e2c5f 138 if(!fkPIDResponseObject){
139 AliWarning("Missing reference data. PID calculation not possible.");
140 return kFALSE;
141 }
ffb1ee30 142
db0e2c5f 143 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
144 Double_t prLayer[AliPID::kSPECIES];
145 Double_t dE[10], s(0.);
146 for(Int_t il(kNlayer); il--;){
147 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
bd58d4b9 148 if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
ffb1ee30 149
db0e2c5f 150 s=0.;
151 Bool_t filled=kTRUE;
152 for(Int_t is(AliPID::kSPECIES); is--;){
bd58d4b9 153 if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
154 if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
f2762b1c 155 AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
db0e2c5f 156 if(prLayer[is]<1.e-30){
157 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
158 filled=kFALSE;
159 break;
160 }
161 s+=prLayer[is];
162 }
163 if(!filled){
164 continue;
165 }
166 if(s<1.e-30){
167 AliDebug(2, Form("Null all species prob layer[%d].", il));
168 continue;
169 }
170 for(Int_t is(AliPID::kSPECIES); is--;){
171 if(kNorm) prLayer[is] /= s;
172 prob[is] *= prLayer[is];
173 }
e0de37e9 174 }
db0e2c5f 175 if(!kNorm) return kTRUE;
176
177 s=0.;
178 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
e0de37e9 179 if(s<1.e-30){
db0e2c5f 180 AliDebug(2, "Null total prob.");
181 return kFALSE;
e0de37e9 182 }
db0e2c5f 183 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
184 return kTRUE;
ffb1ee30 185}
186
ffb1ee30 187//____________________________________________________________
bd58d4b9 188Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
e0de37e9 189 //
190 // Get the non-normalized probability for a certain particle species as coming
191 // from the reference histogram
192 // Interpolation between momentum bins
193 //
194 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
db0e2c5f 195
51a0ce25 196 Double_t probLayer = 0.;
db0e2c5f 197 Float_t pLower, pUpper;
198
bd58d4b9 199 switch(PIDmethod){
db0e2c5f 200 case kNN: // NN
201 break;
202 case kLQ2D: // 2D LQ
203 {
f2762b1c 204 if(species==0||species==2){ // references only for electrons and pions
205 Double_t error;
206 Double_t point[kNslicesLQ2D];
207 for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
db0e2c5f 208
f2762b1c 209 AliTRDTKDInterpolator *refLower = dynamic_cast<AliTRDTKDInterpolator*>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
db0e2c5f 210
f2762b1c 211 if(refLower){
212 refLower->Eval(point,probLayer,error);
213 }
214 else {
215 AliError("No references available");
216 }
217 AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
db0e2c5f 218 }
219 }
220 break;
221 case kLQ1D: // 1D LQ
222 {
223 TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
224 *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
225 // Do Interpolation exept for underflow and overflow
226 if(refLower && refUpper){
227 Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
228 Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
229
230 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
231 } else if(refLower){
232 // underflow
233 probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
234 } else if(refUpper){
235 // overflow
236 probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
237 } else {
238 AliError("No references available");
239 }
f2762b1c 240 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
db0e2c5f 241 }
242 break;
243 default:
244 break;
e0de37e9 245 }
51a0ce25 246 return probLayer;
ffb1ee30 247}
248
249//____________________________________________________________
250void AliTRDPIDResponse::SetOwner(){
e0de37e9 251 //
252 // Make Deep Copy of the Reference Histograms
253 //
db0e2c5f 254 if(!fkPIDResponseObject || IsOwner()) return;
255 const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
256 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
257 SetBit(kIsOwner, kTRUE);
ffb1ee30 258}
259
260//____________________________________________________________
bd58d4b9 261Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
ffb1ee30 262{
f2762b1c 263 //
264 // Recalculate dE/dx
265 //
bd58d4b9 266 switch(PIDmethod){
e0de37e9 267 case kNN: // NN
db0e2c5f 268 break;
269 case kLQ2D: // 2D LQ
270 out[0]=0;
271 out[1]=0;
272 for(Int_t islice = 0; islice < nSlice; islice++){
273 if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
274 else out[1]+= in[islice];
275 }
f2762b1c 276 if(out[0] < 1e-6) return kFALSE;
277 AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
db0e2c5f 278 break;
279 case kLQ1D: // 1D LQ
280 out[0]= 0.;
281 for(Int_t islice = 0; islice < nSlice; islice++)
282 if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
283 if(out[0] < 1e-6) return kFALSE;
f2762b1c 284 AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
db0e2c5f 285 break;
286
e0de37e9 287 default:
db0e2c5f 288 return kFALSE;
e0de37e9 289 }
290 return kTRUE;
ffb1ee30 291}
292
ea235c90 293//____________________________________________________________
bd58d4b9 294Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
db0e2c5f 295 //
296 // Check whether particle is identified as electron assuming a certain electron efficiency level
297 // Only electron and pion hypothesis is taken into account
298 //
299 // Inputs:
300 // Number of tracklets
301 // Likelihood values
302 // Momentum
303 // Electron efficiency level
304 //
305 // If the function fails when the params are not accessible, the function returns true
306 //
307 if(!fkPIDResponseObject){
ea235c90 308 AliError("No PID Param object available");
309 return kTRUE;
310 }
311 Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
ce487a7f 312 Double_t params[4];
bd58d4b9 313 if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
ea235c90 314 AliError("No Params found for the given configuration");
315 return kTRUE;
316 }
ce487a7f 317 Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
ea235c90 318 if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
319 return kFALSE;
320}
321
f2762b1c 322//____________________________________________________________
323Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
324
325 fkPIDResponseObject = obj;
326 if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
327 return kTRUE;
328}