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