[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDCombined.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


//-----------------------------------------------------------------
//        Base class for combining PID of different detectors    //
//        (user selected) and compute Bayesian probabilities     //
//                                                               //
//                                                               //
//   Origin: Pietro Antonioli, INFN-BO Pietro.Antonioli@cern.ch  //
//                                                               //
//-----------------------------------------------------------------

#include <TH1.h>

#include <AliVTrack.h>
#include <AliLog.h>
#include <AliPID.h>
#include <AliPIDResponse.h>

#include "AliPIDCombined.h"

#include "TMath.h"
#include "TFile.h"

#include "AliOADBContainer.h"

// initialize static members
TH2F* AliPIDCombined::fDefaultPriorsTPC[]={0x0};

ClassImp(AliPIDCombined);

AliPIDCombined::AliPIDCombined():
	TNamed("CombinedPID","CombinedPID"),
	fDetectorMask(0),
	fRejectMismatchMask(0x7F),
	fEnablePriors(kTRUE),
	fSelectedSpecies(AliPID::kSPECIES),
	fUseDefaultTPCPriors(kFALSE)
{	
  //
  // default constructor
  //	
  for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
  AliLog::SetClassDebugLevel("AliPIDCombined",10);
}

//-------------------------------------------------------------------------------------------------	
AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
	TNamed(name,title),
	fDetectorMask(0),
	fRejectMismatchMask(0x7F),
	fEnablePriors(kTRUE),
	fSelectedSpecies(AliPID::kSPECIES),
	fUseDefaultTPCPriors(kFALSE)
{
  //
  // default constructor with name and title
  //
  for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
  AliLog::SetClassDebugLevel("AliPIDCombined",10);

}

//-------------------------------------------------------------------------------------------------	
AliPIDCombined::~AliPIDCombined() {

  for(Int_t i=0;i < AliPID::kSPECIESC;i++){
    if(fPriorsDistributions[i])
      delete fPriorsDistributions[i];
  }
}

//-------------------------------------------------------------------------------------------------	
void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
  if ( (type < 0) || ( type >= ((AliPID::EParticleType)AliPID::kSPECIESC) ) ){
    AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
    return;
  }
  if(prior) {
    Int_t i = (Int_t)type;
    if (fPriorsDistributions[i]) {
      delete fPriorsDistributions[i]; 
    }
    fPriorsDistributions[i]=new TH1F(*prior);
  }
}

