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