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 }