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