]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PHOS/AliPHOSPIDv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
New class added.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSPIDv1.cxx
... / ...
CommitLineData
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
146ClassImp( AliPHOSPIDv1)
147
148//____________________________________________________________________________
149AliPHOSPIDv1::AliPHOSPIDv1() :
150 AliPHOSPID(),
151 fBayesian(kFALSE),
152 fDefaultInit(kFALSE),
153 fWrite(kFALSE),
154 fFileNamePrincipalPhoton(),
155 fFileNamePrincipalPi0(),
156 fFileNameParameters(),
157 fPrincipalPhoton(0),
158 fPrincipalPi0(0),
159 fX(0),
160 fPPhoton(0),
161 fPPi0(0),
162 fParameters(0),
163 fVtx(0.),
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)
176{
177 // default ctor
178
179 InitParameters() ;
180 fDefaultInit = kTRUE ;
181}
182
183//____________________________________________________________________________
184AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ) :
185 AliPHOSPID(pid),
186 fBayesian(kFALSE),
187 fDefaultInit(kFALSE),
188 fWrite(kFALSE),
189 fFileNamePrincipalPhoton(),
190 fFileNamePrincipalPi0(),
191 fFileNameParameters(),
192 fPrincipalPhoton(0),
193 fPrincipalPi0(0),
194 fX(0),
195 fPPhoton(0),
196 fPPi0(0),
197 fParameters(0),
198 fVtx(0.),
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
212{
213 // ctor
214 InitParameters() ;
215
216}
217
218//____________________________________________________________________________
219AliPHOSPIDv1::AliPHOSPIDv1(AliPHOSGeometry *geom):
220 AliPHOSPID(geom),
221 fBayesian(kFALSE),
222 fDefaultInit(kFALSE),
223 fWrite(kFALSE),
224 fFileNamePrincipalPhoton(),
225 fFileNamePrincipalPi0(),
226 fFileNameParameters(),
227 fPrincipalPhoton(0),
228 fPrincipalPi0(0),
229 fX(0),
230 fPPhoton(0),
231 fPPi0(0),
232 fParameters(0),
233 fVtx(0.),
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
247{
248 //ctor with the indication on where to look for the track segments
249
250 InitParameters() ;
251 fDefaultInit = kFALSE ;
252}
253
254//____________________________________________________________________________
255AliPHOSPIDv1::~AliPHOSPIDv1()
256{
257 // dtor
258 fPrincipalPhoton = 0;
259 fPrincipalPi0 = 0;
260
261 delete [] fX ; // Principal input
262 delete [] fPPhoton ; // Photon Principal components
263 delete [] fPPi0 ; // Pi0 Principal components
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;
273}
274
275//____________________________________________________________________________
276void AliPHOSPIDv1::InitParameters()
277{
278 // Initialize PID parameters
279 fWrite = kTRUE ;
280 fBayesian = kTRUE ;
281 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 fTphoton[0] = 7.83E8 ;
331 fTphoton[1] = 1.55E-8 ;
332 fTphoton[2] = 5.09E-10 ;
333 fTFphoton = new TFormula("ToF response to photons" , "gaus") ;
334 fTFphoton->SetParameters( fTphoton[0], fTphoton[1], fTphoton[2]) ;
335
336 // Pions
337 //Gaus (0 to max probability)
338 fTpiong[0] = 6.73E8 ;
339 fTpiong[1] = 1.58E-8 ;
340 fTpiong[2] = 5.87E-10 ;
341 fTFpiong = new TFormula("ToF response to pions" , "gaus") ;
342 fTFpiong->SetParameters( fTpiong[0], fTpiong[1], fTpiong[2]) ;
343
344 // Kaons
345 //Gaus (0 to max probability)
346 fTkaong[0] = 3.93E8 ;
347 fTkaong[1] = 1.64E-8 ;
348 fTkaong[2] = 6.07E-10 ;
349 fTFkaong = new TFormula("ToF response to kaon" , "gaus") ;
350 fTFkaong->SetParameters( fTkaong[0], fTkaong[1], fTkaong[2]) ;
351 //Landau (max probability to inf)
352 fTkaonl[0] = 2.0E9 ;
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)
360 fThhadrong[0] = 2.02E8 ;
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)
366 fThhadronl[0] = 1.10E9 ;
367 fThhadronl[1] = 1.74E-8 ;
368 fThhadronl[2] = 1.00E-9 ;
369 fTFhhadronl = new TFormula("ToF response to heavy hadrons" , "landau") ;
370 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 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 ;
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
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 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
440 //charged hadrons gaus
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
447
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
454
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
460 //Threshold to differentiate between charged and neutral
461 fChargedNeutralThreshold = 1e-5;
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
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
476 for (Int_t i =0; i< AliPID::kSPECIESN ; i++)
477 fInitPID[i] = 1.;
478
479}
480
481//________________________________________________________________________
482void 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 if(strstr(option,"tim"))
488 gBenchmark->Start("PHOSPID");
489
490 if(strstr(option,"print")) {
491 Print() ;
492 return ;
493 }
494
495 if(fTrackSegments && //Skip events, where no track segments made
496 fTrackSegments->GetEntriesFast()) {
497
498 GetVertex() ;
499 MakeRecParticles() ;
500
501 if(fBayesian)
502 MakePID() ;
503
504 if(strstr(option,"deb"))
505 PrintRecParticles(option) ;
506 }
507
508 if(strstr(option,"deb"))
509 PrintRecParticles(option);
510 if(strstr(option,"tim")){
511 gBenchmark->Stop("PHOSPID");
512 AliInfo(Form("took %f seconds for PID",
513 gBenchmark->GetCpuTime("PHOSPID")));
514 }
515}
516
517//________________________________________________________________________
518Double_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 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 ;
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] ;
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 if(TMath::Abs(sigma) > 1.e-10){
538 return cnt*TMath::Gaus(y,mean,sigma);
539 }
540 else
541 return 0.;
542
543}
544//________________________________________________________________________
545Double_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 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
555 if(TMath::Abs(sigma) > 1.e-10){
556 return cnt*TMath::Gaus(y,mean,sigma);
557 }
558 else
559 return 0.;
560
561
562
563}
564
565//____________________________________________________________________________
566const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const
567{
568 //Get file name that contains the PCA for a particle ("photon or pi0")
569 particle.ToLower();
570 TString name;
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()));
578 return name;
579}
580
581//____________________________________________________________________________
582Float_t AliPHOSPIDv1::GetParameterCalibration(Int_t i) const
583{
584 // Get the i-th parameter "Calibration"
585 Float_t param = 0.;
586 if (i>2 || i<0) {
587 AliError(Form("Invalid parameter number: %d",i));
588 } else
589 param = (*fParameters)(0,i);
590 return param;
591}
592
593//____________________________________________________________________________
594Float_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 if(i>2 || i<0) {
599 AliError(Form("Invalid parameter number: %d",i));
600 } else {
601 axis.ToLower();
602 if (axis == "x")
603 param = (*fParameters)(1,i);
604 else if (axis == "z")
605 param = (*fParameters)(2,i);
606 else {
607 AliError(Form("Invalid axis name: %s",axis.Data()));
608 }
609 }
610 return param;
611}
612
613//____________________________________________________________________________
614Float_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 axis.ToLower();
620 Float_t p[]={0.,0.,0.};
621 for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
622 Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
623 return sig;
624}
625
626//____________________________________________________________________________
627Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const
628{
629 // Calculates the parameter param of the ellipse
630
631 particle.ToLower();
632 param. ToLower();
633 Float_t p[4]={0.,0.,0.,0.};
634 Float_t value = 0.0;
635 for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
636 if (particle == "photon") {
637 if (param.Contains("a")) e = TMath::Min((Double_t)e,70.);
638 else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.);
639 else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
640 }
641
642 if (particle == "photon")
643 value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
644 else if (particle == "pi0")
645 value = p[0] + p[1]*e + p[2]*e*e;
646
647 return value;
648}
649
650//_____________________________________________________________________________
651Float_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 if (i>3 || i<0) {
657 AliError(Form("Wrong parameter number: %d\n",i));
658 } else
659 param = (*fParameters)(14,i) ;
660 return param;
661}
662
663//____________________________________________________________________________
664Float_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 if (i>2 || i<0) {
670 AliError(Form("Wrong parameter number: %d\n",i));
671 } else
672 param = (*fParameters)(15,i) ;
673 return param;
674}
675
676//____________________________________________________________________________
677Float_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 if(i>2 || i<0) {
683 AliError(Form("Invalid Efficiency-Purity choice %d",i));
684 } else
685 param = (*fParameters)(3,i) ;
686 return param;
687}
688
689//_____________________________________________________________________________
690Float_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 particle.ToLower();
696 param. ToLower();
697 Int_t offset = -1;
698 if (particle == "photon")
699 offset=0;
700 else if (particle == "pi0")
701 offset=5;
702 else
703 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 if (param.Contains("a")) p=4+offset;
710 else if(param.Contains("b")) p=5+offset;
711 else if(param.Contains("c")) p=6+offset;
712 else if(param.Contains("x0"))p=7+offset;
713 else if(param.Contains("y0"))p=8+offset;
714
715 if (i>4 || i<0) {
716 AliError(Form("No parameter with index %d", i)) ;
717 } else if (p==-1) {
718 AliError(Form("No parameter with name %s", param.Data() )) ;
719 } else
720 par = (*fParameters)(p,i) ;
721
722 return par;
723}
724//____________________________________________________________________________
725Int_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 if(effPur>2 || effPur<0)
729 AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
730
731//DP if(ts->GetCpvIndex()<0)
732//DP return 1 ; //no CPV cluster
733
734 Float_t sigX = GetCpv2EmcDistanceCut("X",e);
735 Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
736
737 Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X"));
738 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 if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1)))
744 return 1;//Neutral
745 else
746 return 0;//Charged
747}
748
749//____________________________________________________________________________
750Int_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 particle.ToLower();
755 Int_t prinbit = 0 ;
756 Float_t a = GetEllipseParameter(particle,"a" , e);
757 Float_t b = GetEllipseParameter(particle,"b" , e);
758 Float_t c = GetEllipseParameter(particle,"c" , e);
759 Float_t x0 = GetEllipseParameter(particle,"x0", e);
760 Float_t y0 = GetEllipseParameter(particle,"y0", e);
761
762 Float_t r = TMath::Power((p[0] - x0)/a,2) +
763 TMath::Power((p[1] - y0)/b,2) +
764 c*(p[0] - x0)*(p[1] - y0)/(a*b) ;
765 //3 different ellipses defined
766 if((effPur==2) && (r<1./2.)) prinbit= 1;
767 if((effPur==1) && (r<2. )) prinbit= 1;
768 if((effPur==0) && (r<9./2.)) prinbit= 1;
769
770 if(r<0)
771 AliError("Negative square?") ;
772
773 return prinbit;
774
775}
776//____________________________________________________________________________
777Int_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 Float_t e = emc->GetEnergy();
783 if (e < 30.0) return 0;
784 Float_t m2x = emc->GetM2x();
785 Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
786 TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
787 TMath::Power(GetParameterPhotonBoundary(2),2)) +
788 GetParameterPhotonBoundary(3);
789 AliDebug(1, Form("E=%f, m2x=%f, boundary=%f", e,m2x,m2xBoundary));
790 if (m2x < m2xBoundary)
791 return 1;// A hard photon
792 else
793 return 0;// Not a hard photon
794}
795
796//____________________________________________________________________________
797Int_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 Float_t e = emc->GetEnergy();
803 if (e < 30.0) return 0;
804 Float_t m2x = emc->GetM2x();
805 Float_t m2xBoundary = GetParameterPi0Boundary(0) +
806 e * GetParameterPi0Boundary(1);
807 AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary));
808 if (m2x > m2xBoundary)
809 return 1;// A hard pi0
810 else
811 return 0;// Not a hard pi0
812}
813
814//____________________________________________________________________________
815TVector3 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 TVector3 local ;
824 emc->GetLocalPosition(local) ;
825
826 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 TVector3 vtxOld(0.,0.,0.) ;
834 TVector3 vInc ;
835 Float_t x=local.X() ;
836 Float_t z=local.Z() ;
837 phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ;
838 Float_t depthxOld = 0.;
839 Float_t depthzOld = 0.;
840 Float_t energy = emc->GetEnergy() ;
841 if (energy > 0 && vInc.Y()!=0.) {
842 depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
843 depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
844 }
845 else{
846 AliError("Cluster with zero energy \n");
847 }
848 //Apply Real vertex
849 phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ;
850 Float_t depthx = 0.;
851 Float_t depthz = 0.;
852 if (energy > 0 && vInc.Y()!=0.) {
853 depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
854 depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
855 }
856
857 //Correct for the vertex position and shower depth
858 Double_t xd=x+(depthxOld-depthx) ;
859 Double_t zd=z+(depthzOld-depthz) ;
860 TVector3 dir(0,0,0) ;
861 phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ;
862
863 dir-=fVtx ;
864 dir.SetMag(1.) ;
865
866 return dir ;
867}
868
869//________________________________________________________________________
870Double_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 if (x > par[9]) x = par[9];
877
878 //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ;
879 Double_t cnt = par[0] + par[1] * x + par[2] * x * x ;
880 Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ;
881 Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ;
882
883 if(TMath::Abs(sigma) > 1.e-10){
884 return cnt*TMath::Landau(y,mean,sigma);
885 }
886 else
887 return 0.;
888
889}
890//________________________________________________________________________
891Double_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 Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ;
899 Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ;
900 Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ;
901
902 if(TMath::Abs(sigma) > 1.e-10){
903 return cnt*TMath::Landau(y,mean,sigma);
904 }
905 else
906 return 0.;
907
908
909}
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//____________________________________________________________________________
943void AliPHOSPIDv1::MakePID()
944{
945 // construct the PID weight from a Bayesian Method
946
947 const Int_t kSPECIES = AliPID::kSPECIESN ;
948
949 Int_t nparticles = fRecParticles->GetEntriesFast() ;
950
951 if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
952 AliFatal("RecPoints or TrackSegments not found !") ;
953 }
954
955 TIter next(fTrackSegments) ;
956 AliPHOSTrackSegment * ts ;
957 Int_t index = 0 ;
958
959 Double_t * stof[kSPECIES] ;
960 Double_t * sdp [kSPECIES] ;
961 Double_t * scpv[kSPECIES] ;
962 Double_t * sw [kSPECIES] ;
963 //Info("MakePID","Begin MakePID");
964
965 for (Int_t i =0; i< kSPECIES; i++){
966 stof[i] = new Double_t[nparticles] ;
967 sdp [i] = new Double_t[nparticles] ;
968 scpv[i] = new Double_t[nparticles] ;
969 sw [i] = new Double_t[nparticles] ;
970 }
971
972
973 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
974
975 //cout<<">>>>>> Bayesian Index "<<index<<endl;
976
977 AliPHOSEmcRecPoint * emc = 0 ;
978 if(ts->GetEmcIndex()>=0)
979 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
980
981// AliPHOSCpvRecPoint * cpv = 0 ;
982// if(ts->GetCpvIndex()>=0)
983// cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
984//
985//// Int_t track = 0 ;
986//// track = ts->GetTrackIndex() ; //TPC tracks ?
987
988 if (!emc) {
989 AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
990 }
991
992
993 // ############Tof#############################
994
995 // Info("MakePID", "TOF");
996 Float_t en = emc->GetEnergy();
997 Double_t time = emc->GetTime() ;
998 // cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl;
999
1000 // now get the signals probability
1001 // s(pid) in the Bayesian formulation
1002
1003 stof[AliPID::kPhoton][index] = 1.;
1004 stof[AliPID::kElectron][index] = 1.;
1005 stof[AliPID::kEleCon][index] = 1.;
1006 //We assing the same prob to charged hadrons, sum is 1
1007 stof[AliPID::kPion][index] = 1./3.;
1008 stof[AliPID::kKaon][index] = 1./3.;
1009 stof[AliPID::kProton][index] = 1./3.;
1010 //We assing the same prob to neutral hadrons, sum is 1
1011 stof[AliPID::kNeutron][index] = 1./2.;
1012 stof[AliPID::kKaon0][index] = 1./2.;
1013 stof[AliPID::kMuon][index] = 1.;
1014
1015 if(en < fTOFEnThreshold) {
1016
1017 Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution
1018 Double_t pTofKaon = 0;
1019
1020 if(time < fTkaonl[1])
1021 pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution
1022 else
1023 pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution
1024
1025 Double_t pTofNucleon = 0;
1026
1027 if(time < fThhadronl[1])
1028 pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution
1029 else
1030 pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution
1031 //We assing the same prob to neutral hadrons, sum is the average prob
1032 Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ;
1033 //We assing the same prob to charged hadrons, sum is the average prob
1034 Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ;
1035
1036 stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ;
1037 //gaus distribution
1038 stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ;
1039 //a conversion electron has the photon ToF
1040 stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ;
1041
1042 stof[AliPID::kElectron][index] = pTofPion ;
1043
1044 stof[AliPID::kPion][index] = pTofChHadron ;
1045 stof[AliPID::kKaon][index] = pTofChHadron ;
1046 stof[AliPID::kProton][index] = pTofChHadron ;
1047
1048 stof[AliPID::kKaon0][index] = pTofNeHadron ;
1049 stof[AliPID::kNeutron][index] = pTofNeHadron ;
1050 }
1051
1052 // Info("MakePID", "Dispersion");
1053
1054 // ###########Shower shape: Dispersion####################
1055 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
1061 sdp[AliPID::kPhoton][index] = 1. ;
1062 sdp[AliPID::kElectron][index] = 1. ;
1063 sdp[AliPID::kPion][index] = 1. ;
1064 sdp[AliPID::kKaon][index] = 1. ;
1065 sdp[AliPID::kProton][index] = 1. ;
1066 sdp[AliPID::kNeutron][index] = 1. ;
1067 sdp[AliPID::kEleCon][index] = 1. ;
1068 sdp[AliPID::kKaon0][index] = 1. ;
1069 sdp[AliPID::kMuon][index] = 1. ;
1070
1071 if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){
1072 sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ;
1073 sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ;
1074 sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ;
1075 sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ;
1076 sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ;
1077 sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ;
1078 sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index];
1079 sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ;
1080 sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ;
1081 //landau distribution
1082 }
1083
1084// Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion);
1085// Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]);
1086// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ;
1087// cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl;
1088
1089 //########## CPV-EMC Distance#######################
1090 // Info("MakePID", "Distance");
1091
1092 Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ;
1093 Float_t z = ts->GetCpvDistance("Z") ;
1094
1095 Double_t pcpv = 0 ;
1096 Double_t pcpvneutral = 0. ;
1097
1098 Double_t elprobx = GausF(en , x, fXelectron) ;
1099 Double_t elprobz = GausF(en , z, fZelectron) ;
1100 Double_t chprobx = GausF(en , x, fXcharged) ;
1101 Double_t chprobz = GausF(en , z, fZcharged) ;
1102 Double_t pcpvelectron = elprobx * elprobz;
1103 Double_t pcpvcharged = chprobx * chprobz;
1104
1105// cout<<">>>>energy "<<en<<endl;
1106// cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl;
1107// cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl;
1108// cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1109
1110 // Is neutral or charged?
1111 if(pcpvelectron >= pcpvcharged)
1112 pcpv = pcpvelectron ;
1113 else
1114 pcpv = pcpvcharged ;
1115
1116 if(pcpv < fChargedNeutralThreshold)
1117 {
1118 pcpvneutral = 1. ;
1119 pcpvcharged = 0. ;
1120 pcpvelectron = 0. ;
1121 }
1122 // else
1123 // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl;
1124
1125 scpv[AliPID::kPion][index] = pcpvcharged ;
1126 scpv[AliPID::kKaon][index] = pcpvcharged ;
1127 scpv[AliPID::kProton][index] = pcpvcharged ;
1128
1129 scpv[AliPID::kMuon][index] = pcpvelectron ;
1130 scpv[AliPID::kElectron][index] = pcpvelectron ;
1131 scpv[AliPID::kEleCon][index] = pcpvelectron ;
1132
1133 scpv[AliPID::kPhoton][index] = pcpvneutral ;
1134 scpv[AliPID::kNeutron][index] = pcpvneutral ;
1135 scpv[AliPID::kKaon0][index] = pcpvneutral ;
1136
1137
1138 // Info("MakePID", "CPV passed");
1139
1140 //############## Pi0 #############################
1141 stof[AliPID::kPi0][index] = 0. ;
1142 scpv[AliPID::kPi0][index] = 0. ;
1143 sdp [AliPID::kPi0][index] = 0. ;
1144
1145 if(en > 30.){
1146 // pi0 are detected via decay photon
1147 stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index];
1148 scpv[AliPID::kPi0][index] = pcpvneutral ;
1149 if(emc->GetMultiplicity() > fDispMultThreshold)
1150 sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ;
1151 //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ;
1152// cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = "
1153// <<emc->GetMultiplicity()<<endl;
1154// cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = "
1155// <<sdp [AliPID::kPi0][index]<<endl;
1156 }
1157
1158
1159
1160
1161 //############## muon #############################
1162
1163 if(en > 0.5){
1164 //Muons deposit few energy
1165 scpv[AliPID::kMuon][index] = 0 ;
1166 stof[AliPID::kMuon][index] = 0 ;
1167 sdp [AliPID::kMuon][index] = 0 ;
1168 }
1169
1170 //Weight to apply to hadrons due to energy reconstruction
1171
1172 Float_t weight = fERecWeight ->Eval(en) ;
1173
1174 sw[AliPID::kPhoton][index] = 1. ;
1175 sw[AliPID::kElectron][index] = 1. ;
1176 sw[AliPID::kPion][index] = weight ;
1177 sw[AliPID::kKaon][index] = weight ;
1178 sw[AliPID::kProton][index] = weight ;
1179 sw[AliPID::kNeutron][index] = weight ;
1180 sw[AliPID::kEleCon][index] = 1. ;
1181 sw[AliPID::kKaon0][index] = weight ;
1182 sw[AliPID::kMuon][index] = weight ;
1183 sw[AliPID::kPi0][index] = 1. ;
1184
1185// if(en > 0.5){
1186// cout<<"######################################################"<<endl;
1187// //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ;
1188// cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ;
1189// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl;
1190// cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl;
1191// cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl;
1192// cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1193
1194// cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index]
1195// <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl;
1196// cout<<"EleCon , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index]
1197// <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl;
1198// cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index]
1199// <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl;
1200// cout<<"Muon , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index]
1201// <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl;
1202// cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index]
1203// <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl;
1204// cout<<"Pion , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index]
1205// <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl;
1206// cout<<"Kaon0 , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index]
1207// <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl;
1208// cout<<"Kaon , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index]
1209// <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl;
1210// cout<<"Neutron , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index]
1211// <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl;
1212// cout<<"Proton , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index]
1213// <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl;
1214// cout<<"######################################################"<<endl;
1215// }
1216 index++;
1217 }
1218
1219 //for (index = 0 ; index < kSPECIES ; index++)
1220 // pid[index] /= nparticles ;
1221
1222
1223 // Info("MakePID", "Total Probability calculation");
1224
1225 for(index = 0 ; index < nparticles ; index ++) {
1226
1227 AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
1228
1229 //Conversion electron?
1230
1231 if(recpar->IsEleCon()){
1232 fInitPID[AliPID::kEleCon] = 1. ;
1233 fInitPID[AliPID::kPhoton] = 0. ;
1234 fInitPID[AliPID::kElectron] = 0. ;
1235 }
1236 else{
1237 fInitPID[AliPID::kEleCon] = 0. ;
1238 fInitPID[AliPID::kPhoton] = 1. ;
1239 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 for (jndex = 0 ; jndex < kSPECIES ; jndex++)
1249 wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1250 sw[jndex][index] * fInitPID[jndex] ;
1251
1252 // cout<<"*************wn "<<wn<<endl;
1253 if (TMath::Abs(wn)>0)
1254 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 recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] *
1266 sw[jndex][index] * scpv[jndex][index] *
1267 fInitPID[jndex] / wn) ;
1268 }
1269 }
1270 // Info("MakePID", "Delete");
1271
1272 for (Int_t i =0; i< kSPECIES; i++){
1273 delete [] stof[i];
1274 delete [] sdp [i];
1275 delete [] scpv[i];
1276 delete [] sw [i];
1277 }
1278 // Info("MakePID","End MakePID");
1279}
1280
1281//____________________________________________________________________________
1282void AliPHOSPIDv1::MakeRecParticles()
1283{
1284 // Makes a RecParticle out of a TrackSegment
1285
1286 if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
1287 AliFatal("RecPoints or TrackSegments not found !") ;
1288 }
1289 fRecParticles->Clear();
1290
1291 TIter next(fTrackSegments) ;
1292 AliPHOSTrackSegment * ts ;
1293 Int_t index = 0 ;
1294 AliPHOSRecParticle * rp ;
1295 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1296 // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl;
1297 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
1298 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
1299 rp->SetTrackSegment(index) ;
1300 rp->SetIndexInList(index) ;
1301
1302 AliPHOSEmcRecPoint * emc = 0 ;
1303 if(ts->GetEmcIndex()>=0)
1304 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
1305
1306 AliPHOSCpvRecPoint * cpv = 0 ;
1307 if(ts->GetCpvIndex()>=0)
1308 cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
1309
1310 Int_t track = 0 ;
1311 track = ts->GetTrackIndex() ;
1312
1313 // Now set type (reconstructed) of the particle
1314
1315 // Choose the cluster energy range
1316
1317 if (!emc) {
1318 AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
1319 }
1320
1321 Float_t e = emc->GetEnergy() ;
1322
1323 Float_t lambda[2]={0.,0.} ;
1324 emc->GetElipsAxis(lambda) ;
1325
1326 if((lambda[0]>0.01) && (lambda[1]>0.01)){
1327 // Looking PCA. Define and calculate the data (X),
1328 // introduce in the function X2P that gives the components (P).
1329
1330 Float_t spher = 0. ;
1331 Float_t emaxdtotal = 0. ;
1332
1333 if((lambda[0]+lambda[1])!=0)
1334 spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
1335
1336 emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
1337
1338 fX[0] = lambda[0] ;
1339 fX[1] = lambda[1] ;
1340 fX[2] = emc->GetDispersion() ;
1341 fX[3] = spher ;
1342 fX[4] = emc->GetMultiplicity() ;
1343 fX[5] = emaxdtotal ;
1344 fX[6] = emc->GetCoreEnergy() ;
1345
1346 fPrincipalPhoton->X2P(fX,fPPhoton);
1347 fPrincipalPi0 ->X2P(fX,fPPi0);
1348
1349 }
1350 else{
1351 fPPhoton[0]=-100.0; //We do not accept clusters with
1352 fPPhoton[1]=-100.0; //one cell as a photon-like
1353 fPPi0[0] =-100.0;
1354 fPPi0[1] =-100.0;
1355 }
1356
1357 Float_t time = emc->GetTime() ;
1358 rp->SetTof(time) ;
1359
1360 // Loop of Efficiency-Purity (the 3 points of purity or efficiency
1361 // are taken into account to set the particle identification)
1362 for(Int_t effPur = 0; effPur < 3 ; effPur++){
1363
1364 // Looking at the CPV detector. If RCPV greater than CpvEmcDistance,
1365 // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
1366 // is set to 1
1367 if(GetCPVBit(ts, effPur,e) == 1 ){
1368 rp->SetPIDBit(effPur) ;
1369 //cout<<"CPV bit "<<effPur<<endl;
1370 }
1371 // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th
1372 // bit (depending on the efficiency-purity point )is set to 1
1373 if(time< (*fParameters)(3,effPur))
1374 rp->SetPIDBit(effPur+3) ;
1375
1376 //Photon PCA
1377 //If we are inside the ellipse, 7th, 8th or 9th
1378 // bit (depending on the efficiency-purity point )is set to 1
1379 if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
1380 rp->SetPIDBit(effPur+6) ;
1381
1382 //Pi0 PCA
1383 //If we are inside the ellipse, 10th, 11th or 12th
1384 // bit (depending on the efficiency-purity point )is set to 1
1385 if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1)
1386 rp->SetPIDBit(effPur+9) ;
1387 }
1388 if(GetHardPhotonBit(emc))
1389 rp->SetPIDBit(12) ;
1390 if(GetHardPi0Bit (emc))
1391 rp->SetPIDBit(13) ;
1392
1393 if(track >= 0)
1394 rp->SetPIDBit(14) ;
1395
1396 //Set momentum, energy and other parameters
1397 TVector3 dir = GetMomentumDirection(emc,cpv) ;
1398 dir.SetMag(e) ;
1399 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
1400 rp->SetCalcMass(0);
1401 rp->Name(); //If photon sets the particle pdg name to gamma
1402 rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
1403 rp->SetFirstMother(-1);
1404 rp->SetLastMother(-1);
1405 rp->SetFirstDaughter(-1);
1406 rp->SetLastDaughter(-1);
1407 rp->SetPolarisation(0,0,0);
1408 //Set the position in global coordinate system from the RecPoint
1409 AliPHOSTrackSegment * ts1 = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
1410 AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
1411 TVector3 pos ;
1412 fGeom->GetGlobalPHOS(erp, pos) ;
1413 rp->SetPos(pos);
1414 index++ ;
1415 }
1416}
1417
1418//____________________________________________________________________________
1419void AliPHOSPIDv1::Print(const Option_t *) const
1420{
1421 // Print the parameters used for the particle type identification
1422
1423 AliInfo("=============== AliPHOSPIDv1 ================") ;
1424 printf("Making PID\n") ;
1425 printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ;
1426 printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ;
1427 printf(" Matrix of Parameters: 14x4\n") ;
1428 printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
1429 printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ;
1430 printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
1431 printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
1432 printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
1433 fParameters->Print() ;
1434}
1435
1436
1437
1438//____________________________________________________________________________
1439void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
1440{
1441 // Print table of reconstructed particles
1442
1443 TString message ;
1444 message = " found " ;
1445 message += fRecParticles->GetEntriesFast();
1446 message += " RecParticles\n" ;
1447
1448 if(strstr(option,"all")) { // printing found TS
1449 message += "\n PARTICLE Index \n" ;
1450
1451 Int_t index ;
1452 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
1453 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
1454 message += "\n" ;
1455 message += rp->Name().Data() ;
1456 message += " " ;
1457 message += rp->GetIndexInList() ;
1458 message += " " ;
1459 message += rp->GetType() ;
1460 }
1461 }
1462 AliInfo(message.Data() ) ;
1463}
1464
1465//____________________________________________________________________________
1466void AliPHOSPIDv1::SetParameters()
1467{
1468 // PCA : To do the Principal Components Analysis it is necessary
1469 // the Principal file, which is opened here
1470 fX = new double[7]; // Data for the PCA
1471 fPPhoton = new double[7]; // Eigenvalues of the PCA
1472 fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA
1473
1474 // Read photon principals from the photon file
1475
1476 fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
1477 TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
1478 fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
1479 f.Close() ;
1480
1481 // Read pi0 principals from the pi0 file
1482
1483 fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
1484 TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
1485 fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
1486 fPi0.Close() ;
1487
1488 // Open parameters file and initialization of the Parameters matrix.
1489 // In the File Parameters.dat are all the parameters. These are introduced
1490 // in a matrix of 16x4
1491 //
1492 // All the parameters defined in this file are, in order of row:
1493 // line 0 : calibration
1494 // lines 1,2 : CPV rectangular cat for X and Z
1495 // line 3 : TOF cut
1496 // lines 4-8 : parameters to calculate photon PCA ellipse
1497 // lines 9-13: parameters to calculate pi0 PCA ellipse
1498 // lines 14-15: parameters to calculate border for high-pt photons and pi0
1499
1500 fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
1501 fParameters = new TMatrixF(16,4) ;
1502 const Int_t kMaxLeng=255;
1503 char string[kMaxLeng];
1504
1505 // Open a text file with PID parameters
1506 FILE *fd = fopen(fFileNameParameters.Data(),"r");
1507 if (!fd)
1508 AliFatal(Form("File %s with a PID parameters cannot be opened\n",
1509 fFileNameParameters.Data()));
1510
1511 Int_t i=0;
1512 // Read parameter file line-by-line and skip empty line and comments
1513 while (fgets(string,kMaxLeng,fd) != NULL) {
1514 if (string[0] == '\n' ) continue;
1515 if (string[0] == '!' ) continue;
1516 sscanf(string, "%f %f %f %f",
1517 &(*fParameters)(i,0), &(*fParameters)(i,1),
1518 &(*fParameters)(i,2), &(*fParameters)(i,3));
1519 i++;
1520 AliDebug(1, Form("Line %d: %s",i,string));
1521 }
1522 fclose(fd);
1523}
1524
1525//____________________________________________________________________________
1526void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param)
1527{
1528 // Set parameter "Calibration" i to a value param
1529 if(i>2 || i<0) {
1530 AliError(Form("Invalid parameter number: %d",i));
1531 } else
1532 (*fParameters)(0,i) = param ;
1533}
1534
1535//____________________________________________________________________________
1536void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut)
1537{
1538 // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on
1539 // Purity-Efficiency point i
1540
1541 if(i>2 || i<0) {
1542 AliError(Form("Invalid parameter number: %d",i));
1543 } else {
1544 axis.ToLower();
1545 if (axis == "x") (*fParameters)(1,i) = cut;
1546 else if (axis == "z") (*fParameters)(2,i) = cut;
1547 else {
1548 AliError(Form("Invalid axis name: %s",axis.Data()));
1549 }
1550 }
1551}
1552
1553//____________________________________________________________________________
1554void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param)
1555{
1556 // Set parameter "Hard photon boundary" i to a value param
1557 if(i>4 || i<0) {
1558 AliError(Form("Invalid parameter number: %d",i));
1559 } else
1560 (*fParameters)(14,i) = param ;
1561}
1562
1563//____________________________________________________________________________
1564void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param)
1565{
1566 // Set parameter "Hard pi0 boundary" i to a value param
1567 if(i>1 || i<0) {
1568 AliError(Form("Invalid parameter number: %d",i));
1569 } else
1570 (*fParameters)(15,i) = param ;
1571}
1572
1573//_____________________________________________________________________________
1574void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate)
1575{
1576 // Set the parameter TimeGate depending on Purity-Efficiency point i
1577 if (i>2 || i<0) {
1578 AliError(Form("Invalid Efficiency-Purity choice %d",i));
1579 } else
1580 (*fParameters)(3,i)= gate ;
1581}
1582
1583//_____________________________________________________________________________
1584void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par)
1585{
1586 // Set the parameter "i" that is needed to calculate the ellipse
1587 // parameter "param" for a particle "particle"
1588
1589 particle.ToLower();
1590 param. ToLower();
1591 Int_t p= -1;
1592 Int_t offset=0;
1593
1594 if (particle == "photon") offset=0;
1595 else if (particle == "pi0") offset=5;
1596 else
1597 AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
1598 particle.Data()));
1599
1600 if (param.Contains("a")) p=4+offset;
1601 else if(param.Contains("b")) p=5+offset;
1602 else if(param.Contains("c")) p=6+offset;
1603 else if(param.Contains("x0"))p=7+offset;
1604 else if(param.Contains("y0"))p=8+offset;
1605 if((i>4)||(i<0)) {
1606 AliError(Form("No parameter with index %d", i)) ;
1607 } else if(p==-1) {
1608 AliError(Form("No parameter with name %s", param.Data() )) ;
1609 } else
1610 (*fParameters)(p,i) = par ;
1611}
1612
1613//____________________________________________________________________________
1614void AliPHOSPIDv1::GetVertex(void)
1615{ //extract vertex either using ESD or generator
1616
1617 //Try to extract vertex from data
1618 if(fESD){
1619 const AliESDVertex *esdVtx = fESD->GetVertex() ;
1620 if(esdVtx && esdVtx->GetChi2()!=0.){
1621 fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
1622 return ;
1623 }
1624 }
1625
1626 // Use vertex diamond from CDB GRP folder if the one from ESD is missing
1627 // PLEASE FIX IT
1628 AliWarning("Can not read vertex from data, use fixed \n") ;
1629 fVtx.SetXYZ(0.,0.,0.) ;
1630
1631}
1632//_______________________________________________________________________
1633void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
1634 // Sets values for the initial population of each particle type
1635 for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i];
1636}
1637//_______________________________________________________________________
1638void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
1639 // Gets values for the initial population of each particle type
1640 for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i];
1641}