[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAODPidHF.cxx

/**************************************************************************
 * Copyright(c) 1998-2006, 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.                  *
 * *************************************************************************/

//***********************************************************
// Class AliAODPidHF
// class for PID with AliAODRecoDecayHF
// Authors: D. Caffarri caffarri@pd.infn.it, A.Dainese andrea.dainese@pd.infn.it, S. Dash dash@to.infn.it, F. Prino prino@to.infn.it, R. Romita r.romita@gsi.de, Y. Wang yifei@pi0.physi.uni-heidelberg.de
//***********************************************************
#include "AliAODPidHF.h"
#include "AliAODPid.h"
#include "AliAODTrack.h"
#include "AliPID.h"
#include "AliTPCPIDResponse.h"
#include "AliITSPIDResponse.h"
#include "AliTOFPIDResponse.h"
#include "AliAODpidUtil.h"
#include "AliESDtrack.h"


ClassImp(AliAODPidHF)

//------------------------------
AliAODPidHF::AliAODPidHF():
  AliAODPid(),
  fnNSigma(5),
  fnSigma(),
  fTOFSigma(160.),
  fnPriors(5),
  fPriors(),
  fnPLimit(2),
  fPLimit(),
  fAsym(kFALSE)
{
 //
 // Default constructor
 //
 fPLimit=new Double_t[fnPLimit];
 fnSigma=new Double_t[fnNSigma];
 fPriors=new Double_t[fnPriors];

 for(Int_t i=0;i<fnNSigma;i++){
  fnSigma[i]=0.;
 }
 for(Int_t i=0;i<fnPriors;i++){
  fPriors[i]=0.;
 }
 for(Int_t i=0;i<fnPLimit;i++){
  fPLimit[i]=0.;
 }

}
//----------------------
AliAODPidHF::~AliAODPidHF()
{
      // destructor
 //   if(fPLimit) delete fPLimit;
 //   if(fnSigma)  delete fnSigma;
 //   if(fPriors)  delete fPriors;
}
//------------------------
AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
  AliAODPid(pid),
  fnNSigma(pid.fnNSigma),
  fnSigma(pid.fnSigma),
  fTOFSigma(pid.fTOFSigma),
  fnPriors(pid.fnPriors),
  fPriors(pid.fPriors),
  fnPLimit(pid.fnPLimit),
  fPLimit(pid.fPLimit),
  fAsym(pid.fAsym)
  {
  
  for(Int_t i=0;i<5;i++){
    fPriors[i]=pid.fPriors[i];
  }
  
  }

//----------------------
Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
// raw PID for single detectors, returns the particle type with smaller sigma
   Int_t specie=-1;
   if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
   if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
   if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);

  return specie;

}
//---------------------------
Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
// checks if the track can be a kaon, raw PID applied for single detectors
 Int_t specie=0;

 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);

 if(specie==3) return kTRUE;
 return kFALSE;
}
//---------------------------
Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
// checks if the track can be a pion, raw PID applied for single detectors

 Int_t specie=0;

 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);

 if(specie==2) return kTRUE;
 return kFALSE;
}
//---------------------------
Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
// checks if the track can be a proton raw PID applied for single detectors

 Int_t specie=0;
 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4); 
 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);

 if(specie==4) return kTRUE;

 return kFALSE;
}
//--------------------------
Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
// checks if the track can be an electron raw PID applied for single detectors

 Int_t specie=-1;
 if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
 if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
 if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);

 if(specie==0) return kTRUE;

 return kFALSE;
}
//--------------------------
Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
// n-sigma cut, TPC PID

  if(!CheckStatus(track,"TPC")) return -1;
  AliAODPid *pidObj = track->GetDetPid();
  
  Double_t dedx=pidObj->GetTPCsignal();
  Double_t mom = pidObj->GetTPCmomentum();
  AliTPCPIDResponse tpcResponse;
  Int_t pid=-1;
  if(specie<0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
   Double_t nsigmaMax=fnSigma[0];
   for(Int_t ipart=0;ipart<5;ipart++){
    AliPID::EParticleType type=AliPID::EParticleType(ipart);
    Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
    if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
     pid=ipart;
     nsigmaMax=nsigma;
    }
   }
  }else{ // asks only for one particle specie
   AliPID::EParticleType type=AliPID::EParticleType(specie);
    Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
   if (nsigma>fnSigma[0]) {
    pid=-1; 
   }else{
    pid=specie;
   }
  }

 return pid;

}
//----------------------------
Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
// truncated mean, ITS PID

  if(!CheckStatus(track,"ITS")) return -1;

  Double_t mom=track->P();
  AliAODPid *pidObj = track->GetDetPid();

  Double_t dedx=pidObj->GetITSsignal();
  AliITSPIDResponse itsResponse;
  Int_t pid=-1;
  if(specie<0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
   Double_t nsigmaMax=fnSigma[4];
   for(Int_t ipart=0;ipart<5;ipart++){
    AliPID::EParticleType type=AliPID::EParticleType(ipart);
    Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
    if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
     pid=ipart;
     nsigmaMax=nsigma;
    }
   }
  }else{ // asks only for one particle specie
   AliPID::EParticleType type=AliPID::EParticleType(specie);
    Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
   if (nsigma>fnSigma[4]) {
    pid=-1; 
   }else{
    pid=specie;
   }
  }
 return pid; 
}
//----------------------------
Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
// n-sigma cut, TOF PID

 if(!CheckStatus(track,"TOF")) return -1;

 Double_t time[AliPID::kSPECIESN];
 AliAODPid *pidObj = track->GetDetPid();
 pidObj->GetIntegratedTimes(time);
 Double_t sigTOF=pidObj->GetTOFsignal();
