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