]>
Commit | Line | Data |
---|---|---|
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 | ||
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); | |
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); | |
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 | Double_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 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
259 | const 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 | } |