// AliTOFPIDResponse tofResponse;
 Int_t pid=-1;

  if(specie<0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
   Double_t nsigmaMax=fTOFSigma*fnSigma[3];
   for(Int_t ipart=0;ipart<5;ipart++){
    //AliPID::EParticleType type=AliPID::EParticleType(ipart);
    //Double_t nsigma = tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type));
    Double_t nsigma=TMath::Abs(sigTOF-time[ipart]);
    if((nsigma<nsigmaMax) && (nsigma<fnSigma[3]*fTOFSigma)) {
     pid=ipart;
     nsigmaMax=nsigma;
    }
   }
  }else{ // asks only for one particle specie
   //AliPID::EParticleType type=AliPID::EParticleType(specie);
   //Double_t nsigma = TMath::Abs(tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type)));
    Double_t nsigma=TMath::Abs(sigTOF-time[specie]);
   if (nsigma>fnSigma[3]*fTOFSigma) {
    pid=-1; 
   }else{
    pid=specie;
   }
  }
 return pid; 

}
//------------------------------
void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
// combined PID stored inside the AOD track

 const Double_t *pid=track->PID();
 Float_t max=0.;
 Int_t k=-1;
 for (Int_t i=0; i<10; i++) {
   if (pid[i]>max) {k=i; max=pid[i];}  
 }

 if(k==2) type[0]=kTRUE;
 if(k==3) type[1]=kTRUE;
 if(k==4) type[2]=kTRUE;

 return;
}
//--------------------
void AliAODPidHF::BayesianProbability(AliAODTrack *track,TString detectors,Double_t *pid) const{
// bayesian PID for single detectors or combined

  if(detectors.Contains("ITS")) {BayesianProbabilityITS(track,pid);return;}
  if(detectors.Contains("TPC")) {BayesianProbabilityTPC(track,pid);return;}
  if(detectors.Contains("TOF")) {BayesianProbabilityTOF(track,pid);return;}

  if(detectors.Contains("All")) {
    Double_t probITS[5]={0.,0.,0.,0.,0.};
    Double_t probTPC[5]={0.,0.,0.,0.,0.};
    Double_t probTOF[5]={0.,0.,0.,0.,0.};
    BayesianProbabilityITS(track,probITS);
    BayesianProbabilityTPC(track,probTPC);
    BayesianProbabilityTOF(track,probTOF);
    Double_t probTot[5]={0.,0.,0.,0.,0.};
    for(Int_t i=0;i<5;i++){
     probTot[i]=probITS[i]*probTPC[i]*probTOF[i];
    }
    for(Int_t i2=0;i2<5;i2++){
     pid[i2]=probTot[i2]*fPriors[i2]/(probTot[0]*fPriors[0]+probTot[1]*fPriors[1]+probTot[2]*fPriors[2]+probTot[3]*fPriors[3]+probTot[4]*fPriors[4]);
    }
   }

 return;

}
//------------------------------------
void AliAODPidHF::BayesianProbabilityITS(AliAODTrack *track,Double_t *prob) const{

// bayesian PID for ITS
 AliAODpidUtil pid;
 Double_t itspid[AliPID::kSPECIES];
 pid.MakeITSPID(track,itspid);
 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
  prob[ind]=itspid[ind]*fPriors[ind]/(itspid[0]*fPriors[0]+itspid[1]*fPriors[1]+itspid[2]*fPriors[2]+itspid[3]*fPriors[3]+itspid[4]*fPriors[4]);
 }
 return;

}
//------------------------------------
void AliAODPidHF::BayesianProbabilityTPC(AliAODTrack *track,Double_t *prob) const{
// bayesian PID for TPC

 AliAODpidUtil pid;
 Double_t tpcpid[AliPID::kSPECIES];
 pid.MakeTPCPID(track,tpcpid);
 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
  if(tpcpid[ind]>0.) {
   prob[ind]=tpcpid[ind]*fPriors[ind]/(tpcpid[0]*fPriors[0]+tpcpid[1]*fPriors[1]+tpcpid[2]*fPriors[2]+tpcpid[3]*fPriors[3]+tpcpid[4]*fPriors[4]);
  }else{
   prob[ind]=0.;
  }
 }
 return;

}
//------------------------------------
void AliAODPidHF::BayesianProbabilityTOF(AliAODTrack *track,Double_t *prob) const{
// bayesian PID for TOF

 AliAODpidUtil pid;
 Double_t tofpid[AliPID::kSPECIES];
 pid.MakeTOFPID(track,tofpid);
 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
  prob[ind]=tofpid[ind]*fPriors[ind]/(tofpid[0]*fPriors[0]+tofpid[1]*fPriors[1]+tofpid[2]*fPriors[2]+tofpid[3]*fPriors[3]+tofpid[4]*fPriors[4]);
 }
 return;

}
//---------------------------------
void AliAODPidHF::BayesianProbabilityTRD(AliAODTrack *track,Double_t *prob) const{
// bayesian PID for TRD

 AliAODpidUtil pid;
 Double_t trdpid[AliPID::kSPECIES];
 pid.MakeTRDPID(track,trdpid);
 for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
  if(trdpid[ind]>0.) {
   prob[ind]=trdpid[ind]*fPriors[ind]/(trdpid[0]*fPriors[0]+trdpid[1]*fPriors[1]+trdpid[2]*fPriors[2]+trdpid[3]*fPriors[3]+trdpid[4]*fPriors[4]);
  }else{
   prob[ind]=0.;
  }
 }
  return;

 }
