Including TFile.h
[u/mrichter/AliRoot.git] / PHOS / AliPHOSGammaJet.cxx
index c724b7c2074e310c4d446c6e906f2d13a6ac1176..b04234927419d0e88325bcd9be7fd072c5b48602 100644 (file)
  **************************************************************************/
 /* $Id$ */
 
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.10  2006/01/23 18:04:08  hristov
+ * Removing meaningless const
+ *
+ * Revision 1.9  2006/01/12 16:23:26  schutz
+ * ESD is properly read with methods of macros/AliReadESD.C copied in it
+ *
+ * Revision 1.8  2005/12/20 07:08:32  schutz
+ * corrected error in call AliReadESD
+ *
+ * Revision 1.6  2005/05/28 14:19:04  schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
 //_________________________________________________________________________
 // Class for the analysis of gamma-jet correlations 
-//   
+//  Basically it seaches for a prompt photon in the PHOS acceptance, 
+//  if so we construct a jet around the highest pt particle in the opposite 
+//  side in azimuth, inside the TPC and EMCAL acceptances. First the leading 
+//  particle and then the jet have to fullfill several conditions 
+//  (energy, direction ..) to be accepted. Then the fragmentation function 
+//  of this jet is constructed   
 // 
 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN) 
 //////////////////////////////////////////////////////////////////////////////
 
 // --- ROOT system ---
 
-#include "TString.h"
-#include "TFile.h"
-#include "TLorentzVector.h"
-#include "TList.h"
-#include "TTree.h"
-#include "TNtuple.h"
-#include "TParticle.h"
-#include "TCanvas.h"
-#include "TPaveLabel.h"
-#include "TPad.h"
-#include "TArrayC.h"
-#include "AliRun.h"
+#include <TFile.h>
+#include <TParticle.h>
+#include <TCanvas.h>
+#include <TPaveLabel.h>
+#include <TPad.h>
+#include <TH2.h>
+
 #include "AliPHOSGammaJet.h" 
 #include "AliPHOSGetter.h" 
 #include "AliPHOSGeometry.h"
-#include "TF1.h"
+#include "AliPHOSFastGlobalReconstruction.h"
+#include "AliESD.h"
+#include "AliESDtrack.h"
 #include "Riostream.h"
-#include "../PYTHIA6/AliGenPythia.h"
+
 
 ClassImp(AliPHOSGammaJet)
 
@@ -51,9 +70,10 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
   // ctor
   fAngleMaxParam.Set(4) ;
   fAngleMaxParam.Reset(0.);
-  fAnyConeOrPt      = kFALSE ;
+  fAnyConeOrPt      =  ;
   fEtaCut           = 0. ;
-  fFastRec          = 0 ;
+  fESDdata          = 0  ;
+  fFastRec          = 0  ;
   fInvMassMaxCut    = 0. ;
   fInvMassMinCut    = 0. ;
   fJetRatioMaxCut   = 0. ;
@@ -61,11 +81,11 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
   fJetTPCRatioMaxCut   = 0. ;
   fJetTPCRatioMinCut   = 0. ;
   fMinDistance      = 0. ;
-  fNEvent           = 0 ;
-  fNCone            = 0 ;
-  fNPt              = 0 ;
-  fOnlyCharged      = kFALSE ;
-  fOptFast          = kFALSE ;
+  fNEvent           = 0  ;
+  fNCone            = 0  ;
+  fNPt              = 0  ;
+  fOnlyCharged      =  ;
+  fOptFast          =  ;
   fPhiMaxCut        = 0. ;
   fPhiMinCut        = 0. ;
   fPtCut            = 0. ;
@@ -75,8 +95,13 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
   fRatioMinCut      = 0. ;
   fCone             = 0  ;
   fPtThreshold      = 0  ;
-  fTPCCutsLikeEMCAL = kFALSE ;
-  fSelect           = kFALSE ;
+  fTPCCutsLikeEMCAL = 0  ;
+  fSelect           = 0  ;
+
+  fDirName          = "" ;
+  fESDTree          = "" ;
+  fPattern          = "" ;
+
   for(Int_t i = 0; i<10; i++){
     fCones[i]         = 0.0 ;
     fNameCones[i]     = ""  ;
@@ -100,7 +125,7 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
     }
   }
 
-  fOption         = "" ;
+  fOptionGJ       = "" ;
   fOutputFile     = new TFile(gDirectory->GetName()) ;
   fInputFileName  = gDirectory->GetName() ;
   fOutputFileName = gDirectory->GetName() ;
