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