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