]>
Commit | Line | Data |
---|---|---|
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 | // // | |
00a38d07 | 48 | // NO Parametrization (outside pT range): --> return kFALSE // |
b2138b40 | 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), | |
87da0205 | 64 | fCurrCentrality(-1.), |
b2138b40 | 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), | |
87da0205 | 78 | fCurrCentrality(other.fCurrCentrality), |
b2138b40 | 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; | |
87da0205 | 98 | fCurrCentrality = other.fCurrCentrality; |
b2138b40 | 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 | |
1f631618 | 126 | const TVectorD *params = GetParams(n, pt, charge); |
b2138b40 | 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 | |
1f631618 | 155 | const TVectorD *params = GetParams(n, pt, charge); |
b2138b40 | 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 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
00a38d07 | 185 | Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const { |
b2138b40 | 186 | // |
187 | // | |
188 | Double_t fRange = 5.0; // hardcoded (???) | |
189 | Double_t nsigma = 0.0; | |
b2138b40 | 190 | |
191 | ||
192 | // Check the charge | |
193 | if( charge != -1 && charge != 1){ | |
00a38d07 | 194 | return kFALSE; |
b2138b40 | 195 | } |
196 | ||
197 | ||
198 | // default value (will be returned, if pt below threshold) | |
00a38d07 | 199 | for (Int_t species = 0; species < nSpecies; species++) { |
200 | pEMCAL[species] = 1./nSpecies; | |
b2138b40 | 201 | } |
202 | ||
203 | // set E/p range | |
204 | if(eop < 0.05) eop = 0.05; | |
205 | if(eop > 2.00) eop = 2.00; | |
206 | ||
00a38d07 | 207 | for (Int_t species = 0; species < nSpecies; species++) { |
b2138b40 | 208 | |
209 | AliPID::EParticleType type = AliPID::EParticleType(species); | |
210 | ||
211 | // Get the parameters for this particle type and pt | |
1f631618 | 212 | const TVectorD *params = GetParams(species, pt, charge); |
b2138b40 | 213 | |
00a38d07 | 214 | // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE |
215 | if(!params) return kFALSE; | |
b2138b40 | 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 | } | |
b2138b40 | 251 | } |
252 | ||
00a38d07 | 253 | return kTRUE; |
b2138b40 | 254 | |
255 | } | |
256 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
1f631618 | 257 | const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const { |
b2138b40 | 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 | // | |
87da0205 | 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}; | |
b2138b40 | 276 | |
277 | if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL; | |
1f631618 | 278 | if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons |
b2138b40 | 279 | |
280 | TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle)); | |
281 | if(!particlePar) return NULL; | |
1f631618 | 282 | |
b2138b40 | 283 | const TVectorD *parameters = NULL; |
284 | Double_t momMin = 0.; | |
285 | Double_t momMax = 0.; | |
286 | ||
87da0205 | 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++){ | |
b2138b40 | 294 | |
87da0205 | 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 | } | |
b2138b40 | 316 | |
87da0205 | 317 | // NO centrality dependent parametrization |
318 | else{ | |
b2138b40 | 319 | |
87da0205 | 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 | } | |
b2138b40 | 330 | AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt)); |
331 | ||
332 | return parameters; | |
333 | } |