2dd3b52d790401b001e94ad25b1c08fa5c3984b7
[u/mrichter/AliRoot.git] / PWGHF / 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 "AliAODPid.h"
32 #include "AliAODTrack.h"
33 #include "AliAODMCParticle.h"
34 #include "AliESDtrack.h"
35 #include "AliExternalTrackParam.h"
36 #include "AliLog.h"
37 #include "AliMCParticle.h"
38 #include "AliPID.h"
39 #include "AliPIDResponse.h"
40
41 #include "AliHFEpidTPC.h"
42 #include "AliHFEpidQAmanager.h"
43
44 ClassImp(AliHFEpidTPC)
45
46 //___________________________________________________________________
47 AliHFEpidTPC::AliHFEpidTPC() :
48   // add a list here
49   AliHFEpidBase()
50   , fLineCrossingsEnabled(0)
51   , fkEtaCorrection(NULL)
52   , fHasCutModel(kFALSE)
53   , fNsigmaTPC(3)
54   , fRejectionEnabled(0)
55 {
56   //
57   // default  constructor
58   // 
59
60   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
61   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
62
63   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
64   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
65   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
66   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
67
68 }
69
70 //___________________________________________________________________
71 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
72   // add a list here
73   AliHFEpidBase(name)
74   , fLineCrossingsEnabled(0)
75   , fkEtaCorrection(NULL)
76   , fHasCutModel(kFALSE)
77   , fNsigmaTPC(3)
78   , fRejectionEnabled(0)
79 {
80   //
81   // default  constructor
82   // 
83   //
84   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
85   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
86
87   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
88   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
89   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
90   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
91 }
92
93 //___________________________________________________________________
94 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
95   AliHFEpidBase("")
96   , fLineCrossingsEnabled(0)
97   , fkEtaCorrection(NULL)
98   , fHasCutModel(ref.fHasCutModel)
99   , fNsigmaTPC(2)
100   , fRejectionEnabled(0)
101 {
102   //
103   // Copy constructor
104   //
105   ref.Copy(*this);
106 }
107
108 //___________________________________________________________________
109 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
110   //
111   // Assignment operator
112   //
113   if(this != &ref){
114     ref.Copy(*this);
115   } 
116   return *this;
117 }
118 //___________________________________________________________________
119 void AliHFEpidTPC::Copy(TObject &o) const{
120   //
121   // Copy function 
122   // called in copy constructor and assigment operator
123   //
124   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
125
126   target.fkEtaCorrection =fkEtaCorrection;
127   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
128   target.fHasCutModel = fHasCutModel;
129   target.fNsigmaTPC = fNsigmaTPC;
130   target.fRejectionEnabled = fRejectionEnabled;
131
132   memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
133   memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
134
135   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
136   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
137   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
138  
139   AliHFEpidBase::Copy(target);
140 }
141
142 //___________________________________________________________________
143 AliHFEpidTPC::~AliHFEpidTPC(){
144   //
145   // Destructor
146   //
147 }
148
149 //___________________________________________________________________
150 Bool_t AliHFEpidTPC::InitializePID(Int_t /*run*/){
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(!fkPIDResponse) return 0;
169
170   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
171
172   // Make clone track in order to be able to apply correction
173   const AliVTrack *rectrack;
174   AliESDtrack esdtrack;
175   AliAODTrack aodtrack;
176   if(fkEtaCorrection){
177     // Correction available
178     // apply it on copy
179     if(track->IsESDanalysis()){
180       esdtrack.~AliESDtrack();
181       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
182       ApplyEtaCorrection(&esdtrack, anatype);
183       rectrack = &esdtrack;
184     } else {
185       aodtrack.~AliAODTrack();
186       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
187       ApplyEtaCorrection(&aodtrack, anatype);
188       rectrack = &aodtrack;
189     }
190   } else {
191     // Correction not available - no need to copy
192     rectrack = track->GetRecTrack();
193   }
194   AliHFEpidObject tpctrack(*track);
195   tpctrack.SetRecTrack(rectrack);
196   
197   // QA before selection (after correction)
198   if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
199   AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
200   
201   // make copy of the track in order to allow for applying the correction 
202   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron);
203   AliDebug(1, Form("TPC NSigma: %f", nsigma));
204   // exclude crossing points:
205   // Determine the bethe values for each particle species
206   Bool_t isLineCrossing = kFALSE;
207   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
208     if(ispecies == AliPID::kElectron) continue;
209     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
210     if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
211       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
212       isLineCrossing = kTRUE;      
213       break;
214     }
215   }
216   if(isLineCrossing) return 0;
217
218   // Check particle rejection
219   if(HasParticleRejection()){
220     Int_t reject = Reject(rectrack, anatype);
221     if(reject != 0) return reject;
222   }
223
224   // Check if we have an asymmetric sigma model set
225   Int_t pdg = 0;
226   if(fHasCutModel){
227     pdg = CutSigmaModel(&tpctrack) ? 11 : 0;
228   } else { 
229     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
230     Float_t p = rectrack->P();
231     if(HasAsymmetricSigmaCut()) {
232       
233       //printf("p %f, fPAsigCut[0] %f, fPAsigCut[1] %f\n",p,fPAsigCut[0],fPAsigCut[1]);
234       if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) { 
235         if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; 
236         else pdg = 0;
237       }
238       else pdg = 0;
239     
240     }
241     else {
242       if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
243     }
244   }
245   if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
246   return pdg;
247
248 }
249
250 //___________________________________________________________________
251 Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
252   //
253   // N SigmaCut using parametrization of the cuts
254   //
255   Bool_t isSelected = kTRUE;
256   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
257   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron);
258   Double_t p = GetP(track->GetRecTrack(), anatype);
259   Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
260   AliDebug(2, Form("Centrality: %d\n", centrality));
261   if(centrality > 11) return kFALSE;
262   const TF1 *cutfunction;
263   if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
264   if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
265   return isSelected;
266 }
267
268 //___________________________________________________________________
269 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
270   //
271   // reject particles based on asymmetric sigma cut
272   //
273   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
274   Double_t p = GetP(track, anaType);
275   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
276     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
277     // Particle rejection enabled
278     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
279     Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
280     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
281   }
282   return 0;
283 }
284
285 //___________________________________________________________________
286 void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
287   //
288   // Apply correction for the eta dependence
289   // N.B. This correction has to be applied on a copy track
290   //
291   AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName()));
292   Double_t original = track->GetTPCsignal(),
293            eta = track->Eta();
294   if(anatype == AliHFEpidObject::kESDanalysis){
295     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
296     esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
297   } else {
298     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
299     AliAODPid *pid = aodtrack->GetDetPid();
300     if(pid) pid->SetTPCsignal(original);
301   }
302 }
303
304 //___________________________________________________________________
305 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
306   //
307   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
308   // Stores line center and line sigma
309   //
310   if(species >= AliPID::kSPECIES){
311     AliError("Species doesn't exist");
312     return;
313   }
314   fLineCrossingsEnabled |= 1 << species;
315   fLineCrossingSigma[species] = sigma;
316 }
317
318 //___________________________________________________________________
319 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
320   //
321   // Get the momentum at the inner wall of the TPC
322   //
323   Double_t p = -1;
324   if(anatype == AliHFEpidObject::kESDanalysis){
325     // ESD analysis: Use Inner Params for the momentum estimate
326     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
327     if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
328   } else { 
329     // AOD analysis: Use TPC momentum stored in the AliAODpid object
330     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
331     if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
332   }
333   return p;
334 }