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