c11cb8db4258596959abf68b65d05f638eefa1c2
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEMCALPIDResponse.cxx
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
48 #include <TF1.h>\r
49 #include <TMath.h>\r
50 \r
51 #include "AliEMCALPIDResponse.h"       //class header\r
52 \r
53 #include "AliLog.h"   \r
54 \r
55 ClassImp(AliEMCALPIDResponse)\r
56 \r
57 const Float_t AliEMCALPIDResponse::fLowEoP  = 0.5;   // lower E/p threshold for NON electrons\r
58 const Float_t AliEMCALPIDResponse::fHighEoP = 1.5;   // upper E/p threshold for NON electrons\r
59 \r
60 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
61 AliEMCALPIDResponse::AliEMCALPIDResponse():\r
62   TObject(),\r
63   fNorm(NULL)\r
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
89 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
90 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):\r
91   TObject(other),\r
92   fNorm(other.fNorm)\r
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
113 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
114 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)\r
115 {\r
116   //\r
117   //  The assignment operator\r
118   //\r
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
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
137   \r
138   fPtCutMin[fNptBins] = 0.0;\r
139   SetPtBoundary();\r
140   SetParametrizations();\r
141 \r
142   return *this;\r
143 }\r
144 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
145 AliEMCALPIDResponse::~AliEMCALPIDResponse() {\r
146 \r
147   delete fNorm;\r
148 \r
149 }\r
150 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
151 void 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
164 void 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
240 Int_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
252 Double_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
281 Double_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
311 Double_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
333 Double_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
355 Double_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
382 Double_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
409 Double_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
415   Double_t proba   = 999.;\r
416   \r
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
427     pEMCAL[species] = 999.;\r
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