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