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