@@ -130,7 +155,7 @@ AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) : 
   TTask("GammaJet","Analysis of gamma-jet correlations")
 {
-  // ctor 
+  // ctor
   fInputFileName = inputfilename;
   fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);  
   AliPHOSGetter *  gime = AliPHOSGetter::Instance(fInputFileName) ;
@@ -145,11 +170,12 @@ AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj)
   // cpy ctor
   fAngleMaxParam     = gj.fAngleMaxParam;
   fAnyConeOrPt       = gj.fAnyConeOrPt;
+  fESDdata           = gj.fESDdata;
   fEtaCut            = gj.fEtaCut ;
   fInvMassMaxCut     = gj.fInvMassMaxCut ;
   fInvMassMinCut     = gj.fInvMassMinCut ;
   fFastRec           = gj.fFastRec ;
-  fOption            = gj.fOption ;
+  fOptionGJ          = gj.fOptionGJ ;
   fMinDistance       = gj.fMinDistance ;
   fOptFast           = gj.fOptFast ;
   fOnlyCharged       = gj.fOnlyCharged ;
@@ -183,6 +209,10 @@ AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj)
   fCone              = gj.fCone ;
   fPtThreshold       = gj.fPtThreshold  ;
 
+  fDirName           = gj.fDirName ;
+  fESDTree           = gj.fESDTree ;
+  fPattern           = gj.fPattern ;
+
   SetName (gj.GetName()) ; 
   SetTitle(gj.GetTitle()) ; 
 
@@ -220,8 +250,7 @@ AliPHOSGammaJet::~AliPHOSGammaJet()
 void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList, 
                                      TClonesArray * plCh, 
                                      TClonesArray * plNe,  
