]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEpidTPC.cxx
Updates + addition of EMCal
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidTPC.cxx
CommitLineData
809a4336 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
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**************************************************************************/
50685501 15//
16// Class for TPC PID
17// Implements the abstract base class AliHFEpidBase
18//
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
22//
23// Authors:
24//
25// Markus Fasel <M.Fasel@gsi.de>
26// Markus Heide <mheide@uni-muenster.de>
27//
faee3b18 28#include <TF1.h>
809a4336 29#include <TMath.h>
30
3a72645a 31#include "AliAODpidUtil.h"
6555e2ad 32#include "AliAODPid.h"
722347d8 33#include "AliAODTrack.h"
34#include "AliAODMCParticle.h"
809a4336 35#include "AliESDtrack.h"
36#include "AliExternalTrackParam.h"
37#include "AliLog.h"
722347d8 38#include "AliMCParticle.h"
809a4336 39#include "AliPID.h"
10d100d4 40#include "AliESDpid.h"
809a4336 41
42#include "AliHFEpidTPC.h"
3a72645a 43#include "AliHFEpidQAmanager.h"
809a4336 44
faee3b18 45ClassImp(AliHFEpidTPC)
809a4336 46
3a72645a 47//___________________________________________________________________
48AliHFEpidTPC::AliHFEpidTPC() :
49 // add a list here
50 AliHFEpidBase()
3a72645a 51 , fLineCrossingsEnabled(0)
52 , fUpperSigmaCut(NULL)
53 , fLowerSigmaCut(NULL)
bf892a6a 54 , fElectronMeanCorrection(NULL)
3a72645a 55 , fNsigmaTPC(3)
56 , fRejectionEnabled(0)
57 , fPID(NULL)
58{
59 //
60 // default constructor
61 //
bf892a6a 62 memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
63 memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
64 memset(fPAsigCut, 0, sizeof(Float_t) * 2);
65 memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
66
3a72645a 67}
68
809a4336 69//___________________________________________________________________
70AliHFEpidTPC::AliHFEpidTPC(const char* name) :
71 // add a list here
9bcfd1ab 72 AliHFEpidBase(name)
809a4336 73 , fLineCrossingsEnabled(0)
faee3b18 74 , fUpperSigmaCut(NULL)
75 , fLowerSigmaCut(NULL)
bf892a6a 76 , fElectronMeanCorrection(NULL)
75d81601 77 , fNsigmaTPC(3)
0792aa82 78 , fRejectionEnabled(0)
75d81601 79 , fPID(NULL)
809a4336 80{
81 //
82 // default constructor
83 //
bf892a6a 84 //
85 memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
809a4336 86 memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
722347d8 87 memset(fPAsigCut, 0, sizeof(Float_t) * 2);
88 memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
809a4336 89 fPID = new AliPID;
809a4336 90}
91
92//___________________________________________________________________
93AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
9bcfd1ab 94 AliHFEpidBase("")
809a4336 95 , fLineCrossingsEnabled(0)
faee3b18 96 , fUpperSigmaCut(NULL)
97 , fLowerSigmaCut(NULL)
bf892a6a 98 , fElectronMeanCorrection(NULL)
809a4336 99 , fNsigmaTPC(2)
0792aa82 100 , fRejectionEnabled(0)
75d81601 101 , fPID(NULL)
809a4336 102{
103 //
104 // Copy constructor
105 //
106 ref.Copy(*this);
107}
108
109//___________________________________________________________________
110AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
111 //
112 // Assignment operator
113 //
114 if(this != &ref){
115 ref.Copy(*this);
116 }
117 return *this;
118}
75d81601 119//___________________________________________________________________
809a4336 120void AliHFEpidTPC::Copy(TObject &o) const{
121 //
122 // Copy function
123 // called in copy constructor and assigment operator
124 //
125 AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
126
127 target.fLineCrossingsEnabled = fLineCrossingsEnabled;
faee3b18 128 target.fUpperSigmaCut = fUpperSigmaCut;
129 target.fLowerSigmaCut = fLowerSigmaCut;
bf892a6a 130 target.fElectronMeanCorrection = fElectronMeanCorrection;
809a4336 131 target.fNsigmaTPC = fNsigmaTPC;
0792aa82 132 target.fRejectionEnabled = fRejectionEnabled;
809a4336 133 target.fPID = new AliPID(*fPID);
75d81601 134 memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
722347d8 135 memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
136 memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
137
809a4336 138 AliHFEpidBase::Copy(target);
139}
140
141//___________________________________________________________________
142AliHFEpidTPC::~AliHFEpidTPC(){
143 //
144 // Destructor
145 //
146 if(fPID) delete fPID;
809a4336 147}
148
149//___________________________________________________________________
150Bool_t AliHFEpidTPC::InitializePID(){
151 //
152 // Add TPC dE/dx Line crossings
153 //
75d81601 154 //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
155 //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
809a4336 156 return kTRUE;
157}
158
159//___________________________________________________________________
6555e2ad 160Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
809a4336 161{
162 //
163 // For the TPC pid we use the 2-sigma band around the bethe bloch curve
164 // for electrons
165 // exclusion of the crossing points
166 //
3a72645a 167
168 if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())) return 0;
169
170 // QA before selection
171 if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
e3fc062d 172 AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
3a72645a 173 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
174 Float_t nsigma = NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
175 AliDebug(1, Form("TPC NSigma: %f", nsigma));
809a4336 176 // exclude crossing points:
177 // Determine the bethe values for each particle species
809a4336 178 Bool_t isLineCrossing = kFALSE;
179 for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
75d81601 180 if(ispecies == AliPID::kElectron) continue;
809a4336 181 if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
3a72645a 182 if(TMath::Abs(NumberOfSigmas(track->GetRecTrack(), (AliPID::EParticleType)ispecies, anatype)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
9bcfd1ab 183 // Point in a line crossing region, no PID possible, but !PID still possible ;-)
70da6c5a 184 isLineCrossing = kTRUE;
809a4336 185 break;
186 }
187 }
188 if(isLineCrossing) return 0;
0792aa82 189
190 // Check particle rejection
191 if(HasParticleRejection()){
3a72645a 192 Int_t reject = Reject(track->GetRecTrack(), anatype);
0792aa82 193 if(reject != 0) return reject;
194 }
722347d8 195
faee3b18 196 // Check if we have an asymmetric sigma model set
0792aa82 197 Int_t pdg = 0;
faee3b18 198 if(fUpperSigmaCut || fLowerSigmaCut){
3a72645a 199 pdg = CutSigmaModel(track->GetRecTrack(), anatype) ? 11 : 0;
faee3b18 200 } else {
201 // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
202 Float_t p = 0.;
3a72645a 203 if(HasAsymmetricSigmaCut() && (p = track->GetRecTrack()->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
faee3b18 204 if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
205 } else {
206 if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
207 }
70da6c5a 208 }
3a72645a 209 if(pidqa && pdg != 0) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
0792aa82 210 return pdg;
722347d8 211
809a4336 212}
213
faee3b18 214//___________________________________________________________________
6555e2ad 215Bool_t AliHFEpidTPC::CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
faee3b18 216 //
217 // N SigmaCut using parametrization of the cuts
218 //
219 Bool_t isSelected = kTRUE;
3a72645a 220 Float_t nsigma = NumberOfSigmas(track, AliPID::kElectron, anaType);
6555e2ad 221 Double_t p = GetP(track, anaType);
faee3b18 222 if(fUpperSigmaCut && nsigma > fUpperSigmaCut->Eval(p)) isSelected = kFALSE;
223 if(fLowerSigmaCut && nsigma < fLowerSigmaCut->Eval(p)) isSelected = kFALSE;
224 return isSelected;
225}
226
0792aa82 227//___________________________________________________________________
6555e2ad 228Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
0792aa82 229 //
230 // reject particles based on asymmetric sigma cut
231 //
232 Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
6555e2ad 233 Double_t p = GetP(track, anaType);
0792aa82 234 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
235 if(!TESTBIT(fRejectionEnabled, ispec)) continue;
0792aa82 236 // Particle rejection enabled
237 if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
3a72645a 238 Double_t sigma = NumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec), anaType);
0792aa82 239 if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
240 }
241 return 0;
242}
243
809a4336 244//___________________________________________________________________
75d81601 245void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
809a4336 246 //
247 // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
248 // Stores line center and line sigma
249 //
250 if(species >= AliPID::kSPECIES){
251 AliError("Species doesn't exist");
252 return;
253 }
254 fLineCrossingsEnabled |= 1 << species;
75d81601 255 fLineCrossingSigma[species] = sigma;
809a4336 256}
257
258//___________________________________________________________________
6555e2ad 259Double_t AliHFEpidTPC::NumberOfSigmas(const AliVParticle *track, AliPID::EParticleType species, AliHFEpidObject::AnalysisType_t anaType) const {
3a72645a 260 //
261 // Get the number of sigmas
809a4336 262 //
3a72645a 263 Double_t nSigmas = 100;
264 if(anaType == AliHFEpidObject::kESDanalysis){
265 // ESD analysis
266 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
267 if(esdtrack && fESDpid) nSigmas = fESDpid->NumberOfSigmasTPC(esdtrack, species);
268 } else {
269 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
270 if(aodtrack && fAODpid) nSigmas = fAODpid->NumberOfSigmasTPC(aodtrack, species);
809a4336 271 }
bf892a6a 272 // Correct for the mean o
273 if(fElectronMeanCorrection)
274 nSigmas -= fElectronMeanCorrection->Eval(GetP(track, anaType));
3a72645a 275 return nSigmas;
809a4336 276}
6555e2ad 277
278//___________________________________________________________________
279Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
280 //
281 // Get the momentum at the inner wall of the TPC
282 //
283 Double_t p = -1;
284 if(anatype == AliHFEpidObject::kESDanalysis){
285 // ESD analysis: Use Inner Params for the momentum estimate
286 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
bf892a6a 287 if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
6555e2ad 288 } else {
289 // AOD analysis: Use TPC momentum stored in the AliAODpid object
290 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
bf892a6a 291 if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
6555e2ad 292 }
293 return p;
294}