Update hfe package
[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 "AliTPCdEdxInfo.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 "AliPIDResponse.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   , fkEtaCorrection(NULL)
53   , fHasCutModel(kFALSE)
54   , fUseOnlyOROC(kFALSE)
55   , fNsigmaTPC(3)
56   , fRejectionEnabled(0)
57 {
58   //
59   // default  constructor
60   // 
61
62   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
63   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
64
65   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
66   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
67   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
68   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
69
70 }
71
72 //___________________________________________________________________
73 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
74   // add a list here
75   AliHFEpidBase(name)
76   , fLineCrossingsEnabled(0)
77   , fkEtaCorrection(NULL)
78   , fHasCutModel(kFALSE)
79   , fUseOnlyOROC(kFALSE)
80   , fNsigmaTPC(3)
81   , fRejectionEnabled(0)
82 {
83   //
84   // default  constructor
85   // 
86   //
87   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
88   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
89
90   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
91   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
92   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
93   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
94 }
95
96 //___________________________________________________________________
97 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
98   AliHFEpidBase("")
99   , fLineCrossingsEnabled(0)
100   , fkEtaCorrection(NULL)
101   , fHasCutModel(ref.fHasCutModel)
102   , fUseOnlyOROC(ref.fUseOnlyOROC)
103   , fNsigmaTPC(2)
104   , fRejectionEnabled(0)
105 {
106   //
107   // Copy constructor
108   //
109   ref.Copy(*this);
110 }
111
112 //___________________________________________________________________
113 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
114   //
115   // Assignment operator
116   //
117   if(this != &ref){
118     ref.Copy(*this);
119   } 
120   return *this;
121 }
122 //___________________________________________________________________
123 void AliHFEpidTPC::Copy(TObject &o) const{
124   //
125   // Copy function 
126   // called in copy constructor and assigment operator
127   //
128   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
129
130   target.fkEtaCorrection =fkEtaCorrection;
131   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
132   target.fHasCutModel = fHasCutModel;
133   target.fUseOnlyOROC = fUseOnlyOROC;
134   target.fNsigmaTPC = fNsigmaTPC;
135   target.fRejectionEnabled = fRejectionEnabled;
136
137   memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
138   memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
139
140   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
141   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
142   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
143  
144   AliHFEpidBase::Copy(target);
145 }
146
147 //___________________________________________________________________
148 AliHFEpidTPC::~AliHFEpidTPC(){
149   //
150   // Destructor
151   //
152 }
153
154 //___________________________________________________________________
155 Bool_t AliHFEpidTPC::InitializePID(Int_t /*run*/){
156   //
157   // Add TPC dE/dx Line crossings
158   //
159   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
160   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
161   return kTRUE;
162 }
163
164 //___________________________________________________________________
165 Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
166 {
167   //
168   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
169   // for electrons
170   // exclusion of the crossing points
171   //
172
173   if(!fkPIDResponse) return 0;
174
175   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
176
177   // Make clone track in order to be able to apply correction
178   const AliVTrack *rectrack;
179   AliESDtrack esdtrack;
180   AliAODTrack aodtrack;
181   if(fUseOnlyOROC && !fkEtaCorrection) {
182     if(track->IsESDanalysis()){
183       esdtrack.~AliESDtrack();
184       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
185       UseOROC(&esdtrack, anatype);
186       rectrack = &esdtrack;
187     } else {
188       aodtrack.~AliAODTrack();
189       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
190       UseOROC(&aodtrack, anatype);
191       rectrack = &aodtrack;
192     }
193   }
194   else if(fkEtaCorrection){
195     // Correction available
196     // apply it on copy
197     if(track->IsESDanalysis()){
198       esdtrack.~AliESDtrack();
199       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
200       ApplyEtaCorrection(&esdtrack, anatype);
201       rectrack = &esdtrack;
202     } else {
203       aodtrack.~AliAODTrack();
204       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
205       ApplyEtaCorrection(&aodtrack, anatype);
206       rectrack = &aodtrack;
207     }
208   } else {
209     // Correction not available - no need to copy
210     rectrack = track->GetRecTrack();
211   }
212   AliHFEpidObject tpctrack(*track);
213   tpctrack.SetRecTrack(rectrack);
214   
215   // QA before selection (after correction)
216   if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
217   AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
218   
219   // make copy of the track in order to allow for applying the correction 
220   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron);
221   AliDebug(1, Form("TPC NSigma: %f", nsigma));
222   // exclude crossing points:
223   // Determine the bethe values for each particle species
224   Bool_t isLineCrossing = kFALSE;
225   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
226     if(ispecies == AliPID::kElectron) continue;
227     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
228     if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
229       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
230       isLineCrossing = kTRUE;      
231       break;
232     }
233   }
234   if(isLineCrossing) return 0;
235
236   // Check particle rejection
237   if(HasParticleRejection()){
238     Int_t reject = Reject(rectrack, anatype);
239     if(reject != 0) return reject;
240   }
241
242   // Check if we have an asymmetric sigma model set
243   Int_t pdg = 0;
244   if(fHasCutModel){
245     pdg = CutSigmaModel(&tpctrack) ? 11 : 0;
246   } else { 
247     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
248     Float_t p = rectrack->P();
249     if(HasAsymmetricSigmaCut()) {
250       
251       //printf("p %f, fPAsigCut[0] %f, fPAsigCut[1] %f\n",p,fPAsigCut[0],fPAsigCut[1]);
252       if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) { 
253         if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; 
254         else pdg = 0;
255       }
256       else pdg = 0;
257     
258     }
259     else {
260       if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
261     }
262   }
263   if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
264   return pdg;
265
266 }
267
268 //___________________________________________________________________
269 Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
270   //
271   // N SigmaCut using parametrization of the cuts
272   //
273   Bool_t isSelected = kTRUE;
274   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
275   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron);
276   Double_t p = GetP(track->GetRecTrack(), anatype);
277   Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
278   AliDebug(2, Form("Centrality: %d\n", centrality));
279   if(centrality > 11) return kFALSE;
280   const TF1 *cutfunction;
281   if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
282   if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
283   return isSelected;
284 }
285
286 //___________________________________________________________________
287 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
288   //
289   // reject particles based on asymmetric sigma cut
290   //
291   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
292   Double_t p = GetP(track, anaType);
293   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
294     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
295     // Particle rejection enabled
296     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
297     Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
298     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
299   }
300   return 0;
301 }
302
303 //___________________________________________________________________
304 void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
305   //
306   // Apply correction for the eta dependence
307   // N.B. This correction has to be applied on a copy track
308   //
309   AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName()));
310   Double_t original = track->GetTPCsignal(),
311            eta = track->Eta();
312   if(anatype == AliHFEpidObject::kESDanalysis){
313     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
314     esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
315   } else {
316     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
317     AliAODPid *pid = aodtrack->GetDetPid();
318     if(pid) pid->SetTPCsignal(original);
319   }
320 }
321 //___________________________________________________________________
322 void AliHFEpidTPC::UseOROC(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
323   //
324   // Use TPC signal from the OROC
325   // N.B. This correction has to be applied on a copy track
326   //
327   //Double_t original = track->GetTPCsignal();
328   
329   if(anatype == AliHFEpidObject::kESDanalysis){
330     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
331     AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
332     Double32_t  TPCsignalRegion[4]; // TPC dEdx signal in 4 different regions - 0 - IROC, 1- OROC medium, 2 - OROC long, 3- OROC all, (default truncation used)
333     Char_t      TPCsignalNRegion[3]; // number of clusters above threshold used in the dEdx calculation
334     Char_t      TPCsignalNRowRegion[3]; // number of crosed rows used in the dEdx calculation - signal below threshold included
335     dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
336     esdtrack->SetTPCsignal(TPCsignalRegion[3],esdtrack->GetTPCsignalSigma(),(TPCsignalNRegion[1]+TPCsignalNRegion[2])); // the two last are not ok
337   } else {
338     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
339     AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
340     Double32_t  TPCsignalRegion[4]; // TPC dEdx signal in 4 different regions - 0 - IROC, 1- OROC medium, 2 - OROC long, 3- OROC all, (default truncation used)
341     Char_t      TPCsignalNRegion[3]; // number of clusters above threshold used in the dEdx calculation
342     Char_t      TPCsignalNRowRegion[3]; // number of crosed rows used in the dEdx calculation - signal below threshold included
343     dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
344     AliAODPid *pid = aodtrack->GetDetPid();
345     if(pid) pid->SetTPCsignal(TPCsignalRegion[3]);
346     if(pid) pid->SetTPCsignalN((TPCsignalNRegion[1]+TPCsignalNRegion[2]));
347   }
348
349 }
350 //___________________________________________________________________
351 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
352   //
353   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
354   // Stores line center and line sigma
355   //
356   if(species >= AliPID::kSPECIES){
357     AliError("Species doesn't exist");
358     return;
359   }
360   fLineCrossingsEnabled |= 1 << species;
361   fLineCrossingSigma[species] = sigma;
362 }
363
364 //___________________________________________________________________
365 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
366   //
367   // Get the momentum at the inner wall of the TPC
368   //
369   Double_t p = -1;
370   if(anatype == AliHFEpidObject::kESDanalysis){
371     // ESD analysis: Use Inner Params for the momentum estimate
372     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
373     if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
374   } else { 
375     // AOD analysis: Use TPC momentum stored in the AliAODpid object
376     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
377     if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
378   }
379   return p;
380 }