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