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