New classes: AliPHOSRecParticle, AliPHOSParticleGuesser, AliPHOSAnalyze
[u/mrichter/AliRoot.git] / PHOS / AliPHOSAnalyze.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 //_________________________________________________________________________
17 // Algorythm class to analyze PHOS events
18 //*-- Y. Schutz :   SUBATECH 
19 //////////////////////////////////////////////////////////////////////////////
20
21 // --- ROOT system ---
22
23 #include "TH2.h"
24 #include "TParticle.h"
25 #include "TClonesArray.h"
26 #include "TTree.h"
27 #include "TMath.h"
28 #include "TCanvas.h" 
29
30 // --- Standard library ---
31
32 // --- AliRoot header files ---
33
34 #include "AliRun.h"
35 #include "AliPHOSAnalyze.h"
36 #include "AliPHOSClusterizerv1.h"
37 #include "AliPHOSTrackSegmentMakerv1.h"
38 #include "AliPHOSParticleGuesserv1.h"
39 #include "AliPHOSReconstructioner.h"
40 #include "AliPHOSDigit.h"
41 #include "AliPHOSTrackSegment.h"
42 #include "AliPHOSRecParticle.h"
43
44 ClassImp(AliPHOSAnalyze)
45
46
47 //____________________________________________________________________________
48   AliPHOSAnalyze::AliPHOSAnalyze()
49 {
50   // ctor
51   
52   fRootFile = 0 ; 
53 }
54
55 //____________________________________________________________________________
56 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
57 {
58   // ctor
59   
60   Bool_t OK = OpenRootFile(name)  ; 
61   if ( !OK ) {
62     cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
63   }
64   else { 
65     gAlice = (AliRun*) fRootFile->Get("gAlice");
66     fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
67     fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
68     fEvt = -999 ; 
69   }
70 }
71
72 //____________________________________________________________________________
73 AliPHOSAnalyze::~AliPHOSAnalyze()
74 {
75   // dtor
76
77   fRootFile->Close() ; 
78   delete fRootFile ; 
79   fRootFile = 0 ; 
80
81   delete fPHOS ; 
82   fPHOS = 0 ; 
83
84   delete fClu ; 
85   fClu = 0 ;
86
87   delete fPag ; 
88   fPag = 0 ;
89
90   delete fRec ; 
91   fRec = 0 ;
92
93   delete fTrs ; 
94   fTrs = 0 ;
95
96 }
97
98 //____________________________________________________________________________
99 void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
100 {
101   Bool_t OK = Init(evt) ; 
102   
103   if ( OK ) {
104     //=========== Get the number of entries in the Digits array
105     
106     Int_t nId = fPHOS->Digits()->GetEntries();    
107     printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
108     
109     //=========== Do the reconstruction
110     
111     cout << "AnalyzeOneEvent > Found  " << nId << "  digits in PHOS"   << endl ;  
112     
113     fPHOS->Reconstruction(fRec);  
114     
115     // =========== End of reconstruction
116     
117     cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;   
118   } // OK
119   else
120     cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;   
121
122 }
123
124 //____________________________________________________________________________
125 Bool_t AliPHOSAnalyze::Init(Int_t evt)
126 {
127
128   Bool_t OK = kTRUE ; 
129   
130    //========== Open galice root file  
131
132   if ( fRootFile == 0 ) {
133     Text_t * name  = new Text_t[80] ; 
134     cout << "AnalyzeOneEvent > Enter file root file name : " ;  
135     cin >> name ; 
136     Bool_t OK = OpenRootFile(name) ; 
137     if ( !OK )
138       cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
139     else { 
140       //========== Get AliRun object from file 
141       
142       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
143       
144       //=========== Get the PHOS object and associated geometry from the file 
145       
146       fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
147       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
148     } // else !OK
149   } // if fRootFile
150   
151   if ( OK ) {
152     
153     //========== Create the Clusterizer
154
155     fClu = new AliPHOSClusterizerv1() ; 
156     fClu->SetEmcEnergyThreshold(0.01) ; 
157     fClu->SetEmcClusteringThreshold(0.1) ; 
158     fClu->SetPpsdEnergyThreshold(0.0000001) ; 
159     fClu->SetPpsdClusteringThreshold(0.0000002) ; 
160     fClu->SetLocalMaxCut(0.03) ;
161     fClu->SetCalibrationParameters(0., 0.0000001) ;  
162     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ; 
163     fClu->PrintParameters() ; 
164     
165     //========== Creates the track segment maker
166     
167     fTrs = new AliPHOSTrackSegmentMakerv1() ;
168     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; 
169     
170     //========== Creates the particle guesser
171     
172     fPag = new AliPHOSParticleGuesserv1() ;
173     cout <<  "AnalyzeOneEvent > using particle guess " << fPag->GetName() << endl ; 
174     
175     //========== Creates the Reconstructioner  
176     
177     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPag) ;     
178     
179     //=========== Connect the various Tree's for evt
180     
181     if ( evt == -999 ) {
182       cout <<  "AnalyzeOneEvent > Enter event number : " ; 
183       cin >> evt ;  
184       cout << evt << endl ; 
185     }
186     fEvt = evt ; 
187     gAlice->GetEvent(evt);
188     
189     //=========== Get the Digit TTree
190     
191     gAlice->TreeD()->GetEvent(0) ;     
192     
193   } // OK
194   
195   return OK ; 
196 }
197
198 //____________________________________________________________________________
199 void AliPHOSAnalyze:: DisplayKineEvent(Int_t evt)
200 {
201   if (evt == -999) 
202     evt = fEvt ;
203
204   Int_t Module ; 
205   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ; 
206   cin >> Module ; cout << Module << endl ; 
207
208   Int_t TestParticle ; 
209   cout << " 22      : PHOTON " << endl 
210        << " (-)11   : (POSITRON)ELECTRON " << endl 
211        << " (-)2112 : (ANTI)NEUTRON " << endl  
212        << " -999    : Everything else " << endl ; 
213   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ; 
214   cin >> TestParticle ; cout << TestParticle << endl ; 
215
216   Text_t HistoName[80] ;
217   sprintf(HistoName,"Event %d: Incident particles in module %d", evt, Module) ; 
218
219   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by Module   
220   fGeom->EmcModuleCoverage(Module, tm, tM, pm, pM, kDegre) ;
221
222   Double_t theta, phi ; 
223   fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
224
225   Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
226   Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
227
228   tm -= theta ; 
229   tM += theta ; 
230   pm -= phi ; 
231   pM += phi ; 
232
233   TH2F * HistoParticle = new TH2F("HistoParticle",  HistoName, 
234                                           pdim, pm, pM, tdim, tm, tM) ; 
235   HistoParticle->SetStats(kFALSE) ;
236
237   // Get pointers to Alice Particle TClonesArray
238
239   TParticle * Particle;
240   TClonesArray * ArrayOfParticles  = gAlice->Particles();    
241
242   Text_t CanvasName[80];
243   sprintf(CanvasName,"Particles incident in PHOS/EMC module # %d",Module) ;
244   TCanvas * KineCanvas = new TCanvas("KineCanvas", CanvasName, 650, 500) ; 
245
246   // get the KINE Tree
247
248   TTree * Kine =  gAlice->TreeK() ; 
249   Stat_t NumberOfParticles =  Kine->GetEntries() ; 
250   cout << "DisplayKineEvent > events in Kine " << NumberOfParticles << endl ; 
251   
252   // loop over particles
253
254   Double_t raddeg = 180. / TMath::Pi() ; 
255   Int_t index1 ; 
256   Int_t nparticlein = 0 ; 
257   for (index1 = 0 ; index1 < NumberOfParticles ; index1++){
258     Int_t nparticle = ArrayOfParticles->GetEntriesFast() ;
259     Int_t index2 ; 
260     for( index2 = 0 ; index2 < nparticle ; index2++) {         
261       Particle            = (TParticle*)ArrayOfParticles->UncheckedAt(index2) ;
262       Int_t  ParticleType = Particle->GetPdgCode() ;
263       if (TestParticle == -999 || TestParticle == ParticleType) { 
264         Double_t Phi        = Particle->Phi() ;
265         Double_t Theta      = Particle->Theta() ;
266         Int_t mod ; 
267         Double_t x, z ; 
268         fGeom->ImpactOnEmc(Theta, Phi, mod, z, x) ;
269         if ( mod == Module ) {
270           nparticlein++ ; 
271           HistoParticle->Fill(Phi*raddeg, Theta*raddeg, Particle->Energy() ) ; 
272         } 
273       } 
274     }
275   }
276   KineCanvas->Draw() ; 
277   HistoParticle->Draw("color") ; 
278   TPaveText *  PaveText = new TPaveText(294, 100, 300, 101); 
279   Text_t text[40] ; 
280   sprintf(text, "Particles: %d ", nparticlein) ;
281   PaveText->AddText(text) ; 
282   PaveText->Draw() ; 
283   KineCanvas->Update(); 
284
285 }
286 //____________________________________________________________________________
287 void AliPHOSAnalyze::DisplayRecParticles()
288 {
289   if (fEvt == -999) {
290     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
291     Text_t answer[1] ; 
292     cin >> answer ; cout << answer ; 
293     if ( answer == "y" ) 
294       AnalyzeOneEvent() ;
295   } 
296     if (fEvt != -999) {
297       
298       Int_t Module ; 
299       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ; 
300       cin >> Module ; cout << Module << endl ;
301       Text_t HistoName[80] ; 
302       sprintf(HistoName,"Event %d: Reconstructed particles in module %d", fEvt, Module) ; 
303       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by Module   
304       fGeom->EmcModuleCoverage(Module, tm, tM, pm, pM, kDegre) ;
305       Double_t theta, phi ; 
306       fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
307       Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
308       Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
309       tm -= theta ; 
310       tM += theta ; 
311       pm -= phi ; 
312       TH2F * HistoRParticle = new TH2F("HistoRParticle",  HistoName, 
313                                        pdim, pm, pM, tdim, tm, tM) ; 
314       HistoRParticle->SetStats(kFALSE) ;
315       Text_t CanvasName[80] ; 
316       sprintf(CanvasName, "Reconstructed particles in PHOSmodule # %d", Module) ;
317       TCanvas * RParticleCanvas = new TCanvas("RparticleCanvas", CanvasName, 650, 500) ; 
318       RecParticlesList * rpl = fPHOS->RecParticles() ; 
319       Int_t NRecParticles = rpl->GetEntries() ; 
320       Int_t NRecParticlesInModule = 0 ; 
321       TIter nextRecPart(rpl) ; 
322       AliPHOSRecParticle * rp ; 
323       cout << "DisplayRecParticles > " << NRecParticles << " reconstructed particles " << endl ; 
324       Double_t raddeg = 180. / TMath::Pi() ; 
325       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
326         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
327         if ( ts->GetPHOSMod() == Module ) {  
328           NRecParticlesInModule++ ; 
329           Double_t theta = rp->Theta() * raddeg ;
330           Double_t phi   = rp->Phi() * raddeg ;
331           Double_t energy = rp->Energy() ; 
332           HistoRParticle->Fill(phi, theta, energy) ;
333         }
334       }
335       HistoRParticle->Draw("color") ; 
336       Text_t text[80] ; 
337       sprintf(text, "reconstructed particles: %d", NRecParticlesInModule) ;
338       TPaveText *  PaveText = new TPaveText(292, 100, 300, 101); 
339       PaveText->AddText(text) ; 
340       PaveText->Draw() ; 
341       RParticleCanvas->Update() ; 
342     }
343 }
344
345 //____________________________________________________________________________
346 void AliPHOSAnalyze::DisplayRecPoints()
347 {
348   if (fEvt == -999) {
349     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
350     Text_t answer[1] ; 
351     cin >> answer ; cout << answer ; 
352     if ( answer == "y" ) 
353       AnalyzeOneEvent() ;
354   } 
355     if (fEvt != -999) {
356       
357       Int_t Module ; 
358       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ; 
359       cin >> Module ; cout << Module << endl ; 
360
361       Text_t CanvasName[80];
362       sprintf(CanvasName,"Digits in PHOS/EMC module # %d",Module) ;
363       TCanvas * ModuleCanvas = new TCanvas("Module", CanvasName, 650, 500) ; 
364       ModuleCanvas->Draw() ;
365
366       //=========== Creating 2d-histogram of the PHOS Module
367       // a little bit junkie but is used to test Geom functinalities
368
369       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by Module   
370       
371       fGeom->EmcModuleCoverage(Module, tm, tM, pm, pM); 
372       // convert angles into coordinates local to the EMC module of interest
373
374       Int_t EmcModuleNumber ;
375       Double_t EmcModulexm, EmcModulezm ; // minimum local coordinate in a given EMCA module
376       Double_t EmcModulexM, EmcModulezM ; // maximum local coordinate in a given EMCA module
377       fGeom->ImpactOnEmc(tm, pm, EmcModuleNumber, EmcModulezm, EmcModulexm) ;
378       fGeom->ImpactOnEmc(tM, pM, EmcModuleNumber, EmcModulezM, EmcModulexM) ;
379       Int_t xdim = (Int_t)( ( EmcModulexM - EmcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
380       Int_t zdim = (Int_t)( ( EmcModulezM - EmcModulezm ) / fGeom->GetCrystalSize(2) ) ;
381       Float_t xmin = EmcModulexm - fGeom->GetCrystalSize(0) ; 
382       Float_t xMax = EmcModulexM + fGeom->GetCrystalSize(0) ; 
383       Float_t zmin = EmcModulezm - fGeom->GetCrystalSize(2) ; 
384       Float_t zMax = EmcModulezM + fGeom->GetCrystalSize(2) ;     
385       Text_t HistoName[80];
386       sprintf(HistoName,"Event %d: Digits and RecPoints in module %d", fEvt, Module) ;
387       TH2F * hModule = new TH2F("HistoReconstructed", HistoName,
388                                 xdim, xmin, xMax, zdim, zmin, zMax) ;  
389       hModule->SetMaximum(2.0);
390       hModule->SetMinimum(0.0);
391       hModule->SetStats(kFALSE); 
392
393       TIter next(fPHOS->Digits()) ;
394       Float_t Energy, y, z;
395       Int_t RelId[4]; Int_t NumberOfDigits = 0 ;
396       AliPHOSDigit * digit ; 
397       Float_t Etot ; 
398       while((digit = (AliPHOSDigit *)next())) 
399         {  
400           fGeom->AbsToRelNumbering(digit->GetId(), RelId) ;
401           if (RelId[0] == Module)  
402             {  
403               NumberOfDigits++ ;
404               Energy = fClu->Calibrate(digit->GetAmp()) ;
405               Etot += Energy ; 
406               fGeom->RelPosInModule(RelId,y,z) ; 
407               if (Energy > 0.01 )  
408                 hModule->Fill(y, z, Energy) ;
409             } 
410         }
411       cout <<"DrawRecPoints >  Found in Module " 
412            << Module << " " << NumberOfDigits << "  digits with total energy " << Etot << endl ;
413       hModule->Draw("col2") ;
414
415       //=========== Cluster in Module
416
417       TClonesArray * EmcRP = fPHOS->EmcClusters() ;
418       Etot = 0.; 
419       Int_t TotalNumberOfClusters = 0 ; 
420       Int_t NumberOfClusters = 0 ;
421       TIter nextemc(EmcRP) ;
422       AliPHOSEmcRecPoint * emc ;
423       while((emc = (AliPHOSEmcRecPoint *)nextemc())) 
424         {
425           TotalNumberOfClusters++ ;
426           if ( emc->GetPHOSMod() == Module )
427             { 
428               NumberOfClusters++ ; 
429               Energy = emc->GetTotalEnergy() ;   
430               Etot+= Energy ;  
431               emc->Draw("P") ;
432             }
433         }
434       cout << "DrawRecPoints > Found " << TotalNumberOfClusters << " EMC Clusters in PHOS" << endl ; 
435       cout << "DrawRecPoints > Found in Module " << Module << "  " << NumberOfClusters << " EMC Clusters " << endl ;
436       cout << "DrawRecPoints > Total energy  " << Etot << endl ; 
437
438       TPaveText *  PaveText = new TPaveText(22, 80, 83, 90); 
439       Text_t text[40] ; 
440       sprintf(text, "digits: %d;  clusters: %d", NumberOfDigits, NumberOfClusters) ;
441       PaveText->AddText(text) ; 
442       PaveText->Draw() ; 
443       ModuleCanvas->Update(); 
444  
445       //=========== Cluster in Module PPSD Down
446
447       TClonesArray * PpsdRP = fPHOS->PpsdClusters() ;
448       Etot = 0.; 
449       TIter nextPpsd(PpsdRP) ;
450       AliPHOSPpsdRecPoint * Ppsd ;
451       while((Ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
452         {
453           TotalNumberOfClusters++ ;
454           if ( Ppsd->GetPHOSMod() == Module )
455             { 
456               NumberOfClusters++ ; 
457               Energy = Ppsd->GetEnergy() ;   
458               Etot+=Energy ;  
459               if (!Ppsd->GetUp()) Ppsd->Draw("P") ;
460             }
461         }
462       cout << "DrawRecPoints > Found " << TotalNumberOfClusters << " Ppsd Down Clusters in PHOS" << endl ; 
463       cout << "DrawRecPoints > Found in Module " << Module << "  " << NumberOfClusters << " Ppsd Down Clusters " << endl ;
464       cout << "DrawRecPoints > Total energy  " << Etot << endl ; 
465
466       //=========== Cluster in Module PPSD Up
467   
468       PpsdRP = fPHOS->PpsdClusters() ;
469       Etot = 0.; 
470       TIter nextPpsdUp(PpsdRP) ;
471       while((Ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) 
472         {
473           TotalNumberOfClusters++ ;
474           if ( Ppsd->GetPHOSMod() == Module )
475             { 
476               NumberOfClusters++ ; 
477               Energy = Ppsd->GetEnergy() ;   
478               Etot+=Energy ;  
479               if (Ppsd->GetUp()) Ppsd->Draw("P") ;
480             }
481         }
482   cout << "DrawRecPoints > Found " << TotalNumberOfClusters << " Ppsd Up Clusters in PHOS" << endl ; 
483   cout << "DrawRecPoints > Found in Module " << Module << "  " << NumberOfClusters << " Ppsd Up Clusters " << endl ;
484   cout << "DrawRecPoints > Total energy  " << Etot << endl ; 
485     
486     } // if !-999
487 }
488
489 //____________________________________________________________________________
490 void AliPHOSAnalyze::DisplayTrackSegments()
491 {
492   if (fEvt == -999) {
493     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; 
494     Text_t answer[1] ; 
495     cin >> answer ; cout << answer ; 
496     if ( answer == "y" ) 
497       AnalyzeOneEvent() ;
498   } 
499     if (fEvt != -999) {
500
501       Int_t Module ; 
502       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ; 
503       cin >> Module ; cout << Module << endl ; 
504       //=========== Creating 2d-histogram of the PHOS Module
505       // a little bit junkie but is used to test Geom functinalities
506       
507       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by Module   
508       
509       fGeom->EmcModuleCoverage(Module, tm, tM, pm, pM); 
510       // convert angles into coordinates local to the EMC module of interest
511       
512       Int_t EmcModuleNumber ;
513       Double_t EmcModulexm, EmcModulezm ; // minimum local coordinate in a given EMCA module
514       Double_t EmcModulexM, EmcModulezM ; // maximum local coordinate in a given EMCA module
515       fGeom->ImpactOnEmc(tm, pm, EmcModuleNumber, EmcModulezm, EmcModulexm) ;
516       fGeom->ImpactOnEmc(tM, pM, EmcModuleNumber, EmcModulezM, EmcModulexM) ;
517       Int_t xdim = (Int_t)( ( EmcModulexM - EmcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
518       Int_t zdim = (Int_t)( ( EmcModulezM - EmcModulezm ) / fGeom->GetCrystalSize(2) ) ;
519       Float_t xmin = EmcModulexm - fGeom->GetCrystalSize(0) ; 
520       Float_t xMax = EmcModulexM + fGeom->GetCrystalSize(0) ; 
521       Float_t zmin = EmcModulezm - fGeom->GetCrystalSize(2) ; 
522       Float_t zMax = EmcModulezM + fGeom->GetCrystalSize(2) ;     
523       Text_t HistoName[80];
524       sprintf(HistoName,"Event %d: Track Segments in module %d", fEvt, Module) ; 
525       TH2F * HistoTrack = new TH2F("HistoTrack",  HistoName, 
526                                    xdim, xmin, xMax, zdim, zmin, zMax) ;  
527       HistoTrack->SetStats(kFALSE); 
528       Text_t CanvasName[80];
529       sprintf(CanvasName,"Track segments in PHOS/EMC-PPSD module # %d", Module) ;
530       TCanvas * TrackCanvas = new TCanvas("TrackSegmentCanvas", CanvasName, 650, 500) ; 
531       HistoTrack->Draw() ; 
532
533       TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
534       AliPHOSTrackSegment * trseg ;
535  
536       Int_t NTrackSegments = trsegl->GetEntries() ;
537       Int_t index ;
538       Float_t Etot = 0 ;
539       Int_t NTrackSegmentsInModule = 0 ; 
540       for(index = 0; index < NTrackSegments ; index++){
541         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
542         Etot+= trseg->GetEnergy() ;
543         if ( trseg->GetPHOSMod() == Module ) { 
544           NTrackSegmentsInModule++ ; 
545           trseg->Draw("P");
546         }
547       } 
548       Text_t text[80] ; 
549       sprintf(text, "track segments: %d", NTrackSegmentsInModule) ;
550       TPaveText *  PaveText = new TPaveText(22, 80, 83, 90); 
551       PaveText->AddText(text) ; 
552       PaveText->Draw() ; 
553       TrackCanvas->Update() ; 
554       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< Etot << endl ;
555     
556    }
557 }
558 //____________________________________________________________________________
559 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
560 {
561   fRootFile   = new TFile(name) ;
562   return  fRootFile->IsOpen() ; 
563 }