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 **************************************************************************/
17 //-----------------------------------------------------------------
18 // Base class for combining PID of different detectors //
19 // (user selected) and compute Bayesian probabilities //
22 // Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch //
24 //-----------------------------------------------------------------
28 #include <AliVTrack.h>
31 #include <AliPIDResponse.h>
33 #include "AliPIDCombined.h"
38 #include "AliOADBContainer.h"
40 // initialize static members
41 TH2F* AliPIDCombined::fDefaultPriorsTPC[5];
43 ClassImp(AliPIDCombined);
45 AliPIDCombined::AliPIDCombined():
46 TNamed("CombinedPID","CombinedPID"),
48 fRejectMismatchMask(0x7F),
50 fSelectedSpecies(AliPID::kSPECIES),
51 fUseDefaultTPCPriors(kFALSE)
54 // default constructor
56 for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
57 AliLog::SetClassDebugLevel("AliPIDCombined",10);
60 //-------------------------------------------------------------------------------------------------
61 AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
64 fRejectMismatchMask(0x7F),
66 fSelectedSpecies(AliPID::kSPECIES),
67 fUseDefaultTPCPriors(kFALSE)
70 // default constructor with name and title
72 for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
73 AliLog::SetClassDebugLevel("AliPIDCombined",10);
77 //-------------------------------------------------------------------------------------------------
78 AliPIDCombined::~AliPIDCombined() {
80 for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
81 if(fPriorsDistributions[i])
82 delete fPriorsDistributions[i];
86 //-------------------------------------------------------------------------------------------------
87 void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
88 if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
89 AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
93 Int_t i = (Int_t)type;
94 if (i >= AliPID::kSPECIES ) {
95 if (i < (Int_t)AliPID::kDeuteron) {
96 AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
98 } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
100 if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) { // coverity fix (to put it mildly....)
101 AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
104 if (fPriorsDistributions[i]) {
105 delete fPriorsDistributions[i];
107 fPriorsDistributions[i]=new TH1F(*prior);
111 //-------------------------------------------------------------------------------------------------
112 UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
113 Double_t pITS[fSelectedSpecies];
114 Double_t pTPC[fSelectedSpecies];
115 Double_t pTRD[fSelectedSpecies];
116 Double_t pTOF[fSelectedSpecies];
117 Double_t pHMPID[fSelectedSpecies];
118 Double_t pEMCAL[fSelectedSpecies];
119 Double_t pPHOS[fSelectedSpecies];
120 for (Int_t i=0;i<fSelectedSpecies;i++) {
130 AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
131 Double_t p[fSelectedSpecies];
132 memset(p,0,sizeof(Double_t)*fSelectedSpecies);
134 // getting probability distributions for selected detectors only
135 if (fDetectorMask & AliPIDResponse::kDetITS) {
136 status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
137 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
140 if (fDetectorMask & AliPIDResponse::kDetTPC) {
141 status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
142 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
146 if (fDetectorMask & AliPIDResponse::kDetTRD) {
147 status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
148 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
151 if (fDetectorMask & AliPIDResponse::kDetTOF) {
152 status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
153 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
156 if (fDetectorMask & AliPIDResponse::kDetHMPID) {
157 status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
158 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
162 if (fDetectorMask & AliPIDResponse::kDetEMCAL) {
163 status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
164 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
168 if (fDetectorMask & AliPIDResponse::kDetPHOS) {
169 status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
170 SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
174 for (Int_t i =0; i<fSelectedSpecies; i++){
175 p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
177 Double_t priors[fSelectedSpecies];
178 memset(priors,0,fSelectedSpecies*sizeof(Double_t));
180 GetPriors(track,priors,response->GetCurrentCentrality());
182 // apply tof matching efficiency to priors if TOF joined PID for this track
183 if(fUseDefaultTPCPriors && (usedMask & AliPIDResponse::kDetTOF)){
184 Double_t pt=TMath::Abs(track->Pt());
185 Float_t kaonTOFfactor = 0.1;
186 if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
187 Float_t protonTOFfactor = 0.1;
188 if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
190 if(track->Charge() < 0){ // for negative tracks
191 kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
192 protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
195 p[3] *= kaonTOFfactor;
196 p[4] *= protonTOFfactor;
199 else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
200 ComputeBayesProbabilities(bayesProbabilities,p,priors);
205 //-------------------------------------------------------------------------------------------------
206 void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
209 // get priors from distributions
212 Double_t pt=TMath::Abs(track->Pt());
214 if(fUseDefaultTPCPriors){ // use default priors if requested
215 Float_t usedCentr = centrality+5*(centrality>0);
216 if(usedCentr < -0.99) usedCentr = -0.99;
217 else if(usedCentr > 99) usedCentr = 99;
218 if(pt > 9.99) pt = 9.99;
219 else if(pt < 0.01) pt = 0.01;
221 for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
226 Double_t sumPriors = 0;
227 for (Int_t i=0;i<fSelectedSpecies;++i){
228 if (fPriorsDistributions[i]){
229 p[i]=fPriorsDistributions[i]->Interpolate(pt);
237 // normalizing priors
238 if (sumPriors == 0) return;
239 for (Int_t i=0;i<fSelectedSpecies;++i){
244 //-------------------------------------------------------------------------------------------------
245 void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
248 // get priors from distributions
250 Double_t pt=TMath::Abs(track->Pt());
252 if(fUseDefaultTPCPriors){ // use default priors if requested
253 Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
254 if(usedCentr < -0.99) usedCentr = -0.99;
255 else if(usedCentr > 99) usedCentr = 99;
256 if(pt > 9.99) pt = 9.99;
257 else if(pt < 0.01) pt = 0.01;
259 for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
261 // Extra factor if TOF matching was required
262 if(detUsed & AliPIDResponse::kDetTOF){
263 Float_t kaonTOFfactor = 0.1;
264 if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
265 Float_t protonTOFfactor = 0.1;
266 if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
268 if(track->Charge() < 0){ // for negative tracks
269 kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
270 protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
273 p[3] *= kaonTOFfactor;
274 p[4] *= protonTOFfactor;
281 Double_t sumPriors = 0;
282 for (Int_t i=0;i<fSelectedSpecies;++i){
283 if (fPriorsDistributions[i]){
284 p[i]=fPriorsDistributions[i]->Interpolate(pt);
292 // normalizing priors
293 if (sumPriors == 0) return;
294 for (Int_t i=0;i<fSelectedSpecies;++i){
299 //-------------------------------------------------------------------------------------------------
300 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
304 // calculate Bayesian probabilities
307 for (Int_t i = 0; i < fSelectedSpecies; i++) {
308 sum += probDensity[i] * prior[i];
312 AliError("Invalid probability densities or priors");
313 for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
316 for (Int_t i = 0; i < fSelectedSpecies; i++) {
317 probabilities[i] = probDensity[i] * prior[i] / sum;
323 //----------------------------------------------------------------------------------------
324 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
326 case AliPIDResponse::kDetNoSignal:
328 case AliPIDResponse::kDetPidOk:
329 *maskDetIn = *maskDetIn | bit;
331 case AliPIDResponse::kDetMismatch:
332 if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
333 else *maskDetIn = *maskDetIn | bit;
339 //----------------------------------------------------------------------------------------
340 void AliPIDCombined::SetDefaultTPCPriors(){
342 fUseDefaultTPCPriors = kTRUE;
344 TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
345 oadbfilename += "/PIDdefaultPriors.root";
346 TFile * foadb = TFile::Open(oadbfilename.Data());
347 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
349 AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
350 if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
352 const char *namespecies[5] = {"El","Mu","Pi","Ka","Pr"};
355 for(Int_t i=0;i < 5;i++){
356 snprintf(name,99,"hDefault%sPriors",namespecies[i]);
357 fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
358 if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));