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