Since it contains fixes of coding rule violations, all classes are involved. Further...
[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**************************************************************************/
15/************************************************************************
16 * *
17 * Class for TPC PID *
18 * Implements the abstract base class AliHFEpidBase *
19 * *
20 * Class contains TPC specific cuts and QA histograms *
21 * Two selection strategies are offered: Selection of certain value *
22 * regions in the TPC dE/dx (by IsSelected), and likelihoods *
23 * *
24 * Authors: *
25 * *
26 * Markus Fasel <M.Fasel@gsi.de> *
27 * Markus Heide <mheide@uni-muenster.de> *
28 * *
29 * *
30 ************************************************************************/
31#include <TH2I.h>
32#include <TList.h>
33#include <TMath.h>
34
35#include "AliESDtrack.h"
36#include "AliExternalTrackParam.h"
37#include "AliLog.h"
38#include "AliPID.h"
39#include "AliTPCpidESD.h"
40#include "AliVParticle.h"
41
42#include "AliHFEpidTPC.h"
43
44
45
46//___________________________________________________________________
47AliHFEpidTPC::AliHFEpidTPC(const char* name) :
48 // add a list here
49 AliHFEpidBase(name)
50 , fLineCrossingsEnabled(0)
75d81601 51 , fNsigmaTPC(3)
52 , fPID(NULL)
53 , fPIDtpcESD(NULL)
54 , fQAList(NULL)
809a4336 55{
56 //
57 // default constructor
58 //
809a4336 59 memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
60 fPID = new AliPID;
61 fPIDtpcESD = new AliTPCpidESD;
62}
63
64//___________________________________________________________________
65AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
66 AliHFEpidBase("")
67 , fLineCrossingsEnabled(0)
68 , fNsigmaTPC(2)
75d81601 69 , fPID(NULL)
70 , fPIDtpcESD(NULL)
71 , fQAList(NULL)
809a4336 72{
73 //
74 // Copy constructor
75 //
76 ref.Copy(*this);
77}
78
79//___________________________________________________________________
80AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
81 //
82 // Assignment operator
83 //
84 if(this != &ref){
85 ref.Copy(*this);
86 }
87 return *this;
88}
89
75d81601 90//___________________________________________________________________
809a4336 91void AliHFEpidTPC::Copy(TObject &o) const{
92 //
93 // Copy function
94 // called in copy constructor and assigment operator
95 //
96 AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
97
98 target.fLineCrossingsEnabled = fLineCrossingsEnabled;
99 target.fNsigmaTPC = fNsigmaTPC;
100 target.fPID = new AliPID(*fPID);
101 target.fPIDtpcESD = new AliTPCpidESD(*fPIDtpcESD);
102 target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
75d81601 103 memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
809a4336 104
105 AliHFEpidBase::Copy(target);
106}
107
108//___________________________________________________________________
109AliHFEpidTPC::~AliHFEpidTPC(){
110 //
111 // Destructor
112 //
113 if(fPID) delete fPID;
114 if(fPIDtpcESD) delete fPIDtpcESD;
115 if(fQAList){
116 fQAList->Delete();
117 delete fQAList;
118 }
119}
120
121//___________________________________________________________________
122Bool_t AliHFEpidTPC::InitializePID(){
123 //
124 // Add TPC dE/dx Line crossings
125 //
75d81601 126 //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
127 //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
809a4336 128 return kTRUE;
129}
130
131//___________________________________________________________________
132Int_t AliHFEpidTPC::IsSelected(AliVParticle *track)
133{
134 //
135 // For the TPC pid we use the 2-sigma band around the bethe bloch curve
136 // for electrons
137 // exclusion of the crossing points
138 //
75d81601 139 AliESDtrack *esdTrack = NULL;
809a4336 140 if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
141 if(IsQAon()) FillTPChistograms(esdTrack);
809a4336 142 // exclude crossing points:
143 // Determine the bethe values for each particle species
809a4336 144 Bool_t isLineCrossing = kFALSE;
145 for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
75d81601 146 if(ispecies == AliPID::kElectron) continue;
809a4336 147 if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
75d81601 148 if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies) < fLineCrossingSigma[ispecies]){
809a4336 149 // Point in a line crossing region, no PID possible
150 isLineCrossing = kTRUE;
151 break;
152 }
153 }
154 if(isLineCrossing) return 0;
155 // Check whether distance from the electron line is smaller than n-sigma
75d81601 156 if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron) < fNsigmaTPC ) return 11;
809a4336 157 return 0;
158}
159
160//___________________________________________________________________
75d81601 161void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
809a4336 162 //
163 // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
164 // Stores line center and line sigma
165 //
166 if(species >= AliPID::kSPECIES){
167 AliError("Species doesn't exist");
168 return;
169 }
170 fLineCrossingsEnabled |= 1 << species;
75d81601 171 fLineCrossingSigma[species] = sigma;
809a4336 172}
173
174//___________________________________________________________________
175Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
176{
177 //gives probability for track to come from a certain particle species;
178 //no priors considered -> probabilities for equal abundances of all species!
179 //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
180
181 //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
182
183 if(!track) return -1.;
809a4336 184 Bool_t outlier = kTRUE;
185 // Check whether distance from the respective particle line is smaller than r sigma
75d81601 186 for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
187 if(fPIDtpcESD->GetNumberOfSigmas(track, (AliPID::EParticleType)hypo) > rsig)
188 outlier = kTRUE;
189 else {
190 outlier = kFALSE;
191 break;
192 }
193 }
809a4336 194 if(outlier)
195 return -2.;
196
75d81601 197 Double_t tpcProb[5];
809a4336 198
75d81601 199 track->GetTPCpid(tpcProb);
809a4336 200
75d81601 201 return tpcProb[species];
809a4336 202}
809a4336 203
809a4336 204//___________________________________________________________________
205Double_t AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
206{
207 //ratio of likelihoods to be whatever species/to be an electron;
208 //as a cross-check for possible particle type suppression compared to electrons
209 if(!track) return -20;
210 if((Likelihood(track,species) == -2.)||(Likelihood(track,0)== -2.))
211 return -30;
212 if(Likelihood(track,species) == 0.)
213 return -10;
214 if (Likelihood(track,0) == 0.)
215 return 10.;
216 else
217 return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
218
219
220}
75d81601 221
809a4336 222//___________________________________________________________________
223void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track){
224 //
225 if(!track)
226 return;
227
75d81601 228 Double_t tpcSignal = track->GetTPCsignal();
809a4336 229 Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
230 if(HasMCData()){
231 switch(TMath::Abs(GetPdgCode(const_cast<AliESDtrack *>(track)))){
75d81601 232 case 11: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
809a4336 233 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
234 //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
235 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
236 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
237 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
238 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
239 //___________________________________________________________________________________________
240 //Likelihoods for electrons to be other particle species
241 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
242 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
243 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
244 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
245 break;
246 //___________________________________________________________________________________________
75d81601 247 case 13: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
809a4336 248 //Likelihood of muon to be an electron
249 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
250 //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
251 //below functions are the same for other species
252 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
253 break;
75d81601 254 case 211: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
809a4336 255 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
256 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
257 break;
75d81601 258 case 321: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
809a4336 259 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
260 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
261 break;
75d81601 262 case 2212: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
809a4336 263 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
264 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
265 break;
75d81601 266 default: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
809a4336 267 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
268
269 break;
270 }
271 }
272 //TPC signal and Likelihood to be electron for all tracks (independent of MC information)
75d81601 273 (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
809a4336 274 (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
809a4336 275}
276
277//___________________________________________________________________
278void AliHFEpidTPC::AddQAhistograms(TList *qaList){
279 fQAList = new TList;
280 fQAList->SetName("fTPCqaHistos");
281
282 fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron);
283 fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon);
284 fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion);
285 fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon);
286 fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
287 fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
288 fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
289
290 fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
291 fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
292 fQAList->AddAt(new TH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobMu);
293 fQAList->AddAt(new TH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobKa);
294 fQAList->AddAt(new TH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPro);
295 fQAList->AddAt(new TH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobOth);
296 fQAList->AddAt(new TH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobAll);
297
298
299 fQAList->AddAt(new TH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
300 fQAList->AddAt(new TH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
301 fQAList->AddAt(new TH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
302 fQAList->AddAt(new TH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
303
304 fQAList->AddAt(new TH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
305 fQAList->AddAt(new TH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
306 fQAList->AddAt(new TH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
307 fQAList->AddAt(new TH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
308
309 fQAList->AddAt(new TH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPi);
310 fQAList->AddAt(new TH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobMu);
311 fQAList->AddAt(new TH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobKa);
312 fQAList->AddAt(new TH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPro);
313
314
315 qaList->AddLast(fQAList);
316}