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