//--------------------------------
Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{

// Quality cuts on the tracks, detector by detector

 if(detectors.Contains("ITS")){
  if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
  UChar_t clumap=track->GetITSClusterMap();
  Int_t nPointsForPid=0;
  for(Int_t i=2; i<6; i++){
   if(clumap&(1<<i)) ++nPointsForPid;
  }
  if(nPointsForPid<3) return kFALSE;// track not to be used for PID purposes
 }

 if(detectors.Contains("TPC")){
   if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
   UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
   if (nTPCClus<70) return kFALSE;
 }

 if(detectors.Contains("TOF")){
   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return kFALSE;
   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return kFALSE;
   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return kFALSE;
 }


 if(detectors.Contains("TRD")){
  if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return kFALSE;
  AliAODPid *pidObj = track->GetDetPid();
  Float_t *mom=pidObj->GetTRDmomentum();
  Int_t ntracklets=0;
  for(Int_t iPl=0;iPl<6;iPl++){
   if(mom[iPl]>0.) ntracklets++;
  }
   if(ntracklets<4) return kFALSE;
 }

 return kTRUE;
}
//--------------------------------------------
Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
// TPC nsigma cut PID, different sigmas in different p bins

  if(!CheckStatus(track,"TPC")) return kFALSE;
  AliAODPid *pidObj = track->GetDetPid();
  Double_t mom = pidObj->GetTPCmomentum();
  Double_t dedx=pidObj->GetTPCsignal();
  
  AliTPCPIDResponse tpcResponse;
  AliPID::EParticleType type=AliPID::EParticleType(specie);
  Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type)); 

  if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
  if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
  if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;

 return kFALSE;
}
//------------------
Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track,Int_t mode,Int_t specie,Bool_t compat){
// combination of the PID info coming from TPC and TOF
 if(mode==1){
  //TOF || TPC (a la' Andrea R.)
 // convention (temporary): 
 // for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
 // the method returns the sum of the response of the 2 detectors
  if(!CheckStatus(track,"TPC") && !CheckStatus(track,"TOF")) return 0;

  Int_t TPCinfo=0;
  if(CheckStatus(track,"TPC")) {
   if(fAsym) {
    if(TPCRawAsym(track,specie)) TPCinfo=1;
   }else{
    if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=1;
    if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=1;
    if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=1;
   }


  if(compat && TPCinfo<0){
   Double_t sig0tmp=fnSigma[0];
   SetSigma(0,3.);
   if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=0;
   if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=0;
   if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=0;
   SetSigma(0,sig0tmp);
  }

 }

 if(!CheckStatus(track,"TOF")) return TPCinfo;

 Int_t TOFinfo=-1;
 
  if(specie==2 && IsPionRaw(track,"TOF")) TOFinfo=1;
  if(specie==3 && IsKaonRaw(track,"TOF")) TOFinfo=1;
  if(specie==4 && IsProtonRaw(track,"TOF")) TOFinfo=1;

 if(compat && TOFinfo>0){
  Double_t ptrack=track->P();
  if(ptrack>1.5) TOFinfo=0;
 }

 return TPCinfo+TOFinfo;
 }
 if(mode==2){
  //TPC & TOF (a la' Yifei)
 // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
  Int_t TPCinfo=1; 
  if(CheckStatus(track,"TPC")) {
   if(fAsym){
    if(!TPCRawAsym(track,specie)) TPCinfo=-1;
   }else{
    if(specie==2 && !IsPionRaw(track,"TPC")) TPCinfo=-1;
    if(specie==3 && !IsKaonRaw(track,"TPC")) TPCinfo=-1;
    if(specie==4 && !IsProtonRaw(track,"TPC")) TPCinfo=-1;
   }
  }

  Int_t TOFinfo=1;
  if(!CheckStatus(track,"TOF")) return TPCinfo;

   if(specie==2 && !IsPionRaw(track,"TOF")) TOFinfo=-1;
   if(specie==3 && !IsKaonRaw(track,"TOF")) TOFinfo=-1;
   if(specie==4 && !IsProtonRaw(track,"TOF")) TOFinfo=-1;

   if(TOFinfo==1 && TPCinfo==1) return 1;
   return -1;

 }

 if(mode==3){
 //TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
 // convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
  
  if(!CheckStatus(track,"TPC") && !CheckStatus(track,"TOF")) return 0;

  Double_t ptrack=track->P();
  

  Int_t TPCinfo=-1;

   if(ptrack<fPLimit[0]) {  
    if(!CheckStatus(track,"TPC")) return 0;
    if(fAsym) {
     if(TPCRawAsym(track,specie)) TPCinfo=1;
    }else{
     if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=1;
     if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=1;
     if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=1;
    } 
    return TPCinfo;
   }

   Int_t TOFinfo=-1;
   if(ptrack>=fPLimit[0]){
    if(!CheckStatus(track,"TOF")) return 0;
    if(specie==2 && IsPionRaw(track,"TOF")) TOFinfo=1;
    if(specie==3 && IsKaonRaw(track,"TOF")) TOFinfo=1;
    if(specie==4 && IsProtonRaw(track,"TOF")) TOFinfo=1;
    return TOFinfo;
   }

 }

 return -1;
}
Commit	Line	Data
7ce8802c	1	/**************************************************************************
	2	* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	*
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	* *************************************************************************/
	16
	17	//***********************************************************
	18	// Class AliAODPidHF
	19	// class for PID with AliAODRecoDecayHF
