TRD dEdx OADB framework: Nch and Ncls dependence added in the OADB object
[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"
5cd0300d 38
39#include "AliVTrack.h"
ffb1ee30 40
db0e2c5f 41#include "AliTRDPIDResponseObject.h"
ffb1ee30 42#include "AliTRDPIDResponse.h"
239fe91c 43//#include "AliTRDTKDInterpolator.h"
44#include "AliTRDNDFast.h"
9c499471 45#include "AliTRDdEdxParams.h"
ffb1ee30 46
47ClassImp(AliTRDPIDResponse)
48
ffb1ee30 49//____________________________________________________________
e0de37e9 50AliTRDPIDResponse::AliTRDPIDResponse():
51 TObject()
db0e2c5f 52 ,fkPIDResponseObject(NULL)
9c499471 53 ,fkTRDdEdxParams(NULL)
e0de37e9 54 ,fGainNormalisationFactor(1.)
ffb1ee30 55{
e0de37e9 56 //
57 // Default constructor
58 //
ffb1ee30 59}
60
61//____________________________________________________________
62AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9 63 TObject(ref)
db0e2c5f 64 ,fkPIDResponseObject(NULL)
9c499471 65 ,fkTRDdEdxParams(NULL)
e0de37e9 66 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
ffb1ee30 67{
e0de37e9 68 //
69 // Copy constructor
e0de37e9 70 //
ffb1ee30 71}
72
73//____________________________________________________________
74AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9 75 //
76 // Assignment operator
e0de37e9 77 //
78 if(this == &ref) return *this;
79
80 // Make copy
81 TObject::operator=(ref);
51a0ce25 82 fGainNormalisationFactor = ref.fGainNormalisationFactor;
db0e2c5f 83 fkPIDResponseObject = ref.fkPIDResponseObject;
9c499471 84 fkTRDdEdxParams = ref.fkTRDdEdxParams;
85
e0de37e9 86 return *this;
ffb1ee30 87}
88
89//____________________________________________________________
90AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9 91 //
92 // Destructor
93 //
9c499471 94 if(IsOwner()) {
95 delete fkPIDResponseObject;
96 delete fkTRDdEdxParams;
97 }
ffb1ee30 98}
99
100//____________________________________________________________
db0e2c5f 101Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
e0de37e9 102 //
103 // Load References into the toolkit
104 //
105 AliDebug(1, "Loading reference histos from root file");
2ad33025 106 TDirectory *owd = gDirectory;// store old working directory
e0de37e9 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 gROOT->cd();
db0e2c5f 117 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
51a0ce25 118 in->Close(); delete in;
e0de37e9 119 owd->cd();
e0de37e9 120 SetBit(kIsOwner, kTRUE);
db0e2c5f 121 AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
e0de37e9 122 return kTRUE;
ffb1ee30 123}
124
5cd0300d 125
126
127//____________________________________________________________
128Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
129{
130 //
131 //calculate the TRD nSigma
132 //
133
134 const Double_t badval = -9999;
135 Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
136
137 const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
138
139 const Double_t mean = info[0];
140 const Double_t res = info[1];
141 if(res<0){
142 return badval;
143 }
144
145 const Double_t sigma = mean*res;
bebe75fc 146 const Double_t eps = 1e-12;
147 return delta/(sigma + eps);
5cd0300d 148}
149
150//____________________________________________________________
151Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
152{
153 //
154 //calculate the TRD signal difference w.r.t. the expected
155 //output other information in info[]
156 //
157
9c499471 158 const Double_t badval = -9999;
5cd0300d 159 Double_t pTRD = -999;
160 for(Int_t ich=0; ich<6; ich++){
161 pTRD = track->GetTRDmomentum(ich);
162 if(pTRD>0)
163 break;
164 }
165 if(pTRD<0){
166 return badval;
167 }
168
9c499471 169 if(!fkTRDdEdxParams){
2ddf4abb 170 AliError("fkTRDdEdxParams null");
9c499471 171 return -99999;
172 }
173
a8da1749 174 const Double_t nch = track->GetTRDNchamberdEdx();
175 const Double_t ncls = track->GetTRDNclusterdEdx();
176
177 const TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type, nch, ncls);
178 if(meanvec.GetNrows()==0){
179 return badval;
180 }
181
182 const TVectorF resvec = fkTRDdEdxParams->GetSigmaParameter(type, nch, ncls);
183 if(resvec.GetNrows()==0){
184 return badval;
185 }
186
187 const Float_t *meanpar = meanvec.GetMatrixArray();
188 const Float_t *respar = resvec.GetMatrixArray();
9c499471 189
190 //============================================================================================<<<<<<<<<<<<<
191
5cd0300d 192 const Double_t bg = pTRD/AliPID::ParticleMass(type);
9c499471 193 const Double_t expsig = MeandEdxTR(&bg, meanpar);
5cd0300d 194
195 if(info){
196 info[0]= expsig;
9c499471 197 info[1]= ResolutiondEdxTR(&ncls, respar);
5cd0300d 198 }
199
bebe75fc 200 const Double_t eps = 1e-10;
201
5cd0300d 202 if(ratio){
bebe75fc 203 return track->GetTRDsignal()/(expsig + eps);
5cd0300d 204 }
205 else{
206 return track->GetTRDsignal() - expsig;
207 }
208}
209
210
2ddf4abb 211Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Float_t * par)
5cd0300d 212{
213 //
214 //resolution
215 //npar=3
216 //
217
218 const Double_t ncls = xx[0];
219
220 return par[0]+par[1]*TMath::Power(ncls, par[2]);
221}
222
2ddf4abb 223Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Float_t * pin)
5cd0300d 224{
225 //
226 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
227 //npar = 8 = 3+5
228 //
229
2ddf4abb 230 Float_t ptr[4]={0};
5cd0300d 231 for(int ii=0; ii<3; ii++){
232 ptr[ii+1]=pin[ii];
233 }
234 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
235}
236
2ddf4abb 237Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Float_t * par)
5cd0300d 238{
239 //
240 //LOGISTIC parametrization for TR, in unit of MIP
241 //npar = 4
242 //
243
244 const Double_t bg = xx[0];
245 const Double_t gamma = sqrt(1+bg*bg);
246
247 const Double_t p0 = TMath::Abs(par[1]);
248 const Double_t p1 = TMath::Abs(par[2]);
249 const Double_t p2 = TMath::Abs(par[3]);
250
251 const Double_t zz = TMath::Log(gamma);
252 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
253
254 return par[0]+tryield;
255}
256
2ddf4abb 257Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Float_t * par)
5cd0300d 258{
259 //
260 //ALEPH parametrization for dEdx
261 //npar = 5
262 //
263
264 const Double_t bg = xx[0];
265 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
266
267 const Double_t p0 = TMath::Abs(par[0]);
268 const Double_t p1 = TMath::Abs(par[1]);
269
270 const Double_t p2 = TMath::Abs(par[2]);
271
272 const Double_t p3 = TMath::Abs(par[3]);
273 const Double_t p4 = TMath::Abs(par[4]);
274
275 const Double_t aa = TMath::Power(beta, p3);
276
277 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
278
279 return (p1-aa-bb)*p0/aa;
280
281}
282
283
ffb1ee30 284//____________________________________________________________
239fe91c 285Int_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 286{
e0de37e9 287 //
288 // Calculate TRD likelihood values for the track based on dedx and
289 // momentum values. The likelihoods are calculated by query the
290 // reference data depending on the PID method selected
291 //
292 // Input parameter :
293 // n - number of dedx slices/chamber
294 // dedx - array of dedx slices organized layer wise
295 // p - array of momentum measurements organized layer wise
296 //
297 // Return parameters
298 // prob - probabilities allocated by TRD for particle specis
299 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
300 //
301 // Return value
239fe91c 302 // number of tracklets used for PID, 0 if no PID
bd58d4b9 303 //
304 AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
305
ffb1ee30 306
db0e2c5f 307 if(!fkPIDResponseObject){
31746593 308 AliDebug(3,"Missing reference data. PID calculation not possible.");
239fe91c 309 return 0;
db0e2c5f 310 }
ffb1ee30 311
db0e2c5f 312 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
313 Double_t prLayer[AliPID::kSPECIES];
314 Double_t dE[10], s(0.);
239fe91c 315 Int_t ntrackletsPID=0;
db0e2c5f 316 for(Int_t il(kNlayer); il--;){
317 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
bd58d4b9 318 if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
239fe91c 319 s=0.;
db0e2c5f 320 Bool_t filled=kTRUE;
321 for(Int_t is(AliPID::kSPECIES); is--;){
239fe91c 322 //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
bd58d4b9 323 if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
f2762b1c 324 AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
db0e2c5f 325 if(prLayer[is]<1.e-30){
326 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
327 filled=kFALSE;
328 break;
329 }
330 s+=prLayer[is];
331 }
332 if(!filled){
333 continue;
334 }
335 if(s<1.e-30){
336 AliDebug(2, Form("Null all species prob layer[%d].", il));
337 continue;
338 }
339 for(Int_t is(AliPID::kSPECIES); is--;){
340 if(kNorm) prLayer[is] /= s;
341 prob[is] *= prLayer[is];
342 }
239fe91c 343 ntrackletsPID++;
e0de37e9 344 }
239fe91c 345 if(!kNorm) return ntrackletsPID;
db0e2c5f 346
347 s=0.;
348 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
e0de37e9 349 if(s<1.e-30){
db0e2c5f 350 AliDebug(2, "Null total prob.");
239fe91c 351 return 0;
e0de37e9 352 }
db0e2c5f 353 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
239fe91c 354 return ntrackletsPID;
ffb1ee30 355}
356
ffb1ee30 357//____________________________________________________________
bd58d4b9 358Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
e0de37e9 359 //
360 // Get the non-normalized probability for a certain particle species as coming
361 // from the reference histogram
362 // Interpolation between momentum bins
363 //
364 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
db0e2c5f 365
51a0ce25 366 Double_t probLayer = 0.;
239fe91c 367
db0e2c5f 368 Float_t pLower, pUpper;
369
239fe91c 370 AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
371 *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
372
373 // Do Interpolation exept for underflow and overflow
374 if(refLower && refUpper){
375 Double_t probLower = refLower->Eval(dEdx);
376 Double_t probUpper = refUpper->Eval(dEdx);
377
378 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
379 } else if(refLower){
380 // underflow
381 probLayer = refLower->Eval(dEdx);
382 } else if(refUpper){
383 // overflow
384 probLayer = refUpper->Eval(dEdx);
385 } else {
386 AliError("No references available");
387 }
388 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
389
390 return probLayer;
391
392/* old implementation
393
394switch(PIDmethod){
395case kNN: // NN
db0e2c5f 396 break;
397 case kLQ2D: // 2D LQ
398 {
f2762b1c 399 if(species==0||species==2){ // references only for electrons and pions
2285168d 400 Double_t error = 0.;
f2762b1c 401 Double_t point[kNslicesLQ2D];
402 for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
db0e2c5f 403
239fe91c 404 AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
405 *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
406 // Do Interpolation exept for underflow and overflow
407 if(refLower && refUpper){
408 Double_t probLower=0,probUpper=0;
409 refLower->Eval(point,probLower,error);
410 refUpper->Eval(point,probUpper,error);
411 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
412 } else if(refLower){
413 // underflow
f2762b1c 414 refLower->Eval(point,probLayer,error);
239fe91c 415 } else if(refUpper){
416 // overflow
417 refUpper->Eval(point,probLayer,error);
418 } else {
f2762b1c 419 AliError("No references available");
420 }
421 AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
db0e2c5f 422 }
423 }
424 break;
425 case kLQ1D: // 1D LQ
426 {
427 TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
428 *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
429 // Do Interpolation exept for underflow and overflow
430 if(refLower && refUpper){
431 Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
432 Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
433
434 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
435 } else if(refLower){
436 // underflow
437 probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
438 } else if(refUpper){
439 // overflow
440 probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
441 } else {
442 AliError("No references available");
443 }
f2762b1c 444 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
db0e2c5f 445 }
446 break;
447 default:
448 break;
239fe91c 449 }
450 return probLayer;
451 */
452
ffb1ee30 453}
454
455//____________________________________________________________
456void AliTRDPIDResponse::SetOwner(){
e0de37e9 457 //
458 // Make Deep Copy of the Reference Histograms
459 //
db0e2c5f 460 if(!fkPIDResponseObject || IsOwner()) return;
461 const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
462 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
463 SetBit(kIsOwner, kTRUE);
ffb1ee30 464}
465
466//____________________________________________________________
bd58d4b9 467Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
ffb1ee30 468{
f2762b1c 469 //
470 // Recalculate dE/dx
471 //
bd58d4b9 472 switch(PIDmethod){
e0de37e9 473 case kNN: // NN
db0e2c5f 474 break;
475 case kLQ2D: // 2D LQ
476 out[0]=0;
477 out[1]=0;
478 for(Int_t islice = 0; islice < nSlice; islice++){
239fe91c 479 if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
db0e2c5f 480 if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
481 else out[1]+= in[islice];
482 }
239fe91c 483 // normalize signal to number of slices
484 out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
485 out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
f2762b1c 486 if(out[0] < 1e-6) return kFALSE;
487 AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
db0e2c5f 488 break;
489 case kLQ1D: // 1D LQ
490 out[0]= 0.;
491 for(Int_t islice = 0; islice < nSlice; islice++)
492 if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
239fe91c 493 out[0]*=1./Double_t(nSlice);
db0e2c5f 494 if(out[0] < 1e-6) return kFALSE;
f2762b1c 495 AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
db0e2c5f 496 break;
497
e0de37e9 498 default:
db0e2c5f 499 return kFALSE;
e0de37e9 500 }
501 return kTRUE;
ffb1ee30 502}
503
ea235c90 504//____________________________________________________________
bd58d4b9 505Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
db0e2c5f 506 //
507 // Check whether particle is identified as electron assuming a certain electron efficiency level
508 // Only electron and pion hypothesis is taken into account
509 //
510 // Inputs:
511 // Number of tracklets
512 // Likelihood values
513 // Momentum
514 // Electron efficiency level
515 //
516 // If the function fails when the params are not accessible, the function returns true
517 //
518 if(!fkPIDResponseObject){
31746593 519 AliDebug(3,"No PID Param object available");
ea235c90 520 return kTRUE;
521 }
522 Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
ce487a7f 523 Double_t params[4];
bd58d4b9 524 if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
ea235c90 525 AliError("No Params found for the given configuration");
526 return kTRUE;
527 }
ce487a7f 528 Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
ea235c90 529 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
530 return kFALSE;
531}
532
f2762b1c 533//____________________________________________________________
534Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
535
536 fkPIDResponseObject = obj;
537 if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
538 return kTRUE;
539}
9c499471 540
541
542//____________________________________________________________
543Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
544{
545 fkTRDdEdxParams = par;
546 return kTRUE;
547}