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