Update master to aliroot
[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   fCurrCentrality(-1.),
65   fkPIDParams(NULL)
66 {
67   //
68   //  The default constructor
69   //
70
71
72   fNorm = new TF1("fNorm","gaus",-20,20); 
73 }
74 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
75 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
76   TObject(other),
77   fNorm(other.fNorm),
78   fCurrCentrality(other.fCurrCentrality),
79   fkPIDParams(other.fkPIDParams)
80 {
81   //
82   //  The copy constructor
83   //
84
85 }
86 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
88 {
89   //
90   //  The assignment operator
91   //
92
93   if(this == &other) return *this;
94   
95   // Make copy
96   TObject::operator=(other);
97   fNorm = other.fNorm;
98   fCurrCentrality = other.fCurrCentrality;
99   fkPIDParams = other.fkPIDParams;
100
101  
102   return *this;
103 }
104 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
105 AliEMCALPIDResponse::~AliEMCALPIDResponse() {
106
107   delete fNorm;
108
109 }
110 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {
112   //
113   // Calculates the expected sigma of the PID signal as the function of 
114   // the information stored in the track, for the specified particle type 
115   //  
116   //
117   
118   Double_t norm = 1.;
119
120   // Check the charge
121   if( charge != -1 && charge != 1){
122     return norm;
123   }
124
125   // Get the parameters for this particle type and pt
126   const TVectorD *params = GetParams(n, pt, charge);
127
128   // IF not in momentum range, NULL is returned --> return default value
129   if(!params) return norm;
130
131   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
132   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
133   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
134   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
135   Double_t probLow  = (*params)[6];   // probability to be below eopMin
136   Double_t probHigh = (*params)[7];   // probability to be above eopMax
137
138   // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
139   fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
140   norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
141
142   return norm;
143 }
144 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
145 Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {
146       
147   Double_t nsigma = -999.;
148
149   // Check the charge
150   if( charge != -1 && charge != 1){
151     return nsigma;
152   }
153
154   // Get the parameters for this particle type and pt
155   const TVectorD *params = GetParams(n, pt, charge);
156
157   // IF not in momentum range, NULL is returned --> return default value
158   if(!params) return nsigma;
159
160   Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
161   Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
162   Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
163   Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
164
165   // if electron
166   if(n == AliPID::kElectron){
167     if(sigma != 0) nsigma = (eop - mean) / sigma;
168   }
169
170   // if NON electron
171   else{
172     if ( eop < eopMin )
173       nsigma = -99;    // not parametrized (smaller than eopMin)
174     else if ( eop > eopMax )
175       nsigma = 99.;     // not parametrized (bigger than eopMax)
176     else{
177       if(sigma != 0) nsigma = (eop - mean) / sigma; 
178     }
179   }
180
181   return nsigma;
182
183 }
184 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
185 Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
186   //
187   //
188   Double_t fRange  = 5.0;   // hardcoded (???)
189   Double_t nsigma  = 0.0;
190   
191
192   // Check the charge
193   if( charge != -1 && charge != 1){
194     return kFALSE;
195   }
196
197  
198   // default value (will be returned, if pt below threshold)
199   for (Int_t species = 0; species < nSpecies; species++) {
200     pEMCAL[species] = 1./nSpecies;
201   }
202
203   // set E/p range
204   if(eop < 0.05) eop = 0.05;
205   if(eop > 2.00) eop = 2.00;
206   
207   for (Int_t species = 0; species < nSpecies; species++) {
208     
209     AliPID::EParticleType type = AliPID::EParticleType(species);
210
211     // Get the parameters for this particle type and pt
212     const TVectorD *params = GetParams(species, pt, charge);
213     
214     // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
215     if(!params) return kFALSE;
216
217     Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
218     Double_t probLow  = (*params)[6];   // probability to be below eopMin
219     Double_t probHigh = (*params)[7];   // probability to be above eopMax
220
221     // get nsigma value for each particle type at this E/p value
222     nsigma = GetNumberOfSigmas(pt,eop,type,charge);
223
224     // electrons (standard Gaussian calculation of probabilities)
225     if(type == AliPID::kElectron){
226       if (TMath::Abs(nsigma) > fRange) {
227         pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
228       }
229       else{
230         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
231       }
232     }
233     //NON electrons
234     else{
235       // E/p < eopMin  -->  return probability below E/p = eopMin
236       if ( nsigma == -99){
237         pEMCAL[species] = probLow;
238       }
239       // E/p > eopMax  -->  return probability above E/p = eopMax
240       else if ( nsigma == 99){
241         pEMCAL[species] = probHigh;
242       }
243       // in parametrized region --> calculate probability for corresponding Gauss curve
244       else{
245         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
246         
247         // normalize to total probability == 1
248         pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
249       }
250     }
251   }
252
253   return kTRUE;
254
255 }
256 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
257 const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
258   //
259   // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
260   //
261   // 0 = momMin
262   // 1 = momMax
263   // 2 = mean of Gaus
264   // 3 = sigma of Gaus
265   // 4 = eopLow   
266   // 5 = eopHig   
267   // 6 = probLow  (not used for electrons)
268   // 7 = probHigh (not used for electrons)
269   //
270   // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
271   // so first the correct centrality bin has to be found
272
273   // **** Centrality bins (hard coded for the moment)
274   const Int_t nCent = 7;
275   Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
276
277   if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
278   if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
279
280   TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
281   if(!particlePar) return NULL;
282
283   const TVectorD *parameters = NULL;
284   Double_t momMin = 0.;
285   Double_t momMax = 0.;
286
287   // is the centrality dependent parametrization used
288   TString arrayName = particlePar->GetName();
289
290   // centrality dependent parametrization
291   if(arrayName.Contains("Centrality")){
292     
293     for(Int_t iCent = 0; iCent < nCent; iCent++){
294
295       if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
296         
297         TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
298         if(!centPar) return NULL;
299         
300         TIter centIter(centPar);
301         parameters = NULL;
302         momMin = 0.;
303         momMax = 0.;
304         
305         while((parameters = static_cast<const TVectorD *>(centIter()))){
306           
307           momMin = (*parameters)[0];
308           momMax = (*parameters)[1];
309           
310           if( fPt > momMin && fPt < momMax ) return parameters;
311          
312         } 
313       }
314     }
315   }
316
317   // NO centrality dependent parametrization
318   else{
319
320     TIter parIter(particlePar);
321     while((parameters = static_cast<const TVectorD *>(parIter()))){
322       
323       momMin = (*parameters)[0];
324       momMax = (*parameters)[1];
325       
326       if( fPt > momMin && fPt < momMax ) return parameters;
327       
328     }  
329   }
330   AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
331
332   return parameters;
333 }