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