Since it contains fixes of coding rule violations, all classes are involved. Further...
[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  *                                                                      *
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 //___________________________________________________________________
47 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
48   // add a list here
49     AliHFEpidBase(name)
50   , fLineCrossingsEnabled(0)
51   , fNsigmaTPC(3)
52   , fPID(NULL)
53   , fPIDtpcESD(NULL)
54   , fQAList(NULL)
55 {
56   //
57   // default  constructor
58   // 
59   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
60   fPID = new AliPID;
61   fPIDtpcESD = new AliTPCpidESD;
62 }
63
64 //___________________________________________________________________
65 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
66     AliHFEpidBase("")
67   , fLineCrossingsEnabled(0)
68   , fNsigmaTPC(2)
69   , fPID(NULL)
70   , fPIDtpcESD(NULL)
71   , fQAList(NULL)
72 {
73   //
74   // Copy constructor
75   //
76   ref.Copy(*this);
77 }
78
79 //___________________________________________________________________
80 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
81   //
82   // Assignment operator
83   //
84   if(this != &ref){
85     ref.Copy(*this);
86   } 
87   return *this;
88 }
89
90 //___________________________________________________________________
91 void 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());
103   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
104
105   AliHFEpidBase::Copy(target);
106 }
107
108 //___________________________________________________________________
109 AliHFEpidTPC::~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 //___________________________________________________________________
122 Bool_t AliHFEpidTPC::InitializePID(){
123   //
124   // Add TPC dE/dx Line crossings
125   //
126   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
127   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
128   return kTRUE;
129 }
130
131 //___________________________________________________________________
132 Int_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   //
139   AliESDtrack *esdTrack = NULL;
140   if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
141   if(IsQAon()) FillTPChistograms(esdTrack);
142   // exclude crossing points:
143   // Determine the bethe values for each particle species
144   Bool_t isLineCrossing = kFALSE;
145   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
146     if(ispecies == AliPID::kElectron) continue;
147     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
148     if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies) < fLineCrossingSigma[ispecies]){
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
156   if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron) < fNsigmaTPC ) return 11;
157   return 0;
158 }
159
160 //___________________________________________________________________
161 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
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;
171   fLineCrossingSigma[species] = sigma;
172 }
173
174 //___________________________________________________________________
175 Double_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.;
184   Bool_t outlier = kTRUE;
185   // Check whether distance from the respective particle line is smaller than r sigma
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   }
194   if(outlier)
195     return -2.;
196
197   Double_t tpcProb[5];
198
199   track->GetTPCpid(tpcProb);
200
201   return tpcProb[species];
202 }
203
204 //___________________________________________________________________
205 Double_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 }
221
222 //___________________________________________________________________
223 void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track){
224   //
225   if(!track)
226     return;
227  
228   Double_t tpcSignal = track->GetTPCsignal();
229   Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
230   if(HasMCData()){
231     switch(TMath::Abs(GetPdgCode(const_cast<AliESDtrack *>(track)))){
232       case 11:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
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         //___________________________________________________________________________________________
247     case 13: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
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;
254       case 211:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
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;
258       case 321:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
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;
262       case 2212: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
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;
266       default: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
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)
273   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
274  (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
275 }
276
277 //___________________________________________________________________
278 void 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 }