Merge branch 'master' into LocalDev
[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.,0.,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.,0.,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.,0.,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::kSPECIESCN ; 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::kSPECIESCN ;
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 //Initialize anused species
1004 for(Int_t iii=0; iii<kSPECIES; iii++)stof[iii][index]=0. ;
1005
1006 stof[AliPID::kPhoton][index] = 1.;
1007 stof[AliPID::kElectron][index] = 1.;
1008 stof[AliPID::kEleCon][index] = 1.;
1009 //We assing the same prob to charged hadrons, sum is 1
1010 stof[AliPID::kPion][index] = 1./3.;
1011 stof[AliPID::kKaon][index] = 1./3.;
1012 stof[AliPID::kProton][index] = 1./3.;
1013 //We assing the same prob to neutral hadrons, sum is 1
1014 stof[AliPID::kNeutron][index] = 1./2.;
1015 stof[AliPID::kKaon0][index] = 1./2.;
1016 stof[AliPID::kMuon][index] = 1.;
1017
1018 if(en < fTOFEnThreshold) {
1019
1020 Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution
1021 Double_t pTofKaon = 0;
1022
1023 if(time < fTkaonl[1])
1024 pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution
1025 else
1026 pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution
1027
1028 Double_t pTofNucleon = 0;
1029
1030 if(time < fThhadronl[1])
1031 pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution
1032 else
1033 pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution
1034 //We assing the same prob to neutral hadrons, sum is the average prob
1035 Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ;
1036 //We assing the same prob to charged hadrons, sum is the average prob
1037 Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ;
1038
1039 stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ;
1040 //gaus distribution
1041 stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ;
1042 //a conversion electron has the photon ToF
1043 stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ;
1044
1045 stof[AliPID::kElectron][index] = pTofPion ;
1046
1047 stof[AliPID::kPion][index] = pTofChHadron ;
1048 stof[AliPID::kKaon][index] = pTofChHadron ;
1049 stof[AliPID::kProton][index] = pTofChHadron ;
1050
1051 stof[AliPID::kKaon0][index] = pTofNeHadron ;
1052 stof[AliPID::kNeutron][index] = pTofNeHadron ;
1053 }
1054
1055 // Info("MakePID", "Dispersion");
1056
1057 // ###########Shower shape: Dispersion####################
1058 Float_t dispersion = emc->GetDispersion();
1059 //DP: Correct for non-perpendicular incidence
1060 //DP: still to be done
1061
1062 //dispersion is not well defined if the cluster is only in few crystals
1063 //Initialize anused species
1064 for(Int_t iii=0; iii<kSPECIES; iii++)sdp[iii][index]=0. ;
1065
1066 sdp[AliPID::kPhoton][index] = 1. ;
1067 sdp[AliPID::kElectron][index] = 1. ;
1068 sdp[AliPID::kPion][index] = 1. ;
1069 sdp[AliPID::kKaon][index] = 1. ;
1070 sdp[AliPID::kProton][index] = 1. ;
1071 sdp[AliPID::kNeutron][index] = 1. ;
1072 sdp[AliPID::kEleCon][index] = 1. ;
1073 sdp[AliPID::kKaon0][index] = 1. ;
1074 sdp[AliPID::kMuon][index] = 1. ;
1075
1076 if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){
1077 sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ;
1078 sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ;
1079 sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ;
1080 sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ;
1081 sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ;
1082 sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ;
1083 sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index];
1084 sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ;
1085 sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ;
1086 //landau distribution
1087 }
1088
1089// Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion);
1090// Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]);
1091// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ;
1092// cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl;
1093
1094 //########## CPV-EMC Distance#######################
1095 // Info("MakePID", "Distance");
1096
1097 Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ;
1098 Float_t z = ts->GetCpvDistance("Z") ;
1099
1100 Double_t pcpv = 0 ;
1101 Double_t pcpvneutral = 0. ;
1102
1103 Double_t elprobx = GausF(en , x, fXelectron) ;
1104 Double_t elprobz = GausF(en , z, fZelectron) ;
1105 Double_t chprobx = GausF(en , x, fXcharged) ;
1106 Double_t chprobz = GausF(en , z, fZcharged) ;
1107 Double_t pcpvelectron = elprobx * elprobz;
1108 Double_t pcpvcharged = chprobx * chprobz;
1109
1110// cout<<">>>>energy "<<en<<endl;
1111// cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl;
1112// cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl;
1113// cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1114
1115 // Is neutral or charged?
1116 if(pcpvelectron >= pcpvcharged)
1117 pcpv = pcpvelectron ;
1118 else
1119 pcpv = pcpvcharged ;
1120
1121 if(pcpv < fChargedNeutralThreshold)
1122 {
1123 pcpvneutral = 1. ;
1124 pcpvcharged = 0. ;
1125 pcpvelectron = 0. ;
1126 }
1127 // else
1128 // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl;
1129 //Initialize anused species
1130 for(Int_t iii=0; iii<kSPECIES; iii++)scpv[iii][index]=0. ;
1131
1132 scpv[AliPID::kPion][index] = pcpvcharged ;
1133 scpv[AliPID::kKaon][index] = pcpvcharged ;
1134 scpv[AliPID::kProton][index] = pcpvcharged ;
1135
1136 scpv[AliPID::kMuon][index] = pcpvelectron ;
1137 scpv[AliPID::kElectron][index] = pcpvelectron ;
1138 scpv[AliPID::kEleCon][index] = pcpvelectron ;
1139
1140 scpv[AliPID::kPhoton][index] = pcpvneutral ;
1141 scpv[AliPID::kNeutron][index] = pcpvneutral ;
1142 scpv[AliPID::kKaon0][index] = pcpvneutral ;
1143
1144
1145 // Info("MakePID", "CPV passed");
1146
1147 //############## Pi0 #############################
1148 stof[AliPID::kPi0][index] = 0. ;
1149 scpv[AliPID::kPi0][index] = 0. ;
1150 sdp [AliPID::kPi0][index] = 0. ;
1151
1152 if(en > 30.){
1153 // pi0 are detected via decay photon
1154 stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index];
1155 scpv[AliPID::kPi0][index] = pcpvneutral ;
1156 if(emc->GetMultiplicity() > fDispMultThreshold)
1157 sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ;
1158 //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ;
1159// cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = "
1160// <<emc->GetMultiplicity()<<endl;
1161// cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = "
1162// <<sdp [AliPID::kPi0][index]<<endl;
1163 }
1164
1165
1166
1167
1168 //############## muon #############################
1169
1170 if(en > 0.5){
1171 //Muons deposit few energy
1172 scpv[AliPID::kMuon][index] = 0 ;
1173 stof[AliPID::kMuon][index] = 0 ;
1174 sdp [AliPID::kMuon][index] = 0 ;
1175 }
1176
1177 //Weight to apply to hadrons due to energy reconstruction
1178 //Initialize anused species
1179 for(Int_t iii=0; iii<kSPECIES; iii++)sw[iii][index]=1. ;
1180
1181 Float_t weight = fERecWeight ->Eval(en) ;
1182
1183 sw[AliPID::kPhoton][index] = 1. ;
1184 sw[AliPID::kElectron][index] = 1. ;
1185 sw[AliPID::kPion][index] = weight ;
1186 sw[AliPID::kKaon][index] = weight ;
1187 sw[AliPID::kProton][index] = weight ;
1188 sw[AliPID::kNeutron][index] = weight ;
1189 sw[AliPID::kEleCon][index] = 1. ;
1190 sw[AliPID::kKaon0][index] = weight ;
1191 sw[AliPID::kMuon][index] = weight ;
1192 sw[AliPID::kPi0][index] = 1. ;
1193
1194// if(en > 0.5){
1195// cout<<"######################################################"<<endl;
1196// //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ;
1197// cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ;
1198// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl;
1199// cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl;
1200// cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl;
1201// cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1202
1203// cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index]
1204// <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl;
1205// cout<<"EleCon , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index]
1206// <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl;
1207// cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index]
1208// <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl;
1209// cout<<"Muon , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index]
1210// <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl;
1211// cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index]
1212// <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl;
1213// cout<<"Pion , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index]
1214// <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl;
1215// cout<<"Kaon0 , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index]
1216// <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl;
1217// cout<<"Kaon , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index]
1218// <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl;
1219// cout<<"Neutron , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index]
1220// <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl;
1221// cout<<"Proton , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index]
1222// <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl;
1223// cout<<"######################################################"<<endl;
1224// }
1225 index++;
1226 }
1227
1228 //for (index = 0 ; index < kSPECIES ; index++)
1229 // pid[index] /= nparticles ;
1230
1231
1232 // Info("MakePID", "Total Probability calculation");
1233
1234 for(index = 0 ; index < nparticles ; index ++) {
1235
1236 AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
1237
1238 //Conversion electron?
1239
1240 if(recpar->IsEleCon()){
1241 fInitPID[AliPID::kEleCon] = 1. ;
1242 fInitPID[AliPID::kPhoton] = 0. ;
1243 fInitPID[AliPID::kElectron] = 0. ;
1244 }
1245 else{
1246 fInitPID[AliPID::kEleCon] = 0. ;
1247 fInitPID[AliPID::kPhoton] = 1. ;
1248 fInitPID[AliPID::kElectron] = 1. ;
1249 }
1250 // fInitPID[AliPID::kEleCon] = 0. ;
1251
1252
1253 // calculates the Bayesian weight
1254
1255 Int_t jndex ;
1256 Double_t wn = 0.0 ;
1257 for (jndex = 0 ; jndex < kSPECIES ; jndex++)
1258 wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1259 sw[jndex][index] * fInitPID[jndex] ;
1260
1261 // cout<<"*************wn "<<wn<<endl;
1262 if (TMath::Abs(wn)>0)
1263 for (jndex = 0 ; jndex < kSPECIES ; jndex++) {
1264 //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn"
1265 //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex] << endl;
1266 //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl;
1267 // if(jndex == AliPID::kPi0 || jndex == AliPID::kPhoton){
1268 // cout<<"Particle "<<jndex<<" final prob * wn "
1269 // <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1270 // fInitPID[jndex] <<" wn "<< wn<<endl;
1271 // cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index]
1272 // <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl;
1273 // }
1274 recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] *
1275 sw[jndex][index] * scpv[jndex][index] *
1276 fInitPID[jndex] / wn) ;
1277 }
1278 }
1279 // Info("MakePID", "Delete");
1280
1281 for (Int_t i =0; i< kSPECIES; i++){
1282 delete [] stof[i];
1283 delete [] sdp [i];
1284 delete [] scpv[i];
1285 delete [] sw [i];
1286 }
1287 // Info("MakePID","End MakePID");
1288}
1289
1290//____________________________________________________________________________
1291void AliPHOSPIDv1::MakeRecParticles()
1292{
1293 // Makes a RecParticle out of a TrackSegment
1294
1295 if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
1296 AliFatal("RecPoints or TrackSegments not found !") ;
1297 }
1298 fRecParticles->Clear();
1299
1300 TIter next(fTrackSegments) ;
1301 AliPHOSTrackSegment * ts ;
1302 Int_t index = 0 ;
1303 AliPHOSRecParticle * rp ;
1304 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1305 // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl;
1306 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
1307 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
1308 rp->SetTrackSegment(index) ;
1309 rp->SetIndexInList(index) ;
1310
1311 AliPHOSEmcRecPoint * emc = 0 ;
1312 if(ts->GetEmcIndex()>=0)
1313 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
1314
1315 AliPHOSCpvRecPoint * cpv = 0 ;
1316 if(ts->GetCpvIndex()>=0)
1317 cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
1318
1319 Int_t track = 0 ;
1320 track = ts->GetTrackIndex() ;
1321
1322 // Now set type (reconstructed) of the particle
1323
1324 // Choose the cluster energy range
1325
1326 if (!emc) {
1327 AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
1328 }
1329
1330 Float_t e = emc->GetEnergy() ;
1331
1332 Float_t lambda[2]={0.,0.} ;
1333 emc->GetElipsAxis(lambda) ;
1334
1335 if((lambda[0]>0.01) && (lambda[1]>0.01)){
1336 // Looking PCA. Define and calculate the data (X),
1337 // introduce in the function X2P that gives the components (P).
1338
1339 Float_t spher = 0. ;
1340 Float_t emaxdtotal = 0. ;
1341
1342 if((lambda[0]+lambda[1])!=0)
1343 spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
1344
1345 emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
1346
1347 fX[0] = lambda[0] ;
1348 fX[1] = lambda[1] ;
1349 fX[2] = emc->GetDispersion() ;
1350 fX[3] = spher ;
1351 fX[4] = emc->GetMultiplicity() ;
1352 fX[5] = emaxdtotal ;
1353 fX[6] = emc->GetCoreEnergy() ;
1354
1355 fPrincipalPhoton->X2P(fX,fPPhoton);
1356 fPrincipalPi0 ->X2P(fX,fPPi0);
1357
1358 }
1359 else{
1360 fPPhoton[0]=-100.0; //We do not accept clusters with
1361 fPPhoton[1]=-100.0; //one cell as a photon-like
1362 fPPi0[0] =-100.0;
1363 fPPi0[1] =-100.0;
1364 }
1365
1366 Float_t time = emc->GetTime() ;
1367 rp->SetTof(time) ;
1368
1369 // Loop of Efficiency-Purity (the 3 points of purity or efficiency
1370 // are taken into account to set the particle identification)
1371 for(Int_t effPur = 0; effPur < 3 ; effPur++){
1372
1373 // Looking at the CPV detector. If RCPV greater than CpvEmcDistance,
1374 // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
1375 // is set to 1
1376 if(GetCPVBit(ts, effPur,e) == 1 ){
1377 rp->SetPIDBit(effPur) ;
1378 //cout<<"CPV bit "<<effPur<<endl;
1379 }
1380 // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th
1381 // bit (depending on the efficiency-purity point )is set to 1
1382 if(time< (*fParameters)(3,effPur))
1383 rp->SetPIDBit(effPur+3) ;
1384
1385 //Photon PCA
1386 //If we are inside the ellipse, 7th, 8th or 9th
1387 // bit (depending on the efficiency-purity point )is set to 1
1388 if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
1389 rp->SetPIDBit(effPur+6) ;
1390
1391 //Pi0 PCA
1392 //If we are inside the ellipse, 10th, 11th or 12th
1393 // bit (depending on the efficiency-purity point )is set to 1
1394 if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1)
1395 rp->SetPIDBit(effPur+9) ;
1396 }
1397 if(GetHardPhotonBit(emc))
1398 rp->SetPIDBit(12) ;
1399 if(GetHardPi0Bit (emc))
1400 rp->SetPIDBit(13) ;
1401
1402 if(track >= 0)
1403 rp->SetPIDBit(14) ;
1404
1405 //Set momentum, energy and other parameters
1406 TVector3 dir = GetMomentumDirection(emc,cpv) ;
1407 dir.SetMag(e) ;
1408 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
1409 rp->SetCalcMass(0);
1410 rp->Name(); //If photon sets the particle pdg name to gamma
1411 rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
1412 rp->SetFirstMother(-1);
1413 rp->SetLastMother(-1);
1414 rp->SetFirstDaughter(-1);
1415 rp->SetLastDaughter(-1);
1416 rp->SetPolarisation(0,0,0);
1417 //Set the position in global coordinate system from the RecPoint
1418 AliPHOSTrackSegment * ts1 = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
1419 AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
1420 TVector3 pos ;
1421 fGeom->GetGlobalPHOS(erp, pos) ;
1422 rp->SetPos(pos);
1423 index++ ;
1424 }
1425}
1426
1427//____________________________________________________________________________
1428void AliPHOSPIDv1::Print(const Option_t *) const
1429{
1430 // Print the parameters used for the particle type identification
1431
1432 AliInfo("=============== AliPHOSPIDv1 ================") ;
1433 printf("Making PID\n") ;
1434 printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ;
1435 printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ;
1436 printf(" Matrix of Parameters: 14x4\n") ;
1437 printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
1438 printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ;
1439 printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
1440 printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
1441 printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
1442 fParameters->Print() ;
1443}
1444
1445
1446
1447//____________________________________________________________________________
1448void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
1449{
1450 // Print table of reconstructed particles
1451
1452 TString message ;
1453 message = " found " ;
1454 message += fRecParticles->GetEntriesFast();
1455 message += " RecParticles\n" ;
1456
1457 if(strstr(option,"all")) { // printing found TS
1458 message += "\n PARTICLE Index \n" ;
1459
1460 Int_t index ;
1461 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
1462 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
1463 message += "\n" ;
1464 message += rp->Name().Data() ;
1465 message += " " ;
1466 message += rp->GetIndexInList() ;
1467 message += " " ;
1468 message += rp->GetType() ;
1469 }
1470 }
1471 AliInfo(message.Data() ) ;
1472}
1473
1474//____________________________________________________________________________
1475void AliPHOSPIDv1::SetParameters()
1476{
1477 // PCA : To do the Principal Components Analysis it is necessary
1478 // the Principal file, which is opened here
1479 fX = new double[7]; // Data for the PCA
1480 fPPhoton = new double[7]; // Eigenvalues of the PCA
1481 fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA
1482
1483 // Read photon principals from the photon file
1484
1485 fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
1486 TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
1487 fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
1488 f.Close() ;
1489
1490 // Read pi0 principals from the pi0 file
1491
1492 fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
1493 TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
1494 fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
1495 fPi0.Close() ;
1496
1497 // Open parameters file and initialization of the Parameters matrix.
1498 // In the File Parameters.dat are all the parameters. These are introduced
1499 // in a matrix of 16x4
1500 //
1501 // All the parameters defined in this file are, in order of row:
1502 // line 0 : calibration
1503 // lines 1,2 : CPV rectangular cat for X and Z
1504 // line 3 : TOF cut
1505 // lines 4-8 : parameters to calculate photon PCA ellipse
1506 // lines 9-13: parameters to calculate pi0 PCA ellipse
1507 // lines 14-15: parameters to calculate border for high-pt photons and pi0
1508
1509 fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
1510 fParameters = new TMatrixF(16,4) ;
1511 const Int_t kMaxLeng=255;
1512 char string[kMaxLeng];
1513
1514 // Open a text file with PID parameters
1515 FILE *fd = fopen(fFileNameParameters.Data(),"r");
1516 if (!fd)
1517 AliFatal(Form("File %s with a PID parameters cannot be opened\n",
1518 fFileNameParameters.Data()));
1519
1520 Int_t i=0;
1521 // Read parameter file line-by-line and skip empty line and comments
1522 while (fgets(string,kMaxLeng,fd) != NULL) {
1523 if (string[0] == '\n' ) continue;
1524 if (string[0] == '!' ) continue;
1525 sscanf(string, "%f %f %f %f",
1526 &(*fParameters)(i,0), &(*fParameters)(i,1),
1527 &(*fParameters)(i,2), &(*fParameters)(i,3));
1528 i++;
1529 AliDebug(1, Form("Line %d: %s",i,string));
1530 }
1531 fclose(fd);
1532}
1533
1534//____________________________________________________________________________
1535void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param)
1536{
1537 // Set parameter "Calibration" i to a value param
1538 if(i>2 || i<0) {
1539 AliError(Form("Invalid parameter number: %d",i));
1540 } else
1541 (*fParameters)(0,i) = param ;
1542}
1543
1544//____________________________________________________________________________
1545void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut)
1546{
1547 // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on
1548 // Purity-Efficiency point i
1549
1550 if(i>2 || i<0) {
1551 AliError(Form("Invalid parameter number: %d",i));
1552 } else {
1553 axis.ToLower();
1554 if (axis == "x") (*fParameters)(1,i) = cut;
1555 else if (axis == "z") (*fParameters)(2,i) = cut;
1556 else {
1557 AliError(Form("Invalid axis name: %s",axis.Data()));
1558 }
1559 }
1560}
1561
1562//____________________________________________________________________________
1563void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param)
1564{
1565 // Set parameter "Hard photon boundary" i to a value param
1566 if(i>4 || i<0) {
1567 AliError(Form("Invalid parameter number: %d",i));
1568 } else
1569 (*fParameters)(14,i) = param ;
1570}
1571
1572//____________________________________________________________________________
1573void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param)
1574{
1575 // Set parameter "Hard pi0 boundary" i to a value param
1576 if(i>1 || i<0) {
1577 AliError(Form("Invalid parameter number: %d",i));
1578 } else
1579 (*fParameters)(15,i) = param ;
1580}
1581
1582//_____________________________________________________________________________
1583void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate)
1584{
1585 // Set the parameter TimeGate depending on Purity-Efficiency point i
1586 if (i>2 || i<0) {
1587 AliError(Form("Invalid Efficiency-Purity choice %d",i));
1588 } else
1589 (*fParameters)(3,i)= gate ;
1590}
1591
1592//_____________________________________________________________________________
1593void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par)
1594{
1595 // Set the parameter "i" that is needed to calculate the ellipse
1596 // parameter "param" for a particle "particle"
1597
1598 particle.ToLower();
1599 param. ToLower();
1600 Int_t p= -1;
1601 Int_t offset=0;
1602
1603 if (particle == "photon") offset=0;
1604 else if (particle == "pi0") offset=5;
1605 else
1606 AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
1607 particle.Data()));
1608
1609 if (param.Contains("a")) p=4+offset;
1610 else if(param.Contains("b")) p=5+offset;
1611 else if(param.Contains("c")) p=6+offset;
1612 else if(param.Contains("x0"))p=7+offset;
1613 else if(param.Contains("y0"))p=8+offset;
1614 if((i>4)||(i<0)) {
1615 AliError(Form("No parameter with index %d", i)) ;
1616 } else if(p==-1) {
1617 AliError(Form("No parameter with name %s", param.Data() )) ;
1618 } else
1619 (*fParameters)(p,i) = par ;
1620}
1621
1622//____________________________________________________________________________
1623void AliPHOSPIDv1::GetVertex(void)
1624{ //extract vertex either using ESD or generator
1625
1626 //Try to extract vertex from data
1627 if(fESD){
1628 const AliESDVertex *esdVtx = fESD->GetVertex() ;
1629 if(esdVtx && esdVtx->GetChi2()!=0.){
1630 fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
1631 return ;
1632 }
1633 }
1634
1635 // Use vertex diamond from CDB GRP folder if the one from ESD is missing
1636 // PLEASE FIX IT
1637 AliWarning("Can not read vertex from data, use fixed \n") ;
1638 fVtx.SetXYZ(0.,0.,0.) ;
1639
1640}
1641//_______________________________________________________________________
1642void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
1643 // Sets values for the initial population of each particle type
1644 for (Int_t i=0; i<AliPID::kSPECIESCN; i++) fInitPID[i] = p[i];
1645}
1646//_______________________________________________________________________
1647void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
1648 // Gets values for the initial population of each particle type
1649 for (Int_t i=0; i<AliPID::kSPECIESCN; i++) p[i] = fInitPID[i];
1650}