228f3aba	20	// Authors: D. Caffarri caffarri@pd.infn.it, A.Dainese andrea.dainese@pd.infn.it, S. Dash dash@to.infn.it, F. Prino prino@to.infn.it, R. Romita r.romita@gsi.de, Y. Wang yifei@pi0.physi.uni-heidelberg.de
7ce8802c	21	//***********************************************************
	22	#include "AliAODPidHF.h"
	23	#include "AliAODPid.h"
	24	#include "AliAODTrack.h"
	25	#include "AliPID.h"
	26	#include "AliTPCPIDResponse.h"
	27	#include "AliITSPIDResponse.h"
	28	#include "AliTOFPIDResponse.h"
	29	#include "AliAODpidUtil.h"
f9e3cf8e	30	#include "AliESDtrack.h"
7ce8802c	31
	32
	33	ClassImp(AliAODPidHF)
	34
	35	//------------------------------
	36	AliAODPidHF::AliAODPidHF():
228f3aba	37	AliAODPid(),
	38	fnNSigma(5),
	39	fnSigma(),
	40	fTOFSigma(160.),
	41	fnPriors(5),
	42	fPriors(),
	43	fnPLimit(2),
	44	fPLimit(),
	45	fAsym(kFALSE)
7ce8802c	46	{
	47	//
	48	// Default constructor
	49	//
228f3aba	50	fPLimit=new Double_t[fnPLimit];
	51	fnSigma=new Double_t[fnNSigma];
	52	fPriors=new Double_t[fnPriors];
	53
	54	for(Int_t i=0;i<fnNSigma;i++){
	55	fnSigma[i]=0.;
	56	}
	57	for(Int_t i=0;i<fnPriors;i++){
	58	fPriors[i]=0.;
	59	}
	60	for(Int_t i=0;i<fnPLimit;i++){
	61	fPLimit[i]=0.;
	62	}
7ce8802c	63
	64	}
	65	//----------------------
	66	AliAODPidHF::~AliAODPidHF()
	67	{
	68	// destructor
228f3aba	69	// if(fPLimit) delete fPLimit;
	70	// if(fnSigma) delete fnSigma;
	71	// if(fPriors) delete fPriors;
7ce8802c	72	}
	73	//------------------------
	74	AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
	75	AliAODPid(pid),
228f3aba	76	fnNSigma(pid.fnNSigma),
	77	fnSigma(pid.fnSigma),
	78	fTOFSigma(pid.fTOFSigma),
	79	fnPriors(pid.fnPriors),
	80	fPriors(pid.fPriors),
	81	fnPLimit(pid.fnPLimit),
	82	fPLimit(pid.fPLimit),
	83	fAsym(pid.fAsym)
