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