//-------------------------------------------------------------------------------------------------	
UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const {
  //
  // (1) Get raw probabilities of requested detectors and combine
  // (2) Get priors and propagate depending on detectors used
  // (3) Compute Bayes probabilities
  //


  // (1) Get raw probabilities of selected detectors and combine
  UInt_t usedMask=0;             // detectors actually used for track
  AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
  Double_t p[fSelectedSpecies];  // combined probabilities of selected detectors
  for (Int_t i=0;i<fSelectedSpecies;i++) p[i]=1.; // no decision
  for (Int_t ibit = 0; ibit < 7 ; ibit++) {
    AliPIDResponse::EDetCode detBit = (AliPIDResponse::EDetCode)(1<<ibit);
    if (fDetectorMask & detBit) {  	    // getting probabilities for requested detectors only
      Double_t detProb[fSelectedSpecies];
      status = response->ComputePIDProbability(detBit,track,fSelectedSpecies,detProb);
      if (status == AliPIDResponse::kDetPidOk) {
	if (fRejectMismatchMask & detBit) {         // mismatch check (currently just for TOF)
	  if (detBit == AliPIDResponse::kDetTOF) {
	    Float_t probMis = response->GetTOFMismatchProbability(track);
	    SetCombinedStatus(status,&usedMask,detBit,detProb,probMis);
	  } else {
	    SetCombinedStatus(status,&usedMask,detBit);
	  }
	} else {
	  SetCombinedStatus(status,&usedMask,detBit);
	}
	for (Int_t i=0;i<fSelectedSpecies;i++) p[i] *= detProb[i];
      }
    }
  }
  // if no detectors available there is no point to go further
  if (usedMask == 0) return usedMask;

  // (2) Get priors and propagate depending on detectors used
  Double_t priors[fSelectedSpecies];
  memset(priors,0,fSelectedSpecies*sizeof(Double_t));
  if (fEnablePriors){
    GetPriors(track,priors,response->GetCurrentCentrality());
    
    // for the moment we have three cases
    // TPC+EMCAL    --> apply EMCAL propagation factors (TRD and TOF may be present)
    // TPC+TOF      --> apply TOF propagation factors (TRD may be present)
    // TPC+TRD      --> apply TRD propagation factors
    if(fUseDefaultTPCPriors) {
      Double_t pt=TMath::Abs(track->Pt());
      if ( (usedMask & AliPIDResponse::kDetEMCAL)==AliPIDResponse::kDetEMCAL ) { // EMCAL is the outer having prop. factors for the moment
	// EMCal case (for the moment only in combination with TPC)
	// propagation factors determined from LHC11d MC (LHC12a15f)
	// v2 clusterizer, dEta < 0.015, dPhi < 0.03, NonLinearityFunction = 6
	       
	Float_t electronEMCALfactor=0.1;
	Float_t kaonEMCALfactor = 0.1;
	Float_t protonEMCALfactor = 0.1;
	       
	if(track->Charge() > 0){
	  // positiv charge (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
	  if(pt > 0.75 && pt < 20.0){
	    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)) );
	    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)) );
	    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)) );
		   
	  }
	}
	       
	if(track->Charge() < 0){
	  // negative charge  (start parametrization at 0.75 GeV/c and stop at 20 GeV/c for the moment)
	  if(pt > 0.75 && pt < 20.0){ 		   
	    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)) );
	    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)) );
	    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)) );
	  }
	}
	priors[0] *= electronEMCALfactor;
	priors[3] *= kaonEMCALfactor;
	priors[4] *= protonEMCALfactor;	      
      } // end of EMCAL case
      else if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){
	Float_t kaonTOFfactor = 0.1;
	if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
	Float_t protonTOFfactor = 0.1;
	if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
	       
	if(track->Charge() < 0){ // for negative tracks
	  kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
	  protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
	}
	// TODO: we may need an electron factor for TOF as well, especially if TRD is present!
	priors[3] *= kaonTOFfactor;
	priors[4] *= protonTOFfactor;
      }  // end of TOF case
      else if ( (usedMask & AliPIDResponse::kDetTRD)==AliPIDResponse::kDetTRD ) {
	Float_t electronTRDfactor=0.1;
	Float_t kaonTRDfactor = 0.1;
	Float_t protonTRDfactor = 0.1;
		 
	if(track->Charge() > 0){
	  // positiv charge
	  if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,5.13315e-03)/2.11145e-01)*TMath::Power(pt,-2.97659e+00);
	  if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,-4.29549e-02)/4.87989e-01)*TMath::Power(pt,-1.54270e+00);
	  if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,2.81238e+00)/7.57082e-02)*TMath::Power(pt,-8.12595e-01);
	}
		 
	if(track->Charge() < 0){
	  // negative charge
	  if(pt > 0.25) electronTRDfactor =  1 - TMath::Exp(-TMath::Power(pt,2.45537e-02)/1.90397e-01)*TMath::Power(pt,-3.33121e+00);
	  if(pt > 0.35) kaonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt, -3.42831e-03)/5.57013e-01)*TMath::Power(pt,-1.39202e+00);
	  if(pt > 0.35) protonTRDfactor = 1 - TMath::Exp(-TMath::Power(pt,3.36631e+00)/7.18819e-02)*TMath::Power(pt,-2.00577e-01);
	}
	// what about electrons
	priors[0] *= electronTRDfactor;
	priors[3] *= kaonTRDfactor;
	priors[4] *= protonTRDfactor;	      
      } // end of TRD case
    } // end of fUseDefaultTPCPriors
  }   // end of use priors
  else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}

  // (3) Compute Bayes probabilities
  ComputeBayesProbabilities(bayesProbabilities,p,priors);
  return usedMask;
}


