Compilation warnings fixed
[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
1cb7c1ee 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
581354c5 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
7acf6008 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
8d0f3f77 275//____________________________________________________________________________
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
69183710 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
88cb7938 593//____________________________________________________________________________
fc7e2f43 594Float_t AliPHOSPIDv1::GetCalibratedEnergy(Float_t e) const
88cb7938 595{
596// It calibrates Energy depending on the recpoint energy.
597// The energy of the reconstructed cluster is corrected with
598// the formula A + B* E + C* E^2, whose parameters where obtained
599// through the study of the reconstructed energy distribution of
600// monoenergetic photons.
71994f35 601
602 if(!fEnergyCorrectionOn) return e;
603
88cb7938 604 Float_t p[]={0.,0.,0.};
605 for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
606 Float_t enerec = p[0] + p[1]*e + p[2]*e*e;
607 return enerec ;
608
609}
610
e3817e5f 611//____________________________________________________________________________
fc7e2f43 612Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const
e3817e5f 613{
614 // Get the i-th parameter "CPV-EMC distance" for the specified axis
615 Float_t param = 0.;
351dd634 616 if(i>2 || i<0) {
617 AliError(Form("Invalid parameter number: %d",i));
618 } else {
e3817e5f 619 axis.ToLower();
351dd634 620 if (axis == "x")
621 param = (*fParameters)(1,i);
622 else if (axis == "z")
623 param = (*fParameters)(2,i);
624 else {
625 AliError(Form("Invalid axis name: %s",axis.Data()));
626 }
e3817e5f 627 }
628 return param;
629}
630
631//____________________________________________________________________________
fc7e2f43 632Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const
e3817e5f 633{
88cb7938 634 // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and
635 // Purity-Efficiency point
636
637 axis.ToLower();
638 Float_t p[]={0.,0.,0.};
639 for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
640 Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
641 return sig;
e3817e5f 642}
643
88cb7938 644//____________________________________________________________________________
fc7e2f43 645Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const
88cb7938 646{
647 // Calculates the parameter param of the ellipse
e3817e5f 648
649 particle.ToLower();
650 param. ToLower();
88cb7938 651 Float_t p[4]={0.,0.,0.,0.};
652 Float_t value = 0.0;
653 for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
654 if (particle == "photon") {
655 if (param.Contains("a")) e = TMath::Min((Double_t)e,70.);
656 else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.);
657 else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
658 }
e3817e5f 659
443caba9 660 if (particle == "photon")
661 value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
662 else if (particle == "pi0")
663 value = p[0] + p[1]*e + p[2]*e*e;
664
88cb7938 665 return value;
e3817e5f 666}
667
668//_____________________________________________________________________________
fc7e2f43 669Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const
e3817e5f 670{
671 // Get the parameter "i" to calculate the boundary on the moment M2x
672 // for photons at high p_T
673 Float_t param = 0;
351dd634 674 if (i>3 || i<0) {
675 AliError(Form("Wrong parameter number: %d\n",i));
676 } else
e3817e5f 677 param = (*fParameters)(14,i) ;
678 return param;
148b2bba 679}
e3817e5f 680
148b2bba 681//____________________________________________________________________________
fc7e2f43 682Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const
e3817e5f 683{
684 // Get the parameter "i" to calculate the boundary on the moment M2x
685 // for pi0 at high p_T
686 Float_t param = 0;
351dd634 687 if (i>2 || i<0) {
688 AliError(Form("Wrong parameter number: %d\n",i));
689 } else
e3817e5f 690 param = (*fParameters)(15,i) ;
691 return param;
692}
148b2bba 693
e3817e5f 694//____________________________________________________________________________
fc7e2f43 695Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const
e3817e5f 696{
88cb7938 697 // Get TimeGate parameter depending on Purity-Efficiency i:
698 // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity
699 Float_t param = 0.;
351dd634 700 if(i>2 || i<0) {
701 AliError(Form("Invalid Efficiency-Purity choice %d",i));
702 } else
88cb7938 703 param = (*fParameters)(3,i) ;
704 return param;
e3817e5f 705}
706
148b2bba 707//_____________________________________________________________________________
fc7e2f43 708Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const
88cb7938 709{
710 // Get the parameter "i" that is needed to calculate the ellipse
711 // parameter "param" for the particle "particle" ("photon" or "pi0")
712
e3817e5f 713 particle.ToLower();
714 param. ToLower();
88cb7938 715 Int_t offset = -1;
351dd634 716 if (particle == "photon")
717 offset=0;
718 else if (particle == "pi0")
719 offset=5;
e3817e5f 720 else
351dd634 721 AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
722 particle.Data()));
88cb7938 723
724 Int_t p= -1;
725 Float_t par = 0;
e3817e5f 726
727 if (param.Contains("a")) p=4+offset;
728 else if(param.Contains("b")) p=5+offset;
729 else if(param.Contains("c")) p=6+offset;
730 else if(param.Contains("x0"))p=7+offset;
731 else if(param.Contains("y0"))p=8+offset;
12022e83 732
351dd634 733 if (i>4 || i<0) {
734 AliError(Form("No parameter with index %d", i)) ;
735 } else if (p==-1) {
736 AliError(Form("No parameter with name %s", param.Data() )) ;
737 } else
88cb7938 738 par = (*fParameters)(p,i) ;
739
740 return par;
12022e83 741}
742
12022e83 743
26aa7e4a 744//DP____________________________________________________________________________
745//Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const
746//{
747// // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
748//
7fb9892d 749// AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance();
26aa7e4a 750// TVector3 vecEmc ;
751// TVector3 vecCpv ;
752// if(cpv){
753// emc->GetLocalPosition(vecEmc) ;
754// cpv->GetLocalPosition(vecCpv) ;
755//
756// if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
757// // Correct to difference in CPV and EMC position due to different distance to center.
758// // we assume, that particle moves from center
759// Float_t dCPV = geom->GetIPtoOuterCoverDistance();
760// Float_t dEMC = geom->GetIPtoCrystalSurface() ;
761// dEMC = dEMC / dCPV ;
762// vecCpv = dEMC * vecCpv - vecEmc ;
763// if (axis == "X") return vecCpv.X();
764// if (axis == "Y") return vecCpv.Y();
765// if (axis == "Z") return vecCpv.Z();
766// if (axis == "R") return vecCpv.Mag();
767// }
768// return 100000000 ;
769// }
770// return 100000000 ;
771//}
12022e83 772//____________________________________________________________________________
26aa7e4a 773Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const
bc0c084c 774{
c947e71a 775 //Calculates the pid bit for the CPV selection per each purity.
e3817e5f 776 if(effPur>2 || effPur<0)
351dd634 777 AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
26aa7e4a 778
407d15b3 779//DP if(ts->GetCpvIndex()<0)
780//DP return 1 ; //no CPV cluster
bc0c084c 781
e3817e5f 782 Float_t sigX = GetCpv2EmcDistanceCut("X",e);
783 Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
bc0c084c 784
26aa7e4a 785 Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X"));
786 Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z"));
407d15b3 787// Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ;
7fb46731 788
789 //if(deltaX>sigX*(effPur+1))
790 //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1)))
791 if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1)))
bc0c084c 792 return 1;//Neutral
793 else
794 return 0;//Charged
bc0c084c 795}
69183710 796
6ad0bfa0 797//____________________________________________________________________________
fc7e2f43 798Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const
148b2bba 799{
50739f15 800 //Is the particle inside de PCA ellipse?
581354c5 801
e3817e5f 802 particle.ToLower();
803 Int_t prinbit = 0 ;
7fb46731 804 Float_t a = GetEllipseParameter(particle,"a" , e);
805 Float_t b = GetEllipseParameter(particle,"b" , e);
806 Float_t c = GetEllipseParameter(particle,"c" , e);
e3817e5f 807 Float_t x0 = GetEllipseParameter(particle,"x0", e);
808 Float_t y0 = GetEllipseParameter(particle,"y0", e);
809
810 Float_t r = TMath::Power((p[0] - x0)/a,2) +
811 TMath::Power((p[1] - y0)/b,2) +
812 c*(p[0] - x0)*(p[1] - y0)/(a*b) ;
50739f15 813 //3 different ellipses defined
e3817e5f 814 if((effPur==2) && (r<1./2.)) prinbit= 1;
815 if((effPur==1) && (r<2. )) prinbit= 1;
816 if((effPur==0) && (r<9./2.)) prinbit= 1;
50739f15 817
581354c5 818 if(r<0)
351dd634 819 AliError("Negative square?") ;
1f0e7ccd 820
821 return prinbit;
148b2bba 822
148b2bba 823}
1f0e7ccd 824//____________________________________________________________________________
fc7e2f43 825Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
1f0e7ccd 826{
e3817e5f 827 // Set bit for identified hard photons (E > 30 GeV)
828 // if the second moment M2x is below the boundary
829
830 Float_t e = emc->GetEnergy();
831 if (e < 30.0) return 0;
832 Float_t m2x = emc->GetM2x();
833 Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
834 TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
835 TMath::Power(GetParameterPhotonBoundary(2),2)) +
836 GetParameterPhotonBoundary(3);
a675b8d6 837 AliDebug(1, Form("E=%f, m2x=%f, boundary=%f", e,m2x,m2xBoundary));
e3817e5f 838 if (m2x < m2xBoundary)
839 return 1;// A hard photon
840 else
841 return 0;// Not a hard photon
1f0e7ccd 842}
92f521a9 843
e3817e5f 844//____________________________________________________________________________
fc7e2f43 845Int_t AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const
e3817e5f 846{
847 // Set bit for identified hard pi0 (E > 30 GeV)
848 // if the second moment M2x is above the boundary
849
850 Float_t e = emc->GetEnergy();
851 if (e < 30.0) return 0;
852 Float_t m2x = emc->GetM2x();
853 Float_t m2xBoundary = GetParameterPi0Boundary(0) +
854 e * GetParameterPi0Boundary(1);
351dd634 855 AliDebug(1,Form("E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary));
e3817e5f 856 if (m2x > m2xBoundary)
857 return 1;// A hard pi0
bc0c084c 858 else
e3817e5f 859 return 0;// Not a hard pi0
f0a4c9e9 860}
e3817e5f 861
862//____________________________________________________________________________
8d4608b5 863TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSCpvRecPoint * )const
88cb7938 864{
865 // Calculates the momentum direction:
866 // 1. if only a EMC RecPoint, direction is given by IP and this RecPoint
867 // 2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints
868 // However because of the poor position resolution of PPSD the direction is always taken as if we were
869 // in case 1.
f0a4c9e9 870
407d15b3 871 TVector3 local ;
872 emc->GetLocalPosition(local) ;
873
874 AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
875 //Correct for the non-perpendicular incidence
876 // Correction for the depth of the shower starting point (TDR p 127)
877 Float_t para = 0.925 ;
878 Float_t parb = 6.52 ;
879
880 //Remove Old correction (vertex at 0,0,0)
881 TVector3 vtxOld(0.,0.,0.) ;
882 TVector3 vInc ;
883 Float_t x=local.X() ;
884 Float_t z=local.Z() ;
885 phosgeom->GetIncidentVector(vtxOld,emc->GetPHOSMod(),x,z,vInc) ;
886 Float_t depthxOld = 0.;
887 Float_t depthzOld = 0.;
888 Float_t energy = emc->GetEnergy() ;
889 if (energy > 0 && vInc.Y()!=0.) {
753b19cd 890 depthxOld = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
891 depthzOld = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
7fb46731 892 }
adcca1e6 893 else{
407d15b3 894 AliError("Cluster with zero energy \n");
2c06dc7a 895 }
407d15b3 896 //Apply Real vertex
897 phosgeom->GetIncidentVector(fVtx,emc->GetPHOSMod(),x,z,vInc) ;
898 Float_t depthx = 0.;
899 Float_t depthz = 0.;
900 if (energy > 0 && vInc.Y()!=0.) {
e5b7b511 901 depthx = ( para * TMath::Log(energy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
902 depthz = ( para * TMath::Log(energy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
adcca1e6 903 }
407d15b3 904
e5b7b511 905 //Correct for the vertex position and shower depth
906 Double_t xd=x+(depthxOld-depthx) ;
907 Double_t zd=z+(depthzOld-depthz) ;
908 TVector3 dir(0,0,0) ;
909 phosgeom->Local2Global(emc->GetPHOSMod(),xd,zd,dir) ;
407d15b3 910
e5b7b511 911 dir-=fVtx ;
7fb46731 912 dir.SetMag(1.) ;
e3817e5f 913
88cb7938 914 return dir ;
7b7c1533 915}
916
35adb638 917//________________________________________________________________________
17323043 918Double_t AliPHOSPIDv1::LandauF(Double_t x, Double_t y, Double_t * par)
35adb638 919{
c947e71a 920 //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
921 //this method returns a density probability of this parameter, given by a landau
922 //function whose parameters depend with the energy with a function: a/(x*x)+b/x+b
923
2924941c 924 if (x > par[9]) x = par[9];
925
926 //Double_t cnt = par[1] / (x*x) + par[2] / x + par[0] ;
927 Double_t cnt = par[0] + par[1] * x + par[2] * x * x ;
35adb638 928 Double_t mean = par[4] / (x*x) + par[5] / x + par[3] ;
929 Double_t sigma = par[7] / (x*x) + par[8] / x + par[6] ;
c947e71a 930
931 if(TMath::Abs(sigma) > 1.e-10){
acb5beb7 932 return cnt*TMath::Landau(y,mean,sigma);
35adb638 933 }
934 else
935 return 0.;
936
937}
938//________________________________________________________________________
17323043 939Double_t AliPHOSPIDv1::LandauPol2(Double_t x, Double_t y, Double_t * par)
35adb638 940{
c947e71a 941
942 //Given the energy x and the parameter y (tof, shower dispersion or cpv-emc distance),
943 //this method returns a density probability of this parameter, given by a landau
944 //function whose parameters depend with the energy like second order polinomial
945
35adb638 946 Double_t cnt = par[2] * (x*x) + par[1] * x + par[0] ;
c947e71a 947 Double_t mean = par[5] * (x*x) + par[4] * x + par[3] ;
948 Double_t sigma = par[8] * (x*x) + par[7] * x + par[6] ;
35adb638 949
c947e71a 950 if(TMath::Abs(sigma) > 1.e-10){
acb5beb7 951 return cnt*TMath::Landau(y,mean,sigma);
35adb638 952 }
953 else
954 return 0.;
c947e71a 955
956
35adb638 957}
958// //________________________________________________________________________
959// Double_t AliPHOSPIDv1::ChargedHadronDistProb(Double_t x, Double_t y, Double_t * parg, Double_t * parl)
960// {
961// Double_t cnt = 0.0 ;
962// Double_t mean = 0.0 ;
963// Double_t sigma = 0.0 ;
964// Double_t arg = 0.0 ;
965// if (y < parl[4] / (x*x) + parl[5] / x + parl[3]){
966// cnt = parg[1] / (x*x) + parg[2] / x + parg[0] ;
967// mean = parg[4] / (x*x) + parg[5] / x + parg[3] ;
968// sigma = parg[7] / (x*x) + parg[8] / x + parg[6] ;
969// TF1 * f = new TF1("gaus","gaus",0.,100.);
970// f->SetParameters(cnt,mean,sigma);
971// arg = f->Eval(y) ;
972// }
973// else{
974// cnt = parl[1] / (x*x) + parl[2] / x + parl[0] ;
975// mean = parl[4] / (x*x) + parl[5] / x + parl[3] ;
976// sigma = parl[7] / (x*x) + parl[8] / x + parl[6] ;
977// TF1 * f = new TF1("landau","landau",0.,100.);
978// f->SetParameters(cnt,mean,sigma);
979// arg = f->Eval(y) ;
980// }
981// // Double_t mean = par[3] + par[4] * x + par[5] * x * x ;
982// // Double_t sigma = par[6] + par[7] * x + par[8] * x * x ;
983
984// //Double_t arg = -(y-mean)*(y-mean)/(2*sigma*sigma) ;
985// //return cnt * TMath::Exp(arg) ;
986
987// return arg;
988
989// }
2cc71c1e 990//____________________________________________________________________________
991void AliPHOSPIDv1::MakePID()
992{
993 // construct the PID weight from a Bayesian Method
c947e71a 994
304864ab 995 const Int_t kSPECIES = AliPID::kSPECIESN ;
7fb46731 996
9a2cdbdf 997 Int_t nparticles = fRecParticles->GetEntriesFast() ;
7fb46731 998
9a2cdbdf 999 if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
7fb46731 1000 AliFatal("RecPoints or TrackSegments not found !") ;
1001 }
9a2cdbdf 1002
1003 TIter next(fTrackSegments) ;
7fb46731 1004 AliPHOSTrackSegment * ts ;
1005 Int_t index = 0 ;
1006
35adb638 1007 Double_t * stof[kSPECIES] ;
1008 Double_t * sdp [kSPECIES] ;
1009 Double_t * scpv[kSPECIES] ;
fb7b51ad 1010 Double_t * sw [kSPECIES] ;
35adb638 1011 //Info("MakePID","Begin MakePID");
1012
1013 for (Int_t i =0; i< kSPECIES; i++){
1014 stof[i] = new Double_t[nparticles] ;
1015 sdp [i] = new Double_t[nparticles] ;
1016 scpv[i] = new Double_t[nparticles] ;
fb7b51ad 1017 sw [i] = new Double_t[nparticles] ;
35adb638 1018 }
1019
7fb46731 1020
1021 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
1022
1023 //cout<<">>>>>> Bayesian Index "<<index<<endl;
1024
1025 AliPHOSEmcRecPoint * emc = 0 ;
1026 if(ts->GetEmcIndex()>=0)
9a2cdbdf 1027 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
7fb46731 1028
407d15b3 1029// AliPHOSCpvRecPoint * cpv = 0 ;
1030// if(ts->GetCpvIndex()>=0)
1031// cpv = (AliPHOSCpvRecPoint *) cpvRecPoints->At(ts->GetCpvIndex()) ;
1032//
1033//// Int_t track = 0 ;
1034//// track = ts->GetTrackIndex() ; //TPC tracks ?
cc1fe362 1035
7fb46731 1036 if (!emc) {
a675b8d6 1037 AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
7fb46731 1038 }
c947e71a 1039
2924941c 1040
7fb46731 1041 // ############Tof#############################
c947e71a 1042
7fb46731 1043 // Info("MakePID", "TOF");
fb7b51ad 1044 Float_t en = emc->GetEnergy();
7fb46731 1045 Double_t time = emc->GetTime() ;
1046 // cout<<">>>>>>>Energy "<<en<<"Time "<<time<<endl;
fb7b51ad 1047
35adb638 1048 // now get the signals probability
1049 // s(pid) in the Bayesian formulation
cc1fe362 1050
304864ab 1051 stof[AliPID::kPhoton][index] = 1.;
1052 stof[AliPID::kElectron][index] = 1.;
304864ab 1053 stof[AliPID::kEleCon][index] = 1.;
7fb46731 1054 //We assing the same prob to charged hadrons, sum is 1
1055 stof[AliPID::kPion][index] = 1./3.;
1056 stof[AliPID::kKaon][index] = 1./3.;
1057 stof[AliPID::kProton][index] = 1./3.;
1058 //We assing the same prob to neutral hadrons, sum is 1
1059 stof[AliPID::kNeutron][index] = 1./2.;
1060 stof[AliPID::kKaon0][index] = 1./2.;
304864ab 1061 stof[AliPID::kMuon][index] = 1.;
2924941c 1062
1063 if(en < fTOFEnThreshold) {
7fb46731 1064
1065 Double_t pTofPion = fTFpiong ->Eval(time) ; //gaus distribution
1066 Double_t pTofKaon = 0;
1067
35adb638 1068 if(time < fTkaonl[1])
fb7b51ad 1069 pTofKaon = fTFkaong ->Eval(time) ; //gaus distribution
35adb638 1070 else
fb7b51ad 1071 pTofKaon = fTFkaonl ->Eval(time) ; //landau distribution
7fb46731 1072
1073 Double_t pTofNucleon = 0;
1074
35adb638 1075 if(time < fThhadronl[1])
7fb46731 1076 pTofNucleon = fTFhhadrong ->Eval(time) ; //gaus distribution
35adb638 1077 else
7fb46731 1078 pTofNucleon = fTFhhadronl ->Eval(time) ; //landau distribution
7fb46731 1079 //We assing the same prob to neutral hadrons, sum is the average prob
fb7b51ad 1080 Double_t pTofNeHadron = (pTofKaon + pTofNucleon)/2. ;
1081 //We assing the same prob to charged hadrons, sum is the average prob
1082 Double_t pTofChHadron = (pTofPion + pTofKaon + pTofNucleon)/3. ;
7fb46731 1083
fb7b51ad 1084 stof[AliPID::kPhoton][index] = fTFphoton ->Eval(time) ;
1085 //gaus distribution
1086 stof[AliPID::kEleCon][index] = stof[AliPID::kPhoton][index] ;
1087 //a conversion electron has the photon ToF
304864ab 1088 stof[AliPID::kMuon][index] = stof[AliPID::kPhoton][index] ;
7fb46731 1089
fb7b51ad 1090 stof[AliPID::kElectron][index] = pTofPion ;
7fb46731 1091
1092 stof[AliPID::kPion][index] = pTofChHadron ;
1093 stof[AliPID::kKaon][index] = pTofChHadron ;
1094 stof[AliPID::kProton][index] = pTofChHadron ;
1095
1096 stof[AliPID::kKaon0][index] = pTofNeHadron ;
1097 stof[AliPID::kNeutron][index] = pTofNeHadron ;
cc1fe362 1098 }
c947e71a 1099
1100 // Info("MakePID", "Dispersion");
cc1fe362 1101
7fb46731 1102 // ###########Shower shape: Dispersion####################
35adb638 1103 Float_t dispersion = emc->GetDispersion();
407d15b3 1104 //DP: Correct for non-perpendicular incidence
1105 //DP: still to be done
1106
35adb638 1107 //dispersion is not well defined if the cluster is only in few crystals
cc1fe362 1108
304864ab 1109 sdp[AliPID::kPhoton][index] = 1. ;
1110 sdp[AliPID::kElectron][index] = 1. ;
1111 sdp[AliPID::kPion][index] = 1. ;
1112 sdp[AliPID::kKaon][index] = 1. ;
1113 sdp[AliPID::kProton][index] = 1. ;
1114 sdp[AliPID::kNeutron][index] = 1. ;
1115 sdp[AliPID::kEleCon][index] = 1. ;
1116 sdp[AliPID::kKaon0][index] = 1. ;
1117 sdp[AliPID::kMuon][index] = 1. ;
c947e71a 1118
2924941c 1119 if(en > fDispEnThreshold && emc->GetMultiplicity() > fDispMultThreshold){
7fb46731 1120 sdp[AliPID::kPhoton][index] = GausF(en , dispersion, fDphoton) ;
304864ab 1121 sdp[AliPID::kElectron][index] = sdp[AliPID::kPhoton][index] ;
1122 sdp[AliPID::kPion][index] = LandauF(en , dispersion, fDhadron ) ;
1123 sdp[AliPID::kKaon][index] = sdp[AliPID::kPion][index] ;
1124 sdp[AliPID::kProton][index] = sdp[AliPID::kPion][index] ;
1125 sdp[AliPID::kNeutron][index] = sdp[AliPID::kPion][index] ;
1126 sdp[AliPID::kEleCon][index] = sdp[AliPID::kPhoton][index];
1127 sdp[AliPID::kKaon0][index] = sdp[AliPID::kPion][index] ;
fb7b51ad 1128 sdp[AliPID::kMuon][index] = fDFmuon ->Eval(dispersion) ;
1129 //landau distribution
35adb638 1130 }
cc1fe362 1131
7fb46731 1132// Info("MakePID","multiplicity %d, dispersion %f", emc->GetMultiplicity(), dispersion);
1133// Info("MakePID","ss: photon %f, hadron %f ", sdp[AliPID::kPhoton][index], sdp[AliPID::kPion][index]);
c947e71a 1134// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<", dispersion "<< dispersion<<endl ;
304864ab 1135// cout<<"<<<<<ss: photon "<<sdp[AliPID::kPhoton][index]<<", hadron "<<sdp[AliPID::kPion][index]<<endl;
c947e71a 1136
7fb46731 1137 //########## CPV-EMC Distance#######################
1138 // Info("MakePID", "Distance");
fb7b51ad 1139
26aa7e4a 1140 Float_t x = TMath::Abs(ts->GetCpvDistance("X")) ;
1141 Float_t z = ts->GetCpvDistance("Z") ;
fb7b51ad 1142
7fb46731 1143 Double_t pcpv = 0 ;
c947e71a 1144 Double_t pcpvneutral = 0. ;
7fb46731 1145
1146 Double_t elprobx = GausF(en , x, fXelectron) ;
1147 Double_t elprobz = GausF(en , z, fZelectron) ;
1148 Double_t chprobx = GausF(en , x, fXcharged) ;
1149 Double_t chprobz = GausF(en , z, fZcharged) ;
c947e71a 1150 Double_t pcpvelectron = elprobx * elprobz;
1151 Double_t pcpvcharged = chprobx * chprobz;
7fb46731 1152
1153// cout<<">>>>energy "<<en<<endl;
c947e71a 1154// cout<<">>>>electron : x "<<x<<" xprob "<<elprobx<<" z "<<z<<" zprob "<<elprobz<<endl;
1155// cout<<">>>>hadron : x "<<x<<" xprob "<<chprobx<<" z "<<z<<" zprob "<<chprobz<<endl;
1156// cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
1157
2924941c 1158 // Is neutral or charged?
35adb638 1159 if(pcpvelectron >= pcpvcharged)
1160 pcpv = pcpvelectron ;
1161 else
1162 pcpv = pcpvcharged ;
1163
fb7b51ad 1164 if(pcpv < fChargedNeutralThreshold)
35adb638 1165 {
1166 pcpvneutral = 1. ;
1167 pcpvcharged = 0. ;
1168 pcpvelectron = 0. ;
1169 }
c947e71a 1170 // else
1171 // cout<<">>>>>>>>>>>CHARGED>>>>>>>>>>>"<<endl;
35adb638 1172
304864ab 1173 scpv[AliPID::kPion][index] = pcpvcharged ;
1174 scpv[AliPID::kKaon][index] = pcpvcharged ;
1175 scpv[AliPID::kProton][index] = pcpvcharged ;
7fb46731 1176
1177 scpv[AliPID::kMuon][index] = pcpvelectron ;
304864ab 1178 scpv[AliPID::kElectron][index] = pcpvelectron ;
304864ab 1179 scpv[AliPID::kEleCon][index] = pcpvelectron ;
7fb46731 1180
1181 scpv[AliPID::kPhoton][index] = pcpvneutral ;
1182 scpv[AliPID::kNeutron][index] = pcpvneutral ;
304864ab 1183 scpv[AliPID::kKaon0][index] = pcpvneutral ;
7fb46731 1184
35adb638 1185
1186 // Info("MakePID", "CPV passed");
c947e71a 1187
7fb46731 1188 //############## Pi0 #############################
304864ab 1189 stof[AliPID::kPi0][index] = 0. ;
1190 scpv[AliPID::kPi0][index] = 0. ;
1191 sdp [AliPID::kPi0][index] = 0. ;
c947e71a 1192
35adb638 1193 if(en > 30.){
1194 // pi0 are detected via decay photon
2924941c 1195 stof[AliPID::kPi0][index] = stof[AliPID::kPhoton][index];
304864ab 1196 scpv[AliPID::kPi0][index] = pcpvneutral ;
2924941c 1197 if(emc->GetMultiplicity() > fDispMultThreshold)
1198 sdp [AliPID::kPi0][index] = GausF(en , dispersion, fDpi0) ;
1199 //sdp [AliPID::kPi0][index] = GausPol2(en , dispersion, fDpi0) ;
1200// cout<<"E = "<<en<<" GeV; disp = "<<dispersion<<"; mult = "
1201// <<emc->GetMultiplicity()<<endl;
1202// cout<<"PDF: photon = "<<sdp [AliPID::kPhoton][index]<<"; pi0 = "
1203// <<sdp [AliPID::kPi0][index]<<endl;
35adb638 1204 }
1205
2924941c 1206
7fb46731 1207
2924941c 1208
7fb46731 1209 //############## muon #############################
1210
35adb638 1211 if(en > 0.5){
1212 //Muons deposit few energy
304864ab 1213 scpv[AliPID::kMuon][index] = 0 ;
1214 stof[AliPID::kMuon][index] = 0 ;
1215 sdp [AliPID::kMuon][index] = 0 ;
c947e71a 1216 }
1217
fb7b51ad 1218 //Weight to apply to hadrons due to energy reconstruction
1219
1220 Float_t weight = fERecWeight ->Eval(en) ;
1221
1222 sw[AliPID::kPhoton][index] = 1. ;
1223 sw[AliPID::kElectron][index] = 1. ;
1224 sw[AliPID::kPion][index] = weight ;
1225 sw[AliPID::kKaon][index] = weight ;
1226 sw[AliPID::kProton][index] = weight ;
1227 sw[AliPID::kNeutron][index] = weight ;
1228 sw[AliPID::kEleCon][index] = 1. ;
1229 sw[AliPID::kKaon0][index] = weight ;
1230 sw[AliPID::kMuon][index] = weight ;
1231 sw[AliPID::kPi0][index] = 1. ;
1232
7fb46731 1233// if(en > 0.5){
1234// cout<<"######################################################"<<endl;
1235// //cout<<"MakePID: energy "<<en<<", tof "<<time<<", distance "<<distance<<", dispersion "<<dispersion<<endl ;
1236// cout<<"MakePID: energy "<<en<<", tof "<<time<<", dispersion "<<dispersion<<", x "<<x<<", z "<<z<<endl ;
1237// cout<<">>>>>multiplicity "<<emc->GetMultiplicity()<<endl;
1238// cout<<">>>>electron : xprob "<<elprobx<<" zprob "<<elprobz<<endl;
1239// cout<<">>>>hadron : xprob "<<chprobx<<" zprob "<<chprobz<<endl;
1240// cout<<">>>>electron : px*pz "<<pcpvelectron <<" hadron: px*pz "<<pcpvcharged<<endl;
c947e71a 1241
2924941c 1242// cout<<"Photon , pid "<< fInitPID[AliPID::kPhoton]<<" tof "<<stof[AliPID::kPhoton][index]
1243// <<", cpv "<<scpv[AliPID::kPhoton][index]<<", ss "<<sdp[AliPID::kPhoton][index]<<endl;
7fb46731 1244// cout<<"EleCon , pid "<< fInitPID[AliPID::kEleCon]<<", tof "<<stof[AliPID::kEleCon][index]
1245// <<", cpv "<<scpv[AliPID::kEleCon][index]<<" ss "<<sdp[AliPID::kEleCon][index]<<endl;
1246// cout<<"Electron , pid "<< fInitPID[AliPID::kElectron]<<", tof "<<stof[AliPID::kElectron][index]
1247// <<", cpv "<<scpv[AliPID::kElectron][index]<<" ss "<<sdp[AliPID::kElectron][index]<<endl;
1248// cout<<"Muon , pid "<< fInitPID[AliPID::kMuon]<<", tof "<<stof[AliPID::kMuon][index]
1249// <<", cpv "<<scpv[AliPID::kMuon][index]<<" ss "<<sdp[AliPID::kMuon][index]<<endl;
2924941c 1250// cout<<"Pi0 , pid "<< fInitPID[AliPID::kPi0]<<", tof "<<stof[AliPID::kPi0][index]
1251// <<", cpv "<<scpv[AliPID::kPi0][index]<<" ss "<<sdp[AliPID::kPi0][index]<<endl;
7fb46731 1252// cout<<"Pion , pid "<< fInitPID[AliPID::kPion]<<", tof "<<stof[AliPID::kPion][index]
1253// <<", cpv "<<scpv[AliPID::kPion][index]<<" ss "<<sdp[AliPID::kPion][index]<<endl;
1254// cout<<"Kaon0 , pid "<< fInitPID[AliPID::kKaon0]<<", tof "<<stof[AliPID::kKaon0][index]
1255// <<", cpv "<<scpv[AliPID::kKaon0][index]<<" ss "<<sdp[AliPID::kKaon0][index]<<endl;
1256// cout<<"Kaon , pid "<< fInitPID[AliPID::kKaon]<<", tof "<<stof[AliPID::kKaon][index]
1257// <<", cpv "<<scpv[AliPID::kKaon][index]<<" ss "<<sdp[AliPID::kKaon][index]<<endl;
1258// cout<<"Neutron , pid "<< fInitPID[AliPID::kNeutron]<<", tof "<<stof[AliPID::kNeutron][index]
1259// <<", cpv "<<scpv[AliPID::kNeutron][index]<<" ss "<<sdp[AliPID::kNeutron][index]<<endl;
1260// cout<<"Proton , pid "<< fInitPID[AliPID::kProton]<<", tof "<<stof[AliPID::kProton][index]
1261// <<", cpv "<<scpv[AliPID::kProton][index]<<" ss "<<sdp[AliPID::kProton][index]<<endl;
1262// cout<<"######################################################"<<endl;
1263// }
1264 index++;
cc1fe362 1265 }
35adb638 1266
1267 //for (index = 0 ; index < kSPECIES ; index++)
1268 // pid[index] /= nparticles ;
1269
7fb46731 1270
35adb638 1271 // Info("MakePID", "Total Probability calculation");
1272
cc1fe362 1273 for(index = 0 ; index < nparticles ; index ++) {
2924941c 1274
9a2cdbdf 1275 AliPHOSRecParticle * recpar = static_cast<AliPHOSRecParticle *>(fRecParticles->At(index));
2924941c 1276
1277 //Conversion electron?
1278
1279 if(recpar->IsEleCon()){
1280 fInitPID[AliPID::kEleCon] = 1. ;
1281 fInitPID[AliPID::kPhoton] = 0. ;
1282 fInitPID[AliPID::kElectron] = 0. ;
1283 }
1284 else{
1285 fInitPID[AliPID::kEleCon] = 0. ;
1286 fInitPID[AliPID::kPhoton] = 1. ;
1287 fInitPID[AliPID::kElectron] = 1. ;
1288 }
1289 // fInitPID[AliPID::kEleCon] = 0. ;
1290
1291
cc1fe362 1292 // calculates the Bayesian weight
7fb46731 1293
cc1fe362 1294 Int_t jndex ;
1295 Double_t wn = 0.0 ;
1296 for (jndex = 0 ; jndex < kSPECIES ; jndex++)
2924941c 1297 wn += stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1298 sw[jndex][index] * fInitPID[jndex] ;
1299
7fb46731 1300 // cout<<"*************wn "<<wn<<endl;
e74ea0e9 1301 if (TMath::Abs(wn)>0)
1302 for (jndex = 0 ; jndex < kSPECIES ; jndex++) {
35adb638 1303 //cout<<"jndex "<<jndex<<" wn "<<wn<<" SetPID * wn"
1304 //<<stof[jndex][index] * sdp[jndex][index] * pid[jndex] << endl;
1305 //cout<<" tof "<<stof[jndex][index] << " disp " <<sdp[jndex][index] << " pid "<< fInitPID[jndex] << endl;
2924941c 1306 // if(jndex == AliPID::kPi0 || jndex == AliPID::kPhoton){
1307 // cout<<"Particle "<<jndex<<" final prob * wn "
1308 // <<stof[jndex][index] * sdp[jndex][index] * scpv[jndex][index] *
1309 // fInitPID[jndex] <<" wn "<< wn<<endl;
1310 // cout<<"pid "<< fInitPID[jndex]<<", tof "<<stof[jndex][index]
1311 // <<", cpv "<<scpv[jndex][index]<<" ss "<<sdp[jndex][index]<<endl;
1312 // }
1313 recpar->SetPID(jndex, stof[jndex][index] * sdp[jndex][index] *
1314 sw[jndex][index] * scpv[jndex][index] *
1315 fInitPID[jndex] / wn) ;
e74ea0e9 1316 }
2cc71c1e 1317 }
35adb638 1318 // Info("MakePID", "Delete");
1319
2924941c 1320 for (Int_t i =0; i< kSPECIES; i++){
1321 delete [] stof[i];
1322 delete [] sdp [i];
1323 delete [] scpv[i];
1324 delete [] sw [i];
1325 }
35adb638 1326 // Info("MakePID","End MakePID");
2cc71c1e 1327}
1328
7acf6008 1329//____________________________________________________________________________
e3817e5f 1330void AliPHOSPIDv1::MakeRecParticles()
1331{
b2a60966 1332 // Makes a RecParticle out of a TrackSegment
148b2bba 1333
9a2cdbdf 1334 if ( !fEMCRecPoints || !fCPVRecPoints || !fTrackSegments ) {
351dd634 1335 AliFatal("RecPoints or TrackSegments not found !") ;
148b2bba 1336 }
9a2cdbdf 1337 fRecParticles->Clear();
148b2bba 1338
9a2cdbdf 1339 TIter next(fTrackSegments) ;
7acf6008 1340 AliPHOSTrackSegment * ts ;
6ad0bfa0 1341 Int_t index = 0 ;
09fc14a0 1342 AliPHOSRecParticle * rp ;
7acf6008 1343 while ( (ts = (AliPHOSTrackSegment *)next()) ) {
7fb46731 1344 // cout<<">>>>>>>>>>>>>>>PCA Index "<<index<<endl;
9a2cdbdf 1345 new( (*fRecParticles)[index] ) AliPHOSRecParticle() ;
1346 rp = (AliPHOSRecParticle *)fRecParticles->At(index) ;
f0a4c9e9 1347 rp->SetTrackSegment(index) ;
9688c1dd 1348 rp->SetIndexInList(index) ;
148b2bba 1349
7acf6008 1350 AliPHOSEmcRecPoint * emc = 0 ;
1351 if(ts->GetEmcIndex()>=0)
9a2cdbdf 1352 emc = (AliPHOSEmcRecPoint *) fEMCRecPoints->At(ts->GetEmcIndex()) ;
fad3e5b9 1353
8d4608b5 1354 AliPHOSCpvRecPoint * cpv = 0 ;
7acf6008 1355 if(ts->GetCpvIndex()>=0)
9a2cdbdf 1356 cpv = (AliPHOSCpvRecPoint *) fCPVRecPoints->At(ts->GetCpvIndex()) ;
fad3e5b9 1357
bd76890a 1358 Int_t track = 0 ;
1359 track = ts->GetTrackIndex() ;
1360
148b2bba 1361 // Now set type (reconstructed) of the particle
1362
1363 // Choose the cluster energy range
9fa5f1d0 1364
fbf811ec 1365 if (!emc) {
a675b8d6 1366 AliFatal(Form("-> emc(%d)", ts->GetEmcIndex())) ;
fbf811ec 1367 }
50739f15 1368
e3817e5f 1369 Float_t e = emc->GetEnergy() ;
bc0c084c 1370
6f969528 1371 Float_t lambda[2] ;
1372 emc->GetElipsAxis(lambda) ;
407d15b3 1373
50739f15 1374 if((lambda[0]>0.01) && (lambda[1]>0.01)){
1375 // Looking PCA. Define and calculate the data (X),
bc0c084c 1376 // introduce in the function X2P that gives the components (P).
1377
c947e71a 1378 Float_t spher = 0. ;
1379 Float_t emaxdtotal = 0. ;
50739f15 1380
bc0c084c 1381 if((lambda[0]+lambda[1])!=0)
e4df4b30 1382 spher=TMath::Abs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]);
50739f15 1383
c947e71a 1384 emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy();
50739f15 1385
1386 fX[0] = lambda[0] ;
1387 fX[1] = lambda[1] ;
1388 fX[2] = emc->GetDispersion() ;
c947e71a 1389 fX[3] = spher ;
50739f15 1390 fX[4] = emc->GetMultiplicity() ;
c947e71a 1391 fX[5] = emaxdtotal ;
50739f15 1392 fX[6] = emc->GetCoreEnergy() ;
1393
e3817e5f 1394 fPrincipalPhoton->X2P(fX,fPPhoton);
1395 fPrincipalPi0 ->X2P(fX,fPPi0);
1f0e7ccd 1396
50739f15 1397 }
1398 else{
e3817e5f 1399 fPPhoton[0]=-100.0; //We do not accept clusters with
1400 fPPhoton[1]=-100.0; //one cell as a photon-like
1401 fPPi0[0] =-100.0;
1402 fPPi0[1] =-100.0;
50739f15 1403 }
1404
2cc71c1e 1405 Float_t time = emc->GetTime() ;
1406 rp->SetTof(time) ;
9fa5f1d0 1407
bc0c084c 1408 // Loop of Efficiency-Purity (the 3 points of purity or efficiency
1409 // are taken into account to set the particle identification)
e3817e5f 1410 for(Int_t effPur = 0; effPur < 3 ; effPur++){
50739f15 1411
bc0c084c 1412 // Looking at the CPV detector. If RCPV greater than CpvEmcDistance,
1413 // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
1414 // is set to 1
26aa7e4a 1415 if(GetCPVBit(ts, effPur,e) == 1 ){
e3817e5f 1416 rp->SetPIDBit(effPur) ;
35adb638 1417 //cout<<"CPV bit "<<effPur<<endl;
1418 }
50739f15 1419 // Looking the TOF. If TOF smaller than gate, 4th, 5th or 6th
1420 // bit (depending on the efficiency-purity point )is set to 1
2cc71c1e 1421 if(time< (*fParameters)(3,effPur))
e3817e5f 1422 rp->SetPIDBit(effPur+3) ;
2cc71c1e 1423
e3817e5f 1424 //Photon PCA
50739f15 1425 //If we are inside the ellipse, 7th, 8th or 9th
1426 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 1427 if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1)
1428 rp->SetPIDBit(effPur+6) ;
1f0e7ccd 1429
e3817e5f 1430 //Pi0 PCA
1f0e7ccd 1431 //If we are inside the ellipse, 10th, 11th or 12th
1432 // bit (depending on the efficiency-purity point )is set to 1
e3817e5f 1433 if(GetPrincipalBit("pi0" ,fPPi0 ,effPur,e) == 1)
1434 rp->SetPIDBit(effPur+9) ;
f0a4c9e9 1435 }
e3817e5f 1436 if(GetHardPhotonBit(emc))
1437 rp->SetPIDBit(12) ;
1438 if(GetHardPi0Bit (emc))
1439 rp->SetPIDBit(13) ;
1f0e7ccd 1440
bd76890a 1441 if(track >= 0)
1442 rp->SetPIDBit(14) ;
1443
9fa5f1d0 1444 //Set momentum, energy and other parameters
50739f15 1445 Float_t encal = GetCalibratedEnergy(e);
9fa5f1d0 1446 TVector3 dir = GetMomentumDirection(emc,cpv) ;
1447 dir.SetMag(encal) ;
1448 rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
1449 rp->SetCalcMass(0);
e0ed2e49 1450 rp->Name(); //If photon sets the particle pdg name to gamma
407d15b3 1451 rp->SetProductionVertex(fVtx.X(),fVtx.Y(),fVtx.Z(),0);
e747b8da 1452 rp->SetFirstMother(-1);
1453 rp->SetLastMother(-1);
1454 rp->SetFirstDaughter(-1);
1455 rp->SetLastDaughter(-1);
1456 rp->SetPolarisation(0,0,0);
d956e9b7 1457 //Set the position in global coordinate system from the RecPoint
605cb8bb 1458 AliPHOSTrackSegment * ts1 = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
1459 AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts1->GetEmcIndex()));
d956e9b7 1460 TVector3 pos ;
9a2cdbdf 1461 fGeom->GetGlobalPHOS(erp, pos) ;
d956e9b7 1462 rp->SetPos(pos);
6ad0bfa0 1463 index++ ;
1464 }
6ad0bfa0 1465}
e3817e5f 1466
09fc14a0 1467//____________________________________________________________________________
702ab87e 1468void AliPHOSPIDv1::Print(const Option_t *) const
09fc14a0 1469{
b2a60966 1470 // Print the parameters used for the particle type identification
bc0c084c 1471
351dd634 1472 AliInfo("=============== AliPHOSPIDv1 ================") ;
88cb7938 1473 printf("Making PID\n") ;
1474 printf(" Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() ) ;
1475 printf(" Name of parameters file %s\n", fFileNameParameters.Data() ) ;
1476 printf(" Matrix of Parameters: 14x4\n") ;
1477 printf(" Energy Calibration 1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
1478 printf(" RCPV 2x3 rows x and z, columns function cut parameters\n") ;
1479 printf(" TOF 1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
1480 printf(" PCA 5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
407d15b3 1481 printf(" Pi0 PCA 5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
50739f15 1482 fParameters->Print() ;
09fc14a0 1483}
1484
a496c46c 1485
69183710 1486
7acf6008 1487//____________________________________________________________________________
a4e98857 1488void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
1489{
dd5c4038 1490 // Print table of reconstructed particles
1491
21cd0c07 1492 TString message ;
9a2cdbdf 1493 message = " found " ;
1494 message += fRecParticles->GetEntriesFast();
3bf72d32 1495 message += " RecParticles\n" ;
1496
7acf6008 1497 if(strstr(option,"all")) { // printing found TS
3bf72d32 1498 message += "\n PARTICLE Index \n" ;
7acf6008 1499
1500 Int_t index ;
9a2cdbdf 1501 for (index = 0 ; index < fRecParticles->GetEntries() ; index++) {
1502 AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) fRecParticles->At(index) ;
3bf72d32 1503 message += "\n" ;
1504 message += rp->Name().Data() ;
1505 message += " " ;
1506 message += rp->GetIndexInList() ;
1507 message += " " ;
1508 message += rp->GetType() ;
7acf6008 1509 }
3bf72d32 1510 }
351dd634 1511 AliInfo(message.Data() ) ;
69183710 1512}
88cb7938 1513
1514//____________________________________________________________________________
1515void AliPHOSPIDv1::SetParameters()
1516{
1517 // PCA : To do the Principal Components Analysis it is necessary
1518 // the Principal file, which is opened here
1519 fX = new double[7]; // Data for the PCA
1520 fPPhoton = new double[7]; // Eigenvalues of the PCA
1521 fPPi0 = new double[7]; // Eigenvalues of the Pi0 PCA
1522
1523 // Read photon principals from the photon file
1524
1525 fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
1526 TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
1527 fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ;
1528 f.Close() ;
1529
1530 // Read pi0 principals from the pi0 file
1531
1532 fFileNamePrincipalPi0 = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
1533 TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
1534 fPrincipalPi0 = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ;
1535 fPi0.Close() ;
1536
1537 // Open parameters file and initialization of the Parameters matrix.
1538 // In the File Parameters.dat are all the parameters. These are introduced
1539 // in a matrix of 16x4
1540 //
1541 // All the parameters defined in this file are, in order of row:
1542 // line 0 : calibration
1543 // lines 1,2 : CPV rectangular cat for X and Z
1544 // line 3 : TOF cut
1545 // lines 4-8 : parameters to calculate photon PCA ellipse
1546 // lines 9-13: parameters to calculate pi0 PCA ellipse
1547 // lines 14-15: parameters to calculate border for high-pt photons and pi0
1548
1549 fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
e8d02863 1550 fParameters = new TMatrixF(16,4) ;
c947e71a 1551 const Int_t kMaxLeng=255;
1552 char string[kMaxLeng];
88cb7938 1553
1554 // Open a text file with PID parameters
1555 FILE *fd = fopen(fFileNameParameters.Data(),"r");
1556 if (!fd)
351dd634 1557 AliFatal(Form("File %s with a PID parameters cannot be opened\n",
1558 fFileNameParameters.Data()));
88cb7938 1559
1560 Int_t i=0;
1561 // Read parameter file line-by-line and skip empty line and comments
c947e71a 1562 while (fgets(string,kMaxLeng,fd) != NULL) {
88cb7938 1563 if (string[0] == '\n' ) continue;
1564 if (string[0] == '!' ) continue;
1565 sscanf(string, "%f %f %f %f",
1566 &(*fParameters)(i,0), &(*fParameters)(i,1),
1567 &(*fParameters)(i,2), &(*fParameters)(i,3));
1568 i++;
a675b8d6 1569 AliDebug(1, Form("Line %d: %s",i,string));
88cb7938 1570 }
1571 fclose(fd);
1572}
1573
1574//____________________________________________________________________________
1575void AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param)
1576{
1577 // Set parameter "Calibration" i to a value param
351dd634 1578 if(i>2 || i<0) {
1579 AliError(Form("Invalid parameter number: %d",i));
1580 } else
88cb7938 1581 (*fParameters)(0,i) = param ;
1582}
1583
1584//____________________________________________________________________________
1585void AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut)
1586{
1587 // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on
1588 // Purity-Efficiency point i
1589
351dd634 1590 if(i>2 || i<0) {
1591 AliError(Form("Invalid parameter number: %d",i));
1592 } else {
88cb7938 1593 axis.ToLower();
1594 if (axis == "x") (*fParameters)(1,i) = cut;
1595 else if (axis == "z") (*fParameters)(2,i) = cut;
351dd634 1596 else {
1597 AliError(Form("Invalid axis name: %s",axis.Data()));
1598 }
88cb7938 1599 }
1600}
1601
1602//____________________________________________________________________________
1603void AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param)
1604{
1605 // Set parameter "Hard photon boundary" i to a value param
351dd634 1606 if(i>4 || i<0) {
1607 AliError(Form("Invalid parameter number: %d",i));
1608 } else
88cb7938 1609 (*fParameters)(14,i) = param ;
1610}
1611
1612//____________________________________________________________________________
1613void AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param)
1614{
1615 // Set parameter "Hard pi0 boundary" i to a value param
351dd634 1616 if(i>1 || i<0) {
1617 AliError(Form("Invalid parameter number: %d",i));
1618 } else
88cb7938 1619 (*fParameters)(15,i) = param ;
1620}
1621
1622//_____________________________________________________________________________
1623void AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate)
1624{
1625 // Set the parameter TimeGate depending on Purity-Efficiency point i
351dd634 1626 if (i>2 || i<0) {
1627 AliError(Form("Invalid Efficiency-Purity choice %d",i));
1628 } else
88cb7938 1629 (*fParameters)(3,i)= gate ;
1630}
1631
1632//_____________________________________________________________________________
1633void AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par)
1634{
1635 // Set the parameter "i" that is needed to calculate the ellipse
1636 // parameter "param" for a particle "particle"
1637
1638 particle.ToLower();
1639 param. ToLower();
1640 Int_t p= -1;
1641 Int_t offset=0;
1642
1643 if (particle == "photon") offset=0;
1644 else if (particle == "pi0") offset=5;
1645 else
351dd634 1646 AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
1647 particle.Data()));
88cb7938 1648
1649 if (param.Contains("a")) p=4+offset;
1650 else if(param.Contains("b")) p=5+offset;
1651 else if(param.Contains("c")) p=6+offset;
1652 else if(param.Contains("x0"))p=7+offset;
1653 else if(param.Contains("y0"))p=8+offset;
351dd634 1654 if((i>4)||(i<0)) {
1655 AliError(Form("No parameter with index %d", i)) ;
1656 } else if(p==-1) {
1657 AliError(Form("No parameter with name %s", param.Data() )) ;
1658 } else
88cb7938 1659 (*fParameters)(p,i) = par ;
1660}
1661
407d15b3 1662//____________________________________________________________________________
1663void AliPHOSPIDv1::GetVertex(void)
1664{ //extract vertex either using ESD or generator
1665
1666 //Try to extract vertex from data
1667 if(fESD){
1668 const AliESDVertex *esdVtx = fESD->GetVertex() ;
e5b7b511 1669 if(esdVtx && esdVtx->GetChi2()!=0.){
407d15b3 1670 fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
1671 return ;
1672 }
1673 }
9a2cdbdf 1674
1675 // Use vertex diamond from CDB GRP folder if the one from ESD is missing
1676 // PLEASE FIX IT
407d15b3 1677 AliWarning("Can not read vertex from data, use fixed \n") ;
1678 fVtx.SetXYZ(0.,0.,0.) ;
1679
1680}
35adb638 1681//_______________________________________________________________________
1682void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
1683 // Sets values for the initial population of each particle type
304864ab 1684 for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i];
35adb638 1685}
1686//_______________________________________________________________________
1687void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
1688 // Gets values for the initial population of each particle type
304864ab 1689 for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i];
35adb638 1690}