9fa2be88448c5f5753058ebd98013a1080f74aed
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDCombined.cxx
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
35 #include "TMath.h"
36 #include "TFile.h"
37
38 #include "AliOADBContainer.h"
39
40 // initialize static members
41 TH2F* AliPIDCombined::fDefaultPriorsTPC[]={0x0};
42 Float_t AliPIDCombined::fTOFmismProb = 0;
43
44 ClassImp(AliPIDCombined);
45
46 AliPIDCombined::AliPIDCombined():
47         TNamed("CombinedPID","CombinedPID"),
48         fDetectorMask(0),
49         fRejectMismatchMask(0x7F),
50         fEnablePriors(kTRUE),
51         fSelectedSpecies(AliPID::kSPECIES),
52         fUseDefaultTPCPriors(kFALSE)
53 {       
54   //
55   // default constructor
56   //    
57   for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
58   AliLog::SetClassDebugLevel("AliPIDCombined",10);
59 }
60
61 //-------------------------------------------------------------------------------------------------     
62 AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
63         TNamed(name,title),
64         fDetectorMask(0),
65         fRejectMismatchMask(0x7F),
66         fEnablePriors(kTRUE),
67         fSelectedSpecies(AliPID::kSPECIES),
68         fUseDefaultTPCPriors(kFALSE)
69 {
70   //
71   // default constructor with name and title
72   //
73   for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
74   AliLog::SetClassDebugLevel("AliPIDCombined",10);
75
76 }
77
78 //-------------------------------------------------------------------------------------------------     
79 AliPIDCombined::~AliPIDCombined() {
80
81   for(Int_t i=0;i < AliPID::kSPECIESC;i++){
82     if(fPriorsDistributions[i])
83       delete fPriorsDistributions[i];
84   }
85 }
86
87 //-------------------------------------------------------------------------------------------------     
88 void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
89   if ( (type < 0) || ( type >= ((AliPID::EParticleType)AliPID::kSPECIESC) ) ){
90     AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
91     return;
92   }
93   if(prior) {
94     Int_t i = (Int_t)type;
95     if (fPriorsDistributions[i]) {
96       delete fPriorsDistributions[i]; 
97     }
98     fPriorsDistributions[i]=new TH1F(*prior);
99   }
100 }
101
102 //-------------------------------------------------------------------------------------------------     
103 UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
104   //
105   // (1) Get raw probabilities of requested detectors and combine
106   // (2) Get priors and propagate depending on detectors used
107   // (3) Compute Bayes probabilities
108   //
109
110
111   // (1) Get raw probabilities of selected detectors and combine
112   UInt_t usedMask=0;             // detectors actually used for track
113   fTOFmismProb = 0; // reset TOF mismatch weights
114
115   AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
116   Double_t p[fSelectedSpecies];  // combined probabilities of selected detectors
117   Double_t pMismTOF[fSelectedSpecies];  // combined TOF mismatch probabilities using selected detectors
118   for (Int_t i=0;i<fSelectedSpecies;i++){ p[i]=1.;pMismTOF[i]=1.;} // no decision
119   for (Int_t ibit = 0; ibit < 7 ; ibit++) {
120     AliPIDResponse::EDetCode detBit = (AliPIDResponse::EDetCode)(1<<ibit);
121     if (fDetectorMask & detBit) {           // getting probabilities for requested detectors only
122       Double_t detProb[fSelectedSpecies];
123       status = response->ComputePIDProbability(detBit,track,fSelectedSpecies,detProb);
124       if (status == AliPIDResponse::kDetPidOk) {
125         if (fRejectMismatchMask & detBit) {         // mismatch check (currently just for TOF)
126           if (detBit == AliPIDResponse::kDetTOF) {
127             fTOFmismProb = response->GetTOFMismatchProbability(); // mismatch weights computed with TOF probs (no arguments)
128             Float_t probMis = response->GetTOFMismatchProbability(track); // mismatch compatibility TPC-TOF cut
129             SetCombinedStatus(status,&usedMask,detBit,detProb,probMis);
130           } else {
131             SetCombinedStatus(status,&usedMask,detBit);
132           }
133         } else {
134           SetCombinedStatus(status,&usedMask,detBit);
135         }
136         for (Int_t i=0;i<fSelectedSpecies;i++){
137           p[i] *= detProb[i];
138           if(detBit == AliPIDResponse::kDetTOF) pMismTOF[i] *= fTOFmismProb;
139           else pMismTOF[i] *= detProb[i];
140         }
141       }
142     }
143   }
144   // if no detectors available there is no point to go further
145   if (usedMask == 0) return usedMask;
146
147   // (2) Get priors and propagate depending on detectors used
148   Double_t priors[fSelectedSpecies];
149   memset(priors,0,fSelectedSpecies*sizeof(Double_t));
150   if (fEnablePriors){
151     GetPriors(track,priors,response->GetCurrentCentrality());
152
153     // We apply the propagation factors of the more external detector
154     // 
155     // TPC+HMPID    --> apply HMPID propagation factors (TRD and TOF may be present)
156     // TPC+EMCAL    --> apply EMCAL propagation factors (TRD and TOF may be present)
157     // TPC+TOF      --> apply TOF propagation factors (TRD may be present, HMPID and EMCAL not (if requested))
158     // TPC+TRD      --> apply TRD propagation factors (TOF, HMPID and EMCAL not present (if requested) )
159     // 
160     if(fUseDefaultTPCPriors) {
161
162       Double_t pt=TMath::Abs(track->Pt());
163       if ( ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL) || ( (usedMask & AliPIDResponse::kDetHMPID)==AliPIDResponse::kDetHMPID) ) {
164       // we assume EMCAL and HMPID cannot be simultaneously present
165       if ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL ) {
166         // EMCal case (for the moment only in combination with TPC)
167         // propagation factors determined from LHC11d MC (LHC12a15f)
168         // v2 clusterizer, dEta < 0.015, dPhi < 0.03, NonLinearityFunction = 6
169                
170         Float_t electronEMCALfactor=0.1;
171         Float_t kaonEMCALfactor = 0.1;
172         Float_t protonEMCALfactor = 0.1;
173                
174         if(track->Charge() > 0){
175           // positiv charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
176           if(pt > 0.75 && pt < 20.0){
177             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)) );
178             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)) );
179             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)) );
180                    
181           }
182         }
183                
184         if(track->Charge() < 0){
185           // negative charge  (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
186           if(pt > 0.75 && pt < 20.0){              
187             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)) );
188             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)) );
189             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)) );
190           }
191         }
192         priors[0] *= electronEMCALfactor;
193         priors[3] *= kaonEMCALfactor;
194         priors[4] *= protonEMCALfactor;       
195       } // end of EMCAL case
196
197       else if ( (usedMask & AliPIDResponse::kDetHMPID)==AliPIDResponse::kDetHMPID ) {  // HMPID case
198
199         Float_t kaonHMPIDfactor   = 0.;
200         Float_t protonHMPIDfactor = 0.;                     
201         if(pt>1. && pt<6.)  kaonHMPIDfactor   = (-0.0729337 + pt*0.0999531 - pt*pt*0.0371803 + pt*pt*pt*0.00706436 - pt*pt*pt*pt*0.000643619 + pt*pt*pt*pt*pt*2.21853e-05)/(-0.00896231+ pt*0.0330702 - pt*pt*0.0109562+ pt*pt*pt*0.00232895 - pt*pt*pt*pt*0.000246143 + pt*pt*pt*pt*pt*9.59812e-06); 
202         if(pt>1.4 && pt<6.) protonHMPIDfactor = (-0.0444188 + pt*0.0681506 - pt*pt*0.0231819 + pt*pt*pt*0.00400771 - pt*pt*pt*pt*0.000339315 + pt*pt*pt*pt*pt*1.12616e-05)/(-0.00896231+ pt*0.0330702 - pt*pt*0.0109562+ pt*pt*pt*0.00232895 - pt*pt*pt*pt*0.000246143 + pt*pt*pt*pt*pt*9.59812e-06);        
203         if(pt>6. && pt<8.)  kaonHMPIDfactor   = (-0.0729337 + pt*0.0999531 - pt*pt*0.0371803 + pt*pt*pt*0.00706436 - pt*pt*pt*pt*0.000643619 + pt*pt*pt*pt*pt*2.21853e-05)/0.0530456;        
204         if(pt>8.)           kaonHMPIDfactor   = 0.0550432/0.0530456; 
205         if(pt>6. && pt<8.5) protonHMPIDfactor = (-0.0444188 + pt*0.0681506 - pt*pt*0.0231819 + pt*pt*pt*0.00400771 - pt*pt*pt*pt*0.000339315 + pt*pt*pt*pt*pt*1.12616e-05)/0.0530456;      
206         if(pt>8.5)          protonHMPIDfactor = 0.0530071/0.0530456;       
207                                  
208         if(track->Charge() < 0){ 
209          if(pt>0.4 && pt<6.) protonHMPIDfactor = (-0.0351485 + pt*0.0473821 - pt*pt*0.0147947 + pt*pt*pt*0.00254811- pt*pt*pt*pt*0.000224724 + pt*pt*pt*pt*pt*7.9303e-06)/(-0.00896231+ pt*0.0330702 - pt*pt*0.0109562+ pt*pt*pt*0.00232895 - pt*pt*pt*pt*0.000246143 + pt*pt*pt*pt*pt*9.59812e-06);
210          if(pt>6. && pt<8.5) protonHMPIDfactor = (-0.0351485 + pt*0.0473821 - pt*pt*0.0147947 + pt*pt*pt*0.00254811- pt*pt*pt*pt*0.000224724 + pt*pt*pt*pt*pt*7.9303e-06)/0.0530456; 
211          if(pt>8.5)          protonHMPIDfactor = 0.0457756/0.0530456; 
212         } 
213       
214         priors[3] *= kaonHMPIDfactor;
215         priors[4] *= protonHMPIDfactor;
216                
217       }
218
219     } // end of outer cases: EMCAL/HMPID
220       else if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){
221         Float_t kaonTOFfactor = 0.1;
222         if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
223         Float_t protonTOFfactor = 0.1;
224         if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
225                
226         if(track->Charge() < 0){ // for negative tracks
227           kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
228           protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
229         }
230         // TODO: we may need an electron factor for TOF as well, especially if TRD is present!
231         priors[3] *= kaonTOFfactor;
232         priors[4] *= protonTOFfactor;
233       }  // end of TOF case
234       else if ( (usedMask & AliPIDResponse::kDetTRD)==AliPIDResponse::kDetTRD ) {
235         Float_t electronTRDfactor=0.1;
236         Float_t kaonTRDfactor = 0.1;
237         Float_t protonTRDfactor = 0.1;
238                  
239         if(track->Charge() > 0){
240           // positiv charge
241           if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,5.13315e-03)/2.11145e-01)*TMath::Power(pt,-2.97659e+00);
242           if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,-4.29549e-02)/4.87989e-01)*TMath::Power(pt,-1.54270e+00);
243           if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.81238e+00)/7.57082e-02)*TMath::Power(pt,-8.12595e-01);
244         }
245                  
246         if(track->Charge() < 0){
247           // negative charge
248           if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,2.45537e-02)/1.90397e-01)*TMath::Power(pt,-3.33121e+00);
249           if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt, -3.42831e-03)/5.57013e-01)*TMath::Power(pt,-1.39202e+00);
250           if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,3.36631e+00)/7.18819e-02)*TMath::Power(pt,-2.00577e-01);
251         }
252         // what about electrons
253         priors[0] *= electronTRDfactor;
254         priors[3] *= kaonTRDfactor;
255         priors[4] *= protonTRDfactor;         
256       } // end of TRD case
257     } // end of fUseDefaultTPCPriors
258   }   // end of use priors
259   else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
260
261   // (3) Compute Bayes probabilities
262   ComputeBayesProbabilities(bayesProbabilities,p,priors);
263
264   // compute TOF probability contribution from mismatch
265   fTOFmismProb = 0; 
266   for (Int_t i=0;i<fSelectedSpecies;i++) fTOFmismProb += pMismTOF[i];
267
268   return usedMask;
269 }
270
271
272 //-------------------------------------------------------------------------------------------------
273 void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
274         
275         //
276         // get priors from distributions
277         //
278         
279         Double_t pt=TMath::Abs(track->Pt());
280
281         if(fUseDefaultTPCPriors){ // use default priors if requested
282           Float_t usedCentr = centrality+5*(centrality>0);
283           if(usedCentr < -0.99) usedCentr = -0.99;
284           else if(usedCentr > 99) usedCentr = 99;
285           if(pt > 9.99) pt = 9.99;
286           else if(pt < 0.01)  pt = 0.01;
287           
288           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
289
290           return;
291         }
292         
293         Double_t sumPriors = 0;
294         for (Int_t i=0;i<fSelectedSpecies;++i){
295                 if (fPriorsDistributions[i]){
296                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
297                 }
298                 else {
299                         p[i]=0.;
300                 }
301                 sumPriors+=p[i];                
302         }
303
304         // normalizing priors
305         if (sumPriors == 0) return;
306         for (Int_t i=0;i<fSelectedSpecies;++i){
307            p[i]=p[i]/sumPriors;
308         }
309         return;
310 }
311 //-------------------------------------------------------------------------------------------------
312 void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
313         
314         //
315         // get priors from distributions
316         //
317         Double_t pt=TMath::Abs(track->Pt());
318         
319         if(fUseDefaultTPCPriors){ // use default priors if requested
320           Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
321           if(usedCentr < -0.99) usedCentr = -0.99;
322           else if(usedCentr > 99) usedCentr = 99;
323           if(pt > 9.99) pt = 9.99;
324           else if(pt < 0.01)  pt = 0.01;
325           
326           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
327
328           // Extra factor if TOF matching was required
329           if(detUsed & AliPIDResponse::kDetTOF){
330             Float_t kaonTOFfactor = 0.1;
331             if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
332             Float_t protonTOFfactor = 0.1;
333             if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
334             
335             if(track->Charge() < 0){ // for negative tracks
336               kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
337               protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
338             }
339             
340             p[3] *= kaonTOFfactor;
341             p[4] *= protonTOFfactor;
342           }
343
344           return;
345         }
346
347
348         Double_t sumPriors = 0;
349         for (Int_t i=0;i<fSelectedSpecies;++i){
350                 if (fPriorsDistributions[i]){
351                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
352                 }
353                 else {
354                         p[i]=0.;
355                 }
356                 sumPriors+=p[i];                
357         }
358
359         // normalizing priors
360         if (sumPriors == 0) return;
361         for (Int_t i=0;i<fSelectedSpecies;++i){
362            p[i]=p[i]/sumPriors;
363         }
364         return;
365 }
366 //-------------------------------------------------------------------------------------------------     
367 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior, Double_t* probDensityMism) const {
368
369
370   //
371   // calculate Bayesian probabilities
372   //
373   Double_t sum = 0.;
374   for (Int_t i = 0; i < fSelectedSpecies; i++) {
375     sum += probDensity[i] * prior[i];
376   }
377   if (sum <= 0) {
378
379     AliError("Invalid probability densities or priors");
380     for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
381     return;
382   }
383   for (Int_t i = 0; i < fSelectedSpecies; i++) {
384     probabilities[i] = probDensity[i] * prior[i] / sum;
385     if(probDensityMism) probDensityMism[i] *= prior[i] / sum;
386   }
387
388
389 }
390
391
392 //----------------------------------------------------------------------------------------
393 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit) const {
394   switch (status) {
395   case AliPIDResponse::kDetNoParams:
396   case AliPIDResponse::kDetNoSignal:
397   case AliPIDResponse::kDetMismatch: // for backward compatibilty, we need then to remove kDetMismatch from AliPIDResponse
398     break;
399   case AliPIDResponse::kDetPidOk:
400     *maskDetIn = *maskDetIn | bit;
401     break;
402   }
403
404 }
405
406 //----------------------------------------------------------------------------------------
407 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p, const Float_t probMis) const {
408   switch (status) {
409   case AliPIDResponse::kDetNoParams:
410   case AliPIDResponse::kDetNoSignal:
411   case AliPIDResponse::kDetMismatch: // for backward compatibility, we need then to remove kDeteMismatch from AliPIDResponse
412     break;
413   case AliPIDResponse::kDetPidOk:
414       if (probMis > 0.5) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
415       else *maskDetIn = *maskDetIn | bit;
416     break;
417   }
418
419 }
420
421
422
423 //----------------------------------------------------------------------------------------
424 void AliPIDCombined::SetDefaultTPCPriors(){
425   fEnablePriors=kTRUE;
426   fUseDefaultTPCPriors = kTRUE;
427
428   //check if priors are already initialized
429   if (fDefaultPriorsTPC[0]) return;
430   
431   TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
432   oadbfilename += "/PIDdefaultPriors.root";
433   TFile * foadb = TFile::Open(oadbfilename.Data());
434   if(!foadb || !foadb->IsOpen()) {
435     delete foadb;
436     AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
437     return;
438   }
439   
440   AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
441   if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
442   
443   const char *namespecies[AliPID::kSPECIES] = {"El","Mu","Pi","Ka","Pr"};
444   char name[100];
445
446   for(Int_t i=0;i < AliPID::kSPECIES;i++){
447     snprintf(name,99,"hDefault%sPriors",namespecies[i]);
448     fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
449     if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
450   }
451
452   delete foadb;
453 }