//-------------------------------------------------------------------------------------------------
void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
	
	//
	// get priors from distributions
	//
	
	Double_t pt=TMath::Abs(track->Pt());

        if(fUseDefaultTPCPriors){ // use default priors if requested
	  Float_t usedCentr = centrality+5*(centrality>0);
	  if(usedCentr < -0.99) usedCentr = -0.99;
	  else if(usedCentr > 99) usedCentr = 99;
	  if(pt > 9.99) pt = 9.99;
	  else if(pt < 0.01)  pt = 0.01;
	  
	  for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);

	  return;
	}
	
	Double_t sumPriors = 0;
	for (Int_t i=0;i<fSelectedSpecies;++i){
		if (fPriorsDistributions[i]){
			p[i]=fPriorsDistributions[i]->Interpolate(pt);
		}
		else {
			p[i]=0.;
		}
		sumPriors+=p[i];		
	}

	// normalizing priors
	if (sumPriors == 0) return;
	for (Int_t i=0;i<fSelectedSpecies;++i){
	   p[i]=p[i]/sumPriors;
	}
	return;
}
//-------------------------------------------------------------------------------------------------
void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
	
	//
	// get priors from distributions
	//
        Double_t pt=TMath::Abs(track->Pt());
	
        if(fUseDefaultTPCPriors){ // use default priors if requested
	  Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
	  if(usedCentr < -0.99) usedCentr = -0.99;
	  else if(usedCentr > 99) usedCentr = 99;
	  if(pt > 9.99) pt = 9.99;
	  else if(pt < 0.01)  pt = 0.01;
	  
	  for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);

	  // Extra factor if TOF matching was required
	  if(detUsed & AliPIDResponse::kDetTOF){
	    Float_t kaonTOFfactor = 0.1;
	    if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
	    Float_t protonTOFfactor = 0.1;
	    if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
	    
	    if(track->Charge() < 0){ // for negative tracks
	      kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
	      protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
	    }
	    
	    p[3] *= kaonTOFfactor;
	    p[4] *= protonTOFfactor;
	  }

	  return;
	}


	Double_t sumPriors = 0;
	for (Int_t i=0;i<fSelectedSpecies;++i){
		if (fPriorsDistributions[i]){
			p[i]=fPriorsDistributions[i]->Interpolate(pt);
		}
		else {
			p[i]=0.;
		}
		sumPriors+=p[i];		
	}

	// normalizing priors
	if (sumPriors == 0) return;
	for (Int_t i=0;i<fSelectedSpecies;++i){
	   p[i]=p[i]/sumPriors;
	}
	return;
}
//-------------------------------------------------------------------------------------------------	
void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {


  //
  // calculate Bayesian probabilities
  //
  Double_t sum = 0.;
  for (Int_t i = 0; i < fSelectedSpecies; i++) {
    sum += probDensity[i] * prior[i];
  }
  if (sum <= 0) {

    AliError("Invalid probability densities or priors");
    for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
    return;
  }
  for (Int_t i = 0; i < fSelectedSpecies; i++) {
    probabilities[i] = probDensity[i] * prior[i] / sum;
  }


}


//----------------------------------------------------------------------------------------
void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit) const {
  switch (status) {
  case AliPIDResponse::kDetNoParams:
  case AliPIDResponse::kDetNoSignal:
  case AliPIDResponse::kDetMismatch: // for backward compatibilty, we need then to remove kDetMismatch from AliPIDResponse
    break;
  case AliPIDResponse::kDetPidOk:
    *maskDetIn = *maskDetIn | bit;
    break;
  }

}

//----------------------------------------------------------------------------------------
void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit, Double_t* p, const Float_t probMis) const {
  switch (status) {
  case AliPIDResponse::kDetNoParams:
  case AliPIDResponse::kDetNoSignal:
  case AliPIDResponse::kDetMismatch: // for backward compatibility, we need then to remove kDeteMismatch from AliPIDResponse
    break;
  case AliPIDResponse::kDetPidOk:
      if (probMis > 0.5) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
      else *maskDetIn = *maskDetIn | bit;
    break;
  }

}


//----------------------------------------------------------------------------------------
void AliPIDCombined::SetDefaultTPCPriors(){
  fEnablePriors=kTRUE;
  fUseDefaultTPCPriors = kTRUE;

  //check if priors are already initialized
  if (fDefaultPriorsTPC[0]) return;
  
  TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
  oadbfilename += "/PIDdefaultPriors.root";
  TFile * foadb = TFile::Open(oadbfilename.Data());
  if(!foadb || !foadb->IsOpen()) {
    delete foadb;
    AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
    return;
  }
  
  AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
  if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
  
  const char *namespecies[AliPID::kSPECIES] = {"El","Mu","Pi","Ka","Pr"};
  char name[100];

  for(Int_t i=0;i < AliPID::kSPECIES;i++){
    snprintf(name,99,"hDefault%sPriors",namespecies[i]);
    fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
    if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
  }

  delete foadb;
}
Commit	Line	Data
b8bfee30	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
80f28562	35	#include "TMath.h"
	36	#include "TFile.h"
	37
	38	#include "AliOADBContainer.h"
	39
	40	// initialize static members
00a38d07	41	TH2F* AliPIDCombined::fDefaultPriorsTPC[]={0x0};
80f28562	42
	43	ClassImp(AliPIDCombined);
	44
b8bfee30	45	AliPIDCombined::AliPIDCombined():
	46	TNamed("CombinedPID","CombinedPID"),
	47	fDetectorMask(0),
	48	fRejectMismatchMask(0x7F),
	49	fEnablePriors(kTRUE),
80f28562	50	fSelectedSpecies(AliPID::kSPECIES),
80f28562	51	fUseDefaultTPCPriors(kFALSE)
b8bfee30	52	{
	53	//
	54	// default constructor
	55	//
00a38d07	56	for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
b8bfee30	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),
80f28562	66	fSelectedSpecies(AliPID::kSPECIES),
80f28562	67	fUseDefaultTPCPriors(kFALSE)
b8bfee30	68	{
	69	//
	70	// default constructor with name and title
	71	//
00a38d07	72	for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
b8bfee30	73	AliLog::SetClassDebugLevel("AliPIDCombined",10);
	74
	75	}
	76
	77	//-------------------------------------------------------------------------------------------------
	78	AliPIDCombined::~AliPIDCombined() {
	79
00a38d07	80	for(Int_t i=0;i < AliPID::kSPECIESC;i++){
b8bfee30	81	if(fPriorsDistributions[i])
	82	delete fPriorsDistributions[i];
	83	}
	84	}
	85
	86	//-------------------------------------------------------------------------------------------------
	87	void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
00a38d07	88	if ( (type < 0) \|\| ( type >= ((AliPID::EParticleType)AliPID::kSPECIESC) ) ){
b8bfee30	89	AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
	90	return;
	91	}
	92	if(prior) {
8de3b790	93	Int_t i = (Int_t)type;
b8bfee30	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 {
3ed61c90	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	//
b8bfee30	108
b8bfee30	109
3ed61c90	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);
b8bfee30	130	}
3ed61c90	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
44bea91d	154
3ed61c90	155	Float_t electronEMCALfactor=0.1;
	156	Float_t kaonEMCALfactor = 0.1;
	157	Float_t protonEMCALfactor = 0.1;
44bea91d	158
3ed61c90	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)) );
44bea91d	165
3ed61c90	166	}
3ed61c90	167	}
44bea91d	168
3ed61c90	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);
df0293b3	186
3ed61c90	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;
df0293b3	199
3ed61c90	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	}
df0293b3	206
3ed61c90	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;
b8bfee30	225	}
	226
	227
	228	//-------------------------------------------------------------------------------------------------
80f28562	229	void AliPIDCombined::GetPriors(const AliVTrack track, Double_t p,Float_t centrality) const {
b8bfee30	230
	231	//
	232	// get priors from distributions
	233	//
	234
	235	Double_t pt=TMath::Abs(track->Pt());
80f28562	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
b8bfee30	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	}
80f28562	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	}
b8bfee30	314
80f28562	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	}
b8bfee30	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) {
80f28562	334
b8bfee30	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
3ed61c90	346
b8bfee30	347	//----------------------------------------------------------------------------------------
3ed61c90	348	void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UInt_t *maskDetIn, AliPIDResponse::EDetCode bit) const {
b8bfee30	349	switch (status) {
3ed61c90	350	case AliPIDResponse::kDetNoParams:
	351	case AliPIDResponse::kDetNoSignal:
	352	case AliPIDResponse::kDetMismatch: // for backward compatibilty, we need then to remove kDetMismatch from AliPIDResponse
b8bfee30	353	break;
	354	case AliPIDResponse::kDetPidOk:
	355	maskDetIn = maskDetIn \| bit;
	356	break;
3ed61c90	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;
b8bfee30	371	break;
	372	}
	373
	374	}
	375
3ed61c90	376
3ed61c90	377
80f28562	378	//----------------------------------------------------------------------------------------
	379	void AliPIDCombined::SetDefaultTPCPriors(){
	380	fEnablePriors=kTRUE;
	381	fUseDefaultTPCPriors = kTRUE;
	382
00a38d07	383	//check if priors are already initialized
	384	if (fDefaultPriorsTPC[0]) return;
	385
80f28562	386	TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
	387	oadbfilename += "/PIDdefaultPriors.root";
	388	TFile * foadb = TFile::Open(oadbfilename.Data());
00a38d07	389	if(!foadb \|\| !foadb->IsOpen()) {
	390	delete foadb;
	391	AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
	392	return;
	393	}
80f28562	394
	395	AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
	396	if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
	397
00a38d07	398	const char *namespecies[AliPID::kSPECIES] = {"El","Mu","Pi","Ka","Pr"};
80f28562	399	char name[100];
80f28562	400
00a38d07	401	for(Int_t i=0;i < AliPID::kSPECIES;i++){
4c234df8	402	snprintf(name,99,"hDefault%sPriors",namespecies[i]);
80f28562	403	fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
	404	if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
	405	}
00a38d07	406
00a38d07	407	delete foadb;
80f28562	408	}