* about the suitability of this software for any purpose. It is *
* provided "as is" without express or implied warranty. *
**************************************************************************/
-/************************************************************************
- * *
- * Class for TPC PID *
- * Implements the abstract base class AliHFEpidBase *
- * *
- * Class contains TPC specific cuts and QA histograms *
- * Two selection strategies are offered: Selection of certain value *
- * regions in the TPC dE/dx (by IsSelected), and likelihoods *
- * *
- * Authors: *
- * *
- * Markus Fasel <M.Fasel@gsi.de> *
- * Markus Heide <mheide@uni-muenster.de> *
- * *
- * *
- ************************************************************************/
+//
+// Class for TPC PID
+// Implements the abstract base class AliHFEpidBase
+//
+// Class contains TPC specific cuts and QA histograms
+// Two selection strategies are offered: Selection of certain value
+// regions in the TPC dE/dx (by IsSelected), and likelihoods
+//
+// Authors:
+//
+// Markus Fasel <M.Fasel@gsi.de>
+// Markus Heide <mheide@uni-muenster.de>
+//
#include <TH2I.h>
#include <TList.h>
#include <TMath.h>
+//#include <TParticle.h>
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
#include "AliESDtrack.h"
#include "AliExternalTrackParam.h"
#include "AliLog.h"
+#include "AliMCParticle.h"
#include "AliPID.h"
-#include "AliTPCpidESD.h"
-#include "AliVParticle.h"
+#include "AliESDpid.h"
+//#include "AliVParticle.h"
#include "AliHFEpidTPC.h"
//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const char* name) :
// add a list here
- AliHFEpidBase(name)
+ AliHFEpidBase(name)
+ , fLineCrossingType(0)
, fLineCrossingsEnabled(0)
- , fNsigmaTPC(2)
- , fPID(0x0)
- , fPIDtpcESD(0x0)
- , fQAList(0x0)
+ , fNsigmaTPC(3)
+ , fRejectionEnabled(0)
+ , fPID(NULL)
+ , fESDpid(NULL)
+ , fQAList(NULL)
{
//
// default constructor
//
- memset(fLineCrossingCenter, 0, sizeof(Double_t) * AliPID::kSPECIES);
memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
+ memset(fPAsigCut, 0, sizeof(Float_t) * 2);
+ memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
fPID = new AliPID;
- fPIDtpcESD = new AliTPCpidESD;
+ fESDpid = new AliESDpid;
}
//___________________________________________________________________
AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
- AliHFEpidBase("")
+ AliHFEpidBase("")
+ , fLineCrossingType(0)
, fLineCrossingsEnabled(0)
, fNsigmaTPC(2)
- , fPID(0x0)
- , fPIDtpcESD(0x0)
- , fQAList(0x0)
+ , fRejectionEnabled(0)
+ , fPID(NULL)
+ , fESDpid(NULL)
+ , fQAList(NULL)
{
//
// Copy constructor
return *this;
}
+//___________________________________________________________________
void AliHFEpidTPC::Copy(TObject &o) const{
//
// Copy function
target.fLineCrossingsEnabled = fLineCrossingsEnabled;
target.fNsigmaTPC = fNsigmaTPC;
+ target.fRejectionEnabled = fRejectionEnabled;
target.fPID = new AliPID(*fPID);
- target.fPIDtpcESD = new AliTPCpidESD(*fPIDtpcESD);
+ target.fESDpid = new AliESDpid(*fESDpid);
target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
-
+ memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
+ memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
+ memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
+
AliHFEpidBase::Copy(target);
}
// Destructor
//
if(fPID) delete fPID;
- if(fPIDtpcESD) delete fPIDtpcESD;
+ if(fESDpid) delete fESDpid;
if(fQAList){
fQAList->Delete();
delete fQAList;
//
// Add TPC dE/dx Line crossings
//
- AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
- AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
+ //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
+ //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
return kTRUE;
}
//___________________________________________________________________
-Int_t AliHFEpidTPC::IsSelected(AliVParticle *track)
+Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
{
//
// For the TPC pid we use the 2-sigma band around the bethe bloch curve
// for electrons
// exclusion of the crossing points
//
- AliESDtrack *esdTrack = 0x0;
- if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
- if(IsQAon()) FillTPChistograms(esdTrack);
- Double_t TPCsignal = esdTrack->GetTPCsignal();
+ if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
+ AliESDtrack *esdTrack = NULL;
+ if(!(esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack))) return 0;
+ AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track->fMCtrack);
+ return MakePIDesd(esdTrack, mctrack);
+ }else{
+ AliAODTrack *aodtrack = NULL;
+ if(!(aodtrack = dynamic_cast<AliAODTrack *>(track->fRecTrack))) return 0;
+ AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
+ return MakePIDaod(aodtrack, aodmctrack);
+ }
+}
+//___________________________________________________________________
+Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
+ //
+ // Doing TPC PID as explained in IsSelected for ESD tracks
+ //
+ if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
+ Float_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron);
// exclude crossing points:
// Determine the bethe values for each particle species
- Double_t p = esdTrack->GetInnerParam()->P();
Bool_t isLineCrossing = kFALSE;
+ fLineCrossingType = 0; // default value
for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
+ if(ispecies == AliPID::kElectron) continue;
if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
- if(TMath::Abs(p - fLineCrossingCenter[ispecies]) < fLineCrossingSigma[ispecies]){
- // Point in a line crossing region, no PID possible
+ if(TMath::Abs(fESDpid->NumberOfSigmasTPC(esdTrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
+ // Point in a line crossing region, no PID possible, but !PID still possible ;-)
isLineCrossing = kTRUE;
+ fLineCrossingType = ispecies;
break;
}
}
if(isLineCrossing) return 0;
+
+ // Check particle rejection
+ if(HasParticleRejection()){
+ Int_t reject = Reject(esdTrack);
+ if(reject != 0) return reject;
+ }
// Check whether distance from the electron line is smaller than n-sigma
- Double_t beta = p/fPID->ParticleMass(AliPID::kElectron);
- if(TMath::Abs(TPCsignal - 50*fPIDtpcESD->Bethe(beta)) < GetTPCsigma(p,0)) return 11;
+
+ // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
+ Float_t p = 0.;
+ Int_t pdg = 0;
+ if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
+ if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
+ } else {
+ if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
+ }
+ if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
+
+ return pdg;
+}
+
+//___________________________________________________________________
+Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mctrack*/){
+ AliError("AOD PID not yet implemented");
return 0;
}
//___________________________________________________________________
-void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t p, Double_t sigma_p){
+Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
+ //
+ // reject particles based on asymmetric sigma cut
+ //
+ Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
+ Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
+ for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
+ if(!TESTBIT(fRejectionEnabled, ispec)) continue;
+ // Particle rejection enabled
+ if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
+ Double_t sigma = fESDpid->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
+ if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
+ }
+ return 0;
+}
+
+//___________________________________________________________________
+void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
//
// Add exclusion point for the TPC PID where a dEdx line crosses the electron line
// Stores line center and line sigma
return;
}
fLineCrossingsEnabled |= 1 << species;
- fLineCrossingCenter[species] = p;
- fLineCrossingSigma[species] = sigma_p;
-}
-
-//___________________________________________________________________
-Double_t AliHFEpidTPC::GetTPCsigma(Double_t p, Int_t species){
- //
- // return the TPC sigma, momentum dependent
- //
- if(p < 0.1 || p > 20.) return 0.;
- Double_t beta = p/fPID->ParticleMass(species);
-
-
- return 50*fPIDtpcESD->Bethe(beta) * 0.06;
+ fLineCrossingSigma[species] = sigma;
}
//___________________________________________________________________
//IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
if(!track) return -1.;
- Int_t hypo; //marks particle hypotheses for 2-sigma bands
- Double_t beta;//well o.k., it corresponds to gamma * beta
- Double_t p = track->GetInnerParam()->P();
- Double_t TPCsignal = track->GetTPCsignal();
Bool_t outlier = kTRUE;
// Check whether distance from the respective particle line is smaller than r sigma
- for(hypo = 0; hypo < 5; hypo++)
- {
- beta = p/fPID->ParticleMass(hypo);
- if(TMath::Abs(TPCsignal - (GetTPCsigma(p, hypo))/0.06) > (rsig * GetTPCsigma(p,hypo)))
- outlier = kTRUE;
- else
- {
- outlier = kFALSE;
- break;
- }
- }
+ for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
+ if(TMath::Abs(fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)hypo)) > rsig)
+ outlier = kTRUE;
+ else {
+ outlier = kFALSE;
+ break;
+ }
+ }
if(outlier)
return -2.;
- Double_t TPCprob[5];
+ Double_t tpcProb[5];
- track->GetTPCpid(TPCprob);
+ track->GetTPCpid(tpcProb);
- return TPCprob[species];
+ return tpcProb[species];
}
-//___________________________________________________________________
-Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species)
-{
- //default: rsig = 2.
- // for everything else, see above!
- if(!track) return -1.;
- Int_t hypo; //marks particle hypotheses for 2-sigma bands
- Double_t beta;
- Double_t p = track->GetInnerParam()->P();
- Double_t TPCsignal = track->GetTPCsignal();
- Bool_t outlier = kTRUE;
- // Check whether distance from the respective particle line is smaller than 2 sigma
- for(hypo = 0; hypo < 5; hypo++)
- {
- beta = p/fPID->ParticleMass(hypo);
-
- if(TMath::Abs(TPCsignal - (GetTPCsigma(p, hypo))/0.06) > (2. * GetTPCsigma(p,hypo)))
- outlier = kTRUE;
- else
- {
- outlier = kFALSE;
- break;
- }
- }
- if(outlier == kTRUE)
- return -2.;
-
- Double_t TPCprob[5];
-
- track->GetTPCpid(TPCprob);
-
- return TPCprob[species];
-}
//___________________________________________________________________
Double_t AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
{
//ratio of likelihoods to be whatever species/to be an electron;
//as a cross-check for possible particle type suppression compared to electrons
+ const Double_t kVerySmall = 1e-10;
if(!track) return -20;
- if((Likelihood(track,species) == -2.)||(Likelihood(track,0)== -2.))
+ if((TMath::Abs(Likelihood(track,species) + 2) < kVerySmall)||(TMath::Abs(Likelihood(track,0) + 2 ) < kVerySmall))
return -30;
- if(Likelihood(track,species) == 0.)
+ if(TMath::Abs(Likelihood(track,species)) < kVerySmall)
return -10;
- if (Likelihood(track,0) == 0.)
+ if (TMath::Abs(Likelihood(track,0)) < kVerySmall)
return 10.;
else
return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
}
+
//___________________________________________________________________
-void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track){
+void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track, const AliMCParticle *mctrack){
+ //
+ // Fill the QA histogtrams
//
if(!track)
return;
- Double_t tpc_signal = track->GetTPCsignal();
+ Double_t tpcSignal = track->GetTPCsignal();
Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
if(HasMCData()){
- switch(TMath::Abs(GetPdgCode(const_cast<AliESDtrack *>(track)))){
- case 11: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpc_signal);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
- //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
- //___________________________________________________________________________________________
- //Likelihoods for electrons to be other particle species
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
- break;
- //___________________________________________________________________________________________
- case 13: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpc_signal);
- //Likelihood of muon to be an electron
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
- //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
- //below functions are the same for other species
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
- break;
- case 211: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpc_signal);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
- break;
- case 321: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpc_signal);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
- break;
- case 2212: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpc_signal);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
- break;
- default: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpc_signal);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
-
- break;
+ switch(TMath::Abs(mctrack->Particle()->GetPdgCode())){
+ case 11: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
+ //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
+ //___________________________________________________________________________________________
+ //Likelihoods for electrons to be other particle species
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
+ break;
+ //___________________________________________________________________________________________
+ case 13: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
+ //Likelihood of muon to be an electron
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
+ //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
+ //below functions are the same for other species
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
+ break;
+ case 211: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
+ break;
+ case 321: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
+ break;
+ case 2212: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
+ break;
+ default: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
+
+ break;
}
}
//TPC signal and Likelihood to be electron for all tracks (independent of MC information)
- (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpc_signal);
- (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
-
-
-
+ (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
+ (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
}
//___________________________________________________________________
void AliHFEpidTPC::AddQAhistograms(TList *qaList){
+ //
+ // Create QA histograms for TPC PID
+ //
fQAList = new TList;
fQAList->SetName("fTPCqaHistos");
fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
+ fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 60, 0, 600), kHistTPCselected);
fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff