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 "AliTPCdEdxInfo.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 "AliPIDResponse.h"
42 #include "AliHFEpidTPC.h"
43 #include "AliHFEpidQAmanager.h"
45 ClassImp(AliHFEpidTPC)
47 //___________________________________________________________________
48 AliHFEpidTPC::AliHFEpidTPC() :
51 , fLineCrossingsEnabled(0)
52 , fkEtaCorrection(NULL)
53 , fkCentralityCorrection(NULL)
54 , fHasCutModel(kFALSE)
55 , fUseOnlyOROC(kFALSE)
57 , fRejectionEnabled(0)
61 // default constructor
64 memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
65 memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
67 memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
68 memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
69 memset(fPAsigCut, 0, sizeof(Float_t) * 2);
70 memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
74 //___________________________________________________________________
75 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
78 , fLineCrossingsEnabled(0)
79 , fkEtaCorrection(NULL)
80 , fkCentralityCorrection(NULL)
81 , fHasCutModel(kFALSE)
82 , fUseOnlyOROC(kFALSE)
84 , fRejectionEnabled(0)
88 // default constructor
91 memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
92 memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
94 memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
95 memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
96 memset(fPAsigCut, 0, sizeof(Float_t) * 2);
97 memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
100 //___________________________________________________________________
101 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
103 , fLineCrossingsEnabled(0)
104 , fkEtaCorrection(NULL)
105 , fkCentralityCorrection(NULL)
106 , fHasCutModel(ref.fHasCutModel)
107 , fUseOnlyOROC(ref.fUseOnlyOROC)
109 , fRejectionEnabled(0)
118 //___________________________________________________________________
119 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
121 // Assignment operator
128 //___________________________________________________________________
129 void AliHFEpidTPC::Copy(TObject &o) const{
132 // called in copy constructor and assigment operator
134 AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
136 target.fkEtaCorrection = fkEtaCorrection;
137 target.fkCentralityCorrection = fkCentralityCorrection;
138 target.fLineCrossingsEnabled = fLineCrossingsEnabled;
139 target.fHasCutModel = fHasCutModel;
140 target.fUseOnlyOROC = fUseOnlyOROC;
141 target.fNsigmaTPC = fNsigmaTPC;
142 target.fRejectionEnabled = fRejectionEnabled;
144 memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
145 memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
147 memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
148 memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
149 memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
151 AliHFEpidBase::Copy(target);
154 //___________________________________________________________________
155 AliHFEpidTPC::~AliHFEpidTPC(){
161 //___________________________________________________________________
162 Bool_t AliHFEpidTPC::InitializePID(Int_t /*run*/){
164 // Add TPC dE/dx Line crossings
166 //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
167 //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
171 //___________________________________________________________________
172 Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
175 // For the TPC pid we use the 2-sigma band around the bethe bloch curve
177 // exclusion of the crossing points
180 if(!fkPIDResponse) return 0;
182 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
184 // Make clone track in order to be able to apply correction
185 const AliVTrack *rectrack;
186 AliESDtrack esdtrack;
187 AliAODTrack aodtrack;
188 /*if(fUseOnlyOROC && !(fkEtaCorrection || fkCentralityCorrection)) {
189 if(track->IsESDanalysis()){
190 esdtrack.~AliESDtrack();
191 new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
192 UseOROC(&esdtrack, anatype);
193 rectrack = &esdtrack;
195 aodtrack.~AliAODTrack();
196 new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
197 UseOROC(&aodtrack, anatype);
198 rectrack = &aodtrack;
201 else if(fkEtaCorrection || fkCentralityCorrection){*/
202 if(fkEtaCorrection || fkCentralityCorrection){
203 // Correction available
205 if(track->IsESDanalysis()){
206 esdtrack.~AliESDtrack();
207 new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
208 if(track->IsPbPb() && HasCentralityCorrection())
209 ApplyCentralityCorrection(&esdtrack, track->GetMultiplicity(), anatype);
210 if(HasEtaCorrection())
211 ApplyEtaCorrection(&esdtrack, anatype);
212 rectrack = &esdtrack;
214 aodtrack.~AliAODTrack();
215 new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
216 if(track->IsPbPb() && HasCentralityCorrection())
217 ApplyCentralityCorrection(&aodtrack, track->GetMultiplicity(), anatype);
218 if(HasEtaCorrection())
219 ApplyEtaCorrection(&aodtrack, anatype);
220 rectrack = &aodtrack;
223 // Correction not available - no need to copy
224 rectrack = track->GetRecTrack();
226 AliHFEpidObject tpctrack(*track);
227 tpctrack.SetRecTrack(rectrack);
229 // QA before selection (after correction)
230 if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
231 AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
233 // make copy of the track in order to allow for applying the correction
234 Float_t nsigma = fUsedEdx ? rectrack->GetTPCsignal() : fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron);
235 AliDebug(1, Form("TPC NSigma: %f", nsigma));
236 // exclude crossing points:
237 // Determine the bethe values for each particle species
238 Bool_t isLineCrossing = kFALSE;
239 for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
240 if(ispecies == AliPID::kElectron) continue;
241 if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
242 if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
243 // Point in a line crossing region, no PID possible, but !PID still possible ;-)
244 isLineCrossing = kTRUE;
248 if(isLineCrossing) return 0;
250 // Check particle rejection
251 if(HasParticleRejection()){
252 Int_t reject = Reject(rectrack, anatype);
253 if(reject != 0) return reject;
256 // Check if we have an asymmetric sigma model set
259 pdg = CutSigmaModel(&tpctrack) ? 11 : 0;
261 // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
262 Float_t p = rectrack->P();
263 if(HasAsymmetricSigmaCut()) {
265 //printf("p %f, fPAsigCut[0] %f, fPAsigCut[1] %f\n",p,fPAsigCut[0],fPAsigCut[1]);
266 if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) {
267 if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
274 if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
277 if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
282 //___________________________________________________________________
283 Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
285 // N SigmaCut using parametrization of the cuts
287 Bool_t isSelected = kTRUE;
288 AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
289 Float_t nsigma = fUsedEdx ? track->GetRecTrack()->GetTPCsignal() : fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron);
290 Double_t p = GetP(track->GetRecTrack(), anatype);
291 Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
292 AliDebug(2, Form("Centrality: %d\n", centrality));
293 if(centrality > 11) return kFALSE;
294 const TF1 *cutfunction;
295 if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
296 if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
300 //___________________________________________________________________
301 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
303 // reject particles based on asymmetric sigma cut
305 Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
306 Double_t p = GetP(track, anaType);
307 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
308 if(!TESTBIT(fRejectionEnabled, ispec)) continue;
309 // Particle rejection enabled
310 if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
311 Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
312 if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
317 //___________________________________________________________________
318 void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
320 // Apply correction for the eta dependence
321 // N.B. This correction has to be applied on a copy track
323 AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName()));
324 Double_t original = track->GetTPCsignal(),
326 if(anatype == AliHFEpidObject::kESDanalysis){
327 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
328 if(fkEtaCorrection->Eval(eta)>0.0) esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
330 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
331 AliAODPid *pid = aodtrack->GetDetPid();
332 if(pid && fkEtaCorrection->Eval(eta)>0.0) pid->SetTPCsignal(original/fkEtaCorrection->Eval(eta));
336 //___________________________________________________________________
337 void AliHFEpidTPC::ApplyCentralityCorrection(AliVTrack *track, Double_t centralityEstimator, AliHFEpidObject::AnalysisType_t anatype) const{
339 // Apply correction for the eta dependence
340 // N.B. This correction has to be applied on a copy track
342 AliDebug(1, Form("Applying correction function %s\n", fkCentralityCorrection->GetName()));
343 Double_t original = track->GetTPCsignal();
344 if(anatype == AliHFEpidObject::kESDanalysis){
345 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
346 if(fkCentralityCorrection->Eval(centralityEstimator)>0.0) esdtrack->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
348 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
349 AliAODPid *pid = aodtrack->GetDetPid();
350 if(pid && fkCentralityCorrection->Eval(centralityEstimator)>0.0) pid->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator));
354 //___________________________________________________________________
355 void AliHFEpidTPC::UseOROC(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
357 // Use TPC signal from the OROC
358 // N.B. This correction has to be applied on a copy track
360 //Double_t original = track->GetTPCsignal();
362 if(anatype == AliHFEpidObject::kESDanalysis){
363 AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
364 AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
365 Double32_t TPCsignalRegion[4]; // TPC dEdx signal in 4 different regions - 0 - IROC, 1- OROC medium, 2 - OROC long, 3- OROC all, (default truncation used)
366 Char_t TPCsignalNRegion[3]; // number of clusters above threshold used in the dEdx calculation
367 Char_t TPCsignalNRowRegion[3]; // number of crosed rows used in the dEdx calculation - signal below threshold included
368 dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
369 esdtrack->SetTPCsignal(TPCsignalRegion[3],esdtrack->GetTPCsignalSigma(),(TPCsignalNRegion[1]+TPCsignalNRegion[2])); // the two last are not ok
371 AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
372 AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
373 Double32_t TPCsignalRegion[4]; // TPC dEdx signal in 4 different regions - 0 - IROC, 1- OROC medium, 2 - OROC long, 3- OROC all, (default truncation used)
374 Char_t TPCsignalNRegion[3]; // number of clusters above threshold used in the dEdx calculation
375 Char_t TPCsignalNRowRegion[3]; // number of crosed rows used in the dEdx calculation - signal below threshold included
376 dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
377 AliAODPid *pid = aodtrack->GetDetPid();
378 if(pid) pid->SetTPCsignal(TPCsignalRegion[3]);
379 if(pid) pid->SetTPCsignalN((TPCsignalNRegion[1]+TPCsignalNRegion[2]));
383 //___________________________________________________________________
384 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
386 // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
387 // Stores line center and line sigma
389 if(species >= AliPID::kSPECIES){
390 AliError("Species doesn't exist");
393 fLineCrossingsEnabled |= 1 << species;
394 fLineCrossingSigma[species] = sigma;
397 //___________________________________________________________________
398 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
400 // Get the momentum at the inner wall of the TPC
403 if(anatype == AliHFEpidObject::kESDanalysis){
404 // ESD analysis: Use Inner Params for the momentum estimate
405 const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
406 if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
408 // AOD analysis: Use TPC momentum stored in the AliAODpid object
409 const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
410 if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();