1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////
20 // Implementation of the class to store parameters of ITS //
21 // response funcions for dE/dx based PID //
22 // Origin: F.Prino, Torino, prino@to.infn.it //
24 ///////////////////////////////////////////////////////////////////
29 #include "AliITSPidParams.h"
31 ClassImp(AliITSPidParams)
33 //______________________________________________________________________
34 AliITSPidParams::AliITSPidParams(Bool_t isMC):
36 fSDDElecLandauWidth(0),
37 fSDDElecGaussWidth(0),
38 fSSDElecLandauWidth(0),
39 fSSDElecGaussWidth(0),
40 fSDDPionLandauWidth(0),
41 fSDDPionGaussWidth(0),
42 fSSDPionLandauWidth(0),
43 fSSDPionGaussWidth(0),
44 fSDDKaonLandauWidth(0),
45 fSDDKaonGaussWidth(0),
46 fSSDKaonLandauWidth(0),
47 fSSDKaonGaussWidth(0),
48 fSDDProtLandauWidth(0),
49 fSDDProtGaussWidth(0),
50 fSSDProtLandauWidth(0),
53 // default constructor
57 //______________________________________________________________________
58 AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
60 fSDDElecLandauWidth(0),
61 fSDDElecGaussWidth(0),
62 fSSDElecLandauWidth(0),
63 fSSDElecGaussWidth(0),
64 fSDDPionLandauWidth(0),
65 fSDDPionGaussWidth(0),
66 fSSDPionLandauWidth(0),
67 fSSDPionGaussWidth(0),
68 fSDDKaonLandauWidth(0),
69 fSDDKaonGaussWidth(0),
70 fSSDKaonLandauWidth(0),
71 fSSDKaonGaussWidth(0),
72 fSDDProtLandauWidth(0),
73 fSDDProtGaussWidth(0),
74 fSSDProtLandauWidth(0),
77 // standard constructor
81 //______________________________________________________________________
82 AliITSPidParams::~AliITSPidParams(){
84 if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
85 if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
87 if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
88 if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
90 if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
91 if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
93 if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
94 if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
96 if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
97 if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
99 if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
100 if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
102 if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
103 if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
105 if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
106 if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
108 //______________________________________________________________________
109 Double_t AliITSPidParams::BetheBloch(Double_t mom, Double_t mass, const Double_t *p) const{
111 // This is the empirical parameterization of the Bethe-Bloch formula.
112 // It is normalized to 1 at the minimum.
116 Double_t bg = mom/mass;
117 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
118 Double_t aa = TMath::Power(beta,p[3]);
119 Double_t bb = TMath::Power(1./bg,p[4]);
120 bb=TMath::Log(p[2]+bb);
121 return (p[0]-aa-bb)*p[1]/aa;
123 //______________________________________________________________________
124 Double_t AliITSPidParams::ExtrapolateWidth(Double_t mom, Double_t x1, Double_t y1,Double_t x2, Double_t y2) const{
126 // This is a linear extrapolation of Landau width and Gaussian width
129 Double_t slope = (y2-y1)/(x2-x1);
130 return slope*mom+(y1-slope*x1);
131 }//______________________________________________________________________
132 void AliITSPidParams::InitMC(){
133 // initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
134 // parameter values from LHC10d1
135 // MPV BetheBloch parameters;
137 //sdd MC electrons parameters
138 fSDDElecMPVBetheParams[0] = -0.0931934;
139 fSDDElecMPVBetheParams[1] = 77.8422;
140 fSDDElecMPVBetheParams[2] = -0.889085;
141 fSDDElecMPVBetheParams[3] = -154.455;
142 fSDDElecMPVBetheParams[4] = -0.000256076;
144 //ssd MC electrons parameters
145 fSSDElecMPVBetheParams[0] = -0.0989358;
146 fSSDElecMPVBetheParams[1] = 77.8271;
147 fSSDElecMPVBetheParams[2] = -0.900887;
148 fSSDElecMPVBetheParams[3] =-1241.14;
149 fSSDElecMPVBetheParams[4] = -0.0014204;
152 if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
153 fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
154 fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);
156 if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
157 fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
158 fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);
160 if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
161 fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
162 fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);
164 if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
165 fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
166 fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
168 //sdd MC hadrons parameters
169 fSDDHadronMPVBetheParams[0] = 1.13531;
170 fSDDHadronMPVBetheParams[1] = -156.651;
171 fSDDHadronMPVBetheParams[2] = 1.87562;
172 fSDDHadronMPVBetheParams[3] = 0.45819;
173 fSDDHadronMPVBetheParams[4] = 2.26506;
175 //ssd MC hadrons parameters
176 fSSDHadronMPVBetheParams[0] = -0.451908;
177 fSSDHadronMPVBetheParams[1] = -55.4368;
178 fSSDHadronMPVBetheParams[2] = 0.984636;
179 fSSDHadronMPVBetheParams[3] = 0.97078;
180 fSSDHadronMPVBetheParams[4] = 2.50883;
183 if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
184 fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
185 fSDDPionLandauWidth->SetParameters(0.08026, 5.87922);
187 if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
188 fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
189 fSDDPionGaussWidth->SetParameters(-0.090495, 7.705286);
191 if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
192 fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
193 fSSDPionLandauWidth->SetParameters(0.083882, 5.823419);
195 if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
196 fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
197 fSSDPionGaussWidth->SetParameters(-0.105218, 5.650956);
200 if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
201 fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
202 fSDDKaonLandauWidth->SetParameters(1.430692, 5.581389);
204 if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
205 fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
206 fSDDKaonGaussWidth->SetParameters(-1.6537, 8.071832);
208 if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
209 fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
210 fSSDKaonLandauWidth->SetParameters(1.368824, 5.639291);
212 if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
213 fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
214 fSSDKaonGaussWidth->SetParameters(-1.901858, 6.0932281);
217 if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
218 fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
219 fSDDProtLandauWidth->SetParameters(6.529418, 5.049098);
221 if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
222 fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
223 fSDDProtGaussWidth->SetParameters(-9.360599, 11.318026);
225 if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
226 fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
227 fSSDProtLandauWidth->SetParameters(6.419493, 4.925070);
229 if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
230 fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
231 fSSDProtGaussWidth->SetParameters(-9.599064, 9.358656);
233 //______________________________________________________________________
234 void AliITSPidParams::InitData(){
235 // initialize TFormulas to Real Data values (=p-p ALICE Experiment)
236 // parameter values from LHC10b
237 // MPV BetheBloch parameters;
239 //sdd data electrons parameters
240 fSDDElecMPVBetheParams[0] = -0.130204;
241 fSDDElecMPVBetheParams[1] = 75.1267;
242 fSDDElecMPVBetheParams[2] = -0.889337;
243 fSDDElecMPVBetheParams[3] = 388.372;
244 fSDDElecMPVBetheParams[4] = 0.00134649;
246 //ssd data electrons parameters
247 fSSDElecMPVBetheParams[0] = -0.162773;
248 fSSDElecMPVBetheParams[1] = 72.9393;
249 fSSDElecMPVBetheParams[2] = -0.896944;
250 fSSDElecMPVBetheParams[3] = 3233.02;
251 fSSDElecMPVBetheParams[4] = 0.00146896;
254 if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
255 fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
256 fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);
258 if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
259 fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
260 fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);
262 if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
263 fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
264 fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);
266 if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
267 fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
268 fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
270 //sdd data hadrons parameters
271 fSDDHadronMPVBetheParams[0] = -18.1867;
272 fSDDHadronMPVBetheParams[1] = -3.45806;
273 fSDDHadronMPVBetheParams[2] = 2.23635;
274 fSDDHadronMPVBetheParams[3] = 2.08328;
275 fSDDHadronMPVBetheParams[4] = -1.92331;
277 //ssd data hadrons parameters
278 fSSDHadronMPVBetheParams[0] = -12.6459;
279 fSSDHadronMPVBetheParams[1] = -4.84598;
280 fSSDHadronMPVBetheParams[2] = 1.50253;
281 fSSDHadronMPVBetheParams[3] = 2.08328;
282 fSSDHadronMPVBetheParams[4] = -1.3719;
285 if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
286 fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
287 fSDDPionLandauWidth->SetParameters(0.07694, 5.468728);
289 if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
290 fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
291 fSDDPionGaussWidth->SetParameters(-0.098209, 10.409441);
293 if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
294 fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
295 fSSDPionLandauWidth->SetParameters(0.071602, 5.365442);
297 if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
298 fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
299 fSSDPionGaussWidth->SetParameters(-0.098045, 7.617583);
302 if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
303 fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
304 fSDDKaonLandauWidth->SetParameters(0.998191, 5.461668);
306 if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
307 fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
308 fSDDKaonGaussWidth->SetParameters(1.629308, 9.851873);
310 if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
311 fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
312 fSSDKaonLandauWidth->SetParameters(0.773113, 5.618683);
314 if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
315 fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
316 fSSDKaonGaussWidth->SetParameters(1.510713, 6.862774);
319 if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
320 fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
321 fSDDProtLandauWidth->SetParameters(3.561429, 5.372105);
323 if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
324 fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
325 fSDDProtGaussWidth->SetParameters(5.395926, 10.044613);
327 if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
328 fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
329 fSSDProtLandauWidth->SetParameters(2.647428, 5.678460);
331 if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
332 fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
333 fSSDProtGaussWidth->SetParameters(4.91025, 6.779763);
335 //_______________________________________________________________________
336 Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
337 // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
338 if(TMath::Abs(pdgCode)==11) return GetLandauGausNorm(dedx,AliPID::kElectron,mom,lay);
339 else if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
340 else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
341 else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
344 //_______________________________________________________________________
345 Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
346 // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
350 if(partType==AliPID::kElectron){
351 if(lay==3 || lay==4){
352 par[0]=GetSDDElecLandauWidth(mom);
353 par[1]=GetSDDElecMPV(mom);
354 par[2]=GetSDDElecGaussWidth(mom);
357 else if(lay==5 || lay==6){
358 par[0]=GetSSDElecLandauWidth(mom);
359 par[1]=GetSSDElecMPV(mom);
360 par[2]=GetSSDElecGaussWidth(mom);
363 }else if(partType==AliPID::kPion){
364 if(lay==3 || lay==4){
365 par[0]=GetSDDPionLandauWidth(mom);
366 par[1]=GetSDDPionMPV(mom);
367 par[2]=GetSDDPionGaussWidth(mom);
370 else if(lay==5 || lay==6){
371 par[0]=GetSSDPionLandauWidth(mom);
372 par[1]=GetSSDPionMPV(mom);
373 par[2]=GetSSDPionGaussWidth(mom);
376 }else if(partType==AliPID::kKaon){
377 if(lay==3 || lay==4){
378 par[0]=GetSDDKaonLandauWidth(mom);
379 par[1]=GetSDDKaonMPV(mom);
380 par[2]=GetSDDKaonGaussWidth(mom);
383 else if(lay==5 || lay==6){
384 par[0]=GetSSDKaonLandauWidth(mom);
385 par[1]=GetSSDKaonMPV(mom);
386 par[2]=GetSSDKaonGaussWidth(mom);
389 }else if(partType==AliPID::kProton){
390 if(lay==3 || lay==4){
391 par[0]=GetSDDProtLandauWidth(mom);
392 par[1]=GetSDDProtMPV(mom);
393 par[2]=GetSDDProtGaussWidth(mom);
396 else if(lay==5 || lay==6){
397 par[0]=GetSSDProtLandauWidth(mom);
398 par[1]=GetSSDProtMPV(mom);
399 par[2]=GetSSDProtGaussWidth(mom);
403 if(!isSet) return 0.;
405 const Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
406 const Double_t mpshift = -0.22278298; // Landau maximum location
408 const Double_t np = 100.0; // number of convolution steps
409 const Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
419 // MP shift correction
420 mpc = par[1] - mpshift * par[0];
421 // Range of convolution integral
422 xlow = dedx - sc * par[2];
423 xupp = dedx + sc * par[2];
424 step = (xupp-xlow) / np;
426 // Convolution integral of Landau and Gaussian by sum
427 for(i=1.0; i<=np/2; i++) {
428 xx = xlow + (i-.5) * step;
430 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
431 sum += fland * TMath::Gaus(dedx,xx,par[2]);
433 xx = xupp - (i-.5) * step;
434 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
435 sum += fland * TMath::Gaus(dedx,xx,par[2]);
438 return (step * sum * invsq2pi / par[2]);