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