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