]>
Commit | Line | Data |
---|---|---|
dc293ae9 | 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 | /* $Id$ */ | |
17 | /* History of cvs commits: | |
18 | * | |
19 | * $Log$ | |
4ce73a00 | 20 | * Revision 1.10 2007/03/09 14:34:11 gustavo |
21 | * Correct probability calculation, added missing initialization of data members | |
22 | * | |
4cfd8ef9 | 23 | * Revision 1.9 2007/02/20 20:17:43 hristov |
24 | * Corrected array size, removed warnings (icc) | |
25 | * | |
d69ab345 | 26 | * Revision 1.8 2006/12/19 08:49:35 gustavo |
27 | * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD() | |
28 | * | |
dc293ae9 | 29 | * |
30 | */ | |
31 | // to compute PID for all the clusters in ESDs.root file | |
32 | // the ESDs.root have to be in the same directory as the class | |
33 | // | |
34 | // | |
35 | // | |
36 | // | |
37 | // | |
38 | // AliEMCALPID::CalculPID(Energy,Lambda0) | |
39 | // Calcul PID for all clusters in AliESDs.root file | |
40 | // keep this function for the moment for a simple verification, could be removed | |
41 | // | |
42 | // | |
43 | // | |
44 | // AliEMCALPID::CalculPID(Energy,Lambda0) | |
45 | // calcul PID Weght for a cluster with Energy, Lambda0 . | |
46 | // Double_t PIDFinal[AliPID::kSPECIESN] is the standard PID for : | |
47 | // | |
48 | // | |
49 | // | |
50 | // kElectron : fPIDFinal[0] | |
51 | // kMuon : fPIDFinal[1] | |
52 | // kPion : fPIDFinal[2] | |
53 | // kKaon : fPIDFinal[3] | |
54 | // kProton : fPIDFinal[4] | |
55 | // kPhoton : fPIDFinal[5] | |
56 | // kPi0 : fPIDFinal[6] | |
57 | // kNeutron : fPIDFinal[7] | |
58 | // kKaon0 : fPIDFinal[8] | |
59 | // kEleCon : fPIDFinal[9] | |
60 | // kUnknown : fPIDFinal[10] | |
61 | // | |
62 | // | |
63 | // PID[3] is a simple PID for | |
64 | // Electron & Photon PID[0] | |
65 | // Pi0 PID[1] | |
66 | // Hadron PID[2] | |
67 | // | |
68 | // | |
69 | // | |
70 | // | |
71 | // | |
72 | // --- ROOT system --- | |
73 | ||
74 | // standard C++ includes | |
75 | #include <Riostream.h> | |
76 | // #include <cstdlib> | |
77 | // #include <cmath> | |
78 | ||
79 | // ROOT includes | |
80 | #include "TTree.h" | |
81 | #include "TStyle.h" | |
82 | #include "TVector3.h" | |
83 | #include "TBranch.h" | |
84 | #include "TClonesArray.h" | |
85 | #include "TCanvas.h" | |
86 | #include "TLorentzVector.h" | |
87 | #include "TMath.h" | |
88 | #include "TFile.h" | |
89 | #include "TH1.h" | |
90 | #include "TH2.h" | |
91 | #include "TParticle.h" | |
92 | ||
93 | // // STEER includes | |
94 | // #include "AliRun.h" | |
95 | // #include "AliRunLoader.h" | |
96 | // #include "AliHeader.h" | |
97 | // #include "AliLoader.h" | |
98 | // #include "AliStack.h" | |
99 | // #include "AliESDtrack.h" | |
100 | // #include "AliESD.h" | |
101 | #include "AliLog.h" | |
102 | #include "AliEMCALPID.h" | |
103 | ||
104 | ClassImp(AliEMCALPID) | |
105 | ||
4cfd8ef9 | 106 | //______________________________________________ |
107 | AliEMCALPID::AliEMCALPID(): | |
108 | fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE) | |
dc293ae9 | 109 | { |
4cfd8ef9 | 110 | // |
111 | // Constructor. | |
112 | // Initialize all constant values which have to be used | |
113 | // during PID algorithm execution | |
114 | // | |
115 | ||
116 | // set flag for printing to FALSE by default | |
117 | fPrintInfo = kFALSE; | |
118 | ||
119 | // as a first step, all array elements are initialized to 0.0 | |
120 | Int_t i, j; | |
121 | for (i = 0; i < 6; i++) { | |
122 | for (j = 0; j < 6; j++) { | |
123 | fGamma[i][j] = fHadron[i][j] = fPiZero5to10[i][j] = fPiZero10to60[i][j] = 0.; | |
124 | } | |
125 | } | |
126 | ||
127 | // then, only the ones which must be not zero are initialized | |
128 | // while the others will remain to the value 0.0 | |
129 | ||
130 | fGamma[0][0] = 0.038022; | |
131 | fGamma[0][1] = -0.0001883; | |
132 | fGamma[0][2] = 5.449e-06; | |
133 | ||
134 | fGamma[1][0] = 0.207313; | |
135 | fGamma[1][1] = -0.000978; | |
136 | fGamma[1][2] = 0.00001634; | |
137 | ||
138 | fGamma[2][0] = 0.043364; | |
139 | fGamma[2][1] = -0.0002048; | |
140 | fGamma[2][2] = 8.661e-06; | |
141 | fGamma[2][3] = -1.353e-07; | |
142 | ||
143 | fGamma[3][0] = 0.265004; | |
144 | fGamma[3][1] = 0.061298; | |
145 | fGamma[3][2] = -0.003203; | |
146 | fGamma[3][3] = 4.73e-05; | |
147 | ||
148 | fGamma[4][0] = 0.243579; | |
149 | fGamma[4][1] = -1.614e-05; | |
150 | ||
151 | fGamma[5][0] = 0.002942; | |
152 | fGamma[5][1] = -3.976e-05; | |
153 | ||
154 | fHadron[0][0] = 0.011945 / 3.; | |
155 | fHadron[0][1] = 0.000386 / 3.; | |
156 | fHadron[0][2] = -0.000014 / 3.; | |
157 | fHadron[0][3] = 1.336e-07 / 3.; | |
158 | ||
159 | fHadron[1][0] = 0.496544; | |
160 | fHadron[1][1] = -0.003226; | |
161 | fHadron[1][2] = 0.00001678; | |
162 | ||
163 | fHadron[2][0] = 0.144838; | |
164 | fHadron[2][1] = -0.002954; | |
165 | fHadron[2][2] = 0.00008754; | |
166 | fHadron[2][3] = -7.587e-07; | |
167 | ||
168 | fHadron[3][0] = 1.264461 / 7.; | |
169 | fHadron[3][1] = 0.002097 / 7.; | |
170 | ||
171 | fHadron[4][0] = 0.261950; | |
172 | fHadron[4][1] = -0.001078; | |
173 | fHadron[4][2] = 0.00003237; | |
174 | fHadron[4][3] = -3.241e-07; | |
175 | fHadron[4][4] = 0.; | |
176 | fHadron[4][5] = 0.; | |
177 | fHadron[5][0] = 0.010317; | |
178 | fHadron[5][1] = 0.; | |
179 | fHadron[5][2] = 0.; | |
180 | fHadron[5][3] = 0.; | |
181 | fHadron[5][4] = 0.; | |
182 | fHadron[5][5] = 0.; | |
183 | ||
184 | fPiZero5to10[0][0] = 0.009138; | |
185 | fPiZero5to10[0][1] = 0.0006377; | |
186 | ||
187 | fPiZero5to10[1][0] = 0.08; | |
188 | ||
189 | fPiZero5to10[2][0] = -0.061119; | |
190 | fPiZero5to10[2][1] = 0.019013; | |
191 | ||
192 | fPiZero5to10[3][0] = 0.2; | |
193 | ||
194 | fPiZero5to10[4][0] = 0.252044; | |
195 | fPiZero5to10[4][1] = -0.002315; | |
196 | ||
197 | fPiZero5to10[5][0] = 0.002942; | |
198 | fPiZero5to10[5][1] = -3.976e-05; | |
199 | ||
200 | fPiZero10to60[0][0] = 0.009138; | |
201 | fPiZero10to60[0][1] = 0.0006377; | |
202 | ||
203 | fPiZero10to60[1][0] = 1.272837; | |
204 | fPiZero10to60[1][1] = -0.069708; | |
205 | fPiZero10to60[1][2] = 0.001568; | |
206 | fPiZero10to60[1][3] = -1.162e-05; | |
207 | ||
208 | fPiZero10to60[2][0] = 0.139703; | |
209 | fPiZero10to60[2][1] = 0.003687; | |
210 | fPiZero10to60[2][2] = -0.000568; | |
211 | fPiZero10to60[2][3] = 1.498e-05; | |
212 | fPiZero10to60[2][4] = -1.174e-07; | |
213 | ||
214 | fPiZero10to60[3][0] = -0.826367; | |
215 | fPiZero10to60[3][1] = 0.096951; | |
216 | fPiZero10to60[3][2] = -0.002215; | |
217 | fPiZero10to60[3][3] = 2.523e-05; | |
218 | ||
219 | fPiZero10to60[4][0] = 0.249890; | |
220 | fPiZero10to60[4][1] = -0.000063; | |
221 | ||
222 | fPiZero10to60[5][0] = 0.002942; | |
223 | fPiZero10to60[5][1] = -3.976e-05; | |
224 | ||
225 | fPIDWeight[0] = -1; | |
226 | fPIDWeight[1] = -1; | |
227 | fPIDWeight[2] = -1; | |
228 | fReconstructor = kFALSE; | |
dc293ae9 | 229 | } |
4cfd8ef9 | 230 | |
231 | //______________________________________________ | |
dc293ae9 | 232 | void AliEMCALPID::RunPID(AliESD *esd) |
233 | { | |
234 | // | |
235 | // Make the PID for all the EMCAL clusters containedin the ESDs File | |
236 | // but just gamma/PiO/Hadron | |
237 | // | |
238 | // trivial check against NULL object passed | |
4cfd8ef9 | 239 | |
dc293ae9 | 240 | if (esd == 0x0) { |
4ce73a00 | 241 | AliInfo("NULL ESD object passed !!" ); |
dc293ae9 | 242 | return ; |
243 | } | |
4cfd8ef9 | 244 | |
dc293ae9 | 245 | Int_t nClusters = esd->GetNumberOfEMCALClusters(); |
246 | Int_t firstCluster = esd->GetFirstEMCALCluster(); | |
247 | Double_t energy, lambda0; | |
248 | for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) { | |
249 | ||
250 | AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster); | |
251 | energy = clust->GetClusterEnergy(); | |
252 | lambda0 = clust->GetM02(); | |
253 | // verify cluster type | |
254 | Int_t clusterType= clust->GetClusterType(); | |
4ce73a00 | 255 | if (clusterType == AliESDCaloCluster::kClusterv1 && lambda0 != 0 && energy < 1000) { |
256 | ||
257 | ||
dc293ae9 | 258 | // reject clusters with lambda0 = 0 |
4ce73a00 | 259 | |
260 | ||
dc293ae9 | 261 | ComputePID(energy, lambda0); |
4ce73a00 | 262 | |
263 | ||
dc293ae9 | 264 | if (fPrintInfo) { |
265 | AliInfo("___________________________________________________"); | |
266 | AliInfo(Form( "Particle Energy = %f",energy)); | |
267 | AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) ); | |
268 | AliInfo("PIDWeight of the particle :" ); | |
269 | AliInfo(Form( " GAMMA : %f",fPID[0] )); | |
270 | AliInfo(Form( " PiZero : %f",fPID[1] )); | |
271 | AliInfo(Form( " HADRON : %f", fPID[2] )); | |
272 | AliInfo("_________________________________________"); | |
273 | AliInfo(Form( " kElectron : %f", fPIDFinal[0]) ); | |
274 | AliInfo(Form( " kMuon : %f", fPIDFinal[1] )); | |
275 | AliInfo(Form( " kPion : %f", fPIDFinal[2] )); | |
276 | AliInfo(Form( " kKaon : %f", fPIDFinal[3] )); | |
277 | AliInfo(Form( " kProton : %f", fPIDFinal[4] )); | |
278 | AliInfo(Form( " kPhoton : %f", fPIDFinal[5] )); | |
279 | AliInfo(Form( " kPi0 : %f", fPIDFinal[6] )); | |
280 | AliInfo(Form( " kNeutron : %f", fPIDFinal[7] )); | |
281 | AliInfo(Form( " kKaon0 : %f", fPIDFinal[8] )); | |
282 | AliInfo(Form( " kEleCon : %f", fPIDFinal[9] )); | |
283 | AliInfo(Form( " kUnknown : %f", fPIDFinal[10] )); | |
284 | AliInfo("___________________________________________________"); | |
285 | } | |
4cfd8ef9 | 286 | |
dc293ae9 | 287 | if(fReconstructor) // In case it is called during reconstruction. |
288 | clust->SetPid(fPIDFinal); | |
289 | } // end if (clusterType...) | |
290 | } // end for (iCluster...) | |
291 | } | |
4cfd8ef9 | 292 | |
293 | //__________________________________________________________ | |
dc293ae9 | 294 | void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0) |
295 | { | |
296 | // | |
297 | // This is the main command, which uses the distributions computed and parametrised, | |
298 | // and gives the PID by the bayesian method. | |
299 | // | |
4ce73a00 | 300 | |
301 | if (energy<5){energy =6;} | |
302 | ||
303 | ||
4cfd8ef9 | 304 | TArrayD paramDistribGamma = DistLambda0(energy, 1); |
305 | TArrayD paramDistribPiZero = DistLambda0(energy, 2); | |
306 | TArrayD paramDistribHadron = DistLambda0(energy, 3); | |
307 | ||
308 | Bool_t norm = kFALSE; | |
309 | ||
310 | fProbGamma = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0]; | |
311 | fProbGamma += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3]; | |
312 | fProbPiZero = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0]; | |
313 | fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3]; | |
314 | fProbHadron = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0]; | |
315 | fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3]; | |
316 | ||
317 | // compute PID Weight | |
318 | fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron); | |
319 | fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron); | |
320 | fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron); | |
321 | ||
322 | SetPID(fPIDWeight[0], 0); | |
323 | SetPID(fPIDWeight[1], 1); | |
324 | SetPID(fPIDWeight[2], 2); | |
325 | ||
326 | // sortie ecran pid Weight only for control (= in english ???) | |
327 | if (fPrintInfo) { | |
328 | AliInfo(Form( "Energy in loop = %f", energy) ); | |
329 | AliInfo(Form( "Lambda0 in loop = %f", lambda0) ); | |
330 | AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) ); | |
331 | // AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) ); | |
332 | AliInfo(Form( "fProbaPiZero = %f", fProbPiZero )); | |
333 | AliInfo(Form( "fProbaHadron = %f", fProbHadron) ); | |
334 | AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f", fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) ); | |
335 | AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] )); | |
336 | AliInfo("********************************************************" ); | |
337 | } | |
338 | ||
339 | fPIDFinal[0] = fPIDWeight[0]/2; | |
340 | fPIDFinal[1] = fPIDWeight[2]/8; | |
341 | fPIDFinal[2] = fPIDWeight[2]/8; | |
342 | fPIDFinal[3] = fPIDWeight[2]/8; | |
343 | fPIDFinal[4] = fPIDWeight[2]/8; | |
344 | fPIDFinal[5] = fPIDWeight[0]/2; | |
345 | fPIDFinal[6] = fPIDWeight[1] ; | |
346 | fPIDFinal[7] = fPIDWeight[2]/8; | |
347 | fPIDFinal[8] = fPIDWeight[2]/8; | |
348 | fPIDFinal[9] = fPIDWeight[2]/8; | |
349 | fPIDFinal[10] = fPIDWeight[2]/8; | |
dc293ae9 | 350 | } |
4cfd8ef9 | 351 | |
352 | //________________________________________________________ | |
dc293ae9 | 353 | TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type) |
354 | { | |
4cfd8ef9 | 355 | // |
356 | // Compute the values of the parametrised distributions using the data initialised before. | |
357 | // | |
358 | Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.; | |
359 | Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.; | |
360 | TArrayD distributionParam(6); | |
361 | ||
362 | switch (type) { | |
363 | case 1: | |
364 | constGauss = Polynomial(energy, fGamma[0]); | |
365 | meanGauss = Polynomial(energy, fGamma[1]); | |
366 | sigmaGauss = Polynomial(energy, fGamma[2]); | |
367 | constLandau = Polynomial(energy, fGamma[3]); | |
368 | mpvLandau = Polynomial(energy, fGamma[4]); | |
369 | sigmaLandau = Polynomial(energy, fGamma[5]); | |
370 | break; | |
371 | case 2: | |
372 | if (energy < 10) { | |
373 | constGauss = Polynomial(energy, fPiZero5to10[0]); | |
374 | meanGauss = Polynomial(energy, fPiZero5to10[1]); | |
375 | sigmaGauss = Polynomial(energy, fPiZero5to10[2]); | |
376 | constLandau = Polynomial(energy, fPiZero5to10[3]); | |
377 | mpvLandau = Polynomial(energy, fPiZero5to10[4]); | |
378 | sigmaLandau = Polynomial(energy, fPiZero5to10[5]); | |
379 | } | |
380 | else { | |
381 | constGauss = Polynomial(energy, fPiZero10to60[0]); | |
382 | meanGauss = Polynomial(energy, fPiZero10to60[1]); | |
383 | sigmaGauss = Polynomial(energy, fPiZero10to60[2]); | |
384 | constLandau = Polynomial(energy, fPiZero10to60[3]); | |
385 | mpvLandau = Polynomial(energy, fPiZero10to60[4]); | |
386 | sigmaLandau = Polynomial(energy, fPiZero10to60[5]); | |
387 | } | |
388 | break; | |
389 | case 3: | |
390 | constGauss = Polynomial(energy, fHadron[0]); | |
391 | meanGauss = Polynomial(energy, fHadron[1]); | |
392 | sigmaGauss = Polynomial(energy, fHadron[2]); | |
393 | constLandau = Polynomial(energy, fHadron[3]); | |
394 | mpvLandau = Polynomial(energy, fHadron[4]); | |
395 | sigmaLandau = Polynomial(energy, fHadron[5]); | |
396 | break; | |
397 | } | |
398 | ||
399 | distributionParam[0] = constGauss; | |
400 | distributionParam[1] = meanGauss; | |
401 | distributionParam[2] = sigmaGauss; | |
402 | distributionParam[3] = constLandau; | |
403 | distributionParam[4] = mpvLandau; | |
404 | distributionParam[5] = sigmaLandau; | |
405 | ||
406 | return distributionParam; | |
dc293ae9 | 407 | } |
4cfd8ef9 | 408 | |
409 | //_______________________________________________________ | |
dc293ae9 | 410 | Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params) |
411 | { | |
4cfd8ef9 | 412 | // |
413 | // Compute a polynomial for a given value of 'x' | |
414 | // with the array of parameters passed as the second arg | |
415 | // | |
416 | Double_t y; | |
417 | y = params[0]; | |
418 | y += params[1] * x; | |
419 | y += params[2] * x * x; | |
420 | y += params[3] * x * x * x; | |
421 | y += params[4] * x * x * x * x; | |
422 | y += params[5] * x * x * x * x * x; | |
423 | ||
424 | return y; | |
dc293ae9 | 425 | } |