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[5];
  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::kSPECIES+AliPID::kSPECIESLN;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::kSPECIES+AliPID::kSPECIESLN;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::kSPECIES+AliPID::kSPECIESLN;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::kSPECIESN+AliPID::kSPECIESLN) ) ){
  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 (i >= AliPID::kSPECIES ) {
  95       if (i < (Int_t)AliPID::kDeuteron) {
  96         AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
  97         return;
  98       } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
  99     }
 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));
 102       return;
 103     }
 104     if (fPriorsDistributions[i]) {
 105       delete fPriorsDistributions[i];
 106     }
 107     fPriorsDistributions[i]=new TH1F(*prior);
 108   }
 109 }
 110
 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++) {
 121          pITS[i]=1.;
 122          pTPC[i]=1.;
 123          pTRD[i]=1.;
 124          pTOF[i]=1.;
 125          pHMPID[i]=1.;
 126          pEMCAL[i]=1.;
 127          pPHOS[i]=1.;
 128         }
 129         UInt_t usedMask=0;
 130         AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
 131         Double_t p[fSelectedSpecies];
 132         memset(p,0,sizeof(Double_t)*fSelectedSpecies);
 133
 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);
 138         }
 139
 140         if (fDetectorMask & AliPIDResponse::kDetTPC) {
 141           status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
 142           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
 143         }
 144
 145
 146         if (fDetectorMask & AliPIDResponse::kDetTRD) {
 147           status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
 148           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
 149         }
 150
 151         if (fDetectorMask &  AliPIDResponse::kDetTOF) {
 152           status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
 153           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
 154         }
 155
 156         if (fDetectorMask & AliPIDResponse::kDetHMPID) {
 157           status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
 158           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
 159         }
 160
 161
 162         if (fDetectorMask & AliPIDResponse::kDetEMCAL) {
 163           status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
 164           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
 165         }
 166
 167
 168         if (fDetectorMask & AliPIDResponse::kDetPHOS) {
 169           status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
 170           SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
 171         }
 172
 173
 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];
 176         }
 177         Double_t priors[fSelectedSpecies];
 178         memset(priors,0,fSelectedSpecies*sizeof(Double_t));
 179         if (fEnablePriors){
 180           GetPriors(track,priors,response->GetCurrentCentrality());
 181
 182           // for the moment we have three cases
 183           // TPC+TRD      --> apply TRD propagation factors
 184           // TPC+TOF      --> apply TOF propagation factors
 185           // TPC+TRD+TOF  --> apply TOF propagation factors
 186           // apply tof matching efficiency to priors if TOF joined PID for this track
 187           if(fUseDefaultTPCPriors) {
 188              Double_t pt=TMath::Abs(track->Pt());
 189              if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){ // TOF is the outer having prop. factors for the moment
 190                Float_t kaonTOFfactor = 0.1;
 191                if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
 192                Float_t protonTOFfactor = 0.1;
 193                if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
 194
 195                if(track->Charge() < 0){ // for negative tracks
 196                  kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
 197                  protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
 198                }
 199                p[3] *= kaonTOFfactor;
 200                p[4] *= protonTOFfactor;
 201              }
 202              else {
 203                if ( (usedMask & AliPIDResponse::kDetTRD)==AliPIDResponse::kDetTRD ) {
 204                  Float_t electronTRDfactor=0.1;
 205                  Float_t kaonTRDfactor = 0.1;
 206                  Float_t protonTRDfactor = 0.1;
 207
 208                  if(track->Charge() > 0){
 209                    // positiv charge
 210                    if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,5.13315e-03)/2.11145e-01)*TMath::Power(pt,-2.97659e+00);
 211                    if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,-4.29549e-02)/4.87989e-01)*TMath::Power(pt,-1.54270e+00);
 212                    if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.81238e+00)/7.57082e-02)*TMath::Power(pt,-8.12595e-01);
 213                  }
 214
 215                  if(track->Charge() < 0){
 216                    // negative charge
 217                    if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,2.45537e-02)/1.90397e-01)*TMath::Power(pt,-3.33121e+00);
 218                    if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt, -3.42831e-03)/5.57013e-01)*TMath::Power(pt,-1.39202e+00);
 219                    if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,3.36631e+00)/7.18819e-02)*TMath::Power(pt,-2.00577e-01);
 220                  }
 221                  // what about electrons
 222                  p[0] *= electronTRDfactor;
 223                  p[3] *= kaonTRDfactor;
 224                  p[4] *= protonTRDfactor;
 225                } // end of TRD case
 226              } // end of detectors inner than TOF
 227           } // end of fUseDefaultTPCPriors
 228         }   // end of use priors
 229         else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
 230         ComputeBayesProbabilities(bayesProbabilities,p,priors);
 231         return usedMask;
 232 }
 233
 234
 235 //-------------------------------------------------------------------------------------------------
 236 void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
 237
 238         //
 239         // get priors from distributions
 240         //
 241
 242         Double_t pt=TMath::Abs(track->Pt());
 243
 244         if(fUseDefaultTPCPriors){ // use default priors if requested
 245           Float_t usedCentr = centrality+5*(centrality>0);
 246           if(usedCentr < -0.99) usedCentr = -0.99;
 247           else if(usedCentr > 99) usedCentr = 99;
 248           if(pt > 9.99) pt = 9.99;
 249           else if(pt < 0.01)  pt = 0.01;
 250
 251           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
 252
 253           return;
 254         }
 255
 256         Double_t sumPriors = 0;
 257         for (Int_t i=0;i<fSelectedSpecies;++i){
 258                 if (fPriorsDistributions[i]){
 259                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
 260                 }
 261                 else {
 262                         p[i]=0.;
 263                 }
 264                 sumPriors+=p[i];
 265         }
 266
 267         // normalizing priors
 268         if (sumPriors == 0) return;
 269         for (Int_t i=0;i<fSelectedSpecies;++i){
 270            p[i]=p[i]/sumPriors;
 271         }
 272         return;
 273 }
 274 //-------------------------------------------------------------------------------------------------
 275 void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
 276
 277         //
 278         // get priors from distributions
 279         //
 280         Double_t pt=TMath::Abs(track->Pt());
 281
 282         if(fUseDefaultTPCPriors){ // use default priors if requested
 283           Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
 284           if(usedCentr < -0.99) usedCentr = -0.99;
 285           else if(usedCentr > 99) usedCentr = 99;
 286           if(pt > 9.99) pt = 9.99;
 287           else if(pt < 0.01)  pt = 0.01;
 288
 289           for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
 290
 291           // Extra factor if TOF matching was required
 292           if(detUsed & AliPIDResponse::kDetTOF){
 293             Float_t kaonTOFfactor = 0.1;
 294             if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
 295             Float_t protonTOFfactor = 0.1;
 296             if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
 297
 298             if(track->Charge() < 0){ // for negative tracks
 299               kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
 300               protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
 301             }
 302
 303             p[3] *= kaonTOFfactor;
 304             p[4] *= protonTOFfactor;
 305           }
 306
 307           return;
 308         }
 309
 310
 311         Double_t sumPriors = 0;
 312         for (Int_t i=0;i<fSelectedSpecies;++i){
 313                 if (fPriorsDistributions[i]){
 314                         p[i]=fPriorsDistributions[i]->Interpolate(pt);
 315                 }
 316                 else {
 317                         p[i]=0.;
 318                 }
 319                 sumPriors+=p[i];
 320         }
 321
 322         // normalizing priors
 323         if (sumPriors == 0) return;
 324         for (Int_t i=0;i<fSelectedSpecies;++i){
 325            p[i]=p[i]/sumPriors;
 326         }
 327         return;
 328 }
 329 //-------------------------------------------------------------------------------------------------
 330 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
 331
 332
 333   //
 334   // calculate Bayesian probabilities
 335   //
 336   Double_t sum = 0.;
 337   for (Int_t i = 0; i < fSelectedSpecies; i++) {
 338     sum += probDensity[i] * prior[i];
 339   }
 340   if (sum <= 0) {
 341
 342     AliError("Invalid probability densities or priors");
 343     for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
 344     return;
 345   }
 346   for (Int_t i = 0; i < fSelectedSpecies; i++) {
 347     probabilities[i] = probDensity[i] * prior[i] / sum;
 348   }
 349
 350
 351 }
 352
 353 //----------------------------------------------------------------------------------------
 354 void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p) const {
 355   switch (status) {
 356   case AliPIDResponse::kDetNoSignal:
 357     break;
 358   case AliPIDResponse::kDetPidOk:
 359     *maskDetIn = *maskDetIn | bit;
 360     break;
 361   case AliPIDResponse::kDetMismatch:
 362     if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
 363     else *maskDetIn = *maskDetIn | bit;
 364     break;
 365   }
 366
 367 }
 368
 369 //----------------------------------------------------------------------------------------
 370 void AliPIDCombined::SetDefaultTPCPriors(){
 371   fEnablePriors=kTRUE;
 372   fUseDefaultTPCPriors = kTRUE;
 373
 374   TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
 375   oadbfilename += "/PIDdefaultPriors.root";
 376   TFile * foadb = TFile::Open(oadbfilename.Data());
 377   if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 378
 379   AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
 380   if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
 381
 382   const char *namespecies[5] = {"El","Mu","Pi","Ka","Pr"};
 383   char name[100];
 384
 385   for(Int_t i=0;i < 5;i++){
 386     snprintf(name,99,"hDefault%sPriors",namespecies[i]);
 387     fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
 388     if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
 389   }
 390 }