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