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