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