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