7ce8802c	84	{
	85
	86	for(Int_t i=0;i<5;i++){
	87	fPriors[i]=pid.fPriors[i];
	88	}
	89
	90	}
	91
	92	//----------------------
228f3aba	93	Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector) const{
228f3aba	94	// raw PID for single detectors, returns the particle type with smaller sigma
f9e3cf8e	95	Int_t specie=-1;
	96	if(detector.Contains("ITS")) return ApplyPidITSRaw(track,specie);
	97	if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,specie);
	98	if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,specie);
7ce8802c	99
f9e3cf8e	100	return specie;
7ce8802c	101
	102	}
	103	//---------------------------
228f3aba	104	Bool_t AliAODPidHF::IsKaonRaw(AliAODTrack *track, TString detector) const{
228f3aba	105	// checks if the track can be a kaon, raw PID applied for single detectors
7ce8802c	106	Int_t specie=0;
	107
	108	if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
	109	if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
	110	if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
	111
	112	if(specie==3) return kTRUE;
	113	return kFALSE;
	114	}
	115	//---------------------------
228f3aba	116	Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector) const{
228f3aba	117	// checks if the track can be a pion, raw PID applied for single detectors
7ce8802c	118
	119	Int_t specie=0;
	120
	121	if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
	122	if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
	123	if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
	124
	125	if(specie==2) return kTRUE;
	126	return kFALSE;
	127	}
	128	//---------------------------
228f3aba	129	Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector) const{
228f3aba	130	// checks if the track can be a proton raw PID applied for single detectors
7ce8802c	131
	132	Int_t specie=0;
	133	if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
	134	if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4);
	135	if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
	136
	137	if(specie==4) return kTRUE;
	138
	139	return kFALSE;
	140	}
f9e3cf8e	141	//--------------------------
228f3aba	142	Bool_t AliAODPidHF::IsElectronRaw(AliAODTrack *track, TString detector) const{
228f3aba	143	// checks if the track can be an electron raw PID applied for single detectors
f9e3cf8e	144
	145	Int_t specie=-1;
	146	if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,0);
	147	if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,0);
	148	if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,0);
	149
	150	if(specie==0) return kTRUE;
	151
	152	return kFALSE;
	153	}
	154	//--------------------------
228f3aba	155	Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
228f3aba	156	// n-sigma cut, TPC PID
7ce8802c	157
f9e3cf8e	158	if(!CheckStatus(track,"TPC")) return -1;
7ce8802c	159	AliAODPid *pidObj = track->GetDetPid();
	160
	161	Double_t dedx=pidObj->GetTPCsignal();
	162	Double_t mom = pidObj->GetTPCmomentum();
	163	AliTPCPIDResponse tpcResponse;
f9e3cf8e	164	Int_t pid=-1;
f9e3cf8e	165	if(specie<0){ // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
228f3aba	166	Double_t nsigmaMax=fnSigma[0];
f9e3cf8e	167	for(Int_t ipart=0;ipart<5;ipart++){
7ce8802c	168	AliPID::EParticleType type=AliPID::EParticleType(ipart);
f9e3cf8e	169	Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
228f3aba	170	if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
7ce8802c	171	pid=ipart;
	172	nsigmaMax=nsigma;
	173	}
	174	}
	175	}else{ // asks only for one particle specie
	176	AliPID::EParticleType type=AliPID::EParticleType(specie);
f9e3cf8e	177	Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
228f3aba	178	if (nsigma>fnSigma[0]) {
f9e3cf8e	179	pid=-1;
7ce8802c	180	}else{
	181	pid=specie;
	182	}
	183	}
	184
	185	return pid;
	186
	187	}
	188	//----------------------------
228f3aba	189	Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
228f3aba	190	// truncated mean, ITS PID
f9e3cf8e	191
	192	if(!CheckStatus(track,"ITS")) return -1;
	193
7ce8802c	194	Double_t mom=track->P();
	195	AliAODPid *pidObj = track->GetDetPid();
	196
	197	Double_t dedx=pidObj->GetITSsignal();
	198	AliITSPIDResponse itsResponse;
f9e3cf8e	199	Int_t pid=-1;
f9e3cf8e	200	if(specie<0){ // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
228f3aba	201	Double_t nsigmaMax=fnSigma[4];
f9e3cf8e	202	for(Int_t ipart=0;ipart<5;ipart++){
7ce8802c	203	AliPID::EParticleType type=AliPID::EParticleType(ipart);
7ce8802c	204	Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
228f3aba	205	if((nsigma<nsigmaMax) && (nsigma<fnSigma[4])) {
7ce8802c	206	pid=ipart;
	207	nsigmaMax=nsigma;
	208	}
	209	}
	210	}else{ // asks only for one particle specie
	211	AliPID::EParticleType type=AliPID::EParticleType(specie);
	212	Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
228f3aba	213	if (nsigma>fnSigma[4]) {
f9e3cf8e	214	pid=-1;
7ce8802c	215	}else{
	216	pid=specie;
	217	}
	218	}
	219	return pid;
	220	}
	221	//----------------------------