-                                     TClonesArray * plNePHOS,
-                                     const AliPHOSGeometry * geom )
+                                     TClonesArray * plNePHOS)
 {
 
   // List of particles copied to a root file.
@@ -242,11 +271,11 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
   Char_t fi[100];
   sprintf(fi,"bgrd/%d/list.root",iEvent);
   TFile * f = TFile::Open(fi) ;
-  cout<<f->GetName()<<endl;
+  //cout<<f->GetName()<<endl;
 
   TParticle *particle = new TParticle();
-  TTree *T = (TTree*) f->Get("T");
-  TBranch *branch = T->GetBranch("primaries");
+  TTree *t = (TTree*) f->Get("T");
+  TBranch *branch = t->GetBranch("primaries");
   branch->SetAddress(&particle);
 
   Int_t index   = particleList->GetEntries() ; 
@@ -257,10 +286,14 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
   Int_t iParticle = 0 ; 
   Int_t m = 0;
   Double_t x = 0., z = 0.;
-  //  cout<<"bkg entries "<<T->GetEntries()<<endl;
+  //  cout<<"bkg entries "<<t->GetEntries()<<endl;
+
+  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
+
   if(!fOptFast){
-    for (iParticle=0 ; iParticle < T->GetEntries() ; iParticle++) {
-      T->GetEvent(iParticle) ;
+    for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
+      t->GetEvent(iParticle) ;
  
       m = 0 ;
       x = 0. ;
@@ -280,7 +313,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
          (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
          indexCh++ ;
          
-         if(strstr(fOption,"deb all")||strstr(fOption,"deb")){
+         if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){
            dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
              Fill(particle->Pt());
          }
@@ -290,7 +323,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
         
          if(m != 0)
            {//Is in PHOS
-             if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+             if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                  Fill(particle->Pt());
        
@@ -302,7 +335,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
          if((particle->Phi()>fPhiEMCALCut[0]) && 
             (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
            {//Is in EMCAL 
-             if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+             if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                  Fill(particle->Pt());
 
@@ -328,8 +361,8 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
     //     fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName);
     //     fFastRec->FastReconstruction(iiEvent);
     
-    for (iParticle=0 ; iParticle <   T->GetEntries() ; iParticle++) {
-      T->GetEvent(iParticle) ;
+    for (iParticle=0 ; iParticle <   t->GetEntries() ; iParticle++) {
+      t->GetEvent(iParticle) ;
       m = 0 ;
       x = 0. ;
       z = 0. ;
@@ -397,7 +430,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
            //Check if decay photons are too close for PHOS
            cellDistance = angle*460; //cm
            if (cellDistance < fMinDistance) {
-             if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+             if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                  Fill(particle->Pt());
              
@@ -423,7 +456,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
                  indexNe++ ; 
                }
                if(m != 0){
-                 if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+                 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                    dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                      Fill(particle->Pt());
                  new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
@@ -452,7 +485,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
                  geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
                  if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
                      && m == 0){
-                   if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                   if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                      dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                        Fill(photon1->Pt());
                    new((*particleList)[index]) TParticle(*photon1) ;
@@ -465,7 +498,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
                    p1 = kTRUE;
                  }
                  if(m != 0){
-                   if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                   if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                      dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                        Fill(photon1->Pt());
                    new((*plNePHOS)[indexNePHOS])       TParticle(*photon1) ;
@@ -486,7 +519,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
                  geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
                  if(photon2->Phi()>fPhiEMCALCut[0] && 
                     photon2->Phi()<fPhiEMCALCut[1] && m == 0){
-                   if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                   if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                      dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                        Fill(photon2->Pt());
                    new((*particleList)[index]) TParticle(*photon2) ;
@@ -498,7 +531,7 @@ void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
                    indexNe++ ; 
                  }
                  if(m != 0){
-                   if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                   if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                      dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                        Fill(photon2->Pt());
                    new((*plNePHOS)[indexNePHOS])       TParticle(*photon2) ;
@@ -529,14 +562,14 @@ Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
                                                 const Double_t *x) {
 
   //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]);
-  Double_t Epp = par[0] + par[1] * ptg ;
-  Double_t Spp = par[2] + par[3] * ptg ;
+  Double_t epp = par[0] + par[1] * ptg ;
+  Double_t spp = par[2] + par[3] * ptg ;
   Double_t f   = x[0]   + x[1]   * ptg ;
-  Double_t Epb = Epp + par[4] ;
-  Double_t Spb = TMath::Sqrt(Spp*Spp+ par[5]*par[5]) ;
-  Double_t rat = (Epb - Spb * f) / ptg ;
-  //Info("CalculateLimit","Epp %f, Spp %f, f %f", Epp, Spp, f);
-  //Info("CalculateLimit","Epb %f, Spb %f, rat %f", Epb, Spb, rat);
+  Double_t epb = epp + par[4] ;
+  Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
+  Double_t rat = (epb - spb * f) / ptg ;
+  //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
+  //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
   return rat ;
 }
 
@@ -545,13 +578,14 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
                                         TClonesArray * particleList, 
                                         TClonesArray * plCh, 
                                         TClonesArray * plNe,  
-                                        TClonesArray * plNePHOS,
-                                        const AliPHOSGeometry * geom )
+                                        TClonesArray * plNePHOS)
 {
   //Info("CreateParticleList","Inside");
-  AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ; 
+  AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
+  const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
   gime->Event(iEvent, "X") ;
 
+
   Int_t index = particleList->GetEntries() ; 
   Int_t indexCh     = plCh->GetEntries() ;
   Int_t indexNe     = plNe->GetEntries() ;
@@ -579,7 +613,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
            ->GetParticle(particle->GetPdgCode())->Charge();
          if((charge != 0) && (particle->Pt() > fChargedPtCut)){
            
-           if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+           if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
              dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                Fill(particle->Pt());
            new((*plCh)[indexCh++])       TParticle(*particle) ;
@@ -589,7 +623,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
            geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
            if(m != 0)
              {//Is in PHOS
-               if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+               if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                  dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                    Fill(particle->Pt());
                
@@ -598,7 +632,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
            if((particle->Phi()>fPhiEMCALCut[0]) && 
               (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
              {//Is in EMCAL 
-               if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+               if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                  dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                    Fill(particle->Pt());
                new((*plNe)[indexNe++])       TParticle(*particle) ;
@@ -693,7 +727,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
                    if(part.Pt() > fNeutralPtCut){
                      if(particle->Phi()>fPhiEMCALCut[0] && 
                         particle->Phi()<fPhiEMCALCut[1] && m == 0){
-                       if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+                       if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                          dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                            Fill(particle->Pt());
                        
@@ -707,7 +741,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
                        indexNe++;
                      }//InEMCAL
                      if(m != 0){
-                       if(strstr(fOption,"deb all")|| strstr(fOption,"deb"))
+                       if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
                          dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                            Fill(particle->Pt());
                        new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
@@ -734,7 +768,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
                      geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
                      if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
                          && m == 0){
-                     if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                     if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                        dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                          Fill(photon1->Pt());
                      new((*plNe)[indexNe++])       TParticle(*photon1) ;
@@ -743,7 +777,7 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
                      p1 = kTRUE;
                      }
                      if(m != 0){
-                       if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                       if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                          dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                            Fill(photon1->Pt());
                        new((*plNePHOS)[indexNePHOS++])       TParticle(*photon1) ;
@@ -761,14 +795,14 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
                      geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
                      if(photon2->Phi()>fPhiEMCALCut[0] && 
                         photon2->Phi()<fPhiEMCALCut[1] && m == 0){
-                       if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                       if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                          dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                            Fill(photon2->Pt());
                        new((*plNe)[indexNe++])       TParticle(*photon2) ;
                        new((*particleList)[index++]) TParticle(*photon2) ;
                      }
                      if(m != 0){
-                       if(strstr(fOption,"deb all") || strstr(fOption,"deb"))
+                       if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
                          dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
                            Fill(photon2->Pt());
                        new((*plNePHOS)[indexNePHOS++])       TParticle(*photon2) ;
@@ -791,6 +825,134 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
     }//OptFast
   //gime->Delete() ; 
 }
+//____________________________________________________________________________
+void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl, 
+                                               TClonesArray * plCh, 
+                                               TClonesArray * plNe,  
+                                               TClonesArray * plNePHOS,
+                                               const AliESD * esd){
+  
+  //Create a list of particles from the ESD. These particles have been measured 
+  //by the Central Tracking system (TPC), PHOS and EMCAL 
+  //(EMCAL not available for the moment). 
+  //Info("CreateParticleListFromESD","Inside");
+
+  Int_t index = pl->GetEntries() ; 
+  Int_t npar  = 0 ;
+  Double_t pid[AliPID::kSPECIESN];  
+
+  //########### PHOS ##############
+  //Info("CreateParticleListFromESD","Fill ESD PHOS list");
+  Int_t begphos = esd->GetFirstPHOSParticle();  
+  Int_t endphos = esd->GetFirstPHOSParticle() + 
+    esd->GetNumberOfPHOSParticles() ;  
+  Int_t indexNePHOS = plNePHOS->GetEntries() ;
+  if(strstr(fOptionGJ,"deb all"))
+    Info("CreateParticleListFromESD","PHOS: first particle %d, last particle %d",
+        begphos,endphos);
+
+  for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
+    AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
+   
+    //Create a TParticle to fill the particle list
+
+    Double_t en = track->GetPHOSsignal() ;
+    Double_t * p = new Double_t();
+    track->GetPHOSposition(p) ;
+    TVector3 pos(p[0],p[1],p[2]) ; 
+    Double_t phi  = pos.Phi();
+    Double_t theta= pos.Theta();
+    Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
+    Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
+    Double_t pz = en*TMath::Cos(theta);
+
+    TParticle * particle = new TParticle() ;
+    particle->SetMomentum(px,py,pz,en) ;
+
+    //Select only photons
+    
+    track->GetPHOSpid(pid);
+    //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
+    if( pid[AliPID::kPhoton] > 0.75)
+      new((*plNePHOS)[indexNePHOS++])   TParticle(*particle) ;
+  }
+
+  //########### TPC #####################
+  //Info("CreateParticleListFromESD","Fill ESD TPC list");
+  Int_t begtpc = 0 ;  
+  Int_t endtpc = esd->GetNumberOfTracks() ;
+  Int_t indexCh     = plCh->GetEntries() ;
+  if(strstr(fOptionGJ,"deb all"))
+    Info("CreateParticleListFromESD","TPC: first particle %d, last particle %d",
+        begtpc,endtpc);
+
+  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
+    AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
+    
+    Double_t en = track ->GetTPCsignal() ;
+    TVector3 mom = track->P3() ;
+    Double_t px = mom.Px();
+    Double_t py = mom.Py();
+    Double_t pz = mom.Pz(); //Check with TPC people if this is correct.
+
+    //cout<<"TPC signal "<<en<<endl;
+    //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
+    TParticle * particle = new TParticle() ;
+    particle->SetMomentum(px,py,pz,en) ;
+
+    new((*plCh)[indexCh++])       TParticle(*particle) ;    
+    new((*pl)[index++])           TParticle(*particle) ;
+
+  }
+
+  //################ EMCAL ##############
+  Double_t v[3] ; //vertex ;
+  esd->GetVertex()->GetXYZ(v) ; 
+  //##########Uncomment when ESD for EMCAL works ##########  
+  //Info("CreateParticleListFromESD","Fill ESD EMCAL list");
+  Int_t begem = esd->GetFirstEMCALParticle();  
+  Int_t endem = esd->GetFirstEMCALParticle() + 
+    esd->GetNumberOfEMCALParticles() ;  
+  Int_t indexNe  = plNe->GetEntries() ; 
+  if(strstr(fOptionGJ,"deb all"))
+    Info("CreateParticleListFromESD","EMCAL: first particle %d, last particle %d",
+        begem,endem);
+   
+  for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
+     AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
+  
+    Double_t en = track->GetEMCALsignal() ;
+    Double_t *p = new Double_t();
+    track->GetEMCALposition(p) ;
+    TVector3 pos(p[0],p[1],p[2]) ;
+    Double_t phi  = pos.Phi();
+    Double_t theta= pos.Theta();
+    Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
+    Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
+    Double_t pz = en*TMath::Cos(theta);
+    //cout<<"EMCAL signal "<<en<<endl;
+    //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
+    //TParticle * particle = new TParticle() ;
+    //particle->SetMomentum(px,py,pz,en) ;
+  
+    Int_t pdg = 0;
+    //  //Uncomment if PID IS WORKING, photon and pi0 idenitification.
+    //  //if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
+    //  //pdg = 22;
+    //  //else if( pid[AliPID::kPi0] > 0.75)
+    //  //pdg = 111;
+    pdg = 22; //No PID, assume all photons
+    TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                        px, py, pz, en, v[0], v[1], v[2], 0);
+
+    new((*plNe)[indexNe++])       TParticle(*particle) ; 
+    new((*pl)[index++])           TParticle(*particle) ;
+  }
+  
+  //  Info("CreateParticleListFromESD","End Inside");
+  
+}
 
 
 
@@ -798,12 +960,26 @@ void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
 void AliPHOSGammaJet::Exec(Option_t *option) 
 {
   // does the job
-  fOption = option;
+  fOptionGJ = option;
   MakeHistos() ; 
+  
+  AliESD * esd = 0;
+  TChain * t = 0 ;
+
+  if(fESDdata){
+    // Create chain of esd trees
+    const UInt_t kNevent = static_cast<UInt_t>(GetNEvent()) ; 
+    t = ReadESD(kNevent, fDirName, fESDTree, fPattern) ; 
+    if(!t) {
+      AliError("Could not create the TChain") ; 
+      //break ;
+    }
+    
+    // ESD  
+    t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree  
 
-  AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
-  const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
-
+  }
 
 //   AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator();
 //   pyth->Init();
@@ -814,7 +990,7 @@ void AliPHOSGammaJet::Exec(Option_t *option)
   TClonesArray * plNePHOS     = new TClonesArray("TParticle",1000);
 
   for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) {
-    if(strstr(fOption,"deb")||strstr(fOption,"deb all"))
+    if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
       Info("Exec", "Event %d", iEvent) ;
 
     fRan.SetSeed(0);
@@ -825,45 +1001,56 @@ void AliPHOSGammaJet::Exec(Option_t *option)
  
     TLorentzVector jet   (0,0,0,0);
     TLorentzVector jettpc(0,0,0,0);
-    
-    CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS, geom ); 
 
-//    TLorentzVector pyjet(0,0,0,0);
+    if(fESDdata){
 
-//     Int_t nJ, nJT;
-//     Float_t jets[4][10];
-//     pyth->SetJetReconstructionMode(1);
-//     pyth->LoadEvent();
-//     pyth->GetJets(nJ, nJT, jets);
-    
-//     Float_t pxJ = jets[0][0];
-//     Float_t pyJ = jets[1][0];
-//     Float_t pzJ = jets[2][0];
-//     Float_t eJ  = jets[3][0];
-//     pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
-    
-//     if(nJT > 1){
-//       //Info("Exec",">>>>>>>>>>Number of jets !!!!   %d",nJT);
-//       for (Int_t iJ = 1; iJ < nJT; iJ++) {
-//     Float_t pxJ = jets[0][iJ];
-//     Float_t pyJ = jets[1][iJ];
-//     Float_t pzJ = jets[2][iJ];
-//     Float_t eJ  = jets[3][iJ];
-//     pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
-//     //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
-//     //           iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
-//       }
+      Int_t iNbytes = t->GetEntry(iEvent); // store event in esd
+      //cout<<"nbytes "<<iNbytes<<endl;
+      if ( iNbytes == 0 ) {
+       AliError("Empty TChain") ; 
+       break ;
+      }      
+      CreateParticleListFromESD(particleList, plCh,plNe,plNePHOS, esd); //,iEvent);
+    }
+    else{   
+      CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS); 
+      
+      //    TLorentzVector pyjet(0,0,0,0);
+      
+      //     Int_t nJ, nJT;
+      //     Float_t jets[4][10];
+      //     pyth->SetJetReconstructionMode(1);
+      //     pyth->LoadEvent();
+      //     pyth->GetJets(nJ, nJT, jets);
+      
+      //     Float_t pxJ = jets[0][0];
+      //     Float_t pyJ = jets[1][0];
+      //     Float_t pzJ = jets[2][0];
+      //     Float_t eJ  = jets[3][0];
+      //     pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
       
-//     }
+      //     if(nJT > 1){
+      //       //Info("Exec",">>>>>>>>>>Number of jets !!!!   %d",nJT);
+      //       for (Int_t iJ = 1; iJ < nJT; iJ++) {
+      //       Float_t pxJ = jets[0][iJ];
+      //       Float_t pyJ = jets[1][iJ];
+      //       Float_t pzJ = jets[2][iJ];
+      //       Float_t eJ  = jets[3][iJ];
+      //       pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
+      //       //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
+      //       //           iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
+      //       }
+      
+      //     }
+      
+      if(fHIJING)
+       AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS);
+    }
     
-    if(fHIJING)
-      AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS, geom);
-
-
-    Bool_t IsInPHOS = kFALSE ;
-    GetGammaJet(plNePHOS, ptg, phig, etag, IsInPHOS) ; 
+    Bool_t iIsInPHOS = kFALSE ;
+    GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ; 
 
-    if(IsInPHOS){
+    if(iIsInPHOS){
 
       //Info("Exec"," In PHOS") ;
       dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
@@ -871,7 +1058,7 @@ void AliPHOSGammaJet::Exec(Option_t *option)
        ->Fill(ptg,phig);
       dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
        ->Fill(ptg,etag);
-      if(strstr(fOption,"deb")||strstr(fOption,"deb all"))
+      if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
        Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
             phig,etag) ;
 
@@ -912,7 +1099,7 @@ void AliPHOSGammaJet::Exec(Option_t *option)
            ->Fill(ptg,phig-phich);
          dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
            ->Fill(ptg,etag-etach);
-         if(strstr(fOption,"deb"))
+         if(strstr(fOptionGJ,"deb"))
          Info("Exec"," Charged Leading") ;
        }
        if((ptpi > ptch) && insidepi){
@@ -927,18 +1114,18 @@ void AliPHOSGammaJet::Exec(Option_t *option)
          dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
            ->Fill(ptg,etag-etapi);
          
-         if(ptpi > 0. && strstr(fOption,"deb"))
+         if(ptpi > 0. && strstr(fOptionGJ,"deb"))
            Info("Exec"," Pi0 Leading") ;
        }
                
-       if(strstr(fOption,"deb"))
+       if(strstr(fOptionGJ,"deb"))
          Info("Exec","Leading pt %f, phi %f",ptl,phil);
        if(insidech || insidepi){
          if(!fAnyConeOrPt){
           
            MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet);
 
-           if(strstr(fOption,"deb")){
+           if(strstr(fOptionGJ,"deb")){
 //           Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
 //                pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
              Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f",
@@ -958,7 +1145,7 @@ void AliPHOSGammaJet::Exec(Option_t *option)
        //TPC
        if(fOnlyCharged && ptch > 0.)
          {
-           if(strstr(fOption,"deb"))
+           if(strstr(fOptionGJ,"deb"))
              Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
 
            dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
@@ -972,7 +1159,7 @@ void AliPHOSGammaJet::Exec(Option_t *option)
           
              MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
 
-             if(strstr(fOption,"deb")){
+             if(strstr(fOptionGJ,"deb")){
 //             Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
 //                  pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
                Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f",
@@ -1011,7 +1198,8 @@ void AliPHOSGammaJet::Exec(Option_t *option)
 void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, 
                                    TString conf, TString type)
 {
-  
+  //Fill jet fragmentation histograms if !fAnyCone, 
+  //only for fCone and fPtThres 
   TParticle * particle = 0 ;
   Int_t ipr = -1 ;
   Float_t  charge = 0;
@@ -1042,7 +1230,8 @@ void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
                                               TString conf, TString type,
                                               TString cone, TString ptcut)
 {
-  
+   //Fill jet fragmentation histograms if fAnyCone, 
+   //for several cones and pt thresholds  
    TParticle *particle = 0;
    Int_t ipr=-1;
    Float_t  charge = 0;
@@ -1072,9 +1261,9 @@ void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
 
 //____________________________________________________________________________
 void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt, 
-                                 Double_t &phi, Double_t &eta, Bool_t &Is)
-
+                                 Double_t &phi, Double_t &eta, Bool_t &Is) const 
 {
+  //Search for the prompt photon in PHOS with pt > fPtCut
   pt  = -10.;
   eta = -10.;
   phi = -10.;
@@ -1082,15 +1271,12 @@ void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
 
-    if( (particle->Pt() > fPtCut ) && ( particle->GetStatusCode() == 1) 
-       && ( particle->GetPdgCode() == 22 || particle->GetPdgCode() == 111) ){
-       
-      if (particle->Pt() > pt) {
-       pt  = particle->Pt();          
-       phi = particle->Phi() ;
-       eta = particle->Eta() ;
-       Is  = kTRUE;
-      }
+    if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){
+
+      pt  = particle->Pt();          
+      phi = particle->Phi() ;
+      eta = particle->Eta() ;
+      Is  = kTRUE;
     }
   }
 }
@@ -1098,9 +1284,10 @@ void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
 //____________________________________________________________________________
 void  AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl, 
                                        Double_t ptg, Double_t phig, 
-                                       Double_t &pt, Double_t &eta, Double_t &phi)
+                                       Double_t &pt, Double_t &eta, Double_t &phi) const 
 {  
-
+  //Search for the charged particle with highest with 
+  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
   pt  = -100.;
   eta = -100;
   phi = -100;
@@ -1129,8 +1316,11 @@ void  AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
 //____________________________________________________________________________
 void  AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl, 
                                     Double_t ptg, Double_t phig, 
-                                    Double_t &pt,  Double_t &eta, Double_t &phi)
+                                    Double_t &pt,  Double_t &eta, Double_t &phi)  
 {  
+
+  //Search for the neutral pion with highest with 
+  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
   pt  = -100.;
   eta = -100.;
   phi = -100.;
@@ -1209,8 +1399,7 @@ void  AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
 
     while ( (particle = (TParticle*)next()) ) {
       iPrimary++;        
-      //     if(particle->GetStatusCode() == 1){
-
+    
       ksPdg = particle->GetPdgCode();
       ptl  = particle->Pt();
       if(ksPdg == 111){ //2 gamma
@@ -1314,6 +1503,9 @@ void  AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
 //____________________________________________________________________________
 void AliPHOSGammaJet::InitParameters()
 {
+
+  //Initialize the parameters of the analysis.
+
   fAngleMaxParam.Set(4) ;
   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
   fAngleMaxParam.AddAt(-0.25,1) ;
@@ -1321,7 +1513,7 @@ void AliPHOSGammaJet::InitParameters()
   fAngleMaxParam.AddAt(-2e-4,3) ;
   fAnyConeOrPt    = kFALSE ;
   fOutputFileName = "GammaJet.root" ;
-  fOption         = "";
+  fOptionGJ         = "";
   fHIJINGFileName = "galice.root" ;
   fHIJING         = kFALSE ;
   fMinDistance    = 3.6 ;
@@ -1330,6 +1522,7 @@ void AliPHOSGammaJet::InitParameters()
   fInvMassMinCut  = 0.12 ;
   fOnlyCharged    = kFALSE ;
   fOptFast        = kFALSE ;
+  fESDdata        = kTRUE ;
   fPhiEMCALCut[0] = 60. *TMath::Pi()/180.;
   fPhiEMCALCut[1] = 180.*TMath::Pi()/180.;
   fPhiMaxCut      = 3.4 ;
@@ -1348,6 +1541,10 @@ void AliPHOSGammaJet::InitParameters()
   fJetTPCRatioMinCut = 0.3 ;
   fSelect         = kFALSE  ;
 
+  fDirName      = "./" ;
+  fESDTree      = "esdTree" ;
+  fPattern      = "." ;
+
   //Cut depending on gamma energy
 
   fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed
@@ -1408,7 +1605,9 @@ void AliPHOSGammaJet::InitParameters()
 }
 
 //__________________________________________________________________________-
-Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e){
+Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
+  //Check if the opening angle of the candidate pairs is inside 
+  //our selection windowd
   Bool_t result = kFALSE;
   Double_t mpi0 = 0.1349766;
   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
@@ -1427,6 +1626,7 @@ Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e){
 //__________________________________________________________________________-
 Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj, 
                                      const TString type ){
+  //Check if the energy of the reconstructed jet is within an energy window
 
   Double_t par[6];
   Double_t xmax[2];
@@ -2184,7 +2384,9 @@ void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
                              Double_t ptl, Double_t  phil, Double_t  etal,
                              TString conf, TLorentzVector & jet)
 {
-  
+  //Fill the jet with the particles around the leading particle with 
+  //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and 
+  //check if we select it. Fill jet histograms
   Float_t ptcut = 0. ;
   if(fHIJING){
     if(ptg > fPtJetSelectionCut)  ptcut = 2. ;
@@ -2239,7 +2441,7 @@ void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
   ptjet = jet.Pt();
   ptbkg = bkg.Pt();
 
-  if(strstr(fOption,"deb") || strstr(fOption,"deb all"))
+  if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
     Info("MakeJet","Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg);
  
 
@@ -2259,7 +2461,7 @@ void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
 
 
   if(IsJetSelected(ptg,ptjet,conf) || fSelect){
-    if(strstr(fOption,"deb") || strstr(fOption,"deb all"))
+    if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
       Info("MakeJet","JetSelected");
     dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
       Fill(ptg);
@@ -2278,6 +2480,10 @@ void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
                                         Double_t  phig, Double_t ptl, 
                                         Double_t  phil, Double_t  etal, 
                                         TString conf){
+
+  //Fill the jet with the particles around the leading particle with 
+  //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and 
+  //check if we select it. Fill jet i histograms
   
   TClonesArray * jetList = new TClonesArray("TParticle",1000);
   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
@@ -2337,7 +2543,7 @@ void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
       //Fill histograms
       if(ptjet > 0.) {
 
-       if(strstr(fOption,"deb")){
+       if(strstr(fOptionGJ,"deb")){
          Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]);
          Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg  %f",ptg,ptjet,ptbkg);
        }
@@ -2387,7 +2593,7 @@ void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
 //____________________________________________________________________________
 void  AliPHOSGammaJet::MakePhoton(TLorentzVector & particle)
 {
-  
+  //Fast reconstruction for photons
   Double_t energy = particle.E()  ;
   Double_t modenergy = MakeEnergy(energy) ;
   //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy);
@@ -2475,6 +2681,7 @@ void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
 //____________________________________________________________________________
 void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
 {
+  //Plot some relevant histograms of the analysis
   TH2F * h = dynamic_cast<TH2F*>(fOutputFile->Get(what));
   if(h){
     h->Draw();
@@ -2520,12 +2727,13 @@ void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
 }
 
 //____________________________________________________________________________
-void AliPHOSGammaJet::Print(char * opt) 
+void AliPHOSGammaJet::Print(const Option_t * opt) const
 {
+
+  //Print some relevant parameters set for the analysis
   if(! opt)
     return;
 
-  // Print 
   Info("Print", "%s %s", GetName(), GetTitle() ) ;
  
   printf("Eta cut           : %f\n", fEtaCut) ;
@@ -2543,11 +2751,73 @@ void AliPHOSGammaJet::Print(char * opt)
  
 } 
 
+//__________________________________________________________________________
+TChain * AliPHOSGammaJet::ReadESDfromdisk(const UInt_t eventsToRead, 
+                        const TString dirName, 
+                        const TString esdTreeName, 
+                        const char *  pattern) 
+{
+  // Reads ESDs from Disk
+ TChain *  rv = 0; 
+
+  AliInfo( Form("\nReading files in %s \nESD tree name is %s \nReading %d events", 
+               dirName.Data(), esdTreeName.Data(), eventsToRead) ) ;
+  
+  // create a TChain of all the files 
+  TChain * cESDTree = new TChain(esdTreeName) ; 
+
+  // read from the directory file until the require number of events are collected
+  void * from = gSystem->OpenDirectory(dirName) ;
+  if (!from) {
+    AliError( Form("Directory %s does not exist") ) ;
+    rv = 0 ;
+  }
+  else{ // reading file names from directory
+    const char * subdir ; 
+    // search all subdirectories witch matching pattern
+    while( (subdir = gSystem->GetDirEntry(from))  && 
+          (cESDTree->GetEntries() < eventsToRead)) {
+      if ( strstr(subdir, pattern) != 0 ) { 
+       char file[200] ; 
+        sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);    
+       AliInfo( Form("Adding %s\n", file) );
+       cESDTree->Add(file) ;
+      }
+    } // while file
+  
+    AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ;
+    rv = cESDTree ; 
+    
+  } // reading file names from directory
+  return rv ; 
+}
+
+//__________________________________________________________________________
+TChain * AliPHOSGammaJet::ReadESD(const UInt_t eventsToRead,
+                 const TString dirName, 
+                 const TString esdTreeName, 
+                 const char *  pattern)  
+{
+  // Read AliESDs files and return a Chain of events
+  if ( dirName == "" ) {
+    AliError("Give the name of the DIR where to find files") ; 
+    return 0 ; 
+  }
+  if ( esdTreeName == "" ) 
+    return ReadESDfromdisk(eventsToRead, dirName) ;
+  else if ( strcmp(pattern, "") == 0 )
+    return ReadESDfromdisk(eventsToRead, dirName, esdTreeName) ;
+  else 
+    return ReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;          
+}
+
 //___________________________________________________________________
 void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone, 
                             Double_t eta, Double_t phi)
 {
-  
+
+  //Check if the particle is inside the cone defined by the leading particle
   b = kFALSE;
   
   if(phi > TMath::TwoPi())