Merge branch 'master' into LocalDev
[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),
8f356520 163 fVtx(0.,0.,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),
8f356520 198 fVtx(0.,0.,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),
8f356520 233 fVtx(0.,0.,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
00a38d07 476 for (Int_t i =0; i< AliPID::kSPECIESCN ; 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
00a38d07 947 const Int_t kSPECIES = AliPID::kSPECIESCN ;
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
11d62687 965 for (Int_t i =0; i< kSPECIES; i++){
966 stof[i] = new Double_t[nparticles] ;
35adb638 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() ;
10d07555 998 // cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl;
fb7b51ad 999
35adb638 1000 // now get the signals probability
1001 // s(pid) in the Bayesian formulation
10d07555 1002
1003 //Initialize anused species
11d62687 1004 for(Int_t iii=0; iii<kSPECIES; iii++)stof[iii][index]=0. ;
cc1fe362 1005
304864ab 1006 stof[AliPID::kPhoton][index] = 1.;
1007 stof[AliPID::kElectron][index] = 1.;
304864ab 1008 stof[AliPID::kEleCon][index] = 1.;
7fb46731 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.;
304864ab 1016 stof[AliPID::kMuon][index] = 1.;
2924941c 1017
1018 if(en < fTOFEnThreshold) {
7fb46731 1019
1020 Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution
1021 Double_t pTofKaon = 0;
1022
35adb638 1023 if(time < fTkaonl[1])
fb7b51ad 1024 pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution
35adb638 1025 else
fb7b51ad 1026 pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution
7fb46731 1027
1028 Double_t pTofNucleon = 0;
1029
35adb638 1030 if(time < fThhadronl[1])
7fb46731 1031 pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution
35adb638 1032 else
7fb46731 1033 pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution
7fb46731 1034 //We assing the same prob to neutral hadrons, sum is the average prob
fb7b51ad 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. ;
7fb46731 1038
fb7b51ad 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
304864ab 1043 stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ;
7fb46731 1044
fb7b51ad 1045 stof[AliPID::kElectron][index] = pTofPion ;
7fb46731 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 ;
cc1fe362 1053 }
c947e71a 1054
1055 // Info("MakePID", "Dispersion");
cc1fe362 1056
7fb46731 1057 // ###########Shower shape: Dispersion####################
35adb638 1058 Float_t dispersion = emc->GetDispersion();
407d15b3 1059 //DP: Correct for non-perpendicular incidence
1060 //DP: still to be done
1061
35adb638 1062 //dispersion is not well defined if the cluster is only in few crystals
10d07555 1063 //Initialize anused species
11d62687 1064 for(Int_t iii=0; iii<kSPECIES; iii++)sdp[iii][index]=0. ;
cc1fe362 1065
304864ab 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. ;
c947e71a 1075
2924941c 1076 if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){
7fb46731 1077 sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ;
304864ab 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] ;
fb7b51ad 1085 sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ;
1086 //landau distribution
35adb638 1087 }
cc1fe362 1088
7fb46731 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]);
c947e71a 1091// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ;
304864ab 1092// cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl;
c947e71a 1093
7fb46731 1094 //########## CPV-EMC Distance#######################
1095 // Info("MakePID", "Distance");
fb7b51ad 1096
26aa7e4a 1097 Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ;
1098 Float_t z = ts->GetCpvDistance("Z") ;
fb7b51ad 1099
7fb46731 1100 Double_t pcpv = 0 ;
c947e71a 1101 Double_t pcpvneutral = 0. ;
7fb46731 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) ;
c947e71a 1107 Double_t pcpvelectron = elprobx * elprobz;
1108 Double_t pcpvcharged = chprobx * chprobz;
7fb46731 1109
1110// cout<<">>>>energy "<<en<<endl;
c947e71a 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
2924941c 1115 // Is neutral or charged?
35adb638 1116 if(pcpvelectron >= pcpvcharged)
1117 pcpv = pcpvelectron ;
1118 else
1119 pcpv = pcpvcharged ;
1120
fb7b51ad 1121 if(pcpv < fChargedNeutralThreshold)
35adb638 1122 {
1123 pcpvneutral = 1. ;
1124 pcpvcharged = 0. ;
1125 pcpvelectron = 0. ;
1126 }
c947e71a 1127 // else
1128 // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl;
10d07555 1129 //Initialize anused species
11d62687 1130 for(Int_t iii=0; iii<kSPECIES; iii++)scpv[iii][index]=0. ;
35adb638 1131
304864ab 1132 scpv[AliPID::kPion][index] = pcpvcharged ;
1133 scpv[AliPID::kKaon][index] = pcpvcharged ;
1134 scpv[AliPID::kProton][index] = pcpvcharged ;
7fb46731 1135
1136 scpv[AliPID::kMuon][index] = pcpvelectron ;
304864ab 1137 scpv[AliPID::kElectron][index] = pcpvelectron ;
304864ab 1138 scpv[AliPID::kEleCon][index] = pcpvelectron ;
7fb46731 1139
1140 scpv[AliPID::kPhoton][index] = pcpvneutral ;
1141 scpv[AliPID::kNeutron][index] = pcpvneutral ;
304864ab 1142 scpv[AliPID::kKaon0][index] = pcpvneutral ;
7fb46731 1143
35adb638 1144
1145 // Info("MakePID", "CPV passed");
c947e71a 1146
7fb46731 1147 //############## Pi0 #############################
304864ab 1148 stof[AliPID::kPi0][index] = 0. ;
1149 scpv[AliPID::kPi0][index] = 0. ;
1150 sdp [AliPID::kPi0][index] = 0. ;
c947e71a 1151
35adb638 1152 if(en > 30.){
1153 // pi0 are detected via decay photon
2924941c 1154 stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index];
304864ab 1155 scpv[AliPID::kPi0][index] = pcpvneutral ;
2924941c 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;
35adb638 1163 }
1164
2924941c 1165
7fb46731 1166
2924941c 1167
7fb46731 1168 //############## muon #############################
1169
35adb638 1170 if(en > 0.5){
1171 //Muons deposit few energy
304864ab 1172 scpv[AliPID::kMuon][index] = 0 ;
1173 stof[AliPID::kMuon][index] = 0 ;
1174 sdp [AliPID::kMuon][index] = 0 ;
c947e71a 1175 }
1176
fb7b51ad 1177 //Weight to apply to hadrons due to energy reconstruction
10d07555 1178 //Initialize anused species
11d62687 1179 for(Int_t iii=0; iii<kSPECIES; iii++)sw[iii][index]=1. ;
fb7b51ad 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
7fb46731 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;
c947e71a 1202
2924941c 1203// cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index]
1204// <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl;
7fb46731 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;
2924941c 1211// cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index]
1212// <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl;
7fb46731 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++;
cc1fe362 1226 }
35adb638 1227
1228 //for (index = 0 ; index < kSPECIES ; index++)
1229 // pid[index] /= nparticles ;
1230
7fb46731 1231
35adb638 1232 // Info("MakePID", "Total Probability calculation");
1233
cc1fe362 1234 for(index = 0 ; index < nparticles ; index ++) {
2924941c 1235
9a2cdbdf 1236 AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
2924941c 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
cc1fe362 1253 // calculates the Bayesian weight
7fb46731 1254
cc1fe362 1255 Int_t jndex ;
1256 Double_t wn = 0.0 ;
1257 for (jndex = 0 ; jndex < kSPECIES ; jndex++)
2924941c 1258 wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1259 sw[jndex][index] * fInitPID[jndex] ;
1260
7fb46731 1261 // cout<<"*************wn "<<wn<<endl;
e74ea0e9 1262 if (TMath::Abs(wn)>0)
1263 for (jndex = 0 ; jndex < kSPECIES ; jndex++) {
35adb638 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;
2924941c 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) ;
e74ea0e9 1277 }
2cc71c1e 1278 }
35adb638 1279 // Info("MakePID", "Delete");
1280
2924941c 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 }
35adb638 1287 // Info("MakePID","End MakePID");
2cc71c1e 1288}
1289
1290//____________________________________________________________________________
e3817e5f 1291void AliPHOSPIDv1::MakeRecParticles()
1292{
b2a60966 1293 // Makes a RecParticle out of a TrackSegment
148b2bba 1294
9a2cdbdf 1295 if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
351dd634 1296 AliFatal("RecPoints or TrackSegments not found !") ;
148b2bba 1297 }
9a2cdbdf 1298 fRecParticles->Clear();
148b2bba 1299
9a2cdbdf 1300 TIter next(fTrackSegments) ;
7acf6008 1301 AliPHOSTrackSegment * ts ;
6ad0bfa0 1302 Int_t index = 0 ;
09fc14a0 1303 AliPHOSRecParticle * rp ;
7acf6008 1304 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
7fb46731 1305 // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl;
9a2cdbdf 1306 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
1307 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
f0a4c9e9 1308 rp->SetTrackSegment(index) ;
9688c1dd 1309 rp->SetIndexInList(index) ;
148b2bba 1310
7acf6008 1311 AliPHOSEmcRecPoint * emc = 0 ;
1312 if(ts->GetEmcIndex()>=0)
9a2cdbdf 1313 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
fad3e5b9 1314
8d4608b5 1315 AliPHOSCpvRecPoint * cpv = 0 ;
7acf6008 1316 if(ts->GetCpvIndex()>=0)
9a2cdbdf 1317 cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
fad3e5b9 1318
bd76890a 1319 Int_t track = 0 ;
1320 track = ts->GetTrackIndex() ;
1321
148b2bba 1322 // Now set type (reconstructed) of the particle
1323
1324 // Choose the cluster energy range
9fa5f1d0 1325
fbf811ec 1326 if (!emc) {
a675b8d6 1327 AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
fbf811ec 1328 }
50739f15 1329
e3817e5f 1330 Float_t e = emc->GetEnergy() ;
bc0c084c 1331
25312a8e 1332 Float_t lambda[2]={0.,0.} ;
6f969528 1333 emc->GetElipsAxis(lambda) ;
407d15b3 1334
50739f15 1335 if((lambda[0]>0.01) && (lambda[1]>0.01)){
1336 // Looking PCA. Define and calculate the data (X),
bc0c084c 1337 // introduce in the function X2P that gives the components (P).
1338
c947e71a 1339 Float_t spher = 0. ;
1340 Float_t emaxdtotal = 0. ;
50739f15 1341
bc0c084c 1342 if((lambda[0]+lambda[1])!=0)
e4df4b30 1343 spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
50739f15 1344
c947e71a 1345 emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
50739f15 1346
1347 fX[0] = lambda[0] ;
1348 fX[1] = lambda[1] ;
1349 fX[2] = emc->GetDispersion() ;
c947e71a 1350 fX[3] = spher ;
50739f15 1351 fX[4] = emc->GetMultiplicity() ;
c947e71a 1352 fX[5] = emaxdtotal ;
50739f15 1353 fX[6] = emc->GetCoreEnergy() ;
1354
e3817e5f 1355 fPrincipalPhoton->X2P(fX,fPPhoton);
1356 fPrincipalPi0 ->X2P(fX,fPPi0);
1f0e7ccd 1357
50739f15 1358 }
1359 else{
e3817e5f 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;
50739f15 1364 }
1365
2cc71c1e 1366 Float_t time = emc->GetTime() ;
1367 rp->SetTof(time) ;
9fa5f1d0 1368
bc0c084c 1369 // Loop of Efficiency-Purity (the 3 points of purity or efficiency
1370 // are taken into account to set the particle identification)
e3817e5f 1371 for(Int_t effPur = 0; effPur < 3 ; effPur++){
50739f15 1372
bc0c084c 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
26aa7e4a 1376 if(GetCPVBit(ts, effPur,e) == 1 ){
e3817e5f 1377 rp->SetPIDBit(effPur) ;
35adb638 1378 //cout<<"CPV bit "<<effPur<<endl;
1379 }
50739f15 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
2cc71c1e 1382 if(time< (*fParameters)(3,effPur))
e3817e5f 1383 rp->SetPIDBit(effPur+3) ;
2cc71c1e 1384
e3817e5f 1385 //Photon PCA
50739f15 1386 //If we are inside the ellipse, 7th, 8th or 9th
1387 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 1388 if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
1389 rp->SetPIDBit(effPur+6) ;
1f0e7ccd 1390
e3817e5f 1391 //Pi0 PCA
1f0e7ccd 1392 //If we are inside the ellipse, 10th, 11th or 12th
1393 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 1394 if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1)
1395 rp->SetPIDBit(effPur+9) ;
f0a4c9e9 1396 }
e3817e5f 1397 if(GetHardPhotonBit(emc))
1398 rp->SetPIDBit(12) ;
1399 if(GetHardPi0Bit (emc))
1400 rp->SetPIDBit(13) ;
1f0e7ccd 1401
bd76890a 1402 if(track >= 0)
1403 rp->SetPIDBit(14) ;
1404
9fa5f1d0 1405 //Set momentum, energy and other parameters
9fa5f1d0 1406 TVector3 dir = GetMomentumDirection(emc,cpv) ;
25312a8e 1407 dir.SetMag(e) ;
1408 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
9fa5f1d0 1409 rp->SetCalcMass(0);
e0ed2e49 1410 rp->Name(); //If photon sets the particle pdg name to gamma
407d15b3 1411 rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
e747b8da 1412 rp->SetFirstMother(-1);
1413 rp->SetLastMother(-1);
1414 rp->SetFirstDaughter(-1);
1415 rp->SetLastDaughter(-1);
1416 rp->SetPolarisation(0,0,0);
d956e9b7 1417 //Set the position in global coordinate system from the RecPoint
605cb8bb 1418 AliPHOSTrackSegment * ts1 = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
1419 AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
d956e9b7 1420 TVector3 pos ;
9a2cdbdf 1421 fGeom->GetGlobalPHOS(erp, pos) ;
d956e9b7 1422 rp->SetPos(pos);
6ad0bfa0 1423 index++ ;
1424 }
6ad0bfa0 1425}
e3817e5f 1426
09fc14a0 1427//____________________________________________________________________________
702ab87e 1428void AliPHOSPIDv1::Print(const Option_t *) const
09fc14a0 1429{
b2a60966 1430 // Print the parameters used for the particle type identification
bc0c084c 1431
351dd634 1432 AliInfo("=============== AliPHOSPIDv1 ================") ;
88cb7938 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") ;
407d15b3 1441 printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
50739f15 1442 fParameters->Print() ;
09fc14a0 1443}
1444
8d0f3f77 1445
69183710 1446
7acf6008 1447//____________________________________________________________________________
a4e98857 1448void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
1449{
dd5c4038 1450 // Print table of reconstructed particles
1451
21cd0c07 1452 TString message ;
9a2cdbdf 1453 message = " found " ;
1454 message += fRecParticles->GetEntriesFast();
3bf72d32 1455 message += " RecParticles\n" ;
1456
7acf6008 1457 if(strstr(option,"all")) { // printing found TS
3bf72d32 1458 message += "\n PARTICLE Index \n" ;
7acf6008 1459
1460 Int_t index ;
9a2cdbdf 1461 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
1462 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
3bf72d32 1463 message += "\n" ;
1464 message += rp->Name().Data() ;
1465 message += " " ;
1466 message += rp->GetIndexInList() ;
1467 message += " " ;
1468 message += rp->GetType() ;
7acf6008 1469 }
3bf72d32 1470 }
351dd634 1471 AliInfo(message.Data() ) ;
69183710 1472}
88cb7938 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");
e8d02863 1510 fParameters = new TMatrixF(16,4) ;
c947e71a 1511 const Int_t kMaxLeng=255;
1512 char string[kMaxLeng];
88cb7938 1513
1514 // Open a text file with PID parameters
1515 FILE *fd = fopen(fFileNameParameters.Data(),"r");
1516 if (!fd)
351dd634 1517 AliFatal(Form("File %s with a PID parameters cannot be opened\n",
1518 fFileNameParameters.Data()));
88cb7938 1519
1520 Int_t i=0;
1521 // Read parameter file line-by-line and skip empty line and comments
c947e71a 1522 while (fgets(string,kMaxLeng,fd) != NULL) {
88cb7938 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++;
a675b8d6 1529 AliDebug(1, Form("Line %d: %s",i,string));
88cb7938 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
351dd634 1538 if(i>2 || i<0) {
1539 AliError(Form("Invalid parameter number: %d",i));
1540 } else
88cb7938 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
351dd634 1550 if(i>2 || i<0) {
1551 AliError(Form("Invalid parameter number: %d",i));
1552 } else {
88cb7938 1553 axis.ToLower();
1554 if (axis == "x") (*fParameters)(1,i) = cut;
1555 else if (axis == "z") (*fParameters)(2,i) = cut;
351dd634 1556 else {
1557 AliError(Form("Invalid axis name: %s",axis.Data()));
1558 }
88cb7938 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
351dd634 1566 if(i>4 || i<0) {
1567 AliError(Form("Invalid parameter number: %d",i));
1568 } else
88cb7938 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
351dd634 1576 if(i>1 || i<0) {
1577 AliError(Form("Invalid parameter number: %d",i));
1578 } else
88cb7938 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
351dd634 1586 if (i>2 || i<0) {
1587 AliError(Form("Invalid Efficiency-Purity choice %d",i));
1588 } else
88cb7938 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
351dd634 1606 AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
1607 particle.Data()));
88cb7938 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;
351dd634 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
88cb7938 1619 (*fParameters)(p,i) = par ;
1620}
1621
1622//____________________________________________________________________________
407d15b3 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() ;
e5b7b511 1629 if(esdVtx && esdVtx->GetChi2()!=0.){
407d15b3 1630 fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
1631 return ;
1632 }
1633 }
9a2cdbdf 1634
1635 // Use vertex diamond from CDB GRP folder if the one from ESD is missing
1636 // PLEASE FIX IT
407d15b3 1637 AliWarning("Can not read vertex from data, use fixed \n") ;
1638 fVtx.SetXYZ(0.,0.,0.) ;
1639
1640}
35adb638 1641//_______________________________________________________________________
1642void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
1643 // Sets values for the initial population of each particle type
00a38d07 1644 for (Int_t i=0; i<AliPID::kSPECIESCN; i++) fInitPID[i] = p[i];
35adb638 1645}
1646//_______________________________________________________________________
1647void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
1648 // Gets values for the initial population of each particle type
00a38d07 1649 for (Int_t i=0; i<AliPID::kSPECIESCN; i++) p[i] = fInitPID[i];
35adb638 1650}