c220e3dd5d7fd2b8d82b4ea02e7e7e5a1709dd42
[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 /* $Id$ */
17
18 //_________________________________________________________________________
19 // Algorythm class to analyze PHOSv0 events:
20 // Construct histograms and displays them.
21 // Use the macro EditorBar.C for best access to the functionnalities
22 //
23 //*-- Author: Y. Schutz (SUBATECH) 
24 //////////////////////////////////////////////////////////////////////////////
25
26 // --- ROOT system ---
27
28 #include "TFile.h"
29 #include "TH1.h"
30 #include "TPad.h"
31 #include "TH2.h"
32 #include "TParticle.h"
33 #include "TClonesArray.h"
34 #include "TTree.h"
35 #include "TMath.h"
36 #include "TCanvas.h" 
37
38 // --- Standard library ---
39
40 #include <iostream.h>
41 #include <stdio.h>
42
43 // --- AliRoot header files ---
44
45 #include "AliRun.h"
46 #include "AliPHOSAnalyze.h"
47 #include "AliPHOSClusterizerv1.h"
48 #include "AliPHOSTrackSegmentMakerv1.h"
49 #include "AliPHOSPIDv1.h"
50 #include "AliPHOSReconstructioner.h"
51 #include "AliPHOSDigit.h"
52 #include "AliPHOSTrackSegment.h"
53 #include "AliPHOSRecParticle.h"
54
55 ClassImp(AliPHOSAnalyze)
56
57
58 //____________________________________________________________________________
59   AliPHOSAnalyze::AliPHOSAnalyze()
60 {
61   // default ctor (useless)
62   
63   fRootFile = 0 ; 
64 }
65
66 //____________________________________________________________________________
67 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * name)
68 {
69   // ctor: analyze events from root file: name
70   
71   Bool_t ok = OpenRootFile(name)  ; 
72   if ( !ok ) {
73     cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
74   }
75   else { 
76     gAlice = (AliRun*) fRootFile->Get("gAlice");
77     fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
78     fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() ) ;
79     fEvt = -999 ; 
80   }
81 }
82
83 //____________________________________________________________________________
84 AliPHOSAnalyze::~AliPHOSAnalyze()
85 {
86   // dtor
87
88   fRootFile->Close() ; 
89   delete fRootFile ; 
90   fRootFile = 0 ; 
91
92   delete fPHOS ; 
93   fPHOS = 0 ; 
94
95   delete fClu ; 
96   fClu = 0 ;
97
98   delete fPID ; 
99   fPID = 0 ;
100
101   delete fRec ; 
102   fRec = 0 ;
103
104   delete fTrs ; 
105   fTrs = 0 ;
106
107 }
108
109 //____________________________________________________________________________
110 void AliPHOSAnalyze::AnalyzeOneEvent(Int_t evt)
111 {
112   // analyze one single event with id=evt
113
114   TStopwatch ts ; 
115   ts.Start() ; 
116   Bool_t ok = Init(evt) ; 
117   
118   if ( ok ) {
119     //=========== Get the number of entries in the Digits array
120     
121     Int_t nId = fPHOS->Digits()->GetEntries();    
122     printf("AnalyzeOneEvent > Number of entries in the Digit array is %d \n",nId);
123     
124     //=========== Do the reconstruction
125     
126     cout << "AnalyzeOneEvent > Found  " << nId << "  digits in PHOS"   << endl ;  
127     
128    
129     fPHOS->Reconstruction(fRec);  
130     
131     // =========== End of reconstruction
132
133     // =========== Write the root file
134
135     fRootFile->Write() ;
136     
137     // =========== Finish
138
139     cout << "AnalyzeOneEvent > event # " << fEvt << " processed" << endl ;   
140   } // ok
141   else
142     cout << "AnalyzeOneEvent > filed to process event # " << evt << endl ;   
143   
144   ts.Stop() ;   cout << "CPU  time = " << ts.CpuTime()   << endl ; 
145   cout << "Real time = " << ts.RealTime()  << endl ; 
146 }
147
148 //____________________________________________________________________________
149  void AliPHOSAnalyze::AnalyzeManyEvents(Int_t Nevents, Int_t module)    
150 {
151   // analyzes Nevents events in a single PHOS module  
152
153   if ( fRootFile == 0 ) 
154     cout << "AnalyzeManyEvents > " << "Root File not openned" << endl ;  
155   else
156     {
157       //========== Get AliRun object from file 
158       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
159       //=========== Get the PHOS object and associated geometry from the file      
160       fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
161       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
162       //========== Booking Histograms
163       cout << "AnalyzeManyEvents > " << "Booking Histograms" << endl ; 
164       BookingHistograms();
165       Int_t ievent;
166       Int_t relid[4] ; 
167       AliPHOSDigit * digit ;
168       AliPHOSEmcRecPoint * emc ;
169       AliPHOSPpsdRecPoint * ppsd ;
170       //      AliPHOSTrackSegment * tracksegment ;
171       AliPHOSRecParticle * recparticle;
172       for ( ievent=0; ievent<Nevents; ievent++)
173         {  
174           if (ievent==0)  cout << "AnalyzeManyEvents > " << "Starting Analyzing " << endl ; 
175           //========== Create the Clusterizer
176           fClu = new AliPHOSClusterizerv1() ; 
177           fClu->SetEmcEnergyThreshold(0.025) ; 
178           fClu->SetEmcClusteringThreshold(0.50) ; 
179           fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
180           fClu->SetPpsdClusteringThreshold(0.0000001) ; 
181           fClu->SetLocalMaxCut(0.03) ;
182           fClu->SetCalibrationParameters(0., 0.00000001) ;  
183           //========== Creates the track segment maker
184           fTrs = new AliPHOSTrackSegmentMakerv1()  ;
185           fTrs->UnsetUnfoldFlag() ; 
186           //========== Creates the particle identifier
187           fPID = new AliPHOSPIDv1() ;
188           fPID->SetShowerProfileCuts(0.3, 1.8, 0.3, 1.8 ) ; 
189           fPID->Print() ;           
190           //========== Creates the Reconstructioner  
191           fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ; 
192           //========== Event Number>         
193           if ( ( log10((Float_t)(ievent+1)) - (Int_t)(log10((Float_t)(ievent+1))) ) == 0. ) 
194             cout <<  "AnalyzeManyEvents > " << "Event is " << ievent << endl ;  
195           //=========== Connects the various Tree's for evt
196           gAlice->GetEvent(ievent);
197           //=========== Gets the Digit TTree
198           gAlice->TreeD()->GetEvent(0) ;
199           //=========== Gets the number of entries in the Digits array 
200           TIter nextdigit(fPHOS->Digits()) ;
201           while( ( digit = (AliPHOSDigit *)nextdigit() ) )
202             {
203               fGeom->AbsToRelNumbering(digit->GetId(), relid) ;         
204               if (fClu->IsInEmc(digit)) fhEmcDigit->Fill(fClu->Calibrate(digit->GetAmp())) ; 
205               else    
206                 {  
207                   if (relid[1]<17) fhVetoDigit->Fill(fClu->Calibrate(digit->GetAmp())); 
208                   if (relid[1]>16) fhConvertorDigit->Fill(fClu->Calibrate(digit->GetAmp()));
209                 }
210             }
211           //=========== Do the reconstruction
212           fPHOS->Reconstruction(fRec);
213           //=========== Cluster in module
214           TIter nextEmc(fPHOS->EmcClusters()  ) ;
215           while((emc = (AliPHOSEmcRecPoint *)nextEmc())) 
216             {
217               if ( emc->GetPHOSMod() == module )
218                 {  
219                   fhEmcCluster->Fill(  emc->GetTotalEnergy()  ); 
220                   TIter nextPpsd( fPHOS->PpsdClusters()) ;
221                   while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
222                     {
223                       if ( ppsd->GetPHOSMod() == module )
224                         {                         
225                           if (!ppsd->GetUp()) fhConvertorEmc->Fill(emc->GetTotalEnergy(),ppsd->GetTotalEnergy()) ;
226                         }
227                     } 
228                 }
229             }
230           //=========== Cluster in module PPSD Down
231           TIter nextPpsd(fPHOS->PpsdClusters() ) ;
232           while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
233             {
234               if ( ppsd->GetPHOSMod() == module )
235                 { 
236                   if (!ppsd->GetUp()) fhConvertorCluster->Fill(ppsd->GetTotalEnergy()) ;
237                   if (ppsd->GetUp())  fhVetoCluster     ->Fill(ppsd->GetTotalEnergy()) ;
238                 }
239             }
240           //========== TRackSegments in the event
241           TIter nextRecParticle(fPHOS->RecParticles() ) ; 
242           while((recparticle = (AliPHOSRecParticle *)nextRecParticle())) 
243             {
244               if ( recparticle->GetPHOSTrackSegment()->GetPHOSMod() == module )
245                 { 
246                   cout << "Particle type is " << recparticle->GetType() << endl ; 
247                   Int_t numberofprimaries = 0 ;
248                   Int_t * listofprimaries = recparticle->GetPrimaries(numberofprimaries) ;
249                   cout << "Number of primaries = " << numberofprimaries << endl ; 
250                   Int_t index ;
251                   for ( index = 0 ; index < numberofprimaries ; index++)
252                     cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
253                   switch(recparticle->GetType())
254                     {
255                     case kGAMMA:
256                       fhPhotonEnergy->Fill(recparticle->Energy() ) ; 
257                       //fhPhotonPositionX->Fill(recpart. ) ;
258                       //fhPhotonPositionY->Fill(recpart. ) ;                 
259                       cout << "PHOTON" << endl;
260                       break;
261                     case kELECTRON:
262                       fhElectronEnergy->Fill(recparticle->Energy() ) ; 
263                       //fhElectronPositionX->Fill(recpart. ) ;
264                       //fhElectronPositionY->Fill(recpart. ) ; 
265                       cout << "ELECTRON" << endl;
266                       break;
267                     case kNEUTRALHADRON:
268                       fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ; 
269                       //fhNeutralHadronPositionX->Fill(recpart. ) ;
270                       //fhNeutralHadronPositionY->Fill(recpart. ) ; 
271                       cout << "NEUTRAl HADRON" << endl;
272                       break ;
273                     case kNEUTRALEM:
274                       fhNeutralEMEnergy->Fill(recparticle->Energy() ) ; 
275                       //fhNeutralEMPositionX->Fill(recpart. ) ;
276                       //fhNeutralEMPositionY->Fill(recpart. ) ; 
277                       //cout << "NEUTRAL EM" << endl;
278                       break ;
279                     case kCHARGEDHADRON:
280                       fhChargedHadronEnergy->Fill(recparticle->Energy() ) ; 
281                       //fhChargedHadronPositionX->Fill(recpart. ) ;
282                       //fhChargedHadronPositionY->Fill(recpart. ) ; 
283                       cout << "CHARGED HADRON" << endl;
284                       break ;
285                     case kGAMMAHADRON:
286                       fhPhotonHadronEnergy->Fill(recparticle->Energy() ) ; 
287                       //fhPhotonHadronPositionX->Fill(recpart. ) ;
288                       //fhPhotonHadronPositionY->Fill(recpart. ) ; 
289                       cout << "PHOTON HADRON" << endl;
290                       break ;                 
291                     }
292                 }
293             }
294           // Deleting fClu, fTrs, fPID et fRec
295           fClu->Delete();
296           fTrs->Delete();
297           fPID->Delete();
298           fRec->Delete();
299
300         }   // endfor
301       SavingHistograms();
302     }       // endif
303 }           // endfunction
304
305
306 //____________________________________________________________________________
307 void  AliPHOSAnalyze::BookingHistograms()
308 {
309   // Books the histograms where the results of the analysis are stored (to be changed)
310
311   if (fhEmcDigit )         
312     delete fhEmcDigit  ;
313   if (fhVetoDigit )     
314     delete fhVetoDigit  ;
315   if (fhConvertorDigit ) 
316     delete fhConvertorDigit   ;
317   if (fhEmcCluster   )  
318     delete  fhEmcCluster   ;
319   if (fhVetoCluster )     
320     delete fhVetoCluster   ;
321   if (fhConvertorCluster )
322     delete fhConvertorCluster  ;
323   if (fhConvertorEmc )    
324     delete fhConvertorEmc  ;
325  
326   fhEmcDigit                = new TH1F("hEmcDigit",      "hEmcDigit",         1000,  0. ,  25.);
327   fhVetoDigit               = new TH1F("hVetoDigit",     "hVetoDigit",         500,  0. ,  3.e-5);
328   fhConvertorDigit          = new TH1F("hConvertorDigit","hConvertorDigit",    500,  0. ,  3.e-5);
329   fhEmcCluster              = new TH1F("hEmcCluster",    "hEmcCluster",       1000,  0. ,  30.);
330   fhVetoCluster             = new TH1F("hVetoCluster",   "hVetoCluster",       500,  0. ,  3.e-5);
331   fhConvertorCluster        = new TH1F("hConvertorCluster","hConvertorCluster",500,  0. ,  3.e-5);
332   fhConvertorEmc            = new TH2F("hConvertorEmc",  "hConvertorEmc",      200,  1. ,  3., 200, 0., 3.e-5);
333   fhPhotonEnergy            = new TH1F("hPhotonEnergy",  "hPhotonEnergy",     1000,  0. ,  30.);
334   fhElectronEnergy          = new TH1F("hElectronEnergy","hElectronEnergy",   1000,  0. ,  30.);
335   fhNeutralHadronEnergy     = new TH1F("hNeutralHadronEnergy", "hNeutralHadronEnergy",    1000,  0. ,  30.);
336   fhNeutralEMEnergy         = new TH1F("hNeutralEMEnergy", "hNeutralEMEnergy",    1000,  0. ,  30.);
337   fhChargedHadronEnergy     = new TH1F("hChargedHadronEnergy", "hChargedHadronEnergy",    1000,  0. ,  30.);
338   fhPhotonHadronEnergy      = new TH1F("hPhotonHadronEnergy","hPhotonHadronEnergy",500,-80. , 80.);
339   fhPhotonPositionX         = new TH1F("hPhotonPositionX","hPhotonPositionX",   500,-80. , 80.);
340   fhElectronPositionX       = new TH1F("hElectronPositionX","hElectronPositionX",500,-80. , 80.);
341   fhNeutralHadronPositionX  = new TH1F("hNeutralHadronPositionX","hNeutralHadronPositionX",500,-80. , 80.);
342   fhNeutralEMPositionX      = new TH1F("hNeutralEMPositionX","hNeutralEMPositionX",500,-80. , 80.);
343   fhChargedHadronPositionX  = new TH1F("hChargedHadronPositionX","hChargedHadronPositionX",500,-80. , 80.);
344   fhPhotonHadronPositionX   = new TH1F("hPhotonHadronPositionX","hPhotonHadronPositionX",500,-80. , 80.);
345   fhPhotonPositionY         = new TH1F("hPhotonPositionY","hPhotonPositionY",   500,-80. , 80.);
346   fhElectronPositionY       = new TH1F("hElectronPositionY","hElectronPositionY",500,-80. , 80.);
347   fhNeutralHadronPositionY  = new TH1F("hNeutralHadronPositionY","hNeutralHadronPositionY",500,-80. , 80.);
348   fhNeutralEMPositionY      = new TH1F("hNeutralEMPositionY","hNeutralEMPositionY",500,-80. , 80.);
349   fhChargedHadronPositionY  = new TH1F("hChargedHadronPositionY","hChargedHadronPositionY",500,-80. , 80.);
350   fhPhotonHadronPositionY   = new TH1F("hPhotonHadronPositionY","hPhotonHadronPositionY",500,-80. , 80.);
351
352
353 }
354 //____________________________________________________________________________
355 Bool_t AliPHOSAnalyze::Init(Int_t evt)
356 {
357   // Do a few initializations: open the root file
358   //                           get the AliRun object
359   //                           defines the clusterizer, tracksegment maker and particle identifier
360   //                           sets the associated parameters
361
362   Bool_t ok = kTRUE ; 
363   
364    //========== Open galice root file  
365
366   if ( fRootFile == 0 ) {
367     Text_t * name  = new Text_t[80] ; 
368     cout << "AnalyzeOneEvent > Enter file root file name : " ;  
369     cin >> name ; 
370     Bool_t ok = OpenRootFile(name) ; 
371     if ( !ok )
372       cout << " AliPHOSAnalyze > Error opening " << name << endl ; 
373     else { 
374       //========== Get AliRun object from file 
375       
376       gAlice = (AliRun*) fRootFile->Get("gAlice") ;
377       
378       //=========== Get the PHOS object and associated geometry from the file 
379       
380       fPHOS  = (AliPHOSv0 *)gAlice->GetDetector("PHOS") ;
381       fGeom  = AliPHOSGeometry::GetInstance( fPHOS->GetGeometry()->GetName(), fPHOS->GetGeometry()->GetTitle() );
382     } // else !ok
383   } // if fRootFile
384   
385   if ( ok ) {
386     
387     //========== Create the Clusterizer
388
389     fClu =  new AliPHOSClusterizerv1() ; 
390     fClu->SetEmcEnergyThreshold(0.030) ; 
391     fClu->SetEmcClusteringThreshold(1.0) ; 
392     fClu->SetPpsdEnergyThreshold    (0.0000002) ; 
393     fClu->SetPpsdClusteringThreshold(0.0000001) ; 
394     fClu->SetLocalMaxCut(0.03) ;
395     fClu->SetCalibrationParameters(0., 0.00000001) ;  
396     cout <<  "AnalyzeOneEvent > using clusterizer " << fClu->GetName() << endl ; 
397     fClu->PrintParameters() ; 
398     
399     //========== Creates the track segment maker
400     
401     fTrs = new AliPHOSTrackSegmentMakerv1() ;
402     cout <<  "AnalyzeOneEvent > using tack segment maker " << fTrs->GetName() << endl ; 
403     fTrs->UnsetUnfoldFlag() ;
404     
405     //========== Creates the particle identifier
406     
407     fPID = new AliPHOSPIDv1() ;
408     cout <<  "AnalyzeOneEvent > using particle identifier " << fPID->GetName() << endl ; 
409     
410     //========== Creates the Reconstructioner  
411     
412     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;     
413     
414     //=========== Connect the various Tree's for evt
415     
416     if ( evt == -999 ) {
417       cout <<  "AnalyzeOneEvent > Enter event number : " ; 
418       cin >> evt ;  
419       cout << evt << endl ; 
420     }
421     fEvt = evt ; 
422     gAlice->GetEvent(evt);
423     
424     //=========== Get the Digit TTree
425     
426     gAlice->TreeD()->GetEvent(0) ;     
427     
428   } // ok
429   
430   return ok ; 
431 }
432
433
434 //____________________________________________________________________________
435 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
436 {
437   // Display particles from the Kine Tree in global Alice (theta, phi) coordinates. 
438   // One PHOS module at the time.
439   // The particle type can be selected.
440   
441   if (evt == -999) 
442     evt = fEvt ;
443
444   Int_t module ; 
445   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ; 
446   cin >> module ; cout << module << endl ; 
447
448   Int_t testparticle ; 
449   cout << " 22      : PHOTON " << endl 
450        << " (-)11   : (POSITRON)ELECTRON " << endl 
451        << " (-)2112 : (ANTI)NEUTRON " << endl  
452        << " -999    : Everything else " << endl ; 
453   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ; 
454   cin >> testparticle ; cout << testparticle << endl ; 
455
456   Text_t histoname[80] ;
457   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ; 
458
459   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
460   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
461
462   Double_t theta, phi ; 
463   fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
464
465   Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
466   Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
467
468   tm -= theta ; 
469   tM += theta ; 
470   pm -= phi ; 
471   pM += phi ; 
472
473   TH2F * histoparticle = new TH2F("histoparticle",  histoname, 
474                                           pdim, pm, pM, tdim, tm, tM) ; 
475   histoparticle->SetStats(kFALSE) ;
476
477   // Get pointers to Alice Particle TClonesArray
478
479   TParticle * particle;
480   TClonesArray * particlearray  = gAlice->Particles();    
481
482   Text_t canvasname[80];
483   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
484   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ; 
485
486   // get the KINE Tree
487
488   TTree * kine =  gAlice->TreeK() ; 
489   Stat_t nParticles =  kine->GetEntries() ; 
490   cout << "DisplayKineEvent > events in kine " << nParticles << endl ; 
491   
492   // loop over particles
493
494   Double_t kRADDEG = 180. / TMath::Pi() ; 
495   Int_t index1 ; 
496   Int_t nparticlein = 0 ; 
497   for (index1 = 0 ; index1 < nParticles ; index1++){
498     Int_t nparticle = particlearray->GetEntriesFast() ;
499     Int_t index2 ; 
500     for( index2 = 0 ; index2 < nparticle ; index2++) {         
501       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
502       Int_t  particletype = particle->GetPdgCode() ;
503       if (testparticle == -999 || testparticle == particletype) { 
504         Double_t phi        = particle->Phi() ;
505         Double_t theta      = particle->Theta() ;
506         Int_t mod ; 
507         Double_t x, z ; 
508         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
509         if ( mod == module ) {
510           nparticlein++ ; 
511           if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
512             histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ; 
513         } 
514       } 
515     }
516   }
517   kinecanvas->Draw() ; 
518   histoparticle->Draw("color") ; 
519   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101); 
520   Text_t text[40] ; 
521   sprintf(text, "Particles: %d ", nparticlein) ;
522   pavetext->AddText(text) ; 
523   pavetext->Draw() ; 
524   kinecanvas->Update(); 
525
526 }
527 //____________________________________________________________________________
528 void AliPHOSAnalyze::DisplayRecParticles()
529 {
530   // Display reconstructed particles in global Alice(theta, phi) coordinates. 
531   // One PHOS module at the time.
532   // Click on symbols indicate the reconstructed particle type. 
533
534   if (fEvt == -999) {
535     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; 
536     Text_t answer[1] ; 
537     cin >> answer ; cout << answer ; 
538     if ( answer == "y" ) 
539       AnalyzeOneEvent() ;
540   } 
541     if (fEvt != -999) {
542       
543       Int_t module ; 
544       cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ; 
545       cin >> module ; cout << module << endl ;
546       Text_t histoname[80] ; 
547       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ; 
548       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
549       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
550       Double_t theta, phi ; 
551       fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
552       Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
553       Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
554       tm -= theta ; 
555       tM += theta ; 
556       pm -= phi ; 
557       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname, 
558                                        pdim, pm, pM, tdim, tm, tM) ; 
559       histoRparticle->SetStats(kFALSE) ;
560       Text_t canvasname[80] ; 
561       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
562       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; 
563       RecParticlesList * rpl = fPHOS->RecParticles() ; 
564       Int_t nRecParticles = rpl->GetEntries() ; 
565       Int_t nRecParticlesInModule = 0 ; 
566       TIter nextRecPart(rpl) ; 
567       AliPHOSRecParticle * rp ; 
568       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ; 
569       Double_t kRADDEG = 180. / TMath::Pi() ; 
570       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
571         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
572         if ( ts->GetPHOSMod() == module ) {
573           Int_t numberofprimaries = 0 ;
574           Int_t * listofprimaries = rp->GetPrimaries(numberofprimaries) ;
575           cout << "Number of primaries = " << numberofprimaries << endl ; 
576           Int_t index ;
577           for ( index = 0 ; index < numberofprimaries ; index++)
578             cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
579           
580           nRecParticlesInModule++ ; 
581           Double_t theta = rp->Theta() * kRADDEG ;
582           Double_t phi   = rp->Phi() * kRADDEG ;
583           Double_t energy = rp->Energy() ; 
584           histoRparticle->Fill(phi, theta, energy) ;
585         }
586       }
587       histoRparticle->Draw("color") ; 
588
589       nextRecPart.Reset() ; 
590       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
591         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
592         if ( ts->GetPHOSMod() == module )  
593           rp->Draw("P") ; 
594       }
595
596       Text_t text[80] ; 
597       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
598       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101); 
599       pavetext->AddText(text) ; 
600       pavetext->Draw() ; 
601       rparticlecanvas->Update() ; 
602     }
603 }
604
605 //____________________________________________________________________________
606 void AliPHOSAnalyze::DisplayRecPoints()
607 {
608   // Display reconstructed points in local PHOS-module (x, z) coordinates. 
609   // One PHOS module at the time.
610   // Click on symbols displays the EMC cluster, or PPSD information.
611
612   if (fEvt == -999) {
613     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
614     Text_t answer[1] ; 
615     cin >> answer ; cout << answer ; 
616     if ( answer == "y" ) 
617       AnalyzeOneEvent() ;
618   } 
619     if (fEvt != -999) {
620       
621       Int_t module ; 
622       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ; 
623       cin >> module ; cout << module << endl ; 
624
625       Text_t canvasname[80];
626       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
627       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ; 
628       modulecanvas->Draw() ;
629
630       //=========== Creating 2d-histogram of the PHOS module
631       // a little bit junkie but is used to test Geom functinalities
632
633       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
634       
635       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
636       // convert angles into coordinates local to the EMC module of interest
637
638       Int_t emcModuleNumber ;
639       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
640       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
641       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
642       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
643       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
644       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
645       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
646       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
647       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
648       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
649       Text_t histoname[80];
650       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
651       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
652                                 xdim, xmin, xMax, zdim, zmin, zMax) ;  
653       hModule->SetMaximum(2.0);
654       hModule->SetMinimum(0.0);
655       hModule->SetStats(kFALSE); 
656
657       TIter next(fPHOS->Digits()) ;
658       Float_t energy, y, z;
659       Float_t etot=0.;
660       Int_t relid[4]; Int_t nDigits = 0 ;
661       AliPHOSDigit * digit ; 
662
663       // Making 2D histogram of the EMC module
664       while((digit = (AliPHOSDigit *)next())) 
665         {  
666           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
667           if (relid[0] == module && relid[1] == 0)  
668             {  
669               energy = fClu->Calibrate(digit->GetAmp()) ;
670               cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl; 
671               if (energy >  fClu->GetEmcEnergyThreshold()  ){
672                 nDigits++ ;
673                 etot += energy ; 
674                 fGeom->RelPosInModule(relid,y,z) ;   
675                 hModule->Fill(y, z, energy) ;
676               }
677             } 
678         }
679       cout <<"DrawRecPoints >  Found in module " 
680            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
681       hModule->Draw("col2") ;
682
683       //=========== Cluster in module
684
685       TClonesArray * emcRP = fPHOS->EmcClusters() ;
686       etot = 0.; 
687       Int_t totalnClusters = 0 ; 
688       Int_t nClusters = 0 ;
689       TIter nextemc(emcRP) ;
690       AliPHOSEmcRecPoint * emc ;
691       while((emc = (AliPHOSEmcRecPoint *)nextemc())) 
692         {
693           //      Int_t numberofprimaries ;
694           //      Int_t * primariesarray = new Int_t[10] ;
695           //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
696           totalnClusters++ ;
697           if ( emc->GetPHOSMod() == module )
698             { 
699               nClusters++ ; 
700               energy = emc->GetTotalEnergy() ;   
701               etot+= energy ;  
702               emc->Draw("M") ;
703             }
704         }
705       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; 
706       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
707       cout << "DrawRecPoints > total energy  " << etot << endl ; 
708
709       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
710       Text_t text[40] ; 
711       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
712       pavetext->AddText(text) ; 
713       pavetext->Draw() ; 
714       modulecanvas->Update(); 
715  
716       //=========== Cluster in module PPSD Down
717
718       TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
719       etot = 0.; 
720       TIter nextPpsd(ppsdRP) ;
721       AliPHOSPpsdRecPoint * ppsd ;
722       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
723         {
724           totalnClusters++ ;
725           if ( ppsd->GetPHOSMod() == module )
726             { 
727               nClusters++ ; 
728               energy = ppsd->GetEnergy() ;   
729               etot+=energy ;  
730               if (!ppsd->GetUp()) ppsd->Draw("P") ;
731             }
732         }
733       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ; 
734       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
735       cout << "DrawRecPoints > total energy  " << etot << endl ; 
736
737       //=========== Cluster in module PPSD Up
738   
739       ppsdRP = fPHOS->PpsdClusters() ;
740       etot = 0.; 
741       TIter nextPpsdUp(ppsdRP) ;
742       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) 
743         {
744           totalnClusters++ ;
745           if ( ppsd->GetPHOSMod() == module )
746             { 
747               nClusters++ ; 
748               energy = ppsd->GetEnergy() ;   
749               etot+=energy ;  
750               if (ppsd->GetUp()) ppsd->Draw("P") ;
751             }
752         }
753   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ; 
754   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
755   cout << "DrawRecPoints > total energy  " << etot << endl ; 
756     
757     } // if !-999
758 }
759
760 //____________________________________________________________________________
761 void AliPHOSAnalyze::DisplayTrackSegments()
762 {
763   // Display track segments in local PHOS-module (x, z) coordinates. 
764   // One PHOS module at the time.
765   // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
766
767   if (fEvt == -999) {
768     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; 
769     Text_t answer[1] ; 
770     cin >> answer ; cout << answer ; 
771     if ( answer == "y" ) 
772       AnalyzeOneEvent() ;
773   } 
774     if (fEvt != -999) {
775
776       Int_t module ; 
777       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ; 
778       cin >> module ; cout << module << endl ; 
779       //=========== Creating 2d-histogram of the PHOS module
780       // a little bit junkie but is used to test Geom functinalities
781       
782       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
783       
784       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
785       // convert angles into coordinates local to the EMC module of interest
786       
787       Int_t emcModuleNumber ;
788       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
789       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
790       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
791       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
792       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
793       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
794       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
795       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
796       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
797       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
798       Text_t histoname[80];
799       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ; 
800       TH2F * histotrack = new TH2F("histotrack",  histoname, 
801                                    xdim, xmin, xMax, zdim, zmin, zMax) ;  
802       histotrack->SetStats(kFALSE); 
803       Text_t canvasname[80];
804       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
805       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; 
806       histotrack->Draw() ; 
807
808       TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
809       AliPHOSTrackSegment * trseg ;
810  
811       Int_t nTrackSegments = trsegl->GetEntries() ;
812       Int_t index ;
813       Float_t etot = 0 ;
814       Int_t nTrackSegmentsInModule = 0 ; 
815       for(index = 0; index < nTrackSegments ; index++){
816         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
817         etot+= trseg->GetEnergy() ;
818         if ( trseg->GetPHOSMod() == module ) { 
819           nTrackSegmentsInModule++ ; 
820           trseg->Draw("P");
821         }
822       } 
823       Text_t text[80] ; 
824       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
825       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
826       pavetext->AddText(text) ; 
827       pavetext->Draw() ; 
828       trackcanvas->Update() ; 
829       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
830     
831    }
832 }
833 //____________________________________________________________________________
834 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
835 {
836   // Open the root file named "name"
837   
838   fRootFile   = new TFile(name, "update") ;
839   return  fRootFile->IsOpen() ; 
840 }
841 //____________________________________________________________________________
842 void AliPHOSAnalyze::SavingHistograms()
843 {
844   // Saves the histograms in a root file named "name.analyzed" 
845
846   Text_t outputname[80] ;
847   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
848   TFile output(outputname,"RECREATE");
849   output.cd();
850   if (fhEmcDigit )         
851     fhEmcDigit->Write()  ;
852   if (fhVetoDigit )  
853     fhVetoDigit->Write()  ;
854   if (fhConvertorDigit ) 
855     fhConvertorDigit->Write()   ;
856   if (fhEmcCluster   )
857     fhEmcCluster->Write()   ;
858   if (fhVetoCluster ) 
859     fhVetoCluster->Write()   ;
860   if (fhConvertorCluster )
861     fhConvertorCluster->Write()  ;
862   if (fhConvertorEmc ) 
863     fhConvertorEmc->Write()  ;
864   if (fhPhotonEnergy)    
865     fhPhotonEnergy->Write() ;
866   if (fhPhotonPositionX)  
867     fhPhotonPositionX->Write() ;
868   if (fhPhotonPositionY)  
869     fhPhotonPositionX->Write() ;
870   if (fhElectronEnergy)  
871     fhElectronEnergy->Write() ;
872   if (fhElectronPositionX)
873     fhElectronPositionX->Write() ;
874   if (fhElectronPositionY) 
875     fhElectronPositionX->Write() ;
876   if (fhNeutralHadronEnergy) 
877     fhNeutralHadronEnergy->Write() ;
878   if (fhNeutralHadronPositionX)
879     fhNeutralHadronPositionX->Write() ;
880   if (fhNeutralHadronPositionY) 
881     fhNeutralHadronPositionX->Write() ;
882   if (fhNeutralEMEnergy)   
883     fhNeutralEMEnergy->Write() ;
884   if (fhNeutralEMPositionX)
885     fhNeutralEMPositionX->Write() ;
886   if (fhNeutralEMPositionY) 
887     fhNeutralEMPositionX->Write() ;
888   if (fhChargedHadronEnergy) 
889     fhChargedHadronEnergy->Write() ;
890   if (fhChargedHadronPositionX) 
891     fhChargedHadronPositionX->Write() ;
892   if (fhChargedHadronPositionY)
893     fhChargedHadronPositionX->Write() ;
894   if (fhPhotonHadronEnergy) 
895     fhPhotonHadronEnergy->Write() ;
896   if (fhPhotonHadronPositionX) 
897     fhPhotonHadronPositionX->Write() ;
898   if (fhPhotonHadronPositionY)
899     fhPhotonHadronPositionX->Write() ;
900
901   output.Write();
902   output.Close();
903 }