EMCAL Response in OADB
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEMCALPIDResponse.cxx
CommitLineData
b2138b40 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 999 //
49//////////////////////////////////////////////////////////////////////////
50
51#include <TF1.h>
52#include <TMath.h>
53
54#include "AliEMCALPIDResponse.h" //class header
55
56#include "AliLog.h"
57
58ClassImp(AliEMCALPIDResponse)
59
60//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61AliEMCALPIDResponse::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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
74AliEMCALPIDResponse::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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85AliEMCALPIDResponse & 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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102AliEMCALPIDResponse::~AliEMCALPIDResponse() {
103
104 delete fNorm;
105
106}
107//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
108Double_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);
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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
142Double_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);
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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
182Double_t AliEMCALPIDResponse::ComputeEMCALProbability(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 Double_t proba = 999.;
188
189
190 // Check the charge
191 if( charge != -1 && charge != 1){
192 return proba;
193 }
194
195
196 // default value (will be returned, if pt below threshold)
197 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
198 pEMCAL[species] = 999.;
199 }
200
201 // set E/p range
202 if(eop < 0.05) eop = 0.05;
203 if(eop > 2.00) eop = 2.00;
204
205 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
206
207 AliPID::EParticleType type = AliPID::EParticleType(species);
208
209 // Get the parameters for this particle type and pt
210 const TVectorD *params = GetParams(species, pt);
211
212 // IF not in momentum range, NULL is returned --> return default value
213 if(!params) return proba;
214
215 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
216 Double_t probLow = (*params)[6]; // probability to be below eopMin
217 Double_t probHigh = (*params)[7]; // probability to be above eopMax
218
219 // get nsigma value for each particle type at this E/p value
220 nsigma = GetNumberOfSigmas(pt,eop,type,charge);
221
222 // electrons (standard Gaussian calculation of probabilities)
223 if(type == AliPID::kElectron){
224 if (TMath::Abs(nsigma) > fRange) {
225 pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
226 }
227 else{
228 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
229 }
230 }
231 //NON electrons
232 else{
233 // E/p < eopMin --> return probability below E/p = eopMin
234 if ( nsigma == -99){
235 pEMCAL[species] = probLow;
236 }
237 // E/p > eopMax --> return probability above E/p = eopMax
238 else if ( nsigma == 99){
239 pEMCAL[species] = probHigh;
240 }
241 // in parametrized region --> calculate probability for corresponding Gauss curve
242 else{
243 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
244
245 // normalize to total probability == 1
246 pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
247 }
248 }
249
250 // return the electron probability
251 proba = pEMCAL[AliPID::kElectron];
252
253 }
254
255 return proba;
256
257}
258//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
259const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt) const {
260 //
261 // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
262 //
263 // 0 = momMin
264 // 1 = momMax
265 // 2 = mean of Gaus
266 // 3 = sigma of Gaus
267 // 4 = eopLow
268 // 5 = eopHig
269 // 6 = probLow (not used for electrons)
270 // 7 = probHigh (not used for electrons)
271 //
272
273 if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
274
275 TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
276 if(!particlePar) return NULL;
277
278 TIter parIter(particlePar);
279 const TVectorD *parameters = NULL;
280 Double_t momMin = 0.;
281 Double_t momMax = 0.;
282
283 while((parameters = static_cast<const TVectorD *>(parIter()))){
284
285 momMin = (*parameters)[0];
286 momMax = (*parameters)[1];
287
288 if( fPt > momMin && fPt < momMax ) return parameters;
289
290 }
291 AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
292
293 return parameters;
294}