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