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