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