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