228f3aba	222	Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
228f3aba	223	// n-sigma cut, TOF PID
f9e3cf8e	224
	225	if(!CheckStatus(track,"TOF")) return -1;
	226
7ce8802c	227	Double_t time[AliPID::kSPECIESN];
	228	AliAODPid *pidObj = track->GetDetPid();
	229	pidObj->GetIntegratedTimes(time);
228f3aba	230	Double_t sigTOF=pidObj->GetTOFsignal();
228f3aba	231	// AliTOFPIDResponse tofResponse;
f9e3cf8e	232	Int_t pid=-1;
7ce8802c	233
f9e3cf8e	234	if(specie<0){ // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
228f3aba	235	Double_t nsigmaMax=fTOFSigma*fnSigma[3];
f9e3cf8e	236	for(Int_t ipart=0;ipart<5;ipart++){
228f3aba	237	//AliPID::EParticleType type=AliPID::EParticleType(ipart);
	238	//Double_t nsigma = tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type));
	239	Double_t nsigma=TMath::Abs(sigTOF-time[ipart]);
	240	if((nsigma<nsigmaMax) && (nsigma<fnSigma[3]*fTOFSigma)) {
7ce8802c	241	pid=ipart;
	242	nsigmaMax=nsigma;
	243	}
	244	}
	245	}else{ // asks only for one particle specie
228f3aba	246	//AliPID::EParticleType type=AliPID::EParticleType(specie);
	247	//Double_t nsigma = TMath::Abs(tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type)));
	248	Double_t nsigma=TMath::Abs(sigTOF-time[specie]);
	249	if (nsigma>fnSigma[3]*fTOFSigma) {
f9e3cf8e	250	pid=-1;
7ce8802c	251	}else{
	252	pid=specie;
	253	}
	254	}
	255	return pid;
	256
	257	}
	258	//------------------------------
228f3aba	259	void AliAODPidHF::CombinedProbability(AliAODTrack track,Bool_t type) const{
228f3aba	260	// combined PID stored inside the AOD track
7ce8802c	261
	262	const Double_t *pid=track->PID();
	263	Float_t max=0.;
	264	Int_t k=-1;
	265	for (Int_t i=0; i<10; i++) {
	266	if (pid[i]>max) {k=i; max=pid[i];}
	267	}
	268
	269	if(k==2) type[0]=kTRUE;
	270	if(k==3) type[1]=kTRUE;
	271	if(k==4) type[2]=kTRUE;
	272
	273	return;
	274	}
	275	//--------------------
228f3aba	276	void AliAODPidHF::BayesianProbability(AliAODTrack track,TString detectors,Double_t pid) const{
228f3aba	277	// bayesian PID for single detectors or combined
7ce8802c	278
	279	if(detectors.Contains("ITS")) {BayesianProbabilityITS(track,pid);return;}
	280	if(detectors.Contains("TPC")) {BayesianProbabilityTPC(track,pid);return;}
	281	if(detectors.Contains("TOF")) {BayesianProbabilityTOF(track,pid);return;}
	282
	283	if(detectors.Contains("All")) {
	284	Double_t probITS[5]={0.,0.,0.,0.,0.};
	285	Double_t probTPC[5]={0.,0.,0.,0.,0.};
	286	Double_t probTOF[5]={0.,0.,0.,0.,0.};
	287	BayesianProbabilityITS(track,probITS);
	288	BayesianProbabilityTPC(track,probTPC);
	289	BayesianProbabilityTOF(track,probTOF);
	290	Double_t probTot[5]={0.,0.,0.,0.,0.};
	291	for(Int_t i=0;i<5;i++){
	292	probTot[i]=probITS[i]probTPC[i]probTOF[i];
	293	}
	294	for(Int_t i2=0;i2<5;i2++){
	295	pid[i2]=probTot[i2]fPriors[i2]/(probTot[0]fPriors[0]+probTot[1]fPriors[1]+probTot[2]fPriors[2]+probTot[3]fPriors[3]+probTot[4]fPriors[4]);
	296	}
	297	}
	298
	299	return;
	300
	301	}
	302	//------------------------------------
228f3aba	303	void AliAODPidHF::BayesianProbabilityITS(AliAODTrack track,Double_t prob) const{
7ce8802c	304
228f3aba	305	// bayesian PID for ITS
7ce8802c	306	AliAODpidUtil pid;
	307	Double_t itspid[AliPID::kSPECIES];
	308	pid.MakeITSPID(track,itspid);
	309	for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
	310	prob[ind]=itspid[ind]fPriors[ind]/(itspid[0]fPriors[0]+itspid[1]fPriors[1]+itspid[2]fPriors[2]+itspid[3]fPriors[3]+itspid[4]fPriors[4]);
	311	}
	312	return;
	313
	314	}
	315	//------------------------------------
