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