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