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