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