1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
22 // Markus Fasel <M.Fasel@gsi.de>
23 // Anton Andronic <A.Andronic@gsi.de>
27 #include <TDirectory.h>
32 #include <TObjArray.h>
39 #include "AliVTrack.h"
41 #include "AliTRDPIDResponseObject.h"
42 #include "AliTRDPIDResponse.h"
43 //#include "AliTRDTKDInterpolator.h"
44 #include "AliTRDNDFast.h"
45 #include "AliTRDdEdxParams.h"
47 ClassImp(AliTRDPIDResponse)
49 //____________________________________________________________
50 AliTRDPIDResponse::AliTRDPIDResponse():
52 ,fkPIDResponseObject(NULL)
53 ,fkTRDdEdxParams(NULL)
54 ,fGainNormalisationFactor(1.)
57 // Default constructor
61 //____________________________________________________________
62 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
64 ,fkPIDResponseObject(NULL)
65 ,fkTRDdEdxParams(NULL)
66 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
73 //____________________________________________________________
74 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
76 // Assignment operator
78 if(this == &ref) return *this;
81 TObject::operator=(ref);
82 fGainNormalisationFactor = ref.fGainNormalisationFactor;
83 fkPIDResponseObject = ref.fkPIDResponseObject;
84 fkTRDdEdxParams = ref.fkTRDdEdxParams;
89 //____________________________________________________________
90 AliTRDPIDResponse::~AliTRDPIDResponse(){
95 delete fkPIDResponseObject;
96 delete fkTRDdEdxParams;
100 //____________________________________________________________
101 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
103 // Load References into the toolkit
105 AliDebug(1, "Loading reference histos from root file");
106 TDirectory *owd = gDirectory;// store old working directory
109 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
110 TFile *in = TFile::Open(filename);
112 AliError("Ref file not available.");
117 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
118 in->Close(); delete in;
120 SetBit(kIsOwner, kTRUE);
121 AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
127 //____________________________________________________________
128 Double_t AliTRDPIDResponse::GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType type) const
131 //calculate the TRD nSigma
134 const Double_t badval = -9999;
135 Double_t info[5]; for(int i=0; i<5; i++){info[i]=badval;}
137 const Double_t delta = GetSignalDelta(track, type, kFALSE, info);
139 const Double_t mean = info[0];
140 const Double_t res = info[1];
145 const Double_t sigma = mean*res;
146 const Double_t eps = 1e-12;
147 return delta/(sigma + eps);
150 //____________________________________________________________
151 Double_t AliTRDPIDResponse::GetSignalDelta( const AliVTrack* track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/, Double_t *info/*=0x0*/) const
154 //calculate the TRD signal difference w.r.t. the expected
155 //output other information in info[]
158 const Double_t badval = -9999;
160 //cut on number of chambers
161 const Double_t cutch = 6;
163 //cut on mean number of clusters per chamber
164 const Double_t cutclsperch = 18;
166 const Double_t nch = track->GetTRDNchamberdEdx();
171 const Double_t ncls = track->GetTRDNclusterdEdx();
172 if( ncls/nch < cutclsperch){
176 Double_t pTRD = -999;
177 for(Int_t ich=0; ich<6; ich++){
178 pTRD = track->GetTRDmomentum(ich);
186 if(!fkTRDdEdxParams){
187 AliError("AliTRDPIDResponse::GetSignalDelta fkTRDdEdxParams null");
191 TVectorF meanvec = fkTRDdEdxParams->GetMeanParameter(type);
192 TVectorF resvec = fkTRDdEdxParams->GetSigmaParameter(type);
194 Double_t meanpar[meanvec.GetNrows()];
195 Double_t respar[resvec.GetNrows()];
197 for(Int_t ii=0; ii<meanvec.GetNrows(); ii++){
198 meanpar[ii]=meanvec[ii];
200 for(Int_t ii=0; ii<resvec.GetNrows(); ii++){
201 respar[ii]=resvec[ii];
204 //============================================================================================<<<<<<<<<<<<<
206 const Double_t bg = pTRD/AliPID::ParticleMass(type);
207 const Double_t expsig = MeandEdxTR(&bg, meanpar);
211 info[1]= ResolutiondEdxTR(&ncls, respar);
214 const Double_t eps = 1e-10;
217 return track->GetTRDsignal()/(expsig + eps);
220 return track->GetTRDsignal() - expsig;
225 Double_t AliTRDPIDResponse::ResolutiondEdxTR(const Double_t * xx, const Double_t * par)
232 const Double_t ncls = xx[0];
234 return par[0]+par[1]*TMath::Power(ncls, par[2]);
237 Double_t AliTRDPIDResponse::MeandEdxTR(const Double_t * xx, const Double_t * pin)
240 //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
245 for(int ii=0; ii<3; ii++){
248 return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
251 Double_t AliTRDPIDResponse::MeanTR(const Double_t * xx, const Double_t * par)
254 //LOGISTIC parametrization for TR, in unit of MIP
258 const Double_t bg = xx[0];
259 const Double_t gamma = sqrt(1+bg*bg);
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]);
265 const Double_t zz = TMath::Log(gamma);
266 const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
268 return par[0]+tryield;
271 Double_t AliTRDPIDResponse::MeandEdx(const Double_t * xx, const Double_t * par)
274 //ALEPH parametrization for dEdx
278 const Double_t bg = xx[0];
279 const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
281 const Double_t p0 = TMath::Abs(par[0]);
282 const Double_t p1 = TMath::Abs(par[1]);
284 const Double_t p2 = TMath::Abs(par[2]);
286 const Double_t p3 = TMath::Abs(par[3]);
287 const Double_t p4 = TMath::Abs(par[4]);
289 const Double_t aa = TMath::Power(beta, p3);
291 const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
293 return (p1-aa-bb)*p0/aa;
298 //____________________________________________________________
299 Int_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
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
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
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.
316 // number of tracklets used for PID, 0 if no PID
318 AliDebug(3,Form(" Response for PID method: %d",PIDmethod));
321 if(!fkPIDResponseObject){
322 AliDebug(3,"Missing reference data. PID calculation not possible.");
326 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
327 Double_t prLayer[AliPID::kSPECIES];
328 Double_t dE[10], s(0.);
329 Int_t ntrackletsPID=0;
330 for(Int_t il(kNlayer); il--;){
331 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
332 if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue;
335 for(Int_t is(AliPID::kSPECIES); is--;){
336 //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue;
337 if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod);
338 AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is]));
339 if(prLayer[is]<1.e-30){
340 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
350 AliDebug(2, Form("Null all species prob layer[%d].", il));
353 for(Int_t is(AliPID::kSPECIES); is--;){
354 if(kNorm) prLayer[is] /= s;
355 prob[is] *= prLayer[is];
359 if(!kNorm) return ntrackletsPID;
362 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
364 AliDebug(2, "Null total prob.");
367 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
368 return ntrackletsPID;
371 //____________________________________________________________
372 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) const {
374 // Get the non-normalized probability for a certain particle species as coming
375 // from the reference histogram
376 // Interpolation between momentum bins
378 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
380 Double_t probLayer = 0.;
382 Float_t pLower, pUpper;
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));
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);
392 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
395 probLayer = refLower->Eval(dEdx);
398 probLayer = refUpper->Eval(dEdx);
400 AliError("No references available");
402 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
406 /* old implementation
413 if(species==0||species==2){ // references only for electrons and pions
415 Double_t point[kNslicesLQ2D];
416 for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
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);
428 refLower->Eval(point,probLayer,error);
431 refUpper->Eval(point,probLayer,error);
433 AliError("No references available");
435 AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error));
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]));
448 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
451 probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
454 probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
456 AliError("No references available");
458 AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer));
469 //____________________________________________________________
470 void AliTRDPIDResponse::SetOwner(){
472 // Make Deep Copy of the Reference Histograms
474 if(!fkPIDResponseObject || IsOwner()) return;
475 const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
476 fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
477 SetBit(kIsOwner, kTRUE);
480 //____________________________________________________________
481 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const
492 for(Int_t islice = 0; islice < nSlice; islice++){
493 if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled
494 if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
495 else out[1]+= in[islice];
497 // normalize signal to number of slices
498 out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0());
499 out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0());
500 if(out[0] < 1e-6) return kFALSE;
501 AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1]));
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
507 out[0]*=1./Double_t(nSlice);
508 if(out[0] < 1e-6) return kFALSE;
509 AliDebug(3,Form("CookdEdx dEdx %f",out[0]));
518 //____________________________________________________________
519 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const {
521 // Check whether particle is identified as electron assuming a certain electron efficiency level
522 // Only electron and pion hypothesis is taken into account
525 // Number of tracklets
528 // Electron efficiency level
530 // If the function fails when the params are not accessible, the function returns true
532 if(!fkPIDResponseObject){
533 AliDebug(3,"No PID Param object available");
536 Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
538 if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){
539 AliError("No Params found for the given configuration");
542 Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
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
547 //____________________________________________________________
548 Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){
550 fkPIDResponseObject = obj;
551 if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print("");
556 //____________________________________________________________
557 Bool_t AliTRDPIDResponse::SetdEdxParams(const AliTRDdEdxParams * par)
559 fkTRDdEdxParams = par;