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