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