228f3aba	316	void AliAODPidHF::BayesianProbabilityTPC(AliAODTrack track,Double_t prob) const{
228f3aba	317	// bayesian PID for TPC
7ce8802c	318
	319	AliAODpidUtil pid;
	320	Double_t tpcpid[AliPID::kSPECIES];
	321	pid.MakeTPCPID(track,tpcpid);
	322	for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
	323	if(tpcpid[ind]>0.) {
	324	prob[ind]=tpcpid[ind]fPriors[ind]/(tpcpid[0]fPriors[0]+tpcpid[1]fPriors[1]+tpcpid[2]fPriors[2]+tpcpid[3]fPriors[3]+tpcpid[4]fPriors[4]);
	325	}else{
	326	prob[ind]=0.;
	327	}
	328	}
	329	return;
	330
	331	}
	332	//------------------------------------
228f3aba	333	void AliAODPidHF::BayesianProbabilityTOF(AliAODTrack track,Double_t prob) const{
228f3aba	334	// bayesian PID for TOF
7ce8802c	335
	336	AliAODpidUtil pid;
	337	Double_t tofpid[AliPID::kSPECIES];
228f3aba	338	pid.MakeTOFPID(track,tofpid);
7ce8802c	339	for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
	340	prob[ind]=tofpid[ind]fPriors[ind]/(tofpid[0]fPriors[0]+tofpid[1]fPriors[1]+tofpid[2]fPriors[2]+tofpid[3]fPriors[3]+tofpid[4]fPriors[4]);
	341	}
	342	return;
	343
	344	}
f9e3cf8e	345	//---------------------------------
228f3aba	346	void AliAODPidHF::BayesianProbabilityTRD(AliAODTrack track,Double_t prob) const{
228f3aba	347	// bayesian PID for TRD
f9e3cf8e	348
	349	AliAODpidUtil pid;
	350	Double_t trdpid[AliPID::kSPECIES];
	351	pid.MakeTRDPID(track,trdpid);
	352	for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
	353	if(trdpid[ind]>0.) {
	354	prob[ind]=trdpid[ind]fPriors[ind]/(trdpid[0]fPriors[0]+trdpid[1]fPriors[1]+trdpid[2]fPriors[2]+trdpid[3]fPriors[3]+trdpid[4]fPriors[4]);
	355	}else{
	356	prob[ind]=0.;
	357	}
	358	}
	359	return;
	360
	361	}
	362	//--------------------------------
228f3aba	363	Bool_t AliAODPidHF::CheckStatus(AliAODTrack *track,TString detectors) const{
f9e3cf8e	364
228f3aba	365	// Quality cuts on the tracks, detector by detector
f9e3cf8e	366
	367	if(detectors.Contains("ITS")){
	368	if ((track->GetStatus()&AliESDtrack::kITSin)==0) return kFALSE;
	369	UChar_t clumap=track->GetITSClusterMap();
	370	Int_t nPointsForPid=0;
	371	for(Int_t i=2; i<6; i++){
	372	if(clumap&(1<<i)) ++nPointsForPid;
	373	}
	374	if(nPointsForPid<3) return kFALSE;// track not to be used for PID purposes
	375	}
	376
	377	if(detectors.Contains("TPC")){
	378	if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return kFALSE;
	379	UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
	380	if (nTPCClus<70) return kFALSE;
	381	}
	382
	383	if(detectors.Contains("TOF")){
	384	if ((track->GetStatus()&AliESDtrack::kTOFout )==0) return kFALSE;
	385	if ((track->GetStatus()&AliESDtrack::kTIME )==0) return kFALSE;
	386	if ((track->GetStatus()&AliESDtrack::kTOFpid )==0) return kFALSE;
	387	}
	388
	389
	390	if(detectors.Contains("TRD")){
	391	if ((track->GetStatus()&AliESDtrack::kTRDout )==0) return kFALSE;
	392	AliAODPid *pidObj = track->GetDetPid();
	393	Float_t *mom=pidObj->GetTRDmomentum();
	394	Int_t ntracklets=0;
	395	for(Int_t iPl=0;iPl<6;iPl++){
	396	if(mom[iPl]>0.) ntracklets++;
	397	}
	398	if(ntracklets<4) return kFALSE;
	399	}
	400
	401	return kTRUE;
	402	}
