Coding rule violations corrected
[u/mrichter/AliRoot.git] / STEER / AliAODpidUtil.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 /* $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
25 #include "AliLog.h"
26 #include "AliPID.h"
27 #include "AliAODpidUtil.h"
28 #include "AliAODEvent.h"
29 #include "AliAODTrack.h"
30 #include "AliAODPid.h"
31 #include "AliTRDPIDResponse.h"
32 #include "AliESDtrack.h"
33
34 ClassImp(AliAODpidUtil)
35
36   Int_t AliAODpidUtil::MakePID(AliAODTrack *track,Double_t *p) const {
37   //
38   //  Calculate probabilities for all detectors, except if TPConly==kTRUE
39   //  and combine PID
40   //  
41   //   Option TPConly==kTRUE is used during reconstruction, 
42   //  because ITS tracking uses TPC pid
43   //  HMPID and TRD pid are done in detector reconstructors
44   //
45
46   /*
47     Float_t TimeZeroTOF = 0;
48     if (subtractT0) 
49     TimeZeroTOF = event->GetT0();
50   */
51   Int_t ns=AliPID::kSPECIES;
52   Double_t tpcPid[AliPID::kSPECIES];
53   MakeTPCPID(track,tpcPid);
54   Double_t itsPid[AliPID::kSPECIES];
55   Double_t tofPid[AliPID::kSPECIES];
56   Double_t trdPid[AliPID::kSPECIES];
57   MakeITSPID(track,itsPid);
58   MakeTOFPID(track,tofPid);
59   //MakeHMPIDPID(track);
60   MakeTRDPID(track,trdPid);
61   for (Int_t j=0; j<ns; j++) {
62     p[j]=tpcPid[j]*itsPid[j]*tofPid[j]*trdPid[j];
63   }
64
65   return 0;
66 }
67 //_________________________________________________________________________
68 void AliAODpidUtil::MakeTPCPID(const AliAODTrack *track,Double_t *p) const
69 {
70   //
71   //  TPC pid using bethe-bloch and gaussian response
72   //
73
74   if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return;
75
76   AliAODPid *pidObj = track->GetDetPid();
77   Double_t mom      = track->P();
78   Double_t dedx     = 0.;  
79   UShort_t nTPCClus = 0;
80   if (pidObj) {
81       nTPCClus = pidObj->GetTPCsignalN();
82       dedx     = pidObj->GetTPCsignal();
83       mom      = pidObj->GetTPCmomentum();
84   }
85   
86   Bool_t mismatch=kTRUE;
87
88   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
89     AliPID::EParticleType type=AliPID::EParticleType(j);
90     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type); 
91     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,nTPCClus,type);
92     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
93       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
94     } else {
95       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
96       mismatch=kFALSE;
97     }
98
99   }
100
101   if (mismatch)
102     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
103
104
105   return;
106 }
107 //_________________________________________________________________________
108 void AliAODpidUtil::MakeITSPID(const AliAODTrack *track,Double_t *p) const
109 {
110   //
111   // ITS PID
112   //  1) Truncated mean method
113   //
114
115
116   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
117   UChar_t clumap=track->GetITSClusterMap();
118   Int_t nPointsForPid=0;
119   for(Int_t i=2; i<6; i++){
120    if(clumap&(1<<i)) ++nPointsForPid;
121   }
122   if(nPointsForPid<3) { // track not to be used for combined PID purposes
123     for (Int_t j=0; j<AliPID::kSPECIES; j++) 
124       p[j] = 1./AliPID::kSPECIES;
125     return;
126   }
127   Double_t mom=track->P();  
128   AliAODPid *pidObj = track->GetDetPid();
129
130   Double_t dedx = 0.;
131   if (pidObj) {
132       dedx = pidObj->GetITSsignal();
133   }
134   
135   Bool_t mismatch = kTRUE;
136   Bool_t isSA = kTRUE;
137   if(track->GetStatus() & AliESDtrack::kTPCin){
138     isSA = kFALSE;
139     if (pidObj)
140       mom = pidObj->GetTPCmomentum();
141   }
142   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
143     Double_t mass = AliPID::ParticleMass(j);//GeV/c^2
144     Double_t bethe = fITSResponse.Bethe(mom,mass);
145     Double_t sigma = fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
146     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
147       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
148     } else {
149       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
150       mismatch=kFALSE;
151     }
152
153     // Check for particles heavier than (AliPID::kSPECIES - 1)
154
155   }
156
157   if (mismatch)
158     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
159
160   return;
161
162 }
163 //_________________________________________________________________________
164 void AliAODpidUtil::MakeTOFPID(const AliAODTrack *track, Double_t *p) const
165 {
166   //
167   //   TOF PID using gaussian response
168   //
169   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return;
170   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return;
171   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return;
172
173   Double_t time[AliPID::kSPECIESN];
174   Double_t sigma[AliPID::kSPECIESN];
175   AliAODPid *pidObj = track->GetDetPid();
176   pidObj->GetIntegratedTimes(time);
177
178   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
179     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
180   }
181
182   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
183                   Form("Expected TOF signals [ps]: %f %f %f %f %f",
184                        time[AliPID::kElectron],
185                        time[AliPID::kMuon],
186                        time[AliPID::kPion],
187                        time[AliPID::kKaon],
188                        time[AliPID::kProton]));
189
190   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
191                   Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
192                        sigma[AliPID::kElectron],
193                        sigma[AliPID::kMuon],
194                        sigma[AliPID::kPion],
195                        sigma[AliPID::kKaon],
196                        sigma[AliPID::kProton]
197                        ));
198
199   Double_t tof = pidObj->GetTOFsignal();
200
201   Bool_t mismatch = kTRUE;
202   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
203     Double_t sig = sigma[j];
204     if (TMath::Abs(tof-time[j]) > fRange*sig) {
205       p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
206     } else
207       p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
208
209     // Check the mismatching
210     Double_t mass = AliPID::ParticleMass(j);
211     Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
212     if (p[j]>pm) mismatch = kFALSE;
213
214     // Check for particles heavier than (AliPID::kSPECIES - 1)
215
216   }
217
218   if (mismatch)
219     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
220
221   return;
222 }
223 //_________________________________________________________________________
224 void AliAODpidUtil::MakeTRDPID(const AliAODTrack *track,Double_t *p) const
225 {
226   
227   // Method to recalculate the TRD PID probabilities
228   if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return;
229
230   AliAODPid *pidObj = track->GetDetPid();
231   Float_t *mom=pidObj->GetTRDmomentum();
232   Int_t ntracklets=0;
233   for(Int_t iPl=0;iPl<6;iPl++){
234    if(mom[iPl]>0.) ntracklets++;
235   }
236    if(ntracklets<4) return;
237
238   Double_t* dedx=pidObj->GetTRDsignal();
239   Bool_t norm=kTRUE;
240   fTRDResponse.GetResponse(pidObj->GetTRDnSlices(),dedx,mom,p,norm);
241   return;
242 }
243
244 //_________________________________________________________________________
245 Float_t AliAODpidUtil::NumberOfSigmasTPC(const AliAODTrack *track, AliPID::EParticleType type) const {
246   
247   Double_t mom = 0.0;
248   AliAODPid *pidObj = 0x0;
249   if (track) {
250     mom = track->P();
251     pidObj = track->GetDetPid();
252   }
253   UShort_t nTPCClus=0;
254   Double_t tpcSignal=0.0;
255   if (pidObj) {
256     nTPCClus=pidObj->GetTPCsignalN();
257     mom = pidObj->GetTPCmomentum();
258     tpcSignal = pidObj->GetTPCsignal();
259   }
260   return fTPCResponse.GetNumberOfSigmas(mom,tpcSignal,nTPCClus,type); 
261 }
262
263 //_________________________________________________________________________
264 Float_t AliAODpidUtil::NumberOfSigmasTOF(const AliAODTrack *track, AliPID::EParticleType type) const {
265   Double_t times[AliPID::kSPECIES]={0.};
266   Double_t sigmaTOFPid[AliPID::kSPECIES]={0.};
267   AliAODPid *pidObj = 0x0;
268   Double_t mom = 0.0;
269   if (track) {
270     pidObj = track->GetDetPid();
271     mom = track->P();
272   }
273   Double_t tofSignal = 0.0;
274   if (pidObj) {
275     pidObj->GetIntegratedTimes(times);
276     pidObj->GetTOFpidResolution(sigmaTOFPid);
277     tofSignal = pidObj->GetTOFsignal();
278   }
279   if (sigmaTOFPid[type]>0) return (tofSignal - times[type])/sigmaTOFPid[type];
280   else return (tofSignal - times[type])/fTOFResponse.GetExpectedSigma(mom,times[type],AliPID::ParticleMass(type));
281 }
282
283 //_________________________________________________________________________
284 Float_t AliAODpidUtil::NumberOfSigmasITS(const AliAODTrack *track, AliPID::EParticleType type) const {
285   AliAODPid *pidObj = 0x0;
286   UChar_t clumap=0;
287   Float_t mom=0.0;
288   UShort_t nTPCClus=0;
289   if (track) {
290     pidObj = track->GetDetPid();
291     clumap=track->GetITSClusterMap();
292     mom=track->P();
293     nTPCClus=track->GetTPCNcls();
294   }
295   Int_t nPointsForPid=0;
296   for(Int_t i=2; i<6; i++){
297     if(clumap&(1<<i)) ++nPointsForPid;
298   }
299   Float_t dEdx=0.0;
300   if (pidObj) dEdx=pidObj->GetITSsignal();
301   Bool_t isSA=kTRUE;  
302   if(nTPCClus>0) isSA=kFALSE;
303   return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA);
304 }