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