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