New PID in AliPHOSPIDv1
[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()
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.05) ; 
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 kNEUTRAL_HA:
270                       fhNeutralHadronEnergy->Fill(recparticle->Energy() ) ; 
271                       //fhNeutralHadronPositionX->Fill(recpart. ) ;
272                       //fhNeutralHadronPositionY->Fill(recpart. ) ; 
273                       cout << "NEUTRAl HADRON" << endl;
274                       break ;
275                     case kNEUTRAL_EM:
276                       fhNeutralEMEnergy->Fill(recparticle->Energy() ) ; 
277                       //fhNeutralEMPositionX->Fill(recpart. ) ;
278                       //fhNeutralEMPositionY->Fill(recpart. ) ; 
279                       //cout << "NEUTRAL EM" << endl;
280                       break ;
281                     case kCHARGED_HA:
282                       fhChargedHadronEnergy->Fill(recparticle->Energy() ) ; 
283                       //fhChargedHadronPositionX->Fill(recpart. ) ;
284                       //fhChargedHadronPositionY->Fill(recpart. ) ; 
285                       cout << "CHARGED HADRON" << endl;
286                       break ;
287                     case kGAMMA_HA:
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(0.50) ; 
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     //fPID->SetShowerProfileCuts(Float_t l1m, Float_t l1M, Float_t l2m, Float_t l2M) ; 
418     fPID->SetShowerProfileCuts(0.7, 2.0 , 0.6 , 1.5) ; 
419
420     //========== Creates the Reconstructioner  
421     
422     fRec = new AliPHOSReconstructioner(fClu, fTrs, fPID) ;
423     fRec -> SetDebugReconstruction(kTRUE);     
424     
425     //=========== Connect the various Tree's for evt
426     
427     if ( evt == -999 ) {
428       cout <<  "AnalyzeOneEvent > Enter event number : " ; 
429       cin >> evt ;  
430       cout << evt << endl ; 
431     }
432     fEvt = evt ; 
433     gAlice->GetEvent(evt);
434     
435     //=========== Get the Digit TTree
436     
437     gAlice->TreeD()->GetEvent(0) ;     
438     
439   } // ok
440   
441   return ok ; 
442 }
443
444
445 //____________________________________________________________________________
446 void AliPHOSAnalyze::DisplayKineEvent(Int_t evt)
447 {
448   // Display particles from the Kine Tree in global Alice (theta, phi) coordinates. 
449   // One PHOS module at the time.
450   // The particle type can be selected.
451   
452   if (evt == -999) 
453     evt = fEvt ;
454
455   Int_t module ; 
456   cout <<  "DisplayKineEvent > which module (1-5,  -1: all) ? " ; 
457   cin >> module ; cout << module << endl ; 
458
459   Int_t testparticle ; 
460   cout << " 22      : PHOTON " << endl 
461        << " (-)11   : (POSITRON)ELECTRON " << endl 
462        << " (-)2112 : (ANTI)NEUTRON " << endl  
463        << " -999    : Everything else " << endl ; 
464   cout  <<  "DisplayKineEvent > enter PDG particle code to display " ; 
465   cin >> testparticle ; cout << testparticle << endl ; 
466
467   Text_t histoname[80] ;
468   sprintf(histoname,"Event %d: Incident particles in module %d", evt, module) ; 
469
470   Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
471   fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
472
473   Double_t theta, phi ; 
474   fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
475
476   Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
477   Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
478
479   tm -= theta ; 
480   tM += theta ; 
481   pm -= phi ; 
482   pM += phi ; 
483
484   TH2F * histoparticle = new TH2F("histoparticle",  histoname, 
485                                           pdim, pm, pM, tdim, tm, tM) ; 
486   histoparticle->SetStats(kFALSE) ;
487
488   // Get pointers to Alice Particle TClonesArray
489
490   TParticle * particle;
491   TClonesArray * particlearray  = gAlice->Particles();    
492
493   Text_t canvasname[80];
494   sprintf(canvasname,"Particles incident in PHOS/EMC module # %d",module) ;
495   TCanvas * kinecanvas = new TCanvas("kinecanvas", canvasname, 650, 500) ; 
496
497   // get the KINE Tree
498
499   TTree * kine =  gAlice->TreeK() ; 
500   Stat_t nParticles =  kine->GetEntries() ; 
501   cout << "DisplayKineEvent > events in kine " << nParticles << endl ; 
502   
503   // loop over particles
504
505   Double_t kRADDEG = 180. / TMath::Pi() ; 
506   Int_t index1 ; 
507   Int_t nparticlein = 0 ; 
508   for (index1 = 0 ; index1 < nParticles ; index1++){
509     Int_t nparticle = particlearray->GetEntriesFast() ;
510     Int_t index2 ; 
511     for( index2 = 0 ; index2 < nparticle ; index2++) {         
512       particle            = (TParticle*)particlearray->UncheckedAt(index2) ;
513       Int_t  particletype = particle->GetPdgCode() ;
514       if (testparticle == -999 || testparticle == particletype) { 
515         Double_t phi        = particle->Phi() ;
516         Double_t theta      = particle->Theta() ;
517         Int_t mod ; 
518         Double_t x, z ; 
519         fGeom->ImpactOnEmc(theta, phi, mod, z, x) ;
520         if ( mod == module ) {
521           nparticlein++ ; 
522           if (particle->Energy() >  fClu->GetEmcClusteringThreshold()  )
523             histoparticle->Fill(phi*kRADDEG, theta*kRADDEG, particle->Energy() ) ; 
524         } 
525       } 
526     }
527   }
528   kinecanvas->Draw() ; 
529   histoparticle->Draw("color") ; 
530   TPaveText *  pavetext = new TPaveText(294, 100, 300, 101); 
531   Text_t text[40] ; 
532   sprintf(text, "Particles: %d ", nparticlein) ;
533   pavetext->AddText(text) ; 
534   pavetext->Draw() ; 
535   kinecanvas->Update(); 
536
537 }
538 //____________________________________________________________________________
539 void AliPHOSAnalyze::DisplayRecParticles()
540 {
541   // Display reconstructed particles in global Alice(theta, phi) coordinates. 
542   // One PHOS module at the time.
543   // Click on symbols indicate the reconstructed particle type. 
544
545   if (fEvt == -999) {
546     cout << "DisplayRecParticles > Analyze an event first ... (y/n) " ; 
547     Text_t answer[1] ; 
548     cin >> answer ; cout << answer ; 
549     if ( answer == "y" ) 
550       AnalyzeOneEvent() ;
551   } 
552     if (fEvt != -999) {
553       
554       Int_t module ; 
555       cout <<  "DisplayRecParticles > which module (1-5,  -1: all) ? " ; 
556       cin >> module ; cout << module << endl ;
557       Text_t histoname[80] ; 
558       sprintf(histoname,"Event %d: Reconstructed particles in module %d", fEvt, module) ; 
559       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
560       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM, kDegre) ;
561       Double_t theta, phi ; 
562       fGeom->EmcXtalCoverage(theta, phi, kDegre) ;
563       Int_t tdim = (Int_t)( (tM - tm) / theta ) ; 
564       Int_t pdim = (Int_t)( (pM - pm) / phi ) ; 
565       tm -= theta ; 
566       tM += theta ; 
567       pm -= phi ; 
568       TH2F * histoRparticle = new TH2F("histoRparticle",  histoname, 
569                                        pdim, pm, pM, tdim, tm, tM) ; 
570       histoRparticle->SetStats(kFALSE) ;
571       Text_t canvasname[80] ; 
572       sprintf(canvasname, "Reconstructed particles in PHOSmodule # %d", module) ;
573       TCanvas * rparticlecanvas = new TCanvas("RparticleCanvas", canvasname, 650, 500) ; 
574       RecParticlesList * rpl = fPHOS->RecParticles() ; 
575       Int_t nRecParticles = rpl->GetEntries() ; 
576       Int_t nRecParticlesInModule = 0 ; 
577       TIter nextRecPart(rpl) ; 
578       AliPHOSRecParticle * rp ; 
579       cout << "DisplayRecParticles > " << nRecParticles << " reconstructed particles " << endl ; 
580       Double_t kRADDEG = 180. / TMath::Pi() ; 
581       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
582         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
583         if ( ts->GetPHOSMod() == module ) {
584           Int_t numberofprimaries = 0 ;
585           Int_t * listofprimaries = rp->GetPrimaries(numberofprimaries) ;
586           cout << "Number of primaries = " << numberofprimaries << endl ; 
587           Int_t index ;
588           for ( index = 0 ; index < numberofprimaries ; index++)
589             cout << "    primary # " << index << " =  " << listofprimaries[index] << endl ;  
590           
591           nRecParticlesInModule++ ; 
592           Double_t theta = rp->Theta() * kRADDEG ;
593           Double_t phi   = rp->Phi() * kRADDEG ;
594           Double_t energy = rp->Energy() ; 
595           histoRparticle->Fill(phi, theta, energy) ;
596         }
597       }
598       histoRparticle->Draw("color") ; 
599
600       nextRecPart.Reset() ; 
601       while ( (rp = (AliPHOSRecParticle *)nextRecPart() ) ) {
602         AliPHOSTrackSegment * ts = rp->GetPHOSTrackSegment() ; 
603         if ( ts->GetPHOSMod() == module )  
604           rp->Draw("P") ; 
605       }
606
607       Text_t text[80] ; 
608       sprintf(text, "reconstructed particles: %d", nRecParticlesInModule) ;
609       TPaveText *  pavetext = new TPaveText(292, 100, 300, 101); 
610       pavetext->AddText(text) ; 
611       pavetext->Draw() ; 
612       rparticlecanvas->Update() ; 
613     }
614 }
615
616 //____________________________________________________________________________
617 void AliPHOSAnalyze::DisplayRecPoints()
618 {
619   // Display reconstructed points in local PHOS-module (x, z) coordinates. 
620   // One PHOS module at the time.
621   // Click on symbols displays the EMC cluster, or PPSD information.
622
623   if (fEvt == -999) {
624     cout << "DisplayRecPoints > Analyze an event first ... (y/n) " ; 
625     Text_t answer[1] ; 
626     cin >> answer ; cout << answer ; 
627     if ( answer == "y" ) 
628       AnalyzeOneEvent() ;
629   } 
630     if (fEvt != -999) {
631       
632       Int_t module ; 
633       cout <<  "DisplayRecPoints > which module (1-5,  -1: all) ? " ; 
634       cin >> module ; cout << module << endl ; 
635
636       Text_t canvasname[80];
637       sprintf(canvasname,"Digits in PHOS/EMC module # %d",module) ;
638       TCanvas * modulecanvas = new TCanvas("module", canvasname, 650, 500) ; 
639       modulecanvas->Draw() ;
640
641       //=========== Creating 2d-histogram of the PHOS module
642       // a little bit junkie but is used to test Geom functinalities
643
644       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
645       
646       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
647       // convert angles into coordinates local to the EMC module of interest
648
649       Int_t emcModuleNumber ;
650       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
651       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
652       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
653       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
654       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
655       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
656       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
657       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
658       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
659       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
660       Text_t histoname[80];
661       sprintf(histoname,"Event %d: Digits and RecPoints in module %d", fEvt, module) ;
662       TH2F * hModule = new TH2F("HistoReconstructed", histoname,
663                                 xdim, xmin, xMax, zdim, zmin, zMax) ;  
664       hModule->SetMaximum(2.0);
665       hModule->SetMinimum(0.0);
666       hModule->SetStats(kFALSE); 
667
668       TIter next(fPHOS->Digits()) ;
669       Float_t energy, y, z;
670       Float_t etot=0.;
671       Int_t relid[4]; Int_t nDigits = 0 ;
672       AliPHOSDigit * digit ; 
673
674       // Making 2D histogram of the EMC module
675       while((digit = (AliPHOSDigit *)next())) 
676         {  
677           fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
678           if (relid[0] == module && relid[1] == 0)  
679             {  
680               energy = fClu->Calibrate(digit->GetAmp()) ;
681               cout << "Energy is " << energy << " and threshold is " << fClu->GetEmcEnergyThreshold() << endl; 
682               if (energy >  fClu->GetEmcEnergyThreshold()  ){
683                 nDigits++ ;
684                 etot += energy ; 
685                 fGeom->RelPosInModule(relid,y,z) ;   
686                 hModule->Fill(y, z, energy) ;
687               }
688             } 
689         }
690       cout <<"DrawRecPoints >  Found in module " 
691            << module << " " << nDigits << "  digits with total energy " << etot << endl ;
692       hModule->Draw("col2") ;
693
694       //=========== Cluster in module
695
696       //      TClonesArray * emcRP = fPHOS->EmcClusters() ; 
697       TObjArray * emcRP = fPHOS->EmcRecPoints() ; 
698       
699       etot = 0.; 
700       Int_t totalnClusters = 0 ; 
701       Int_t nClusters = 0 ;
702       TIter nextemc(emcRP) ;
703       AliPHOSEmcRecPoint * emc ;
704       while((emc = (AliPHOSEmcRecPoint *)nextemc())) 
705         {
706           //      Int_t numberofprimaries ;
707           //      Int_t * primariesarray = new Int_t[10] ;
708           //      emc->GetPrimaries(numberofprimaries, primariesarray) ;
709           totalnClusters++ ;
710           if ( emc->GetPHOSMod() == module )
711             { 
712               nClusters++ ; 
713               energy = emc->GetTotalEnergy() ;   
714               etot+= energy ;  
715               emc->Draw("M") ;
716             }
717         }
718       cout << "DrawRecPoints > Found " << totalnClusters << " EMC Clusters in PHOS" << endl ; 
719       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " EMC Clusters " << endl ;
720       cout << "DrawRecPoints > total energy  " << etot << endl ; 
721
722       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
723       Text_t text[40] ; 
724       sprintf(text, "digits: %d;  clusters: %d", nDigits, nClusters) ;
725       pavetext->AddText(text) ; 
726       pavetext->Draw() ; 
727       modulecanvas->Update(); 
728  
729       //=========== Cluster in module PPSD Down
730
731       //      TClonesArray * ppsdRP = fPHOS->PpsdClusters() ;
732       TObjArray * ppsdRP = fPHOS->PpsdRecPoints() ;
733  
734       etot = 0.; 
735       TIter nextPpsd(ppsdRP) ;
736       AliPHOSPpsdRecPoint * ppsd ;
737       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsd())) 
738         {
739           totalnClusters++ ;
740           if ( ppsd->GetPHOSMod() == module )
741             { 
742               nClusters++ ; 
743               energy = ppsd->GetEnergy() ;   
744               etot+=energy ;  
745               if (!ppsd->GetUp()) ppsd->Draw("P") ;
746             }
747         }
748       cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Down Clusters in PHOS" << endl ; 
749       cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Down Clusters " << endl ;
750       cout << "DrawRecPoints > total energy  " << etot << endl ; 
751
752       //=========== Cluster in module PPSD Up
753   
754       ppsdRP = fPHOS->PpsdRecPoints() ;
755      
756       etot = 0.; 
757       TIter nextPpsdUp(ppsdRP) ;
758       while((ppsd = (AliPHOSPpsdRecPoint *)nextPpsdUp())) 
759         {
760           totalnClusters++ ;
761           if ( ppsd->GetPHOSMod() == module )
762             { 
763               nClusters++ ; 
764               energy = ppsd->GetEnergy() ;   
765               etot+=energy ;  
766               if (ppsd->GetUp()) ppsd->Draw("P") ;
767             }
768         }
769   cout << "DrawRecPoints > Found " << totalnClusters << " Ppsd Up Clusters in PHOS" << endl ; 
770   cout << "DrawRecPoints > Found in module " << module << "  " << nClusters << " Ppsd Up Clusters " << endl ;
771   cout << "DrawRecPoints > total energy  " << etot << endl ; 
772     
773     } // if !-999
774 }
775
776 //____________________________________________________________________________
777 void AliPHOSAnalyze::DisplayTrackSegments()
778 {
779   // Display track segments in local PHOS-module (x, z) coordinates. 
780   // One PHOS module at the time.
781   // One symbol per PHOS subsystem: EMC, upper PPSD, lower PPSD.
782
783   if (fEvt == -999) {
784     cout << "DisplayTrackSegments > Analyze an event first ... (y/n) " ; 
785     Text_t answer[1] ; 
786     cin >> answer ; cout << answer ; 
787     if ( answer == "y" ) 
788       AnalyzeOneEvent() ;
789   } 
790     if (fEvt != -999) {
791
792       Int_t module ; 
793       cout <<  "DisplayTrackSegments > which module (1-5,  -1: all) ? " ; 
794       cin >> module ; cout << module << endl ; 
795       //=========== Creating 2d-histogram of the PHOS module
796       // a little bit junkie but is used to test Geom functinalities
797       
798       Double_t tm, tM, pm, pM ; // min and Max theta and phi covered by module   
799       
800       fGeom->EmcModuleCoverage(module, tm, tM, pm, pM); 
801       // convert angles into coordinates local to the EMC module of interest
802       
803       Int_t emcModuleNumber ;
804       Double_t emcModulexm, emcModulezm ; // minimum local coordinate in a given EMCA module
805       Double_t emcModulexM, emcModulezM ; // maximum local coordinate in a given EMCA module
806       fGeom->ImpactOnEmc(tm, pm, emcModuleNumber, emcModulezm, emcModulexm) ;
807       fGeom->ImpactOnEmc(tM, pM, emcModuleNumber, emcModulezM, emcModulexM) ;
808       Int_t xdim = (Int_t)( ( emcModulexM - emcModulexm ) / fGeom->GetCrystalSize(0) ) ;  
809       Int_t zdim = (Int_t)( ( emcModulezM - emcModulezm ) / fGeom->GetCrystalSize(2) ) ;
810       Float_t xmin = emcModulexm - fGeom->GetCrystalSize(0) ; 
811       Float_t xMax = emcModulexM + fGeom->GetCrystalSize(0) ; 
812       Float_t zmin = emcModulezm - fGeom->GetCrystalSize(2) ; 
813       Float_t zMax = emcModulezM + fGeom->GetCrystalSize(2) ;     
814       Text_t histoname[80];
815       sprintf(histoname,"Event %d: Track Segments in module %d", fEvt, module) ; 
816       TH2F * histotrack = new TH2F("histotrack",  histoname, 
817                                    xdim, xmin, xMax, zdim, zmin, zMax) ;  
818       histotrack->SetStats(kFALSE); 
819       Text_t canvasname[80];
820       sprintf(canvasname,"Track segments in PHOS/EMC-PPSD module # %d", module) ;
821       TCanvas * trackcanvas = new TCanvas("TrackSegmentCanvas", canvasname, 650, 500) ; 
822       histotrack->Draw() ; 
823
824       TrackSegmentsList * trsegl = fPHOS->TrackSegments() ;
825       AliPHOSTrackSegment * trseg ;
826  
827       Int_t nTrackSegments = trsegl->GetEntries() ;
828       Int_t index ;
829       Float_t etot = 0 ;
830       Int_t nTrackSegmentsInModule = 0 ; 
831       for(index = 0; index < nTrackSegments ; index++){
832         trseg = (AliPHOSTrackSegment * )trsegl->At(index) ;
833         etot+= trseg->GetEnergy() ;
834         if ( trseg->GetPHOSMod() == module ) { 
835           nTrackSegmentsInModule++ ; 
836           trseg->Draw("P");
837         }
838       } 
839       Text_t text[80] ; 
840       sprintf(text, "track segments: %d", nTrackSegmentsInModule) ;
841       TPaveText *  pavetext = new TPaveText(22, 80, 83, 90); 
842       pavetext->AddText(text) ; 
843       pavetext->Draw() ; 
844       trackcanvas->Update() ; 
845       cout << "DisplayTrackSegments > Found " << trsegl->GetEntries() << " Track segments with total energy "<< etot << endl ;
846     
847    }
848 }
849 //____________________________________________________________________________
850 Bool_t AliPHOSAnalyze::OpenRootFile(Text_t * name)
851 {
852   // Open the root file named "name"
853   
854   fRootFile   = new TFile(name, "update") ;
855   return  fRootFile->IsOpen() ; 
856 }
857 //____________________________________________________________________________
858 void AliPHOSAnalyze::SavingHistograms()
859 {
860   // Saves the histograms in a root file named "name.analyzed" 
861
862   Text_t outputname[80] ;
863   sprintf(outputname,"%s.analyzed",fRootFile->GetName());
864   TFile output(outputname,"RECREATE");
865   output.cd();
866   if (fhEmcDigit )         
867     fhEmcDigit->Write()  ;
868   if (fhVetoDigit )  
869     fhVetoDigit->Write()  ;
870   if (fhConvertorDigit ) 
871     fhConvertorDigit->Write()   ;
872   if (fhEmcCluster   )
873     fhEmcCluster->Write()   ;
874   if (fhVetoCluster ) 
875     fhVetoCluster->Write()   ;
876   if (fhConvertorCluster )
877     fhConvertorCluster->Write()  ;
878   if (fhConvertorEmc ) 
879     fhConvertorEmc->Write()  ;
880   if (fhPhotonEnergy)    
881     fhPhotonEnergy->Write() ;
882   if (fhPhotonPositionX)  
883     fhPhotonPositionX->Write() ;
884   if (fhPhotonPositionY)  
885     fhPhotonPositionX->Write() ;
886   if (fhElectronEnergy)  
887     fhElectronEnergy->Write() ;
888   if (fhElectronPositionX)
889     fhElectronPositionX->Write() ;
890   if (fhElectronPositionY) 
891     fhElectronPositionX->Write() ;
892   if (fhNeutralHadronEnergy) 
893     fhNeutralHadronEnergy->Write() ;
894   if (fhNeutralHadronPositionX)
895     fhNeutralHadronPositionX->Write() ;
896   if (fhNeutralHadronPositionY) 
897     fhNeutralHadronPositionX->Write() ;
898   if (fhNeutralEMEnergy)   
899     fhNeutralEMEnergy->Write() ;
900   if (fhNeutralEMPositionX)
901     fhNeutralEMPositionX->Write() ;
902   if (fhNeutralEMPositionY) 
903     fhNeutralEMPositionX->Write() ;
904   if (fhChargedHadronEnergy) 
905     fhChargedHadronEnergy->Write() ;
906   if (fhChargedHadronPositionX) 
907     fhChargedHadronPositionX->Write() ;
908   if (fhChargedHadronPositionY)
909     fhChargedHadronPositionX->Write() ;
910   if (fhPhotonHadronEnergy) 
911     fhPhotonHadronEnergy->Write() ;
912   if (fhPhotonHadronPositionX) 
913     fhPhotonHadronPositionX->Write() ;
914   if (fhPhotonHadronPositionY)
915     fhPhotonHadronPositionX->Write() ;
916
917   output.Write();
918   output.Close();
919 }