#include "TParticle.h"
#include "AliAODMCParticle.h"
#include "AliPIDResponse.h"
+#include "AliPIDCombined.h"
#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
using namespace std;
const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
-const char * kDetectorName[]={"TPC","TOF"} ;
+const char * kDetectorName[]={"ITS","TPC","TOF"} ;
const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
ClassImp(AliHelperPID)
-AliHelperPID::AliHelperPID() : TNamed("HelperPID", "PID object"),fisMC(0), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fPIDResponse(0),fOutputList(0),fRequestTOFPID(1),fPtTOFPID(.6),fHasTOFPID(0){
+AliHelperPID::AliHelperPID() : TNamed("HelperPID", "PID object"),fisMC(0), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fBayesCut(0.8), fPIDResponse(0x0), fPIDCombined(0x0),fOutputList(0x0),fRequestTOFPID(1),fRemoveTracksT0Fill(0),fUseExclusiveNSigma(0),fPtTOFPID(.6),fHasTOFPID(0){
for(Int_t ipart=0;ipart<kNSpecies;ipart++)
for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++)
Double_t miny=-30;
Double_t maxy=30;
if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
- TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,200,miny,maxy);
+ TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
fOutputList->Add(fHistoNSigma);
Double_t maxy=10;
if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
- Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
+ Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
fOutputList->Add(fHistoNSigma);
}
}
+ //BayesRec plot
+ for(Int_t ipart=0;ipart<kNSpecies;ipart++){
+ Double_t miny=0.;
+ Double_t maxy=1;
+ TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
+ Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
+ fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
+ fOutputList->Add(fHistoBayes);
+ }
+
//nsigmaDC plot
for(Int_t ipart=0;ipart<kNSpecies;ipart++){
for(Int_t ipid=0;ipid<=kNSigmaPIDType;ipid++){
Double_t maxy=10;
if(ipid==kNSigmaTPCTOF){miny=0;maxy=20;}
TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
- Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
+ Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
fOutputList->Add(fHistoNSigma);
Double_t maxy=30;
if(ipid==kNSigmaTPCTOF){miny=0;maxy=50;}
TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
- Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,100,miny,maxy);
+ Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
fOutputList->Add(fHistoNSigma);
//PID signal plot
for(Int_t idet=0;idet<kNDetectors;idet++){
- Double_t maxy=1000;
- if(idet==kTOF)maxy=200000;
- TH2F *fHistoPID=new TH2F(Form("PID_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,200,-maxy,maxy);
+ for(Int_t ipart=0;ipart<kNSpecies;ipart++){
+ Double_t maxy=500;
+ if(idet==kTOF)maxy=1.1;
+ TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
+ fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
+ fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
+ fOutputList->Add(fHistoPID);
+ }
+ }
+ //PID signal plot, before PID cut
+ for(Int_t idet=0;idet<kNDetectors;idet++){
+ Double_t maxy=500;
+ if(idet==kTOF)maxy=1.1;
+ TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
fOutputList->Add(fHistoPID);
Int_t AliHelperPID::GetParticleSpecies(AliVTrack * trk, Bool_t FIllQAHistos){
//return the specie according to the minimum nsigma value
//no double counting, this has to be evaluated using CheckDoubleCounting()
- //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);
- CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
-
- return FindMinNSigma(trk,FIllQAHistos);//NSigmaRec distr filled here
+ Int_t ID=kSpUndefined;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////
-
-void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
- //defines data member fnsigmas
+ //get the PID response
if(!fPIDResponse) {
AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
AliFatal("Cannot get pid response");
}
+ //calculate nsigmas (used also by the bayesian)
+ CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
+
+ //Do PID
+ if(fPIDType==kBayes){//use bayesianPID
+
+ if(!fPIDCombined) {
+ // ------- setup PIDCombined
+ fPIDCombined=new AliPIDCombined;
+ fPIDCombined->SetDefaultTPCPriors();
+ fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
+ }
+ if(!fPIDCombined) {
+ AliFatal("PIDCombined object not found");
+ }
+
+ ID = GetIDBayes(trk,FIllQAHistos);
+
+ }else{ //use nsigma PID
+
+ ID=FindMinNSigma(trk,FIllQAHistos);
+
+ if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
+ Bool_t *HasDC;
+ HasDC=GetDoubleCounting(trk,FIllQAHistos);
+ for(Int_t ipart=0;ipart<kNSpecies;ipart++){
+ if(HasDC[ipart]==kTRUE) ID = kSpUndefined;
+ }
+ }
+ }
+
+ //Fill PID signal plot
+ if(ID != kSpUndefined){
+ for(Int_t idet=0;idet<kNDetectors;idet++){
+ TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
+ if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
+ if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
+ if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
+ }
+ }
+ //Fill PID signal plot without cuts
+ for(Int_t idet=0;idet<kNDetectors;idet++){
+ TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
+ if(idet==kITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
+ if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
+ if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),TOFBetaCalc(trk)*trk->Charge());
+ }
+
+ return ID;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+Int_t AliHelperPID::GetParticleSpecies(AliVParticle * part) {
+ // return PID according to MC truth
+ switch(TMath::Abs(part->PdgCode())){
+ case 2212:
+ return kSpProton;
+ break;
+ case 321:
+ return kSpKaon;
+ break;
+ case 211:
+ return kSpPion;
+ break;
+ default:
+ return kSpUndefined;
+ }
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+Int_t AliHelperPID::GetIDBayes(AliVTrack * trk, Bool_t FIllQAHistos){
+
+ Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
+
+ Double_t probBayes[AliPID::kSPECIES];
+
+ UInt_t detUsed= 0;
+ CheckTOF(trk);
+ if(fHasTOFPID && trk->Pt()>fPtTOFPID){//use TOF information
+ detUsed = CalcPIDCombined(trk, fPIDResponse, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
+ if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return kSpUndefined;//check that TPC and TOF are used
+ }else{
+ detUsed = CalcPIDCombined(trk, fPIDResponse,AliPIDResponse::kDetTPC, probBayes);
+ if(detUsed != AliPIDResponse::kDetTPC)return kSpUndefined;//check that TPC is used
+ }
+
+ //the probability has to be normalized to one, we check it
+ Double_t sump=0.;
+ for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
+ if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
+ AliFatal("Bayesian probability not normalized to one");
+ }
+
+ //probabilities are normalized to one, if the cut is above .5 there is no problem
+ if(probBayes[AliPID::kPion]>fBayesCut && IDs[kSpPion]==1){
+ TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpPion));
+ h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
+ return kSpPion;
+ }
+ else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[kSpKaon]==1){
+ TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpKaon));
+ h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
+ return kSpKaon;
+ }
+ else if(probBayes[AliPID::kProton]>fBayesCut && IDs[kSpProton]==1){
+ TH2F *h=GetHistogram2D(Form("BayesRec_%d",kSpProton));
+ h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
+ return kSpProton;
+ }
+ else{
+ return kSpUndefined;
+ }
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+UInt_t AliHelperPID::CalcPIDCombined(const AliVTrack *track,const AliPIDResponse *PIDResponse, Int_t detMask, Double_t* prob) const{
+ //
+ // Bayesian PID calculation
+ //
+ for(Int_t i=0;i<AliPID::kSPECIES;i++)
+ {
+ prob[i]=0.;
+ }
+ fPIDCombined->SetDetectorMask(detMask);
+
+ return fPIDCombined->ComputeProbabilities((AliVTrack*)track, PIDResponse, prob);
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+void AliHelperPID::CalculateNSigmas(AliVTrack * trk, Bool_t FIllQAHistos){
+ //defines data member fnsigmas
+
// Compute nsigma for each hypthesis
AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
// --- TPC
nsigmaTOFkProton = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton);
nsigmaTOFkKaon = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon);
nsigmaTOFkPion = fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion);
+ Double_t d2Proton=nsigmaTPCkProton * nsigmaTPCkProton + nsigmaTOFkProton * nsigmaTOFkProton;
+ Double_t d2Kaon=nsigmaTPCkKaon * nsigmaTPCkKaon + nsigmaTOFkKaon * nsigmaTOFkKaon;
+ Double_t d2Pion=nsigmaTPCkPion * nsigmaTPCkPion + nsigmaTOFkPion * nsigmaTOFkPion;
+ //commented, this is done in analogy with AliESDTrackCuts, nsigma combind according to the probability
// --- combined
- nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
- nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
- nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
+ // -----------------------------------
+ // How to get to a n-sigma cut?
+ //
+ // The accumulated statistics from 0 to d is
+ //
+ // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
+ // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
+ // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
+ //
+ // work around precision problem
+ // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
+ //if(TMath::Exp(- d2Proton / 2) > 1e-15)nsigmaTPCTOFkProton = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Proton / 2));
+ //if(TMath::Exp(- d2Kaon / 2) > 1e-15)nsigmaTPCTOFkKaon = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Kaon / 2));
+ //if(TMath::Exp(- d2Pion / 2) > 1e-15)nsigmaTPCTOFkPion = TMath::Sqrt(2)*TMath::ErfInverse(1 - TMath::Exp(- d2Pion / 2));
+
+ //used for the 2PC PID paper (circular cut)
+ nsigmaTPCTOFkProton = TMath::Sqrt(d2Proton);
+ nsigmaTPCTOFkKaon = TMath::Sqrt(d2Kaon);
+ nsigmaTPCTOFkPion = TMath::Sqrt(d2Pion);
}else{
// --- combined
// if TOF is missing and below fPtTOFPID only the TPC information is used
h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
}
}
- //Fill PID signal plot
- for(Int_t idet=0;idet<kNDetectors;idet++){
- TH2F *h=GetHistogram2D(Form("PID_%d",idet));
- if(idet==kTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
- if(idet==kTOF && fHasTOFPID)h->Fill(trk->P(),trk->GetTOFsignal()*trk->Charge());
- }
}
}
nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
break;
+ case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
+ nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
+ nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
+ break;
}
// guess the particle based on the smaller nsigma (within fNSigmaPID)
Bool_t* AliHelperPID::GetDoubleCounting(AliVTrack * trk,Bool_t FIllQAHistos){
//if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
//fill DC histos
- for(Int_t ipart=0;ipart<=kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
+ for(Int_t ipart=0;ipart<kNSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
break;
+ case kBayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
+ nsigmaProton = TMath::Abs(fnsigmas[kSpProton][kNSigmaTPCTOF]);
+ nsigmaKaon = TMath::Abs(fnsigmas[kSpKaon][kNSigmaTPCTOF]) ;
+ nsigmaPion = TMath::Abs(fnsigmas[kSpPion][kNSigmaTPCTOF]) ;
+ break;
}
if(nsigmaPion<fNSigmaPID && MinNSigma!=kSpPion)fHasDoubleCounting[kSpPion]=kTRUE;
if(nsigmaKaon<fNSigmaPID && MinNSigma!=kSpKaon)fHasDoubleCounting[kSpKaon]=kTRUE;
//////////////////////////////////////////////////////////////////////////////////////////////////
+Bool_t* AliHelperPID::GetAllCompatibleIdentitiesNSigma(AliVTrack * trk,Bool_t FIllQAHistos){
+
+ Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
+ IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
+ return IDs;
+
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
Int_t AliHelperPID::GetMCParticleSpecie(AliVEvent* event, AliVTrack * trk, Bool_t FillQAHistos){
//return the specie according to the MC truth
CheckTOF(trk);
void AliHelperPID::CheckTOF(AliVTrack * trk)
{
//check if the particle has TOF Matching
- UInt_t status;
- status=trk->GetStatus();
- if((status&AliVTrack::kTOFout)==0 || (status&AliVTrack::kTIME)==0)fHasTOFPID=kFALSE;
+
+ //get the PIDResponse
+ if(fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF,trk)==0)fHasTOFPID=kFALSE;
else fHasTOFPID=kTRUE;
- //in addition to KTOFout and kTIME we look at the pt
+
+ //in addition to TOF status we look at the pt
if(trk->Pt()<fPtTOFPID)fHasTOFPID=kFALSE;
+
+ if(fRemoveTracksT0Fill)
+ {
+ Int_t startTimeMask = fPIDResponse->GetTOFResponse().GetStartTimeMask(trk->P());
+ if (startTimeMask < 0)fHasTOFPID=kFALSE;
+ }
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+Double_t AliHelperPID::TOFBetaCalc(AliVTrack *track) const{
+ //TOF beta calculation
+ Double_t tofTime=track->GetTOFsignal();
+
+ Double_t c=TMath::C()*1.E-9;// m/ns
+ Float_t startTime = fPIDResponse->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
+ Double_t length= fPIDResponse->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
+ tofTime -= startTime; // subtract startTime to the signal
+ Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
+ tof=tof*c;
+ return length/tof;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////
+
+Double_t AliHelperPID::GetMass(AliHelperParticleSpecies_t id) const{
+ //return Mass according to AliHelperParticleSpecies_t. If undefined return -999.
+ Double_t mass=-999.;
+ if (id == kSpProton) { mass=9.38271999999999995e-01; }
+ if (id == kSpKaon) { mass=4.93676999999999977e-01; }
+ if (id == kSpPion) { mass=1.39570000000000000e-01; }
+ return mass;
}
//////////////////////////////////////////////////////////////////////////////////////////////////