[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[5];

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::kSPECIES+AliPID::kSPECIESLN;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::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
  AliLog::SetClassDebugLevel("AliPIDCombined",10);

}

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

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

//-------------------------------------------------------------------------------------------------	
void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
  if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
    AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
    return;
  }
  if(prior) {
    Int_t i = (Int_t)type;
    if (i >= AliPID::kSPECIES ) {
      if (i < (Int_t)AliPID::kDeuteron) {
	AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
	return;
      } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
    }
    if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) {             // coverity fix (to put it mildly....)
      AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
      return;
    }
    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 {
	Double_t pITS[fSelectedSpecies];
	Double_t pTPC[fSelectedSpecies];
	Double_t pTRD[fSelectedSpecies];
	Double_t pTOF[fSelectedSpecies];
	Double_t pHMPID[fSelectedSpecies];
	Double_t pEMCAL[fSelectedSpecies];
	Double_t pPHOS[fSelectedSpecies];
	for (Int_t i=0;i<fSelectedSpecies;i++) {
	 pITS[i]=1.;
	 pTPC[i]=1.;
	 pTRD[i]=1.;
	 pTOF[i]=1.;
	 pHMPID[i]=1.;
	 pEMCAL[i]=1.;
	 pPHOS[i]=1.;
	}
	UInt_t usedMask=0;
	AliPIDResponse::EDetPidStatus status=AliPIDResponse::kDetNoSignal;
	Double_t p[fSelectedSpecies];
	memset(p,0,sizeof(Double_t)*fSelectedSpecies);

	// getting probability distributions for selected detectors only
	if (fDetectorMask & AliPIDResponse::kDetITS) {
	  status = response->ComputeITSProbability(track, fSelectedSpecies, pITS);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetITS,pITS);
	}

	if (fDetectorMask & AliPIDResponse::kDetTPC) { 
	  status = response->ComputeTPCProbability(track, fSelectedSpecies, pTPC);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTPC,pTPC);
	}


	if (fDetectorMask & AliPIDResponse::kDetTRD) { 
	  status = response->ComputeTRDProbability(track, fSelectedSpecies, pTRD);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTRD,pTRD);
	}

	if (fDetectorMask &  AliPIDResponse::kDetTOF) { 
	  status = response->ComputeTOFProbability(track, fSelectedSpecies, pTOF);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetTOF,pTOF);
	}

	if (fDetectorMask & AliPIDResponse::kDetHMPID) { 
	  status = response->ComputeHMPIDProbability(track, fSelectedSpecies, pHMPID);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetHMPID,pHMPID);
	}


	if (fDetectorMask & AliPIDResponse::kDetEMCAL) { 
	  status = response->ComputeEMCALProbability(track, fSelectedSpecies, pEMCAL);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetEMCAL,pEMCAL);
	}


	if (fDetectorMask & AliPIDResponse::kDetPHOS) { 
	  status = response->ComputePHOSProbability(track, fSelectedSpecies, pPHOS);
	  SetCombinedStatus(status,&usedMask,AliPIDResponse::kDetPHOS,pPHOS);
	}


	for (Int_t i =0; i<fSelectedSpecies; i++){
	  p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
	}
	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+TRD      --> apply TRD propagation factors
	  // TPC+TOF      --> apply TOF propagation factors
	  // TPC+TRD+TOF  --> apply TOF propagation factors
	  // apply tof matching efficiency to priors if TOF joined PID for this track
	  if(fUseDefaultTPCPriors) {
	     Double_t pt=TMath::Abs(track->Pt());
	     if ( (usedMask & AliPIDResponse::kDetTOF) == AliPIDResponse::kDetTOF ){ // TOF is the outer having prop. factors for the moment
	       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;
	     }
	     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
		 p[0] *= electronTRDfactor;
		 p[3] *= kaonTRDfactor;
		 p[4] *= protonTRDfactor;	      
	       } // end of TRD case
	     } // end of detectors inner than TOF
	  } // end of fUseDefaultTPCPriors
	}   // end of use priors
	else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
	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, Double_t* p) const {
  switch (status) {
  case AliPIDResponse::kDetNoSignal:
    break;
  case AliPIDResponse::kDetPidOk:
    *maskDetIn = *maskDetIn | bit;
    break;
  case AliPIDResponse::kDetMismatch:
    if ( fRejectMismatchMask & bit) for (Int_t j=0;j<fSelectedSpecies;j++) p[j]=1./fSelectedSpecies;
    else *maskDetIn = *maskDetIn | bit;
    break;
  }

}

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

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

  for(Int_t i=0;i < 5;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]));
  }
}
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
	41	TH2F* AliPIDCombined::fDefaultPriorsTPC[5];
	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	//
	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),
80f28562	66	fSelectedSpecies(AliPID::kSPECIES),
80f28562	67	fUseDefaultTPCPriors(kFALSE)
b8bfee30	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) {
8de3b790	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));
b8bfee30	97	return;
8de3b790	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;
b8bfee30	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];
4c234df8	132	memset(p,0,sizeof(Double_t)*fSelectedSpecies);
b8bfee30	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];
4c234df8	178	memset(priors,0,fSelectedSpecies*sizeof(Double_t));
80f28562	179	if (fEnablePriors){
	180	GetPriors(track,priors,response->GetCurrentCentrality());
	181
df0293b3	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
80f28562	186	// apply tof matching efficiency to priors if TOF joined PID for this track
df0293b3	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
b8bfee30	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	//-------------------------------------------------------------------------------------------------
80f28562	236	void AliPIDCombined::GetPriors(const AliVTrack track, Double_t p,Float_t centrality) const {
b8bfee30	237
	238	//
	239	// get priors from distributions
	240	//
	241
	242	Double_t pt=TMath::Abs(track->Pt());
80f28562	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
b8bfee30	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	}
80f28562	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	}
b8bfee30	321
80f28562	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	}
b8bfee30	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) {
80f28562	341
b8bfee30	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
80f28562	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++){
4c234df8	386	snprintf(name,99,"hDefault%sPriors",namespecies[i]);
80f28562	387	fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
	388	if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
	389	}
	390	}