Bug 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(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   Double_t mom = track->P();
77   AliAODPid *pidObj = track->GetDetPid();
78   if (pidObj) mom = pidObj->GetTPCmomentum();
79   UShort_t nTPCClus=pidObj->GetTPCsignalN();
80    
81   Double_t dedx=pidObj->GetTPCsignal(); 
82   Bool_t mismatch=kTRUE;
83
84   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
85     AliPID::EParticleType type=AliPID::EParticleType(j);
86     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type); 
87     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,nTPCClus,type);
88     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
89       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
90     } else {
91       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
92       mismatch=kFALSE;
93     }
94
95   }
96
97   if (mismatch)
98     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
99
100
101   return;
102 }
103 //_________________________________________________________________________
104 void AliAODpidUtil::MakeITSPID(AliAODTrack *track,Double_t *p) const
105 {
106   //
107   // ITS PID
108   //  1) Truncated mean method
109   //
110
111
112   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
113   UChar_t clumap=track->GetITSClusterMap();
114   Int_t nPointsForPid=0;
115   for(Int_t i=2; i<6; i++){
116    if(clumap&(1<<i)) ++nPointsForPid;
117   }
118   if(nPointsForPid<3) { // track not to be used for combined PID purposes
119     for (Int_t j=0; j<AliPID::kSPECIES; j++) 
120       p[j] = 1./AliPID::kSPECIES;
121     return;
122   }
123   Double_t mom=track->P();  
124   AliAODPid *pidObj = track->GetDetPid();
125
126   Double_t dedx = pidObj->GetITSsignal();
127   Bool_t mismatch = kTRUE;
128   Bool_t isSA = kTRUE;
129   if(track->GetStatus() & AliESDtrack::kTPCin){
130     isSA = kFALSE;
131     if (pidObj)
132       mom = pidObj->GetTPCmomentum();
133   }
134   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
135     Double_t mass = AliPID::ParticleMass(j);//GeV/c^2
136     Double_t bethe = fITSResponse.Bethe(mom,mass);
137     Double_t sigma = fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
138     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
139       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
140     } else {
141       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
142       mismatch=kFALSE;
143     }
144
145     // Check for particles heavier than (AliPID::kSPECIES - 1)
146
147   }
148
149   if (mismatch)
150     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
151
152   return;
153
154 }
155 //_________________________________________________________________________
156 void AliAODpidUtil::MakeTOFPID(AliAODTrack *track, Double_t *p) const
157 {
158   //
159   //   TOF PID using gaussian response
160   //
161   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return;
162   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return;
163   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return;
164
165   Double_t time[AliPID::kSPECIESN];
166   Double_t sigma[AliPID::kSPECIESN];
167   AliAODPid *pidObj = track->GetDetPid();
168   pidObj->GetIntegratedTimes(time);
169
170   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
171     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
172   }
173
174   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
175                   Form("Expected TOF signals [ps]: %f %f %f %f %f",
176                        time[AliPID::kElectron],
177                        time[AliPID::kMuon],
178                        time[AliPID::kPion],
179                        time[AliPID::kKaon],
180                        time[AliPID::kProton]));
181
182   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
183                   Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
184                        sigma[AliPID::kElectron],
185                        sigma[AliPID::kMuon],
186                        sigma[AliPID::kPion],
187                        sigma[AliPID::kKaon],
188                        sigma[AliPID::kProton]
189                        ));
190
191   Double_t tof = pidObj->GetTOFsignal();
192
193   Bool_t mismatch = kTRUE;
194   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
195     Double_t sig = sigma[j];
196     if (TMath::Abs(tof-time[j]) > fRange*sig) {
197       p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
198     } else
199       p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
200
201     // Check the mismatching
202     Double_t mass = AliPID::ParticleMass(j);
203     Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
204     if (p[j]>pm) mismatch = kFALSE;
205
206     // Check for particles heavier than (AliPID::kSPECIES - 1)
207
208   }
209
210   if (mismatch)
211     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
212
213   return;
214 }
215 //_________________________________________________________________________
216 void AliAODpidUtil::MakeTRDPID(AliAODTrack *track,Double_t *p) const
217 {
218   
219   // Method to recalculate the TRD PID probabilities
220   if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return;
221
222   AliAODPid *pidObj = track->GetDetPid();
223   Float_t *mom=pidObj->GetTRDmomentum();
224   Int_t ntracklets=0;
225   for(Int_t iPl=0;iPl<6;iPl++){
226    if(mom[iPl]>0.) ntracklets++;
227   }
228    if(ntracklets<4) return;
229
230   Double_t* dedx=pidObj->GetTRDsignal();
231   Bool_t norm=kTRUE;
232   fTRDResponse.GetResponse(pidObj->GetTRDnSlices(),dedx,mom,p,norm);
233   return;
234 }