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