]>
Commit | Line | Data |
---|---|---|
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 | ||
6555e2ad | 31 | #include "AliAODPid.h" |
722347d8 | 32 | #include "AliAODTrack.h" |
33 | #include "AliAODMCParticle.h" | |
809a4336 | 34 | #include "AliESDtrack.h" |
35 | #include "AliExternalTrackParam.h" | |
36 | #include "AliLog.h" | |
722347d8 | 37 | #include "AliMCParticle.h" |
809a4336 | 38 | #include "AliPID.h" |
8c1c76e9 | 39 | #include "AliPIDResponse.h" |
809a4336 | 40 | |
41 | #include "AliHFEpidTPC.h" | |
3a72645a | 42 | #include "AliHFEpidQAmanager.h" |
809a4336 | 43 | |
faee3b18 | 44 | ClassImp(AliHFEpidTPC) |
809a4336 | 45 | |
3a72645a | 46 | //___________________________________________________________________ |
47 | AliHFEpidTPC::AliHFEpidTPC() : | |
48 | // add a list here | |
49 | AliHFEpidBase() | |
3a72645a | 50 | , fLineCrossingsEnabled(0) |
11ff28c5 | 51 | , fkEtaCorrection(NULL) |
e156c3bb | 52 | , fHasCutModel(kFALSE) |
3a72645a | 53 | , fNsigmaTPC(3) |
54 | , fRejectionEnabled(0) | |
3a72645a | 55 | { |
56 | // | |
57 | // default constructor | |
58 | // | |
e156c3bb | 59 | |
60 | memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12); | |
61 | memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12); | |
62 | ||
bf892a6a | 63 | memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES); |
64 | memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES); | |
65 | memset(fPAsigCut, 0, sizeof(Float_t) * 2); | |
66 | memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2); | |
67 | ||
3a72645a | 68 | } |
69 | ||
809a4336 | 70 | //___________________________________________________________________ |
71 | AliHFEpidTPC::AliHFEpidTPC(const char* name) : | |
72 | // add a list here | |
9bcfd1ab | 73 | AliHFEpidBase(name) |
809a4336 | 74 | , fLineCrossingsEnabled(0) |
11ff28c5 | 75 | , fkEtaCorrection(NULL) |
e156c3bb | 76 | , fHasCutModel(kFALSE) |
75d81601 | 77 | , fNsigmaTPC(3) |
0792aa82 | 78 | , fRejectionEnabled(0) |
809a4336 | 79 | { |
80 | // | |
81 | // default constructor | |
82 | // | |
bf892a6a | 83 | // |
e156c3bb | 84 | memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12); |
85 | memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12); | |
86 | ||
bf892a6a | 87 | memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES); |
809a4336 | 88 | memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES); |
722347d8 | 89 | memset(fPAsigCut, 0, sizeof(Float_t) * 2); |
90 | memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2); | |
809a4336 | 91 | } |
92 | ||
93 | //___________________________________________________________________ | |
94 | AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) : | |
9bcfd1ab | 95 | AliHFEpidBase("") |
809a4336 | 96 | , fLineCrossingsEnabled(0) |
11ff28c5 | 97 | , fkEtaCorrection(NULL) |
e156c3bb | 98 | , fHasCutModel(ref.fHasCutModel) |
809a4336 | 99 | , fNsigmaTPC(2) |
0792aa82 | 100 | , fRejectionEnabled(0) |
809a4336 | 101 | { |
102 | // | |
103 | // Copy constructor | |
104 | // | |
105 | ref.Copy(*this); | |
106 | } | |
107 | ||
108 | //___________________________________________________________________ | |
109 | AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){ | |
110 | // | |
111 | // Assignment operator | |
112 | // | |
113 | if(this != &ref){ | |
114 | ref.Copy(*this); | |
115 | } | |
116 | return *this; | |
117 | } | |
75d81601 | 118 | //___________________________________________________________________ |
809a4336 | 119 | void AliHFEpidTPC::Copy(TObject &o) const{ |
120 | // | |
121 | // Copy function | |
122 | // called in copy constructor and assigment operator | |
123 | // | |
124 | AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o); | |
125 | ||
11ff28c5 | 126 | target.fkEtaCorrection =fkEtaCorrection; |
809a4336 | 127 | target.fLineCrossingsEnabled = fLineCrossingsEnabled; |
e156c3bb | 128 | target.fHasCutModel = fHasCutModel; |
809a4336 | 129 | target.fNsigmaTPC = fNsigmaTPC; |
0792aa82 | 130 | target.fRejectionEnabled = fRejectionEnabled; |
e156c3bb | 131 | |
132 | memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12); | |
133 | memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12); | |
134 | ||
75d81601 | 135 | memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES); |
722347d8 | 136 | memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2); |
137 | memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2); | |
138 | ||
809a4336 | 139 | AliHFEpidBase::Copy(target); |
140 | } | |
141 | ||
142 | //___________________________________________________________________ | |
143 | AliHFEpidTPC::~AliHFEpidTPC(){ | |
144 | // | |
145 | // Destructor | |
146 | // | |
809a4336 | 147 | } |
148 | ||
149 | //___________________________________________________________________ | |
8c1c76e9 | 150 | Bool_t AliHFEpidTPC::InitializePID(Int_t /*run*/){ |
809a4336 | 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 | 160 | Int_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 | |
8c1c76e9 | 168 | if(!fkPIDResponse) return 0; |
11ff28c5 | 169 | |
170 | AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis; | |
171 | ||
172 | // Make clone track in order to be able to apply correction | |
173 | const AliVTrack *rectrack; | |
174 | AliESDtrack esdtrack; | |
175 | AliAODTrack aodtrack; | |
176 | if(fkEtaCorrection){ | |
177 | // Correction available | |
178 | // apply it on copy | |
179 | if(track->IsESDanalysis()){ | |
180 | esdtrack.~AliESDtrack(); | |
181 | new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack()))); | |
182 | ApplyEtaCorrection(&esdtrack, anatype); | |
183 | rectrack = &esdtrack; | |
184 | } else { | |
185 | aodtrack.~AliAODTrack(); | |
186 | new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack()))); | |
187 | ApplyEtaCorrection(&aodtrack, anatype); | |
188 | rectrack = &aodtrack; | |
189 | } | |
190 | } else { | |
191 | // Correction not available - no need to copy | |
192 | rectrack = track->GetRecTrack(); | |
193 | } | |
194 | AliHFEpidObject tpctrack(*track); | |
195 | tpctrack.SetRecTrack(rectrack); | |
3a72645a | 196 | |
11ff28c5 | 197 | // QA before selection (after correction) |
198 | if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID); | |
e3fc062d | 199 | AliDebug(1, "Doing TPC PID based on n-Sigma cut approach"); |
11ff28c5 | 200 | |
201 | // make copy of the track in order to allow for applying the correction | |
202 | Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron); | |
3a72645a | 203 | AliDebug(1, Form("TPC NSigma: %f", nsigma)); |
809a4336 | 204 | // exclude crossing points: |
205 | // Determine the bethe values for each particle species | |
809a4336 | 206 | Bool_t isLineCrossing = kFALSE; |
207 | for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){ | |
75d81601 | 208 | if(ispecies == AliPID::kElectron) continue; |
809a4336 | 209 | if(!(fLineCrossingsEnabled & 1 << ispecies)) continue; |
11ff28c5 | 210 | if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){ |
9bcfd1ab | 211 | // Point in a line crossing region, no PID possible, but !PID still possible ;-) |
70da6c5a | 212 | isLineCrossing = kTRUE; |
809a4336 | 213 | break; |
214 | } | |
215 | } | |
216 | if(isLineCrossing) return 0; | |
0792aa82 | 217 | |
218 | // Check particle rejection | |
219 | if(HasParticleRejection()){ | |
11ff28c5 | 220 | Int_t reject = Reject(rectrack, anatype); |
0792aa82 | 221 | if(reject != 0) return reject; |
222 | } | |
722347d8 | 223 | |
faee3b18 | 224 | // Check if we have an asymmetric sigma model set |
0792aa82 | 225 | Int_t pdg = 0; |
e156c3bb | 226 | if(fHasCutModel){ |
11ff28c5 | 227 | pdg = CutSigmaModel(&tpctrack) ? 11 : 0; |
faee3b18 | 228 | } else { |
229 | // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut | |
5cd679b7 | 230 | Float_t p = rectrack->P(); |
231 | if(HasAsymmetricSigmaCut()) { | |
232 | ||
233 | //printf("p %f, fPAsigCut[0] %f, fPAsigCut[1] %f\n",p,fPAsigCut[0],fPAsigCut[1]); | |
234 | if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) { | |
235 | if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; | |
236 | else pdg = 0; | |
237 | } | |
238 | else pdg = 0; | |
239 | ||
240 | } | |
241 | else { | |
faee3b18 | 242 | if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11; |
243 | } | |
70da6c5a | 244 | } |
11ff28c5 | 245 | if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID); |
0792aa82 | 246 | return pdg; |
722347d8 | 247 | |
809a4336 | 248 | } |
249 | ||
faee3b18 | 250 | //___________________________________________________________________ |
e156c3bb | 251 | Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const { |
faee3b18 | 252 | // |
253 | // N SigmaCut using parametrization of the cuts | |
254 | // | |
255 | Bool_t isSelected = kTRUE; | |
e156c3bb | 256 | AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis; |
8c1c76e9 | 257 | Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron); |
e156c3bb | 258 | Double_t p = GetP(track->GetRecTrack(), anatype); |
259 | Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0; | |
260 | AliDebug(2, Form("Centrality: %d\n", centrality)); | |
8c1c76e9 | 261 | if(centrality > 11) return kFALSE; |
e156c3bb | 262 | const TF1 *cutfunction; |
263 | if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE; | |
264 | if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE; | |
faee3b18 | 265 | return isSelected; |
266 | } | |
267 | ||
0792aa82 | 268 | //___________________________________________________________________ |
6555e2ad | 269 | Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{ |
0792aa82 | 270 | // |
271 | // reject particles based on asymmetric sigma cut | |
272 | // | |
273 | Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212}; | |
6555e2ad | 274 | Double_t p = GetP(track, anaType); |
0792aa82 | 275 | for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){ |
276 | if(!TESTBIT(fRejectionEnabled, ispec)) continue; | |
0792aa82 | 277 | // Particle rejection enabled |
278 | if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue; | |
8c1c76e9 | 279 | Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec)); |
0792aa82 | 280 | if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge(); |
281 | } | |
282 | return 0; | |
283 | } | |
284 | ||
11ff28c5 | 285 | //___________________________________________________________________ |
286 | void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{ | |
287 | // | |
288 | // Apply correction for the eta dependence | |
289 | // N.B. This correction has to be applied on a copy track | |
290 | // | |
291 | AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName())); | |
292 | Double_t original = track->GetTPCsignal(), | |
293 | eta = track->Eta(); | |
294 | if(anatype == AliHFEpidObject::kESDanalysis){ | |
295 | AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track); | |
296 | esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN()); | |
297 | } else { | |
298 | AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track); | |
299 | AliAODPid *pid = aodtrack->GetDetPid(); | |
300 | if(pid) pid->SetTPCsignal(original); | |
301 | } | |
302 | } | |
303 | ||
809a4336 | 304 | //___________________________________________________________________ |
75d81601 | 305 | void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){ |
809a4336 | 306 | // |
307 | // Add exclusion point for the TPC PID where a dEdx line crosses the electron line | |
308 | // Stores line center and line sigma | |
309 | // | |
310 | if(species >= AliPID::kSPECIES){ | |
311 | AliError("Species doesn't exist"); | |
312 | return; | |
313 | } | |
314 | fLineCrossingsEnabled |= 1 << species; | |
75d81601 | 315 | fLineCrossingSigma[species] = sigma; |
809a4336 | 316 | } |
317 | ||
6555e2ad | 318 | //___________________________________________________________________ |
319 | Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const { | |
320 | // | |
321 | // Get the momentum at the inner wall of the TPC | |
322 | // | |
323 | Double_t p = -1; | |
324 | if(anatype == AliHFEpidObject::kESDanalysis){ | |
325 | // ESD analysis: Use Inner Params for the momentum estimate | |
326 | const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track); | |
bf892a6a | 327 | if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P(); |
6555e2ad | 328 | } else { |
329 | // AOD analysis: Use TPC momentum stored in the AliAODpid object | |
330 | const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track); | |
bf892a6a | 331 | if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P(); |
6555e2ad | 332 | } |
333 | return p; | |
334 | } |