]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEERBase/AliEMCALPIDResponse.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New directory
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEMCALPIDResponse.cxx
CommitLineData
fdb4683b 1/**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15//////////////////////////////////////////////////////////////////////////\r
16// //\r
17// AliEMCALPIDResponse //\r
18// //\r
19// EMCAL class to perfom PID //\r
20// This is a prototype and still under development //\r
21// //\r
22// ---------------------------------------------------------------------//\r
23// GetNumberOfSigmas(): //\r
24// //\r
25// Electrons: Number of Sigmas for E/p value //\r
26// Parametrization of LHC11a (after recalibration) // \r
27// //\r
28// NON electrons: //\r
29// Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //\r
30// --> return +/- 999 //\r
31// Otherwise //\r
32// --> return nsigma (parametrization of LHC10e) // \r
33// //\r
34// ---------------------------------------------------------------------//\r
35// ComputeEMCALProbability(): //\r
36// //\r
37// Electrons: Probability from Gaussian distribution //\r
38// //\r
39// NON electrons: //\r
40// Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //\r
41// --> probability to find particles below or above thr. //\r
42// Otherwise //\r
43// --> Probability from Gaussian distribution // \r
44// (proper normalization to each other?) //\r
45// //\r
46//////////////////////////////////////////////////////////////////////////\r
47\r
24ed4170 48#include <TF1.h>\r
49#include <TMath.h>\r
fdb4683b 50\r
51#include "AliEMCALPIDResponse.h" //class header\r
52\r
fdb4683b 53#include "AliLog.h" \r
54\r
55ClassImp(AliEMCALPIDResponse)\r
9ec6e1da 56\r
57const Float_t AliEMCALPIDResponse::fLowEoP = 0.5; // lower E/p threshold for NON electrons\r
58const Float_t AliEMCALPIDResponse::fHighEoP = 1.5; // upper E/p threshold for NON electrons\r
59\r
fdb4683b 60//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
61AliEMCALPIDResponse::AliEMCALPIDResponse():\r
7cd9f4fb 62 TObject(),\r
63 fNorm(NULL)\r
fdb4683b 64{\r
65 //\r
66 // The default constructor\r
67 //\r
68\r
69 for(Int_t i = 0; i < fNptBins; i++){\r
70\r
71 fPtCutMin[i] = 0.0;\r
72\r
73 for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++){\r
74\r
75 fMeanEoP[j][i] = 0.0;\r
76 fSigmaEoP[j][i] = 0.0;\r
77 fProbLow[j][i] = 0.0;\r
78 fProbHigh[j][i] = 0.0;\r
79\r
80 }\r
81 } \r
82 fPtCutMin[fNptBins] = 0.0;\r
83\r
84 fNorm = new TF1("fNorm","gaus",-20,20); \r
85\r
86 SetPtBoundary();\r
87 SetParametrizations();\r
88}\r
9a2d18d5 89//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
e96b9916 90AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):\r
7cd9f4fb 91 TObject(other),\r
92 fNorm(other.fNorm)\r
e96b9916 93{\r
94 //\r
95 // The copy constructor\r
96 //\r
97 for(Int_t i = 0; i < fNptBins; i++)\r
98 {\r
99 fPtCutMin[i] = 0.0;\r
100 for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++)\r
101 {\r
102 fMeanEoP[j][i] = 0.0;\r
103 fSigmaEoP[j][i] = 0.0;\r
104 fProbLow[j][i] = 0.0;\r
105 fProbHigh[j][i] = 0.0;\r
106 }\r
107 }\r
108 \r
109 fPtCutMin[fNptBins] = 0.0;\r
110 SetPtBoundary();\r
111 SetParametrizations();\r
112}\r
9a2d18d5 113//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
e96b9916 114AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)\r
115{\r
116 //\r
117 // The assignment operator\r
118 //\r
7cd9f4fb 119\r
120 if(this == &other) return *this;\r
121 \r
122 // Make copy\r
123 TObject::operator=(other);\r
124 fNorm = other.fNorm;\r
125\r
126 for(Int_t i = 0; i < fNptBins; i++)\r
e96b9916 127 {\r
128 fPtCutMin[i] = 0.0;\r
129 for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++)\r
130 {\r
131 fMeanEoP[j][i] = 0.0;\r
132 fSigmaEoP[j][i] = 0.0;\r
133 fProbLow[j][i] = 0.0;\r
134 fProbHigh[j][i] = 0.0;\r
135 }\r
136 }\r
7cd9f4fb 137 \r
138 fPtCutMin[fNptBins] = 0.0;\r
139 SetPtBoundary();\r
140 SetParametrizations();\r
141\r
142 return *this;\r
e96b9916 143}\r
9a2d18d5 144//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
145AliEMCALPIDResponse::~AliEMCALPIDResponse() {\r
e96b9916 146\r
9a2d18d5 147 delete fNorm;\r
148\r
149}\r
fdb4683b 150//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
151void AliEMCALPIDResponse::SetPtBoundary(){\r
152 //\r
153 // Set boundaries for momentum bins\r
154 //\r
155 fPtCutMin[0] = 1.5;\r
156 fPtCutMin[1] = 2.5;\r
157 fPtCutMin[2] = 3.5;\r
158 fPtCutMin[3] = 4.5;\r
159 fPtCutMin[4] = 5.5;\r
160 fPtCutMin[5] = 6.5;\r
161\r
162}\r
163//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
164void AliEMCALPIDResponse::SetParametrizations(){\r
165\r
166// This are the preliminary parametrizations (hard coded)\r
167// For electrons from LHC11a (Deepa Thomas)\r
168// For NON-electrons from LHC10e (TOF/TPC analysis)\r
169 \r
170 // Gaussian mean\r
171 Float_t mean[4][6] = {\r
172 { 0.932, 0.997, 0.998, 1.001, 1.011, 1.011 }, // electrons\r
173 { 0.227804, 0.34839, 0.404077, -0.107795, -4.14584, 0.5 }, // NON electrons\r
174 { -2.10377, 0.0582898, 0.0582898, 0.0582898, 0.0582898, 0.0582898 }, // protons\r
175 { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5} // anti-protons\r
176 }; \r
177 \r
178 // Gaussian sigma\r
179 Float_t sigma[4][6]= {\r
180 { 0.0866, 0.0693, 0.0664, 0.0583, 0.0488, 0.0515 }, // electrons\r
181 { 0.310831, 0.267586, 0.404077, 0.381968, 1.46183, 0.314687 }, // NON electrons\r
182 { 0.603209, 0.255332, 0.255332, 0.255332, 0.255332, 0.255332}, // protons\r
183 { 0.516837, 0.351516,0.351516,0.351516,0.351516,0.351516 } // anti - protons\r
184 };\r
185\r
186 // lower probability\r
187 Float_t probL[3][6] = {\r
188 { 0.928689, 0.938455, 0.940448, 0.948496, 0.955981, 0.951923 }, // NON electrons\r
189 { 0.974518, 0.978088, 0.975089, 0.975089, 0.975089,0.975089}, // protons\r
190 { 0.824037, 0.861149, 0.898734, 0.898734, 0.898734, 0.898734}, // anti - protons\r
191 };\r
192\r
193 // upper probability\r
194 Float_t probH[3][6] = {\r
195 { 0.00030227, 4.04106e-05, 0.000147406, 0., 0.000956938, 0.00106838 }, // NON electrons\r
196 { 0.000157945, 0., 0., 0., 0., 0. }, // protons\r
197 { 0.00343237, 0., 0., 0., 0., 0.} // anti - protons\r
198 };\r
199 \r
200 \r
201 // set parametrizations\r
202 Int_t spec = 0;\r
203 for (Int_t species = 0; species < 2*AliPID::kSPECIES; species++) { // first negative particles and then positive\r
204 for (Int_t pt = 0; pt < fNptBins; pt++){\r
205\r
206 switch(species){\r
207 case 0: // electrons\r
208 spec = 0;\r
209 break;\r
210 case 4: // anti - protons\r
211 spec = 3;\r
212 break;\r
213 case 5: // positrons\r
214 spec = 0;\r
215 break;\r
216 case 9: // protons\r
217 spec = 2;\r
218 break;\r
219 default: // NON electrons\r
220 spec = 1;\r
221 break;\r
222 }\r
223 \r
224\r
225 fMeanEoP[species][pt] = mean[spec][pt]; \r
226 fSigmaEoP[species][pt] = sigma[spec][pt]; \r
227 if( spec == 0) { // electrons have NO lower and upper probability thresholds --> set to 0\r
228 fProbLow[species][pt] = 0.;\r
229 fProbHigh[species][pt] = 0.; \r
230 }\r
231 else{\r
232 fProbLow[species][pt] = probL[spec-1][pt];\r
233 fProbHigh[species][pt] = probH[spec-1][pt]; \r
234 } \r
235 \r
236 }//loop pt bins\r
237 }//loop species\r
238}\r
239//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
240Int_t AliEMCALPIDResponse::GetPtBin(Float_t pt) const {\r
241 //\r
242 // Returns the momentum bin index\r
243 //\r
244\r
245 Int_t i = -1;\r
246 while(pt > fPtCutMin[i+1] && i+1 < fNptBins) i++;\r
247\r
248 return i;\r
249}\r
250\r
251//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
252Double_t AliEMCALPIDResponse::GetExpectedSignal( Float_t pt, AliPID::EParticleType n, Int_t charge) const {\r
253 //\r
254 // Calculates the expected PID signal as the function of \r
255 // the information stored in the track, for the specified particle type \r
256 // \r
257\r
258 Double_t signal = 0.;\r
259\r
260 // Check the charge\r
261 if( charge != -1 && charge != 1){\r
262 \r
263 return signal;\r
264 }\r
265 \r
266 // Get the pt bin\r
267 Int_t ptBin = GetPtBin(pt);\r
268 \r
269 // Get the species (first negative , then positive)\r
270 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
271\r
272 // Get the signal\r
273 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
274 signal = fMeanEoP[species][ptBin];\r
275 }\r
276\r
277 return signal;\r
278\r
279}\r
280//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
281Double_t AliEMCALPIDResponse::GetExpectedSigma( Float_t pt, AliPID::EParticleType n, Int_t charge) const {\r
282 //\r
283 // Calculates the expected sigma of the PID signal as the function of \r
284 // the information stored in the track, for the specified particle type \r
285 // \r
286 //\r
287 \r
288 Double_t sigma = 999.;\r
289\r
290 // Check the charge\r
291 if( charge != -1 && charge != 1){\r
292 \r
293 return sigma;\r
294 } \r
295\r
296 // Get the pt bin\r
297 Int_t ptBin = GetPtBin(pt);\r
298 \r
299 // Get the species (first negative , then positive)\r
300 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
301\r
302 // Get the sigma\r
303 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
304 sigma = fSigmaEoP[species][ptBin];\r
305 }\r
306\r
307 return sigma;\r
308 \r
309}\r
310//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
311Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const {\r
312 //\r
313 // Calculates the expected sigma of the PID signal as the function of \r
314 // the information stored in the track, for the specified particle type \r
315 // \r
316 //\r
317 \r
318 Double_t norm = 1.;\r
319\r
320 // Check the charge\r
321 if( charge != -1 && charge != 1){\r
322 \r
323 return norm;\r
324 }\r
325\r
326 // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )\r
327 fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*GetExpectedSigma(pt,n,charge)*GetExpectedSigma(pt,n,charge)),GetExpectedSignal(pt,n,charge),GetExpectedSigma(pt,n,charge));\r
328 norm = 1./fNorm->Integral(fLowEoP,fHighEoP)*(1-GetLowProb(pt,n,charge)-GetHighProb(pt,n,charge));\r
329\r
330 return norm;\r
331}\r
332//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
333Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const {\r
334 \r
335 Double_t mean = GetExpectedSignal(pt,n,charge);\r
336 Double_t sigma = GetExpectedSigma(pt,n,charge);\r
337\r
338 // if electron\r
339 if(n == AliPID::kElectron){\r
340 return (eop - mean) / sigma;\r
341 }\r
342\r
343 // if NON electron\r
344 else{\r
345 if ( eop < fLowEoP )\r
346 return -999.; // not parametrized \r
347 else if ( eop > fHighEoP )\r
348 return 999.; // not parametrized \r
349 else{\r
350 return (eop - mean) / sigma; \r
351 }\r
352 }\r
353}\r
354//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
355Double_t AliEMCALPIDResponse::GetLowProb( Float_t pt, AliPID::EParticleType n, Int_t charge) const {\r
356 //\r
357 //\r
358 \r
359 Double_t prob = 0.;\r
360\r
361 // Check the charge\r
362 if( charge != -1 && charge != 1){\r
363 \r
364 return prob;\r
365 }\r
366 \r
367 // Get the pt bin\r
368 Int_t ptBin = GetPtBin(pt);\r
369 \r
370 // Get the species (first negative , then positive)\r
371 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
372\r
373 // Get the probability\r
374 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
375 prob = fProbLow[species][ptBin];\r
376 }\r
377\r
378 return prob;\r
379 \r
380}\r
381//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
382Double_t AliEMCALPIDResponse::GetHighProb( Float_t pt, AliPID::EParticleType n, Int_t charge) const {\r
383 //\r
384 //\r
385 \r
386 Double_t prob = 0.;\r
387\r
388 // Check the charge\r
389 if( charge != -1 && charge != 1){\r
390 \r
391 return prob;\r
392 }\r
393 \r
394 // Get the pt bin\r
395 Int_t ptBin = GetPtBin(pt);\r
396 \r
397 // Get the species (first negative , then positive)\r
398 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
399\r
400 // Get the probability\r
401 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
402 prob = fProbHigh[species][ptBin];\r
403 }\r
404\r
405 return prob;\r
406 \r
407}\r
408//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
409Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {\r
410 //\r
411 //\r
412 Double_t fRange = 5.0; // hardcoded \r
413 Double_t nsigma = 0.0;\r
414 Double_t sigma = 0.0;\r
e96b9916 415 Double_t proba = 999.;\r
416 \r
fdb4683b 417\r
418 // Check the charge\r
419 if( charge != -1 && charge != 1){\r
420 \r
421 return proba;\r
422 }\r
423\r
424 \r
425 // default value (will be returned, if pt below threshold)\r
426 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
e96b9916 427 pEMCAL[species] = 999.;\r
fdb4683b 428 }\r
429\r
430 if( GetPtBin(pt) > -1 ){\r
431\r
432 // set E/p range\r
433 if(eop < 0.05) eop = 0.05;\r
434 if(eop > 2.00) eop = 2.00;\r
435\r
436 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
437\r
438 AliPID::EParticleType type = AliPID::EParticleType(species);\r
439\r
440 // get nsigma value for each particle type at this E/p value\r
441 nsigma = GetNumberOfSigmas(pt,eop,type,charge);\r
442 sigma = GetExpectedSigma(pt,type,charge);\r
443\r
444 // electrons (standard Gaussian calculation of probabilities)\r
445 if(type == AliPID::kElectron){\r
446 if (TMath::Abs(nsigma) > fRange) {\r
447 pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
448 }\r
449 else{\r
450 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
451 }\r
452 }\r
453 //NON electrons\r
454 else{\r
455 // E/p < 0.5 --> return probability below E/p = 0.5\r
456 if ( nsigma == -999){\r
457 pEMCAL[species] = GetLowProb(pt,type,charge);\r
458 }\r
459 // E/p > 1.5 --> return probability above E/p = 1.5\r
460 else if ( nsigma == 999){\r
461 pEMCAL[species] = GetHighProb(pt,type,charge);\r
462 }\r
463 // in parametrized region --> calculate probability for corresponding Gauss curve\r
464 else{\r
465 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
466 \r
467 // normalize to total probability == 1\r
468 pEMCAL[species]*=GetExpectedNorm(pt,type,charge);\r
469 }\r
470 }\r
471 }\r
472\r
473 // return the electron probability\r
474 proba = pEMCAL[AliPID::kElectron]; \r
475\r
476 }\r
477\r
478 return proba;\r
479\r
480}\r