Particle identification improvment by Gustavo
[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     cerr << "ERROR: AliPHOSPIDv1::Init -> Could not obtain the Getter object !" << endl ; 
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)cout<<"Error: Negative square?"<<endl;
353   return prinsign;
354
355 }
356
357 //_____________________________________________________________________________
358 void  AliPHOSPIDv1::SetCpvtoEmcDistanceCut(Float_t Cluster_En, TString Eff_Pur, Float_t cut) 
359 {
360
361   // Set the parameter Cpvto EmcDistanceCut depending on the cluster energy and 
362   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
363   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
364   // EFFICIENCY by PURITY)
365
366
367   Int_t eff_pur = GetEffPurOption(Eff_Pur);
368   Int_t cluster = GetClusterOption(Cluster_En) ;
369   if((cluster!= -1)&&(eff_pur != -1))
370     (*fParameters)(cluster,eff_pur) = cut ;
371 }
372 //_____________________________________________________________________________
373 void  AliPHOSPIDv1::SetTimeGate(Float_t Cluster_En, TString Eff_Pur, Float_t gate) 
374 {
375
376   // Set the parameter TimeGate depending on the cluster energy and 
377   // Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
378   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
379   // EFFICIENCY by PURITY)
380     
381   Int_t eff_pur = GetEffPurOption(Eff_Pur);
382   Int_t cluster = GetClusterOption(Cluster_En) ;
383   if((cluster!= -1)&&(eff_pur != -1))
384     (*fParameters)(cluster+6,eff_pur) = gate ;
385
386 //_____________________________________________________________________________
387 void  AliPHOSPIDv1::SetParameters() 
388 {
389   // PCA : To do the Principal Components Analysis it is necessary 
390   // the Principal file, which is opened here
391   fX         = new double[7]; // Data for the PCA 
392   fP         = new double[7]; // Eigenvalues of the PCA
393   
394   // Open principal and parameters files to be used
395   
396   fFileName  = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ;
397   fFileNamePar = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat"); 
398   TFile f( fFileName.Data(), "read" ) ;
399   fPrincipal = dynamic_cast<TPrincipal*> (f.Get("principal")) ; 
400   f.Close() ; 
401   
402   // Initialization of the Parameters matrix. In the File Parameters.dat
403   // are all the parameters. These are introduced in a matrix of 18x4  
404   // 
405   // All the parameters defined in this file are, in order of row: 
406   // CpvtoEmcDistanceCut (6 rows, each one depends on the energy range of the 
407   // particle, and 3 columns, each one depending on the efficiency-purity 
408   // point), TimeGate (the same) and the parameters of the functions that 
409   // calculate the ellipse parameters, x0,y0,a, b, c. These 5 parameters 
410   // (5 rows) depend on 4 parameters (columns). Finally there is a row with 
411   // the energy calibration parameters, 3 parameters. 
412   
413   fParameters = new TMatrixD(18,4) ;
414   
415   ifstream paramFile(fFileNamePar) ; 
416   Int_t h,n;
417   for(h = 0; h< 18; h++){
418     for(n = 0; n< 4; n++){
419       paramFile >> (*fParameters)(h,n) ;
420     }
421   }
422   paramFile.close();
423 }
424 //_____________________________________________________________________________
425 const Int_t  AliPHOSPIDv1::GetClusterOption(const Float_t Cluster_En) const
426 {
427   // Gives the cluster energy range, for each range there is associated a TOF or RCPV
428   // parameter.
429   Int_t cluster = -1;
430   if((Cluster_En > 0.0 )&&(Cluster_En <= 2.0 )) cluster = 0 ;
431   if((Cluster_En > 2.0 )&&(Cluster_En <= 5.0 )) cluster = 1 ;
432   if((Cluster_En > 5.0 )&&(Cluster_En <= 10.0)) cluster = 2 ;
433   if((Cluster_En > 10.0)&&(Cluster_En <= 20.0)) cluster = 3 ;
434   if((Cluster_En > 20.0)&&(Cluster_En <= 30.0)) cluster = 4 ;
435   if( Cluster_En > 30.0)                        cluster = 5 ;
436   
437   return cluster;
438 }
439 //____________________________________________________________________________
440 const Int_t  AliPHOSPIDv1::GetEffPurOption(const TString Eff_Pur) const
441 {
442
443   // Looks for the Purity-Efficiency point (possible options "HIGH EFFICIENCY" 
444   // "MEDIUM EFFICIENCY" "LOW EFFICIENCY" and 3 more options changing 
445   // EFFICIENCY by PURITY)
446
447   Int_t eff_pur = -1 ;
448
449   if(Eff_Pur.Contains("HIGH EFFICIENCY") ||Eff_Pur.Contains("LOW PURITY") )
450     eff_pur = 0 ;
451   else if(Eff_Pur.Contains("MEDIUM EFFICIENCY") ||Eff_Pur.Contains("MEDIUM PURITY") ) 
452     eff_pur = 1 ;
453   else if(Eff_Pur.Contains("LOW EFFICIENCY")||Eff_Pur.Contains("HIGH PURITY") ) 
454     eff_pur = 2 ;
455   else{
456     eff_pur = -1;
457     cout<<"Invalid Efficiency-Purity option"<<endl;
458     cout<<"Possible options: HIGH EFFICIENCY =    LOW PURITY"<<endl;
459     cout<<"                MEDIUM EFFICIENCY = MEDIUM PURITY"<<endl;
460     cout<<"                   LOW EFFICIENCY =   HIGH PURITY"<<endl;
461   }
462
463   return eff_pur;
464 }
465 //________________________________________________________________________
466 void  AliPHOSPIDv1::SetEllipseParameter(TString Param, Int_t i, Double_t par) 
467 {  
468   // Set the parameter "i" that is needed to calculate the ellipse 
469   // parameter "Param".
470
471   Int_t p= -1;
472   
473   if(Param.Contains("a"))p=12; 
474   if(Param.Contains("b"))p=13; 
475   if(Param.Contains("c"))p=14; 
476   if(Param.Contains("x0"))p=15; 
477   if(Param.Contains("y0"))p=16;
478   if((i>4)||(i<0))
479     cout<<"Error:: No parameter with index "<<i<<endl; 
480   else if(p==-1)
481     cout<<"Error:: No parameter with name "<<Param<<endl; 
482   else
483     (*fParameters)(p,i) = par ;
484
485 //________________________________________________________________________
486 const Double_t  AliPHOSPIDv1::GetParameterToCalculateEllipse(const TString Param, const Int_t i) const
487
488   // Get the parameter "i" that is needed to calculate the ellipse 
489   // parameter "Param".
490
491   Int_t p= -1;
492   Double_t par = -1;
493
494   if(Param.Contains("a"))p=12; 
495   if(Param.Contains("b"))p=13; 
496   if(Param.Contains("c"))p=14; 
497   if(Param.Contains("x0"))p=15; 
498   if(Param.Contains("y0"))p=16;
499
500   if((i>4)||(i<0))
501     cout<<"Error:: No parameter with index "<<i<<endl; 
502   else if(p==-1)
503     cout<<"Error:: No parameter with name "<<Param<<endl; 
504   else
505     par = (*fParameters)(p,i) ;
506   
507   return par;
508
509
510 //____________________________________________________________________________
511 void  AliPHOSPIDv1::SetCalibrationParameter(Int_t i,Double_t param) 
512 {
513   (*fParameters)(17,i) = param ;
514 }
515 //____________________________________________________________________________
516 const Double_t  AliPHOSPIDv1::GetCalibrationParameter(const Int_t i) const 
517 {
518   Float_t param = (*fParameters)(17,i);
519   return param;
520 }
521 //____________________________________________________________________________
522 const Double_t  AliPHOSPIDv1::GetEllipseParameter(const TString Param,Float_t E) const 
523 {
524   Double_t p[4]={0.,0.,0.,0.};
525   Double_t value = 0.0;
526   Int_t i;
527
528   if(Param.Contains("a")){
529     for(i=0;i<4;i++)p[i]=(*fParameters)(12,i);
530     if(E>70.)E=70.;
531   }
532   
533   if(Param.Contains("b")){
534     for(i=0;i<4;i++)p[i]=(*fParameters)(13,i);
535     if(E>70.)E=70.;
536   }
537   
538   if(Param.Contains("c"))
539     for(i=0;i<4;i++)p[i]=(*fParameters)(14,i);
540   
541   if(Param.Contains("x0")){
542     for(i=0;i<4;i++)p[i]=(*fParameters)(15,i);
543     if(E<1.)E=1.1;
544   }
545   if(Param.Contains("y0"))
546     for(i=0;i<4;i++)p[i]=(*fParameters)(16,i);
547   
548   value = p[0]/TMath::Sqrt(E)+p[1]*E+p[2]*E*E+p[3];
549   return value;
550 }
551 //____________________________________________________________________________
552
553 void  AliPHOSPIDv1::Exec(Option_t * option) 
554 {
555   //Steering method
556   
557   if( strcmp(GetName(), "")== 0 ) 
558     Init() ;
559   
560   if(strstr(option,"tim"))
561     gBenchmark->Start("PHOSPID");
562   
563   if(strstr(option,"print")) {
564     Print("") ; 
565     return ; 
566   }
567
568
569 //   gAlice->GetEvent(0) ;
570
571 //   //check, if the branch with name of this" already exits?
572 //   if (gAlice->TreeR()) {
573 //     TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
574 //     TIter next(lob) ; 
575 //     TBranch * branch = 0 ;  
576 //     Bool_t phospidfound = kFALSE, pidfound = kFALSE ; 
577     
578 //     TString taskName(GetName()) ; 
579 //     taskName.Remove(taskName.Index(Version())-1) ;
580     
581 //     while ( (branch = (TBranch*)next()) && (!phospidfound || !pidfound) ) {
582 //       if ( (strcmp(branch->GetName(), "PHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) 
583 //      phospidfound = kTRUE ;
584       
585 //       else if ( (strcmp(branch->GetName(), "AliPHOSPID")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) 
586 //      pidfound = kTRUE ; 
587 //     }
588     
589 //     if ( phospidfound || pidfound ) {
590 //       cerr << "WARNING: AliPHOSPIDv1::Exec -> RecParticles and/or PIDtMaker branch with name " 
591 //         << taskName.Data() << " already exits" << endl ;
592 //       return ; 
593 //     }       
594 //   }
595
596 //   Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
597 //   Int_t ievent ;
598 //   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;  
599
600   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
601   if(gime->BranchExists("RecParticles") )
602     return ;
603   Int_t nevents = gime->MaxEvent() ;       //(Int_t) gAlice->TreeE()->GetEntries() ;
604   Int_t ievent ;
605
606
607   for(ievent = 0; ievent < nevents; ievent++){
608     gime->Event(ievent,"R") ;
609  
610     MakeRecParticles() ;
611     
612     WriteRecParticles(ievent);
613     
614     if(strstr(option,"deb"))
615       PrintRecParticles(option) ;
616
617     //increment the total number of rec particles per run 
618     fRecParticlesInRun += gime->RecParticles(BranchName())->GetEntriesFast() ; 
619
620   }
621   
622   if(strstr(option,"tim")){
623     gBenchmark->Stop("PHOSPID");
624     cout << "AliPHOSPID:" << endl ;
625     cout << "  took " << gBenchmark->GetCpuTime("PHOSPID") << " seconds for PID " 
626          <<  gBenchmark->GetCpuTime("PHOSPID")/nevents << " seconds per event " << endl ;
627     cout << endl ;
628   }
629   
630 }
631
632 //____________________________________________________________________________
633 void  AliPHOSPIDv1::MakeRecParticles(){
634
635   // Makes a RecParticle out of a TrackSegment
636   
637   AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; 
638   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
639   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
640   TClonesArray * trackSegments = gime->TrackSegments() ; 
641   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
642     cerr << "ERROR:  AliPHOSPIDv1::MakeRecParticles -> RecPoints or TrackSegments not found ! " << endl ; 
643     abort() ; 
644   }
645   TClonesArray * recParticles  = gime->RecParticles() ; 
646   recParticles->Clear();
647  
648
649   TIter next(trackSegments) ; 
650   AliPHOSTrackSegment * ts ; 
651   Int_t index = 0 ; 
652   AliPHOSRecParticle * rp ; 
653   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
654     
655     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
656     rp = (AliPHOSRecParticle *)recParticles->At(index) ; 
657     rp->SetTrackSegment(index) ;
658     rp->SetIndexInList(index) ;
659         
660     AliPHOSEmcRecPoint * emc = 0 ;
661     if(ts->GetEmcIndex()>=0)
662       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
663     
664     AliPHOSRecPoint    * cpv = 0 ;
665     if(ts->GetCpvIndex()>=0)
666       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
667     
668     // Now set type (reconstructed) of the particle
669
670     // Choose the cluster energy range
671     
672     // YK: check if (emc != 0) !!!
673     if (!emc) {
674       cerr << "ERROR:  AliPHOSPIDv1::MakeRecParticles -> emc("
675            <<ts->GetEmcIndex()<<") = "         <<emc<< endl;
676       abort();
677     }
678     Float_t    e = emc->GetEnergy() ;   
679     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.
680     if(cluster== -1) continue ;
681
682     Float_t  lambda[2] ;
683     emc->GetElipsAxis(lambda) ;
684     
685     if((lambda[0]>0.01) && (lambda[1]>0.01)){
686       // Looking PCA. Define and calculate the data (X),
687       // introduce in the function 
688       // X2P that gives the components (P).  
689       Float_t  Spher = 0. ;
690       Float_t  Emaxdtotal = 0. ; 
691       
692       if((lambda[0]+lambda[1])!=0) Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
693       
694       Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); 
695       
696       fX[0] = lambda[0] ;  
697       fX[1] = lambda[1] ; 
698       fX[2] = emc->GetDispersion() ; 
699       fX[3] = Spher ; 
700       fX[4] = emc->GetMultiplicity() ;  
701       fX[5] = Emaxdtotal ;  
702       fX[6] = emc->GetCoreEnergy() ;  
703       
704       fPrincipal->X2P(fX,fP);
705     }
706     else{
707       fP[0]=-100.0;  //We do not accept clusters with 
708       fP[1]=-100.0;  //one cell as a photon-like
709     }
710     
711     Float_t time =emc->GetTime() ;
712     
713     // Loop of Efficiency-Purity (the 3 points of purity or efficiency are taken 
714     // into account to set the particle identification)
715     for(Int_t eff_pur = 0; eff_pur < 3 ; eff_pur++){
716       
717       // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, 1st, 
718       // 2nd or 3rd bit (depending on the efficiency-purity point )is set to 1 . 
719       
720       if(GetDistance(emc, cpv,  "R") > (*fParameters)(cluster,eff_pur) )  
721         rp->SetPIDBit(eff_pur) ;
722       
723       // Looking the TOF. If TOF smaller than gate,  4th, 5th or 6th 
724       // bit (depending on the efficiency-purity point )is set to 1             
725       if(time< (*fParameters)(cluster+6,eff_pur)) {
726         rp->SetPIDBit(eff_pur+3) ;                  
727       }
728       
729       //If we are inside the ellipse, 7th, 8th or 9th 
730       // bit (depending on the efficiency-purity point )is set to 1 
731       if(GetPrincipalSign(fP,eff_pur,e) == 1) 
732         rp->SetPIDBit(eff_pur+6) ;
733     }
734       
735     //Set momentum, energy and other parameters 
736     Float_t  encal = GetCalibratedEnergy(e);
737     TVector3 dir   = GetMomentumDirection(emc,cpv) ; 
738     dir.SetMag(encal) ;
739     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
740     rp->SetCalcMass(0);
741     rp->Name(); //If photon sets the particle pdg name to gamma
742     rp->SetProductionVertex(0,0,0,0);
743     rp->SetFirstMother(-1);
744     rp->SetLastMother(-1);
745     rp->SetFirstDaughter(-1);
746     rp->SetLastDaughter(-1);
747     rp->SetPolarisation(0,0,0);
748     index++ ; 
749   }
750   
751 }
752
753 //____________________________________________________________________________
754 void  AliPHOSPIDv1:: Print()
755 {
756   // Print the parameters used for the particle type identification
757     cout <<  "=============== AliPHOSPID1 ================" << endl ;
758     cout <<  "Making PID "<< endl ;
759 //     cout <<  "    Headers file:               " << fHeaderFileName.Data() << endl ;
760 //     cout <<  "    RecPoints branch title:     " << fRecPointsTitle.Data() << endl ;
761 //     cout <<  "    TrackSegments Branch title: " << fTrackSegmentsTitle.Data() << endl ;
762 //     cout <<  "    RecParticles Branch title   " << fRecParticlesTitle.Data() << endl;
763     cout <<  "    Pricipal analysis file from 0.5 to 100 " << fFileName.Data() << endl;
764     cout <<  "    Name of parameters file     "<<fFileNamePar.Data() << endl ;
765     cout <<  "    Matrix of Parameters: 18x4"<<endl;
766     cout <<  "        RCPV 6x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]"<<endl;
767     cout <<  "        TOF  6x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]"<<endl;
768     cout <<  "        PCA  5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]"<<endl;
769     cout <<  "        Energy Calibration  1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]"<<endl;
770     fParameters->Print() ;
771     cout <<  "============================================" << endl ;
772 }
773
774 //____________________________________________________________________________
775 void  AliPHOSPIDv1::WriteRecParticles(Int_t event)
776 {
777  
778   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; 
779
780   TClonesArray * recParticles = gime->RecParticles() ; 
781   recParticles->Expand(recParticles->GetEntriesFast() ) ;
782   TTree * treeR ;
783
784   if(fToSplit){
785     if(!fSplitFile)
786       return ;
787     fSplitFile->cd() ;
788     char name[10] ;
789     sprintf(name,"%s%d", "TreeR",event) ;
790     treeR = dynamic_cast<TTree*>(fSplitFile->Get(name)); 
791   }
792   else{
793     treeR = gAlice->TreeR();
794   }
795   
796   if(!treeR){
797     gAlice->MakeTree("R", fSplitFile);
798     treeR = gAlice->TreeR() ;
799   }
800   
801   //First rp
802   Int_t bufferSize = 32000 ;    
803   TBranch * rpBranch = treeR->Branch("PHOSRP",&recParticles,bufferSize);
804   rpBranch->SetTitle(BranchName());
805
806   
807   //second, pid
808   Int_t splitlevel = 0 ; 
809   AliPHOSPIDv1 * pid = this ;
810   TBranch * pidBranch = treeR->Branch("AliPHOSPID","AliPHOSPIDv1",&pid,bufferSize,splitlevel);
811   pidBranch->SetTitle(BranchName());
812   
813   rpBranch->Fill() ;
814   pidBranch->Fill() ; 
815   
816   treeR->AutoSave() ; //Write(0,kOverwrite) ;  
817   if(gAlice->TreeR()!=treeR){
818     treeR->Delete();
819   }
820 }
821
822 //____________________________________________________________________________
823 TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const 
824
825   // Calculates the momentum direction:
826   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
827   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints 
828   //  However because of the poor position resolution of PPSD the direction is always taken as if we were 
829   //  in case 1.
830
831   TVector3 dir(0,0,0) ; 
832   
833   TVector3 emcglobalpos ;
834   TMatrix  dummy ;
835   
836   emc->GetGlobalPosition(emcglobalpos, dummy) ;
837   
838
839   dir = emcglobalpos ;  
840   dir.SetZ( -dir.Z() ) ;   // why ?  
841   dir.SetMag(1.) ;
842
843   //account correction to the position of IP
844   Float_t xo,yo,zo ; //Coordinates of the origin
845   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
846   TVector3 origin(xo,yo,zo);
847   dir = dir - origin ;
848
849   return dir ;  
850 }
851 //____________________________________________________________________________
852 void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
853 {
854   // Print table of reconstructed particles
855
856   AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ; 
857
858   TClonesArray * recParticles = gime->RecParticles(BranchName()) ; 
859   
860   cout << "AliPHOSPIDv1: event "<<gAlice->GetEvNumber()  << endl ;
861   cout << "       found " << recParticles->GetEntriesFast() << " RecParticles " << endl ;
862   
863   if(strstr(option,"all")) {  // printing found TS
864     
865     cout << "  PARTICLE       "   
866          << "  Index    "  << endl ;
867     
868     Int_t index ;
869     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
870        AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;       
871
872        cout << setw(10) << rp->Name() << "  "
873             << setw(5) <<  rp->GetIndexInList() << " " <<endl;
874        cout << "Type "<<  rp->GetType() << endl;
875     }
876     cout << "-------------------------------------------" << endl ;
877   }
878   
879 }
880
881
882