In AddTaskPHOSPi0Flow.C set Cent. Bin past event buffers/lists,
[u/mrichter/AliRoot.git] / STEER / AOD / AliAODpidUtil.cxx
CommitLineData
1423bac9 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/* $Id: AliAODpidUtil.cxx 38329 2010-01-17 19:17:24Z hristov $ */
17
18//-----------------------------------------------------------------
19// Implementation of the combined PID class
20// For the AOD Class
21// containing information on the particle identification
22// Origin: Rosa Romita, GSI, r.romita@gsi.de
23//-----------------------------------------------------------------
24
539a5a59 25#include "TRandom.h"
1423bac9 26#include "AliLog.h"
27#include "AliPID.h"
28#include "AliAODpidUtil.h"
29#include "AliAODEvent.h"
30#include "AliAODTrack.h"
31#include "AliAODPid.h"
32#include "AliTRDPIDResponse.h"
1c2fc6e0 33#include "AliESDtrack.h"
539a5a59 34#include "AliAODMCHeader.h"
35#include "AliAODMCParticle.h"
1423bac9 36
37ClassImp(AliAODpidUtil)
38
03945936 39 Int_t AliAODpidUtil::MakePID(const AliAODTrack *track,Double_t *p) const {
1423bac9 40 //
41 // Calculate probabilities for all detectors, except if TPConly==kTRUE
42 // and combine PID
43 //
44 // Option TPConly==kTRUE is used during reconstruction,
45 // because ITS tracking uses TPC pid
46 // HMPID and TRD pid are done in detector reconstructors
47 //
48
49 /*
50 Float_t TimeZeroTOF = 0;
51 if (subtractT0)
52 TimeZeroTOF = event->GetT0();
53 */
54 Int_t ns=AliPID::kSPECIES;
55 Double_t tpcPid[AliPID::kSPECIES];
56 MakeTPCPID(track,tpcPid);
57 Double_t itsPid[AliPID::kSPECIES];
58 Double_t tofPid[AliPID::kSPECIES];
59 Double_t trdPid[AliPID::kSPECIES];
60 MakeITSPID(track,itsPid);
1c2fc6e0 61 MakeTOFPID(track,tofPid);
1423bac9 62 //MakeHMPIDPID(track);
63 MakeTRDPID(track,trdPid);
64 for (Int_t j=0; j<ns; j++) {
65 p[j]=tpcPid[j]*itsPid[j]*tofPid[j]*trdPid[j];
66 }
67
68 return 0;
69}
70//_________________________________________________________________________
539a5a59 71Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
72 AliAODTrack *track = (AliAODTrack *) t;
73 Float_t dedx = track->GetTPCsignalTunedOnData();
74 if(dedx > 0) return dedx;
75
76 Double_t mom = t->GetTPCmomentum();
77
78 dedx = t->GetTPCsignal();
79 track->SetTPCsignalTunedOnData(dedx);
80
81 if(dedx < 20) return dedx;
82
83
84 AliPID::EParticleType type = AliPID::kPion;
85
86 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(track->GetAODEvent()->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
87 if (mcHeader) {
88 TClonesArray *mcArray = (TClonesArray*)track->GetAODEvent()->GetList()->FindObject(AliAODMCParticle::StdBranchName());
89
90 Bool_t kGood = kTRUE;
91
92 Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(t->GetLabel())))->GetPdgCode());
93 if(iS==AliPID::ParticleCode(AliPID::kElectron)){
94 type = AliPID::kElectron;
95 }
96 else if(iS==AliPID::ParticleCode(AliPID::kMuon)){
97 type = AliPID::kMuon;
98 }
99 else if(iS==AliPID::ParticleCode(AliPID::kPion)){
100 type = AliPID::kPion;
101 }
102 else if(iS==AliPID::ParticleCode(AliPID::kKaon)){
103 type = AliPID::kKaon;
104 }
105 else if(iS==AliPID::ParticleCode(AliPID::kProton)){
106 type = AliPID::kProton;
107 }
108 else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d
109 type = AliPID::kDeuteron;
110 }
111 else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t
112 type = AliPID::kTriton;
113 }
114 else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He
115 type = AliPID::kHe3;
116 }
117 else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He
118 type = AliPID::kAlpha;
119 }
120 else
121 kGood = kFALSE;
122
123 if(kGood){
124 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
125 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
126 dedx = gRandom->Gaus(bethe,sigma);
127
128 if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
129 }
130
131 }
132
133 track->SetTPCsignalTunedOnData(dedx);
134 return dedx;
135}
136//_________________________________________________________________________
5e6a3170 137void AliAODpidUtil::MakeTPCPID(const AliAODTrack *track,Double_t *p) const
1423bac9 138{
139 //
140 // TPC pid using bethe-bloch and gaussian response
141 //
142
1c2fc6e0 143 if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return;
1c2fc6e0 144
1423bac9 145 AliAODPid *pidObj = track->GetDetPid();
a67eebca 146 Double_t mom = track->P();
147 Double_t dedx = 0.;
148 UShort_t nTPCClus = 0;
149 if (pidObj) {
150 nTPCClus = pidObj->GetTPCsignalN();
151 dedx = pidObj->GetTPCsignal();
152 mom = pidObj->GetTPCmomentum();
153 }
154
1423bac9 155 Bool_t mismatch=kTRUE;
156
157 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
158 AliPID::EParticleType type=AliPID::EParticleType(j);
159 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
1c2fc6e0 160 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,nTPCClus,type);
1423bac9 161 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
162 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
163 } else {
164 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
165 mismatch=kFALSE;
166 }
167
168 }
169
170 if (mismatch)
77638021 171 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
1423bac9 172
173
174 return;
175}
176//_________________________________________________________________________
5e6a3170 177void AliAODpidUtil::MakeITSPID(const AliAODTrack *track,Double_t *p) const
1423bac9 178{
179 //
180 // ITS PID
181 // 1) Truncated mean method
182 //
183
184
1c2fc6e0 185 if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
186 UChar_t clumap=track->GetITSClusterMap();
187 Int_t nPointsForPid=0;
188 for(Int_t i=2; i<6; i++){
189 if(clumap&(1<<i)) ++nPointsForPid;
190 }
15e979c9 191 if(nPointsForPid<3) { // track not to be used for combined PID purposes
192 for (Int_t j=0; j<AliPID::kSPECIES; j++)
193 p[j] = 1./AliPID::kSPECIES;
194 return;
195 }
1423bac9 196 Double_t mom=track->P();
197 AliAODPid *pidObj = track->GetDetPid();
198
0f9da745 199 Double_t dedx = 0.;
200 if (pidObj) {
201 dedx = pidObj->GetITSsignal();
202 }
203
15e979c9 204 Bool_t mismatch = kTRUE;
205 Bool_t isSA = kTRUE;
206 if(track->GetStatus() & AliESDtrack::kTPCin){
207 isSA = kFALSE;
208 if (pidObj)
209 mom = pidObj->GetTPCmomentum();
210 }
1423bac9 211 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
15e979c9 212 Double_t mass = AliPID::ParticleMass(j);//GeV/c^2
213 Double_t bethe = fITSResponse.Bethe(mom,mass);
214 Double_t sigma = fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
1423bac9 215 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
216 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
217 } else {
218 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
219 mismatch=kFALSE;
220 }
221
222 // Check for particles heavier than (AliPID::kSPECIES - 1)
223
224 }
225
226 if (mismatch)
227 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
228
229 return;
230
231}
232//_________________________________________________________________________
5e6a3170 233void AliAODpidUtil::MakeTOFPID(const AliAODTrack *track, Double_t *p) const
1423bac9 234{
235 //
236 // TOF PID using gaussian response
237 //
1c2fc6e0 238 if ((track->GetStatus()&AliESDtrack::kTOFout )==0) return;
239 if ((track->GetStatus()&AliESDtrack::kTIME )==0) return;
240 if ((track->GetStatus()&AliESDtrack::kTOFpid )==0) return;
1423bac9 241
242 Double_t time[AliPID::kSPECIESN];
243 Double_t sigma[AliPID::kSPECIESN];
244 AliAODPid *pidObj = track->GetDetPid();
245 pidObj->GetIntegratedTimes(time);
246
247 for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
248 sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
249 }
250
251 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
252 Form("Expected TOF signals [ps]: %f %f %f %f %f",
253 time[AliPID::kElectron],
254 time[AliPID::kMuon],
255 time[AliPID::kPion],
256 time[AliPID::kKaon],
257 time[AliPID::kProton]));
258
259 AliDebugGeneral("AliESDpid::MakeTOFPID",2,
260 Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
261 sigma[AliPID::kElectron],
262 sigma[AliPID::kMuon],
263 sigma[AliPID::kPion],
264 sigma[AliPID::kKaon],
265 sigma[AliPID::kProton]
266 ));
267
1c2fc6e0 268 Double_t tof = pidObj->GetTOFsignal();
1423bac9 269
270 Bool_t mismatch = kTRUE;
271 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
272 Double_t sig = sigma[j];
273 if (TMath::Abs(tof-time[j]) > fRange*sig) {
274 p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
275 } else
276 p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
277
278 // Check the mismatching
279 Double_t mass = AliPID::ParticleMass(j);
280 Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
281 if (p[j]>pm) mismatch = kFALSE;
282
283 // Check for particles heavier than (AliPID::kSPECIES - 1)
284
285 }
286
287 if (mismatch)
77638021 288 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
1423bac9 289
290 return;
291}
292//_________________________________________________________________________
5e6a3170 293void AliAODpidUtil::MakeTRDPID(const AliAODTrack *track,Double_t *p) const
1423bac9 294{
ea235c90 295 ComputeTRDProbability(track, AliPID::kSPECIES, p);
1423bac9 296 return;
297}
ea235c90 298