// 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 P. Antonioli pietro.antonioli@bo.infn.it
//***********************************************************
#include <TCanvas.h>
+#include <TString.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TFile.h>
#include "AliAODPidHF.h"
#include "AliAODPid.h"
fppLowEn2011(kFALSE),
fPbPb(kFALSE),
fTOFdecide(kFALSE),
- fOldPid(kTRUE),
+ fOldPid(kFALSE),
fPtThresholdTPC(999999.),
+ fMaxTrackMomForCombinedPID(999999.),
fPidResponse(0),
fPidCombined(new AliPIDCombined()),
- fTPCResponse(new AliTPCPIDResponse())
+ fTPCResponse(new AliTPCPIDResponse()),
+ fPriorsH(),
+ fCombDetectors(kTPCTOF),
+ fUseCombined(kFALSE)
{
//
// Default constructor
// if(fnSigma) delete fnSigma;
// if(fPriors) delete fPriors;
delete fTPCResponse;
+ for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
+ delete fPriorsH[ispecies];
+ }
}
//------------------------
AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
fTOFdecide(pid.fTOFdecide),
fOldPid(pid.fOldPid),
fPtThresholdTPC(pid.fPtThresholdTPC),
+ fMaxTrackMomForCombinedPID(pid.fMaxTrackMomForCombinedPID),
fPidResponse(pid.fPidResponse),
fPidCombined(pid.fPidCombined),
- fTPCResponse(pid.fTPCResponse)
+ fTPCResponse(0x0),
+ fCombDetectors(pid.fCombDetectors),
+ fUseCombined(pid.fUseCombined)
{
fnSigma = new Double_t[fnNSigma];
for(Int_t i=0;i<fnPLimit;i++){
fPLimit[i]=pid.fPLimit[i];
}
+ fPriors = new Double_t[fnPriors];
+ for(Int_t i=0;i<fnPriors;i++){
+ fPriors[i]=pid.fPriors[i];
+ }
+ for(Int_t i=0;i<AliPID::kSPECIES;i++){
+ fPriorsH[i]=pid.fPriorsH[i];
+ }
-
+ if(pid.fTPCResponse) fTPCResponse = new AliTPCPIDResponse(*(pid.fTPCResponse));
//fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
//fPidCombined = new AliPIDCombined(*(pid.fPidCombined));
Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const{
// n-sigma cut, TPC PID
- if(!CheckTPCPIDStatus(track)) return 0;
+ Double_t nsigma;
Int_t pid=-1;
- if(fOldPid){
- AliAODPid *pidObj = track->GetDetPid();
-
- Double_t dedx=pidObj->GetTPCsignal();
- Double_t mom = pidObj->GetTPCmomentum();
- if(mom>fPtThresholdTPC) return 0;
- UShort_t nTPCClus=pidObj->GetTPCsignalN();
- if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
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(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
- if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
- pid=ipart;
- nsigmaMax=nsigma;
+ Double_t nsigmaMin=999.;
+ for(Int_t ipart=0;ipart<5;ipart++){
+ if(GetnSigmaTPC(track,ipart,nsigma)==1){
+ nsigma=TMath::Abs(nsigma);
+ if((nsigma<nsigmaMin) && (nsigma<fnSigma[0])) {
+ pid=ipart;
+ nsigmaMin=nsigma;
+ }
+ }
}
- }
}else{ // asks only for one particle specie
- AliPID::EParticleType type=AliPID::EParticleType(specie);
- Double_t nsigma = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
- if (nsigma>fnSigma[0]) {
- pid=-1;
- }else{
- pid=specie;
- }
- }
- }else{ //old pid
- 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(fPidResponse->NumberOfSigmasTPC(track,type));
- if((nsigma<nsigmaMax) && (nsigma<fnSigma[0])) {
- pid=ipart;
- nsigmaMax=nsigma;
+ if(GetnSigmaTPC(track,specie,nsigma)==1){
+ nsigma=TMath::Abs(nsigma);
+ if (nsigma>fnSigma[0]) pid=-1;
+ else pid=specie;
}
- }
- }else{ // asks only for one particle specie
- AliPID::EParticleType type=AliPID::EParticleType(specie);
- Double_t nsigma = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
- if (nsigma>fnSigma[0]) {
- pid=-1;
- }else{
- pid=specie;
- }
}
- } //new pid
-
- return pid;
-
+ return pid;
}
//----------------------------
Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie) const{
Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const{
// n-sigma cut, TOF PID
- if(!CheckTOFPIDStatus(track)) return 0;
+ Double_t nsigma;
+ Int_t pid=-1;
+ if(specie<0){
+ Double_t nsigmaMin=999.;
+ for(Int_t ipart=0;ipart<5;ipart++){
+ if(GetnSigmaTOF(track,ipart,nsigma)==1){
+ nsigma=TMath::Abs(nsigma);
+ if((nsigma<nsigmaMin)&& (nsigma<fnSigma[3])){
+ pid=ipart;
+ nsigmaMin=nsigma;
+ }
+ }
+ }
+ }else{ // asks only for one particle specie
+ if(GetnSigmaTOF(track,specie,nsigma)==1){
+ nsigma=TMath::Abs(nsigma);
+ if (nsigma>fnSigma[3]) pid=-1;
+ else pid=specie;
+ }
+ }
+ return pid;
+ /*
Double_t time[AliPID::kSPECIESN];
Double_t sigmaTOFPid[AliPID::kSPECIES];
AliAODPid *pidObj = track->GetDetPid();
pidObj->GetIntegratedTimes(time);
Double_t sigTOF=pidObj->GetTOFsignal();
- pidObj->GetTOFpidResolution(sigmaTOFPid);
+
+ AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
+ if (event) {
+ AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
+ if (tofH && fPidResponse) { // reading new AOD with new aliroot
+ AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
+ sigTOF -= TOFres.GetStartTime(track->P());
+ if (specie<0) {
+ for (Int_t ipart = 0; ipart<5; ipart++) {
+ sigmaTOFPid[ipart]=TOFres.GetExpectedSigma(track->P(),time[ipart],AliPID::ParticleMass(ipart));
+ }
+ }
+ else sigmaTOFPid[specie]=TOFres.GetExpectedSigma(track->P(),time[specie],AliPID::ParticleMass(specie)); //fTOFResponse is set in InitialiseEvent
+ } else pidObj->GetTOFpidResolution(sigmaTOFPid); // reading old AOD with new aliroot
+ } else pidObj->GetTOFpidResolution(sigmaTOFPid); //reading old AOD with old aliroot
Int_t pid=-1;
}
}
return pid;
-
+ */
}
//------------------------------
void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type) const{
return;
}
-//--------------------
-void AliAODPidHF::BayesianProbability(AliAODTrack *track,Double_t *pid) const{
-// bayesian PID for single detectors or combined
-
- if(fITS && !fTPC && !fTOF) {BayesianProbabilityITS(track,pid);return;}
- if(fTPC && !fITS && !fTOF) {BayesianProbabilityTPC(track,pid);return;}
- if(fTOF && !fITS && !fTPC) {BayesianProbabilityTOF(track,pid);return;}
-
- Double_t probITS[5]={1.,1.,1.,1.,1.};
- Double_t probTPC[5]={1.,1.,1.,1.,1.};
- Double_t probTOF[5]={1.,1.,1.,1.,1.};
- if(fITS) BayesianProbabilityITS(track,probITS);
- if(fTPC) BayesianProbabilityTPC(track,probTPC);
- if(fTOF) 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++){
- if(fTOF || fTPC || fTRD){
- prob[ind]=itspid[ind];
- }else{
- 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(fTOF || fITS || fTRD){
- prob[ind]=tpcpid[ind];
- }else{
- 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]);
- }
-}
- 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++){
- if(fTPC || fITS || fTRD){
- prob[ind]=tofpid[ind];
- }else{
- 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(fTPC || fITS || fTOF){
- prob[ind]=trdpid[ind];
- }else{
- 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]);
- }
-}
- return;
-
- }
//--------------------------------
Bool_t AliAODPidHF::CheckITSPIDStatus(AliAODTrack *track) const{
// ITS PID quality cuts
Bool_t AliAODPidHF::TPCRawAsym(AliAODTrack* track,Int_t specie) const{
// TPC nsigma cut PID, different sigmas in different p bins
- if(!CheckTPCPIDStatus(track)) return kFALSE;
AliAODPid *pidObj = track->GetDetPid();
Double_t mom = pidObj->GetTPCmomentum();
- if(mom>fPtThresholdTPC) return 0;
- Double_t nsigma=999.;
- if(fOldPid){
- Double_t dedx=pidObj->GetTPCsignal();
- UShort_t nTPCClus=pidObj->GetTPCsignalN();
- if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
-
- AliPID::EParticleType type=AliPID::EParticleType(specie);
- nsigma = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
+ if(mom>fPtThresholdTPC) return kTRUE;
+
+ Double_t nsigma;
+ if(GetnSigmaTPC(track,specie,nsigma)!=1) return kFALSE;
+ nsigma=TMath::Abs(nsigma);
- }else{ //old pid
- AliPID::EParticleType type=AliPID::EParticleType(specie);
- nsigma = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
- } //new pid
if(mom<fPLimit[0] && nsigma<fnSigma[0]) return kTRUE;
if(mom<fPLimit[1] && mom>fPLimit[0] && nsigma<fnSigma[1]) return kTRUE;
Int_t AliAODPidHF::MatchTPCTOF(AliAODTrack *track, Int_t specie){
// combination of the PID info coming from TPC and TOF
+ Double_t ptrack=track->P();
+ if(ptrack>fMaxTrackMomForCombinedPID) return 1;
+
Bool_t okTPC=CheckTPCPIDStatus(track);
Bool_t okTOF=CheckTOFPIDStatus(track);
+ if(ptrack>fPtThresholdTPC) okTPC=kFALSE;
if(fMatch==1){
//TOF || TPC (a la' Andrea R.)
tTOFinfo=-1;
if(ApplyPidTOFRaw(track,specie)==specie) tTOFinfo=1;
if(fCompat && tTOFinfo>0){
- Double_t ptrack=track->P();
if(ptrack>fPCompatTOF) {
Double_t sig0tmp=fnSigma[3];
SetSigma(3,fnSigmaCompat[1]);
// convention (temporary): -1 = kFALSE, 1 = kTRUE, 0 = not identified
if(fTPC && fTOF) if(!okTPC && !okTOF) return 0;
- Double_t ptrack=track->P();
-
+
Int_t tTPCinfo=-1;
if(ptrack>=fPLimit[0] && ptrack<fPLimit[1] && fTPC) {
if(!okTPC) return 0;
}
//-----------------------
Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
+ // TOF proton compatibility
+ if(!CheckTOFPIDStatus(track)) return 0;
- if(!CheckTOFPIDStatus(track)) return 0;
-
- Double_t time[AliPID::kSPECIESN];
+ Double_t nsigma;
+ if(GetnSigmaTOF(track,3,nsigma)==1){
+ if(nsigma>nsigmaK) return kTRUE;
+ }
+ return kFALSE;
+ /* Double_t time[AliPID::kSPECIESN];
Double_t sigmaTOFPid[AliPID::kSPECIES];
AliAODPid *pidObj = track->GetDetPid();
pidObj->GetIntegratedTimes(time);
Double_t sigTOF=pidObj->GetTOFsignal();
- pidObj->GetTOFpidResolution(sigmaTOFPid);
+
+ AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
+ if (event) {
+ AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
+ if (tofH && fPidResponse) {
+ AliTOFPIDResponse TOFres = (AliTOFPIDResponse)fPidResponse->GetTOFResponse();
+ sigTOF -= TOFres.GetStartTime(track->P());
+ sigmaTOFPid[3]=TOFres.GetExpectedSigma(track->P(),time[3],AliPID::ParticleMass(3));
+ }
+ else pidObj->GetTOFpidResolution(sigmaTOFPid);
+ } else pidObj->GetTOFpidResolution(sigmaTOFPid);
Double_t sigmaTOFtrack;
if (sigmaTOFPid[3]>0) sigmaTOFtrack=sigmaTOFPid[3];
else sigmaTOFtrack=fTOFSigma; // backward compatibility for old AODs
if((sigTOF-time[3])>nsigmaK*sigmaTOFtrack)return kTRUE;// K, Pi excluded (->LIKELY A PROTON)
return kFALSE;
-
+ */
}
//--------------------------------------------------------------------------
void AliAODPidHF::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior){
}
//--------------------------------------------------------------------------
-Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const{
-
+Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
+ // get n sigma for TPC
+
if(!CheckTPCPIDStatus(track)) return -1;
Double_t nsigmaTPC=-999;
UShort_t nTPCClus=pidObj->GetTPCsignalN();
if(nTPCClus==0) {nTPCClus=track->GetTPCNcls();}
AliPID::EParticleType type=AliPID::EParticleType(species);
- nsigmaTPC = TMath::Abs(fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type));
- sigma=nsigmaTPC;
+ nsigmaTPC = fTPCResponse->GetNumberOfSigmas(mom,dedx,nTPCClus,type);
+ nsigma=nsigmaTPC;
} else{
-
+ if(!fPidResponse) return -1;
AliPID::EParticleType type=AliPID::EParticleType(species);
- nsigmaTPC = TMath::Abs(fPidResponse->NumberOfSigmasTPC(track,type));
- sigma=nsigmaTPC;
+ nsigmaTPC = fPidResponse->NumberOfSigmasTPC(track,type);
+ nsigma=nsigmaTPC;
}
return 1;
}
//-----------------------------
-Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &sigma) const{
+Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
+ // get n sigma for TOF
if(!CheckTOFPIDStatus(track)) return -1;
-
- Double_t time[AliPID::kSPECIESN];
- Double_t sigmaTOFPid[AliPID::kSPECIES];
- AliAODPid *pidObj = track->GetDetPid();
- pidObj->GetIntegratedTimes(time);
- Double_t sigTOF=pidObj->GetTOFsignal();
- pidObj->GetTOFpidResolution(sigmaTOFPid);
-
- if(sigmaTOFPid[species]<1e-99) return -2;
-
- Double_t sigmaTOF=TMath::Abs(sigTOF-time[species])/sigmaTOFPid[species];
- sigma=sigmaTOF;
- return 1;
+
+ if(fPidResponse){
+ nsigma = fPidResponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)species);
+ return 1;
+ }else{
+ AliFatal("To use TOF PID you need to attach AliPIDResponseTask");
+ nsigma=-999.;
+ return -1;
+ }
}
+//-----------------------
+Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detectors) {
+ // Exclude a given hypothesis (labelTracks) in detector
+
+ if (detectors.Contains("ITS")) {
+
+ AliInfo("Nothing to be done");
+ /*
+ Double_t nsigma=0.;
+ if (GetnSigmaITS(track,labelTrack,nsigma)==1){
+ if(nsigma>nsigmaCut) return kTRUE;
+ }
+ */
+ return kFALSE;
+
+ } else if (detectors.Contains("TPC")) {
+
+ Double_t nsigma=0.;
+ if (GetnSigmaTPC(track,labelTrack,nsigma)==1){
+ if(nsigma>nsigmaCut) return kTRUE;
+ }
+ return kFALSE;
+
+ } else if (detectors.Contains("TOF")) {
+
+ if (!(CheckTOFPIDStatus(track))) return kFALSE;
+ Double_t nsigma=0.;
+ if (GetnSigmaTOF(track,labelTrack,nsigma)==1){
+ if(nsigma>nsigmaCut) return kTRUE;
+ }
+ return kFALSE;
+
+ }
+ return kFALSE;
+
+}
//-----------------------------
+void AliAODPidHF::SetPriorsHistos(TString priorFileName){
+ // Set histograms with priors
+
+ for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
+ if(fPriorsH[ispecies]) delete fPriorsH[ispecies];
+ TString nt ="name";
+ nt+="_prior_";
+ nt+=AliPID::ParticleName(ispecies);
+ }
+ TDirectory *current = gDirectory;
+ TFile *priorFile=TFile::Open(priorFileName);
+ if (priorFile) {
+ TH1F* h3=static_cast<TH1F*>(priorFile->Get("priors3step9"));
+ TH1F* h2=static_cast<TH1F*>(priorFile->Get("priors2step9"));
+ TH1F* h1=static_cast<TH1F*>(priorFile->Get("priors1step9"));
+ current->cd();
+ fPriorsH[AliPID::kProton] = new TH1F(*h3);
+ fPriorsH[AliPID::kKaon ] = new TH1F(*h2);
+ fPriorsH[AliPID::kPion ] = new TH1F(*h1);
+ priorFile->Close();
+ delete priorFile;
+ TF1 *salt=new TF1("salt","1.e-10",0,10);
+ fPriorsH[AliPID::kProton]->Add(salt);
+ fPriorsH[AliPID::kKaon ]->Add(salt);
+ fPriorsH[AliPID::kPion ]->Add(salt);
+ delete salt;
+ }
+}
+//----------------------------------
+void AliAODPidHF::SetUpCombinedPID(){
+ // Configuration of combined Bayesian PID
+
+ fPidCombined->SetSelectedSpecies(AliPID::kSPECIES);
+ for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
+ fPidCombined->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),fPriorsH[ispecies]);
+ }
+ switch (fCombDetectors){
+ case kTPCTOF:
+ fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+ break;
+ case kTPCITS:
+ fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetITS);
+ break;
+ case kTPC:
+ fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+ break;
+ case kTOF:
+ fPidCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+ break;
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff