6ad0bfa0 |
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 | |
b2a60966 |
16 | /* $Id$ */ |
17 | |
702ab87e |
18 | /* History of cvs commits: |
19 | * |
20 | * $Log$ |
9a2cdbdf |
21 | * Revision 1.113 2007/08/07 14:12:03 kharlov |
22 | * Quality assurance added (Yves Schutz) |
23 | * |
ddd1a39c |
24 | * Revision 1.112 2007/07/11 13:43:30 hristov |
25 | * New class AliESDEvent, backward compatibility with the old AliESD (Christian) |
26 | * |
af885e0f |
27 | * Revision 1.111 2007/05/04 14:49:29 policheh |
28 | * AliPHOSRecPoint inheritance from AliCluster |
29 | * |
9ee9f78d |
30 | * Revision 1.110 2007/04/24 10:08:03 kharlov |
31 | * Vertex extraction from GenHeader |
32 | * |
2c06dc7a |
33 | * Revision 1.109 2007/04/18 09:34:05 kharlov |
34 | * Geometry bug fixes |
35 | * |
e5b7b511 |
36 | * Revision 1.108 2007/04/16 09:03:37 kharlov |
37 | * Incedent angle correction fixed |
38 | * |
753b19cd |
39 | * Revision 1.107 2007/04/02 15:00:16 cvetan |
40 | * No more calls to gAlice in the reconstruction |
41 | * |
999f9a8f |
42 | * Revision 1.106 2007/04/01 15:40:15 kharlov |
43 | * Correction for actual vertex position implemented |
44 | * |
407d15b3 |
45 | * Revision 1.105 2007/03/06 06:57:46 kharlov |
46 | * DP:calculation of distance to CPV done in TSM |
47 | * |
26aa7e4a |
48 | * Revision 1.104 2006/12/15 10:46:26 hristov |
49 | * Using TMath::Abs instead of fabs |
50 | * |
e4df4b30 |
51 | * Revision 1.103 2006/09/07 18:31:08 kharlov |
52 | * Effective c++ corrections (T.Pocheptsov) |
53 | * |
3f7dbdb7 |
54 | * Revision 1.102 2006/01/23 17:51:48 hristov |
55 | * Using the recommended way of forward declarations for TVector and TMatrix (see v5-08-00 release notes). Additional clean-up |
56 | * |
e8d02863 |
57 | * Revision 1.101 2005/05/28 14:19:04 schutz |
58 | * Compilation warnings fixed by T.P. |
59 | * |
702ab87e |
60 | */ |
61 | |
6ad0bfa0 |
62 | //_________________________________________________________________________ |
b2a60966 |
63 | // Implementation version v1 of the PHOS particle identifier |
7acf6008 |
64 | // Particle identification based on the |
148b2bba |
65 | // - RCPV: distance from CPV recpoint to EMCA recpoint. |
66 | // - TOF |
67 | // - PCA: Principal Components Analysis.. |
68 | // The identified particle has an identification number corresponding |
69 | // to a 9 bits number: |
bc0c084c |
70 | // -Bit 0 to 2: bit set if RCPV > CpvEmcDistance (each bit corresponds |
148b2bba |
71 | // to a different efficiency-purity point of the photon identification) |
bc0c084c |
72 | // -Bit 3 to 5: bit set if TOF < TimeGate (each bit corresponds |
148b2bba |
73 | // to a different efficiency-purity point of the photon identification) |
74 | // -Bit 6 to 9: bit set if Principal Components are |
50739f15 |
75 | // inside an ellipse defined by the parameters a, b, c, x0 and y0. |
148b2bba |
76 | // (each bit corresponds to a different efficiency-purity point of the |
50739f15 |
77 | // photon identification) |
78 | // The PCA (Principal components analysis) needs a file that contains |
79 | // a previous analysis of the correlations between the particles. This |
bc0c084c |
80 | // file is $ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root. Analysis done for |
50739f15 |
81 | // energies between 0.5 and 100 GeV. |
9fa5f1d0 |
82 | // A calibrated energy is calculated. The energy of the reconstructed |
50739f15 |
83 | // cluster is corrected with the formula A + B * E + C * E^2, whose |
bc0c084c |
84 | // parameters where obtained through the study of the reconstructed |
50739f15 |
85 | // energy distribution of monoenergetic photons. |
a4e98857 |
86 | // |
bc0c084c |
87 | // All the parameters (RCPV(2 rows-3 columns),TOF(1r-3c),PCA(5r-4c) |
50739f15 |
88 | // and calibration(1r-3c))are stored in a file called |
89 | // $ALICE_ROOT/PHOS/Parameters.dat. Each time that AliPHOSPIDv1 is |
bc0c084c |
90 | // initialized, this parameters are copied to a Matrix (9,4), a |
50739f15 |
91 | // TMatrixD object. |
7acf6008 |
92 | // |
a4e98857 |
93 | // use case: |
50739f15 |
94 | // root [0] AliPHOSPIDv1 * p = new AliPHOSPIDv1("galice1.root") |
a4e98857 |
95 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
50739f15 |
96 | // // reading headers from file galice1.root and create RecParticles |
97 | // TrackSegments and RecPoints are used |
98 | // // set file name for the branch RecParticles |
f0a4c9e9 |
99 | // root [1] p->ExecuteTask("deb all time") |
50739f15 |
100 | // // available options |
101 | // // "deb" - prints # of reconstructed particles |
102 | // // "deb all" - prints # and list of RecParticles |
103 | // // "time" - prints benchmarking results |
7acf6008 |
104 | // |
50739f15 |
105 | // root [2] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1",kTRUE) |
148b2bba |
106 | // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated |
50739f15 |
107 | // //Split mode. |
f0a4c9e9 |
108 | // root [3] p2->ExecuteTask() |
109 | // |
50739f15 |
110 | |
f0a4c9e9 |
111 | |
7acf6008 |
112 | //*-- Author: Yves Schutz (SUBATECH) & Gines Martinez (SUBATECH) & |
148b2bba |
113 | // Gustavo Conesa April 2002 |
50739f15 |
114 | // PCA redesigned by Gustavo Conesa October 2002: |
115 | // The way of using the PCA has changed. Instead of 2 |
116 | // files with the PCA, each one with different energy ranges |
117 | // of application, we use the wide one (0.5-100 GeV), and instead |
bc0c084c |
118 | // of fixing 3 ellipses for different ranges of energy, it has been |
50739f15 |
119 | // studied the dependency of the ellipses parameters with the |
120 | // energy, and they are implemented in the code as a funtion |
121 | // of the energy. |
122 | // |
123 | // |
124 | // |
6ad0bfa0 |
125 | // --- ROOT system --- |
c947e71a |
126 | |
127 | |
128 | // --- Standard library --- |
e8d02863 |
129 | #include <TMatrixF.h> |
acb5beb7 |
130 | #include "TFormula.h" |
7acf6008 |
131 | #include "TBenchmark.h" |
148b2bba |
132 | #include "TPrincipal.h" |
c947e71a |
133 | #include "TFile.h" |
e3817e5f |
134 | #include "TSystem.h" |
148b2bba |
135 | |
6ad0bfa0 |
136 | // --- AliRoot header files --- |
c947e71a |
137 | //#include "AliLog.h" |
e3817e5f |
138 | #include "AliPHOS.h" |
26d4b141 |
139 | #include "AliPHOSPIDv1.h" |
af885e0f |
140 | #include "AliESDEvent.h" |
407d15b3 |
141 | #include "AliESDVertex.h" |
9a2cdbdf |
142 | #include "AliPHOSTrackSegment.h" |
143 | #include "AliPHOSEmcRecPoint.h" |
144 | #include "AliPHOSRecParticle.h" |
6ad0bfa0 |
145 | |
26d4b141 |
146 | ClassImp( AliPHOSPIDv1) |
6ad0bfa0 |
147 | |
1cb7c1ee |
148 | //____________________________________________________________________________ |
3f7dbdb7 |
149 | AliPHOSPIDv1::AliPHOSPIDv1() : |
ddd1a39c |
150 | AliPHOSPID(), |
3f7dbdb7 |
151 | fBayesian(kFALSE), |
152 | fDefaultInit(kFALSE), |
153 | fWrite(kFALSE), |
3f7dbdb7 |
154 | fFileNamePrincipalPhoton(), |
155 | fFileNamePrincipalPi0(), |
156 | fFileNameParameters(), |
157 | fPrincipalPhoton(0), |
158 | fPrincipalPi0(0), |
159 | fX(0), |
160 | fPPhoton(0), |
161 | fPPi0(0), |
3f7dbdb7 |
162 | fParameters(0), |
ddd1a39c |
163 | fVtx(0.), |
3f7dbdb7 |
164 | fTFphoton(0), |
165 | fTFpiong(0), |
166 | fTFkaong(0), |
167 | fTFkaonl(0), |
168 | fTFhhadrong(0), |
169 | fTFhhadronl(0), |
170 | fDFmuon(0), |
171 | fERecWeight(0), |
172 | fChargedNeutralThreshold(0.), |
173 | fTOFEnThreshold(0), |
174 | fDispEnThreshold(0), |
175 | fDispMultThreshold(0) |
1cb7c1ee |
176 | { |
a4e98857 |
177 | // default ctor |
148b2bba |
178 | |
8d0f3f77 |
179 | InitParameters() ; |
92f521a9 |
180 | fDefaultInit = kTRUE ; |
7acf6008 |
181 | } |
182 | |
581354c5 |
183 | //____________________________________________________________________________ |
3f7dbdb7 |
184 | AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ) : |
185 | AliPHOSPID(pid), |
186 | fBayesian(kFALSE), |
187 | fDefaultInit(kFALSE), |
188 | fWrite(kFALSE), |
3f7dbdb7 |
189 | fFileNamePrincipalPhoton(), |
190 | fFileNamePrincipalPi0(), |
191 | fFileNameParameters(), |
192 | fPrincipalPhoton(0), |
193 | fPrincipalPi0(0), |
194 | fX(0), |
195 | fPPhoton(0), |
196 | fPPi0(0), |
3f7dbdb7 |
197 | fParameters(0), |
ddd1a39c |
198 | fVtx(0.), |
3f7dbdb7 |
199 | fTFphoton(0), |
200 | fTFpiong(0), |
201 | fTFkaong(0), |
202 | fTFkaonl(0), |
203 | fTFhhadrong(0), |
204 | fTFhhadronl(0), |
205 | fDFmuon(0), |
206 | fERecWeight(0), |
207 | fChargedNeutralThreshold(0.), |
208 | fTOFEnThreshold(0), |
209 | fDispEnThreshold(0), |
210 | fDispMultThreshold(0) |
211 | |
581354c5 |
212 | { |
386aef34 |
213 | // ctor |
581354c5 |
214 | InitParameters() ; |
581354c5 |
215 | |
216 | } |
217 | |
7acf6008 |
218 | //____________________________________________________________________________ |
9a2cdbdf |
219 | AliPHOSPIDv1::AliPHOSPIDv1(AliPHOSGeometry *geom): |
220 | AliPHOSPID(geom), |
3f7dbdb7 |
221 | fBayesian(kFALSE), |
222 | fDefaultInit(kFALSE), |
223 | fWrite(kFALSE), |
3f7dbdb7 |
224 | fFileNamePrincipalPhoton(), |
225 | fFileNamePrincipalPi0(), |
226 | fFileNameParameters(), |
227 | fPrincipalPhoton(0), |
228 | fPrincipalPi0(0), |
229 | fX(0), |
230 | fPPhoton(0), |
231 | fPPi0(0), |
3f7dbdb7 |
232 | fParameters(0), |
ddd1a39c |
233 | fVtx(0.), |
3f7dbdb7 |
234 | fTFphoton(0), |
235 | fTFpiong(0), |
236 | fTFkaong(0), |
237 | fTFkaonl(0), |
238 | fTFhhadrong(0), |
239 | fTFhhadronl(0), |
240 | fDFmuon(0), |
241 | fERecWeight(0), |
242 | fChargedNeutralThreshold(0.), |
243 | fTOFEnThreshold(0), |
244 | fDispEnThreshold(0), |
245 | fDispMultThreshold(0) |
246 | |
7acf6008 |
247 | { |
a4e98857 |
248 | //ctor with the indication on where to look for the track segments |
7b7c1533 |
249 | |
8d0f3f77 |
250 | InitParameters() ; |
92f521a9 |
251 | fDefaultInit = kFALSE ; |
7acf6008 |
252 | } |
7b7c1533 |
253 | |
7acf6008 |
254 | //____________________________________________________________________________ |
255 | AliPHOSPIDv1::~AliPHOSPIDv1() |
256 | { |
79bb1b62 |
257 | // dtor |
acb5beb7 |
258 | fPrincipalPhoton = 0; |
259 | fPrincipalPi0 = 0; |
9fa5f1d0 |
260 | |
e3817e5f |
261 | delete [] fX ; // Principal input |
262 | delete [] fPPhoton ; // Photon Principal components |
263 | delete [] fPPi0 ; // Pi0 Principal components |
acb5beb7 |
264 | |
265 | delete fParameters; |
266 | delete fTFphoton; |
267 | delete fTFpiong; |
268 | delete fTFkaong; |
269 | delete fTFkaonl; |
270 | delete fTFhhadrong; |
271 | delete fTFhhadronl; |
272 | delete fDFmuon; |
7acf6008 |
273 | } |
a496c46c |
274 | |
8d0f3f77 |
275 | //____________________________________________________________________________ |
276 | void AliPHOSPIDv1::InitParameters() |
277 | { |
e3817e5f |
278 | // Initialize PID parameters |
adcca1e6 |
279 | fWrite = kTRUE ; |
35adb638 |
280 | fBayesian = kTRUE ; |
9fa5f1d0 |
281 | SetParameters() ; // fill the parameters matrix from parameters file |
35adb638 |
282 | |
cc1fe362 |
283 | // initialisation of response function parameters |
284 | // Tof |
2924941c |
285 | |
286 | // // Photons |
287 | // fTphoton[0] = 0.218 ; |
288 | // fTphoton[1] = 1.55E-8 ; |
289 | // fTphoton[2] = 5.05E-10 ; |
290 | // fTFphoton = new TFormula("ToF response to photons" , "gaus") ; |
291 | // fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; |
292 | |
293 | // // Pions |
294 | // //Gaus (0 to max probability) |
295 | // fTpiong[0] = 0.0971 ; |
296 | // fTpiong[1] = 1.58E-8 ; |
297 | // fTpiong[2] = 5.69E-10 ; |
298 | // fTFpiong = new TFormula("ToF response to pions" , "gaus") ; |
299 | // fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; |
300 | |
301 | // // Kaons |
302 | // //Gaus (0 to max probability) |
303 | // fTkaong[0] = 0.0542 ; |
304 | // fTkaong[1] = 1.64E-8 ; |
305 | // fTkaong[2] = 6.07E-10 ; |
306 | // fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; |
307 | // fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; |
308 | // //Landau (max probability to inf) |
309 | // fTkaonl[0] = 0.264 ; |
310 | // fTkaonl[1] = 1.68E-8 ; |
311 | // fTkaonl[2] = 4.10E-10 ; |
312 | // fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; |
313 | // fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; |
314 | |
315 | // //Heavy Hadrons |
316 | // //Gaus (0 to max probability) |
317 | // fThhadrong[0] = 0.0302 ; |
318 | // fThhadrong[1] = 1.73E-8 ; |
319 | // fThhadrong[2] = 9.52E-10 ; |
320 | // fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; |
321 | // fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; |
322 | // //Landau (max probability to inf) |
323 | // fThhadronl[0] = 0.139 ; |
324 | // fThhadronl[1] = 1.745E-8 ; |
325 | // fThhadronl[2] = 1.00E-9 ; |
326 | // fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; |
327 | // fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; |
328 | |
cc1fe362 |
329 | // Photons |
2924941c |
330 | fTphoton[0] = 7.83E8 ; |
35adb638 |
331 | fTphoton[1] = 1.55E-8 ; |
2924941c |
332 | fTphoton[2] = 5.09E-10 ; |
35adb638 |
333 | fTFphoton = new TFormula("ToF response to photons" , "gaus") ; |
cc1fe362 |
334 | fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ; |
35adb638 |
335 | |
336 | // Pions |
337 | //Gaus (0 to max probability) |
2924941c |
338 | fTpiong[0] = 6.73E8 ; |
35adb638 |
339 | fTpiong[1] = 1.58E-8 ; |
2924941c |
340 | fTpiong[2] = 5.87E-10 ; |
35adb638 |
341 | fTFpiong = new TFormula("ToF response to pions" , "gaus") ; |
342 | fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ; |
35adb638 |
343 | |
344 | // Kaons |
345 | //Gaus (0 to max probability) |
2924941c |
346 | fTkaong[0] = 3.93E8 ; |
35adb638 |
347 | fTkaong[1] = 1.64E-8 ; |
2924941c |
348 | fTkaong[2] = 6.07E-10 ; |
35adb638 |
349 | fTFkaong = new TFormula("ToF response to kaon" , "gaus") ; |
350 | fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ; |
351 | //Landau (max probability to inf) |
2924941c |
352 | fTkaonl[0] = 2.0E9 ; |
35adb638 |
353 | fTkaonl[1] = 1.68E-8 ; |
354 | fTkaonl[2] = 4.10E-10 ; |
355 | fTFkaonl = new TFormula("ToF response to kaon" , "landau") ; |
356 | fTFkaonl->SetParameters( fTkaonl[0], fTkaonl[1], fTkaonl[2]) ; |
357 | |
358 | //Heavy Hadrons |
359 | //Gaus (0 to max probability) |
2924941c |
360 | fThhadrong[0] = 2.02E8 ; |
35adb638 |
361 | fThhadrong[1] = 1.73E-8 ; |
362 | fThhadrong[2] = 9.52E-10 ; |
363 | fTFhhadrong = new TFormula("ToF response to heavy hadrons" , "gaus") ; |
364 | fTFhhadrong->SetParameters( fThhadrong[0], fThhadrong[1], fThhadrong[2]) ; |
365 | //Landau (max probability to inf) |
2924941c |
366 | fThhadronl[0] = 1.10E9 ; |
367 | fThhadronl[1] = 1.74E-8 ; |
368 | fThhadronl[2] = 1.00E-9 ; |
35adb638 |
369 | fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ; |
370 | fTFhhadronl->SetParameters( fThhadronl[0], fThhadronl[1], fThhadronl[2]) ; |
371 | |
35adb638 |
372 | |
373 | |
374 | // Shower shape: dispersion gaussian parameters |
375 | // Photons |
2924941c |
376 | |
377 | // fDphoton[0] = 4.62e-2; fDphoton[1] = 1.39e-2 ; fDphoton[2] = -3.80e-2;//constant |
378 | // fDphoton[3] = 1.53 ; fDphoton[4] =-6.62e-2 ; fDphoton[5] = 0.339 ;//mean |
379 | // fDphoton[6] = 6.89e-2; fDphoton[7] =-6.59e-2 ; fDphoton[8] = 0.194 ;//sigma |
380 | |
381 | // fDpi0[0] = 0.0586 ; fDpi0[1] = 1.06E-3 ; fDpi0[2] = 0. ;//constant |
382 | // fDpi0[3] = 2.67 ; fDpi0[4] =-2.00E-2 ; fDpi0[5] = 9.37E-5 ;//mean |
383 | // fDpi0[6] = 0.153 ; fDpi0[7] = 9.34E-4 ; fDpi0[8] =-1.49E-5 ;//sigma |
384 | |
385 | // fDhadron[0] = 1.61E-2 ; fDhadron[1] = 3.03E-3 ; fDhadron[2] = 1.01E-2 ;//constant |
386 | // fDhadron[3] = 3.81 ; fDhadron[4] = 0.232 ; fDhadron[5] =-1.25 ;//mean |
387 | // fDhadron[6] = 0.897 ; fDhadron[7] = 0.0987 ; fDhadron[8] =-0.534 ;//sigma |
388 | |
389 | fDphoton[0] = 1.5 ; fDphoton[1] = 0.49 ; fDphoton[2] =-1.7E-2 ;//constant |
390 | fDphoton[3] = 1.5 ; fDphoton[4] = 4.0E-2 ; fDphoton[5] = 0.21 ;//mean |
391 | fDphoton[6] = 4.8E-2 ; fDphoton[7] =-0.12 ; fDphoton[8] = 0.27 ;//sigma |
392 | fDphoton[9] = 16.; //for E> fDphoton[9] parameters calculated at fDphoton[9] |
393 | |
394 | fDpi0[0] = 0.25 ; fDpi0[1] = 3.3E-2 ; fDpi0[2] =-1.0e-5 ;//constant |
395 | fDpi0[3] = 1.50 ; fDpi0[4] = 398. ; fDpi0[5] = 12. ;//mean |
396 | fDpi0[6] =-7.0E-2 ; fDpi0[7] =-524. ; fDpi0[8] = 22. ;//sigma |
397 | fDpi0[9] = 110.; //for E> fDpi0[9] parameters calculated at fDpi0[9] |
398 | |
399 | fDhadron[0] = 6.5 ; fDhadron[1] =-5.3 ; fDhadron[2] = 1.5 ;//constant |
400 | fDhadron[3] = 3.8 ; fDhadron[4] = 0.23 ; fDhadron[5] =-1.2 ;//mean |
401 | fDhadron[6] = 0.88 ; fDhadron[7] = 9.3E-2 ; fDhadron[8] =-0.51 ;//sigma |
402 | fDhadron[9] = 2.; //for E> fDhadron[9] parameters calculated at fDhadron[9] |
403 | |
404 | fDmuon[0] = 0.0631 ; |
405 | fDmuon[1] = 1.4 ; |
35adb638 |
406 | fDmuon[2] = 0.0557 ; |
407 | fDFmuon = new TFormula("Shower shape response to muons" , "landau") ; |
408 | fDFmuon->SetParameters( fDmuon[0], fDmuon[1], fDmuon[2]) ; |
409 | |
35adb638 |
410 | |
c947e71a |
411 | // x(CPV-EMC) distance gaussian parameters |
412 | |
2924941c |
413 | // fXelectron[0] = 8.06e-2 ; fXelectron[1] = 1.00e-2; fXelectron[2] =-5.14e-2;//constant |
414 | // fXelectron[3] = 0.202 ; fXelectron[4] = 8.15e-3; fXelectron[5] = 4.55 ;//mean |
415 | // fXelectron[6] = 0.334 ; fXelectron[7] = 0.186 ; fXelectron[8] = 4.32e-2;//sigma |
c947e71a |
416 | |
2924941c |
417 | // //charged hadrons gaus |
418 | // fXcharged[0] = 6.43e-3 ; fXcharged[1] =-4.19e-5; fXcharged[2] = 1.42e-3;//constant |
419 | // fXcharged[3] = 2.75 ; fXcharged[4] =-0.40 ; fXcharged[5] = 1.68 ;//mean |
420 | // fXcharged[6] = 3.135 ; fXcharged[7] =-9.41e-2; fXcharged[8] = 1.31e-2;//sigma |
c947e71a |
421 | |
2924941c |
422 | // // z(CPV-EMC) distance gaussian parameters |
423 | |
424 | // fZelectron[0] = 8.22e-2 ; fZelectron[1] = 5.11e-3; fZelectron[2] =-3.05e-2;//constant |
425 | // fZelectron[3] = 3.09e-2 ; fZelectron[4] = 5.87e-2; fZelectron[5] =-9.49e-2;//mean |
426 | // fZelectron[6] = 0.263 ; fZelectron[7] =-9.02e-3; fZelectron[8] = 0.151 ;//sigma |
c947e71a |
427 | |
2924941c |
428 | // //charged hadrons gaus |
c947e71a |
429 | |
2924941c |
430 | // fZcharged[0] = 1.00e-2 ; fZcharged[1] = 2.82E-4 ; fZcharged[2] = 2.87E-3 ;//constant |
431 | // fZcharged[3] =-4.68e-2 ; fZcharged[4] =-9.21e-3 ; fZcharged[5] = 4.91e-2 ;//mean |
432 | // fZcharged[6] = 1.425 ; fZcharged[7] =-5.90e-2 ; fZcharged[8] = 5.07e-2 ;//sigma |
433 | |
434 | |
435 | fXelectron[0] =-1.6E-2 ; fXelectron[1] = 0.77 ; fXelectron[2] =-0.15 ;//constant |
436 | fXelectron[3] = 0.35 ; fXelectron[4] = 0.25 ; fXelectron[5] = 4.12 ;//mean |
437 | fXelectron[6] = 0.30 ; fXelectron[7] = 0.11 ; fXelectron[8] = 0.16 ;//sigma |
438 | fXelectron[9] = 3.; //for E> fXelectron[9] parameters calculated at fXelectron[9] |
439 | |
c947e71a |
440 | //charged hadrons gaus |
2924941c |
441 | fXcharged[0] = 0.14 ; fXcharged[1] =-3.0E-2 ; fXcharged[2] = 0 ;//constant |
442 | fXcharged[3] = 1.4 ; fXcharged[4] =-9.3E-2 ; fXcharged[5] = 1.4 ;//mean |
443 | fXcharged[6] = 5.7 ; fXcharged[7] = 0.27 ; fXcharged[8] =-1.8 ;//sigma |
444 | fXcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] |
445 | |
446 | // z(CPV-EMC) distance gaussian parameters |
c947e71a |
447 | |
2924941c |
448 | fZelectron[0] = 0.49 ; fZelectron[1] = 0.53 ; fZelectron[2] =-9.8E-2 ;//constant |
449 | fZelectron[3] = 2.8E-2 ; fZelectron[4] = 5.0E-2 ; fZelectron[5] =-8.2E-2 ;//mean |
450 | fZelectron[6] = 0.25 ; fZelectron[7] =-1.7E-2 ; fZelectron[8] = 0.17 ;//sigma |
451 | fZelectron[9] = 3.; //for E> fZelectron[9] parameters calculated at fZelectron[9] |
452 | |
453 | //charged hadrons gaus |
c947e71a |
454 | |
2924941c |
455 | fZcharged[0] = 0.46 ; fZcharged[1] =-0.65 ; fZcharged[2] = 0.52 ;//constant |
456 | fZcharged[3] = 1.1E-2 ; fZcharged[4] = 0. ; fZcharged[5] = 0. ;//mean |
457 | fZcharged[6] = 0.60 ; fZcharged[7] =-8.2E-2 ; fZcharged[8] = 0.45 ;//sigma |
458 | fZcharged[9] = 1.2; //for E> fXcharged[9] parameters calculated at fXcharged[9] |
459 | |
fb7b51ad |
460 | //Threshold to differentiate between charged and neutral |
461 | fChargedNeutralThreshold = 1e-5; |
2924941c |
462 | fTOFEnThreshold = 2; //Maximum energy to use TOF |
463 | fDispEnThreshold = 0.5; //Minimum energy to use shower shape |
464 | fDispMultThreshold = 3; //Minimum multiplicity to use shower shape |
fb7b51ad |
465 | |
466 | //Weight to hadrons recontructed energy |
467 | |
468 | fERecWeightPar[0] = 0.32 ; |
469 | fERecWeightPar[1] = 3.8 ; |
470 | fERecWeightPar[2] = 5.4E-3 ; |
471 | fERecWeightPar[3] = 5.6E-2 ; |
472 | fERecWeight = new TFormula("Weight for hadrons" , "[0]*exp(-x*[1])+[2]*exp(-x*[3])") ; |
473 | fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ; |
474 | |
475 | |
304864ab |
476 | for (Int_t i =0; i< AliPID::kSPECIESN ; i++) |
35adb638 |
477 | fInitPID[i] = 1.; |
fb7b51ad |
478 | |
8d0f3f77 |
479 | } |
480 | |
88cb7938 |
481 | //________________________________________________________________________ |
9a2cdbdf |
482 | void AliPHOSPIDv1::TrackSegments2RecParticles(Option_t *option) |
88cb7938 |
483 | { |
eabde521 |
484 | // Steering method to perform particle reconstruction and identification |
485 | // for the event range from fFirstEvent to fLastEvent. |
61d3d6aa |
486 | |
88cb7938 |
487 | if(strstr(option,"tim")) |
488 | gBenchmark->Start("PHOSPID"); |
489 | |
490 | if(strstr(option,"print")) { |
491 | Print() ; |
492 | return ; |
493 | } |
494 | |
9a2cdbdf |
495 | if(fTrackSegments && //Skip events, where no track segments made |
496 | fTrackSegments->GetEntriesFast()) { |
88cb7938 |
497 | |
9a2cdbdf |
498 | GetVertex() ; |
499 | MakeRecParticles() ; |
500 | |
501 | if(fBayesian) |
502 | MakePID() ; |
35adb638 |
503 | |
9a2cdbdf |
504 | if(strstr(option,"deb")) |
505 | PrintRecParticles(option) ; |
88cb7938 |
506 | } |
9a2cdbdf |
507 | |
ff417097 |
508 | if(strstr(option,"deb")) |
509 | PrintRecParticles(option); |
88cb7938 |
510 | if(strstr(option,"tim")){ |
511 | gBenchmark->Stop("PHOSPID"); |
9a2cdbdf |
512 | AliInfo(Form("took %f seconds for PID", |
513 | gBenchmark->GetCpuTime("PHOSPID"))); |
88cb7938 |
514 | } |
88cb7938 |
515 | } |
516 | |
35adb638 |
517 | //________________________________________________________________________ |
17323043 |
518 | Double_t AliPHOSPIDv1::GausF(Double_t x, Double_t y, Double_t * par) |
35adb638 |
519 | { |
c947e71a |
520 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
521 | //this method returns a density probability of this parameter, given by a gaussian |
522 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b |
2924941c |
523 | //Float_t xorg = x; |
524 | if (x > par[9]) x = par[9]; |
525 | |
526 | //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; |
527 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
35adb638 |
528 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
529 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; |
c947e71a |
530 | |
2924941c |
531 | // if(xorg > 30) |
532 | // cout<<"En_in = "<<xorg<<"; En_out = "<<x<<"; cnt = "<<cnt |
533 | // <<"; mean = "<<mean<<"; sigma = "<<sigma<<endl; |
534 | |
35adb638 |
535 | // Double_t arg = - (y-mean) * (y-mean) / (2*sigma*sigma) ; |
536 | // return cnt * TMath::Exp(arg) ; |
c947e71a |
537 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 |
538 | return cnt*TMath::Gaus(y,mean,sigma); |
35adb638 |
539 | } |
540 | else |
541 | return 0.; |
c947e71a |
542 | |
35adb638 |
543 | } |
544 | //________________________________________________________________________ |
17323043 |
545 | Double_t AliPHOSPIDv1::GausPol2(Double_t x, Double_t y, Double_t * par) |
35adb638 |
546 | { |
c947e71a |
547 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
548 | //this method returns a density probability of this parameter, given by a gaussian |
549 | //function whose parameters depend with the energy like second order polinomial |
550 | |
35adb638 |
551 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
552 | Double_t mean = par[3] + par[4] * x + par[5] * x * x ; |
553 | Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; |
554 | |
c947e71a |
555 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 |
556 | return cnt*TMath::Gaus(y,mean,sigma); |
35adb638 |
557 | } |
558 | else |
559 | return 0.; |
c947e71a |
560 | |
561 | |
562 | |
35adb638 |
563 | } |
564 | |
69183710 |
565 | //____________________________________________________________________________ |
e3817e5f |
566 | const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const |
148b2bba |
567 | { |
e3817e5f |
568 | //Get file name that contains the PCA for a particle ("photon or pi0") |
569 | particle.ToLower(); |
570 | TString name; |
351dd634 |
571 | if (particle=="photon") |
572 | name = fFileNamePrincipalPhoton ; |
573 | else if (particle=="pi0" ) |
574 | name = fFileNamePrincipalPi0 ; |
575 | else |
576 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
577 | particle.Data())); |
e3817e5f |
578 | return name; |
579 | } |
bc0c084c |
580 | |
e3817e5f |
581 | //____________________________________________________________________________ |
fc7e2f43 |
582 | Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const |
e3817e5f |
583 | { |
584 | // Get the i-th parameter "Calibration" |
585 | Float_t param = 0.; |
351dd634 |
586 | if (i>2 || i<0) { |
587 | AliError(Form("Invalid parameter number: %d",i)); |
588 | } else |
e3817e5f |
589 | param = (*fParameters)(0,i); |
590 | return param; |
591 | } |
bc0c084c |
592 | |
88cb7938 |
593 | //____________________________________________________________________________ |
fc7e2f43 |
594 | Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const |
88cb7938 |
595 | { |
596 | // It calibrates Energy depending on the recpoint energy. |
597 | // The energy of the reconstructed cluster is corrected with |
598 | // the formula A + B* E + C* E^2, whose parameters where obtained |
599 | // through the study of the reconstructed energy distribution of |
600 | // monoenergetic photons. |
71994f35 |
601 | |
602 | if(!fEnergyCorrectionOn) return e; |
603 | |
88cb7938 |
604 | Float_t p[]={0.,0.,0.}; |
605 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i); |
606 | Float_t enerec = p[0] + p[1]*e + p[2]*e*e; |
607 | return enerec ; |
608 | |
609 | } |
610 | |
e3817e5f |
611 | //____________________________________________________________________________ |
fc7e2f43 |
612 | Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const |
e3817e5f |
613 | { |
614 | // Get the i-th parameter "CPV-EMC distance" for the specified axis |
615 | Float_t param = 0.; |
351dd634 |
616 | if(i>2 || i<0) { |
617 | AliError(Form("Invalid parameter number: %d",i)); |
618 | } else { |
e3817e5f |
619 | axis.ToLower(); |
351dd634 |
620 | if (axis == "x") |
621 | param = (*fParameters)(1,i); |
622 | else if (axis == "z") |
623 | param = (*fParameters)(2,i); |
624 | else { |
625 | AliError(Form("Invalid axis name: %s",axis.Data())); |
626 | } |
e3817e5f |
627 | } |
628 | return param; |
629 | } |
630 | |
631 | //____________________________________________________________________________ |
fc7e2f43 |
632 | Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const |
e3817e5f |
633 | { |
88cb7938 |
634 | // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and |
635 | // Purity-Efficiency point |
636 | |
637 | axis.ToLower(); |
638 | Float_t p[]={0.,0.,0.}; |
639 | for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis); |
640 | Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e); |
641 | return sig; |
e3817e5f |
642 | } |
643 | |
88cb7938 |
644 | //____________________________________________________________________________ |
fc7e2f43 |
645 | Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const |
88cb7938 |
646 | { |
647 | // Calculates the parameter param of the ellipse |
e3817e5f |
648 | |
649 | particle.ToLower(); |
650 | param. ToLower(); |
88cb7938 |
651 | Float_t p[4]={0.,0.,0.,0.}; |
652 | Float_t value = 0.0; |
653 | for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i); |
654 | if (particle == "photon") { |
655 | if (param.Contains("a")) e = TMath::Min((Double_t)e,70.); |
656 | else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.); |
657 | else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1); |
658 | } |
e3817e5f |
659 | |
443caba9 |
660 | if (particle == "photon") |
661 | value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3]; |
662 | else if (particle == "pi0") |
663 | value = p[0] + p[1]*e + p[2]*e*e; |
664 | |
88cb7938 |
665 | return value; |
e3817e5f |
666 | } |
667 | |
668 | //_____________________________________________________________________________ |
fc7e2f43 |
669 | Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const |
e3817e5f |
670 | { |
671 | // Get the parameter "i" to calculate the boundary on the moment M2x |
672 | // for photons at high p_T |
673 | Float_t param = 0; |
351dd634 |
674 | if (i>3 || i<0) { |
675 | AliError(Form("Wrong parameter number: %d\n",i)); |
676 | } else |
e3817e5f |
677 | param = (*fParameters)(14,i) ; |
678 | return param; |
148b2bba |
679 | } |
e3817e5f |
680 | |
148b2bba |
681 | //____________________________________________________________________________ |
fc7e2f43 |
682 | Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const |
e3817e5f |
683 | { |
684 | // Get the parameter "i" to calculate the boundary on the moment M2x |
685 | // for pi0 at high p_T |
686 | Float_t param = 0; |
351dd634 |
687 | if (i>2 || i<0) { |
688 | AliError(Form("Wrong parameter number: %d\n",i)); |
689 | } else |
e3817e5f |
690 | param = (*fParameters)(15,i) ; |
691 | return param; |
692 | } |
148b2bba |
693 | |
e3817e5f |
694 | //____________________________________________________________________________ |
fc7e2f43 |
695 | Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const |
e3817e5f |
696 | { |
88cb7938 |
697 | // Get TimeGate parameter depending on Purity-Efficiency i: |
698 | // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity |
699 | Float_t param = 0.; |
351dd634 |
700 | if(i>2 || i<0) { |
701 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); |
702 | } else |
88cb7938 |
703 | param = (*fParameters)(3,i) ; |
704 | return param; |
e3817e5f |
705 | } |
706 | |
148b2bba |
707 | //_____________________________________________________________________________ |
fc7e2f43 |
708 | Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const |
88cb7938 |
709 | { |
710 | // Get the parameter "i" that is needed to calculate the ellipse |
711 | // parameter "param" for the particle "particle" ("photon" or "pi0") |
712 | |
e3817e5f |
713 | particle.ToLower(); |
714 | param. ToLower(); |
88cb7938 |
715 | Int_t offset = -1; |
351dd634 |
716 | if (particle == "photon") |
717 | offset=0; |
718 | else if (particle == "pi0") |
719 | offset=5; |
e3817e5f |
720 | else |
351dd634 |
721 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
722 | particle.Data())); |
88cb7938 |
723 | |
724 | Int_t p= -1; |
725 | Float_t par = 0; |
e3817e5f |
726 | |
727 | if (param.Contains("a")) p=4+offset; |
728 | else if(param.Contains("b")) p=5+offset; |
729 | else if(param.Contains("c")) p=6+offset; |
730 | else if(param.Contains("x0"))p=7+offset; |
731 | else if(param.Contains("y0"))p=8+offset; |
12022e83 |
732 | |
351dd634 |
733 | if (i>4 || i<0) { |
734 | AliError(Form("No parameter with index %d", i)) ; |
735 | } else if (p==-1) { |
736 | AliError(Form("No parameter with name %s", param.Data() )) ; |
737 | } else |
88cb7938 |
738 | par = (*fParameters)(p,i) ; |
739 | |
740 | return par; |
12022e83 |
741 | } |
742 | |
12022e83 |
743 | |
26aa7e4a |
744 | //DP____________________________________________________________________________ |
745 | //Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const |
746 | //{ |
747 | // // Calculates the distance between the EMC RecPoint and the PPSD RecPoint |
748 | // |
7fb9892d |
749 | // AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance(); |
26aa7e4a |
750 | // TVector3 vecEmc ; |
751 | // TVector3 vecCpv ; |
752 | // if(cpv){ |
753 | // emc->GetLocalPosition(vecEmc) ; |
754 | // cpv->GetLocalPosition(vecCpv) ; |
755 | // |
756 | // if(emc->GetPHOSMod() == cpv->GetPHOSMod()){ |
757 | // // Correct to difference in CPV and EMC position due to different distance to center. |
758 | // // we assume, that particle moves from center |
759 | // Float_t dCPV = geom->GetIPtoOuterCoverDistance(); |
760 | // Float_t dEMC = geom->GetIPtoCrystalSurface() ; |
761 | // dEMC = dEMC / dCPV ; |
762 | // vecCpv = dEMC * vecCpv - vecEmc ; |
763 | // if (axis == "X") return vecCpv.X(); |
764 | // if (axis == "Y") return vecCpv.Y(); |
765 | // if (axis == "Z") return vecCpv.Z(); |
766 | // if (axis == "R") return vecCpv.Mag(); |
767 | // } |
768 | // return 100000000 ; |
769 | // } |
770 | // return 100000000 ; |
771 | //} |
12022e83 |
772 | //____________________________________________________________________________ |
26aa7e4a |
773 | Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const |
bc0c084c |
774 | { |
c947e71a |
775 | //Calculates the pid bit for the CPV selection per each purity. |
e3817e5f |
776 | if(effPur>2 || effPur<0) |
351dd634 |
777 | AliError(Form("Invalid Efficiency-Purity choice %d",effPur)); |
26aa7e4a |
778 | |
407d15b3 |
779 | //DP if(ts->GetCpvIndex()<0) |
780 | //DP return 1 ; //no CPV cluster |
bc0c084c |
781 | |
e3817e5f |
782 | Float_t sigX = GetCpv2EmcDistanceCut("X",e); |
783 | Float_t sigZ = GetCpv2EmcDistanceCut("Z",e); |
bc0c084c |
784 | |
26aa7e4a |
785 | Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X")); |
786 | Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z")); |
407d15b3 |
787 | // Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ; |
7fb46731 |
788 | |
789 | //if(deltaX>sigX*(effPur+1)) |
790 | //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1))) |
791 | if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1))) |
bc0c084c |
792 | return 1;//Neutral |
793 | else |
794 | return 0;//Charged |
bc0c084c |
795 | } |
69183710 |
796 | |
6ad0bfa0 |
797 | //____________________________________________________________________________ |
fc7e2f43 |
798 | Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const |
148b2bba |
799 | { |
50739f15 |
800 | //Is the particle inside de PCA ellipse? |
581354c5 |
801 | |
e3817e5f |
802 | particle.ToLower(); |
803 | Int_t prinbit = 0 ; |
7fb46731 |
804 | Float_t a = GetEllipseParameter(particle,"a" , e); |
805 | Float_t b = GetEllipseParameter(particle,"b" , e); |
806 | Float_t c = GetEllipseParameter(particle,"c" , e); |
e3817e5f |
807 | Float_t x0 = GetEllipseParameter(particle,"x0", e); |
808 | Float_t y0 = GetEllipseParameter(particle,"y0", e); |
809 | |
810 | Float_t r = TMath::Power((p[0] - x0)/a,2) + |
811 | TMath::Power((p[1] - y0)/b,2) + |
812 | c*(p[0] - x0)*(p[1] - y0)/(a*b) ; |
50739f15 |
813 | //3 different ellipses defined |
e3817e5f |
814 | if((effPur==2) && (r<1./2.)) prinbit= 1; |
815 | if((effPur==1) && (r<2. )) prinbit= 1; |
816 | if((effPur==0) && (r<9./2.)) prinbit= 1; |
50739f15 |
817 | |
581354c5 |
818 | if(r<0) |
351dd634 |
819 | AliError("Negative square?") ; |
1f0e7ccd |
820 | |
821 | return prinbit; |
148b2bba |
822 | |
148b2bba |
823 | } |
1f0e7ccd |
824 | //____________________________________________________________________________ |
fc7e2f43 |
825 | Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const |
1f0e7ccd |
826 | { |
e3817e5f |
827 | // Set bit for identified hard photons (E > 30 GeV) |
828 | // if the second moment M2x is below the boundary |
829 | |
830 | Float_t e = emc->GetEnergy(); |
831 | if (e < 30.0) return 0; |
832 | Float_t m2x = emc->GetM2x(); |
833 | Float_t m2xBoundary = GetParameterPhotonBoundary(0) * |
834 | TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/ |
835 | TMath::Power(GetParameterPhotonBoundary(2),2)) + |
836 | GetParameterPhotonBoundary(3); |
a675b8d6 |
837 | AliDebug(1, Form("E=%f, m2x=%f, boundary=%f", e,m2x,m2xBoundary)); |
e3817e5f |
838 | if (m2x < m2xBoundary) |
839 | return 1;// A hard photon |
840 | else |
841 | return 0;// Not a hard photon |
1f0e7ccd |
842 | } |
92f521a9 |
843 | |
e3817e5f |
844 | //____________________________________________________________________________ |
fc7e2f43 |
845 | Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const |
e3817e5f |
846 | { |
847 | // Set bit for identified hard pi0 (E > 30 GeV) |
848 | // if the second moment M2x is above the boundary |
849 | |
850 | Float_t e = emc->GetEnergy(); |
851 | if (e < 30.0) return 0; |
852 | Float_t m2x = emc->GetM2x(); |
853 | Float_t m2xBoundary = GetParameterPi0Boundary(0) + |
854 | e * GetParameterPi0Boundary(1); |
351dd634 |
855 | AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary)); |
e3817e5f |
856 | if (m2x > m2xBoundary) |
857 | return 1;// A hard pi0 |
bc0c084c |
858 | else |
e3817e5f |
859 | return 0;// Not a hard pi0 |
f0a4c9e9 |
860 | } |
e3817e5f |
861 | |
862 | //____________________________________________________________________________ |
8d4608b5 |
863 | TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const |
88cb7938 |
864 | { |
865 | // Calculates the momentum direction: |
866 | // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint |
867 | // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints |
868 | // However because of the poor position resolution of PPSD the direction is always taken as if we were |
869 | // in case 1. |
f0a4c9e9 |
870 | |
407d15b3 |
871 | TVector3 local ; |
872 | emc->GetLocalPosition(local) ; |
873 | |
874 | AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; |
875 | //Correct for the non-perpendicular incidence |
876 | // Correction for the depth of the shower starting point (TDR p 127) |
877 | Float_t para = 0.925 ; |
878 | Float_t parb = 6.52 ; |
879 | |
880 | //Remove Old correction (vertex at 0,0,0) |
881 | TVector3 vtxOld(0.,0.,0.) ; |
882 | TVector3 vInc ; |
883 | Float_t x=local.X() ; |
884 | Float_t z=local.Z() ; |
885 | phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ; |
886 | Float_t depthxOld = 0.; |
887 | Float_t depthzOld = 0.; |
888 | Float_t energy = emc->GetEnergy() ; |
889 | if (energy > 0 && vInc.Y()!=0.) { |
753b19cd |
890 | depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; |
891 | depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; |
7fb46731 |
892 | } |
adcca1e6 |
893 | else{ |
407d15b3 |
894 | AliError("Cluster with zero energy \n"); |
2c06dc7a |
895 | } |
407d15b3 |
896 | //Apply Real vertex |
897 | phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ; |
898 | Float_t depthx = 0.; |
899 | Float_t depthz = 0.; |
900 | if (energy > 0 && vInc.Y()!=0.) { |
e5b7b511 |
901 | depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; |
902 | depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; |
adcca1e6 |
903 | } |
407d15b3 |
904 | |
e5b7b511 |
905 | //Correct for the vertex position and shower depth |
906 | Double_t xd=x+(depthxOld-depthx) ; |
907 | Double_t zd=z+(depthzOld-depthz) ; |
908 | TVector3 dir(0,0,0) ; |
909 | phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ; |
407d15b3 |
910 | |
e5b7b511 |
911 | dir-=fVtx ; |
7fb46731 |
912 | dir.SetMag(1.) ; |
e3817e5f |
913 | |
88cb7938 |
914 | return dir ; |
7b7c1533 |
915 | } |
916 | |
35adb638 |
917 | //________________________________________________________________________ |
17323043 |
918 | Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par) |
35adb638 |
919 | { |
c947e71a |
920 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
921 | //this method returns a density probability of this parameter, given by a landau |
922 | //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b |
923 | |
2924941c |
924 | if (x > par[9]) x = par[9]; |
925 | |
926 | //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ; |
927 | Double_t cnt = par[0] + par[1] * x + par[2] * x * x ; |
35adb638 |
928 | Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ; |
929 | Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ; |
c947e71a |
930 | |
931 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 |
932 | return cnt*TMath::Landau(y,mean,sigma); |
35adb638 |
933 | } |
934 | else |
935 | return 0.; |
936 | |
937 | } |
938 | //________________________________________________________________________ |
17323043 |
939 | Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par) |
35adb638 |
940 | { |
c947e71a |
941 | |
942 | //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance), |
943 | //this method returns a density probability of this parameter, given by a landau |
944 | //function whose parameters depend with the energy like second order polinomial |
945 | |
35adb638 |
946 | Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ; |
c947e71a |
947 | Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ; |
948 | Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ; |
35adb638 |
949 | |
c947e71a |
950 | if(TMath::Abs(sigma) > 1.e-10){ |
acb5beb7 |
951 | return cnt*TMath::Landau(y,mean,sigma); |
35adb638 |
952 | } |
953 | else |
954 | return 0.; |
c947e71a |
955 | |
956 | |
35adb638 |
957 | } |
958 | // //________________________________________________________________________ |
959 | // Double_t AliPHOSPIDv1::ChargedHadronDistProb(Double_t x, Double_t y, Double_t * parg, Double_t * parl) |
960 | // { |
961 | // Double_t cnt = 0.0 ; |
962 | // Double_t mean = 0.0 ; |
963 | // Double_t sigma = 0.0 ; |
964 | // Double_t arg = 0.0 ; |
965 | // if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){ |
966 | // cnt = parg[1] / (x*x) + parg[2] / x + parg[0] ; |
967 | // mean = parg[4] / (x*x) + parg[5] / x + parg[3] ; |
968 | // sigma = parg[7] / (x*x) + parg[8] / x + parg[6] ; |
969 | // TF1 * f = new TF1("gaus","gaus",0.,100.); |
970 | // f->SetParameters(cnt,mean,sigma); |
971 | // arg = f->Eval(y) ; |
972 | // } |
973 | // else{ |
974 | // cnt = parl[1] / (x*x) + parl[2] / x + parl[0] ; |
975 | // mean = parl[4] / (x*x) + parl[5] / x + parl[3] ; |
976 | // sigma = parl[7] / (x*x) + parl[8] / x + parl[6] ; |
977 | // TF1 * f = new TF1("landau","landau",0.,100.); |
978 | // f->SetParameters(cnt,mean,sigma); |
979 | // arg = f->Eval(y) ; |
980 | // } |
981 | // // Double_t mean = par[3] + par[4] * x + par[5] * x * x ; |
982 | // // Double_t sigma = par[6] + par[7] * x + par[8] * x * x ; |
983 | |
984 | // //Double_t arg = -(y-mean)*(y-mean)/(2*sigma*sigma) ; |
985 | // //return cnt * TMath::Exp(arg) ; |
986 | |
987 | // return arg; |
988 | |
989 | // } |
2cc71c1e |
990 | //____________________________________________________________________________ |
991 | void AliPHOSPIDv1::MakePID() |
992 | { |
993 | // construct the PID weight from a Bayesian Method |
c947e71a |
994 | |
304864ab |
995 | const Int_t kSPECIES = AliPID::kSPECIESN ; |
7fb46731 |
996 | |
9a2cdbdf |
997 | Int_t nparticles = fRecParticles->GetEntriesFast() ; |
7fb46731 |
998 | |
9a2cdbdf |
999 | if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) { |
7fb46731 |
1000 | AliFatal("RecPoints or TrackSegments not found !") ; |
1001 | } |
9a2cdbdf |
1002 | |
1003 | TIter next(fTrackSegments) ; |
7fb46731 |
1004 | AliPHOSTrackSegment * ts ; |
1005 | Int_t index = 0 ; |
1006 | |
35adb638 |
1007 | Double_t * stof[kSPECIES] ; |
1008 | Double_t * sdp [kSPECIES] ; |
1009 | Double_t * scpv[kSPECIES] ; |
fb7b51ad |
1010 | Double_t * sw [kSPECIES] ; |
35adb638 |
1011 | //Info("MakePID","Begin MakePID"); |
1012 | |
1013 | for (Int_t i =0; i< kSPECIES; i++){ |
1014 | stof[i] = new Double_t[nparticles] ; |
1015 | sdp [i] = new Double_t[nparticles] ; |
1016 | scpv[i] = new Double_t[nparticles] ; |
fb7b51ad |
1017 | sw [i] = new Double_t[nparticles] ; |
35adb638 |
1018 | } |
1019 | |
7fb46731 |
1020 | |
1021 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
1022 | |
1023 | //cout<<">>>>>> Bayesian Index "<<index<<endl; |
1024 | |
1025 | AliPHOSEmcRecPoint * emc = 0 ; |
1026 | if(ts->GetEmcIndex()>=0) |
9a2cdbdf |
1027 | emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ; |
7fb46731 |
1028 | |
407d15b3 |
1029 | // AliPHOSCpvRecPoint * cpv = 0 ; |
1030 | // if(ts->GetCpvIndex()>=0) |
1031 | // cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ; |
1032 | // |
1033 | //// Int_t track = 0 ; |
1034 | //// track = ts->GetTrackIndex() ; //TPC tracks ? |
cc1fe362 |
1035 | |
7fb46731 |
1036 | if (!emc) { |
a675b8d6 |
1037 | AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ; |
7fb46731 |
1038 | } |
c947e71a |
1039 | |
2924941c |
1040 | |
7fb46731 |
1041 | // ############Tof############################# |
c947e71a |
1042 | |
7fb46731 |
1043 | // Info("MakePID", "TOF"); |
fb7b51ad |
1044 | Float_t en = emc->GetEnergy(); |
7fb46731 |
1045 | Double_t time = emc->GetTime() ; |
1046 | // cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl; |
fb7b51ad |
1047 | |
35adb638 |
1048 | // now get the signals probability |
1049 | // s(pid) in the Bayesian formulation |
cc1fe362 |
1050 | |
304864ab |
1051 | stof[AliPID::kPhoton][index] = 1.; |
1052 | stof[AliPID::kElectron][index] = 1.; |
304864ab |
1053 | stof[AliPID::kEleCon][index] = 1.; |
7fb46731 |
1054 | //We assing the same prob to charged hadrons, sum is 1 |
1055 | stof[AliPID::kPion][index] = 1./3.; |
1056 | stof[AliPID::kKaon][index] = 1./3.; |
1057 | stof[AliPID::kProton][index] = 1./3.; |
1058 | //We assing the same prob to neutral hadrons, sum is 1 |
1059 | stof[AliPID::kNeutron][index] = 1./2.; |
1060 | stof[AliPID::kKaon0][index] = 1./2.; |
304864ab |
1061 | stof[AliPID::kMuon][index] = 1.; |
2924941c |
1062 | |
1063 | if(en < fTOFEnThreshold) { |
7fb46731 |
1064 | |
1065 | Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution |
1066 | Double_t pTofKaon = 0; |
1067 | |
35adb638 |
1068 | if(time < fTkaonl[1]) |
fb7b51ad |
1069 | pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution |
35adb638 |
1070 | else |
fb7b51ad |
1071 | pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution |
7fb46731 |
1072 | |
1073 | Double_t pTofNucleon = 0; |
1074 | |
35adb638 |
1075 | if(time < fThhadronl[1]) |
7fb46731 |
1076 | pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution |
35adb638 |
1077 | else |
7fb46731 |
1078 | pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution |
7fb46731 |
1079 | //We assing the same prob to neutral hadrons, sum is the average prob |
fb7b51ad |
1080 | Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ; |
1081 | //We assing the same prob to charged hadrons, sum is the average prob |
1082 | Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ; |
7fb46731 |
1083 | |
fb7b51ad |
1084 | stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ; |
1085 | //gaus distribution |
1086 | stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ; |
1087 | //a conversion electron has the photon ToF |
304864ab |
1088 | stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ; |
7fb46731 |
1089 | |
fb7b51ad |
1090 | stof[AliPID::kElectron][index] = pTofPion ; |
7fb46731 |
1091 | |
1092 | stof[AliPID::kPion][index] = pTofChHadron ; |
1093 | stof[AliPID::kKaon][index] = pTofChHadron ; |
1094 | stof[AliPID::kProton][index] = pTofChHadron ; |
1095 | |
1096 | stof[AliPID::kKaon0][index] = pTofNeHadron ; |
1097 | stof[AliPID::kNeutron][index] = pTofNeHadron ; |
cc1fe362 |
1098 | } |
c947e71a |
1099 | |
1100 | // Info("MakePID", "Dispersion"); |
cc1fe362 |
1101 | |
7fb46731 |
1102 | // ###########Shower shape: Dispersion#################### |
35adb638 |
1103 | Float_t dispersion = emc->GetDispersion(); |
407d15b3 |
1104 | //DP: Correct for non-perpendicular incidence |
1105 | //DP: still to be done |
1106 | |
35adb638 |
1107 | //dispersion is not well defined if the cluster is only in few crystals |
cc1fe362 |
1108 | |
304864ab |
1109 | sdp[AliPID::kPhoton][index] = 1. ; |
1110 | sdp[AliPID::kElectron][index] = 1. ; |
1111 | sdp[AliPID::kPion][index] = 1. ; |
1112 | sdp[AliPID::kKaon][index] = 1. ; |
1113 | sdp[AliPID::kProton][index] = 1. ; |
1114 | sdp[AliPID::kNeutron][index] = 1. ; |
1115 | sdp[AliPID::kEleCon][index] = 1. ; |
1116 | sdp[AliPID::kKaon0][index] = 1. ; |
1117 | sdp[AliPID::kMuon][index] = 1. ; |
c947e71a |
1118 | |
2924941c |
1119 | if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){ |
7fb46731 |
1120 | sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ; |
304864ab |
1121 | sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ; |
1122 | sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ; |
1123 | sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ; |
1124 | sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ; |
1125 | sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ; |
1126 | sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index]; |
1127 | sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ; |
fb7b51ad |
1128 | sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ; |
1129 | //landau distribution |
35adb638 |
1130 | } |
cc1fe362 |
1131 | |
7fb46731 |
1132 | // Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion); |
1133 | // Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]); |
c947e71a |
1134 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ; |
304864ab |
1135 | // cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl; |
c947e71a |
1136 | |
7fb46731 |
1137 | //########## CPV-EMC Distance####################### |
1138 | // Info("MakePID", "Distance"); |
fb7b51ad |
1139 | |
26aa7e4a |
1140 | Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ; |
1141 | Float_t z = ts->GetCpvDistance("Z") ; |
fb7b51ad |
1142 | |
7fb46731 |
1143 | Double_t pcpv = 0 ; |
c947e71a |
1144 | Double_t pcpvneutral = 0. ; |
7fb46731 |
1145 | |
1146 | Double_t elprobx = GausF(en , x, fXelectron) ; |
1147 | Double_t elprobz = GausF(en , z, fZelectron) ; |
1148 | Double_t chprobx = GausF(en , x, fXcharged) ; |
1149 | Double_t chprobz = GausF(en , z, fZcharged) ; |
c947e71a |
1150 | Double_t pcpvelectron = elprobx * elprobz; |
1151 | Double_t pcpvcharged = chprobx * chprobz; |
7fb46731 |
1152 | |
1153 | // cout<<">>>>energy "<<en<<endl; |
c947e71a |
1154 | // cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl; |
1155 | // cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl; |
1156 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; |
1157 | |
2924941c |
1158 | // Is neutral or charged? |
35adb638 |
1159 | if(pcpvelectron >= pcpvcharged) |
1160 | pcpv = pcpvelectron ; |
1161 | else |
1162 | pcpv = pcpvcharged ; |
1163 | |
fb7b51ad |
1164 | if(pcpv < fChargedNeutralThreshold) |
35adb638 |
1165 | { |
1166 | pcpvneutral = 1. ; |
1167 | pcpvcharged = 0. ; |
1168 | pcpvelectron = 0. ; |
1169 | } |
c947e71a |
1170 | // else |
1171 | // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl; |
35adb638 |
1172 | |
304864ab |
1173 | scpv[AliPID::kPion][index] = pcpvcharged ; |
1174 | scpv[AliPID::kKaon][index] = pcpvcharged ; |
1175 | scpv[AliPID::kProton][index] = pcpvcharged ; |
7fb46731 |
1176 | |
1177 | scpv[AliPID::kMuon][index] = pcpvelectron ; |
304864ab |
1178 | scpv[AliPID::kElectron][index] = pcpvelectron ; |
304864ab |
1179 | scpv[AliPID::kEleCon][index] = pcpvelectron ; |
7fb46731 |
1180 | |
1181 | scpv[AliPID::kPhoton][index] = pcpvneutral ; |
1182 | scpv[AliPID::kNeutron][index] = pcpvneutral ; |
304864ab |
1183 | scpv[AliPID::kKaon0][index] = pcpvneutral ; |
7fb46731 |
1184 | |
35adb638 |
1185 | |
1186 | // Info("MakePID", "CPV passed"); |
c947e71a |
1187 | |
7fb46731 |
1188 | //############## Pi0 ############################# |
304864ab |
1189 | stof[AliPID::kPi0][index] = 0. ; |
1190 | scpv[AliPID::kPi0][index] = 0. ; |
1191 | sdp [AliPID::kPi0][index] = 0. ; |
c947e71a |
1192 | |
35adb638 |
1193 | if(en > 30.){ |
1194 | // pi0 are detected via decay photon |
2924941c |
1195 | stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index]; |
304864ab |
1196 | scpv[AliPID::kPi0][index] = pcpvneutral ; |
2924941c |
1197 | if(emc->GetMultiplicity() > fDispMultThreshold) |
1198 | sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ; |
1199 | //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ; |
1200 | // cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = " |
1201 | // <<emc->GetMultiplicity()<<endl; |
1202 | // cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = " |
1203 | // <<sdp [AliPID::kPi0][index]<<endl; |
35adb638 |
1204 | } |
1205 | |
2924941c |
1206 | |
7fb46731 |
1207 | |
2924941c |
1208 | |
7fb46731 |
1209 | //############## muon ############################# |
1210 | |
35adb638 |
1211 | if(en > 0.5){ |
1212 | //Muons deposit few energy |
304864ab |
1213 | scpv[AliPID::kMuon][index] = 0 ; |
1214 | stof[AliPID::kMuon][index] = 0 ; |
1215 | sdp [AliPID::kMuon][index] = 0 ; |
c947e71a |
1216 | } |
1217 | |
fb7b51ad |
1218 | //Weight to apply to hadrons due to energy reconstruction |
1219 | |
1220 | Float_t weight = fERecWeight ->Eval(en) ; |
1221 | |
1222 | sw[AliPID::kPhoton][index] = 1. ; |
1223 | sw[AliPID::kElectron][index] = 1. ; |
1224 | sw[AliPID::kPion][index] = weight ; |
1225 | sw[AliPID::kKaon][index] = weight ; |
1226 | sw[AliPID::kProton][index] = weight ; |
1227 | sw[AliPID::kNeutron][index] = weight ; |
1228 | sw[AliPID::kEleCon][index] = 1. ; |
1229 | sw[AliPID::kKaon0][index] = weight ; |
1230 | sw[AliPID::kMuon][index] = weight ; |
1231 | sw[AliPID::kPi0][index] = 1. ; |
1232 | |
7fb46731 |
1233 | // if(en > 0.5){ |
1234 | // cout<<"######################################################"<<endl; |
1235 | // //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ; |
1236 | // cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ; |
1237 | // cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl; |
1238 | // cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl; |
1239 | // cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl; |
1240 | // cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl; |
c947e71a |
1241 | |
2924941c |
1242 | // cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index] |
1243 | // <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl; |
7fb46731 |
1244 | // cout<<"EleCon , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index] |
1245 | // <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl; |
1246 | // cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index] |
1247 | // <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl; |
1248 | // cout<<"Muon , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index] |
1249 | // <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl; |
2924941c |
1250 | // cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index] |
1251 | // <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl; |
7fb46731 |
1252 | // cout<<"Pion , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index] |
1253 | // <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl; |
1254 | // cout<<"Kaon0 , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index] |
1255 | // <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl; |
1256 | // cout<<"Kaon , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index] |
1257 | // <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl; |
1258 | // cout<<"Neutron , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index] |
1259 | // <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl; |
1260 | // cout<<"Proton , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index] |
1261 | // <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl; |
1262 | // cout<<"######################################################"<<endl; |
1263 | // } |
1264 | index++; |
cc1fe362 |
1265 | } |
35adb638 |
1266 | |
1267 | //for (index = 0 ; index < kSPECIES ; index++) |
1268 | // pid[index] /= nparticles ; |
1269 | |
7fb46731 |
1270 | |
35adb638 |
1271 | // Info("MakePID", "Total Probability calculation"); |
1272 | |
cc1fe362 |
1273 | for(index = 0 ; index < nparticles ; index ++) { |
2924941c |
1274 | |
9a2cdbdf |
1275 | AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index)); |
2924941c |
1276 | |
1277 | //Conversion electron? |
1278 | |
1279 | if(recpar->IsEleCon()){ |
1280 | fInitPID[AliPID::kEleCon] = 1. ; |
1281 | fInitPID[AliPID::kPhoton] = 0. ; |
1282 | fInitPID[AliPID::kElectron] = 0. ; |
1283 | } |
1284 | else{ |
1285 | fInitPID[AliPID::kEleCon] = 0. ; |
1286 | fInitPID[AliPID::kPhoton] = 1. ; |
1287 | fInitPID[AliPID::kElectron] = 1. ; |
1288 | } |
1289 | // fInitPID[AliPID::kEleCon] = 0. ; |
1290 | |
1291 | |
cc1fe362 |
1292 | // calculates the Bayesian weight |
7fb46731 |
1293 | |
cc1fe362 |
1294 | Int_t jndex ; |
1295 | Double_t wn = 0.0 ; |
1296 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) |
2924941c |
1297 | wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * |
1298 | sw[jndex][index] * fInitPID[jndex] ; |
1299 | |
7fb46731 |
1300 | // cout<<"*************wn "<<wn<<endl; |
e74ea0e9 |
1301 | if (TMath::Abs(wn)>0) |
1302 | for (jndex = 0 ; jndex < kSPECIES ; jndex++) { |
35adb638 |
1303 | //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn" |
1304 | //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex] << endl; |
1305 | //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl; |
2924941c |
1306 | // if(jndex == AliPID::kPi0 || jndex == AliPID::kPhoton){ |
1307 | // cout<<"Particle "<<jndex<<" final prob * wn " |
1308 | // <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] * |
1309 | // fInitPID[jndex] <<" wn "<< wn<<endl; |
1310 | // cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index] |
1311 | // <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl; |
1312 | // } |
1313 | recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] * |
1314 | sw[jndex][index] * scpv[jndex][index] * |
1315 | fInitPID[jndex] / wn) ; |
e74ea0e9 |
1316 | } |
2cc71c1e |
1317 | } |
35adb638 |
1318 | // Info("MakePID", "Delete"); |
1319 | |
2924941c |
1320 | for (Int_t i =0; i< kSPECIES; i++){ |
1321 | delete [] stof[i]; |
1322 | delete [] sdp [i]; |
1323 | delete [] scpv[i]; |
1324 | delete [] sw [i]; |
1325 | } |
35adb638 |
1326 | // Info("MakePID","End MakePID"); |
2cc71c1e |
1327 | } |
1328 | |
7acf6008 |
1329 | //____________________________________________________________________________ |
e3817e5f |
1330 | void AliPHOSPIDv1::MakeRecParticles() |
1331 | { |
b2a60966 |
1332 | // Makes a RecParticle out of a TrackSegment |
148b2bba |
1333 | |
9a2cdbdf |
1334 | if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) { |
351dd634 |
1335 | AliFatal("RecPoints or TrackSegments not found !") ; |
148b2bba |
1336 | } |
9a2cdbdf |
1337 | fRecParticles->Clear(); |
148b2bba |
1338 | |
9a2cdbdf |
1339 | TIter next(fTrackSegments) ; |
7acf6008 |
1340 | AliPHOSTrackSegment * ts ; |
6ad0bfa0 |
1341 | Int_t index = 0 ; |
09fc14a0 |
1342 | AliPHOSRecParticle * rp ; |
7acf6008 |
1343 | while ( (ts = (AliPHOSTrackSegment *)next()) ) { |
7fb46731 |
1344 | // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl; |
9a2cdbdf |
1345 | new( (*fRecParticles)[index] ) AliPHOSRecParticle() ; |
1346 | rp = (AliPHOSRecParticle *)fRecParticles->At(index) ; |
f0a4c9e9 |
1347 | rp->SetTrackSegment(index) ; |
9688c1dd |
1348 | rp->SetIndexInList(index) ; |
148b2bba |
1349 | |
7acf6008 |
1350 | AliPHOSEmcRecPoint * emc = 0 ; |
1351 | if(ts->GetEmcIndex()>=0) |
9a2cdbdf |
1352 | emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ; |
fad3e5b9 |
1353 | |
8d4608b5 |
1354 | AliPHOSCpvRecPoint * cpv = 0 ; |
7acf6008 |
1355 | if(ts->GetCpvIndex()>=0) |
9a2cdbdf |
1356 | cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ; |
fad3e5b9 |
1357 | |
bd76890a |
1358 | Int_t track = 0 ; |
1359 | track = ts->GetTrackIndex() ; |
1360 | |
148b2bba |
1361 | // Now set type (reconstructed) of the particle |
1362 | |
1363 | // Choose the cluster energy range |
9fa5f1d0 |
1364 | |
fbf811ec |
1365 | if (!emc) { |
a675b8d6 |
1366 | AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ; |
fbf811ec |
1367 | } |
50739f15 |
1368 | |
e3817e5f |
1369 | Float_t e = emc->GetEnergy() ; |
bc0c084c |
1370 | |
6f969528 |
1371 | Float_t lambda[2] ; |
1372 | emc->GetElipsAxis(lambda) ; |
407d15b3 |
1373 | |
50739f15 |
1374 | if((lambda[0]>0.01) && (lambda[1]>0.01)){ |
1375 | // Looking PCA. Define and calculate the data (X), |
bc0c084c |
1376 | // introduce in the function X2P that gives the components (P). |
1377 | |
c947e71a |
1378 | Float_t spher = 0. ; |
1379 | Float_t emaxdtotal = 0. ; |
50739f15 |
1380 | |
bc0c084c |
1381 | if((lambda[0]+lambda[1])!=0) |
e4df4b30 |
1382 | spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); |
50739f15 |
1383 | |
c947e71a |
1384 | emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); |
50739f15 |
1385 | |
1386 | fX[0] = lambda[0] ; |
1387 | fX[1] = lambda[1] ; |
1388 | fX[2] = emc->GetDispersion() ; |
c947e71a |
1389 | fX[3] = spher ; |
50739f15 |
1390 | fX[4] = emc->GetMultiplicity() ; |
c947e71a |
1391 | fX[5] = emaxdtotal ; |
50739f15 |
1392 | fX[6] = emc->GetCoreEnergy() ; |
1393 | |
e3817e5f |
1394 | fPrincipalPhoton->X2P(fX,fPPhoton); |
1395 | fPrincipalPi0 ->X2P(fX,fPPi0); |
1f0e7ccd |
1396 | |
50739f15 |
1397 | } |
1398 | else{ |
e3817e5f |
1399 | fPPhoton[0]=-100.0; //We do not accept clusters with |
1400 | fPPhoton[1]=-100.0; //one cell as a photon-like |
1401 | fPPi0[0] =-100.0; |
1402 | fPPi0[1] =-100.0; |
50739f15 |
1403 | } |
1404 | |
2cc71c1e |
1405 | Float_t time = emc->GetTime() ; |
1406 | rp->SetTof(time) ; |
9fa5f1d0 |
1407 | |
bc0c084c |
1408 | // Loop of Efficiency-Purity (the 3 points of purity or efficiency |
1409 | // are taken into account to set the particle identification) |
e3817e5f |
1410 | for(Int_t effPur = 0; effPur < 3 ; effPur++){ |
50739f15 |
1411 | |
bc0c084c |
1412 | // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, |
1413 | // 1st,2nd or 3rd bit (depending on the efficiency-purity point ) |
1414 | // is set to 1 |
26aa7e4a |
1415 | if(GetCPVBit(ts, effPur,e) == 1 ){ |
e3817e5f |
1416 | rp->SetPIDBit(effPur) ; |
35adb638 |
1417 | //cout<<"CPV bit "<<effPur<<endl; |
1418 | } |
50739f15 |
1419 | // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th |
1420 | // bit (depending on the efficiency-purity point )is set to 1 |
2cc71c1e |
1421 | if(time< (*fParameters)(3,effPur)) |
e3817e5f |
1422 | rp->SetPIDBit(effPur+3) ; |
2cc71c1e |
1423 | |
e3817e5f |
1424 | //Photon PCA |
50739f15 |
1425 | //If we are inside the ellipse, 7th, 8th or 9th |
1426 | // bit (depending on the efficiency-purity point )is set to 1 |
e3817e5f |
1427 | if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) |
1428 | rp->SetPIDBit(effPur+6) ; |
1f0e7ccd |
1429 | |
e3817e5f |
1430 | //Pi0 PCA |
1f0e7ccd |
1431 | //If we are inside the ellipse, 10th, 11th or 12th |
1432 | // bit (depending on the efficiency-purity point )is set to 1 |
e3817e5f |
1433 | if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1) |
1434 | rp->SetPIDBit(effPur+9) ; |
f0a4c9e9 |
1435 | } |
e3817e5f |
1436 | if(GetHardPhotonBit(emc)) |
1437 | rp->SetPIDBit(12) ; |
1438 | if(GetHardPi0Bit (emc)) |
1439 | rp->SetPIDBit(13) ; |
1f0e7ccd |
1440 | |
bd76890a |
1441 | if(track >= 0) |
1442 | rp->SetPIDBit(14) ; |
1443 | |
9fa5f1d0 |
1444 | //Set momentum, energy and other parameters |
50739f15 |
1445 | Float_t encal = GetCalibratedEnergy(e); |
9fa5f1d0 |
1446 | TVector3 dir = GetMomentumDirection(emc,cpv) ; |
1447 | dir.SetMag(encal) ; |
1448 | rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ; |
1449 | rp->SetCalcMass(0); |
e0ed2e49 |
1450 | rp->Name(); //If photon sets the particle pdg name to gamma |
407d15b3 |
1451 | rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0); |
e747b8da |
1452 | rp->SetFirstMother(-1); |
1453 | rp->SetLastMother(-1); |
1454 | rp->SetFirstDaughter(-1); |
1455 | rp->SetLastDaughter(-1); |
1456 | rp->SetPolarisation(0,0,0); |
d956e9b7 |
1457 | //Set the position in global coordinate system from the RecPoint |
605cb8bb |
1458 | AliPHOSTrackSegment * ts1 = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex())); |
1459 | AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex())); |
d956e9b7 |
1460 | TVector3 pos ; |
9a2cdbdf |
1461 | fGeom->GetGlobalPHOS(erp, pos) ; |
d956e9b7 |
1462 | rp->SetPos(pos); |
6ad0bfa0 |
1463 | index++ ; |
1464 | } |
6ad0bfa0 |
1465 | } |
e3817e5f |
1466 | |
09fc14a0 |
1467 | //____________________________________________________________________________ |
702ab87e |
1468 | void AliPHOSPIDv1::Print(const Option_t *) const |
09fc14a0 |
1469 | { |
b2a60966 |
1470 | // Print the parameters used for the particle type identification |
bc0c084c |
1471 | |
351dd634 |
1472 | AliInfo("=============== AliPHOSPIDv1 ================") ; |
88cb7938 |
1473 | printf("Making PID\n") ; |
1474 | printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ; |
1475 | printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ; |
1476 | printf(" Matrix of Parameters: 14x4\n") ; |
1477 | printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ; |
1478 | printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ; |
1479 | printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ; |
1480 | printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ; |
407d15b3 |
1481 | printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ; |
50739f15 |
1482 | fParameters->Print() ; |
09fc14a0 |
1483 | } |
1484 | |
a496c46c |
1485 | |
69183710 |
1486 | |
7acf6008 |
1487 | //____________________________________________________________________________ |
a4e98857 |
1488 | void AliPHOSPIDv1::PrintRecParticles(Option_t * option) |
1489 | { |
dd5c4038 |
1490 | // Print table of reconstructed particles |
1491 | |
21cd0c07 |
1492 | TString message ; |
9a2cdbdf |
1493 | message = " found " ; |
1494 | message += fRecParticles->GetEntriesFast(); |
3bf72d32 |
1495 | message += " RecParticles\n" ; |
1496 | |
7acf6008 |
1497 | if(strstr(option,"all")) { // printing found TS |
3bf72d32 |
1498 | message += "\n PARTICLE Index \n" ; |
7acf6008 |
1499 | |
1500 | Int_t index ; |
9a2cdbdf |
1501 | for (index = 0 ; index < fRecParticles->GetEntries() ; index++) { |
1502 | AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ; |
3bf72d32 |
1503 | message += "\n" ; |
1504 | message += rp->Name().Data() ; |
1505 | message += " " ; |
1506 | message += rp->GetIndexInList() ; |
1507 | message += " " ; |
1508 | message += rp->GetType() ; |
7acf6008 |
1509 | } |
3bf72d32 |
1510 | } |
351dd634 |
1511 | AliInfo(message.Data() ) ; |
69183710 |
1512 | } |
88cb7938 |
1513 | |
1514 | //____________________________________________________________________________ |
1515 | void AliPHOSPIDv1::SetParameters() |
1516 | { |
1517 | // PCA : To do the Principal Components Analysis it is necessary |
1518 | // the Principal file, which is opened here |
1519 | fX = new double[7]; // Data for the PCA |
1520 | fPPhoton = new double[7]; // Eigenvalues of the PCA |
1521 | fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA |
1522 | |
1523 | // Read photon principals from the photon file |
1524 | |
1525 | fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; |
1526 | TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ; |
1527 | fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; |
1528 | f.Close() ; |
1529 | |
1530 | // Read pi0 principals from the pi0 file |
1531 | |
1532 | fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ; |
1533 | TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ; |
1534 | fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; |
1535 | fPi0.Close() ; |
1536 | |
1537 | // Open parameters file and initialization of the Parameters matrix. |
1538 | // In the File Parameters.dat are all the parameters. These are introduced |
1539 | // in a matrix of 16x4 |
1540 | // |
1541 | // All the parameters defined in this file are, in order of row: |
1542 | // line 0 : calibration |
1543 | // lines 1,2 : CPV rectangular cat for X and Z |
1544 | // line 3 : TOF cut |
1545 | // lines 4-8 : parameters to calculate photon PCA ellipse |
1546 | // lines 9-13: parameters to calculate pi0 PCA ellipse |
1547 | // lines 14-15: parameters to calculate border for high-pt photons and pi0 |
1548 | |
1549 | fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); |
e8d02863 |
1550 | fParameters = new TMatrixF(16,4) ; |
c947e71a |
1551 | const Int_t kMaxLeng=255; |
1552 | char string[kMaxLeng]; |
88cb7938 |
1553 | |
1554 | // Open a text file with PID parameters |
1555 | FILE *fd = fopen(fFileNameParameters.Data(),"r"); |
1556 | if (!fd) |
351dd634 |
1557 | AliFatal(Form("File %s with a PID parameters cannot be opened\n", |
1558 | fFileNameParameters.Data())); |
88cb7938 |
1559 | |
1560 | Int_t i=0; |
1561 | // Read parameter file line-by-line and skip empty line and comments |
c947e71a |
1562 | while (fgets(string,kMaxLeng,fd) != NULL) { |
88cb7938 |
1563 | if (string[0] == '\n' ) continue; |
1564 | if (string[0] == '!' ) continue; |
1565 | sscanf(string, "%f %f %f %f", |
1566 | &(*fParameters)(i,0), &(*fParameters)(i,1), |
1567 | &(*fParameters)(i,2), &(*fParameters)(i,3)); |
1568 | i++; |
a675b8d6 |
1569 | AliDebug(1, Form("Line %d: %s",i,string)); |
88cb7938 |
1570 | } |
1571 | fclose(fd); |
1572 | } |
1573 | |
1574 | //____________________________________________________________________________ |
1575 | void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) |
1576 | { |
1577 | // Set parameter "Calibration" i to a value param |
351dd634 |
1578 | if(i>2 || i<0) { |
1579 | AliError(Form("Invalid parameter number: %d",i)); |
1580 | } else |
88cb7938 |
1581 | (*fParameters)(0,i) = param ; |
1582 | } |
1583 | |
1584 | //____________________________________________________________________________ |
1585 | void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) |
1586 | { |
1587 | // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on |
1588 | // Purity-Efficiency point i |
1589 | |
351dd634 |
1590 | if(i>2 || i<0) { |
1591 | AliError(Form("Invalid parameter number: %d",i)); |
1592 | } else { |
88cb7938 |
1593 | axis.ToLower(); |
1594 | if (axis == "x") (*fParameters)(1,i) = cut; |
1595 | else if (axis == "z") (*fParameters)(2,i) = cut; |
351dd634 |
1596 | else { |
1597 | AliError(Form("Invalid axis name: %s",axis.Data())); |
1598 | } |
88cb7938 |
1599 | } |
1600 | } |
1601 | |
1602 | //____________________________________________________________________________ |
1603 | void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) |
1604 | { |
1605 | // Set parameter "Hard photon boundary" i to a value param |
351dd634 |
1606 | if(i>4 || i<0) { |
1607 | AliError(Form("Invalid parameter number: %d",i)); |
1608 | } else |
88cb7938 |
1609 | (*fParameters)(14,i) = param ; |
1610 | } |
1611 | |
1612 | //____________________________________________________________________________ |
1613 | void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) |
1614 | { |
1615 | // Set parameter "Hard pi0 boundary" i to a value param |
351dd634 |
1616 | if(i>1 || i<0) { |
1617 | AliError(Form("Invalid parameter number: %d",i)); |
1618 | } else |
88cb7938 |
1619 | (*fParameters)(15,i) = param ; |
1620 | } |
1621 | |
1622 | //_____________________________________________________________________________ |
1623 | void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) |
1624 | { |
1625 | // Set the parameter TimeGate depending on Purity-Efficiency point i |
351dd634 |
1626 | if (i>2 || i<0) { |
1627 | AliError(Form("Invalid Efficiency-Purity choice %d",i)); |
1628 | } else |
88cb7938 |
1629 | (*fParameters)(3,i)= gate ; |
1630 | } |
1631 | |
1632 | //_____________________________________________________________________________ |
1633 | void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) |
1634 | { |
1635 | // Set the parameter "i" that is needed to calculate the ellipse |
1636 | // parameter "param" for a particle "particle" |
1637 | |
1638 | particle.ToLower(); |
1639 | param. ToLower(); |
1640 | Int_t p= -1; |
1641 | Int_t offset=0; |
1642 | |
1643 | if (particle == "photon") offset=0; |
1644 | else if (particle == "pi0") offset=5; |
1645 | else |
351dd634 |
1646 | AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n", |
1647 | particle.Data())); |
88cb7938 |
1648 | |
1649 | if (param.Contains("a")) p=4+offset; |
1650 | else if(param.Contains("b")) p=5+offset; |
1651 | else if(param.Contains("c")) p=6+offset; |
1652 | else if(param.Contains("x0"))p=7+offset; |
1653 | else if(param.Contains("y0"))p=8+offset; |
351dd634 |
1654 | if((i>4)||(i<0)) { |
1655 | AliError(Form("No parameter with index %d", i)) ; |
1656 | } else if(p==-1) { |
1657 | AliError(Form("No parameter with name %s", param.Data() )) ; |
1658 | } else |
88cb7938 |
1659 | (*fParameters)(p,i) = par ; |
1660 | } |
1661 | |
407d15b3 |
1662 | //____________________________________________________________________________ |
1663 | void AliPHOSPIDv1::GetVertex(void) |
1664 | { //extract vertex either using ESD or generator |
1665 | |
1666 | //Try to extract vertex from data |
1667 | if(fESD){ |
1668 | const AliESDVertex *esdVtx = fESD->GetVertex() ; |
e5b7b511 |
1669 | if(esdVtx && esdVtx->GetChi2()!=0.){ |
407d15b3 |
1670 | fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ; |
1671 | return ; |
1672 | } |
1673 | } |
9a2cdbdf |
1674 | |
1675 | // Use vertex diamond from CDB GRP folder if the one from ESD is missing |
1676 | // PLEASE FIX IT |
407d15b3 |
1677 | AliWarning("Can not read vertex from data, use fixed \n") ; |
1678 | fVtx.SetXYZ(0.,0.,0.) ; |
1679 | |
1680 | } |
35adb638 |
1681 | //_______________________________________________________________________ |
1682 | void AliPHOSPIDv1::SetInitPID(const Double_t *p) { |
1683 | // Sets values for the initial population of each particle type |
304864ab |
1684 | for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i]; |
35adb638 |
1685 | } |
1686 | //_______________________________________________________________________ |
1687 | void AliPHOSPIDv1::GetInitPID(Double_t *p) const { |
1688 | // Gets values for the initial population of each particle type |
304864ab |
1689 | for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i]; |
35adb638 |
1690 | } |