]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEERBase/AliITSPIDResponse.cxx
added a new class/task for light nuclei
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliITSPIDResponse.cxx
CommitLineData
10d100d4 1/**************************************************************************
2 * Copyright(c) 2005-2007, 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$ */
17
18//-----------------------------------------------------------------
19// ITS PID method # 1
20// Implementation of the ITS PID class
21// Very naive one... Should be made better by the detector experts...
22// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
23//-----------------------------------------------------------------
24#include "TMath.h"
25#include "AliITSPIDResponse.h"
26#include "AliITSPidParams.h"
27#include "AliExternalTrackParam.h"
28
10d100d4 29ClassImp(AliITSPIDResponse)
30
15e979c9 31AliITSPIDResponse::AliITSPIDResponse(Bool_t isMC):
10d100d4 32 fRes(0.13),
33 fKp1(15.77),
34 fKp2(4.95),
35 fKp3(0.312),
36 fKp4(2.14),
37 fKp5(0.82)
38{
15e979c9 39 if(!isMC){
40 fBBtpcits[0]=0.73;
41 fBBtpcits[1]=14.68;
42 fBBtpcits[2]=0.905;
43 fBBtpcits[3]=1.2;
44 fBBtpcits[4]=6.6;
62ccfebf 45 fBBsa[0]=2.73198E7; //pure PHOBOS parameterization
88f46717 46 fBBsa[1]=6.92389;
47 fBBsa[2]=1.90088E-6;
48 fBBsa[3]=1.90088E-6;
49 fBBsa[4]=3.40644E-7;
62ccfebf 50 fBBsaHybrid[0]=1.43505E7; //PHOBOS+Polinomial parameterization
51 fBBsaHybrid[1]=49.3402;
52 fBBsaHybrid[2]=1.77741E-7;
53 fBBsaHybrid[3]=1.77741E-7;
54 fBBsaHybrid[4]=1.01311E-7;
55 fBBsaHybrid[5]=77.2777;
56 fBBsaHybrid[6]=33.4099;
57 fBBsaHybrid[7]=46.0089;
58 fBBsaHybrid[8]=-2.26583;
59 fBBsaElectron[0]=4.05799E6; //electrons in the ITS
88f46717 60 fBBsaElectron[1]=38.5713;
61 fBBsaElectron[2]=1.46462E-7;
62 fBBsaElectron[3]=1.46462E-7;
63 fBBsaElectron[4]=4.40284E-7;
8abeb05b 64 fResolSA[0]=1.; // 0 cluster tracks should not be used
88f46717 65 fResolSA[1]=0.25; // rough values for tracks with 1
66 fResolSA[2]=0.131; // value from pp 2010 run (L. Milano, 16-Jun-11)
67 fResolSA[3]=0.113; // value from pp 2010 run
8abeb05b 68 fResolSA[4]=0.104; // value from pp 2010 run
15e979c9 69 for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
70 }else{
99daa709 71 fBBtpcits[0]=1.04;
72 fBBtpcits[1]=27.14;
73 fBBtpcits[2]=1.00;
74 fBBtpcits[3]=0.964;
75 fBBtpcits[4]=2.59;
62ccfebf 76 fBBsa[0]=2.02078E7; //pure PHOBOS parameterization
88f46717 77 fBBsa[1]=14.0724;
78 fBBsa[2]=3.84454E-7;
79 fBBsa[3]=3.84454E-7;
80 fBBsa[4]=2.43913E-7;
62ccfebf 81 fBBsaHybrid[0]=1.05381E7; //PHOBOS+Polinomial parameterization
82 fBBsaHybrid[1]=89.3933;
83 fBBsaHybrid[2]=2.4831E-7;
84 fBBsaHybrid[3]=2.4831E-7;
85 fBBsaHybrid[4]=7.80591E-8;
86 fBBsaHybrid[5]=62.9214;
87 fBBsaHybrid[6]=32.347;
88 fBBsaHybrid[7]=58.7661;
89 fBBsaHybrid[8]=-3.39869;
90 fBBsaElectron[0]=2.26807E6; //electrons in the ITS
88f46717 91 fBBsaElectron[1]=99.985;
92 fBBsaElectron[2]=0.000714841;
93 fBBsaElectron[3]=0.000259585;
94 fBBsaElectron[4]=1.39412E-7;
8abeb05b 95 fResolSA[0]=1.; // 0 cluster tracks should not be used
88f46717 96 fResolSA[1]=0.25; // rough values for tracks with 1
97 fResolSA[2]=0.126; // value from pp 2010 simulations (L. Milano, 16-Jun-11)
98 fResolSA[3]=0.109; // value from pp 2010 simulations
99 fResolSA[4]=0.097; // value from pp 2010 simulations
15e979c9 100 for(Int_t i=0; i<5;i++) fResolTPCITS[i]=0.13;
101 }
10d100d4 102}
103
56576f1e 104/*
10d100d4 105//_________________________________________________________________________
106AliITSPIDResponse::AliITSPIDResponse(Double_t *param):
9ebbddd4 107 fRes(param[0]),
10d100d4 108 fKp1(15.77),
109 fKp2(4.95),
110 fKp3(0.312),
111 fKp4(2.14),
112 fKp5(0.82)
113{
114 //
115 // The main constructor
116 //
6b4634a4 117 for (Int_t i=0; i<5;i++) {
118 fBBsa[i]=0.;
119 fBBtpcits[i]=0.;
120 fResolSA[i]=0.;
121 fResolTPCITS[i]=0.;
122 }
10d100d4 123}
56576f1e 124*/
10d100d4 125
8abeb05b 126//_________________________________________________________________________
15e979c9 127Double_t AliITSPIDResponse::BetheAleph(Double_t p, Double_t mass) const {
10d100d4 128 //
129 // returns AliExternalTrackParam::BetheBloch normalized to
130 // fgMIP at the minimum
131 //
15e979c9 132
10d100d4 133 Double_t bb=
134 AliExternalTrackParam::BetheBlochAleph(p/mass,fKp1,fKp2,fKp3,fKp4,fKp5);
9ebbddd4 135 return bb;
10d100d4 136}
137
8abeb05b 138//_________________________________________________________________________
15e979c9 139Double_t AliITSPIDResponse::Bethe(Double_t p, Double_t mass, Bool_t isSA) const {
140 //
141 // returns AliExternalTrackParam::BetheBloch normalized to
142 // fgMIP at the minimum
143 //
144
145 Double_t bg=p/mass;
146 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
147 Double_t gamma=bg/beta;
62ccfebf 148 Double_t bb=1.;
149
15e979c9 150 Double_t par[5];
151 if(isSA){
88f46717 152 if(mass>0.0005 && mass<0.00052){
153 //if is an electron use a specific BB parameterization
154 //To be used only between 100 and 160 MeV/c
155 for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
156 }else{
157 for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsa[ip];
158 }
15e979c9 159 }else{
160 for(Int_t ip=0; ip<5;ip++) par[ip]=fBBtpcits[ip];
161 }
162 Double_t eff=1.0;
163 if(bg<par[2])
164 eff=(bg-par[3])*(bg-par[3])+par[4];
165 else
166 eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
88f46717 167
15e979c9 168 if(gamma>=0. && beta>0.){
169 bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
170 }
171 return bb;
172}
173
62ccfebf 174//_________________________________________________________________________
175Double_t AliITSPIDResponse::BetheITSsaHybrid(Double_t p, Double_t mass) const {
176 //
177 // returns AliExternalTrackParam::BetheBloch normalized to
178 // fgMIP at the minimum. The PHOBOS parameterization is used for beta*gamma>0.76.
179 // For beta*gamma<0.76 a polinomial function is used
180
181 Double_t bg=p/mass;
182 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
183 Double_t gamma=bg/beta;
184 Double_t bb=1.;
185
186 Double_t par[9];
187 //parameters for pi, K, p
188 for(Int_t ip=0; ip<9;ip++) par[ip]=fBBsaHybrid[ip];
189 //if it is an electron the PHOBOS part of the parameterization is tuned for e
190 //in the range used for identification beta*gamma is >0.76 for electrons
191 //To be used only between 100 and 160 MeV/c
192 if(mass>0.0005 && mass<0.00052)for(Int_t ip=0; ip<5;ip++) par[ip]=fBBsaElectron[ip];
193
194 if(gamma>=0. && beta>0. && bg>0.1){
195 if(bg>0.76){//PHOBOS
196 Double_t eff=1.0;
197 if(bg<par[2])
198 eff=(bg-par[3])*(bg-par[3])+par[4];
199 else
200 eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
201
202 bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
203 }else{//Polinomial
204 bb=par[5] + par[6]/bg + par[7]/(bg*bg) + par[8]/(bg*bg*bg);
205 }
206 }
207 return bb;
208}
209
8abeb05b 210//_________________________________________________________________________
15e979c9 211Double_t AliITSPIDResponse::GetResolution(Double_t bethe,
212 Int_t nPtsForPid,
213 Bool_t isSA) const {
10d100d4 214 //
215 // Calculate expected resolution for truncated mean
216 //
15e979c9 217 Float_t r;
218 if(isSA) r=fResolSA[nPtsForPid];
219 else r=fResolTPCITS[nPtsForPid];
220 return r*bethe;
10d100d4 221}
222
15e979c9 223
224
225
8abeb05b 226//_________________________________________________________________________
b52bfc67 227void AliITSPIDResponse::GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES], Bool_t isMC) const {
10d100d4 228 //
229 // Method to calculate PID probabilities for a single track
230 // using the likelihood method
231 //
232 const Int_t nLay = 4;
2ca1f4ee 233 const Int_t nPart= 4;
10d100d4 234
b52bfc67 235 static AliITSPidParams pars(isMC); // Pid parametrisation parameters
10d100d4 236
2ca1f4ee 237 Double_t itsProb[nPart] = {1,1,1,1}; // e, p, K, pi
10d100d4 238
239 for (Int_t iLay = 0; iLay < nLay; iLay++) {
2ca1f4ee 240 if (qclu[iLay] <= 50.)
10d100d4 241 continue;
242
243 Float_t dedx = qclu[iLay];
244 Float_t layProb = pars.GetLandauGausNorm(dedx,AliPID::kProton,mom,iLay+3);
245 itsProb[0] *= layProb;
246
247 layProb = pars.GetLandauGausNorm(dedx,AliPID::kKaon,mom,iLay+3);
10d100d4 248 itsProb[1] *= layProb;
249
250 layProb = pars.GetLandauGausNorm(dedx,AliPID::kPion,mom,iLay+3);
251 itsProb[2] *= layProb;
2ca1f4ee 252
253 layProb = pars.GetLandauGausNorm(dedx,AliPID::kElectron,mom,iLay+3);
254 itsProb[3] *= layProb;
10d100d4 255 }
256
257 // Normalise probabilities
258 Double_t sumProb = 0;
259 for (Int_t iPart = 0; iPart < nPart; iPart++) {
260 sumProb += itsProb[iPart];
261 }
2ca1f4ee 262 sumProb += itsProb[2]; // muon cannot be distinguished from pions
10d100d4 263
264 for (Int_t iPart = 0; iPart < nPart; iPart++) {
265 itsProb[iPart]/=sumProb;
266 }
2ca1f4ee 267 condprobfun[AliPID::kElectron] = itsProb[3];
897a0e31 268 condprobfun[AliPID::kMuon] = itsProb[2];
269 condprobfun[AliPID::kPion] = itsProb[2];
10d100d4 270 condprobfun[AliPID::kKaon] = itsProb[1];
271 condprobfun[AliPID::kProton] = itsProb[0];
272 return;
273}
15e979c9 274
8abeb05b 275//_________________________________________________________________________
276Int_t AliITSPIDResponse::GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA) const{
277 // method to get particle identity with simple cuts on dE/dx vs. momentum
278
279 Double_t massp=AliPID::ParticleMass(AliPID::kProton);
280 Double_t massk=AliPID::ParticleMass(AliPID::kKaon);
281 Double_t bethep=Bethe(mom,massp,isSA);
282 Double_t bethek=Bethe(mom,massk,isSA);
283 if(signal>(0.5*(bethep+bethek))) return AliPID::kProton;
284 Double_t masspi=AliPID::ParticleMass(AliPID::kPion);
285 Double_t bethepi=Bethe(mom,masspi,isSA);
286 if(signal>(0.5*(bethepi+bethek))) return AliPID::kKaon;
287 return AliPID::kPion;
288
289}