326613ad33bf62c2bd851888cff34901d4a54ce3
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpidTPC.cxx
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 **************************************************************************/
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 //  
28 #include <TF1.h>
29 #include <TMath.h>
30
31 #include "AliAODpidUtil.h"
32 #include "AliAODPid.h"
33 #include "AliAODTrack.h"
34 #include "AliAODMCParticle.h"
35 #include "AliESDtrack.h"
36 #include "AliExternalTrackParam.h"
37 #include "AliLog.h"
38 #include "AliMCParticle.h"
39 #include "AliPID.h"
40 #include "AliESDpid.h"
41
42 #include "AliHFEpidTPC.h"
43 #include "AliHFEpidQAmanager.h"
44
45 ClassImp(AliHFEpidTPC)
46
47 //___________________________________________________________________
48 AliHFEpidTPC::AliHFEpidTPC() :
49   // add a list here
50   AliHFEpidBase()
51   , fLineCrossingsEnabled(0)
52   , fUpperSigmaCut(NULL)
53   , fLowerSigmaCut(NULL)
54   , fElectronMeanCorrection(NULL)
55   , fNsigmaTPC(3)
56   , fRejectionEnabled(0)
57   , fPID(NULL)
58 {
59   //
60   // default  constructor
61   // 
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
67 }
68
69 //___________________________________________________________________
70 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
71   // add a list here
72   AliHFEpidBase(name)
73   , fLineCrossingsEnabled(0)
74   , fUpperSigmaCut(NULL)
75   , fLowerSigmaCut(NULL)
76   , fElectronMeanCorrection(NULL)
77   , fNsigmaTPC(3)
78   , fRejectionEnabled(0)
79   , fPID(NULL)
80 {
81   //
82   // default  constructor
83   // 
84   //
85   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
86   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
87   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
88   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
89   fPID = new AliPID;
90 }
91
92 //___________________________________________________________________
93 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
94   AliHFEpidBase("")
95   , fLineCrossingsEnabled(0)
96   , fUpperSigmaCut(NULL)
97   , fLowerSigmaCut(NULL)
98   , fElectronMeanCorrection(NULL)
99   , fNsigmaTPC(2)
100   , fRejectionEnabled(0)
101   , fPID(NULL)
102 {
103   //
104   // Copy constructor
105   //
106   ref.Copy(*this);
107 }
108
109 //___________________________________________________________________
110 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
111   //
112   // Assignment operator
113   //
114   if(this != &ref){
115     ref.Copy(*this);
116   } 
117   return *this;
118 }
119 //___________________________________________________________________
120 void 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;
128   target.fUpperSigmaCut = fUpperSigmaCut;
129   target.fLowerSigmaCut = fLowerSigmaCut;
130   target.fElectronMeanCorrection = fElectronMeanCorrection;
131   target.fNsigmaTPC = fNsigmaTPC;
132   target.fRejectionEnabled = fRejectionEnabled;
133   target.fPID = new AliPID(*fPID);
134   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
135   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
136   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
137  
138   AliHFEpidBase::Copy(target);
139 }
140
141 //___________________________________________________________________
142 AliHFEpidTPC::~AliHFEpidTPC(){
143   //
144   // Destructor
145   //
146   if(fPID) delete fPID;
147 }
148
149 //___________________________________________________________________
150 Bool_t AliHFEpidTPC::InitializePID(){
151   //
152   // Add TPC dE/dx Line crossings
153   //
154   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
155   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
156   return kTRUE;
157 }
158
159 //___________________________________________________________________
160 Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
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   //
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);
172   AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
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));
176   // exclude crossing points:
177   // Determine the bethe values for each particle species
178   Bool_t isLineCrossing = kFALSE;
179   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
180     if(ispecies == AliPID::kElectron) continue;
181     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
182     if(TMath::Abs(NumberOfSigmas(track->GetRecTrack(), (AliPID::EParticleType)ispecies, anatype)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
183       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
184       isLineCrossing = kTRUE;      
185       break;
186     }
187   }
188   if(isLineCrossing) return 0;
189
190   // Check particle rejection
191   if(HasParticleRejection()){
192     Int_t reject = Reject(track->GetRecTrack(), anatype);
193     if(reject != 0) return reject;
194   }
195
196   // Check if we have an asymmetric sigma model set
197   Int_t pdg = 0;
198   if(fUpperSigmaCut || fLowerSigmaCut){
199     pdg = CutSigmaModel(track->GetRecTrack(), anatype) ? 11 : 0;
200   } else { 
201     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
202     Float_t p = 0.;
203     if(HasAsymmetricSigmaCut() && (p = track->GetRecTrack()->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){ 
204       if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; 
205     } else {
206       if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
207     }
208   }
209   if(pidqa && pdg != 0) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
210   return pdg;
211
212 }
213
214 //___________________________________________________________________
215 Bool_t AliHFEpidTPC::CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
216   //
217   // N SigmaCut using parametrization of the cuts
218   //
219   Bool_t isSelected = kTRUE;
220   Float_t nsigma = NumberOfSigmas(track, AliPID::kElectron, anaType);
221   Double_t p = GetP(track, anaType);
222   if(fUpperSigmaCut && nsigma > fUpperSigmaCut->Eval(p)) isSelected = kFALSE;
223   if(fLowerSigmaCut && nsigma < fLowerSigmaCut->Eval(p)) isSelected = kFALSE;
224   return isSelected;
225 }
226
227 //___________________________________________________________________
228 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
229   //
230   // reject particles based on asymmetric sigma cut
231   //
232   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
233   Double_t p = GetP(track, anaType);
234   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
235     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
236     // Particle rejection enabled
237     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
238     Double_t sigma = NumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec), anaType);
239     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
240   }
241   return 0;
242 }
243
244 //___________________________________________________________________
245 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
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;
255   fLineCrossingSigma[species] = sigma;
256 }
257
258 //___________________________________________________________________
259 Double_t AliHFEpidTPC::NumberOfSigmas(const AliVParticle *track, AliPID::EParticleType species, AliHFEpidObject::AnalysisType_t anaType) const {
260   //    
261   // Get the number of sigmas
262   //
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);
271   }
272   // Correct for the mean o
273   if(fElectronMeanCorrection)
274     nSigmas -= fElectronMeanCorrection->Eval(GetP(track, anaType));   
275   return nSigmas;
276 }
277
278 //___________________________________________________________________
279 Double_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);
287     if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
288   } else { 
289     // AOD analysis: Use TPC momentum stored in the AliAODpid object
290     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
291     if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
292   }
293   return p;
294 }