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