- Reshuffling of the particle codes in AliPID. Now the light nuclei are between the
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEMCALPIDResponse.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 //                                                                      //
17 // AliEMCALPIDResponse                                                  //
18 //                                                                      //
19 // EMCAL class to perfom PID                                            //
20 // This is a prototype and still under development                      //
21 //                                                                      //
22 // ---------------------------------------------------------------------//
23 // GetNumberOfSigmas():                                                 //
24 //                                                                      //
25 // Electrons:  Number of Sigmas for E/p value                           //
26 //             Parametrization of LHC11a (after recalibration)          //
27 //                                                                      //
28 // NON electrons:                                                       //
29 //             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
30 //             --> return +/- 99                                        //
31 //             Otherwise                                                //
32 //             --> return nsigma (parametrization of LHC10e)            //
33 //                                                                      //
34 // NO Parametrization (outside pT range): --> return -999               //
35 //                                                                      //
36 // ---------------------------------------------------------------------//
37 // ComputeEMCALProbability():                                           //
38 //                                                                      //
39 // Electrons:  Probability from Gaussian distribution                   //
40 //                                                                      //
41 // NON electrons:                                                       //
42 //             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
43 //             --> probability to find particles below or above thr.    //
44 //             Otherwise                                                //
45 //             -->  Probability from Gaussian distribution              //
46 //                  (proper normalization to each other?)               //
47 //                                                                      //
48 // NO Parametrization (outside pT range): --> return kFALSE             //
49 //////////////////////////////////////////////////////////////////////////
50
51 #include <TF1.h>
52 #include <TMath.h>
53
54 #include "AliEMCALPIDResponse.h"       //class header
55
56 #include "AliLog.h"   
57
58 ClassImp(AliEMCALPIDResponse)
59
60 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61 AliEMCALPIDResponse::AliEMCALPIDResponse():
62   TObject(),
63   fNorm(NULL),
64   fkPIDParams(NULL)
65 {
66   //
67   //  The default constructor
68   //
69
70
71   fNorm = new TF1("fNorm","gaus",-20,20); 
72 }
73 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
74 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
75   TObject(other),
76   fNorm(other.fNorm),
77   fkPIDParams(other.fkPIDParams)
78 {
79   //
80   //  The copy constructor
81   //
82
83 }
84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
86 {
87   //
88   //  The assignment operator
89   //
90
91   if(this == &other) return *this;
92   
93   // Make copy
94   TObject::operator=(other);
95   fNorm = other.fNorm;
96   fkPIDParams = other.fkPIDParams;
97
98  
99   return *this;
100 }
101 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102 AliEMCALPIDResponse::~AliEMCALPIDResponse() {
103
104   delete fNorm;
105
106 }
107 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
108 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {
109   //
110   // Calculates the expected sigma of the PID signal as the function of 
111   // the information stored in the track, for the specified particle type 
112   //  
113   //
114   
115   Double_t norm = 1.;
116
117   // Check the charge
118   if( charge != -1 && charge != 1){
119     return norm;
120   }
121
122   // Get the parameters for this particle type and pt
123   const TVectorD *params = GetParams(n, pt, charge);
124
125   // IF not in momentum range, NULL is returned --> return default value
126   if(!params) return norm;
127
128   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
129   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
130   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
131   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
132   Double_t probLow  = (*params)[6];   // probability to be below eopMin
133   Double_t probHigh = (*params)[7];   // probability to be above eopMax
134
135   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
136   fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
137   norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
138
139   return norm;
140 }
141 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
142 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {
143       
144   Double_t nsigma = -999.;
145
146   // Check the charge
147   if( charge != -1 && charge != 1){
148     return nsigma;
149   }
150
151   // Get the parameters for this particle type and pt
152   const TVectorD *params = GetParams(n, pt, charge);
153
154   // IF not in momentum range, NULL is returned --> return default value
155   if(!params) return nsigma;
156
157   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
158   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
159   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
160   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
161
162   // if electron
163   if(n == AliPID::kElectron){
164     if(sigma != 0) nsigma = (eop - mean) / sigma;
165   }
166
167   // if NON electron
168   else{
169     if ( eop < eopMin )
170       nsigma = -99;    // not parametrized (smaller than eopMin)
171     else if ( eop > eopMax )
172       nsigma = 99.;     // not parametrized (bigger than eopMax)
173     else{
174       if(sigma != 0) nsigma = (eop - mean) / sigma; 
175     }
176   }
177
178   return nsigma;
179
180 }
181 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
182 Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
183   //
184   //
185   Double_t fRange  = 5.0;   // hardcoded (???)
186   Double_t nsigma  = 0.0;
187   
188
189   // Check the charge
190   if( charge != -1 && charge != 1){
191     return kFALSE;
192   }
193
194  
195   // default value (will be returned, if pt below threshold)
196   for (Int_t species = 0; species < nSpecies; species++) {
197     pEMCAL[species] = 1./nSpecies;
198   }
199
200   // set E/p range
201   if(eop < 0.05) eop = 0.05;
202   if(eop > 2.00) eop = 2.00;
203   
204   for (Int_t species = 0; species < nSpecies; species++) {
205     
206     AliPID::EParticleType type = AliPID::EParticleType(species);
207
208     // Get the parameters for this particle type and pt
209     const TVectorD *params = GetParams(species, pt, charge);
210     
211     // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
212     if(!params) return kFALSE;
213
214     Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
215     Double_t probLow  = (*params)[6];   // probability to be below eopMin
216     Double_t probHigh = (*params)[7];   // probability to be above eopMax
217
218     // get nsigma value for each particle type at this E/p value
219     nsigma = GetNumberOfSigmas(pt,eop,type,charge);
220
221     // electrons (standard Gaussian calculation of probabilities)
222     if(type == AliPID::kElectron){
223       if (TMath::Abs(nsigma) > fRange) {
224         pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
225       }
226       else{
227         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
228       }
229     }
230     //NON electrons
231     else{
232       // E/p < eopMin  -->  return probability below E/p = eopMin
233       if ( nsigma == -99){
234         pEMCAL[species] = probLow;
235       }
236       // E/p > eopMax  -->  return probability above E/p = eopMax
237       else if ( nsigma == 99){
238         pEMCAL[species] = probHigh;
239       }
240       // in parametrized region --> calculate probability for corresponding Gauss curve
241       else{
242         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
243         
244         // normalize to total probability == 1
245         pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
246       }
247     }
248   }
249
250   return kTRUE;
251
252 }
253 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
254 const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
255   //
256   // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
257   //
258   // 0 = momMin
259   // 1 = momMax
260   // 2 = mean of Gaus
261   // 3 = sigma of Gaus
262   // 4 = eopLow   
263   // 5 = eopHig   
264   // 6 = probLow  (not used for electrons)
265   // 7 = probHigh (not used for electrons)
266   //
267
268   if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
269   if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
270
271   TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
272   if(!particlePar) return NULL;
273
274   TIter parIter(particlePar);
275   const TVectorD *parameters = NULL;
276   Double_t momMin = 0.;
277   Double_t momMax = 0.;
278
279   while((parameters = static_cast<const TVectorD *>(parIter()))){
280
281     momMin = (*parameters)[0];
282     momMax = (*parameters)[1];
283
284     if( fPt > momMin && fPt < momMax ) return parameters;
285
286   }  
287   AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
288
289   return parameters;
290 }