take into account changing LHC periods between runs
[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.
602
603 Float_t p[]={0.,0.,0.};
604 for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
605 Float_t enerec = p[0] + p[1]*e + p[2]*e*e;
606 return enerec ;
607
608}
609
610//____________________________________________________________________________
fc7e2f43 611Float_t AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const
e3817e5f 612{
613 // Get the i-th parameter "CPV-EMC distance" for the specified axis
614 Float_t param = 0.;
351dd634 615 if(i>2 || i<0) {
616 AliError(Form("Invalid parameter number: %d",i));
617 } else {
e3817e5f 618 axis.ToLower();
351dd634 619 if (axis == "x")
620 param = (*fParameters)(1,i);
621 else if (axis == "z")
622 param = (*fParameters)(2,i);
623 else {
624 AliError(Form("Invalid axis name: %s",axis.Data()));
625 }
e3817e5f 626 }
627 return param;
628}
629
630//____________________________________________________________________________
fc7e2f43 631Float_t AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const
e3817e5f 632{
88cb7938 633 // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and
634 // Purity-Efficiency point
635
636 axis.ToLower();
637 Float_t p[]={0.,0.,0.};
638 for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
639 Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
640 return sig;
e3817e5f 641}
642
88cb7938 643//____________________________________________________________________________
fc7e2f43 644Float_t AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const
88cb7938 645{
646 // Calculates the parameter param of the ellipse
e3817e5f 647
648 particle.ToLower();
649 param. ToLower();
88cb7938 650 Float_t p[4]={0.,0.,0.,0.};
651 Float_t value = 0.0;
652 for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
653 if (particle == "photon") {
654 if (param.Contains("a")) e = TMath::Min((Double_t)e,70.);
655 else if (param.Contains("b")) e = TMath::Min((Double_t)e,70.);
656 else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
657 }
e3817e5f 658
443caba9 659 if (particle == "photon")
660 value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
661 else if (particle == "pi0")
662 value = p[0] + p[1]*e + p[2]*e*e;
663
88cb7938 664 return value;
e3817e5f 665}
666
667//_____________________________________________________________________________
fc7e2f43 668Float_t AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const
e3817e5f 669{
670 // Get the parameter "i" to calculate the boundary on the moment M2x
671 // for photons at high p_T
672 Float_t param = 0;
351dd634 673 if (i>3 || i<0) {
674 AliError(Form("Wrong parameter number: %d\n",i));
675 } else
e3817e5f 676 param = (*fParameters)(14,i) ;
677 return param;
148b2bba 678}
e3817e5f 679
148b2bba 680//____________________________________________________________________________
fc7e2f43 681Float_t AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const
e3817e5f 682{
683 // Get the parameter "i" to calculate the boundary on the moment M2x
684 // for pi0 at high p_T
685 Float_t param = 0;
351dd634 686 if (i>2 || i<0) {
687 AliError(Form("Wrong parameter number: %d\n",i));
688 } else
e3817e5f 689 param = (*fParameters)(15,i) ;
690 return param;
691}
148b2bba 692
e3817e5f 693//____________________________________________________________________________
fc7e2f43 694Float_t AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const
e3817e5f 695{
88cb7938 696 // Get TimeGate parameter depending on Purity-Efficiency i:
697 // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity
698 Float_t param = 0.;
351dd634 699 if(i>2 || i<0) {
700 AliError(Form("Invalid Efficiency-Purity choice %d",i));
701 } else
88cb7938 702 param = (*fParameters)(3,i) ;
703 return param;
e3817e5f 704}
705
e3817e5f 706//_____________________________________________________________________________
fc7e2f43 707Float_t AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const
88cb7938 708{
709 // Get the parameter "i" that is needed to calculate the ellipse
710 // parameter "param" for the particle "particle" ("photon" or "pi0")
711
e3817e5f 712 particle.ToLower();
713 param. ToLower();
88cb7938 714 Int_t offset = -1;
351dd634 715 if (particle == "photon")
716 offset=0;
717 else if (particle == "pi0")
718 offset=5;
e3817e5f 719 else
351dd634 720 AliError(Form("Wrong particle name: %s (choose from pi0/photon)\n",
721 particle.Data()));
88cb7938 722
723 Int_t p= -1;
724 Float_t par = 0;
e3817e5f 725
726 if (param.Contains("a")) p=4+offset;
727 else if(param.Contains("b")) p=5+offset;
728 else if(param.Contains("c")) p=6+offset;
729 else if(param.Contains("x0"))p=7+offset;
730 else if(param.Contains("y0"))p=8+offset;
12022e83 731
351dd634 732 if (i>4 || i<0) {
733 AliError(Form("No parameter with index %d", i)) ;
734 } else if (p==-1) {
735 AliError(Form("No parameter with name %s", param.Data() )) ;
736 } else
88cb7938 737 par = (*fParameters)(p,i) ;
738
739 return par;
12022e83 740}
741
12022e83 742
26aa7e4a 743//DP____________________________________________________________________________
744//Float_t AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSCpvRecPoint * cpv, Option_t * axis)const
745//{
746// // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
747//
7fb9892d 748// AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance();
26aa7e4a 749// TVector3 vecEmc ;
750// TVector3 vecCpv ;
751// if(cpv){
752// emc->GetLocalPosition(vecEmc) ;
753// cpv->GetLocalPosition(vecCpv) ;
754//
755// if(emc->GetPHOSMod() == cpv->GetPHOSMod()){
756// // Correct to difference in CPV and EMC position due to different distance to center.
757// // we assume, that particle moves from center
758// Float_t dCPV = geom->GetIPtoOuterCoverDistance();
759// Float_t dEMC = geom->GetIPtoCrystalSurface() ;
760// dEMC = dEMC / dCPV ;
761// vecCpv = dEMC * vecCpv - vecEmc ;
762// if (axis == "X") return vecCpv.X();
763// if (axis == "Y") return vecCpv.Y();
764// if (axis == "Z") return vecCpv.Z();
765// if (axis == "R") return vecCpv.Mag();
766// }
767// return 100000000 ;
768// }
769// return 100000000 ;
770//}
12022e83 771//____________________________________________________________________________
26aa7e4a 772Int_t AliPHOSPIDv1::GetCPVBit(AliPHOSTrackSegment * ts, Int_t effPur, Float_t e) const
bc0c084c 773{
c947e71a 774 //Calculates the pid bit for the CPV selection per each purity.
e3817e5f 775 if(effPur>2 || effPur<0)
351dd634 776 AliError(Form("Invalid Efficiency-Purity choice %d",effPur));
26aa7e4a 777
407d15b3 778//DP if(ts->GetCpvIndex()<0)
779//DP return 1 ; //no CPV cluster
bc0c084c 780
e3817e5f 781 Float_t sigX = GetCpv2EmcDistanceCut("X",e);
782 Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
bc0c084c 783
26aa7e4a 784 Float_t deltaX = TMath::Abs(ts->GetCpvDistance("X"));
785 Float_t deltaZ = TMath::Abs(ts->GetCpvDistance("Z"));
407d15b3 786// Info("GetCPVBit"," xdist %f, sigx %f, zdist %f, sigz %f",deltaX, sigX, deltaZ,sigZ) ;
7fb46731 787
788 //if(deltaX>sigX*(effPur+1))
789 //if((deltaX>sigX*(effPur+1)) || (deltaZ>sigZ*(effPur+1)))
790 if((deltaX>sigX*(effPur+1)) && (deltaZ>sigZ*(effPur+1)))
bc0c084c 791 return 1;//Neutral
792 else
793 return 0;//Charged
bc0c084c 794}
69183710 795
69183710 796//____________________________________________________________________________
fc7e2f43 797Int_t AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p, Int_t effPur, Float_t e)const
148b2bba 798{
50739f15 799 //Is the particle inside de PCA ellipse?
581354c5 800
e3817e5f 801 particle.ToLower();
802 Int_t prinbit = 0 ;
7fb46731 803 Float_t a = GetEllipseParameter(particle,"a" , e);
804 Float_t b = GetEllipseParameter(particle,"b" , e);
805 Float_t c = GetEllipseParameter(particle,"c" , e);
e3817e5f 806 Float_t x0 = GetEllipseParameter(particle,"x0", e);
807 Float_t y0 = GetEllipseParameter(particle,"y0", e);
808
809 Float_t r = TMath::Power((p[0] - x0)/a,2) +
810 TMath::Power((p[1] - y0)/b,2) +
811 c*(p[0] - x0)*(p[1] - y0)/(a*b) ;
50739f15 812 //3 different ellipses defined
e3817e5f 813 if((effPur==2) && (r<1./2.)) prinbit= 1;
814 if((effPur==1) && (r<2. )) prinbit= 1;
815 if((effPur==0) && (r<9./2.)) prinbit= 1;
50739f15 816
581354c5 817 if(r<0)
351dd634 818 AliError("Negative square?") ;
1f0e7ccd 819
820 return prinbit;
148b2bba 821
148b2bba 822}
1f0e7ccd 823//____________________________________________________________________________
fc7e2f43 824Int_t AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
1f0e7ccd 825{
e3817e5f 826 // Set bit for identified hard photons (E > 30 GeV)
827 // if the second moment M2x is below the boundary
828
829 Float_t e = emc->GetEnergy();
830 if (e < 30.0) return 0;
831 Float_t m2x = emc->GetM2x();
832 Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
833 TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
834 TMath::Power(GetParameterPhotonBoundary(2),2)) +
835 GetParameterPhotonBoundary(3);
351dd634 836 AliDebug(1, Form("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",
837 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 ;
7acf6008 915}
7b7c1533 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// }
7b7c1533 990//____________________________________________________________________________
2cc71c1e 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) {
1037 AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
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
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) {
351dd634 1366 AliFatal(Form("-> emc(%d) = %d", ts->GetEmcIndex(), emc )) ;
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
9a2cdbdf 1458 AliPHOSTrackSegment * ts = static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(rp->GetPHOSTSIndex()));
1459 AliPHOSEmcRecPoint * erp = static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(ts->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
8d0f3f77 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++;
351dd634 1569 AliDebug(1, Form("SetParameters", "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
1662//____________________________________________________________________________
407d15b3 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}