Line data Source code
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: AliESDpid.cxx 64123 2013-09-05 15:09:53Z morsch $ */
17 :
18 : //-----------------------------------------------------------------
19 : // Implementation of the combined PID class
20 : // For the Event Summary Data Class
21 : // produced by the reconstruction process
22 : // and containing information on the particle identification
23 : // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
24 : //-----------------------------------------------------------------
25 :
26 : #include "TArrayI.h"
27 : #include "TArrayF.h"
28 :
29 : #include "AliLog.h"
30 : #include "AliPID.h"
31 : #include "AliTOFHeader.h"
32 : #include "AliESDpid.h"
33 : #include "AliESDEvent.h"
34 : #include "AliESDtrack.h"
35 : #include "AliMCEvent.h"
36 : #include "AliTOFPIDParams.h"
37 :
38 : #include <AliDetectorPID.h>
39 :
40 172 : ClassImp(AliESDpid)
41 :
42 : Int_t AliESDpid::MakePID(AliESDEvent *event, Bool_t TPConly, Float_t /*timeZeroTOF*/) const {
43 : //
44 : // Calculate probabilities for all detectors, except if TPConly==kTRUE
45 : // and combine PID
46 : //
47 : // Option TPConly==kTRUE is used during reconstruction,
48 : // because ITS tracking uses TPC pid
49 : // HMPID and TRD pid are done in detector reconstructors
50 : //
51 :
52 : /*
53 : Float_t timeZeroTOF = 0;
54 : if (subtractT0)
55 : timeZeroTOF = event->GetT0();
56 : */
57 0 : Int_t nTrk=event->GetNumberOfTracks();
58 0 : for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
59 0 : AliESDtrack *track=event->GetTrack(iTrk);
60 0 : MakeTPCPID(track);
61 0 : if (!TPConly) {
62 0 : MakeITSPID(track);
63 : //MakeTOFPID(track, timeZeroTOF);
64 : //MakeHMPIDPID(track);
65 : //MakeTRDPID(track);
66 0 : }
67 0 : CombinePID(track);
68 : }
69 0 : return 0;
70 : }
71 :
72 : //_________________________________________________________________________
73 : void AliESDpid::MakeTPCPID(AliESDtrack *track) const
74 : {
75 : //
76 : // TPC pid using bethe-bloch and gaussian response
77 : //
78 0 : if ((track->GetStatus()&AliESDtrack::kTPCin )==0)
79 0 : if ((track->GetStatus()&AliESDtrack::kTPCout)==0) return;
80 :
81 0 : Double_t mom = track->GetP();
82 0 : const AliExternalTrackParam *in=track->GetInnerParam();
83 0 : if (in) mom = in->GetP();
84 :
85 0 : Double_t p[AliPID::kSPECIES];
86 0 : Double_t dedx=track->GetTPCsignal();
87 : Bool_t mismatch=kTRUE, heavy=kTRUE;
88 :
89 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) {
90 : AliPID::EParticleType type=AliPID::EParticleType(j);
91 0 : Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
92 0 : Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
93 0 : if (TMath::Abs(dedx-bethe) > fRange*sigma) {
94 0 : p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
95 0 : } else {
96 0 : p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
97 : mismatch=kFALSE;
98 : }
99 :
100 : // Check for particles heavier than (AliPID::kSPECIES - 1)
101 0 : if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
102 :
103 : }
104 :
105 0 : if (mismatch)
106 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
107 :
108 0 : track->SetTPCpid(p);
109 :
110 0 : if (heavy) track->ResetStatus(AliESDtrack::kTPCpid);
111 :
112 0 : }
113 : //_________________________________________________________________________
114 : void AliESDpid::MakeITSPID(AliESDtrack *track) const
115 : {
116 : //
117 : // ITS PID
118 : // Two options, depending on fITSPIDmethod:
119 : // 1) Truncated mean method
120 : // 2) Likelihood, using charges measured in all 4 layers and
121 : // Landau+gaus response functions
122 : //
123 :
124 0 : if ((track->GetStatus()&AliESDtrack::kITSin)==0 &&
125 0 : (track->GetStatus()&AliESDtrack::kITSout)==0) return;
126 :
127 0 : Double_t mom=track->GetP();
128 0 : if (fITSPIDmethod == kITSTruncMean) {
129 0 : Double_t dedx=track->GetITSsignal();
130 : Bool_t isSA=kTRUE;
131 : Double_t momITS=mom;
132 0 : ULong_t trStatus=track->GetStatus();
133 0 : if(trStatus&AliESDtrack::kTPCin) isSA=kFALSE;
134 0 : UChar_t clumap=track->GetITSClusterMap();
135 : Int_t nPointsForPid=0;
136 0 : for(Int_t i=2; i<6; i++){
137 0 : if(clumap&(1<<i)) ++nPointsForPid;
138 : }
139 :
140 0 : if(nPointsForPid<3) { // track not to be used for combined PID purposes
141 0 : track->ResetStatus(AliESDtrack::kITSpid);
142 0 : return;
143 : }
144 :
145 0 : Double_t p[AliPID::kSPECIES];
146 :
147 : Bool_t mismatch=kTRUE, heavy=kTRUE;
148 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) {
149 0 : Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
150 0 : Double_t bethe=fITSResponse.Bethe(momITS,mass);
151 0 : Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
152 0 : if (TMath::Abs(dedx-bethe) > fRange*sigma) {
153 0 : p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
154 0 : } else {
155 0 : p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
156 : mismatch=kFALSE;
157 : }
158 :
159 : // Check for particles heavier than (AliPID::kSPECIES - 1)
160 0 : if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
161 :
162 : }
163 :
164 0 : if (mismatch)
165 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
166 :
167 0 : track->SetITSpid(p);
168 :
169 0 : if (heavy) track->ResetStatus(AliESDtrack::kITSpid);
170 0 : }
171 : else { // Likelihood method
172 0 : Double_t condprobfun[AliPID::kSPECIES];
173 0 : Double_t qclu[4];
174 0 : track->GetITSdEdxSamples(qclu);
175 0 : fITSResponse.GetITSProbabilities(mom,qclu,condprobfun);
176 0 : track->SetITSpid(condprobfun);
177 0 : }
178 :
179 0 : }
180 : //_________________________________________________________________________
181 : void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t /*timeZeroTOF*/) const
182 : {
183 : //
184 : // TOF PID using gaussian response
185 : //
186 :
187 0 : if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return;
188 0 : if ((track->GetStatus()&AliESDtrack::kTIME)==0) return;
189 0 : if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
190 :
191 0 : Int_t ibin = fTOFResponse.GetMomBin(track->GetP());
192 0 : Float_t timezero = fTOFResponse.GetT0bin(ibin);
193 :
194 0 : Double_t time[AliPID::kSPECIESC];
195 0 : track->GetIntegratedTimes(time);
196 :
197 0 : Double_t sigma[AliPID::kSPECIES];
198 0 : for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
199 0 : sigma[iPart] = fTOFResponse.GetExpectedSigma(track->GetP(),time[iPart],AliPID::ParticleMass(iPart));
200 : }
201 :
202 0 : AliDebugGeneral("AliESDpid::MakeTOFPID",2,
203 : Form("Expected TOF signals [ps]: %f %f %f %f %f",
204 : time[AliPID::kElectron],
205 : time[AliPID::kMuon],
206 : time[AliPID::kPion],
207 : time[AliPID::kKaon],
208 : time[AliPID::kProton]));
209 :
210 0 : AliDebugGeneral("AliESDpid::MakeTOFPID",2,
211 : Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
212 : sigma[AliPID::kElectron],
213 : sigma[AliPID::kMuon],
214 : sigma[AliPID::kPion],
215 : sigma[AliPID::kKaon],
216 : sigma[AliPID::kProton]
217 : ));
218 :
219 0 : Double_t tof = track->GetTOFsignal() - timezero;
220 :
221 0 : Double_t p[AliPID::kSPECIES];
222 : // Bool_t mismatch = kTRUE;
223 : Bool_t heavy = kTRUE;
224 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) {
225 0 : Double_t sig = sigma[j];
226 0 : if (TMath::Abs(tof-time[j]) > (fRange+2)*sig) {
227 0 : p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
228 0 : } else
229 0 : p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
230 :
231 : // Check the mismatching
232 :
233 : // Double_t mass = AliPID::ParticleMass(j);
234 : // Double_t pm = fTOFResponse.GetMismatchProbability(track->GetP(),mass);
235 : // if (p[j]>pm) mismatch = kFALSE;
236 :
237 : // Check for particles heavier than (AliPID::kSPECIES - 1)
238 0 : if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
239 :
240 : }
241 :
242 : /*
243 : if (mismatch)
244 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
245 : */
246 :
247 0 : track->SetTOFpid(p);
248 :
249 0 : if (heavy) track->ResetStatus(AliESDtrack::kTOFpid);
250 :
251 : // kTOFmismatch flas is not set because deprecated from 18/02/2013
252 : // if (!CheckTOFMatching(track)) track->SetStatus(AliESDtrack::kTOFmismatch);
253 : // else track->ResetStatus(AliESDtrack::kTOFmismatch);
254 0 : }
255 : //_________________________________________________________________________
256 : void AliESDpid::MakeTRDPID(AliESDtrack *track) const
257 : {
258 : //
259 : // Method to recalculate the TRD PID probabilities
260 : //
261 0 : Double_t prob[AliPID::kSPECIES];
262 0 : GetComputeTRDProbability(track, AliPID::kSPECIES, prob);
263 0 : track->SetTRDpid(prob);
264 0 : }
265 : //_________________________________________________________________________
266 : void AliESDpid::CombinePID(AliESDtrack *track) const
267 : {
268 : //
269 : // Combine the information of various detectors
270 : // to determine the Particle Identification
271 : //
272 0 : Double_t p[AliPID::kSPECIES]={1.};
273 :
274 0 : if (track->IsOn(AliESDtrack::kITSpid)) {
275 0 : Double_t d[AliPID::kSPECIES];
276 0 : track->GetITSpid(d);
277 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
278 0 : }
279 :
280 0 : if (track->IsOn(AliESDtrack::kTPCpid)) {
281 0 : Double_t d[AliPID::kSPECIES];
282 0 : track->GetTPCpid(d);
283 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
284 0 : }
285 :
286 0 : if (track->IsOn(AliESDtrack::kTRDpid)) {
287 0 : Double_t d[AliPID::kSPECIES];
288 0 : track->GetTRDpid(d);
289 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
290 0 : }
291 :
292 0 : if (track->IsOn(AliESDtrack::kTOFpid)) {
293 0 : Double_t d[AliPID::kSPECIES];
294 0 : track->GetTOFpid(d);
295 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
296 0 : }
297 :
298 0 : if (track->IsOn(AliESDtrack::kHMPIDpid)) {
299 0 : Double_t d[AliPID::kSPECIES];
300 0 : track->GetHMPIDpid(d);
301 0 : for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]*=d[j];
302 0 : }
303 :
304 0 : track->SetESDpid(p);
305 0 : }
306 : //_________________________________________________________________________
307 : Bool_t AliESDpid::CheckTOFMatching(AliESDtrack *track) const{
308 : //
309 : // Check pid matching of TOF with TPC as reference
310 : //
311 : Bool_t status = kFALSE;
312 :
313 0 : Double_t exptimes[AliPID::kSPECIESC];
314 0 : track->GetIntegratedTimes(exptimes);
315 :
316 0 : Float_t p = track->P();
317 :
318 0 : Float_t dedx = track->GetTPCsignal();
319 0 : Float_t time = track->GetTOFsignal() - fTOFResponse.GetStartTime(p);
320 :
321 0 : Double_t ptpc[3];
322 0 : track->GetInnerPxPyPz(ptpc);
323 0 : Float_t momtpc=TMath::Sqrt(ptpc[0]*ptpc[0] + ptpc[1]*ptpc[1] + ptpc[2]*ptpc[2]);
324 :
325 0 : for(Int_t i=0;i < AliPID::kSPECIES;i++){
326 : AliPID::EParticleType type=AliPID::EParticleType(i);
327 :
328 0 : Float_t resolutionTOF = fTOFResponse.GetExpectedSigma(p, exptimes[i], AliPID::ParticleMass(i));
329 0 : if(TMath::Abs(exptimes[i] - time) < fRange * resolutionTOF){
330 0 : Float_t dedxExp = fTPCResponse.GetExpectedSignal(momtpc,type);
331 0 : Float_t resolutionTPC = fTPCResponse.GetExpectedSigma(momtpc,track->GetTPCsignalN(),type);
332 :
333 0 : if(TMath::Abs(dedx - dedxExp) < fRangeTOFMismatch * resolutionTPC){
334 : status = kTRUE;
335 0 : }
336 0 : }
337 : }
338 :
339 : // for nuclei
340 0 : Float_t resolutionTOFpr = fTOFResponse.GetExpectedSigma(p, exptimes[4], AliPID::ParticleMass(4));
341 0 : if(!status && (exptimes[4] + fRange*resolutionTOFpr < time)) status = kTRUE;
342 :
343 :
344 0 : return status;
345 0 : }
346 :
347 : //_________________________________________________________________________
348 : Float_t AliESDpid::GetSignalDeltaTOFold(const AliVParticle *track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
349 : {
350 : //
351 : // TOF signal - expected
352 : //
353 0 : AliVTrack *vtrack=(AliVTrack*)track;
354 :
355 0 : const Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
356 : Double_t tofTime = 99999;
357 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) tofTime = (Double_t)this->GetTOFsignalTunedOnData((AliVTrack*)vtrack);
358 0 : else tofTime=vtrack->GetTOFsignal();
359 0 : tofTime = tofTime - fTOFResponse.GetStartTime(vtrack->P());
360 : Double_t delta=-9999.;
361 :
362 0 : if (!ratio) delta=tofTime-expTime;
363 0 : else if (expTime>1.e-20) delta=tofTime/expTime;
364 :
365 0 : return delta;
366 : }
367 :
368 : //_________________________________________________________________________
369 : Float_t AliESDpid::GetNumberOfSigmasTOFold(const AliVParticle *track, AliPID::EParticleType type) const
370 : {
371 : //
372 : // Number of sigma implementation for the TOF
373 : //
374 :
375 0 : AliVTrack *vtrack=(AliVTrack*)track;
376 : Double_t tofTime = 99999;
377 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) tofTime = (Double_t)this->GetTOFsignalTunedOnData((AliVTrack*)vtrack);
378 0 : else tofTime=vtrack->GetTOFsignal();
379 :
380 0 : Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
381 0 : return (tofTime - fTOFResponse.GetStartTime(vtrack->P()) - expTime)/fTOFResponse.GetExpectedSigma(vtrack->P(),expTime,AliPID::ParticleMassZ(type));
382 : }
383 :
384 : //_________________________________________________________________________
385 : void AliESDpid::SetPIDForTracking(AliESDtrack *esdtr) const
386 : {
387 : // assign mass for tracking
388 : //
389 :
390 : // in principle AliPIDCombined could be used to also set priors
391 : //AliPIDCombined pidProb;
392 : //pidProb.SetDetectorMask(kDetTPC); // use only TPC information, couls also be changed
393 : //pidProb.SetSelectedSpecies(AliPID::kSPECIESC); // number of species to use
394 : //pidProb.SetDefaultTPCPriors(); // load default priors
395 :
396 576 : Double_t prob[AliPID::kSPECIESC]={0.};
397 288 : EDetPidStatus pidStatus=ComputePIDProbability(kTPC, esdtr, AliPID::kSPECIESC, prob);
398 : // check if a valid signal was found, otherwise return pion mass
399 288 : if (pidStatus==kDetNoSignal) { //kDetPidOk) {
400 154 : esdtr->SetPIDForTracking(AliPID::kPion);
401 154 : return;
402 : }
403 :
404 : // or with AliPIDCombined
405 : // pidProb.ComputeProbabilities(esdtr, this, p);
406 :
407 : // find max probability
408 : Float_t max=0.,min=1.e9;
409 : Int_t pid=-1;
410 2680 : for (Int_t i=0; i<AliPID::kSPECIESC; ++i) {
411 1512 : if (prob[i]>max) {pid=i; max=prob[i];}
412 1998 : if (prob[i]<min) min=prob[i];
413 : }
414 :
415 : //int pid = AliPID::kPion; // this should be substituted by real most probable TPC pid (e,mu -> pion) or poin if no PID possible
416 :
417 : //
418 268 : if (pid>AliPID::kSPECIESC-1 || (min>=max)) pid = AliPID::kPion;
419 : //
420 134 : esdtr->SetPIDForTracking( pid );
421 : //
422 422 : }
423 :
424 :
425 : //_________________________________________________________________________
426 : void AliESDpid::MakePIDForTracking(AliESDEvent *event) const
427 : {
428 : // assign masses using for tracking
429 32 : Int_t nTrk=event->GetNumberOfTracks();
430 608 : for (Int_t iTrk=nTrk; iTrk--;) {
431 288 : AliESDtrack *track = event->GetTrack(iTrk);
432 288 : SetPIDForTracking(track);
433 : }
434 :
435 16 : }
|
