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