Inheritance from TObject. Automatic streamers.
[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 PHOS events. In this class we demostrate, 
20 // how to handle reconstructed objects with AliPHSOIndexToObject.
21 // As an example we propose sulotions for four most frequently used tasks:
22 //    DrawRecon(...) - to draw reconstructed objects in the PHOS plane,
23 //                     very usefull in the debuging
24 //    InvarianMass(...) - to calculate "REAL" and "MIXED" photon pairs 
25 //                        invariant mass distributions
26 //    EnergyResoluition(...) -\ Energy and position resolutions of the 
27 //    PositionResolution(...)-/ reconstructed photons
28 //    Contamination(...) - calculates contamination of the photon spectrum and 
29 //                         pobability of reconstruction of several primaries as 
30 //                         kGAMMA,kELECTRON etc.
31 //
32 //    User Case:
33 //    root [0] AliPHOSAnalyze * a = new AliPHOSAnalyze("galice.root")
34 //                    // set the file you want to analyse
35 //    root [1] a->DrawRecon(1,3)
36 //                    // plot RecObjects, made in event 1, PHOS module 3 
37 //    root [2] a->DrawRecon(1,3,"PHOSRP","another PID")
38 //                    // plot RecObjets made in the event 1, PHOS module 3,
39 //                    // produced in the another reconstruction pass,
40 //                    // which produced PHOS RecParticles ("PHOSRP") with 
41 //                    // title "another PID".
42 //    root [3] a->InvariantMass()
43 //                    // Calculates "REAL" and "MIXED" invariant mass 
44 //                    // distributions of kGAMMA and (kGAMMA+kNEUTRALEM)
45 //                    // and APPENDS this to the file "invmass.root"
46 //    root [4] a->PositionResolution()
47 //                    // calculates two dimentional histos: energy of the primary
48 //                    // photon vs distance betwin incedence point and reconstructed 
49 //                    // poisition. One can analyse the produced file position.root 
50 //                    // with macro PhotonPosition.C
51 //    root [5] a->EnergyResolution()
52 //                    // calculates two dimentional histos: energy of the primary
53 //                    // photon vs energy of the reconstructed particle. One can 
54 //                    // analyse the produced file energy.root 
55 //                    // with macro PhotonEnergy.C
56 //    root [6] a->Contamination()
57 //                    // fills spectra of primary photons and several kinds of 
58 //                    // reconstructed particles, so that analyzing them one can 
59 //                    // estimate conatmination, efficiency of registration etc.
60 //*--
61 //*-- Author: Dmitri Peressounko (SUBATECH & RRC Kurchatov Institute)
62 //////////////////////////////////////////////////////////////////////////////
63
64
65 // --- ROOT system ---
66
67 #include "TFile.h"
68 #include "TH1.h"
69 #include "TH2.h"
70 #include "TH2.h"
71 #include "TParticle.h"
72 #include "TClonesArray.h"
73 #include "TTree.h"
74 #include "TMath.h"
75 #include "TROOT.h"
76 #include "TFolder.h"
77
78 // --- Standard library ---
79
80 #include <iomanip.h>
81
82 // --- AliRoot header files ---
83
84 #include "AliRun.h"
85 #include "AliPHOSv1.h"
86 #include "AliPHOSAnalyze.h"
87 #include "AliPHOSDigit.h"
88 #include "AliPHOSSDigitizer.h"
89 #include "AliPHOSTrackSegment.h"
90 #include "AliPHOSRecParticle.h"
91 #include "AliPHOSCpvRecPoint.h"
92 #include "AliPHOSGetter.h"
93
94
95 ClassImp(AliPHOSAnalyze)
96
97 //____________________________________________________________________________
98   AliPHOSAnalyze::AliPHOSAnalyze()
99 {
100   // default ctor (useless)
101   fCorrection = 1.2 ;  //Value calculated for default parameters of reconstruction  
102 }
103
104 //____________________________________________________________________________
105 AliPHOSAnalyze::AliPHOSAnalyze(Text_t * fileName)
106 {
107   // ctor: analyze events from root file "name"
108   ffileName = fileName ;
109   fCorrection = 1.05 ;  //Value calculated for default parameters of reconstruction   
110 }
111
112 //____________________________________________________________________________
113 AliPHOSAnalyze::AliPHOSAnalyze(const AliPHOSAnalyze & ana)
114 {
115   // copy ctor
116   ( (AliPHOSAnalyze &)ana ).Copy(*this) ;
117 }
118
119 //____________________________________________________________________________
120 AliPHOSAnalyze::~AliPHOSAnalyze()
121 {
122   // dtor
123
124 }
125 //____________________________________________________________________________
126 void AliPHOSAnalyze::DrawRecon(Int_t Nevent,Int_t Nmod,const char * branchName,const char* branchTitle){
127   //Draws pimary particles and reconstructed 
128   //digits, RecPoints, RecPartices etc 
129   //for event Nevent in the module Nmod.
130
131   //========== Create ObjectGetter
132   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
133   if(Nevent >= gAlice->TreeE()->GetEntries() ) {
134     cout << "There is no event " << Nevent << ", only " << gAlice->TreeE()->GetEntries() << "events available " <<endl ;
135     return ;
136   }
137   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
138   gime->Event(Nevent);
139
140   Int_t nx = phosgeom->GetNPhi() ;
141   Int_t nz = phosgeom->GetNZ() ;
142   Float_t * cri= phosgeom->GetEMCAGeometry()->GetCrystalHalfSize() ;
143   Float_t x = nx*cri[0] ;
144   Float_t z = nz*cri[2] ;
145   Int_t nxCPV = (Int_t) (nx*phosgeom->GetPadSizePhi()/(2.*cri[0])) ;
146   Int_t nzCPV = (Int_t) (nz*phosgeom->GetPadSizeZ()/(2.*cri[2])) ;
147   
148   TH2F * emcDigits = (TH2F*) gROOT->FindObject("emcDigits") ;
149   if(emcDigits)
150     emcDigits->Delete() ;
151   emcDigits = new TH2F("emcDigits","EMC digits",  nx,-x,x,nz,-z,z);
152   TH2F * emcSdigits =(TH2F*) gROOT->FindObject("emcSdigits") ;  
153   if(emcSdigits)
154     emcSdigits->Delete() ;
155   emcSdigits = new TH2F("emcSdigits","EMC sdigits", nx,-x,x,nz,-z,z);
156   TH2F * emcRecPoints = (TH2F*)gROOT->FindObject("emcRecPoints") ; 
157   if(emcRecPoints)
158     emcRecPoints->Delete() ;
159   emcRecPoints = new TH2F("emcRecPoints","EMC RecPoints",nx,-x,x,nz,-z,z);
160   TH2F * cpvSdigits =(TH2F*) gROOT->FindObject("cpvSdigits") ;
161   if(cpvSdigits)
162     cpvSdigits->Delete() ;
163   cpvSdigits = new TH2F("cpvSdigits","CPV sdigits", nx,-x,x,nz,-z,z);
164   TH2F * cpvDigits = (TH2F*)gROOT->FindObject("cpvDigits") ;
165   if(cpvDigits)
166     cpvDigits->Delete() ;
167   cpvDigits = new TH2F("cpvDigits","CPV digits",   nxCPV,-x,x,nzCPV,-z,z) ;
168   TH2F * cpvRecPoints= (TH2F*)gROOT->FindObject("cpvRecPoints") ; 
169   if(cpvRecPoints)
170     cpvRecPoints->Delete() ;
171   cpvRecPoints = new TH2F("cpvRecPoints","CPV RecPoints",    nxCPV,-x,x,nzCPV,-z,z) ;
172
173   TH2F * phot = (TH2F*)gROOT->FindObject("phot") ;
174   if(phot)
175     phot->Delete() ;
176   phot = new TH2F("phot","Primary Photon",  nx,-x,x,nz,-z,z);
177   TH2F * recPhot = (TH2F*)gROOT->FindObject("recPhot") ; 
178   if(recPhot)
179     recPhot->Delete() ;
180   recPhot = new TH2F("recPhot","RecParticles with primary Photon",nx,-x,x,nz,-z,z);
181   
182   
183   //Plot Primary Particles
184   const TParticle * primary ;
185   Int_t iPrimary ;
186   for ( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++)
187     {
188       primary = gime->Primary(iPrimary) ;
189       Int_t primaryType = primary->GetPdgCode() ;
190 //       if( (primaryType == 211)||(primaryType == -211)||(primaryType == 2212)||(primaryType == -2212)
191 //        ||(primaryType == 11)||(primaryType == -11) ) {
192 //         Int_t moduleNumber ;
193 //         Double_t primX, primZ ;
194 //         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
195 //         if(moduleNumber==Nmod)
196 //           charg->Fill(primZ,primX,primary->Energy()) ;
197 //       }
198       if( primaryType == 22 ) {
199         Int_t moduleNumber ;
200         Double_t primX, primZ ;
201         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
202         if(moduleNumber==Nmod) 
203           phot->Fill(primZ,primX,primary->Energy()) ;
204       }
205 //       else{
206 //         if( primaryType == -2112 ) {
207 //           Int_t moduleNumber ;
208 //           Double_t primX, primZ ;
209 //           phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
210 //           if(moduleNumber==Nmod)
211 //             nbar->Fill(primZ,primX,primary->Energy()) ;
212 //         }
213 //       }
214     }  
215
216   
217   Int_t iSDigit ;
218   AliPHOSDigit * sdigit ;
219   const TClonesArray * sdigits = gime->SDigits() ;
220   Int_t nsdig[5] = {0,0,0,0,0} ;
221   if(sdigits){
222     for(iSDigit = 0; iSDigit < sdigits->GetEntriesFast() ; iSDigit++)
223       {
224         sdigit = (AliPHOSDigit *) sdigits->At(iSDigit) ;
225         Int_t relid[4];
226         phosgeom->AbsToRelNumbering(sdigit->GetId(), relid) ;
227         Float_t x,z ;
228         phosgeom->RelPosInModule(relid,x,z) ;
229         Float_t e = gime->SDigitizer()->Calibrate(sdigit->GetAmp()) ;
230         nsdig[relid[0]-1]++ ;
231         if(relid[0]==Nmod){
232           if(relid[1]==0)  //EMC
233             emcSdigits->Fill(x,z,e) ;
234           if( relid[1]!=0 )
235             cpvSdigits->Fill(x,z,e) ;
236         }
237       }
238   }
239   cout << "Number of EMC + CPV SDigits per module: " <<endl ;
240   cout << nsdig[0] << " " << nsdig[1] << " " << nsdig[2] << " " << nsdig[3]<< " " << nsdig[4] << endl ;
241   cout << endl ;
242
243
244   //Plot digits
245   Int_t iDigit ;
246   AliPHOSDigit * digit ;
247   const TClonesArray * digits = gime->Digits(); 
248   if(digits) {
249     for(iDigit = 0; iDigit < digits->GetEntriesFast(); iDigit++)
250       {
251         digit = (AliPHOSDigit *) digits->At(iDigit) ;
252         Int_t relid[4];
253         phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
254         Float_t x,z ;
255         phosgeom->RelPosInModule(relid,x,z) ;
256         Float_t e = gime->SDigitizer()->Calibrate(digit->GetAmp()) ;
257         if(relid[0]==Nmod){
258           if(relid[1]==0)  //EMC
259             emcDigits->Fill(x,z,e) ;
260           if( relid[1]!=0 )
261             cpvDigits->Fill(x,z,e) ;
262         }
263       }
264   }
265   
266   
267   //Plot RecPoints
268   Int_t irecp ;
269   TVector3 pos ;
270   TObjArray * emcrp = gime->EmcRecPoints() ;
271   if(emcrp) {
272     for(irecp = 0; irecp < emcrp->GetEntriesFast() ; irecp ++){
273       AliPHOSEmcRecPoint * emc = (AliPHOSEmcRecPoint *) emcrp->At(irecp) ;
274       if(emc->GetPHOSMod()==Nmod){
275         emc->GetLocalPosition(pos) ;
276         emcRecPoints->Fill(pos.X(),pos.Z(),emc->GetEnergy());
277       }
278     }
279   }
280   
281   TObjArray * cpvrp = gime->CpvRecPoints() ;
282   if(cpvrp) {
283     for(irecp = 0; irecp < cpvrp->GetEntriesFast() ; irecp ++){
284       AliPHOSRecPoint * cpv = (AliPHOSCpvRecPoint *) cpvrp->At(irecp) ;
285       if(cpv->GetPHOSMod()==Nmod){
286         cpv->GetLocalPosition(pos) ;
287         cpvRecPoints->Fill(pos.X(),pos.Z(),cpv->GetEnergy());
288       }
289     }
290   }
291     
292   //Plot RecParticles
293   AliPHOSRecParticle * recParticle ;
294   Int_t iRecParticle ;
295   TClonesArray * rp = gime->RecParticles() ;
296   TClonesArray * ts = gime->TrackSegments() ;
297   if(rp && ts && emcrp) {
298     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ; iRecParticle++ )
299       {
300         recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
301         Int_t moduleNumberRec ;
302         Double_t recX, recZ ;
303         phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
304         if(moduleNumberRec == Nmod){
305           
306           Double_t minDistance = 5. ;
307           Int_t closestPrimary = -1 ;   
308
309           //extract list of primaries: it is stored at EMC RecPoints
310           Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
311           Int_t numberofprimaries ;
312           Int_t * listofprimaries  = ((AliPHOSRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
313           Int_t index ;
314           const TParticle * primary ;
315           Double_t distance = minDistance ;
316           
317           for ( index = 0 ; index < numberofprimaries ; index++){
318             primary = gime->Primary(listofprimaries[index]) ;
319             Int_t moduleNumber ;
320             Double_t primX, primZ ;
321             phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
322             if(moduleNumberRec == moduleNumber)
323               distance = TMath::Sqrt((recX-primX)*(recX-primX)+(recZ-primZ)*(recZ-primZ) ) ;
324             if(minDistance > distance)
325               {
326                 minDistance = distance ;
327                 closestPrimary = listofprimaries[index] ;
328               }
329           }
330           
331           if(closestPrimary >=0 ){
332             
333             Int_t primaryType = gime->Primary(closestPrimary)->GetPdgCode() ;
334             
335             if(primaryType==22)
336               recPhot->Fill(recZ,recX,recParticle->Energy()) ;
337 //          else
338 //            if(primaryType==-2112)
339 //              recNbar->Fill(recZ,recX,recParticle->Energy()) ; 
340           }
341         }
342       }
343
344   }
345   
346   //Plot made histograms
347   emcSdigits->Draw("box") ;
348   emcDigits->SetLineColor(5) ;
349   emcDigits->Draw("boxsame") ;
350   emcRecPoints->SetLineColor(2) ;
351   emcRecPoints->Draw("boxsame") ;
352   cpvSdigits->SetLineColor(1) ;
353   cpvSdigits->Draw("boxsame") ;
354   
355 }
356 //____________________________________________________________________________
357 void AliPHOSAnalyze::Ls(){
358   //lists branches and titles of PHOS-related branches of TreeR, TreeD, TreeS
359   
360   AliPHOSGetter::GetInstance(ffileName.Data()) ;
361
362   Int_t ibranch;
363   TObjArray * branches; 
364   
365   branches = gAlice->TreeS()->GetListOfBranches() ;
366  
367   cout << "TreeS: " << endl ;
368   for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
369     TBranch * branch=(TBranch *) branches->At(ibranch) ;
370     if(strstr(branch->GetName(),"PHOS") )
371       cout << "       " << branch->GetName() << "     " << branch->GetTitle() << endl ;
372   }
373   cout << endl ;
374  
375   branches = gAlice->TreeD()->GetListOfBranches() ;
376   
377   cout << "TreeD: " << endl ;
378   for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
379     TBranch * branch=(TBranch *) branches->At(ibranch) ;
380     if(strstr(branch->GetName(),"PHOS") )
381       cout << "       " << branch->GetName() << "     " << branch->GetTitle() << endl ;
382   }
383   cout << endl ;
384   
385
386   branches = gAlice->TreeR()->GetListOfBranches() ;
387   
388   cout << "TreeR: " << endl ;
389   for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
390     TBranch * branch=(TBranch *) branches->At(ibranch) ;
391     if(strstr(branch->GetName(),"PHOS") )
392       cout << "       " << branch->GetName() << "     " << branch->GetTitle() << endl ;
393   }
394   cout << endl ;
395
396
397 }
398 //____________________________________________________________________________
399  void AliPHOSAnalyze::InvariantMass(const char* branchTitle)    
400 {
401   // Calculates Real and Mixed invariant mass distributions
402   AliPHOSGetter * gime  = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
403
404   Int_t nMixedEvents = 4 ; //# of events used for calculation of 'mixed' distribution 
405
406   //opening file
407   TFile * mfile = new TFile("invmass.root","update");
408   
409   //========== Reading /Booking Histograms
410   TH2D * hRealEM = 0 ;
411   hRealEM = (TH2D*) mfile->Get("hRealEM") ;
412   if(hRealEM == 0) 
413     hRealEM = new TH2D("hRealEM",   "Real for EM particles",      250,0.,1.,40,0.,4.) ;
414   TH2D * hRealPhot = 0 ;
415
416   hRealPhot = (TH2D*)mfile->Get("hRealPhot");
417   if(hRealPhot == 0)
418     hRealPhot = new TH2D("hRealPhot", "Real for kPhoton particles", 250,0.,1.,40,0.,4.) ;
419
420   TH2D * hMixedEM = 0 ;
421   hMixedEM = (TH2D*) mfile->Get("hMixedEM") ;
422   if(hMixedEM == 0)
423     hMixedEM = new TH2D("hMixedEM",  "Mixed for EM particles",     250,0.,1.,40,0.,4.) ;
424
425   TH2D * hMixedPhot = 0 ;
426   hMixedPhot = (TH2D*) mfile->Get("hMixedPhot") ;
427   if(hMixedPhot == 0)
428     hMixedPhot = new TH2D("hMixedPhot","Mixed for kPhoton particles",250,0.,1.,40,0.,4.) ;
429   
430
431   //reading event and copyng it to TConesArray of all photons
432
433   TClonesArray * allRecParticleList  = new TClonesArray("AliPHOSRecParticle", 1000) ;
434   Int_t * nRecParticles = new Int_t[nMixedEvents] ; // to mark boundaries of each event in the total list
435   for(Int_t index = 0; index < nMixedEvents; index ++)
436     nRecParticles[index] = 0 ;
437   Int_t iRecPhot = 0 ;                // number of EM particles in total list
438   
439   //scan over all events
440   Int_t event ;
441   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
442   //  for(event = 0; event < gime->MaxEvent(); event++  ){
443   for(event = 0; event < maxevent; event++  ){
444     gime->Event(event,"R");  //will read only TreeR 
445     
446     //copy EM RecParticles to the "total" list        
447     const AliPHOSRecParticle * recParticle ;
448     Int_t iRecParticle ;
449     TClonesArray * rp = gime->RecParticles() ;
450     if(!rp){
451       cout << "AliPHOSAnalyze::InvariantMass --> Can't find RecParticles " << endl ;
452       return ;
453     }
454
455     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ )
456       {
457         recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
458         if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
459            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW))
460           new( (*allRecParticleList)[iRecPhot++] ) AliPHOSRecParticle(*recParticle) ;
461       }
462     
463     Int_t mevent = event%nMixedEvents ; //event number in the "mixed" cicle
464     nRecParticles[mevent] = iRecPhot-1 ;  
465     
466     //check, if it is time to calculate invariant mass?
467     Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
468     if((mevent == 0) && (event +1 == maxevent)){
469       
470       //   if((mevent == 0) && (event +1 == gime->MaxEvent())){
471       
472       //calculate invariant mass:
473       Int_t irp1,irp2 ;
474       Int_t nCurEvent = 0 ;
475       
476       for(irp1 = 0; irp1 < allRecParticleList->GetEntries()-1; irp1++){
477         AliPHOSRecParticle * rp1 = (AliPHOSRecParticle *)allRecParticleList->At(irp1) ;
478         
479         for(irp2 = irp1+1; irp2 < allRecParticleList->GetEntries(); irp2++){
480           AliPHOSRecParticle * rp2 = (AliPHOSRecParticle *)allRecParticleList->At(irp2) ;
481           
482           Double_t invMass ;
483           invMass = (rp1->Energy()+rp2->Energy())*(rp1->Energy()+rp2->Energy())-
484             (rp1->Px()+rp2->Px())*(rp1->Px()+rp2->Px())-
485             (rp1->Py()+rp2->Py())*(rp1->Py()+rp2->Py())-
486             (rp1->Pz()+rp2->Pz())*(rp1->Pz()+rp2->Pz()) ;
487           
488           if(invMass> 0)
489             invMass = TMath::Sqrt(invMass);
490           
491           Double_t pt ; 
492           pt = TMath::Sqrt((rp1->Px()+rp2->Px() )*( rp1->Px()+rp2->Px() ) + 
493                            (rp1->Py()+rp2->Py() )*( rp1->Py()+rp2->Py() ) );
494           
495           if(irp1 > nRecParticles[nCurEvent])
496             nCurEvent++;
497           
498           if(irp2 <= nRecParticles[nCurEvent]){ //'Real' event
499             hRealEM->Fill(invMass,pt);
500             if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
501                (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
502               hRealPhot->Fill(invMass,pt);
503           }
504           else{
505             hMixedEM->Fill(invMass,pt);
506             if((rp1->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)&&
507                (rp2->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST) )
508               hMixedPhot->Fill(invMass,pt);
509           } //real-mixed
510           
511       } //loop over second rp
512       }//loop over first rp
513       
514       //Make some cleanings 
515       for(Int_t index = 0; index < nMixedEvents; index ++)
516         nRecParticles[index] = 0 ;
517       iRecPhot = 0 ;              
518       allRecParticleList->Clear() ;
519       
520     }
521   }
522   delete allRecParticleList ;
523   
524   //writing output
525   mfile->cd();
526   
527   hRealEM->Write(0,kOverwrite) ;
528   hRealPhot->Write(0,kOverwrite) ;
529   hMixedEM->Write(0,kOverwrite) ;
530   hMixedPhot->Write(0,kOverwrite) ;
531   
532   mfile->Write();
533   mfile->Close();
534   delete mfile ;
535   delete nRecParticles;
536
537 }
538
539 //____________________________________________________________________________
540  void AliPHOSAnalyze::EnergyResolution(const char * branchTitle)    
541 {
542   //fills two dimentional histo: energy of primary vs. energy of reconstructed
543
544   TH2F * hAllEnergy = 0 ;  //all reconstructed with primary photon   
545   TH2F * hPhotEnergy= 0 ;  //kGamma  with primary photon   
546   TH2F * hEMEnergy  = 0 ;  //electromagnetic with primary photon   
547
548   //opening file and reading histograms if any
549   TFile * efile = new TFile("energy.root","update");
550
551   hAllEnergy = (TH2F*)efile->Get("hAllEnergy") ;
552   if(hAllEnergy == 0)
553     hAllEnergy  = new TH2F("hAllEnergy",  "Energy of any RP with primary photon",100, 0., 5., 100, 0., 5.);
554
555   hPhotEnergy =(TH2F*) efile->Get("hPhotEnergy") ;
556   if(hPhotEnergy == 0)
557     hPhotEnergy = new TH2F("hPhotEnergy", "Energy of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
558
559   hEMEnergy =(TH2F*) efile->Get("hEMEnergy");
560   if(hEMEnergy == 0)
561     hEMEnergy   = new TH2F("hEMEnergy",   "Energy of EM with primary photon",    100, 0., 5., 100, 0., 5.);
562
563
564   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
565   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
566
567   Int_t ievent;
568   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
569   for ( ievent=0; ievent < maxevent ; ievent++){
570
571     //read the current event
572     gime->Event(ievent) ;
573
574     const AliPHOSRecParticle * recParticle ;
575     Int_t iRecParticle ;
576     TClonesArray * rp = gime->RecParticles() ;
577     if(!rp) {
578       cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent 
579            << ",  Can't find RecParticles " << endl ;
580       return ;
581     }
582     TClonesArray * ts = gime->TrackSegments() ;
583     if(!ts) {
584       cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent 
585            << ",  Can't find TrackSegments " << endl ;
586       return ;
587     }
588     TObjArray * emcrp = gime->EmcRecPoints() ;
589     if(!emcrp){
590       cout << "AliPHOSAnalyze::EnergyResolution --> Event " <<ievent 
591            << ",  Can't find EmcRecPoints " << endl ;
592       return ;
593     }
594       
595     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast() ;iRecParticle++ ){
596       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
597       
598       //find the closest primary
599       Int_t moduleNumberRec ;
600       Double_t recX, recZ ;
601       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
602       
603       Double_t minDistance  = 100. ;
604       Int_t closestPrimary = -1 ;
605       
606       //extract list of primaries: it is stored at EMC RecPoints
607       Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
608       Int_t numberofprimaries ;
609       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint*) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
610       
611       Int_t index ;
612       const TParticle * primary ;
613       Double_t distance = minDistance ;
614       Double_t dX, dZ; 
615       Double_t dXmin = 0.; 
616       Double_t dZmin = 0. ;
617       for ( index = 0 ; index < numberofprimaries ; index++){
618         primary = gime->Primary(listofprimaries[index]) ;
619         Int_t moduleNumber ;
620         Double_t primX, primZ ;
621         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
622         if(moduleNumberRec == moduleNumber) {
623           dX = recX - primX;
624           dZ = recZ - primZ;
625           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
626           if(minDistance > distance) {
627             minDistance = distance ;
628             dXmin = dX;
629             dZmin = dZ;
630             closestPrimary = listofprimaries[index] ;
631           }
632         }
633       }
634
635       //if found primary, fill histograms
636       if(closestPrimary >=0 ){
637         const TParticle * primary = gime->Primary(closestPrimary) ;
638         if(primary->GetPdgCode() == 22){
639           hAllEnergy->Fill(primary->Energy(), recParticle->Energy()) ;
640           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
641             hPhotEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
642             hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
643           }
644           else
645             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
646               hEMEnergy->Fill(primary->Energy(), recParticle->Energy() ) ; 
647         }
648       }
649     }
650   }
651
652   //write filled histograms
653   efile->cd() ;
654   hAllEnergy->Write(0,kOverwrite) ;
655   hPhotEnergy->Write(0,kOverwrite) ;
656   hEMEnergy->Write(0,kOverwrite)  ;
657   //  efile->Write() ;
658   efile->Close() ;
659   delete efile ;
660
661 }
662 //____________________________________________________________________________
663 void AliPHOSAnalyze::PositionResolution(const char * branchTitle)    
664 {
665   //fills two dimentional histo: energy vs. primary - reconstructed distance  
666
667
668
669   TH2F * hAllPosition  = 0;     // Position of any RP with primary photon
670   TH2F * hPhotPosition = 0;    // Position of kGAMMA with primary photon
671   TH2F * hEMPosition   = 0;      // Position of EM with primary photon
672
673   TH1F * hAllPositionX = 0;    // X-Position Resolution of photons with photon primary
674   TH1F * hAllPositionZ = 0;    // Z-Position Resolution of photons with photon primary
675
676
677   //opening file and reading histograms if any
678   TFile * pfile = new TFile("position.root","update");
679
680   hAllPosition = (TH2F*)pfile->Get("hAllPosition");
681   if(hAllPosition == 0)
682     hAllPosition  = new TH2F("hAllPosition",  
683                              "Position of any RP with primary photon",100, 0., 5., 100, 0., 5.);
684   hPhotPosition= (TH2F*)pfile->Get("hPhotPosition");
685   if(hPhotPosition == 0)
686     hPhotPosition = new TH2F("hPhotPosition", 
687                              "Position of kGAMMA with primary photon",100, 0., 5., 100, 0., 5.);
688   hEMPosition= (TH2F*)pfile->Get("hEMPosition") ;
689   if(hEMPosition == 0)
690     hEMPosition   = new TH2F("hEMPosition",   
691                              "Position of EM with primary photon",    100, 0., 5., 100, 0., 5.);
692   hAllPositionX = (TH1F*)pfile->Get("hAllPositionX") ;                     
693   if(hAllPositionX == 0)
694     hAllPositionX = new TH1F("hAllPositionX", 
695                              "Delta X of any RP with primary photon",100, -2., 2.);
696   hAllPositionZ =(TH1F*) pfile->Get("hAllPositionZ") ;
697   if(hAllPositionZ == 0)
698     hAllPositionZ = new TH1F("hAllPositionZ", 
699                              "Delta X of any RP with primary photon",100, -2., 2.);
700
701
702   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),branchTitle) ;
703   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
704
705   Int_t ievent;
706   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
707   for ( ievent=0; ievent < maxevent ; ievent++){
708     
709     //read the current event
710     gime->Event(ievent) ;
711     TClonesArray * rp = gime->RecParticles() ;
712     if(!rp) {
713       cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent 
714            << ",  Can't find RecParticles " << endl ;
715       return ;
716     }
717     TClonesArray * ts = gime->TrackSegments() ;
718     if(!ts) {
719       cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent 
720            << ",  Can't find TrackSegments " << endl ;
721       return ;
722     }
723     TObjArray * emcrp = gime->EmcRecPoints() ;
724     if(!emcrp){
725       cout << "AliPHOSAnalyze::PositionResolution --> Event " <<ievent 
726            << ",  Can't find EmcRecPoints " << endl ;
727       return ;
728     }
729     
730  
731     const AliPHOSRecParticle * recParticle ;
732     Int_t iRecParticle ;
733     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
734       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
735       
736       //find the closest primary
737       Int_t moduleNumberRec ;
738       Double_t recX, recZ ;
739       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
740       
741       Double_t minDistance  = 100. ;
742       Int_t closestPrimary = -1 ;
743       
744       //extract list of primaries: it is stored at EMC RecPoints
745       Int_t emcIndex = ((AliPHOSTrackSegment*) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
746       Int_t numberofprimaries ;
747       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
748
749       Int_t index ;
750       const TParticle * primary ;
751       Double_t distance = minDistance ;
752       Double_t dX = 1000; // incredible number
753       Double_t dZ = 1000; // for the case if no primary will be found
754       Double_t dXmin = 0.; 
755       Double_t dZmin = 0. ;
756       for ( index = 0 ; index < numberofprimaries ; index++){
757         primary = gime->Primary(listofprimaries[index]) ;
758         Int_t moduleNumber ;
759         Double_t primX, primZ ;
760         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
761         if(moduleNumberRec == moduleNumber) {
762           dX = recX - primX;
763           dZ = recZ - primZ;
764           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
765           if(minDistance > distance) {
766             minDistance = distance ;
767             dXmin = dX;
768             dZmin = dZ;
769             closestPrimary = listofprimaries[index] ;
770           }
771         }
772       }
773       
774       //if found primary, fill histograms
775       if(closestPrimary >=0 ){
776         const TParticle * primary = gime->Primary(closestPrimary) ;
777         if(primary->GetPdgCode() == 22){
778           hAllPosition->Fill(primary->Energy(), minDistance) ;
779           hAllPositionX->Fill(primary->Energy(), dX) ;
780           hAllPositionZ->Fill(primary->Energy(), dZ) ;
781           if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
782             hPhotPosition->Fill(primary->Energy(), minDistance ) ; 
783             hEMPosition->Fill(primary->Energy(), minDistance ) ; 
784           }
785           else
786             if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)
787               hEMPosition->Fill(primary->Energy(), minDistance ) ; 
788         }
789       }
790     }
791   }
792   
793   //Write output histgrams
794   pfile->cd() ;
795   hAllPosition->Write(0,kOverwrite) ;
796   hAllPositionX->Write(0,kOverwrite) ;
797   hAllPositionZ->Write(0,kOverwrite) ;
798   hPhotPosition->Write(0,kOverwrite) ;
799   hEMPosition->Write(0,kOverwrite) ;
800   pfile->Write() ;
801   pfile->Close() ;
802   delete pfile ;
803
804
805 }
806 //____________________________________________________________________________
807 void AliPHOSAnalyze::Contamination(const char* RecPointsTitle){
808 // fills spectra of primary photons and several kinds of 
809 // reconstructed particles, so that analyzing them one can 
810 // estimate conatmination, efficiency of registration etc.
811
812   //define several general histograms
813   TH1F * hPrimary = 0;   //spectrum (P_t distribution) of primary photons         
814   TH1F * hAllRP   = 0;   //spectrum of all RecParticles in PHOS
815   TH1F * hPhot    = 0;   //spectrum of kGAMMA RecParticles
816   TH1F * hShape   = 0;   //spectrum of all EM RecParticles
817   TH1F * hVeto    = 0;   //spectrum of all neutral RecParticles
818
819   //Now separate histograms in accoradance with primary
820   //primary - photon
821   TH1F * hPhotReg = 0;   //Registeres as photon
822   TH1F * hPhotEM  = 0;   //Registered as EM       
823
824   //primary - n
825   TH1F * hNReg = 0;   //Registeres as photon          
826   TH1F * hNEM  = 0;   //Registered as EM            
827
828   //primary - nBar
829   TH1F * hNBarReg = 0;   //Registeres as photon
830   TH1F * hNBarEM  = 0;   //Registered as EM          
831
832   //primary - charged hadron (pBar excluded)
833   TH1F * hChargedReg = 0;   //Registeres as photon  
834   TH1F * hChargedEM  = 0;   //Registered as EM           
835
836   //primary - pBar
837   TH1F * hPbarReg = 0;   //Registeres as photon  
838   TH1F * hPbarEM  = 0;   //Registered as EM 
839
840
841   //Reading histograms from the file
842   TFile * cfile = new TFile("contamination.root","update") ;
843
844   //read general histograms
845   hPrimary = (TH1F*) cfile->Get("hPrimary") ;
846   if(hPrimary == 0)
847     hPrimary= new TH1F("hPrimary", "Primary photon spectrum", 100, 0., 5.);
848   hAllRP = (TH1F*)cfile->Get("hAllRP") ;
849   if(hAllRP == 0)
850     hAllRP = new TH1F("hAllRP","All Reconstructed particles", 100, 0., 5.);
851   hPhot  = (TH1F*)cfile->Get("hPhot") ;
852   if(hPhot == 0)
853     hPhot = new TH1F("hPhot","All kGAMMA RecParticles",100, 0., 5.);
854   hShape = (TH1F*) cfile->Get("hShape") ;
855   if(hShape == 0)
856     hShape = new TH1F("hShape","All particles with EM shower",100, 0., 5.);
857   hVeto= (TH1F*)cfile->Get("hVeto") ;
858   if(hVeto == 0) 
859     hVeto  = new TH1F("hVeto", "All uncharged particles",      100, 0., 5.);
860
861
862   //primary - photon
863   hPhotReg = (TH1F*)cfile->Get("hPhotReg");
864   if(hPhotReg == 0)
865     hPhotReg   = new TH1F("hPhotReg","Photon registered as photon",100, 0., 5.);
866   hPhotEM  =(TH1F*)cfile->Get("hPhotEM");
867   if(hPhotEM== 0)
868     hPhotEM   = new TH1F("hPhotEM",  "Photon registered as EM", 100, 0., 5.);
869
870   //primary - n
871   hNReg = (TH1F*)cfile->Get("hNReg");
872   if(hNReg== 0)
873    hNReg      = new TH1F("hNReg",   "N registered as photon",              100, 0., 5.);
874   hNEM  = (TH1F*)cfile->Get("hNEM"); 
875   if(hNEM== 0)
876     hNEM      = new TH1F("hNEM",    "N registered as EM",      100, 0., 5.);
877
878   //primary - nBar
879   hNBarReg =(TH1F*)cfile->Get("hNBarReg");
880   if(hNBarReg== 0)
881    hNBarReg   = new TH1F("hNBarReg", "NBar registered as photon",           100, 0., 5.);
882   hNBarEM  =(TH1F*)cfile->Get("hNBarEM"); 
883   if(hNBarEM== 0)
884     hNBarEM   = new TH1F("hNBarEM",  "NBar registered as EM",   100, 0., 5.);
885
886   //primary - charged hadron (pBar excluded)
887   hChargedReg = (TH1F*)cfile->Get("hChargedReg");
888   if(hChargedReg== 0)
889     hChargedReg= new TH1F("hChargedReg", "Charged hadron registered as photon",100, 0., 5.);
890   hChargedEM  = (TH1F*)cfile->Get("hChargedEM"); 
891   if(hChargedEM== 0)
892     hChargedEM= new TH1F("hChargedEM","Charged registered as EM",100, 0., 5.);
893  
894   //primary - pBar
895   hPbarReg = (TH1F*)cfile->Get("hPbarReg");
896   if(hPbarReg== 0)
897     hPbarReg= new TH1F("hPbarReg", "pBar registered as photon",100, 0., 5.);
898   hPbarEM  = (TH1F*)cfile->Get("hPbarEM");
899   if(hPbarEM== 0)
900     hPbarEM= new TH1F("hPbarEM","Pbar registered as EM",100, 0., 5.);
901   
902
903   //Now make some initializations
904
905   Int_t counter[8][5] ;      //# of registered particles 
906   Int_t i1,i2 ;
907   for(i1 = 0; i1<8; i1++)
908     for(i2 = 0; i2<5; i2++)
909       counter[i1][i2] = 0 ;
910
911
912
913   AliPHOSGetter * gime = AliPHOSGetter::GetInstance(ffileName.Data(),RecPointsTitle) ;
914   const AliPHOSGeometry * phosgeom = gime->PHOSGeometry() ; 
915   
916   Int_t ievent;
917   Int_t maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
918   for ( ievent=0; ievent < maxevent ; ievent++){
919     
920     gime->Event(ievent) ;
921     
922     TClonesArray * rp = gime->RecParticles() ;
923     if(!rp) {
924       cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent 
925            << ",  Can't find RecParticles " << endl ;
926       return ;
927     }
928     TClonesArray * ts = gime->TrackSegments() ;
929     if(!ts) {
930       cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent 
931            << ",  Can't find TrackSegments " << endl ;
932       return ;
933     }
934     TObjArray * emcrp = gime->EmcRecPoints() ;
935     if(!emcrp){
936       cout << "AliPHOSAnalyze::Contamination --> Event " <<ievent 
937            << ",  Can't find EmcRecPoints " << endl ;
938       return ;
939     }
940     
941     
942     //=========== Make spectrum of the primary photons
943     const TParticle * primary ;
944     Int_t iPrimary ;
945     for( iPrimary = 0 ; iPrimary < gime->NPrimaries() ; iPrimary++){
946       primary = gime->Primary(iPrimary) ;
947       Int_t primaryType = primary->GetPdgCode() ;
948       if( primaryType == 22 ) {
949         //check, if photons folls onto PHOS
950         Int_t moduleNumber ;
951         Double_t primX, primZ ;
952         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
953         if(moduleNumber)
954           hPrimary->Fill(primary->Energy()) ;
955         
956       }
957       
958     }
959     
960     //========== Now scan over RecParticles            
961     const AliPHOSRecParticle * recParticle ;
962     Int_t iRecParticle ;
963     for(iRecParticle = 0; iRecParticle < rp->GetEntriesFast(); iRecParticle++ ){
964       recParticle = (AliPHOSRecParticle *) rp->At(iRecParticle) ;
965       //fill histo spectrum of all RecParticles
966       hAllRP->Fill(CorrectedEnergy(recParticle->Energy())) ;
967       
968       //==========find the closest primary      
969       Int_t moduleNumberRec ;
970       Double_t recX, recZ ;
971       phosgeom->ImpactOnEmc(recParticle->Theta(), recParticle->Phi(), moduleNumberRec, recX, recZ) ;
972       
973       Double_t minDistance  = 100. ;
974       Int_t closestPrimary = -1 ;
975       
976       //extract list of primaries: it is stored at EMC RecPoints
977       Int_t emcIndex = ((AliPHOSTrackSegment *) ts->At(recParticle->GetPHOSTSIndex()))->GetEmcIndex() ;
978       Int_t numberofprimaries ;
979       Int_t * listofprimaries  = ((AliPHOSEmcRecPoint *) emcrp->At(emcIndex))->GetPrimaries(numberofprimaries)  ;
980       Int_t index ;
981       const TParticle * primary ;
982       Double_t distance = minDistance ;
983       Double_t dX, dZ; 
984       Double_t dXmin = 0.; 
985       Double_t dZmin = 0. ;
986       for ( index = 0 ; index < numberofprimaries ; index++){
987         primary = gime->Primary(listofprimaries[index]) ;
988         Int_t moduleNumber ;
989         Double_t primX, primZ ;
990         phosgeom->ImpactOnEmc(primary->Theta(), primary->Phi(), moduleNumber, primX, primZ) ;
991         if(moduleNumberRec == moduleNumber) {
992           dX = recX - primX;
993           dZ = recZ - primZ;
994           distance = TMath::Sqrt(dX*dX + dZ*dZ) ;
995           if(minDistance > distance) {
996             minDistance = distance ;
997             dXmin = dX;
998             dZmin = dZ;
999             closestPrimary = listofprimaries[index] ;
1000           }
1001         }
1002       }
1003       
1004       //===========define the "type" of closest primary
1005       if(closestPrimary >=0 ){
1006         Int_t primaryCode = -1;
1007         const TParticle * primary = gime->Primary(closestPrimary) ;
1008         Int_t primaryType = primary->GetPdgCode() ;
1009         if(primaryType == 22) // photon ?
1010           primaryCode = 0 ;
1011         else
1012           if(primaryType == 2112) // neutron
1013             primaryCode = 1 ; 
1014           else
1015             if(primaryType == -2112) // Anti neutron
1016               primaryCode = 2 ;
1017             else
1018               if(primaryType == -2122) //Anti proton
1019                 primaryCode = 4 ;
1020               else {
1021                 TParticle tempo(*primary) ; 
1022                 if(tempo.GetPDG()->Charge())
1023                   primaryCode = 3 ;
1024               }
1025
1026         //==========Now look at the type of RecParticle
1027         Float_t energy = CorrectedEnergy(recParticle->Energy()) ;
1028         if(recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST){
1029           hPhot->Fill(energy ) ;        
1030           switch(primaryCode){
1031           case 0:
1032             hPhotReg->Fill(energy ) ; 
1033             break ;
1034           case 1:
1035             hNReg->Fill(energy ) ; 
1036             break ;
1037           case 2:
1038             hNBarReg->Fill(energy ) ; 
1039             break ;
1040           case 3:
1041             hChargedReg->Fill(energy ) ;
1042             break ;
1043           case 4:
1044             hPbarReg->Fill(energy ) ;
1045             break ;
1046           default:
1047             break ;
1048           }
1049         }
1050         if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1051            (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMFAST)||
1052            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW)||
1053            (recParticle->GetType() == AliPHOSFastRecParticle::kCHARGEDEMSLOW) ){ //with EM shower
1054           hShape->Fill(energy ) ;
1055           switch(primaryCode){
1056           case 0:
1057             hPhotEM->Fill(energy ) ; 
1058             break ;
1059           case 1:
1060             hNEM->Fill(energy ) ; 
1061             break ;
1062           case 2:
1063             hNBarEM->Fill(energy ) ; 
1064             break ;
1065           case 3:
1066             hChargedEM->Fill(energy ) ; 
1067             break ;
1068           case 4:
1069             hPbarEM->Fill(energy ) ; 
1070             break ;
1071           default:
1072             break ;
1073           }
1074         }
1075         
1076         if((recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMFAST)||
1077            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHAFAST) ||
1078            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALEMSLOW) ||
1079            (recParticle->GetType() == AliPHOSFastRecParticle::kNEUTRALHASLOW) ) //nuetral
1080           hVeto->Fill(energy ) ;
1081         
1082         //fill number of primaries identified as ...
1083         if(primaryCode >= 0) // Primary code defined
1084           counter[recParticle->GetType()][primaryCode]++ ; 
1085         
1086       }
1087       
1088     } // no closest primary found
1089   }     
1090   
1091   
1092   //===================  SaveHistograms
1093   cfile->cd() ;
1094   hPrimary->Write(0,kOverwrite); 
1095   hAllRP->Write(0,kOverwrite);  
1096   hPhot->Write(0,kOverwrite);  
1097   hShape->Write(0,kOverwrite); 
1098   hVeto->Write(0,kOverwrite);  
1099   hPhotReg->Write(0,kOverwrite); 
1100   hPhotEM->Write(0,kOverwrite);   
1101   hNReg ->Write(0,kOverwrite);  
1102   hNEM  ->Write(0,kOverwrite); 
1103   hNBarReg ->Write(0,kOverwrite); 
1104   hNBarEM  ->Write(0,kOverwrite); 
1105   hChargedReg ->Write(0,kOverwrite); 
1106   hChargedEM  ->Write(0,kOverwrite); 
1107   hPbarReg ->Write(0,kOverwrite); 
1108   hPbarEM  ->Write(0,kOverwrite); 
1109   
1110   cfile->Write(0,kOverwrite); 
1111   cfile->Close();
1112   delete cfile ;
1113  
1114   
1115   //print Final Table
1116   maxevent = (Int_t)gAlice->TreeE()->GetEntries() ; 
1117  //  cout << "Resolutions: Analyzed " << gime->MaxEvent() << " event(s)" << endl ;
1118
1119   cout << "Resolutions: Analyzed " << maxevent << " event(s)" << endl ;
1120   cout << endl ;
1121   
1122   cout << "        Primary:    Photon  Neutron  Antineutron  Charged hadron   AntiProton" << endl ; 
1123   cout << "--------------------------------------------------------------------------------" << endl ;
1124   cout << "         kGAMMA: " 
1125        << setw(8) << counter[2][0] << setw(9)  << counter[2][1] << setw(13) << counter[2][2] 
1126        << setw(15)<< counter[2][3] << setw(13) << counter[2][4] << endl ;
1127   cout << "       kGAMMAHA: " 
1128        << setw(8) << counter[3][0] << setw(9)  << counter[3][1] << setw(13) << counter[3][2] 
1129        << setw(15)<< counter[3][3] << setw(13) << counter[3][4] << endl ;
1130   cout << "     kNEUTRALEM: " 
1131        << setw(8) << counter[0][0] << setw(9)  << counter[0][1] << setw(13) << counter[0][2] 
1132        << setw(15)<< counter[0][3] << setw(13) << counter[0][4] << endl ;
1133   cout << "     kNEUTRALHA: " 
1134        << setw(8) << counter[1][0] << setw(9)  << counter[1][1] << setw(13) << counter[1][2] 
1135        << setw(15)<< counter[1][3] << setw(13) << counter[1][4] << endl ;
1136   cout << "      kABSURDEM: " 
1137        << setw(8) << counter[4][0] << setw(9)  << counter[4][1] << setw(13) << counter[4][2] 
1138        << setw(15)<< counter[4][3] << setw(13) << counter[4][4] << endl ;
1139   cout << "      kABSURDHA: " 
1140        << setw(8) << counter[5][0] << setw(9)  << counter[5][1] << setw(13) << counter[5][2] 
1141        << setw(15)<< counter[5][3] << setw(13) << counter[5][4] << endl ;
1142   cout << "      kELECTRON: " 
1143        << setw(8) << counter[6][0] << setw(9)  << counter[6][1] << setw(13) << counter[6][2] 
1144        << setw(15)<< counter[6][3] << setw(13) << counter[6][4] << endl ; 
1145   cout << "     kCHARGEDHA: " 
1146        << setw(8) << counter[7][0] << setw(9)  << counter[7][1] << setw(13) << counter[7][2] 
1147        << setw(15)<< counter[7][3] << setw(13) << counter[7][4] << endl ;
1148   cout << "--------------------------------------------------------------------------------" << endl ;
1149       
1150   Int_t totalInd = 0 ;
1151   for(i1 = 0; i1<8; i1++)
1152     for(i2 = 0; i2<5; i2++)
1153       totalInd+=counter[i1][i2] ;
1154   cout << "Indentified particles: " << totalInd << endl ;
1155   
1156 }
1157
1158
1159