New PID version that uses TPrincipal
[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 //_________________________________________________________________________
19 // Implementation version v1 of the PHOS particle identifier 
20 // Particle identification based on the 
21 //     - RCPV: distance from CPV recpoint to EMCA recpoint.
22 //     - TOF 
23 //     - PCA: Principal Components Analysis..
24 // The identified particle has an identification number corresponding 
25 // to a 9 bits number:
26 //     -Bit 0 to 2: bit set if RCPV > fCpvEmcDistance (each bit corresponds
27 //      to a different efficiency-purity point of the photon identification) 
28 //     -Bit 3 to 5: bit set if TOF  < fTimeGate (each bit corresponds
29 //      to a different efficiency-purity point of the photon identification) 
30 //     -Bit 6 to 9: bit set if Principal Components are 
31 //      inside an ellipse defined by fX_center, fY_center, fA, fB, fAngle
32 //      (each bit corresponds to a different efficiency-purity point of the 
33 //      photon identification) 
34 //
35 //
36 // use case:
37 //  root [0] AliPHOSPIDv1 * p1 = new AliPHOSPIDv1("galice.root")
38 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
39 //  root [1] p1->SetIdentificationMethod("disp ellipse")
40 //  root [2] p1->ExecuteTask()
41 //  root [3] AliPHOSPIDv1 * p2 = new AliPHOSPIDv1("galice1.root","v1")
42 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
43 //                // reading headers from file galice1.root and create  RecParticles with title v1
44                   // TrackSegments and RecPoints with title "v1" are used 
45 //                // set file name for the branch RecParticles
46 //  root [4] p2->ExecuteTask("deb all time")
47 //                // available options
48 //                // "deb" - prints # of reconstructed particles
49 //                // "deb all" -  prints # and list of RecParticles
50 //                // "time" - prints benchmarking results
51 //                  
52 //  root [5] AliPHOSPIDv1 * p3 = new AliPHOSPIDv1("galice1.root","v1","v0")
53 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
54 //                // reading headers from file galice1.root and create  RecParticles with title v1
55                   // RecPoints and TrackSegments with title "v0" are used 
56 //  root [6] p3->ExecuteTask()
57 //*-- Author: Yves Schutz (SUBATECH)  & Gines Martinez (SUBATECH) & 
58 //            Gustavo Conesa April 2002
59
60 // --- ROOT system ---
61 #include "TROOT.h"
62 #include "TTree.h"
63 #include "TFile.h"
64 #include "TF2.h"
65 #include "TFormula.h"
66 #include "TCanvas.h"
67 #include "TFolder.h"
68 #include "TSystem.h"
69 #include "TBenchmark.h"
70 #include "TMatrixD.h"
71 #include "TPrincipal.h"
72 #include "TSystem.h"
73
74 // --- Standard library ---
75
76 #include <iostream.h>
77 #include <fstream.h>
78 #include <iomanip.h>
79
80 // --- AliRoot header files ---
81
82 #include "AliRun.h"
83 #include "AliGenerator.h"
84 #include "AliPHOS.h"
85 #include "AliPHOSPIDv1.h"
86 #include "AliPHOSClusterizerv1.h"
87 #include "AliPHOSTrackSegment.h"
88 #include "AliPHOSTrackSegmentMakerv1.h"
89 #include "AliPHOSRecParticle.h"
90 #include "AliPHOSGeometry.h"
91 #include "AliPHOSGetter.h"
92
93 ClassImp( AliPHOSPIDv1) 
94
95 //____________________________________________________________________________
96 AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID()
97
98   // default ctor
99  
100   fHeaderFileName    = "" ; 
101   fTrackSegmentsTitle= "" ; 
102   fRecPointsTitle    = "" ; 
103   fRecParticlesTitle = "" ; 
104   fRecParticlesInRun = 0 ;
105   fClusterizer       = 0 ; 
106   fTSMaker           = 0 ;
107   fFrom              = "" ;  
108   
109 }
110
111 //____________________________________________________________________________
112 AliPHOSPIDv1::AliPHOSPIDv1(const char * headerFile,const char * name, const char * from) : AliPHOSPID(headerFile, name)
113
114   //ctor with the indication on where to look for the track segments
115  
116   fHeaderFileName     = GetTitle() ; 
117   fTrackSegmentsTitle = GetName() ; 
118   fRecPointsTitle     = GetName() ; 
119   fRecParticlesTitle  = GetName() ;
120   TString tempo(GetName()) ; 
121   tempo.Append(":") ;
122   tempo.Append(Version()) ; 
123   SetName(tempo) ; 
124   fRecParticlesInRun = 0 ; 
125   if ( from == 0 ) 
126     fFrom = name ; 
127   else
128     fFrom = from ; 
129   Init() ;
130
131 }
132
133 //____________________________________________________________________________
134 AliPHOSPIDv1::~AliPHOSPIDv1()
135
136   delete [] fX ; // Principal input
137   delete [] fP ; // Principal components
138   delete fParameters ; // Matrix of Parameters 
139 }
140
141 //____________________________________________________________________________
142 void AliPHOSPIDv1::Init()
143 {
144   // Make all memory allocations that are not possible in default constructor
145   // Add the PID task to the list of PHOS tasks
146   
147   if ( strcmp(GetTitle(), "") == 0 )
148     SetTitle("galice.root") ;
149   
150   TString taskName(GetName()) ; 
151   taskName.Remove(taskName.Index(Version())-1) ;
152
153   // PCA : To do the Principal Components Analysis it is necessary 
154   // the Principal file, which is opened here
155   fX         = new double[7]; // Data for the PCA 
156   fP         = new double[7]; // Eigenvalues of the PCA  
157   fFileName  = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; 
158   TFile f( fFileName.Data(), "read" ) ;
159   fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ; 
160   f.Close() ; 
161   
162   // Initialization of the Parameters matrix. In the File Parameters.dat
163   // are all the parameters. These are introduced in a matrix of 21x3 elements.
164   // All the parameters defined in this file are, in order of row (there are
165   // 3 rows per parameter): CpvtoEmcDistanceCut, TimeGate (and the ellipse 
166   // parameters), X_center, Y_center, a, b, angle. Each row of a given parameter
167   // depends on the cluster energy range (0.3-1,1-2, >2 )
168   // Each column designes the parameters for a point in the Efficiency-Purity
169   // of the photon identification P1(0.959,0.625), P2(0.919,0.835), P3(0.833,0.901).
170  
171   fParameters = new TMatrixD(21,3) ; 
172
173   fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
174   ifstream paramFile(fFileNamePar, ios::in) ; 
175  
176   Int_t i,j ;
177
178   for(i = 0; i< 21; i++){
179     for(j = 0; j< 3; j++){
180       paramFile >> (*fParameters)(i,j) ;
181     }
182   }
183   paramFile.close();
184
185   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), fFrom.Data()) ; 
186
187   gime->SetRecParticlesTitle(taskName) ;
188   if ( gime == 0 ) {
189     cerr << "ERROR: AliPHOSPIDv1::Init -> Could not obtain the Getter object !" << endl ; 
190     return ;
191   } 
192    
193   gime->PostPID(this) ;
194   // create a folder on the white board //YSAlice/WhiteBoard/RecParticles/PHOS/recparticlesName
195   gime->PostRecParticles(taskName.Data() ) ; 
196   
197 }
198 //____________________________________________________________________________
199 Double_t  AliPHOSPIDv1::GetCpvtoEmcDistanceCut(const Float_t Cluster_En, const TString Eff_Pur)const 
200 {
201   // Get CpvtoEmcDistanceCut parameter depending on the cluster energy and 
202   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
203   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
204   // EFFICIENCY by PURITY)
205  
206   Int_t cluster = -1 ;
207   Int_t eff_pur = -1 ;
208
209   // Check the cluster energy range  
210   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
211   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
212   if( Cluster_En > 2.0) cluster = 2 ;
213   
214   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
215   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
216   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
217   
218   if(cluster ==-1){
219     cout<<"Invalid Cluster Energy option"<<endl;
220   }
221   else if(eff_pur ==-1){
222     cout<<"Invalid Efficiency-Purity option"<<endl;
223   }
224   else{
225     return (*fParameters)(cluster,eff_pur) ;
226   }
227 }
228 //____________________________________________________________________________
229
230 Double_t  AliPHOSPIDv1::GetTimeGate(const Float_t Cluster_En, const TString Eff_Pur) const 
231 {
232   // Get TimeGate parameter depending on the cluster energy and 
233   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
234   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
235   // EFFICIENCY by PURITY)
236   Int_t cluster = -1 ;
237   Int_t eff_pur = -1 ;
238   
239   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
240   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
241   if( Cluster_En > 2.0) cluster = 2 ;
242   
243   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
244   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
245   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
246    if(cluster ==-1){
247     cout<<"Invalid Cluster Energy option"<<endl;
248   }
249   else if(eff_pur ==-1){
250     cout<<"Invalid Efficiency-Purity option"<<endl; 
251   }
252   else{
253     return (*fParameters)(cluster+3,eff_pur) ; 
254   }
255 }
256 //_____________________________________________________________________________
257 Float_t  AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  Axis)const
258 {
259   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
260   
261   const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ; 
262   TVector3 vecEmc ;
263   TVector3 vecCpv ;
264   if(cpv){
265     emc->GetLocalPosition(vecEmc) ;
266     cpv->GetLocalPosition(vecCpv) ; 
267     if(emc->GetPHOSMod() == cpv->GetPHOSMod()){      
268       // Correct to difference in CPV and EMC position due to different distance to center.
269       // we assume, that particle moves from center
270       Float_t dCPV = geom->GetIPtoOuterCoverDistance();
271       Float_t dEMC = geom->GetIPtoCrystalSurface() ;
272       dEMC         = dEMC / dCPV ;
273       vecCpv = dEMC * vecCpv  - vecEmc ; 
274       if (Axis == "X") return vecCpv.X();
275       if (Axis == "Y") return vecCpv.Y();
276       if (Axis == "Z") return vecCpv.Z();
277       if (Axis == "R") return vecCpv.Mag();
278   } 
279     
280     return 100000000 ;
281   }
282   return 100000000 ;
283 }
284
285 //____________________________________________________________________________
286 Int_t  AliPHOSPIDv1::GetPrincipalSign(Double_t* P, Int_t cluster, Int_t eff_pur )const
287 {
288   //This method gives if the PCA of the particle are inside a defined ellipse
289   
290   // Get the parameters that define the ellipse stored in the 
291   // fParameters matrix.
292   Double_t X_center = (*fParameters)(cluster+6,eff_pur) ; 
293   Double_t Y_center = (*fParameters)(cluster+9,eff_pur) ; 
294   Double_t A        = (*fParameters)(cluster+12,eff_pur) ; 
295   Double_t B        = (*fParameters)(cluster+15,eff_pur) ; 
296   Double_t Angle    = (*fParameters)(cluster+18,eff_pur) ;
297
298   Int_t      prinsign;
299   Double_t   Dx        = 0. ; 
300   Double_t   Delta     = 0. ; 
301   Double_t   Y         = 0. ; 
302   Double_t   Y_1       = 0. ; 
303   Double_t   Y_2       = 0. ;
304   Double_t   Pi        = TMath::Pi() ;
305   Double_t   Cos_Theta = TMath::Cos(Pi*Angle/180.) ;
306   Double_t   Sin_Theta = TMath::Sin(Pi*Angle/180.) ;   
307
308   Dx = P[0] - X_center ; 
309   Delta = 4.*A*A*A*B* (A*A*Cos_Theta*Cos_Theta 
310                         + B*B*Sin_Theta*Sin_Theta - Dx*Dx) ; 
311   if (Delta < 0.) 
312     {prinsign=0;} 
313   
314   else if (Delta == 0.) 
315     { 
316       Y = Cos_Theta*Sin_Theta*(A*A - B*B)*Dx / 
317         (A*A*Cos_Theta*Cos_Theta + B*B*Sin_Theta*Sin_Theta) ; 
318       Y += Y_center ; 
319       if(P[1]==Y ) 
320         {prinsign=1;} 
321       else 
322         {prinsign=0;} 
323     } 
324   else 
325     { 
326       Y_1 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx +
327              TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta + 
328                                      B*B*Sin_Theta*Sin_Theta) ; 
329       Y_2 = (Cos_Theta*Sin_Theta*(A*A - B*B) *Dx -
330              TMath::Sqrt(Delta)/2.)/(A*A*Cos_Theta*Cos_Theta 
331                                       + B*B*Sin_Theta*Sin_Theta) ; 
332       Y_1 += Y_center ; 
333       Y_2 += Y_center ; 
334       if ((P[1]<=Y_1) && (P[1]>=Y_2)) 
335         {prinsign=1;} 
336       else 
337         {prinsign=0;}  
338     } 
339   return prinsign;
340 }
341
342 //____________________________________________________________________________
343  void  AliPHOSPIDv1::SetEllipseParameters(Float_t Cluster_En, TString Eff_Pur, Float_t x, Float_t y,Float_t a, Float_t b,Float_t angle)
344 {
345
346   // Set all ellipse parameters depending on the cluster energy and 
347   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
348   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
349   // EFFICIENCY by PURITY)
350
351   Int_t cluster = -1 ;
352   Int_t eff_pur = -1 ;
353
354  if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
355  if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
356  if( Cluster_En > 2.0) cluster= 2 ;
357
358  if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
359  if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
360  if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
361
362  if(cluster ==-1){
363     cout<<"Invalid Cluster Energy option"<<endl;
364   }
365   else if(eff_pur ==-1){
366     cout<<"Invalid Efficiency-Purity option"<<endl;
367   }
368   else{    
369     (*fParameters)(cluster+6,eff_pur) = x ;
370     (*fParameters)(cluster+9,eff_pur) = y ;
371     (*fParameters)(cluster+12,eff_pur) = a ;
372     (*fParameters)(cluster+15,eff_pur) = b ;
373     (*fParameters)(cluster+18,eff_pur) = angle ;
374   }
375
376 }
377 //__________________________________________________________________________ 
378 void  AliPHOSPIDv1::SetEllipseXCenter(Float_t Cluster_En, TString Eff_Pur, Float_t x) 
379 {
380   // Set the ellipse parameter x_center depending on the custer energy and 
381   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
382   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
383   // EFFICIENCY by PURITY)
384
385   Int_t cluster = -1 ;
386   Int_t eff_pur = -1 ;
387   
388   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
389   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
390   if( Cluster_En > 2.0) cluster = 2 ;
391   
392   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
393   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
394   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
395    
396   if(cluster ==-1){
397     cout<<"Invalid Cluster Energy option"<<endl;
398   }
399   else if(eff_pur ==-1){
400     cout<<"Invalid Efficiency-Purity option"<<endl;
401   }
402   else{
403     (*fParameters)(cluster+6,eff_pur) = x ; 
404   }
405 }
406 //_________________________________________________________________________    
407 void  AliPHOSPIDv1::SetEllipseYCenter(Float_t Cluster_En, TString Eff_Pur, Float_t y) 
408 {
409
410   // Set the ellipse parameter y_center depending on the cluster energy and 
411   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
412   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
413   // EFFICIENCY by PURITY)
414   
415   Int_t cluster = -1 ;
416   Int_t eff_pur = -1 ;
417   
418   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
419   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
420   if( Cluster_En > 2.0) cluster = 2 ;
421   
422   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
423   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
424   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
425
426   if(cluster ==-1){
427     cout<<"Invalid Cluster Energy option"<<endl;
428   }
429   else if(eff_pur ==-1){
430     cout<<"Invalid Efficiency-Purity option"<<endl ; 
431   }
432   else{
433     (*fParameters)(cluster+9,eff_pur) = y ;
434   }
435 }
436 //_________________________________________________________________________
437 void  AliPHOSPIDv1::SetEllipseAParameter(Float_t Cluster_En, TString Eff_Pur, Float_t a) 
438 {
439   // Set the ellipse parameter a depending on the cluster energy and 
440   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
441   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
442   // EFFICIENCY by PURITY)
443
444   Int_t cluster = -1 ;
445   Int_t eff_pur = -1 ;
446
447   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
448   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
449   if( Cluster_En > 2.0) cluster = 2 ;
450   
451   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
452   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
453   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
454   
455   if(cluster ==-1){
456     cout<<"Invalid Cluster Energy option"<<endl;
457   }
458   else if(eff_pur ==-1){
459     cout<<"Invalid Efficiency-Purity option"<<endl; 
460   }
461   else{ 
462     (*fParameters)(cluster+12,eff_pur) = a ;    
463   }
464 }
465 //________________________________________________________________________
466 void  AliPHOSPIDv1::SetEllipseBParameter(Float_t Cluster_En, TString Eff_Pur, Float_t b) 
467 {
468   // Set the ellipse parameter b depending on the cluster energy and 
469   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
470   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
471   // EFFICIENCY by PURITY)
472
473   Int_t cluster = -1 ;
474   Int_t eff_pur = -1 ;
475
476   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
477   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
478   if( Cluster_En > 2.0) cluster = 2 ;
479   
480   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
481   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
482  if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
483
484  if(cluster ==-1){
485    cout<<"Invalid Cluster Energy option"<<endl;
486  }
487  else if(eff_pur ==-1){
488    cout<<"Invalid Efficiency-Purity option"<<endl; 
489  }
490  else{
491    (*fParameters)(cluster+15,eff_pur) = b ;
492  }
493 }
494 //________________________________________________________________________
495 void  AliPHOSPIDv1::SetEllipseAngle(Float_t Cluster_En, TString Eff_Pur, Float_t angle) 
496 {
497
498   // Set the ellipse parameter angle depending on the cluster energy and 
499   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
500   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
501   // EFFICIENCY by PURITY)
502
503   Int_t cluster = -1 ;
504   Int_t eff_pur = -1 ;
505
506   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
507   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
508   if( Cluster_En > 2.0) cluster = 2 ;
509   
510   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
511   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
512   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
513  
514   if(cluster ==-1){
515     cout<<"Invalid Cluster Energy option"<<endl;
516   }
517   else if(eff_pur ==-1){
518     cout<<"Invalid Efficiency-Purity option"<<endl; 
519   }
520   else{
521     (*fParameters)(cluster+18,eff_pur) = angle ;
522   }
523
524 //_____________________________________________________________________________
525 void  AliPHOSPIDv1::SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut) 
526 {
527
528   // Set the parameter Cpvto EmcDistanceCut depending on the cluster energy and 
529   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
530   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
531   // EFFICIENCY by PURITY)
532
533   Int_t cluster = -1 ;
534   Int_t eff_pur = -1 ;
535
536   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
537   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
538   if( Cluster_En > 2.0) cluster = 2 ;
539   
540   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
541   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
542   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
543   
544   if(cluster ==-1){
545     cout<<"Invalid Cluster Energy option"<<endl;
546   }
547   else if(eff_pur ==-1){
548     cout<<"Invalid Efficiency-Purity option"<<endl; 
549   }
550   else{
551     (*fParameters)(cluster,eff_pur) = cut ;
552   }
553 }  
554 //_____________________________________________________________________________
555 void  AliPHOSPIDv1::SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate) 
556 {
557
558   // Set the parameter TimeGate depending on the cluster energy and 
559   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
560   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
561   // EFFICIENCY by PURITY)
562
563   Int_t cluster = -1 ;
564   Int_t eff_pur = -1 ;
565
566   if((Cluster_En > 0.3)&&(Cluster_En <= 1.0)) cluster = 0 ;
567   if((Cluster_En > 1.0)&&(Cluster_En <= 2.0)) cluster = 1 ;
568   if( Cluster_En > 2.0) cluster = 2 ;
569   
570   if(strstr( Eff_Pur,"HIGH EFFICIENCY" )||strstr( Eff_Pur,"LOW PURITY" ) ) eff_pur = 0 ;
571   if(strstr( Eff_Pur,"MEDIUM EFFICIENCY" )||strstr( Eff_Pur,"MEDIUM PURITY" ) ) eff_pur = 1 ;
572   if(strstr( Eff_Pur,"LOW EFFICIENCY" )||strstr( Eff_Pur,"HIGH PURITY" ) ) eff_pur = 2 ;
573
574   if(cluster ==-1){
575     cout<<"Invalid Cluster Energy option"<<endl;
576   }
577   else if(eff_pur ==-1){
578     cout<<"Invalid Efficiency-Purity option"<<endl; 
579   }
580   else{
581     (*fParameters)(cluster+3,eff_pur) = gate ;
582   }
583
584 //____________________________________________________________________________
585
586 void  AliPHOSPIDv1::Exec(Option_t * option) 
587 {
588   //Steering method
589   
590   if( strcmp(GetName(), "")== 0 ) 
591     Init() ;
592   
593   if(strstr(option,"tim"))
594     gBenchmark->Start("PHOSPID");
595   
596   if(strstr(option,"print")) {
597     Print("") ; 
598     return ; 
599   }
600
601   //cout << gDirectory->GetName() << endl ; 
602
603   gAlice->GetEvent(0) ;
604
605   //check, if the branch with name of this" already exits?
606   TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
607   TIter next(lob) ; 
608   TBranch * branch = 0 ;  
609   Bool_t phospidfound = kFALSE, pidfound = kFALSE ; 
610   
611   TString taskName(GetName()) ; 
612   taskName.Remove(taskName.Index(Version())-1) ;
613
614   while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
615     if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) 
616       phospidfound = kTRUE ;
617     
618     else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) 
619       pidfound = kTRUE ; 
620   }
621
622   if ( phospidfound || pidfound ) {
623     cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name " 
624          << taskName.Data() << " already exits" << endl ;
625     return ; 
626   }       
627   
628   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
629   Int_t ievent ;
630   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;  
631   for(ievent = 0; ievent < nevents; ievent++){
632     gime->Event(ievent,"RA") ;
633  
634     MakeRecParticles() ;
635     
636     WriteRecParticles(ievent);
637     
638     if(strstr(option,"deb"))
639       PrintRecParticles(option) ;
640
641     //increment the total number of rec particles per run 
642     fRecParticlesInRun += gime->RecParticles(taskName)->GetEntriesFast() ; 
643
644   }
645   
646   if(strstr(option,"tim")){
647     gBenchmark->Stop("PHOSPID");
648     cout << "AliPHOSPID:" << endl ;
649     cout << "  took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID " 
650          <<  gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ;
651     cout << endl ;
652   }
653   
654 }
655
656 //____________________________________________________________________________
657 void  AliPHOSPIDv1::MakeRecParticles(){
658
659   // Makes a RecParticle out of a TrackSegment
660
661   TString taskName(GetName()) ; 
662   taskName.Remove(taskName.Index(Version())-1) ;
663   
664   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
665   TObjArray * emcRecPoints = gime->EmcRecPoints(fFrom) ; 
666   TObjArray * cpvRecPoints = gime->CpvRecPoints(fFrom) ; 
667   TClonesArray * trackSegments = gime->TrackSegments(fFrom) ; 
668   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
669     cerr << "ERROR:  AliPHOSPIDv1::MakeRecParticles -> RecPoints or TrackSegments with name " 
670          << fFrom << " not found ! " << endl ; 
671     abort() ; 
672   }
673   TClonesArray * recParticles  = gime->RecParticles(taskName) ; 
674   recParticles->Clear();
675  
676
677   TIter next(trackSegments) ; 
678   AliPHOSTrackSegment * ts ; 
679   Int_t index = 0 ; 
680   AliPHOSRecParticle * rp ; 
681   
682   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
683     
684     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
685     rp = (AliPHOSRecParticle *)recParticles->At(index) ; 
686     rp->SetTraskSegment(index) ;
687     rp->SetIndexInList(index) ;
688         
689     AliPHOSEmcRecPoint * emc = 0 ;
690     if(ts->GetEmcIndex()>=0)
691       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
692     
693     AliPHOSRecPoint    * cpv = 0 ;
694     if(ts->GetCpvIndex()>=0)
695       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
696     
697     //set momentum and energy first
698     Float_t    e = emc->GetEnergy() ;
699     TVector3 dir = GetMomentumDirection(emc,cpv) ; 
700     dir.SetMag(e) ;
701
702     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),e) ;
703     rp->SetCalcMass(0);
704     
705     // Now set type (reconstructed) of the particle
706
707     // Choose the cluster energy range
708     Int_t cluster = 0 ; // Ellipse and rcpv cut in function of the cluster energy
709     if((e > 0.3)&&(e <= 1.0)) cluster = 0 ;
710     if((e > 1.0)&&(e <= 2.0)) cluster = 1 ;
711     if( e > 2.0) cluster = 2 ;
712
713     // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken 
714     // into account to set the particle identification)
715     for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){
716
717     // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, 1st, 
718     // 2nd or 3rd bit (depending on the efficiency-purity point )is set to 1 .  
719       if(GetDistance(emc, cpv,  "R") > (*fParameters)(cluster,eff_pur) )  
720         rp->SetPIDBit(eff_pur) ;
721     
722     // Looking the TOF. If TOF smaller than gate,  4th, 5th or 6th 
723     // bit (depending on the efficiency-purity point )is set to 1             
724       if(emc->GetTime()< (*fParameters)(cluster+3,eff_pur))  
725         rp->SetPIDBit(eff_pur+3) ;                  
726      
727     // Looking PCA. Define and calculate the data (X), introduce in the function 
728     // X2P that gives the components (P).  
729     Float_t    fSpher = 0. ;
730     Float_t    fEmaxdtotal = 0. ; 
731     Float_t lambda[2] ;
732     
733     emc->GetElipsAxis(lambda) ;
734     if((lambda[0]+lambda[1])!=0) fSpher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
735     
736     fEmaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); 
737     
738     fX[0] = lambda[0] ;  
739     fX[1] = lambda[1] ; 
740     fX[2] = emc->GetDispersion() ; 
741     fX[3] = fSpher ; 
742     fX[4] = emc->GetMultiplicity() ;  
743     fX[5] = fEmaxdtotal ;  
744     fX[6] = emc->GetCoreEnergy() ;  
745     
746     fPrincipal->X2P(fX,fP);
747     
748     //If we are inside the ellipse, 7th, 8th or 9th 
749     // bit (depending on the efficiency-purity point )is set to 1 
750     if(GetPrincipalSign(fP,cluster,eff_pur) == 1) 
751       rp->SetPIDBit(eff_pur+6) ;
752  
753   }
754     rp->SetProductionVertex(0,0,0,0);
755     rp->SetFirstMother(-1);
756     rp->SetLastMother(-1);
757     rp->SetFirstDaughter(-1);
758     rp->SetLastDaughter(-1);
759     rp->SetPolarisation(0,0,0);
760     index++ ; 
761   }
762   
763 }
764
765 //____________________________________________________________________________
766 void  AliPHOSPIDv1:: Print(Option_t * option) const
767 {
768   // Print the parameters used for the particle type identification
769     cout <<  "=============== AliPHOSPID1 ================" << endl ;
770     cout <<  "Making PID "<< endl ;
771     cout <<  "    Headers file:               " << fHeaderFileName.Data() << endl ;
772     cout <<  "    RecPoints branch title:     " << fRecPointsTitle.Data() << endl ;
773     cout <<  "    TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ;
774     cout <<  "    RecParticles Branch title   " << fRecParticlesTitle.Data() << endl;
775     cout <<  "with parameters: " << endl ;
776     //    cout <<  "    Maximal EMC - CPV  distance (cm) " << fCpvEmcDistance << endl ;
777     //cout <<  "    Time Gate used:                  " << fTimeGate <<  endl ;
778     //cout <<  "    Principal Ellipse Parameters     " << endl ;
779     //cout <<  "       Ellipse center   (x,y)           (" << fX_center<<","<<fY_center<<")"<< endl;
780     //cout <<  "       Ellipse focus    (a,b)           (" << fA<<","<<fB<<")"<< endl;
781     //cout <<  "       Ellipse angle                     " << fAngle<< endl;        
782     cout <<  "============================================" << endl ;
783 }
784
785 //____________________________________________________________________________
786 void  AliPHOSPIDv1::WriteRecParticles(Int_t event)
787 {
788  
789   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; 
790   TString taskName(GetName()) ; 
791   taskName.Remove(taskName.Index(Version())-1) ;
792   TClonesArray * recParticles = gime->RecParticles(taskName) ; 
793   recParticles->Expand(recParticles->GetEntriesFast() ) ;
794
795   //Make branch in TreeR for RecParticles 
796   char * filename = 0;
797   if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){   //generating file name
798     filename = new char[strlen(gAlice->GetBaseFile())+20] ;
799     sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ; 
800   }
801   
802   TDirectory *cwd = gDirectory;
803   
804   //First rp
805   Int_t bufferSize = 32000 ;    
806   TBranch * rpBranch = gAlice->TreeR()->Branch("PHOSRP",&recParticles,bufferSize);
807   rpBranch->SetTitle(fRecParticlesTitle);
808   if (filename) {
809     rpBranch->SetFile(filename);
810     TIter next( rpBranch->GetListOfBranches());
811     TBranch * sb ;
812     while ((sb=(TBranch*)next())) {
813       sb->SetFile(filename);
814     }   
815     cwd->cd();
816   }
817   
818   //second, pid
819   Int_t splitlevel = 0 ; 
820   AliPHOSPIDv1 * pid = this ;
821   TBranch * pidBranch = gAlice->TreeR()->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel);
822   pidBranch->SetTitle(fRecParticlesTitle.Data());
823   if (filename) {
824     pidBranch->SetFile(filename);
825     TIter next( pidBranch->GetListOfBranches());
826     TBranch * sb ;
827     while ((sb=(TBranch*)next())) {
828       sb->SetFile(filename);
829     }   
830     cwd->cd();
831   }    
832   
833   rpBranch->Fill() ;
834   pidBranch->Fill() ;
835   
836   gAlice->TreeR()->Write(0,kOverwrite) ;  
837   //pidBranch->Write(0,kOverwrite) ;  
838
839   delete [] filename ; 
840 }
841
842 //____________________________________________________________________________
843 TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const 
844
845   // Calculates the momentum direction:
846   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
847   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints 
848   //  However because of the poor position resolution of PPSD the direction is always taken as if we were 
849   //  in case 1.
850
851   TVector3 dir(0,0,0) ; 
852   
853   TVector3 emcglobalpos ;
854   TMatrix  dummy ;
855   
856   emc->GetGlobalPosition(emcglobalpos, dummy) ;
857   
858
859   dir = emcglobalpos ;  
860   dir.SetZ( -dir.Z() ) ;   // why ?  
861   dir.SetMag(1.) ;
862
863   //account correction to the position of IP
864   Float_t xo,yo,zo ; //Coordinates of the origin
865   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
866   TVector3 origin(xo,yo,zo);
867   dir = dir - origin ;
868
869   return dir ;  
870 }
871 //____________________________________________________________________________
872 void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
873 {
874   // Print table of reconstructed particles
875
876   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; 
877
878   TString taskName(GetName()) ; 
879   taskName.Remove(taskName.Index(Version())-1) ;
880   TClonesArray * recParticles = gime->RecParticles(taskName) ; 
881   
882   cout << "AliPHOSPIDv1: event "<<gAlice->GetEvNumber()  << endl ;
883   cout << "       found " << recParticles->GetEntriesFast() << " RecParticles " << endl ;
884   
885   if(strstr(option,"all")) {  // printing found TS
886     
887     cout << "  PARTICLE "   
888          << "  Index    "  << endl ;
889     
890     Int_t index ;
891     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
892        AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;       
893        //cout<<" Type " <<rp->GetType()<<endl;
894        //Text_t particle[11];
895        // switch(rp->GetType()) {
896 //        case  AliPHOSFastRecParticle::kCHARGEDHASLOW:
897 //       strcpy(particle, "CHARGED HA SLOW") ;
898 //       break ; 
899 //        case  AliPHOSFastRecParticle::kNEUTRALHASLOW: 
900 //       strcpy(particle, "NEUTRAL HA SLOW");
901 //       break ;    
902 //        case  AliPHOSFastRecParticle::kCHARGEDHAFAST:
903 //       strcpy(particle, "CHARGED HA FAST") ;
904 //       break ;        
905 //        case  AliPHOSFastRecParticle::kNEUTRALHAFAST:
906 //       strcpy(particle, "NEUTRAL HA FAST");
907 //       break;
908 //        case  AliPHOSFastRecParticle::kCHARGEDEMSLOW:
909 //       strcpy(particle, "CHARGED EM SLOW") ;
910 //       break ;
911 //        case  AliPHOSFastRecParticle::kNEUTRALEMSLOW:
912 //       strcpy(particle, "NEUTRAL EM SLOW");
913 //       break ;
914 //        case  AliPHOSFastRecParticle::kCHARGEDEMFAST:
915 //       strcpy(particle, "CHARGED EM FAST") ;
916 //       break ;
917 //        case  AliPHOSFastRecParticle::kNEUTRALEMFAST:
918 //       strcpy( particle, "NEUTRAL EM FAST");
919 //       break;
920 //       }
921
922       cout << setw(10) << rp->GetType() << "  "
923            << setw(5) <<  rp->GetIndexInList() << " " <<endl;
924    
925     }
926     cout << "-------------------------------------------" << endl ;
927   }
928   
929 }
930
931
932