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