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