228f3aba	403	//--------------------------------------------
	404	Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
	405	// TPC nsigma cut PID, different sigmas in different p bins
	406
	407	if(!CheckStatus(track,"TPC")) return kFALSE;
	408	AliAODPid *pidObj = track->GetDetPid();
	409	Double_t mom = pidObj->GetTPCmomentum();
	410	Double_t dedx=pidObj->GetTPCsignal();
	411
	412	AliTPCPIDResponse tpcResponse;
	413	AliPID::EParticleType type=AliPID::EParticleType(specie);
	414	Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,track->GetTPCNcls(),type));
	415
	416	if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
	417	if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
	418	if(mom>fPLimit[1] && nsigma<fnSigma[2]) return kTRUE;
	419
	420	return kFALSE;
	421	}
	422	//------------------
	423	Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track,Int_t mode,Int_t specie,Bool_t compat){
	424	// combination of the PID info coming from TPC and TOF
	425	if(mode==1){
	426	//TOF \|\| TPC (a la' Andrea R.)
	427	// convention (temporary):
	428	// for the single detectors: -1 = kFALSE, 1 = kTRUE, 0 = compatible
	429	// the method returns the sum of the response of the 2 detectors
	430	if(!CheckStatus(track,"TPC") && !CheckStatus(track,"TOF")) return 0;
	431
	432	Int_t TPCinfo=0;
	433	if(CheckStatus(track,"TPC")) {
	434	if(fAsym) {
	435	if(TPCRawAsym(track,specie)) TPCinfo=1;
	436	}else{
	437	if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=1;
	438	if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=1;
	439	if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=1;
	440	}
	441
	442
	443	if(compat && TPCinfo<0){
53c43440	444	Double_t sig0tmp=fnSigma[0];
228f3aba	445	SetSigma(0,3.);
	446	if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=0;
	447	if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=0;
	448	if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=0;
53c43440	449	SetSigma(0,sig0tmp);
228f3aba	450	}
	451
	452	}
	453
	454	if(!CheckStatus(track,"TOF")) return TPCinfo;
	455
	456	Int_t TOFinfo=-1;
	457
	458	if(specie==2 && IsPionRaw(track,"TOF")) TOFinfo=1;
	459	if(specie==3 && IsKaonRaw(track,"TOF")) TOFinfo=1;
	460	if(specie==4 && IsProtonRaw(track,"TOF")) TOFinfo=1;
	461
	462	if(compat && TOFinfo>0){
	463	Double_t ptrack=track->P();
	464	if(ptrack>1.5) TOFinfo=0;
	465	}
	466
	467	return TPCinfo+TOFinfo;
	468	}
	469	if(mode==2){
	470	//TPC & TOF (a la' Yifei)
	471	// convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
	472	Int_t TPCinfo=1;
	473	if(CheckStatus(track,"TPC")) {
	474	if(fAsym){
	475	if(!TPCRawAsym(track,specie)) TPCinfo=-1;
	476	}else{
	477	if(specie==2 && !IsPionRaw(track,"TPC")) TPCinfo=-1;
	478	if(specie==3 && !IsKaonRaw(track,"TPC")) TPCinfo=-1;
	479	if(specie==4 && !IsProtonRaw(track,"TPC")) TPCinfo=-1;
	480	}
	481	}
	482
	483	Int_t TOFinfo=1;
	484	if(!CheckStatus(track,"TOF")) return TPCinfo;
	485
	486	if(specie==2 && !IsPionRaw(track,"TOF")) TOFinfo=-1;
	487	if(specie==3 && !IsKaonRaw(track,"TOF")) TOFinfo=-1;
	488	if(specie==4 && !IsProtonRaw(track,"TOF")) TOFinfo=-1;
	489
	490	if(TOFinfo==1 && TPCinfo==1) return 1;
	491	return -1;
	492
	493	}
	494
	495	if(mode==3){
	496	//TPC for p<fPLimit[0], TOF for p>=fPLimit[0] (a la' Andrea A.)
	497	// convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
	498
	499	if(!CheckStatus(track,"TPC") && !CheckStatus(track,"TOF")) return 0;
	500
	501	Double_t ptrack=track->P();
	502
	503
	504	Int_t TPCinfo=-1;
	505
	506	if(ptrack<fPLimit[0]) {
	507	if(!CheckStatus(track,"TPC")) return 0;
	508	if(fAsym) {
	509	if(TPCRawAsym(track,specie)) TPCinfo=1;
	510	}else{
	511	if(specie==2 && IsPionRaw(track,"TPC")) TPCinfo=1;
	512	if(specie==3 && IsKaonRaw(track,"TPC")) TPCinfo=1;
	513	if(specie==4 && IsProtonRaw(track,"TPC")) TPCinfo=1;
514	}
515	return TPCinfo;
516	}
517
518	Int_t TOFinfo=-1;
519	if(ptrack>=fPLimit[0]){
520	if(!CheckStatus(track,"TOF")) return 0;
521	if(specie==2 && IsPionRaw(track,"TOF")) TOFinfo=1;
522	if(specie==3 && IsKaonRaw(track,"TOF")) TOFinfo=1;
523	if(specie==4 && IsProtonRaw(track,"TOF")) TOFinfo=1;
524	return TOFinfo;
525	}
526
527	}
528
529	return -1;
530	}
531