Removing useless include
[u/mrichter/AliRoot.git] / STEER / 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 AliEMCALPIDResponse::AliEMCALPIDResponse():\r
58 fNorm(0)\r
59 {\r
60   //\r
61   //  The default constructor\r
62   //\r
63 \r
64   for(Int_t i = 0; i < fNptBins; i++){\r
65 \r
66     fPtCutMin[i] = 0.0;\r
67 \r
68     for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++){\r
69 \r
70       fMeanEoP[j][i]  = 0.0;\r
71       fSigmaEoP[j][i] = 0.0;\r
72       fProbLow[j][i]  = 0.0;\r
73       fProbHigh[j][i] = 0.0;\r
74 \r
75     }\r
76   } \r
77   fPtCutMin[fNptBins] = 0.0;\r
78 \r
79   fNorm = new TF1("fNorm","gaus",-20,20); \r
80 \r
81   SetPtBoundary();\r
82   SetParametrizations();\r
83 }\r
84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
85 void AliEMCALPIDResponse::SetPtBoundary(){\r
86   //\r
87   // Set boundaries for momentum bins\r
88   //\r
89   fPtCutMin[0] = 1.5;\r
90   fPtCutMin[1] = 2.5;\r
91   fPtCutMin[2] = 3.5;\r
92   fPtCutMin[3] = 4.5;\r
93   fPtCutMin[4] = 5.5;\r
94   fPtCutMin[5] = 6.5;\r
95 \r
96 }\r
97 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
98 void AliEMCALPIDResponse::SetParametrizations(){\r
99 \r
100 // This are the preliminary parametrizations (hard coded)\r
101 // For electrons from LHC11a (Deepa Thomas)\r
102 // For NON-electrons from LHC10e (TOF/TPC analysis)\r
103   \r
104   // Gaussian mean\r
105   Float_t mean[4][6]    = {\r
106     { 0.932, 0.997, 0.998, 1.001, 1.011, 1.011 },                              // electrons\r
107     { 0.227804, 0.34839, 0.404077, -0.107795, -4.14584, 0.5 },             // NON electrons\r
108     { -2.10377, 0.0582898, 0.0582898, 0.0582898, 0.0582898, 0.0582898  },   // protons\r
109     { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}                                              // anti-protons\r
110   };        \r
111   \r
112   // Gaussian sigma\r
113   Float_t sigma[4][6]= {\r
114     { 0.0866, 0.0693, 0.0664, 0.0583, 0.0488, 0.0515  },                      // electrons\r
115     { 0.310831, 0.267586, 0.404077, 0.381968, 1.46183, 0.314687 },          // NON electrons\r
116     { 0.603209, 0.255332, 0.255332, 0.255332, 0.255332, 0.255332},          // protons\r
117     { 0.516837, 0.351516,0.351516,0.351516,0.351516,0.351516 }              // anti - protons\r
118   };\r
119 \r
120   // lower probability\r
121   Float_t probL[3][6] = {\r
122     { 0.928689, 0.938455, 0.940448, 0.948496, 0.955981, 0.951923 },         // NON electrons\r
123     { 0.974518, 0.978088, 0.975089, 0.975089, 0.975089,0.975089},           // protons\r
124     { 0.824037, 0.861149, 0.898734, 0.898734, 0.898734, 0.898734},          // anti - protons\r
125   };\r
126 \r
127   // upper probability\r
128   Float_t probH[3][6] = {\r
129     { 0.00030227, 4.04106e-05, 0.000147406, 0., 0.000956938, 0.00106838 },  // NON electrons\r
130     { 0.000157945, 0., 0., 0., 0., 0. },                                    // protons\r
131     { 0.00343237, 0., 0., 0., 0., 0.}                                       // anti - protons\r
132   };\r
133   \r
134   \r
135   // set parametrizations\r
136   Int_t spec = 0;\r
137   for (Int_t species = 0; species < 2*AliPID::kSPECIES; species++) {          // first negative particles and then positive\r
138     for (Int_t pt = 0; pt < fNptBins; pt++){\r
139 \r
140       switch(species){\r
141       case 0:      // electrons\r
142         spec = 0;\r
143         break;\r
144       case 4:      // anti - protons\r
145         spec = 3;\r
146         break;\r
147       case 5:      // positrons\r
148         spec = 0;\r
149         break;\r
150       case 9:      // protons\r
151         spec = 2;\r
152         break;\r
153       default:     // NON electrons\r
154         spec = 1;\r
155         break;\r
156       }\r
157     \r
158 \r
159       fMeanEoP[species][pt]  = mean[spec][pt];    \r
160       fSigmaEoP[species][pt] = sigma[spec][pt];  \r
161       if( spec == 0) {     // electrons have NO lower and upper probability thresholds --> set to 0\r
162         fProbLow[species][pt]  = 0.;\r
163         fProbHigh[species][pt] = 0.;    \r
164       }\r
165       else{\r
166         fProbLow[species][pt]  = probL[spec-1][pt];\r
167         fProbHigh[species][pt] = probH[spec-1][pt];     \r
168       } \r
169     \r
170     }//loop pt bins\r
171   }//loop species\r
172 }\r
173 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
174 Int_t AliEMCALPIDResponse::GetPtBin(Float_t pt) const {\r
175   //\r
176   // Returns the momentum bin index\r
177   //\r
178 \r
179   Int_t i = -1;\r
180   while(pt > fPtCutMin[i+1] && i+1 < fNptBins) i++;\r
181 \r
182   return i;\r
183 }\r
184 \r
185 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
186 Double_t AliEMCALPIDResponse::GetExpectedSignal( Float_t pt, AliPID::EParticleType n, Int_t charge) const  {\r
187   //\r
188   // Calculates the expected PID signal as the function of \r
189   // the information stored in the track, for the specified particle type \r
190   //  \r
191 \r
192   Double_t signal = 0.;\r
193 \r
194   // Check the charge\r
195   if( charge != -1 && charge != 1){\r
196     \r
197     return signal;\r
198   }\r
199   \r
200   // Get the pt bin\r
201   Int_t ptBin = GetPtBin(pt);\r
202   \r
203   // Get the species (first negative , then positive)\r
204   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
205 \r
206   // Get the signal\r
207   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
208     signal = fMeanEoP[species][ptBin];\r
209   }\r
210 \r
211   return signal;\r
212 \r
213 }\r
214 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
215 Double_t AliEMCALPIDResponse::GetExpectedSigma( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
216   //\r
217   // Calculates the expected sigma of the PID signal as the function of \r
218   // the information stored in the track, for the specified particle type \r
219   //  \r
220   //\r
221   \r
222   Double_t sigma = 999.;\r
223 \r
224   // Check the charge\r
225   if( charge != -1 && charge != 1){\r
226     \r
227     return sigma;\r
228   }  \r
229 \r
230   // Get the pt bin\r
231   Int_t ptBin = GetPtBin(pt);\r
232   \r
233   // Get the species (first negative , then positive)\r
234   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
235 \r
236   // Get the sigma\r
237   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
238     sigma = fSigmaEoP[species][ptBin];\r
239   }\r
240 \r
241   return sigma;\r
242  \r
243 }\r
244 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
245 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
246   //\r
247   // Calculates the expected sigma of the PID signal as the function of \r
248   // the information stored in the track, for the specified particle type \r
249   //  \r
250   //\r
251   \r
252   Double_t norm = 1.;\r
253 \r
254   // Check the charge\r
255   if( charge != -1 && charge != 1){\r
256     \r
257     return norm;\r
258   }\r
259 \r
260   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )\r
261   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
262   norm = 1./fNorm->Integral(fLowEoP,fHighEoP)*(1-GetLowProb(pt,n,charge)-GetHighProb(pt,n,charge));\r
263 \r
264   return norm;\r
265 }\r
266 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
267 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {\r
268       \r
269   Double_t mean  = GetExpectedSignal(pt,n,charge);\r
270   Double_t sigma = GetExpectedSigma(pt,n,charge);\r
271 \r
272   // if electron\r
273   if(n == AliPID::kElectron){\r
274     return (eop - mean) / sigma;\r
275   }\r
276 \r
277   // if NON electron\r
278   else{\r
279     if ( eop < fLowEoP )\r
280       return -999.;    // not parametrized \r
281     else if ( eop > fHighEoP )\r
282       return 999.;     // not parametrized \r
283     else{\r
284       return (eop - mean) / sigma; \r
285     }\r
286   }\r
287 }\r
288 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
289 Double_t AliEMCALPIDResponse::GetLowProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
290   //\r
291   //\r
292   \r
293   Double_t prob = 0.;\r
294 \r
295   // Check the charge\r
296   if( charge != -1 && charge != 1){\r
297     \r
298     return prob;\r
299   }\r
300   \r
301   // Get the pt bin\r
302   Int_t ptBin = GetPtBin(pt);\r
303   \r
304   // Get the species (first negative , then positive)\r
305   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
306 \r
307   // Get the probability\r
308   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
309     prob = fProbLow[species][ptBin];\r
310   }\r
311 \r
312   return prob;\r
313  \r
314 }\r
315 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
316 Double_t AliEMCALPIDResponse::GetHighProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
317   //\r
318   //\r
319   \r
320   Double_t prob = 0.;\r
321 \r
322   // Check the charge\r
323   if( charge != -1 && charge != 1){\r
324     \r
325     return prob;\r
326   }\r
327   \r
328   // Get the pt bin\r
329   Int_t ptBin = GetPtBin(pt);\r
330   \r
331   // Get the species (first negative , then positive)\r
332   Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
333 \r
334   // Get the probability\r
335   if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
336     prob = fProbHigh[species][ptBin];\r
337   }\r
338 \r
339   return prob;\r
340   \r
341 }\r
342 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
343 Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {\r
344   //\r
345   //\r
346   Double_t fRange  = 5.0;   // hardcoded \r
347   Double_t nsigma  = 0.0;\r
348   Double_t sigma   = 0.0;\r
349   Double_t proba   = -1.;\r
350 \r
351 \r
352   // Check the charge\r
353   if( charge != -1 && charge != 1){\r
354     \r
355     return proba;\r
356   }\r
357 \r
358  \r
359   // default value (will be returned, if pt below threshold)\r
360   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
361     pEMCAL[species] = -1;\r
362   }\r
363 \r
364   if( GetPtBin(pt) > -1 ){\r
365 \r
366     // set E/p range\r
367     if(eop < 0.05) eop = 0.05;\r
368     if(eop > 2.00) eop = 2.00;\r
369 \r
370     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
371 \r
372       AliPID::EParticleType type = AliPID::EParticleType(species);\r
373 \r
374       // get nsigma value for each particle type at this E/p value\r
375       nsigma = GetNumberOfSigmas(pt,eop,type,charge);\r
376       sigma  = GetExpectedSigma(pt,type,charge);\r
377 \r
378       // electrons (standard Gaussian calculation of probabilities)\r
379       if(type == AliPID::kElectron){\r
380         if (TMath::Abs(nsigma) > fRange) {\r
381           pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
382         }\r
383         else{\r
384           pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
385         }\r
386       }\r
387       //NON electrons\r
388       else{\r
389         // E/p < 0.5  -->  return probability below E/p = 0.5\r
390         if ( nsigma == -999){\r
391           pEMCAL[species] = GetLowProb(pt,type,charge);\r
392         }\r
393         // E/p > 1.5  -->  return probability above E/p = 1.5\r
394         else if ( nsigma == 999){\r
395           pEMCAL[species] = GetHighProb(pt,type,charge);\r
396         }\r
397         // in parametrized region --> calculate probability for corresponding Gauss curve\r
398         else{\r
399           pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
400         \r
401           // normalize to total probability == 1\r
402           pEMCAL[species]*=GetExpectedNorm(pt,type,charge);\r
403         }\r
404       }\r
405     }\r
406 \r
407     // return the electron probability\r
408     proba = pEMCAL[AliPID::kElectron];  \r
409 \r
410   }\r
411 \r
412   return proba;\r
413 \r
414 }\r