Transition to NewIO
[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 > CpvEmcDistance (each bit corresponds
27 //      to a different efficiency-purity point of the photon identification) 
28 //     -Bit 3 to 5: bit set if TOF  < TimeGate (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 done 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 through the study of the reconstructed 
41 //      energy distribution of monoenergetic photons. 
42 //
43 //      All the parameters (RCPV(2 rows-3 columns),TOF(1r-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 (9,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 ellipses 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
98 // --- AliRoot header files ---
99
100 #include "AliGenerator.h"
101 #include "AliPHOS.h"
102 #include "AliPHOSPIDv1.h"
103 #include "AliPHOSClusterizerv1.h"
104 #include "AliPHOSTrackSegment.h"
105 #include "AliPHOSTrackSegmentMakerv1.h"
106 #include "AliPHOSRecParticle.h"
107 #include "AliPHOSGeometry.h"
108 #include "AliPHOSGetter.h"
109
110 ClassImp( AliPHOSPIDv1) 
111
112 //____________________________________________________________________________
113 AliPHOSPIDv1::AliPHOSPIDv1():AliPHOSPID()
114
115   // default ctor
116  
117   InitParameters() ; 
118   fDefaultInit = kTRUE ; 
119 }
120
121 //____________________________________________________________________________
122 AliPHOSPIDv1::AliPHOSPIDv1(const AliPHOSPIDv1 & pid ):AliPHOSPID(pid)
123
124   // ctor
125   InitParameters() ; 
126   Init() ;
127
128 }
129
130 //____________________________________________________________________________
131 AliPHOSPIDv1::AliPHOSPIDv1(const TString alirunFileName, const TString eventFolderName):AliPHOSPID(alirunFileName, eventFolderName)
132
133   //ctor with the indication on where to look for the track segments
134  
135   InitParameters() ; 
136   Init() ;
137   fDefaultInit = kFALSE ; 
138 }
139
140 //____________________________________________________________________________
141 AliPHOSPIDv1::~AliPHOSPIDv1()
142
143   // dtor
144
145   delete [] fX ;       // Principal input 
146   delete [] fPPhoton ; // Photon Principal components
147   delete [] fPPi0 ;    // Pi0 Principal components
148 }
149 //____________________________________________________________________________
150 const TString AliPHOSPIDv1::BranchName() const 
151 {  
152
153   return GetName() ;
154 }
155  
156 //____________________________________________________________________________
157 void AliPHOSPIDv1::Init()
158 {
159   // Make all memory allocations that are not possible in default constructor
160   // Add the PID task to the list of PHOS tasks
161
162   AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data()) ; 
163
164   if ( !gime->PID() ) 
165     gime->PostPID(this) ;
166 }
167
168 //____________________________________________________________________________
169 void AliPHOSPIDv1::InitParameters()
170 {
171   // Initialize PID parameters
172   fRecParticlesInRun = 0 ; 
173   fNEvent            = 0 ;            
174   fRecParticlesInRun = 0 ;
175   SetParameters() ; // fill the parameters matrix from parameters file
176 }
177
178 //________________________________________________________________________
179 void  AliPHOSPIDv1::Exec(Option_t * option) 
180 {
181   //Steering method
182   
183
184   if(strstr(option,"tim"))
185     gBenchmark->Start("PHOSPID");
186   
187   if(strstr(option,"print")) {
188     Print() ; 
189     return ; 
190   }
191
192
193   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
194  
195   Int_t nevents = gime->MaxEvent() ;       
196   Int_t ievent ;
197
198
199   for(ievent = 0; ievent < nevents; ievent++){
200     gime->Event(ievent,"TR") ;
201     if(gime->TrackSegments() && //Skip events, where no track segments made
202        gime->TrackSegments()->GetEntriesFast()) {
203       MakeRecParticles() ;
204       WriteRecParticles(ievent);
205       if(strstr(option,"deb"))
206         PrintRecParticles(option) ;
207       //increment the total number of rec particles per run 
208       fRecParticlesInRun += gime->RecParticles()->GetEntriesFast() ; 
209     }
210   }
211   if(strstr(option,"tim")){
212     gBenchmark->Stop("PHOSPID");
213     Info("Exec", "took %f seconds for PID %f seconds per event", 
214          gBenchmark->GetCpuTime("PHOSPID"),  
215          gBenchmark->GetCpuTime("PHOSPID")/nevents) ;
216   }
217  
218   Unload();
219 }
220
221 //____________________________________________________________________________
222 const TString AliPHOSPIDv1::GetFileNamePrincipal(TString particle) const
223 {
224   //Get file name that contains the PCA for a particle ("photon or pi0")
225   particle.ToLower();
226   TString name;
227   if      (particle=="photon") name = fFileNamePrincipalPhoton ;
228   else if (particle=="pi0"   ) name = fFileNamePrincipalPi0    ;
229   else    Error("GetFileNamePrincipal","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data());
230   return name;
231 }
232
233 //____________________________________________________________________________
234 const Float_t  AliPHOSPIDv1::GetParameterCalibration(Int_t i) const 
235 {
236   // Get the i-th parameter "Calibration"
237   Float_t param = 0.;
238   if (i>2 || i<0)
239     Error("GetParameterCalibration","Invalid parameter number: %d",i);
240   else
241     param = (*fParameters)(0,i);
242   return param;
243 }
244
245 //____________________________________________________________________________
246 const Float_t  AliPHOSPIDv1::GetCalibratedEnergy(const Float_t e) const
247 {
248 //      It calibrates Energy depending on the recpoint energy.
249 //      The energy of the reconstructed cluster is corrected with 
250 //      the formula A + B* E  + C* E^2, whose parameters where obtained 
251 //      through the study of the reconstructed energy distribution of 
252 //      monoenergetic photons.
253  
254   Float_t p[]={0.,0.,0.};
255   for (Int_t i=0; i<3; i++) p[i] = GetParameterCalibration(i);
256   Float_t enerec = p[0] +  p[1]*e + p[2]*e*e;
257   return enerec ;
258
259 }
260
261 //____________________________________________________________________________
262 const Float_t  AliPHOSPIDv1::GetParameterCpv2Emc(Int_t i, TString axis) const 
263 {
264   // Get the i-th parameter "CPV-EMC distance" for the specified axis
265   Float_t param = 0.;
266   if(i>2 || i<0)
267     Error("GetParameterCpv2Emc","Invalid parameter number: %d",i);
268   else {
269     axis.ToLower();
270     if      (axis == "x") param = (*fParameters)(1,i);
271     else if (axis == "z") param = (*fParameters)(2,i);
272     else Error("GetParameterCpv2Emc","Invalid axis name: %s",axis.Data());
273   }
274   return  param;
275 }
276
277 //____________________________________________________________________________
278 const Float_t  AliPHOSPIDv1::GetCpv2EmcDistanceCut(TString axis, Float_t e) const
279 {
280   // Get CpvtoEmcDistance Cut depending on the cluster energy, axis and 
281   // Purity-Efficiency point 
282
283   axis.ToLower();
284   Float_t p[]={0.,0.,0.};
285   for (Int_t i=0; i<3; i++) p[i] = GetParameterCpv2Emc(i,axis);
286   Float_t sig = p[0] + TMath::Exp(p[1] - p[2]*e);
287   return sig;
288 }
289
290 //____________________________________________________________________________
291 const Float_t  AliPHOSPIDv1::GetEllipseParameter(TString particle, TString param, Float_t e) const 
292 {
293   // Calculates the parameter param of the ellipse
294
295   particle.ToLower();
296   param.   ToLower();
297   Float_t p[4]={0.,0.,0.,0.};
298   Float_t value = 0.0;
299   for (Int_t i=0; i<4; i++) p[i] = GetParameterToCalculateEllipse(particle,param,i);
300   if (particle == "photon") {
301     if      (param.Contains("a"))  e = TMath::Min((Double_t)e,70.);
302     else if (param.Contains("b"))  e = TMath::Min((Double_t)e,70.);
303     else if (param.Contains("x0")) e = TMath::Max((Double_t)e,1.1);
304   }
305
306   value = p[0]/TMath::Sqrt(e) + p[1]*e + p[2]*e*e + p[3];
307   return value;
308 }
309
310 //_____________________________________________________________________________
311 const Float_t  AliPHOSPIDv1::GetParameterPhotonBoundary (Int_t i) const
312
313   // Get the parameter "i" to calculate the boundary on the moment M2x
314   // for photons at high p_T
315   Float_t param = 0;
316   if (i>3 || i<0)
317     Error("GetParameterPhotonBoundary","Wrong parameter number: %d\n",i);
318   else
319     param = (*fParameters)(14,i) ;
320   return param;
321 }
322
323 //____________________________________________________________________________
324 const Float_t  AliPHOSPIDv1::GetParameterPi0Boundary (Int_t i) const
325
326   // Get the parameter "i" to calculate the boundary on the moment M2x
327   // for pi0 at high p_T
328   Float_t param = 0;
329   if (i>2 || i<0)
330     Error("GetParameterPi0Boundary","Wrong parameter number: %d\n",i);
331   else
332     param = (*fParameters)(15,i) ;
333   return param;
334 }
335
336 //____________________________________________________________________________
337 const Float_t  AliPHOSPIDv1::GetParameterTimeGate(Int_t i) const
338 {
339   // Get TimeGate parameter depending on Purity-Efficiency i:
340   // i=0 - Low purity, i=1 - Medium purity, i=2 - High purity
341   Float_t param = 0.;
342   if(i>2 || i<0)
343     Error("GetParameterTimeGate","Invalid Efficiency-Purity choice %d",i);
344   else
345     param = (*fParameters)(3,i) ; 
346   return param;
347 }
348
349 //_____________________________________________________________________________
350 const Float_t  AliPHOSPIDv1::GetParameterToCalculateEllipse(TString particle, TString param, Int_t i) const
351
352   // Get the parameter "i" that is needed to calculate the ellipse 
353   // parameter "param" for the particle "particle" ("photon" or "pi0")
354
355   particle.ToLower();
356   param.   ToLower();
357   Int_t offset = -1;
358   if      (particle == "photon") offset=0;
359   else if (particle == "pi0")    offset=5;
360   else
361     Error("GetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data());
362
363   Int_t p= -1;
364   Float_t par = 0;
365
366   if     (param.Contains("a")) p=4+offset; 
367   else if(param.Contains("b")) p=5+offset; 
368   else if(param.Contains("c")) p=6+offset; 
369   else if(param.Contains("x0"))p=7+offset; 
370   else if(param.Contains("y0"))p=8+offset;
371
372   if      (i>4 || i<0)
373     Error("GetParameterToCalculateEllipse", "No parameter with index", i) ; 
374   else if (p==-1)
375     Error("GetParameterToCalculateEllipse", "No parameter with name %s", param.Data() ) ; 
376   else
377     par = (*fParameters)(p,i) ;
378   
379   return par;
380 }
381
382
383 //____________________________________________________________________________
384 const Float_t  AliPHOSPIDv1::GetDistance(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv, Option_t *  axis)const
385 {
386   // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
387   
388   const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; 
389   TVector3 vecEmc ;
390   TVector3 vecCpv ;
391   if(cpv){
392     emc->GetLocalPosition(vecEmc) ;
393     cpv->GetLocalPosition(vecCpv) ; 
394     if(emc->GetPHOSMod() == cpv->GetPHOSMod()){      
395       // Correct to difference in CPV and EMC position due to different distance to center.
396       // we assume, that particle moves from center
397       Float_t dCPV = geom->GetIPtoOuterCoverDistance();
398       Float_t dEMC = geom->GetIPtoCrystalSurface() ;
399       dEMC         = dEMC / dCPV ;
400       vecCpv = dEMC * vecCpv  - vecEmc ; 
401       if (axis == "X") return vecCpv.X();
402       if (axis == "Y") return vecCpv.Y();
403       if (axis == "Z") return vecCpv.Z();
404       if (axis == "R") return vecCpv.Mag();
405     }
406     return 100000000 ;
407   }
408   return 100000000 ;
409 }
410 //____________________________________________________________________________
411 const Int_t  AliPHOSPIDv1::GetCPVBit(AliPHOSEmcRecPoint * emc,AliPHOSRecPoint * cpv,const Int_t effPur, const Float_t e) const
412 {
413   if(effPur>2 || effPur<0)
414     Error("GetCPVBit","Invalid Efficiency-Purity choice %d",effPur);
415   
416   Float_t sigX = GetCpv2EmcDistanceCut("X",e);
417   Float_t sigZ = GetCpv2EmcDistanceCut("Z",e);
418   
419   Float_t deltaX = TMath::Abs(GetDistance(emc, cpv,  "X"));
420   Float_t deltaZ = TMath::Abs(GetDistance(emc, cpv,  "Z"));
421        
422   if((deltaX>sigX*(effPur+1))|(deltaZ>sigZ*(effPur+1)))
423     return 1;//Neutral
424   else
425     return 0;//Charged
426 }
427
428 //____________________________________________________________________________
429 const Int_t  AliPHOSPIDv1::GetPrincipalBit(TString particle, const Double_t* p,const Int_t effPur, const Float_t e)const
430 {
431   //Is the particle inside de PCA ellipse?
432   
433   particle.ToLower();
434   Int_t    prinbit  = 0 ;
435   Float_t a        = GetEllipseParameter(particle,"a" , e); 
436   Float_t b        = GetEllipseParameter(particle,"b" , e);
437   Float_t c        = GetEllipseParameter(particle,"c" , e);
438   Float_t x0 = GetEllipseParameter(particle,"x0", e); 
439   Float_t y0 = GetEllipseParameter(particle,"y0", e);
440   
441   Float_t r = TMath::Power((p[0] - x0)/a,2) + 
442               TMath::Power((p[1] - y0)/b,2) +
443             c*(p[0] -  x0)*(p[1] - y0)/(a*b) ;
444   //3 different ellipses defined
445   if((effPur==2) && (r<1./2.)) prinbit= 1;
446   if((effPur==1) && (r<2.   )) prinbit= 1;
447   if((effPur==0) && (r<9./2.)) prinbit= 1;
448
449   if(r<0)
450     Error("GetPrincipalBit", "Negative square?") ;
451
452   return prinbit;
453
454 }
455 //____________________________________________________________________________
456 const Int_t  AliPHOSPIDv1::GetHardPhotonBit(AliPHOSEmcRecPoint * emc) const
457 {
458   // Set bit for identified hard photons (E > 30 GeV)
459   // if the second moment M2x is below the boundary
460
461   Float_t e   = emc->GetEnergy();
462   if (e < 30.0) return 0;
463   Float_t m2x = emc->GetM2x();
464   Float_t m2xBoundary = GetParameterPhotonBoundary(0) *
465     TMath::Exp(-TMath::Power(e-GetParameterPhotonBoundary(1),2)/2.0/
466                 TMath::Power(GetParameterPhotonBoundary(2),2)) +
467     GetParameterPhotonBoundary(3);
468   Info("GetHardPhotonBit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary);
469   if (m2x < m2xBoundary)
470     return 1;// A hard photon
471   else
472     return 0;// Not a hard photon
473 }
474
475 //____________________________________________________________________________
476 const Int_t  AliPHOSPIDv1::GetHardPi0Bit(AliPHOSEmcRecPoint * emc) const
477 {
478   // Set bit for identified hard pi0  (E > 30 GeV)
479   // if the second moment M2x is above the boundary
480
481   Float_t e   = emc->GetEnergy();
482   if (e < 30.0) return 0;
483   Float_t m2x = emc->GetM2x();
484   Float_t m2xBoundary = GetParameterPi0Boundary(0) +
485                     e * GetParameterPi0Boundary(1);
486   Info("GetHardPi0Bit","E=%f, m2x=%f, boundary=%f",e,m2x,m2xBoundary);
487   if (m2x > m2xBoundary)
488     return 1;// A hard pi0
489   else
490     return 0;// Not a hard pi0
491 }
492
493 //____________________________________________________________________________
494 TVector3 AliPHOSPIDv1::GetMomentumDirection(AliPHOSEmcRecPoint * emc, AliPHOSRecPoint * cpv)const 
495
496   // Calculates the momentum direction:
497   //   1. if only a EMC RecPoint, direction is given by IP and this RecPoint
498   //   2. if a EMC RecPoint and CPV RecPoint, direction is given by the line through the 2 recpoints 
499   //  However because of the poor position resolution of PPSD the direction is always taken as if we were 
500   //  in case 1.
501
502   TVector3 dir(0,0,0) ; 
503   
504   TVector3 emcglobalpos ;
505   TMatrix  dummy ;
506   
507   emc->GetGlobalPosition(emcglobalpos, dummy) ;
508   
509
510   dir = emcglobalpos ;  
511   dir.SetZ( -dir.Z() ) ;   // why ?  
512   dir.SetMag(1.) ;
513
514   //account correction to the position of IP
515   Float_t xo,yo,zo ; //Coordinates of the origin
516   gAlice->Generator()->GetOrigin(xo,yo,zo) ;
517   TVector3 origin(xo,yo,zo);
518   dir = dir - origin ;
519
520   return dir ;  
521 }
522
523 //____________________________________________________________________________
524 void  AliPHOSPIDv1::MakeRecParticles()
525 {
526   // Makes a RecParticle out of a TrackSegment
527   
528   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
529   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
530   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
531   TClonesArray * trackSegments = gime->TrackSegments() ; 
532   if ( !emcRecPoints || !cpvRecPoints || !trackSegments ) {
533     Fatal("MakeRecParticles", "RecPoints or TrackSegments not found !") ;  
534   }
535   TClonesArray * recParticles  = gime->RecParticles() ; 
536   recParticles->Clear();
537
538   TIter next(trackSegments) ; 
539   AliPHOSTrackSegment * ts ; 
540   Int_t index = 0 ; 
541   AliPHOSRecParticle * rp ; 
542   while ( (ts = (AliPHOSTrackSegment *)next()) ) {
543     
544     new( (*recParticles)[index] ) AliPHOSRecParticle() ;
545     rp = (AliPHOSRecParticle *)recParticles->At(index) ; 
546     rp->SetTrackSegment(index) ;
547     rp->SetIndexInList(index) ;
548         
549     AliPHOSEmcRecPoint * emc = 0 ;
550     if(ts->GetEmcIndex()>=0)
551       emc = (AliPHOSEmcRecPoint *) emcRecPoints->At(ts->GetEmcIndex()) ;
552     
553     AliPHOSRecPoint    * cpv = 0 ;
554     if(ts->GetCpvIndex()>=0)
555       cpv = (AliPHOSRecPoint *)   cpvRecPoints->At(ts->GetCpvIndex()) ;
556     
557     // Now set type (reconstructed) of the particle
558
559     // Choose the cluster energy range
560     
561     if (!emc) {
562       Fatal("MakeRecParticles", "-> emc(%d) = %d", ts->GetEmcIndex(), emc ) ;
563     }
564
565     Float_t e = emc->GetEnergy() ;   
566     
567     Float_t  lambda[2] ;
568     emc->GetElipsAxis(lambda) ;
569     
570     if((lambda[0]>0.01) && (lambda[1]>0.01)){
571       // Looking PCA. Define and calculate the data (X),
572       // introduce in the function X2P that gives the components (P).  
573
574       Float_t  Spher = 0. ;
575       Float_t  Emaxdtotal = 0. ; 
576       
577       if((lambda[0]+lambda[1])!=0) 
578         Spher=fabs(lambda[0]-lambda[1])/(lambda[0]+lambda[1]); 
579       
580       Emaxdtotal=emc->GetMaximalEnergy()/emc->GetEnergy(); 
581       
582       fX[0] = lambda[0] ;  
583       fX[1] = lambda[1] ; 
584       fX[2] = emc->GetDispersion() ; 
585       fX[3] = Spher ; 
586       fX[4] = emc->GetMultiplicity() ;  
587       fX[5] = Emaxdtotal ;  
588       fX[6] = emc->GetCoreEnergy() ;  
589       
590       fPrincipalPhoton->X2P(fX,fPPhoton);
591       fPrincipalPi0   ->X2P(fX,fPPi0);
592
593     }
594     else{
595       fPPhoton[0]=-100.0;  //We do not accept clusters with 
596       fPPhoton[1]=-100.0;  //one cell as a photon-like
597       fPPi0[0]   =-100.0;
598       fPPi0[1]   =-100.0;
599     }
600     
601     Float_t time =emc->GetTime() ;
602     
603     // Loop of Efficiency-Purity (the 3 points of purity or efficiency 
604     // are taken into account to set the particle identification)
605     for(Int_t effPur = 0; effPur < 3 ; effPur++){
606       
607       // Looking at the CPV detector. If RCPV greater than CpvEmcDistance, 
608       // 1st,2nd or 3rd bit (depending on the efficiency-purity point )
609       // is set to 1
610       if(GetCPVBit(emc, cpv, effPur,e) == 1 )  
611         rp->SetPIDBit(effPur) ;
612       
613       // Looking the TOF. If TOF smaller than gate,  4th, 5th or 6th 
614       // bit (depending on the efficiency-purity point )is set to 1             
615       if(time< (*fParameters)(2,effPur)) 
616         rp->SetPIDBit(effPur+3) ;                   
617       
618       //Photon PCA
619       //If we are inside the ellipse, 7th, 8th or 9th 
620       // bit (depending on the efficiency-purity point )is set to 1 
621       if(GetPrincipalBit("photon",fPPhoton,effPur,e) == 1) 
622         rp->SetPIDBit(effPur+6) ;
623
624       //Pi0 PCA
625       //If we are inside the ellipse, 10th, 11th or 12th 
626       // bit (depending on the efficiency-purity point )is set to 1 
627       if(GetPrincipalBit("pi0"   ,fPPi0   ,effPur,e) == 1) 
628         rp->SetPIDBit(effPur+9) ;
629     }
630     if(GetHardPhotonBit(emc))
631       rp->SetPIDBit(12) ;
632     if(GetHardPi0Bit   (emc))
633       rp->SetPIDBit(13) ;
634     
635     //Set momentum, energy and other parameters 
636     Float_t  encal = GetCalibratedEnergy(e);
637     TVector3 dir   = GetMomentumDirection(emc,cpv) ; 
638     dir.SetMag(encal) ;
639     rp->SetMomentum(dir.X(),dir.Y(),dir.Z(),encal) ;
640     rp->SetCalcMass(0);
641     rp->Name(); //If photon sets the particle pdg name to gamma
642     rp->SetProductionVertex(0,0,0,0);
643     rp->SetFirstMother(-1);
644     rp->SetLastMother(-1);
645     rp->SetFirstDaughter(-1);
646     rp->SetLastDaughter(-1);
647     rp->SetPolarisation(0,0,0);
648     index++ ; 
649   }
650 }
651   
652 //____________________________________________________________________________
653 void  AliPHOSPIDv1::Print() const
654 {
655   // Print the parameters used for the particle type identification
656
657     Info("Print", "=============== AliPHOSPIDv1 ================") ;
658     printf("Making PID\n") ;
659     printf("    Pricipal analysis file from 0.5 to 100 %s\n", fFileNamePrincipalPhoton.Data() )   ; 
660     printf("    Name of parameters file     %s\n", fFileNameParameters.Data() )  ;
661     printf("    Matrix of Parameters: 14x4\n") ;
662     printf("        Energy Calibration  1x3 [3 parametres to calibrate energy: A + B* E + C * E^2]\n") ;
663     printf("        RCPV 2x3 rows x and z, columns function cut parameters\n") ;
664     printf("        TOF  1x3 [High Eff-Low Pur,Medium Eff-Pur, Low Eff-High Pur]\n") ;
665     printf("        PCA  5x4 [5 ellipse parametres and 4 parametres to calculate them: A/Sqrt(E) + B* E + C * E^2 + D]\n") ;
666     Printf("    Pi0 PCA  5x3 [5 ellipse parametres and 3 parametres to calculate them: A + B* E + C * E^2]\n") ;
667     fParameters->Print() ;
668 }
669
670
671
672 //____________________________________________________________________________
673 void AliPHOSPIDv1::PrintRecParticles(Option_t * option)
674 {
675   // Print table of reconstructed particles
676
677   AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
678
679   TClonesArray * recParticles = gime->RecParticles() ; 
680
681   TString message ; 
682   message  = "\nevent " ;
683   message += gAlice->GetEvNumber() ; 
684   message += "       found " ; 
685   message += recParticles->GetEntriesFast(); 
686   message += " RecParticles\n" ; 
687
688   if(strstr(option,"all")) {  // printing found TS
689     message += "\n  PARTICLE         Index    \n" ; 
690     
691     Int_t index ;
692     for (index = 0 ; index < recParticles->GetEntries() ; index++) {
693       AliPHOSRecParticle * rp = (AliPHOSRecParticle * ) recParticles->At(index) ;       
694       message += "\n" ;
695       message += rp->Name().Data() ;  
696       message += " " ;
697       message += rp->GetIndexInList() ;  
698       message += " " ;
699       message += rp->GetType()  ;
700     }
701   }
702   Info("Print", message.Data() ) ; 
703 }
704
705 //____________________________________________________________________________
706 void  AliPHOSPIDv1::SetParameters() 
707 {
708   // PCA : To do the Principal Components Analysis it is necessary 
709   // the Principal file, which is opened here
710   fX       = new double[7]; // Data for the PCA 
711   fPPhoton = new double[7]; // Eigenvalues of the PCA
712   fPPi0    = new double[7]; // Eigenvalues of the Pi0 PCA
713
714   // Read photon principals from the photon file
715   
716   fFileNamePrincipalPhoton = "$ALICE_ROOT/PHOS/PCA8pa15_0.5-100.root" ; 
717   TFile f( fFileNamePrincipalPhoton.Data(), "read" ) ;
718   fPrincipalPhoton = dynamic_cast<TPrincipal*> (f.Get("principal")) ; 
719   f.Close() ; 
720
721   // Read pi0 principals from the pi0 file
722
723   fFileNamePrincipalPi0    = "$ALICE_ROOT/PHOS/PCA_pi0_40-120.root" ;
724   TFile fPi0( fFileNamePrincipalPi0.Data(), "read" ) ;
725   fPrincipalPi0    = dynamic_cast<TPrincipal*> (fPi0.Get("principal")) ; 
726   fPi0.Close() ;
727
728   // Open parameters file and initialization of the Parameters matrix. 
729   // In the File Parameters.dat are all the parameters. These are introduced 
730   // in a matrix of 16x4  
731   // 
732   // All the parameters defined in this file are, in order of row: 
733   // line   0   : calibration 
734   // lines  1,2 : CPV rectangular cat for X and Z
735   // line   3   : TOF cut
736   // lines  4-8 : parameters to calculate photon PCA ellipse
737   // lines  9-13: parameters to calculate pi0 PCA ellipse
738   // lines 14-15: parameters to calculate border for high-pt photons and pi0
739
740   fFileNameParameters = gSystem->ExpandPathName("$ALICE_ROOT/PHOS/Parameters.dat");
741   fParameters = new TMatrix(16,4) ;
742   const Int_t maxLeng=255;
743   char string[maxLeng];
744
745   // Open a text file with PID parameters
746   FILE *fd = fopen(fFileNameParameters.Data(),"r");
747   if (!fd)
748     Fatal("SetParameter","File %s with a PID parameters cannot be opened\n",
749           fFileNameParameters.Data());
750
751   Int_t i=0;
752   // Read parameter file line-by-line and skip empty line and comments
753   while (fgets(string,maxLeng,fd) != NULL) {
754     if (string[0] == '\n' ) continue;
755     if (string[0] == '!'  ) continue;
756     sscanf(string, "%f %f %f %f",
757            &(*fParameters)(i,0), &(*fParameters)(i,1), 
758            &(*fParameters)(i,2), &(*fParameters)(i,3));
759     i++;
760     //printf("line %d: %s",i,string);
761   }
762   fclose(fd);
763 }
764
765 //____________________________________________________________________________
766 void  AliPHOSPIDv1::SetParameterCalibration(Int_t i,Float_t param) 
767 {
768   // Set parameter "Calibration" i to a value param
769   if(i>2 || i<0)
770     Error("SetParameterCalibration","Invalid parameter number: %d",i);
771   else
772     (*fParameters)(0,i) = param ;
773 }
774
775 //____________________________________________________________________________
776 void  AliPHOSPIDv1::SetParameterCpv2Emc(Int_t i, TString axis, Float_t cut) 
777 {
778   // Set the parameters to calculate Cpv-to-Emc Distance Cut depending on 
779   // Purity-Efficiency point i
780
781   if(i>2 || i<0)
782     Error("SetParameterCpv2Emc","Invalid parameter number: %d",i);
783   else {
784     axis.ToLower();
785     if      (axis == "x") (*fParameters)(1,i) = cut;
786     else if (axis == "z") (*fParameters)(2,i) = cut;
787     else Error("SetParameterCpv2Emc","Invalid axis name: %s",axis.Data());
788   }
789 }
790
791 //____________________________________________________________________________
792 void  AliPHOSPIDv1::SetParameterPhotonBoundary(Int_t i,Float_t param) 
793 {
794   // Set parameter "Hard photon boundary" i to a value param
795   if(i>4 || i<0)
796     Error("SetParameterPhotonBoundary","Invalid parameter number: %d",i);
797   else
798     (*fParameters)(14,i) = param ;
799 }
800
801 //____________________________________________________________________________
802 void  AliPHOSPIDv1::SetParameterPi0Boundary(Int_t i,Float_t param) 
803 {
804   // Set parameter "Hard pi0 boundary" i to a value param
805   if(i>1 || i<0)
806     Error("SetParameterPi0Boundary","Invalid parameter number: %d",i);
807   else
808     (*fParameters)(15,i) = param ;
809 }
810
811 //_____________________________________________________________________________
812 void  AliPHOSPIDv1::SetParameterTimeGate(Int_t i, Float_t gate) 
813 {
814   // Set the parameter TimeGate depending on Purity-Efficiency point i 
815   if (i>2 || i<0)
816     Error("SetParameterTimeGate","Invalid Efficiency-Purity choice %d",i);
817   else
818     (*fParameters)(3,i)= gate ; 
819
820
821 //_____________________________________________________________________________
822 void  AliPHOSPIDv1::SetParameterToCalculateEllipse(TString particle, TString param, Int_t i, Float_t par) 
823 {  
824   // Set the parameter "i" that is needed to calculate the ellipse 
825   // parameter "param" for a particle "particle"
826   
827   particle.ToLower();
828   param.   ToLower();
829   Int_t p= -1;
830   Int_t offset=0;
831
832   if      (particle == "photon") offset=0;
833   else if (particle == "pi0")    offset=5;
834   else
835     Error("SetParameterToCalculateEllipse","Wrong particle name: %s (choose from pi0/photon)\n",particle.Data());
836
837   if     (param.Contains("a")) p=4+offset; 
838   else if(param.Contains("b")) p=5+offset; 
839   else if(param.Contains("c")) p=6+offset; 
840   else if(param.Contains("x0"))p=7+offset; 
841   else if(param.Contains("y0"))p=8+offset;
842   if((i>4)||(i<0))
843     Error("SetEllipseParameter", "No parameter with index %d", i) ; 
844   else if(p==-1)
845     Error("SetEllipseParameter", "No parameter with name %s", param.Data() ) ; 
846   else
847     (*fParameters)(p,i) = par ;
848
849
850 //____________________________________________________________________________
851 void AliPHOSPIDv1::Unload() 
852 {
853   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
854   gime->PhosLoader()->UnloadRecPoints() ;
855   gime->PhosLoader()->UnloadTracks() ;
856   gime->PhosLoader()->UnloadRecParticles() ;
857 }
858
859 //____________________________________________________________________________
860 void  AliPHOSPIDv1::WriteRecParticles(Int_t event)
861 {
862  
863   AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
864
865   TClonesArray * recParticles = gime->RecParticles() ; 
866   recParticles->Expand(recParticles->GetEntriesFast() ) ;
867   TTree * treeP =  gime->TreeP();
868   
869   //First rp
870   Int_t bufferSize = 32000 ;
871   TBranch * rpBranch = treeP->Branch("PHOSRP",&recParticles,bufferSize);
872   rpBranch->SetTitle(BranchName());
873   
874   rpBranch->Fill() ;
875  
876   gime->WriteRecParticles("OVERWRITE");
877   gime->WritePID("OVERWRITE");
878 }
879