1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Implements the abstract base class AliHFEpidBase
19 // Class contains TPC specific cuts and QA histograms
20 // Two selection strategies are offered: Selection of certain value
21 // regions in the TPC dE/dx (by IsSelected), and likelihoods
25 // Markus Fasel <M.Fasel@gsi.de>
26 // Markus Heide <mheide@uni-muenster.de>
31 #include "AliAODpidUtil.h"
32 #include "AliAODPid.h"
33 #include "AliAODTrack.h"
34 #include "AliAODMCParticle.h"
35 #include "AliESDtrack.h"
36 #include "AliExternalTrackParam.h"
38 #include "AliMCParticle.h"
40 #include "AliESDpid.h"
42 #include "AliHFEpidTPC.h"
43 #include "AliHFEpidQAmanager.h"
45 ClassImp(AliHFEpidTPC)
47 //___________________________________________________________________
48 AliHFEpidTPC::AliHFEpidTPC() :
51 , fLineCrossingsEnabled(0)
52 , fUpperSigmaCut(NULL)
53 , fLowerSigmaCut(NULL)
55 , fRejectionEnabled(0)
59 // default constructor
63 //___________________________________________________________________
64 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
67 , fLineCrossingsEnabled(0)
68 , fUpperSigmaCut(NULL)
69 , fLowerSigmaCut(NULL)
71 , fRejectionEnabled(0)
75 // default constructor
77 memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
78 memset(fPAsigCut, 0, sizeof(Float_t) * 2);
79 memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
83 //___________________________________________________________________
84 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
86 , fLineCrossingsEnabled(0)
87 , fUpperSigmaCut(NULL)
88 , fLowerSigmaCut(NULL)
90 , fRejectionEnabled(0)
99 //___________________________________________________________________
100 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
102 // Assignment operator
109 //___________________________________________________________________
110 void AliHFEpidTPC::Copy(TObject &o) const{
113 // called in copy constructor and assigment operator
115 AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
117 target.fLineCrossingsEnabled = fLineCrossingsEnabled;
118 target.fUpperSigmaCut = fUpperSigmaCut;
119 target.fLowerSigmaCut = fLowerSigmaCut;
120 target.fNsigmaTPC = fNsigmaTPC;
121 target.fRejectionEnabled = fRejectionEnabled;
122 target.fPID = new AliPID(*fPID);
123 memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
124 memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
125 memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
127 AliHFEpidBase::Copy(target);
130 //___________________________________________________________________
131 AliHFEpidTPC::~AliHFEpidTPC(){
135 if(fPID) delete fPID;
138 //___________________________________________________________________
139 Bool_t AliHFEpidTPC::InitializePID(){
141 // Add TPC dE/dx Line crossings
143 //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
144 //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
148 //___________________________________________________________________
149 Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
152 // For the TPC pid we use the 2-sigma band around the bethe bloch curve
154 // exclusion of the crossing points
157 if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())) return 0;
159 // QA before selection
160 if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
161 AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
162 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
163 Float_t nsigma = NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
164 AliDebug(1, Form("TPC NSigma: %f", nsigma));
165 // exclude crossing points:
166 // Determine the bethe values for each particle species
167 Bool_t isLineCrossing = kFALSE;
168 for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
169 if(ispecies == AliPID::kElectron) continue;
170 if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
171 if(TMath::Abs(NumberOfSigmas(track->GetRecTrack(), (AliPID::EParticleType)ispecies, anatype)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
172 // Point in a line crossing region, no PID possible, but !PID still possible ;-)
173 isLineCrossing = kTRUE;
177 if(isLineCrossing) return 0;
179 // Check particle rejection
180 if(HasParticleRejection()){
181 Int_t reject = Reject(track->GetRecTrack(), anatype);
182 if(reject != 0) return reject;
185 // Check if we have an asymmetric sigma model set
187 if(fUpperSigmaCut || fLowerSigmaCut){
188 pdg = CutSigmaModel(track->GetRecTrack(), anatype) ? 11 : 0;
190 // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
192 if(HasAsymmetricSigmaCut() && (p = track->GetRecTrack()->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
193 if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
195 if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
198 if(pidqa && pdg != 0) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
203 //___________________________________________________________________
204 Bool_t AliHFEpidTPC::CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
206 // N SigmaCut using parametrization of the cuts
208 Bool_t isSelected = kTRUE;
209 Float_t nsigma = NumberOfSigmas(track, AliPID::kElectron, anaType);
210 Double_t p = GetP(track, anaType);
211 if(fUpperSigmaCut && nsigma > fUpperSigmaCut->Eval(p)) isSelected = kFALSE;
212 if(fLowerSigmaCut && nsigma < fLowerSigmaCut->Eval(p)) isSelected = kFALSE;
216 //___________________________________________________________________
217 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
219 // reject particles based on asymmetric sigma cut
221 Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
222 Double_t p = GetP(track, anaType);
223 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
224 if(!TESTBIT(fRejectionEnabled, ispec)) continue;
225 // Particle rejection enabled
226 if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
227 Double_t sigma = NumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec), anaType);
228 if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
233 //___________________________________________________________________
234 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
236 // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
237 // Stores line center and line sigma
239 if(species >= AliPID::kSPECIES){
240 AliError("Species doesn't exist");
243 fLineCrossingsEnabled |= 1 << species;
244 fLineCrossingSigma[species] = sigma;
247 //___________________________________________________________________
248 Double_t AliHFEpidTPC::NumberOfSigmas(const AliVParticle *track, AliPID::EParticleType species, AliHFEpidObject::AnalysisType_t anaType) const {
250 // Get the number of sigmas
252 Double_t nSigmas = 100;
253 if(anaType == AliHFEpidObject::kESDanalysis){
255 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
256 if(esdtrack && fESDpid) nSigmas = fESDpid->NumberOfSigmasTPC(esdtrack, species);
258 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
259 if(aodtrack && fAODpid) nSigmas = fAODpid->NumberOfSigmasTPC(aodtrack, species);
264 //___________________________________________________________________
265 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
267 // Get the momentum at the inner wall of the TPC
270 if(anatype == AliHFEpidObject::kESDanalysis){
271 // ESD analysis: Use Inner Params for the momentum estimate
272 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
273 p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
275 // AOD analysis: Use TPC momentum stored in the AliAODpid object
276 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
277 p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();