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