TRD nSigma OADB related new codes and modifications and OADB root file -- Xianguo Lu
[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
160 //cut on number of chambers
161 const Double_t cutch = 6;
162
163 //cut on mean number of clusters per chamber
164 const Double_t cutclsperch = 18;
5cd0300d 165
166 const Double_t nch = track->GetTRDNchamberdEdx();
167 if(nch < cutch){
168 return badval;
169 }
170
171 const Double_t ncls = track->GetTRDNclusterdEdx();
172 if( ncls/nch < cutclsperch){
173 return badval;
174 }
175
176 Double_t pTRD = -999;
177 for(Int_t ich=0; ich<6; ich++){
178 pTRD = track->GetTRDmomentum(ich);
179 if(pTRD>0)
180 break;
181 }
182 if(pTRD<0){
183 return badval;
184 }
185
9c499471 186 if(!fkTRDdEdxParams){
187 AliError("AliTRDPIDResponse::GetSignalDelta fkTRDdEdxParams null");
188 return -99999;
189 }
190
191 TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type);
192 TVectorF resvec = fkTRDdEdxParams->GetSigmaParameter(type);
193
194 Double_t meanpar[meanvec.GetNrows()];
195 Double_t respar[resvec.GetNrows()];
196
197 for(Int_t ii=0; ii<meanvec.GetNrows(); ii++){
198 meanpar[ii]=meanvec[ii];
199 }
200 for(Int_t ii=0; ii<resvec.GetNrows(); ii++){
201 respar[ii]=resvec[ii];
202 }
203
204 //============================================================================================<<<<<<<<<<<<<
205
5cd0300d 206 const Double_t bg = pTRD/AliPID::ParticleMass(type);
9c499471 207 const Double_t expsig = MeandEdxTR(&bg, meanpar);
5cd0300d 208
209 if(info){
210 info[0]= expsig;
9c499471 211 info[1]= ResolutiondEdxTR(&ncls, respar);
5cd0300d 212 }
213
bebe75fc 214 const Double_t eps = 1e-10;
215
5cd0300d 216 if(ratio){
bebe75fc 217 return track->GetTRDsignal()/(expsig + eps);
5cd0300d 218 }
219 else{
220 return track->GetTRDsignal() - expsig;
221 }
222}
223
224
225Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Double_t * par)
226{
227 //
228 //resolution
229 //npar=3
230 //
231
232 const Double_t ncls = xx[0];
233
234 return par[0]+par[1]*TMath::Power(ncls, par[2]);
235}
236
237Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Double_t * pin)
238{
239 //
240 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
241 //npar = 8 = 3+5
242 //
243
244 Double_t ptr[4]={0};
245 for(int ii=0; ii<3; ii++){
246 ptr[ii+1]=pin[ii];
247 }
248 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
249}
250
251Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Double_t * par)
252{
253 //
254 //LOGISTIC parametrization for TR, in unit of MIP
255 //npar = 4
256 //
257
258 const Double_t bg = xx[0];
259 const Double_t gamma = sqrt(1+bg*bg);
260
261 const Double_t p0 = TMath::Abs(par[1]);
262 const Double_t p1 = TMath::Abs(par[2]);
263 const Double_t p2 = TMath::Abs(par[3]);
264
265 const Double_t zz = TMath::Log(gamma);
266 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
267
268 return par[0]+tryield;
269}
270
271Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Double_t * par)
272{
273 //
274 //ALEPH parametrization for dEdx
275 //npar = 5
276 //
277
278 const Double_t bg = xx[0];
279 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
280
281 const Double_t p0 = TMath::Abs(par[0]);
282 const Double_t p1 = TMath::Abs(par[1]);
283
284 const Double_t p2 = TMath::Abs(par[2]);
285
286 const Double_t p3 = TMath::Abs(par[3]);
287 const Double_t p4 = TMath::Abs(par[4]);
288
289 const Double_t aa = TMath::Power(beta, p3);
290
291 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
292
293 return (p1-aa-bb)*p0/aa;
294
295}
296
297
ffb1ee30 298//____________________________________________________________
239fe91c 299Int_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 300{
e0de37e9 301 //
302 // Calculate TRD likelihood values for the track based on dedx and
303 // momentum values. The likelihoods are calculated by query the
304 // reference data depending on the PID method selected
305 //
306 // Input parameter :
307 // n - number of dedx slices/chamber
308 // dedx - array of dedx slices organized layer wise
309 // p - array of momentum measurements organized layer wise
310 //
311 // Return parameters
312 // prob - probabilities allocated by TRD for particle specis
313 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
314 //
315 // Return value
239fe91c 316 // number of tracklets used for PID, 0 if no PID
bd58d4b9 317 //
318 AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
319
ffb1ee30 320
db0e2c5f 321 if(!fkPIDResponseObject){
31746593 322 AliDebug(3,"Missing reference data. PID calculation not possible.");
239fe91c 323 return 0;
db0e2c5f 324 }
ffb1ee30 325
db0e2c5f 326 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
327 Double_t prLayer[AliPID::kSPECIES];
328 Double_t dE[10], s(0.);
239fe91c 329 Int_t ntrackletsPID=0;
db0e2c5f 330 for(Int_t il(kNlayer); il--;){
331 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
bd58d4b9 332 if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
239fe91c 333 s=0.;
db0e2c5f 334 Bool_t filled=kTRUE;
335 for(Int_t is(AliPID::kSPECIES); is--;){
239fe91c 336 //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
bd58d4b9 337 if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
f2762b1c 338 AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
db0e2c5f 339 if(prLayer[is]<1.e-30){
340 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
341 filled=kFALSE;
342 break;
343 }
344 s+=prLayer[is];
345 }
346 if(!filled){
347 continue;
348 }
349 if(s<1.e-30){
350 AliDebug(2, Form("Null all species prob layer[%d].", il));
351 continue;
352 }
353 for(Int_t is(AliPID::kSPECIES); is--;){
354 if(kNorm) prLayer[is] /= s;
355 prob[is] *= prLayer[is];
356 }
239fe91c 357 ntrackletsPID++;
e0de37e9 358 }
239fe91c 359 if(!kNorm) return ntrackletsPID;
db0e2c5f 360
361 s=0.;
362 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
e0de37e9 363 if(s<1.e-30){
db0e2c5f 364 AliDebug(2, "Null total prob.");
239fe91c 365 return 0;
e0de37e9 366 }
db0e2c5f 367 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
239fe91c 368 return ntrackletsPID;
ffb1ee30 369}
370
ffb1ee30 371//____________________________________________________________
bd58d4b9 372Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
e0de37e9 373 //
374 // Get the non-normalized probability for a certain particle species as coming
375 // from the reference histogram
376 // Interpolation between momentum bins
377 //
378 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
db0e2c5f 379
51a0ce25 380 Double_t probLayer = 0.;
239fe91c 381
db0e2c5f 382 Float_t pLower, pUpper;
383
239fe91c 384 AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)),
385 *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod));
386
387 // Do Interpolation exept for underflow and overflow
388 if(refLower && refUpper){
389 Double_t probLower = refLower->Eval(dEdx);
390 Double_t probUpper = refUpper->Eval(dEdx);
391
392 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
393 } else if(refLower){
394 // underflow
395 probLayer = refLower->Eval(dEdx);
396 } else if(refUpper){
397 // overflow
398 probLayer = refUpper->Eval(dEdx);
399 } else {
400 AliError("No references available");
401 }
402 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
403
404 return probLayer;
405
406/* old implementation
407
408switch(PIDmethod){
409case kNN: // NN
db0e2c5f 410 break;
411 case kLQ2D: // 2D LQ
412 {
f2762b1c 413 if(species==0||species==2){ // references only for electrons and pions
2285168d 414 Double_t error = 0.;
f2762b1c 415 Double_t point[kNslicesLQ2D];
416 for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
db0e2c5f 417
239fe91c 418 AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)),
419 *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
420 // Do Interpolation exept for underflow and overflow
421 if(refLower && refUpper){
422 Double_t probLower=0,probUpper=0;
423 refLower->Eval(point,probLower,error);
424 refUpper->Eval(point,probUpper,error);
425 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
426 } else if(refLower){
427 // underflow
f2762b1c 428 refLower->Eval(point,probLayer,error);
239fe91c 429 } else if(refUpper){
430 // overflow
431 refUpper->Eval(point,probLayer,error);
432 } else {
f2762b1c 433 AliError("No references available");
434 }
435 AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
db0e2c5f 436 }
437 }
438 break;
439 case kLQ1D: // 1D LQ
440 {
441 TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
442 *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
443 // Do Interpolation exept for underflow and overflow
444 if(refLower && refUpper){
445 Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
446 Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
447
448 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
449 } else if(refLower){
450 // underflow
451 probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
452 } else if(refUpper){
453 // overflow
454 probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
455 } else {
456 AliError("No references available");
457 }
f2762b1c 458 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
db0e2c5f 459 }
460 break;
461 default:
462 break;
239fe91c 463 }
464 return probLayer;
465 */
466
ffb1ee30 467}
468
469//____________________________________________________________
470void AliTRDPIDResponse::SetOwner(){
e0de37e9 471 //
472 // Make Deep Copy of the Reference Histograms
473 //
db0e2c5f 474 if(!fkPIDResponseObject || IsOwner()) return;
475 const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
476 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
477 SetBit(kIsOwner, kTRUE);
ffb1ee30 478}
479
480//____________________________________________________________
bd58d4b9 481Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
ffb1ee30 482{
f2762b1c 483 //
484 // Recalculate dE/dx
485 //
bd58d4b9 486 switch(PIDmethod){
e0de37e9 487 case kNN: // NN
db0e2c5f 488 break;
489 case kLQ2D: // 2D LQ
490 out[0]=0;
491 out[1]=0;
492 for(Int_t islice = 0; islice < nSlice; islice++){
239fe91c 493 if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
db0e2c5f 494 if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
495 else out[1]+= in[islice];
496 }
239fe91c 497 // normalize signal to number of slices
498 out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
499 out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
f2762b1c 500 if(out[0] < 1e-6) return kFALSE;
501 AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
db0e2c5f 502 break;
503 case kLQ1D: // 1D LQ
504 out[0]= 0.;
505 for(Int_t islice = 0; islice < nSlice; islice++)
506 if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
239fe91c 507 out[0]*=1./Double_t(nSlice);
db0e2c5f 508 if(out[0] < 1e-6) return kFALSE;
f2762b1c 509 AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
db0e2c5f 510 break;
511
e0de37e9 512 default:
db0e2c5f 513 return kFALSE;
e0de37e9 514 }
515 return kTRUE;
ffb1ee30 516}
517
ea235c90 518//____________________________________________________________
bd58d4b9 519Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
db0e2c5f 520 //
521 // Check whether particle is identified as electron assuming a certain electron efficiency level
522 // Only electron and pion hypothesis is taken into account
523 //
524 // Inputs:
525 // Number of tracklets
526 // Likelihood values
527 // Momentum
528 // Electron efficiency level
529 //
530 // If the function fails when the params are not accessible, the function returns true
531 //
532 if(!fkPIDResponseObject){
31746593 533 AliDebug(3,"No PID Param object available");
ea235c90 534 return kTRUE;
535 }
536 Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
ce487a7f 537 Double_t params[4];
bd58d4b9 538 if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
ea235c90 539 AliError("No Params found for the given configuration");
540 return kTRUE;
541 }
ce487a7f 542 Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
ea235c90 543 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
544 return kFALSE;
545}
546
f2762b1c 547//____________________________________________________________
548Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
549
550 fkPIDResponseObject = obj;
551 if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
552 return kTRUE;
553}
9c499471 554
555
556//____________________________________________________________
557Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
558{
559 fkTRDdEdxParams = par;
560 return kTRUE;
561}