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