Coding rules (Markus)
[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 <TH2I.h>
29 #include <TList.h>
30 #include <TMath.h>
31 //#include <TParticle.h>
32
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 "AliTPCpidESD.h"
41 //#include "AliVParticle.h"
42
43 #include "AliHFEpidTPC.h"
44
45
46
47 //___________________________________________________________________
48 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
49   // add a list here
50     AliHFEpidBase(name)
51   , fLineCrossingsEnabled(0)
52   , fNsigmaTPC(3)
53   , fRejectionEnabled(0)
54   , fPID(NULL)
55   , fPIDtpcESD(NULL)
56   , fQAList(NULL)
57 {
58   //
59   // default  constructor
60   // 
61   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
62   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
63   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
64   fPID = new AliPID;
65   fPIDtpcESD = new AliTPCpidESD;
66 }
67
68 //___________________________________________________________________
69 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
70     AliHFEpidBase("")
71   , fLineCrossingsEnabled(0)
72   , fNsigmaTPC(2)
73   , fRejectionEnabled(0)
74   , fPID(NULL)
75   , fPIDtpcESD(NULL)
76   , fQAList(NULL)
77 {
78   //
79   // Copy constructor
80   //
81   ref.Copy(*this);
82 }
83
84 //___________________________________________________________________
85 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
86   //
87   // Assignment operator
88   //
89   if(this != &ref){
90     ref.Copy(*this);
91   } 
92   return *this;
93 }
94
95 //___________________________________________________________________
96 void AliHFEpidTPC::Copy(TObject &o) const{
97   //
98   // Copy function 
99   // called in copy constructor and assigment operator
100   //
101   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
102
103   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
104   target.fNsigmaTPC = fNsigmaTPC;
105   target.fRejectionEnabled = fRejectionEnabled;
106   target.fPID = new AliPID(*fPID);
107   target.fPIDtpcESD = new AliTPCpidESD(*fPIDtpcESD);
108   target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
109   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
110   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
111   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
112  
113   AliHFEpidBase::Copy(target);
114 }
115
116 //___________________________________________________________________
117 AliHFEpidTPC::~AliHFEpidTPC(){
118   //
119   // Destructor
120   //
121   if(fPID) delete fPID;
122   if(fPIDtpcESD) delete fPIDtpcESD;
123   if(fQAList){
124     fQAList->Delete();
125     delete fQAList;
126   }
127 }
128
129 //___________________________________________________________________
130 Bool_t AliHFEpidTPC::InitializePID(){
131   //
132   // Add TPC dE/dx Line crossings
133   //
134   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
135   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
136   return kTRUE;
137 }
138
139 //___________________________________________________________________
140 Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
141 {
142   //
143   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
144   // for electrons
145   // exclusion of the crossing points
146   //
147   if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
148     AliESDtrack *esdTrack = NULL;
149     if(!(esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack))) return 0;
150     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track->fMCtrack);
151     return MakePIDesd(esdTrack, mctrack);
152   }else{
153     AliAODTrack *aodtrack = NULL;
154     if(!(aodtrack = dynamic_cast<AliAODTrack *>(track->fRecTrack))) return 0;
155     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
156     return MakePIDaod(aodtrack, aodmctrack);
157   }
158 }
159
160 //___________________________________________________________________
161 Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
162   //
163   //  Doing TPC PID as explained in IsSelected for ESD tracks
164   //
165   if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
166   // exclude crossing points:
167   // Determine the bethe values for each particle species
168   Bool_t isLineCrossing = kFALSE;
169   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
170     if(ispecies == AliPID::kElectron) continue;
171     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
172     if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies) < fLineCrossingSigma[ispecies]){
173       // Point in a line crossing region, no PID possible
174       isLineCrossing = kTRUE;
175       break;
176     }
177   }
178   if(isLineCrossing) return 0;
179
180   // Check particle rejection
181   if(HasParticleRejection()){
182     Int_t reject = Reject(esdTrack);
183     AliDebug(1, Form("PID code from Rejection: %d", reject));
184     if(reject != 0) return reject;
185   }
186   // Check whether distance from the electron line is smaller than n-sigma
187
188   // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
189   Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron);
190   Float_t p = 0.;
191   Int_t pdg = 0;
192   if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){ 
193     if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; 
194   } else {
195     if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
196   }
197   if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
198   return pdg;
199 }
200
201 //___________________________________________________________________
202 Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mctrack*/){
203   AliError("AOD PID not yet implemented");
204   return 0;
205 }
206
207 //___________________________________________________________________
208 Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
209   //
210   // reject particles based on asymmetric sigma cut
211   //
212   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
213   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
214   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
215     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
216     AliDebug(1, Form("Particle Rejection enabled for species %d", ispec));
217     // Particle rejection enabled
218     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
219     Double_t sigma = fPIDtpcESD->GetNumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec));
220     AliDebug(1, Form("Sigma %f, min %f, max %f", sigma, fRejection[4*ispec + 1], fRejection[4*ispec+3]));
221     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
222   }
223   return 0;
224 }
225
226 //___________________________________________________________________
227 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
228   //
229   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
230   // Stores line center and line sigma
231   //
232   if(species >= AliPID::kSPECIES){
233     AliError("Species doesn't exist");
234     return;
235   }
236   fLineCrossingsEnabled |= 1 << species;
237   fLineCrossingSigma[species] = sigma;
238 }
239
240 //___________________________________________________________________
241 Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
242 {
243   //gives probability for track to come from a certain particle species;
244   //no priors considered -> probabilities for equal abundances of all species!
245   //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
246
247   //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
248   
249   if(!track) return -1.;
250   Bool_t outlier = kTRUE;
251   // Check whether distance from the respective particle line is smaller than r sigma
252   for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
253     if(TMath::Abs(fPIDtpcESD->GetNumberOfSigmas(track, (AliPID::EParticleType)hypo)) > rsig)
254       outlier = kTRUE;
255     else {
256             outlier = kFALSE;
257             break;
258           }
259   }
260   if(outlier)
261     return -2.;
262
263   Double_t tpcProb[5];
264
265   track->GetTPCpid(tpcProb);
266
267   return tpcProb[species];
268 }
269
270 //___________________________________________________________________
271 Double_t  AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
272 {
273   //ratio of likelihoods to be whatever species/to be an electron;
274   //as a cross-check for possible particle type suppression compared to electrons
275   const Double_t kVerySmall = 1e-10;
276   if(!track) return -20;
277   if((TMath::Abs(Likelihood(track,species) + 2) < kVerySmall)||(TMath::Abs(Likelihood(track,0) + 2 ) < kVerySmall))
278     return -30;
279   if(TMath::Abs(Likelihood(track,species)) < kVerySmall)
280     return -10;
281   if (TMath::Abs(Likelihood(track,0)) < kVerySmall)
282     return 10.;
283   else
284     return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
285
286
287 }
288
289 //___________________________________________________________________
290 void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track, const AliMCParticle *mctrack){
291   // 
292   // Fill the QA histogtrams
293   //
294   if(!track)
295     return;
296  
297   Double_t tpcSignal = track->GetTPCsignal();
298   Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
299   if(HasMCData()){
300     switch(TMath::Abs(mctrack->Particle()->GetPdgCode())){
301       case 11:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
302                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
303                         //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
304                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
305                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
306                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
307                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
308                         //___________________________________________________________________________________________
309                         //Likelihoods for electrons to be other particle species
310                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
311                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
312                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
313                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
314                         break;
315             //___________________________________________________________________________________________
316       case 13:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
317                   //Likelihood of muon to be an electron
318                   (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
319                   //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
320                   //below functions are the same for other species
321                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
322                   break;
323       case 211:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
324                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
325                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
326                   break;
327       case 321:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
328                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
329                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
330                   break;
331       case 2212:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
332                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
333                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
334                   break;
335       default:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
336                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
337
338                         break;
339     }
340   }
341   //TPC signal and Likelihood to be electron for all tracks (independent of MC information)
342   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
343   (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
344 }
345
346 //___________________________________________________________________
347 void AliHFEpidTPC::AddQAhistograms(TList *qaList){
348   //
349   // Create QA histograms for TPC PID
350   //
351   fQAList = new TList;
352   fQAList->SetName("fTPCqaHistos");
353
354   fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 200, 0, 200), kHistTPCelectron); 
355   fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 200, 0, 200), kHistTPCmuon);
356   fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 200, 0, 200), kHistTPCpion);
357   fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 200, 0, 200), kHistTPCkaon);
358   fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 200, 0, 200), kHistTPCproton);
359   fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 200, 0, 200), kHistTPCothers);
360   fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 200, 0, 200), kHistTPCall);
361   fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 200, 0, 200), kHistTPCselected);
362
363   fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
364   fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
365   fQAList->AddAt(new TH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobMu);
366   fQAList->AddAt(new TH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobKa);
367   fQAList->AddAt(new TH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPro);
368   fQAList->AddAt(new TH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobOth);
369   fQAList->AddAt(new TH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobAll);
370
371
372  fQAList->AddAt(new TH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
373  fQAList->AddAt(new TH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
374  fQAList->AddAt(new TH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
375  fQAList->AddAt(new TH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
376
377  fQAList->AddAt(new TH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
378  fQAList->AddAt(new TH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
379  fQAList->AddAt(new TH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
380  fQAList->AddAt(new TH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
381
382  fQAList->AddAt(new TH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPi);
383  fQAList->AddAt(new TH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobMu);
384  fQAList->AddAt(new TH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobKa);
385  fQAList->AddAt(new TH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPro);
386
387
388   qaList->AddLast(fQAList);
389 }