]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEERBase/AliPIDCombined.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Here the comment:
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDCombined.cxx
CommitLineData
b8bfee30 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
17//-----------------------------------------------------------------
18// Base class for combining PID of different detectors //
19// (user selected) and compute Bayesian probabilities //
20// //
21// //
22// Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch //
23// //
24//-----------------------------------------------------------------
25
26#include <TH1.h>
27
28#include <AliVTrack.h>
29#include <AliLog.h>
30#include <AliPID.h>
31#include <AliPIDResponse.h>
32
33#include "AliPIDCombined.h"
34
80f28562 35#include "TMath.h"
36#include "TFile.h"
37
38#include "AliOADBContainer.h"
39
40// initialize static members
00a38d07 41TH2F* AliPIDCombined::fDefaultPriorsTPC[]={0x0};
80f28562 42
43ClassImp(AliPIDCombined);
44
b8bfee30 45AliPIDCombined::AliPIDCombined():
46 TNamed("CombinedPID","CombinedPID"),
47 fDetectorMask(0),
48 fRejectMismatchMask(0x7F),
49 fEnablePriors(kTRUE),
80f28562 50 fSelectedSpecies(AliPID::kSPECIES),
51 fUseDefaultTPCPriors(kFALSE)
b8bfee30 52{
53 //
54 // default constructor
55 //
00a38d07 56 for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
b8bfee30 57 AliLog::SetClassDebugLevel("AliPIDCombined",10);
58}
59
60//-------------------------------------------------------------------------------------------------
61AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
62 TNamed(name,title),
63 fDetectorMask(0),
64 fRejectMismatchMask(0x7F),
65 fEnablePriors(kTRUE),
80f28562 66 fSelectedSpecies(AliPID::kSPECIES),
67 fUseDefaultTPCPriors(kFALSE)
b8bfee30 68{
69 //
70 // default constructor with name and title
71 //
00a38d07 72 for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
b8bfee30 73 AliLog::SetClassDebugLevel("AliPIDCombined",10);
74
75}
76
77//-------------------------------------------------------------------------------------------------
78AliPIDCombined::~AliPIDCombined() {
79
00a38d07 80 for(Int_t i=0;i < AliPID::kSPECIESC;i++){
b8bfee30 81 if(fPriorsDistributions[i])
82 delete fPriorsDistributions[i];
83 }
84}
85
86//-------------------------------------------------------------------------------------------------
87void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
00a38d07 88 if ( (type < 0) || ( type >= ((AliPID::EParticleType)AliPID::kSPECIESC) ) ){
b8bfee30 89 AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
90 return;
91 }
92 if(prior) {
8de3b790 93 Int_t i = (Int_t)type;
b8bfee30 94 if (fPriorsDistributions[i]) {
95 delete fPriorsDistributions[i];
96 }
97 fPriorsDistributions[i]=new TH1F(*prior);
98 }
99}
100
101//-------------------------------------------------------------------------------------------------
102UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
103 Double_t pITS[fSelectedSpecies];
104 Double_t pTPC[fSelectedSpecies];
105 Double_t pTRD[fSelectedSpecies];
106 Double_t pTOF[fSelectedSpecies];
107 Double_t pHMPID[fSelectedSpecies];
108 Double_t pEMCAL[fSelectedSpecies];
109 Double_t pPHOS[fSelectedSpecies];
110 for (Int_t i=0;i<fSelectedSpecies;i++) {
111 pITS[i]=1.;
112 pTPC[i]=1.;
113 pTRD[i]=1.;
114 pTOF[i]=1.;
115 pHMPID[i]=1.;
116 pEMCAL[i]=1.;
117 pPHOS[i]=1.;
118 }
119 UInt_t usedMask=0;
120 AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
121 Double_t p[fSelectedSpecies];
4c234df8 122 memset(p,0,sizeof(Double_t)*fSelectedSpecies);
b8bfee30 123
124 // getting probability distributions for selected detectors only
125 if (fDetectorMask & AliPIDResponse::kDetITS) {
126 status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
127 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
128 }
129
130 if (fDetectorMask & AliPIDResponse::kDetTPC) {
131 status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
132 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
133 }
134
135
136 if (fDetectorMask & AliPIDResponse::kDetTRD) {
137 status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
138 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
139 }
140
141 if (fDetectorMask & AliPIDResponse::kDetTOF) {
142 status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
143 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
144 }
145
146 if (fDetectorMask & AliPIDResponse::kDetHMPID) {
147 status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
148 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
149 }
150
151
152 if (fDetectorMask & AliPIDResponse::kDetEMCAL) {
153 status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
154 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
155 }
156
157
158 if (fDetectorMask & AliPIDResponse::kDetPHOS) {
159 status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
160 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
161 }
162
163
164 for (Int_t i =0; i<fSelectedSpecies; i++){
165 p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
166 }
167 Double_t priors[fSelectedSpecies];
4c234df8 168 memset(priors,0,fSelectedSpecies*sizeof(Double_t));
80f28562 169 if (fEnablePriors){
170 GetPriors(track,priors,response->GetCurrentCentrality());
171
44bea91d 172 // for the moment we have four cases
df0293b3 173 // TPC+TRD --> apply TRD propagation factors
44bea91d 174 // TPC+TOF --> apply TOF propagation factors (TRD may be present)
175 // TPC+EMCAL --> apply EMCAL propagation factors (TRD and TOF if present)
df0293b3 176 if(fUseDefaultTPCPriors) {
177 Double_t pt=TMath::Abs(track->Pt());
44bea91d 178 if ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL ) { // EMCAL is the outer having prop. factors for the moment
179 // EMCal case (for the moment only in combination with TPC)
180 // propagation factors determined from LHC11d MC (LHC12a15f)
181 // v2 clusterizer, dEta < 0.015, dPhi < 0.03, NonLinearityFunction = 6
182
183 Float_t electronEMCALfactor=0.1;
184 Float_t kaonEMCALfactor = 0.1;
185 Float_t protonEMCALfactor = 0.1;
186
187 if(track->Charge() > 0){
188 // positiv charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
189 if(pt > 0.75 && pt < 20.0){
190 electronEMCALfactor = ( 0.214159 * ( 1 - TMath::Exp(-TMath::Power(pt,0.484512)/0.700499)*TMath::Power(pt,-0.669644)) ) / ( 0.210436 * ( 1 - TMath::Exp(-TMath::Power(pt,-0.219228)/0.947432)*TMath::Power(pt,-0.700792)) );
191 kaonEMCALfactor = ( 0.208686 * ( 1 - TMath::Exp(-TMath::Power(pt,-3.98149e-05)/1.28447)*TMath::Power(pt,-0.629191)) ) / ( 0.210436 * ( 1 - TMath::Exp(-TMath::Power(pt,-0.219228)/0.947432)*TMath::Power(pt,-0.700792)) );
192 protonEMCALfactor = ( 0.27555 * ( 1 - TMath::Exp(-TMath::Power(pt,-1.97226e-05)/1.52719)*TMath::Power(pt,-0.209068)) ) / ( 0.210436 * ( 1 - TMath::Exp(-TMath::Power(pt,-0.219228)/0.947432)*TMath::Power(pt,-0.700792)) );
193
194 }
195 }
196
197 if(track->Charge() < 0){
198 // negative charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
199 if(pt > 0.75 && pt < 20.0){
200 electronEMCALfactor = ( 0.216895 * ( 1 - TMath::Exp(-TMath::Power(pt,0.000105924)/0.865938)*TMath::Power(pt,-1.32787)) ) / ( 0.210385 * ( 1 - TMath::Exp(-TMath::Power(pt,4.41206e-07)/1.08984)*TMath::Power(pt,-0.544375)) );
201 kaonEMCALfactor = ( 0.204117 * ( 1 - TMath::Exp(-TMath::Power(pt,-1.6853e-05)/1.61765)*TMath::Power(pt,-0.738355)) ) / ( 0.210385 * ( 1 - TMath::Exp(-TMath::Power(pt,4.41206e-07)/1.08984)*TMath::Power(pt,-0.544375)) );
202 protonEMCALfactor = ( 0.215679 * ( 1 - TMath::Exp(-TMath::Power(pt,-4.10015e-05)/1.40921)*TMath::Power(pt,-0.533752)) ) / ( 0.210385 * ( 1 - TMath::Exp(-TMath::Power(pt,4.41206e-07)/1.08984)*TMath::Power(pt,-0.544375)) );
203 }
204 }
205 p[0] *= electronEMCALfactor;
206 p[3] *= kaonEMCALfactor;
207 p[4] *= protonEMCALfactor;
208 } // end of EMCAL case
209 else if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){
df0293b3 210 Float_t kaonTOFfactor = 0.1;
211 if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
212 Float_t protonTOFfactor = 0.1;
213 if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
214
215 if(track->Charge() < 0){ // for negative tracks
216 kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
217 protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
218 }
44bea91d 219 // TODO: we may need an electron factor for TOF as well, especially if TRD is present!
df0293b3 220 p[3] *= kaonTOFfactor;
221 p[4] *= protonTOFfactor;
44bea91d 222 } // end of TOF case
223 else if ( (usedMask & AliPIDResponse::kDetTRD)==AliPIDResponse::kDetTRD ) {
df0293b3 224 Float_t electronTRDfactor=0.1;
225 Float_t kaonTRDfactor = 0.1;
226 Float_t protonTRDfactor = 0.1;
227
228 if(track->Charge() > 0){
229 // positiv charge
230 if(pt > 0.25) electronTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,5.13315e-03)/2.11145e-01)*TMath::Power(pt,-2.97659e+00);
231 if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,-4.29549e-02)/4.87989e-01)*TMath::Power(pt,-1.54270e+00);
232 if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.81238e+00)/7.57082e-02)*TMath::Power(pt,-8.12595e-01);
233 }
234
235 if(track->Charge() < 0){
236 // negative charge
237 if(pt > 0.25) electronTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.45537e-02)/1.90397e-01)*TMath::Power(pt,-3.33121e+00);
238 if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt, -3.42831e-03)/5.57013e-01)*TMath::Power(pt,-1.39202e+00);
239 if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,3.36631e+00)/7.18819e-02)*TMath::Power(pt,-2.00577e-01);
240 }
241 // what about electrons
242 p[0] *= electronTRDfactor;
243 p[3] *= kaonTRDfactor;
244 p[4] *= protonTRDfactor;
44bea91d 245 } // end of TRD case
246 } // end of fUseDefaultTPCPriors
df0293b3 247 } // end of use priors
b8bfee30 248 else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
249 ComputeBayesProbabilities(bayesProbabilities,p,priors);
250 return usedMask;
251}
252
253
254//-------------------------------------------------------------------------------------------------
80f28562 255void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
b8bfee30 256
257 //
258 // get priors from distributions
259 //
260
261 Double_t pt=TMath::Abs(track->Pt());
80f28562 262
263 if(fUseDefaultTPCPriors){ // use default priors if requested
264 Float_t usedCentr = centrality+5*(centrality>0);
265 if(usedCentr < -0.99) usedCentr = -0.99;
266 else if(usedCentr > 99) usedCentr = 99;
267 if(pt > 9.99) pt = 9.99;
268 else if(pt < 0.01) pt = 0.01;
269
270 for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
271
272 return;
273 }
274
b8bfee30 275 Double_t sumPriors = 0;
276 for (Int_t i=0;i<fSelectedSpecies;++i){
277 if (fPriorsDistributions[i]){
278 p[i]=fPriorsDistributions[i]->Interpolate(pt);
279 }
280 else {
281 p[i]=0.;
282 }
283 sumPriors+=p[i];
284 }
285
286 // normalizing priors
287 if (sumPriors == 0) return;
288 for (Int_t i=0;i<fSelectedSpecies;++i){
289 p[i]=p[i]/sumPriors;
290 }
291 return;
292}
80f28562 293//-------------------------------------------------------------------------------------------------
294void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
295
296 //
297 // get priors from distributions
298 //
299 Double_t pt=TMath::Abs(track->Pt());
300
301 if(fUseDefaultTPCPriors){ // use default priors if requested
302 Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
303 if(usedCentr < -0.99) usedCentr = -0.99;
304 else if(usedCentr > 99) usedCentr = 99;
305 if(pt > 9.99) pt = 9.99;
306 else if(pt < 0.01) pt = 0.01;
307
308 for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
309
310 // Extra factor if TOF matching was required
311 if(detUsed & AliPIDResponse::kDetTOF){
312 Float_t kaonTOFfactor = 0.1;
313 if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
314 Float_t protonTOFfactor = 0.1;
315 if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
316
317 if(track->Charge() < 0){ // for negative tracks
318 kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
319 protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
320 }
321
322 p[3] *= kaonTOFfactor;
323 p[4] *= protonTOFfactor;
324 }
325
326 return;
327 }
328
329
330 Double_t sumPriors = 0;
331 for (Int_t i=0;i<fSelectedSpecies;++i){
332 if (fPriorsDistributions[i]){
333 p[i]=fPriorsDistributions[i]->Interpolate(pt);
334 }
335 else {
336 p[i]=0.;
337 }
338 sumPriors+=p[i];
339 }
b8bfee30 340
80f28562 341 // normalizing priors
342 if (sumPriors == 0) return;
343 for (Int_t i=0;i<fSelectedSpecies;++i){
344 p[i]=p[i]/sumPriors;
345 }
346 return;
347}
b8bfee30 348//-------------------------------------------------------------------------------------------------
349void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
350
351
352 //
353 // calculate Bayesian probabilities
354 //
355 Double_t sum = 0.;
356 for (Int_t i = 0; i < fSelectedSpecies; i++) {
357 sum += probDensity[i] * prior[i];
358 }
359 if (sum <= 0) {
80f28562 360
b8bfee30 361 AliError("Invalid probability densities or priors");
362 for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
363 return;
364 }
365 for (Int_t i = 0; i < fSelectedSpecies; i++) {
366 probabilities[i] = probDensity[i] * prior[i] / sum;
367 }
368
369
370}
371
372//----------------------------------------------------------------------------------------
373void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
374 switch (status) {
375 case AliPIDResponse::kDetNoSignal:
376 break;
377 case AliPIDResponse::kDetPidOk:
378 *maskDetIn = *maskDetIn | bit;
379 break;
380 case AliPIDResponse::kDetMismatch:
381 if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
382 else *maskDetIn = *maskDetIn | bit;
383 break;
384 }
385
386}
387
80f28562 388//----------------------------------------------------------------------------------------
389void AliPIDCombined::SetDefaultTPCPriors(){
390 fEnablePriors=kTRUE;
391 fUseDefaultTPCPriors = kTRUE;
392
00a38d07 393 //check if priors are already initialized
394 if (fDefaultPriorsTPC[0]) return;
395
80f28562 396 TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
397 oadbfilename += "/PIDdefaultPriors.root";
398 TFile * foadb = TFile::Open(oadbfilename.Data());
00a38d07 399 if(!foadb || !foadb->IsOpen()) {
400 delete foadb;
401 AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
402 return;
403 }
80f28562 404
405 AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
406 if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
407
00a38d07 408 const char *namespecies[AliPID::kSPECIES] = {"El","Mu","Pi","Ka","Pr"};
80f28562 409 char name[100];
410
00a38d07 411 for(Int_t i=0;i < AliPID::kSPECIES;i++){
4c234df8 412 snprintf(name,99,"hDefault%sPriors",namespecies[i]);
80f28562 413 fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
414 if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
415 }
00a38d07 416
417 delete foadb;
80f28562 418}