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