B?\008New version from Gustavo

diff --git a/PWG4/AliAnaGammaDirect.cxx b/PWG4/AliAnaGammaDirect.cxx
index 7c6ee3bf5b2e62c240b3056c5df45873c1a23162..7d0bd7550c7b05ea3a33a1b76b934f3e620d093c 100644 (file)
--- a/PWG4/AliAnaGammaDirect.cxx
+++ b/PWG4/AliAnaGammaDirect.cxx
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/23 17:17:29  schutz
+ * New Gamma package
+ *
  *
  */
 
@@ -41,7 +44,7 @@
 #include <TFile.h>
 #include <TParticle.h>
 #include <TH2.h>
-
+#include <TChain.h>
 #include "AliAnaGammaDirect.h" 
 #include "AliESD.h"
 #include "AliESDtrack.h"
@@ -58,7 +61,7 @@ ClassImp(AliAnaGammaDirect)
     fPrintInfo(0), fMinGammaPt(0.),
     fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
     fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
-    fNCones(0),fNPtThres(0),fMakeICMethod(0)
+    fMakeICMethod(0)
 {
   //Ctor        
   TList * list = gDirectory->GetListOfKeys() ; 
@@ -71,11 +74,6 @@ ClassImp(AliAnaGammaDirect)
     fOutputContainer->Add(h) ; 
   }
   
-  for(Int_t i = 0; i < 10 ; i++){
-    fConeSizes[i]=0;
-    fPtThresholds[i]=0;
-  }
-  
   // Input slot #0 works with an Ntuple
   DefineInput(0, TChain::Class());
   // Output slot #0 writes into a TH1 container
@@ -93,17 +91,11 @@ AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
   fConeSize(g.fConeSize),
   fPtThreshold(g.fPtThreshold),
   fPtSumThreshold(g.fPtSumThreshold), 
-  fNCones(g.fNCones),fNPtThres(g.fNPtThres),
   fMakeICMethod(g.fMakeICMethod)
 {
   // cpy ctor
   SetName (g.GetName()) ; 
   SetTitle(g.GetTitle()) ; 
-
-  for(Int_t i = 0; i < 10 ; i++){
-    fConeSizes[i]=  g.fConeSizes[i];
-    fPtThresholds[i]=   g.fPtThresholds[i];
-  }
 }
 
 //____________________________________________________________________________
@@ -115,9 +107,61 @@ AliAnaGammaDirect::~AliAnaGammaDirect()
   delete fhNGamma    ;  
   delete fhPhiGamma  ; 
   delete fhEtaGamma   ;  
-  delete fhPtCandidate ;
-  delete [] fhPtThresIsolated ;
-  delete [] fhPtSumIsolated ;
+
+}
+
+//______________________________________________________________________________
+void AliAnaGammaDirect::ConnectInputData(const Option_t*)
+{
+  // Initialisation of branch container and histograms 
+    
+  AliInfo(Form("*** Initialization of %s", GetName())) ; 
+  
+  // Get input data
+  fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
+  if (!fChain) {
+    AliError(Form("Input 0 for %s not found\n", GetName()));
+    return ;
+  }
+  
+  // One should first check if the branch address was taken by some other task
+  char ** address = (char **)GetBranchAddress(0, "ESD");
+  if (address) {
+    fESD = (AliESD*)(*address);
+  } else {
+    fESD = new AliESD();
+    SetBranchAddress(0, "ESD", &fESD);
+  }
+
+}
+
+//________________________________________________________________________
+void AliAnaGammaDirect::CreateOutputObjects()
+{  
+
+  // Init parameteres and create histograms to be saved in output file and 
+  // stores them in fOutputContainer
+  InitParameters();
+
+  fOutputContainer = new TObjArray(3) ;
+
+  //Histograms of highest gamma identified in Event
+  fhNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120); 
+  fhNGamma->SetYTitle("N");
+  fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
+  fOutputContainer->Add(fhNGamma) ; 
+  
+  fhPhiGamma  = new TH2F
+    ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
+  fhPhiGamma->SetYTitle("#phi");
+  fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhPhiGamma) ; 
+  
+  fhEtaGamma  = new TH2F
+    ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
+  fhEtaGamma->SetYTitle("#eta");
+  fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhEtaGamma) ;
 
 }
 
@@ -287,7 +331,7 @@ void AliAnaGammaDirect::Exec(Option_t *)
 {
   
   // Processing of one event
-    
+
   //Get ESDs
   Long64_t entry = fChain->GetReadEntry() ;
   
@@ -308,92 +352,53 @@ void AliAnaGammaDirect::Exec(Option_t *)
   TClonesArray * plEMCAL   = new TClonesArray("TParticle",1000);  // All particles measured in EMCAL as Gamma calorimeter
 
   TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
-  TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
-  
+ 
   //Fill lists with photons, neutral particles and charged particles
   //look for the highest energy photon in the event inside fCalorimeter
-
+  
   AliDebug(2, "Fill particle lists, get prompt gamma");
-
+  
   //Fill particle lists 
   CreateParticleList(particleList, plCTS,plEMCAL,plPHOS); 
-
+  
   if(fCalorimeter == "PHOS")
     plNe = plPHOS;
   if(fCalorimeter == "EMCAL")
     plNe = plEMCAL;
-
-
-  //_______________Analysis 1__________________________
-  //Look for the prompt photon in the selected calorimeter
-  if(fMakeICMethod < 3){
-    
-    Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
-    //Search highest energy prompt gamma in calorimeter
-    GetPromptGamma(plNe,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
-    
-    AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
-    
-    //If there is any photon candidate in fCalorimeter
-    if(iIsInPHOSorEMCAL){
-      if (fPrintInfo)
-       AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
-      
-    }//Gamma in Calo
-    
-  }//Analysis 1
-
-
-  //_______________Analysis 2__________________________
-  //Look for the prompt photon in the selected calorimeter
-  //Isolation Cut Analysis for both methods and different pt cuts and cones
-  if(fMakeICMethod == 3){
+  
+  
+  Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
+  //Search highest energy prompt gamma in calorimeter
+  GetPromptGamma(plNe,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
+  
+  AliDebug(1, Form("Is Gamma in %s? %d",fCalorimeter.Data(),iIsInPHOSorEMCAL));
+  
+  //If there is any photon candidate in fCalorimeter
+  if(iIsInPHOSorEMCAL){
+    if (fPrintInfo)
+      AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
     
-    for(Int_t ipr = 0; ipr < plNe->GetEntries() ; ipr ++ ){
-      TParticle * pCandidate = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
-      
-      if(pCandidate->Pt() > fMinGammaPt){
-       
-       Bool_t  icPtThres   = kFALSE;
-       Bool_t  icPtSum     = kFALSE;
-       
-       Float_t ptC             = pCandidate->Pt() ;
-       
-       fhPtCandidate->Fill(ptC);
-       
-       for(Int_t icone = 0; icone<fNCones; icone++){
-         fConeSize  = fConeSizes[icone] ;
-         Float_t coneptsum = 0 ;
-         for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
-           fPtThreshold   =  fPtThresholds[ipt] ;
-           MakeIsolationCut(plCTS,plNe, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
-           AliDebug(4,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
-                           pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
-           
-           fhPtThresIsolated[icone][ipt]->Fill(ptC); 
-         }//pt thresh loop
-         fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
-       }//cone size loop
-      }//min pt candidate
-    }//candidate loop
-  }//Analysis 2
+  }//Gamma in Calo
   
   AliDebug(2, "End of analysis, delete pointers");
   
   particleList->Delete() ; 
   plCTS->Delete() ;
   plNe->Delete() ;
+  plEMCAL->Delete() ;
   plPHOS->Delete() ;
-  pLeading->Delete();
   pGamma->Delete();
+ 
+  PostData(0, fOutputContainer);
 
-  delete plNe ;
+ delete plNe ;
   delete plCTS ;
+  //delete plPHOS ;
+  //delete plEMCAL ;
   delete particleList ;
-  //  delete pLeading;
+
   //  delete pGamma;
 
-  PostData(0, fOutputContainer);
 }    
 
 
@@ -443,40 +448,13 @@ void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS,
 }
 
   //____________________________________________________________________________
-void AliAnaGammaDirect::Init(const Option_t * )
+void AliAnaGammaDirect::InitParameters()
 {
-  // Initialisation of branch container 
  
-  AliDebug(2,Form("*** Initialization of %s", GetName())) ; 
-  
-  // Get input data
-  fChain = dynamic_cast<TChain *>(GetInputData(0)) ;
-  if (!fChain) {
-    AliError(Form("Input 0 for %s not found\n", GetName()));
-    return ;
-  }
-  
-  if (!fESD) {
-    // One should first check if the branch address was taken by some other task
-    char ** address = (char **)GetBranchAddress(0, "ESD") ;
-    if (address) 
-      fESD = (AliESD *)(*address) ; 
-    if (!fESD) 
-      fChain->SetBranchAddress("ESD", &fESD) ;  
-  }
-  // The output objects will be written to 
-  TDirectory * cdir = gDirectory ; 
-  // Open a file for output #0
-  char outputName[1024] ; 
-  sprintf(outputName, "%s.root", GetName() ) ; 
-  OpenFile(0, outputName , "RECREATE") ; 
-  if (cdir) 
-    cdir->cd() ; 
-
- //  //Initialize the parameters of the analysis.
+  //Initialize the parameters of the analysis.
   fCalorimeter="PHOS";
   fPrintInfo           = kTRUE;
-  fMinGammaPt  = 10. ;
+  fMinGammaPt  = 5. ;
 
   //Fill particle lists when PID is ok
   fEMCALPID = kFALSE;
@@ -486,15 +464,7 @@ void AliAnaGammaDirect::Init(const Option_t * )
   fPtThreshold         = 2.0; 
   fPtSumThreshold  = 1.; 
 
-  fNCones           = 4 ; 
-  fNPtThres         = 4 ; 
-  fConeSizes[0] = 0.1; fConeSizes[0] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
-  fPtThresholds[0]=1.; fPtThresholds[0]=2.; fPtThresholds[0]=3.; fPtThresholds[0]=4.;
-
-  fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method, 3 make isolation study
-
-  //Initialization of histograms 
-  MakeHistos() ;
+  fMakeICMethod = 1; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
 }
 
 //__________________________________________________________________
@@ -562,61 +532,6 @@ void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
 
 }
 
-//___________________________________________________________________
-void AliAnaGammaDirect::MakeHistos()
-{
-  // Create histograms to be saved in output file and 
-  // stores them in fOutputContainer
-  
-  fOutputContainer = new TObjArray(10000) ;
-
-  //Histograms of highest gamma identified in Event
-  fhNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120); 
-  fhNGamma->SetYTitle("N");
-  fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
-  fOutputContainer->Add(fhNGamma) ; 
-  
-  fhPhiGamma  = new TH2F
-    ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
-  fhPhiGamma->SetYTitle("#phi");
-  fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhPhiGamma) ; 
-  
-  fhEtaGamma  = new TH2F
-    ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
-  fhEtaGamma->SetYTitle("#eta");
-  fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhEtaGamma) ;
-
-  if( fMakeICMethod== 3 ){
-
-    //Isolation cut histograms
-    fhPtCandidate  = new TH1F
-      ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120); 
-    fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
-    fOutputContainer->Add(fhPtCandidate) ;
-    
-    char name[128];
-    char title[128];
-    for(Int_t icone = 0; icone<fNCones; icone++){
-      sprintf(name,"PtSumIsolated_Cone_%d",icone);
-      sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
-      fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
-      fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
-      fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
-      fOutputContainer->Add(fhPtSumIsolated[icone]) ; 
-      
-      for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
-       sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
-       sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
-       fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
-       fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
-       fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
-      }//icone loop
-  }//ipt loop
-  }
-}
-
 void AliAnaGammaDirect::Print(const Option_t * opt) const
 {
 
@@ -625,10 +540,11 @@ void AliAnaGammaDirect::Print(const Option_t * opt) const
     return;
 
   Info("Print", "%s %s", GetName(), GetTitle() ) ;
+  printf("IC method               =     %d\n", fMakeICMethod) ; 
   printf("Cone Size               =     %f\n", fConeSize) ; 
   printf("pT threshold           =     %f\n", fPtThreshold) ;
   printf("pT sum threshold   =     %f\n", fPtSumThreshold) ; 
-  printf("Min Gamma pT      =     %f\n", fMinGammaPt) ; 
+  printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ; 
   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ; 
 } 
 

diff --git a/PWG4/AliAnaGammaDirect.h b/PWG4/AliAnaGammaDirect.h
index 65a983ded341ebb5e37594e6a958700b0650351d..8afe1040c9162da0a544f22592af550bb6e0b9da 100644 (file)
--- a/PWG4/AliAnaGammaDirect.h
+++ b/PWG4/AliAnaGammaDirect.h
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/23 17:17:29  schutz
+ * New Gamma package
+ *
  *
  */
 
@@ -15,9 +18,10 @@
 // Class for the analysis of gamma  (gamma-jet, 
 // gamma-hadron(Arleo, TODO))
 // This class only contains 3 methods: one to fill lists of particles (ESDs) comming 
-//  from the CTS (ITS+TPC) and the calorimeters: The other search in the 
+//  from the CTS (ITS+TPC) and the calorimeters; the other search in the 
 //  corresponing calorimeter for the highest energy cluster, identify it as 
-//  prompt photon(Shower Shape (done) and Isolation Cut (TODO in new class?)).
+//  prompt photon(Shower Shape and Isolation Cut), the last method does the 
+//  isolation selection
 //
 //  Class created from old AliPHOSGammaJet
 //  (see AliRoot versions previous Release 4-09)
@@ -25,15 +29,11 @@
 //*-- Author: Gustavo Conesa (INFN-LNF)
 
 // --- ROOT system ---
-#include <TROOT.h>
-#include <TChain.h>
-#include "TTask.h"
-#include "TArrayD.h"
-#include "TChain.h"
-#include <TH2F.h>
-#include <TTree.h> 
 #include <TParticle.h> 
+#include <TClonesArray.h> 
+#include <TTree.h> 
 #include "AliAnalysisTask.h" 
+#include <TH2F.h>
 
 class AliESD ; 
  
@@ -44,10 +44,13 @@ public:
   AliAnaGammaDirect(const char *name) ; // default ctor
   AliAnaGammaDirect(const AliAnaGammaDirect & g) ; // cpy ctor
   virtual ~AliAnaGammaDirect() ; //virtual dtor
+
   virtual void Exec(Option_t * opt = "") ;
-  virtual void Init(Option_t * opt = "");
+  virtual void ConnectInputData(Option_t *);
+  virtual void CreateOutputObjects();
   virtual void Terminate(Option_t * opt = "");
   
+  void InitParameters();
   TTree *     GetChain()                const {return fChain ; }
   AliESD *    GetESD()                  const {return fESD ; }
   TObjArray * GetOutputContainer()      const {return fOutputContainer ; }
@@ -72,9 +75,6 @@ public:
   void SetPtSumThreshold(Float_t pt) {fPtSumThreshold = pt; };
   void SetICMethod(Int_t i )          {fMakeICMethod = i ; }
   
-  void SetConeSizes(Int_t i, Float_t r)              {fConeSizes[i] = r ; }
-  void SetPtThresholds(Int_t i, Float_t pt)        {fPtThresholds[i] = pt; };
-
   void SetEMCALPIDOn(Bool_t pid){ fEMCALPID= pid ; }
   void SetPHOSPIDOn(Bool_t pid){ fPHOSPID= pid ; }
 
@@ -105,22 +105,14 @@ public:
   Float_t      fConeSize ; //Size of the isolation cone 
   Float_t      fPtThreshold ; //Mimium pt of the particles in the cone to set isolation
   Float_t      fPtSumThreshold ; //Mimium pt sum of the particles in the cone to set isolation  
-  Int_t         fNCones   ; //Number of cone sizes to test
-  Int_t         fNPtThres ; //Number of ptThres to test
-  Float_t     fConeSizes[10] ; // Arrat with cones to test
-  Float_t     fPtThresholds[10] ; // Array with pt thresholds to test
   Int_t        fMakeICMethod ; //Isolation cut method to be used
                                            // 0: No isolation
                                            // 1: Pt threshold method
                                            // 2: Cone pt sum method
-                                           // 3: Study both methods for several cones and pt.
   //Histograms  
   TH1F * fhNGamma    ; 
   TH2F * fhPhiGamma    ; 
   TH2F * fhEtaGamma    ; 
-  TH1F * fhPtCandidate ;
-  TH1F* fhPtThresIsolated[10][10] ;
-  TH2F* fhPtSumIsolated[10] ;
 
   ClassDef(AliAnaGammaDirect,0)
 } ;

diff --git a/PWG4/AliAnaGammaHadron.cxx b/PWG4/AliAnaGammaHadron.cxx
index 8b0049d9fb9e649607cec132ede28ed9ce407303..b555ee2eef3267676f9526b032b765db08fb245c 100644 (file)
--- a/PWG4/AliAnaGammaHadron.cxx
+++ b/PWG4/AliAnaGammaHadron.cxx
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/23 17:17:29  schutz
+ * New Gamma package
+ *
  *
  */
 
@@ -127,14 +130,181 @@ AliAnaGammaHadron::~AliAnaGammaHadron()
 
 }
 
+//______________________________________________________________________________
+void AliAnaGammaHadron::ConnectInputData(const Option_t*)
+{
+  // Initialisation of branch container and histograms 
+  AliAnaGammaDirect::ConnectInputData("");  
+}
+
+//________________________________________________________________________
+void AliAnaGammaHadron::CreateOutputObjects()
+{  
+
+  // Init parameteres and create histograms to be saved in output file and 
+  // stores them in fOutputContainer
+  InitParameters();
+  AliAnaGammaDirect::CreateOutputObjects();
+
+  fOutputContainer = new TObjArray(100) ;
+  
+  //Use histograms in AliAnaGammaDirect
+  TObjArray  * outputContainer =GetOutputContainer();
+  for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
+    fOutputContainer->Add(outputContainer->At(i)) ;
+  
+  fhPhiCharged  = new TH2F
+    ("PhiCharged","#phi_{#pi^{#pm}}  vs p_{T #gamma}",
+     120,0,120,120,0,7); 
+  fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
+  fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhPhiCharged) ;
+  
+  fhPhiNeutral  = new TH2F
+    ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #gamma}",
+     120,0,120,120,0,7); 
+  fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
+  fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhPhiNeutral) ;  
+  
+  fhEtaCharged  = new TH2F
+    ("EtaCharged","#eta_{#pi^{#pm}}  vs p_{T #gamma}",
+     120,0,120,120,-1,1); 
+  fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
+  fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhEtaCharged) ;
+
+  fhEtaNeutral  = new TH2F
+    ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #gamma}",
+     120,0,120,120,-1,1); 
+  fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
+  fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhEtaNeutral) ;  
 
+  fhDeltaPhiGammaCharged  = new TH2F
+    ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
+     200,0,120,200,0,6.4); 
+  fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
+  fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhDeltaPhiGammaCharged) ; 
+  
+  fhDeltaEtaGammaCharged  = new TH2F
+    ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
+     200,0,120,200,-2,2); 
+  fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
+  fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhDeltaEtaGammaCharged) ; 
+
+  fhDeltaPhiGammaNeutral  = new TH2F
+    ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
+     200,0,120,200,0,6.4); 
+  fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
+  fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhDeltaPhiGammaNeutral) ; 
+  
+  fhDeltaEtaGammaNeutral  = new TH2F
+    ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
+     200,0,120,200,-2,2); 
+  fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
+  fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhDeltaEtaGammaNeutral) ; 
+  
+  //
+  fhAnglePairAccepted  = new TH2F
+    ("AnglePairAccepted",
+     "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
+     200,0,50,200,0,0.2); 
+  fhAnglePairAccepted->SetYTitle("Angle (rad)");
+  fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+  fOutputContainer->Add(fhAnglePairAccepted) ; 
+  
+  fhAnglePairNoCut  = new TH2F
+    ("AnglePairNoCut",
+     "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2); 
+  fhAnglePairNoCut->SetYTitle("Angle (rad)");
+  fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+  fOutputContainer->Add(fhAnglePairNoCut) ; 
+  
+  fhAnglePairAzimuthCut  = new TH2F
+    ("AnglePairAzimuthCut",
+     "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
+     200,0,50,200,0,0.2); 
+  fhAnglePairAzimuthCut->SetYTitle("Angle (rad)");
+  fhAnglePairAzimuthCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+  fOutputContainer->Add(fhAnglePairAzimuthCut) ; 
+  
+    fhAnglePairOpeningAngleCut  = new TH2F
+      ("AnglePairOpeningAngleCut",
+       "Angle between all #gamma pair (opening angle + azimuth cut) vs p_{T  #pi^{0}}"
+       ,200,0,50,200,0,0.2); 
+    fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
+    fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairOpeningAngleCut) ;
+    
+    fhAnglePairAllCut  = new TH2F
+      ("AnglePairAllCut",
+       "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs p_{T  #pi^{0}}"
+       ,200,0,50,200,0,0.2); 
+    fhAnglePairAllCut->SetYTitle("Angle (rad)");
+    fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairAllCut) ; 
+    
+    
+    //
+    fhInvMassPairNoCut  = new TH2F
+      ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairNoCut) ; 
+    
+    fhInvMassPairAzimuthCut  = new TH2F
+      ("InvMassPairAzimuthCut",
+       "Invariant Mass of #gamma pair (azimuth cuts) vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairAzimuthCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairAzimuthCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairAzimuthCut) ; 
+    
+    fhInvMassPairOpeningAngleCut  = new TH2F
+      ("InvMassPairOpeningAngleCut",
+       "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairOpeningAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
+    
+    fhInvMassPairAllCut  = new TH2F
+      ("InvMassPairAllCut",
+       "Invariant Mass of #gamma pair (opening angle+invmass cut+azimuth) vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairAllCut) ; 
+ 
+    //   
+    fhCorrelationGammaCharged  = 
+      new TH2F("CorrelationGammaCharged","z_{#gamma #pi} = p_{T #pi^{#pm}} / p_{T #gamma}",
+              240,0.,120.,1000,0.,1.2); 
+    fhCorrelationGammaCharged->SetYTitle("z_{#gamma #pi}");
+    fhCorrelationGammaCharged->SetXTitle("p_{T #gamma}");
+    fOutputContainer->Add(fhCorrelationGammaCharged) ;
+
+    fhCorrelationGammaNeutral  = 
+      new TH2F("CorrelationGammaNeutral","z_{#gamma #pi} = p_{T #pi^{0}} / p_{T #gamma}",
+              240,0.,120.,1000,0.,1.2); 
+    fhCorrelationGammaNeutral->SetYTitle("z_{#gamma #pi}");
+    fhCorrelationGammaNeutral->SetXTitle("p_{T #gamma}");
+    fOutputContainer->Add(fhCorrelationGammaNeutral) ;
+
+}
 
 //____________________________________________________________________________
 void AliAnaGammaHadron::Exec(Option_t *) 
 {
   
   // Processing of one event
-    
+
   //Get ESDs
   Long64_t entry = GetChain()->GetReadEntry() ;
   
@@ -362,10 +532,10 @@ void  AliAnaGammaHadron::MakeGammaNeutralCorrelation(TClonesArray * pl, TParticl
 }
 
   //____________________________________________________________________________
-void AliAnaGammaHadron::Init(const Option_t * )
+void AliAnaGammaHadron::InitParameters()
 {
   // Initialisation of branch container 
-  AliAnaGammaDirect::Init();
+  //AliAnaGammaDirect::InitParameters();
  
   //Initialize the parameters of the analysis.
   //fCalorimeter="PHOS";
@@ -387,10 +557,6 @@ void AliAnaGammaHadron::Init(const Option_t * )
   // fEMCALPID = kFALSE;
   // fPHOSPID = kFALSE;
 
-  //Initialization of histograms 
-
-  MakeHistos() ; 
-
 }
 
 //__________________________________________________________________________-
@@ -413,165 +579,6 @@ Bool_t AliAnaGammaHadron::IsAngleInWindow(const Float_t angle,const Float_t e) {
   return result;
 }
 
-//____________________________________________________________________________
-void AliAnaGammaHadron::MakeHistos()
-{
-  // Create histograms to be saved in output file and 
-  // stores them in fOutputContainer
-  
-  fOutputContainer = new TObjArray(10000) ;
-
-  //Use histograms in AliAnaGammaDirect
-  TObjArray  * outputContainer =GetOutputContainer();
-  for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
-    fOutputContainer->Add(outputContainer->At(i)) ;
-    
-  fhPhiCharged  = new TH2F
-    ("PhiCharged","#phi_{#pi^{#pm}}  vs p_{T #gamma}",
-     120,0,120,120,0,7); 
-  fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
-  fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhPhiCharged) ;
-  
-  fhPhiNeutral  = new TH2F
-    ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #gamma}",
-     120,0,120,120,0,7); 
-  fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
-  fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhPhiNeutral) ;  
-  
-  fhEtaCharged  = new TH2F
-    ("EtaCharged","#eta_{#pi^{#pm}}  vs p_{T #gamma}",
-     120,0,120,120,-1,1); 
-  fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
-  fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhEtaCharged) ;
-
-  fhEtaNeutral  = new TH2F
-    ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #gamma}",
-     120,0,120,120,-1,1); 
-  fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
-  fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhEtaNeutral) ;  
-
-  fhDeltaPhiGammaCharged  = new TH2F
-    ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
-     200,0,120,200,0,6.4); 
-  fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
-  fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhDeltaPhiGammaCharged) ; 
-  
-  fhDeltaEtaGammaCharged  = new TH2F
-    ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
-     200,0,120,200,-2,2); 
-  fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
-  fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhDeltaEtaGammaCharged) ; 
-
-  fhDeltaPhiGammaNeutral  = new TH2F
-    ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
-     200,0,120,200,0,6.4); 
-  fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
-  fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhDeltaPhiGammaNeutral) ; 
-  
-  fhDeltaEtaGammaNeutral  = new TH2F
-    ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
-     200,0,120,200,-2,2); 
-  fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
-  fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhDeltaEtaGammaNeutral) ; 
-  
-  //
-  fhAnglePairAccepted  = new TH2F
-    ("AnglePairAccepted",
-     "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
-     200,0,50,200,0,0.2); 
-  fhAnglePairAccepted->SetYTitle("Angle (rad)");
-  fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-  fOutputContainer->Add(fhAnglePairAccepted) ; 
-  
-  fhAnglePairNoCut  = new TH2F
-    ("AnglePairNoCut",
-     "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2); 
-  fhAnglePairNoCut->SetYTitle("Angle (rad)");
-  fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-  fOutputContainer->Add(fhAnglePairNoCut) ; 
-  
-  fhAnglePairAzimuthCut  = new TH2F
-    ("AnglePairAzimuthCut",
-     "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
-     200,0,50,200,0,0.2); 
-  fhAnglePairAzimuthCut->SetYTitle("Angle (rad)");
-  fhAnglePairAzimuthCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-  fOutputContainer->Add(fhAnglePairAzimuthCut) ; 
-  
-    fhAnglePairOpeningAngleCut  = new TH2F
-      ("AnglePairOpeningAngleCut",
-       "Angle between all #gamma pair (opening angle + azimuth cut) vs p_{T  #pi^{0}}"
-       ,200,0,50,200,0,0.2); 
-    fhAnglePairOpeningAngleCut->SetYTitle("Angle (rad)");
-    fhAnglePairOpeningAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairOpeningAngleCut) ;
-    
-    fhAnglePairAllCut  = new TH2F
-      ("AnglePairAllCut",
-       "Angle between all #gamma pair (opening angle + inv mass cut+azimuth) vs p_{T  #pi^{0}}"
-       ,200,0,50,200,0,0.2); 
-    fhAnglePairAllCut->SetYTitle("Angle (rad)");
-    fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairAllCut) ; 
-    
-    
-    //
-    fhInvMassPairNoCut  = new TH2F
-      ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairNoCut) ; 
-    
-    fhInvMassPairAzimuthCut  = new TH2F
-      ("InvMassPairAzimuthCut",
-       "Invariant Mass of #gamma pair (azimuth cuts) vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairAzimuthCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairAzimuthCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairAzimuthCut) ; 
-    
-    fhInvMassPairOpeningAngleCut  = new TH2F
-      ("InvMassPairOpeningAngleCut",
-       "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairOpeningAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairOpeningAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairOpeningAngleCut) ; 
-    
-    fhInvMassPairAllCut  = new TH2F
-      ("InvMassPairAllCut",
-       "Invariant Mass of #gamma pair (opening angle+invmass cut+azimuth) vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairAllCut) ; 
- 
-    //   
-    fhCorrelationGammaCharged  = 
-      new TH2F("CorrelationGammaCharged","z_{#gamma #pi} = p_{T #pi^{#pm}} / p_{T #gamma}",
-              240,0.,120.,1000,0.,1.2); 
-    fhCorrelationGammaCharged->SetYTitle("z_{#gamma #pi}");
-    fhCorrelationGammaCharged->SetXTitle("p_{T #gamma}");
-    fOutputContainer->Add(fhCorrelationGammaCharged) ;
-
-    fhCorrelationGammaNeutral  = 
-      new TH2F("CorrelationGammaNeutral","z_{#gamma #pi} = p_{T #pi^{0}} / p_{T #gamma}",
-              240,0.,120.,1000,0.,1.2); 
-    fhCorrelationGammaNeutral->SetYTitle("z_{#gamma #pi}");
-    fhCorrelationGammaNeutral->SetXTitle("p_{T #gamma}");
-    fOutputContainer->Add(fhCorrelationGammaNeutral) ;
-
-}
-
 //____________________________________________________________________________
 void AliAnaGammaHadron::Print(const Option_t * opt) const
 {

diff --git a/PWG4/AliAnaGammaHadron.h b/PWG4/AliAnaGammaHadron.h
index cf334715a893e9630b58d0fa894a0c036e9500aa..80f16eb9c880cf14211710ba93b163061fb018ac 100644 (file)
--- a/PWG4/AliAnaGammaHadron.h
+++ b/PWG4/AliAnaGammaHadron.h
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/23 17:17:29  schutz
+ * New Gamma package
+ *
  *
  */
 
@@ -40,9 +43,11 @@ public:
   AliAnaGammaHadron(const AliAnaGammaHadron & gj) ; // cpy ctor
   virtual ~AliAnaGammaHadron() ; //virtual dtor
   virtual void Exec(Option_t * opt = "") ;
-  virtual void Init(Option_t * opt = "");
+  virtual void ConnectInputData(Option_t *);
+  virtual void CreateOutputObjects();
   virtual void Terminate(Option_t * opt = "");
- 
+
+  void InitParameters(); 
   Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
   Double_t GetInvMassMaxCut() const {return fInvMassMaxCut ; }
   Double_t GetInvMassMinCut() const {return fInvMassMinCut ; }

diff --git a/PWG4/AliAnaGammaIsolCut.cxx b/PWG4/AliAnaGammaIsolCut.cxx
index 480248e874e58067666b988ea71f337c7d2611b8..4741b6f86d96c5a7ccdde38d663bf0d4dad03abc 100644 (file)
--- a/PWG4/AliAnaGammaIsolCut.cxx
+++ b/PWG4/AliAnaGammaIsolCut.cxx
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/25 17:24:20  schutz
+ * new class
+ *
  * Revision 1.1  2007/01/23 17:17:29  schutz
  * New Gamma package
  *
@@ -111,22 +114,63 @@ AliAnaGammaIsolCut::~AliAnaGammaIsolCut()
 
 }
 
+//______________________________________________________________________________
+void AliAnaGammaIsolCut::ConnectInputData(const Option_t*)
+{
+  // Initialisation of branch container and histograms 
+  AliAnaGammaDirect::ConnectInputData("");
+
+}
+
+//____________________________________________________
+void AliAnaGammaIsolCut::CreateOutputObjects()
+{  
+
+  // Init parameteres and create histograms to be saved in output file and 
+  // stores them in fOutputContainer
+  InitParameters();
+
+  fOutputContainer = new TObjArray(100) ; 
+
+  //Isolation cut histograms
+  fhPtCandidate  = new TH1F
+    ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120); 
+  fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
+  fOutputContainer->Add(fhPtCandidate) ;
+  
+  char name[128];
+  char title[128];
+  for(Int_t icone = 0; icone<fNCones; icone++){
+    sprintf(name,"PtSumIsolated_Cone_%d",icone);
+    sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
+    fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
+    fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
+    fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
+    fOutputContainer->Add(fhPtSumIsolated[icone]) ; 
+    
+    for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
+      sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
+      sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
+      fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
+      fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
+      fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
+    }//icone loop
+  }//ipt loop
 
+}
 
 //____________________________________________________________________________
 void AliAnaGammaIsolCut::Exec(Option_t *) 
 {
   
   // Processing of one event
-    
   //Get ESDs
+
   Long64_t entry = GetChain()->GetReadEntry() ;
-  
   if (!GetESD()) {
     AliError("fESD is not connected to the input!") ; 
     return ; 
   }
-  
   if (GetPrintInfo()) 
     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ; 
 
@@ -197,56 +241,17 @@ void AliAnaGammaIsolCut::Exec(Option_t *)
 }    
 
   //____________________________________________________________________________
-void AliAnaGammaIsolCut::Init(const Option_t * )
+void AliAnaGammaIsolCut::InitParameters()
 {
   // Initialisation of branch container 
-  AliAnaGammaDirect::Init();
+  AliAnaGammaDirect::InitParameters();
 
   fNCones           = 4 ; 
   fNPtThres         = 4 ; 
   fConeSizes[0] = 0.1; fConeSizes[0] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
   fPtThresholds[0]=1.; fPtThresholds[0]=2.; fPtThresholds[0]=3.; fPtThresholds[0]=4.;
-
-  //Initialization of histograms 
-  MakeHistos() ;
 }
 
-//___________________________________________________________________
-void AliAnaGammaIsolCut::MakeHistos()
-{
-  // Create histograms to be saved in output file and 
-  // stores them in fOutputContainer
-  
-  fOutputContainer = new TObjArray(10000) ; 
-
-  //Isolation cut histograms
-  fhPtCandidate  = new TH1F
-    ("PtCandidate","p_{T} of candidate particles for isolation",240,0,120); 
-  fhPtCandidate->SetXTitle("p_{T} (GeV/c)");
-  fOutputContainer->Add(fhPtCandidate) ;
-  
-  char name[128];
-  char title[128];
-  for(Int_t icone = 0; icone<fNCones; icone++){
-    sprintf(name,"PtSumIsolated_Cone_%d",icone);
-    sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
-    fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
-    fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
-    fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
-    fOutputContainer->Add(fhPtSumIsolated[icone]) ; 
-    
-    for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
-      sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
-      sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
-      fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
-      fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
-      fOutputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
-    }//icone loop
-  }//ipt loop
-
-}
-
-
 void AliAnaGammaIsolCut::Print(const Option_t * opt) const
 {
 

diff --git a/PWG4/AliAnaGammaIsolCut.h b/PWG4/AliAnaGammaIsolCut.h
index e6020f3f58f200d9dcc7e710aef05c053585d3d2..12a40ba1f0a151608b6c06c89b085acc705ab0c5 100644 (file)
--- a/PWG4/AliAnaGammaIsolCut.h
+++ b/PWG4/AliAnaGammaIsolCut.h
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/25 17:24:20  schutz
+ * new class
+ *
  * Revision 1.1  2007/01/23 17:17:29  schutz
  * New Gamma package
  *
@@ -43,9 +46,11 @@ public:
   AliAnaGammaIsolCut(const AliAnaGammaIsolCut & g) ; // cpy ctor
   virtual ~AliAnaGammaIsolCut() ; //virtual dtor
   virtual void Exec(Option_t * opt = "") ;
-  virtual void Init(Option_t * opt = "");
+  virtual void ConnectInputData(Option_t *);
+  virtual void CreateOutputObjects();
   virtual void Terminate(Option_t * opt = "");
 
+  void InitParameters();
   void Print(const Option_t * opt)const;
 
   Int_t    GetNCones()                  const {return fNCones ; }

diff --git a/PWG4/AliAnaGammaJet.cxx b/PWG4/AliAnaGammaJet.cxx
index 66a4c3a184b6a084ee4667e1dd105388bdeacae9..ef3932a7baa6ff8c147eb2717f4f611a821de647 100644 (file)
--- a/PWG4/AliAnaGammaJet.cxx
+++ b/PWG4/AliAnaGammaJet.cxx
@@ -17,6 +17,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/23 17:17:29  schutz
+ * New Gamma package
+ *
  *
  */
 
@@ -227,932 +230,939 @@ Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg,
   return rat ;
 }
 
-
-//____________________________________________________________________________
-void AliAnaGammaJet::Exec(Option_t *) 
-{
-  
-  // Processing of one event
-  
-  //Get ESDs
-  Long64_t entry = GetChain()->GetReadEntry() ;
-  
-  if (!GetESD()) {
-    AliError("fESD is not connected to the input!") ; 
-    return ; 
-  }
-  
-  if (GetPrintInfo()) 
-    AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ; 
-  
-  //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
-  
-  TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
-  TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
-  TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter
-  TClonesArray * plCalo     = new TClonesArray("TParticle",1000);  // All particles measured in Prompt Gamma calorimeter 
-  
-  
-  TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
-  TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
-  
-  Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
-  
-  //Fill lists with photons, neutral particles and charged particles
-  //look for the highest energy photon in the event inside fCalorimeter
-  AliDebug(2, "Fill particle lists, get prompt gamma");
-  
-  //Fill particle lists 
-  
-  if(GetCalorimeter() == "PHOS")
-    CreateParticleList(particleList, plCTS,plNe,plCalo); 
-  else if(GetCalorimeter() == "EMCAL")
-    CreateParticleList(particleList, plCTS,plCalo,plNe); 
-  else
-    AliError("No calorimeter selected");
-  
-  //Search highest energy prompt gamma in calorimeter
-  GetPromptGamma(plCalo,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
-  
-  AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
-  AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
-                 plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries()));
-
-  //If there is any prompt photon in fCalorimeter, 
-  //search jet leading particle
-  
-  if(iIsInPHOSorEMCAL){
-    if (GetPrintInfo())
-      AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
- 
-    AliDebug(2, "Get Leading Particles, Make jet");
-
-    //Search leading particles in CTS and EMCAL 
-    if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
-      if (GetPrintInfo())
-       AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
-
-      //Search Jet
-      if(!fSeveralConeAndPtCuts)
-       MakeJet(particleList,pGamma,pLeading,"");
-      else{
-       for(Int_t icone = 0; icone<fNCone; icone++) {
-         for(Int_t ipt = 0; ipt<fNPt;ipt++) {  
-           TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt];
-           MakeJet(particleList, pGamma, pLeading,lastname);
-         }//icone
-       }//ipt
-      }//fSeveralConeAndPtCuts
-    }//Leading
-  }//Gamma in Calo
-     
-  AliDebug(2, "End of analysis, delete pointers");
-
-  particleList->Delete() ; 
-  plCTS->Delete() ;
-  plNe->Delete() ;
-  plCalo->Delete() ;
-  pLeading->Delete();
-  pGamma->Delete();
-
-  delete plNe ;
-  delete plCalo ;
-  delete plCTS ;
-  delete particleList ;
-  //  delete pLeading;
-  //  delete pGamma;
-
-  PostData(0, fOutputContainer);
-}    
-
-//____________________________________________________________________________
-void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
+//______________________________________________________________________________
+void AliAnaGammaJet::ConnectInputData(const Option_t*)
 {
-  //Fill histograms wth jet fragmentation 
-  //and number of selected jets and leading particles
-  //and the background multiplicity
-  TParticle * particle = 0 ;
-  Int_t ipr = 0;
-  Float_t  charge = 0;
-
-  TIter next(pl) ; 
-  while ( (particle = dynamic_cast<TParticle*>(next())) ) {
-    ipr++ ;
-    Double_t pt = particle->Pt();
+  // Initialisation of branch container and histograms 
+  AliAnaGammaDirect::ConnectInputData("");    
 
-    charge = TDatabasePDG::Instance()
-      ->GetParticle(particle->GetPdgCode())->Charge();
-    if(charge != 0){//Only jet Charged particles 
-      dynamic_cast<TH2F*>
-       (fOutputContainer->FindObject(type+"Fragment"+lastname))
-       ->Fill(ptg,pt/ptg);
-      dynamic_cast<TH2F*>
-       (fOutputContainer->FindObject(type+"PtDist"+lastname))
-       ->Fill(ptg,pt);
-    }//charged
-
-  }//while
-
-  if(type == "Bkg")
-    dynamic_cast<TH1F*>
-      (fOutputContainer->FindObject("NBkg"+lastname))
-      ->Fill(ipr);
-  else{
-    dynamic_cast<TH1F*>
-      (fOutputContainer->FindObject("NJet"+lastname))->
-      Fill(ptg);
-    dynamic_cast<TH2F*>
-      (fOutputContainer->FindObject("NLeading"+lastname))
-      ->Fill(ptg,ptl);
-  }
-  
 }
 
-//____________________________________________________________________________
-void  AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const 
+//________________________________________________________________________
+void AliAnaGammaJet::CreateOutputObjects()
 {  
-  //Search for the charged particle with highest with 
-  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
-  Double_t pt  = -100.;
-  Double_t phi = -100.;
 
-  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
-
-    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
-
-    Double_t ptl  = particle->Pt();
-    Double_t rat  = ptl/pGamma->Pt() ;
-    Double_t phil = particle->Phi() ;
-    Double_t phig = pGamma->Phi();
-
-    //Selection within angular and energy limits
-    if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
-        (rat > fRatioMinCut) && (rat < fRatioMaxCut)  && (ptl  > pt)) {
-       phi = phil ;
-       pt  = ptl ;
-       pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
-       AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
-     }
-  }
+  // Init parameteres and create histograms to be saved in output file and 
+  // stores them in fOutputContainer
+  InitParameters();
+  AliAnaGammaDirect::CreateOutputObjects();
   
-  AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
-
-}
-
-
-//____________________________________________________________________________
-void  AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)  
-{  
-
-  //Search for the neutral pion with highest with 
-  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
-  Double_t pt  = -100.;
-  Double_t phi = -100.;
-  Double_t ptl = -100.;
-  Double_t rat = -100.; 
-  Double_t phil = -100. ;
-  Double_t ptg  = pGamma->Pt();
-  Double_t phig = pGamma->Phi();
+  fOutputContainer = new TObjArray(100) ;
+ 
+  TObjArray  * outputContainer =GetOutputContainer();
+  for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
+    fOutputContainer->Add(outputContainer->At(i)) ;
 
-  TIter next(pl);
-  TParticle * particle = 0;
+  //
+  fhChargeRatio  = new TH2F
+    ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
+     120,0,120,120,0,1); 
+  fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
+  fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhChargeRatio) ;
   
-  Int_t iPrimary = -1;
-  TLorentzVector gammai,gammaj;
-  Double_t angle = 0., e = 0., invmass = 0.;
-  Double_t anglef = 0., ef = 0., invmassf = 0.;
-  Int_t ksPdg = 0;
-  Int_t jPrimary=-1;
-
-  while ( (particle = (TParticle*)next()) ) {
-    iPrimary++;          
+  fhDeltaPhiCharge  = new TH2F
+    ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
+     200,0,120,200,0,6.4); 
+  fhDeltaPhiCharge->SetYTitle("#Delta #phi");
+  fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhDeltaPhiCharge) ; 
+  
+  fhDeltaEtaCharge  = new TH2F
+    ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
+     200,0,120,200,-2,2); 
+  fhDeltaEtaCharge->SetYTitle("#Delta #eta");
+  fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
+  fOutputContainer->Add(fhDeltaEtaCharge) ; 
+  
+  //
+  if(!fJetsOnlyInCTS){
+    fhPi0Ratio  = new TH2F
+      ("Pi0Ratio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
+       120,0,120,120,0,1); 
+    fhPi0Ratio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
+    fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhPi0Ratio) ; 
     
-    ksPdg = particle->GetPdgCode();
-    ptl  = particle->Pt();
-    if(ksPdg == 111){ //2 gamma overlapped, found with PID
-      rat = ptl/ptg ;
-      phil = particle->Phi() ;
-      //Selection within angular and energy limits
-      if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) && 
-        ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
-       phi = phil ;
-       pt  = ptl ;
-       pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
-       AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
-      }// cuts
-    }// pdg = 111
-    else if(ksPdg == 22){//1 gamma
-      //Search the photon companion in case it comes from  a Pi0 decay
-      //Apply several cuts to select the good pair
-      particle->Momentum(gammai);
-      jPrimary=-1;
-      TIter next2(pl);
-      while ( (particle = (TParticle*)next2()) ) {
-       jPrimary++;
-       if(jPrimary>iPrimary){
-         ksPdg = particle->GetPdgCode();
-
-         if(ksPdg == 22){
-           particle->Momentum(gammaj);
-           //Info("GetLeadingPi0","Egammai %f, Egammaj %f", 
-           //gammai.Pt(),gammaj.Pt());
-           ptl  = (gammai+gammaj).Pt();
-           phil = (gammai+gammaj).Phi();
-           if(phil < 0)
-             phil+=TMath::TwoPi();
-           rat = ptl/ptg ;
-           invmass = (gammai+gammaj).M();
-           angle = gammaj.Angle(gammai.Vect());
-           //Info("GetLeadingPi0","Angle %f", angle);
-           e = (gammai+gammaj).E();
-           //Fill histograms with no cuts applied.
-           fhAnglePairNoCut->Fill(e,angle);
-           fhInvMassPairNoCut->Fill(ptg,invmass);
-           //First cut on the energy and azimuth of the pair
-           if((rat > fRatioMinCut) && (rat < fRatioMaxCut) && 
-              ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
-             
-             fhAnglePairLeadingCut->Fill(e,angle);
-             fhInvMassPairLeadingCut->Fill(ptg,invmass);
-             //Second cut on the aperture of the pair
-             if(IsAngleInWindow(angle,e)){
-               fhAnglePairAngleCut->Fill(e,angle);
-               fhInvMassPairAngleCut->Fill(ptg,invmass);
-               
-               //Info("GetLeadingPi0","InvMass %f", invmass);
-               //Third cut on the invariant mass of the pair
-               if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
-                 fhInvMassPairAllCut->Fill(ptg,invmass);
-                 fhAnglePairAllCut->Fill(e,angle);
-                 if(ptl > pt ){
-                   pt       = ptl;
-                   phi      = phil ;
-                   ef       = e ;
-                   anglef   = angle ;
-                   invmassf = invmass ;
-                   pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
-                   AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
-                 }
-               }//cuts
-             }//(invmass>0.125) && (invmass<0.145)
-           }//gammaj.Angle(gammai.Vect())<0.04
-         }//if pdg = 22
-       }//iprimary<jprimary
-      }//while
-    }// if pdg = 22
-    //     }
-  }//while
-  
-  if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms
-    fhInvMassPairLeading->Fill(ptg,invmassf);
-    fhAnglePairLeading->Fill(ef,anglef);
+    fhDeltaPhiPi0  = new TH2F
+      ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
+       200,0,120,200,0,6.4); 
+    fhDeltaPhiPi0->SetYTitle("#Delta #phi");
+    fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhDeltaPhiPi0) ; 
+    
+    fhDeltaEtaPi0  = new TH2F
+      ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
+       200,0,120,200,-2,2); 
+    fhDeltaEtaPi0->SetYTitle("#Delta #eta");
+    fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhDeltaEtaPi0) ; 
+ 
+    //
+    fhAnglePair  = new TH2F
+      ("AnglePair",
+       "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}",
+       200,0,50,200,0,0.2); 
+    fhAnglePair->SetYTitle("Angle (rad)");
+    fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePair) ; 
+    
+    fhAnglePairAccepted  = new TH2F
+      ("AnglePairAccepted",
+       "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
+       200,0,50,200,0,0.2); 
+    fhAnglePairAccepted->SetYTitle("Angle (rad)");
+    fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairAccepted) ; 
+    
+    fhAnglePairNoCut  = new TH2F
+      ("AnglePairNoCut",
+       "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2); 
+    fhAnglePairNoCut->SetYTitle("Angle (rad)");
+    fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairNoCut) ; 
+    
+    fhAnglePairLeadingCut  = new TH2F
+      ("AnglePairLeadingCut",
+       "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
+       200,0,50,200,0,0.2); 
+    fhAnglePairLeadingCut->SetYTitle("Angle (rad)");
+    fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairLeadingCut) ; 
+    
+    fhAnglePairAngleCut  = new TH2F
+      ("AnglePairAngleCut",
+       "Angle between all #gamma pair (angle + leading cut) vs p_{T  #pi^{0}}"
+       ,200,0,50,200,0,0.2); 
+    fhAnglePairAngleCut->SetYTitle("Angle (rad)");
+    fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairAngleCut) ;
+    
+    fhAnglePairAllCut  = new TH2F
+      ("AnglePairAllCut",
+       "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T  #pi^{0}}"
+       ,200,0,50,200,0,0.2); 
+    fhAnglePairAllCut->SetYTitle("Angle (rad)");
+    fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairAllCut) ; 
+    
+    fhAnglePairLeading  = new TH2F
+      ("AnglePairLeading",
+       "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
+       200,0,50,200,0,0.2); 
+    fhAnglePairLeading->SetYTitle("Angle (rad)");
+    fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
+    fOutputContainer->Add(fhAnglePairLeading) ; 
+    
+    //
+    fhInvMassPairNoCut  = new TH2F
+      ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairNoCut) ; 
+    
+    fhInvMassPairLeadingCut  = new TH2F
+      ("InvMassPairLeadingCut",
+       "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairLeadingCut) ; 
+    
+    fhInvMassPairAngleCut  = new TH2F
+      ("InvMassPairAngleCut",
+       "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairAngleCut) ; 
+    
+    fhInvMassPairAllCut  = new TH2F
+      ("InvMassPairAllCut",
+       "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairAllCut) ; 
+    
+    fhInvMassPairLeading  = new TH2F
+      ("InvMassPairLeading",
+       "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
+       120,0,120,360,0,0.5); 
+    fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
+    fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhInvMassPairLeading) ; 
   }
   
-  AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
-}
-
-//____________________________________________________________________________
-Bool_t  AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,  
-                                        TParticle * pGamma, TParticle * pLeading) 
-{
-  //Search Charged or Neutral leading particle, select the highest one.
-  
-  TParticle * pLeadingCh = new TParticle();
-  TParticle * pLeadingPi0 = new TParticle();
-  
-  Double_t ptg  =  pGamma->Pt(); 
-  Double_t phig = pGamma->Phi(); 
-  Double_t etag = pGamma->Eta(); 
   
-  if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS)
-    {
-      AliDebug(3, "GetLeadingPi0");
-      GetLeadingPi0   (plNe, pGamma, pLeadingPi0) ;
-      AliDebug(3, "GetLeadingCharge");
-      GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
-      
-      Double_t ptch = pLeadingCh->Pt(); 
-      Double_t phich = pLeadingCh->Phi(); 
-      Double_t etach = pLeadingCh->Eta(); 
-      Double_t ptpi = pLeadingPi0->Pt(); 
-      Double_t phipi = pLeadingPi0->Phi(); 
-      Double_t etapi = pLeadingPi0->Eta(); 
-
-      //Is leading cone inside EMCAL acceptance?
-      
-      Bool_t insidech = kFALSE ;
-      if((phich - fCone) >  fPhiEMCALCut[0] && 
-        (phich + fCone) <  fPhiEMCALCut[1] && 
-       (etach-fCone) < fEtaEMCALCut )
-       insidech = kTRUE ;
-      
-      Bool_t insidepi = kFALSE ;
-      if((phipi - fCone) >  fPhiEMCALCut[0] && 
-        (phipi + fCone) <  fPhiEMCALCut[1] &&
-       (etapi-fCone) < fEtaEMCALCut )
-       insidepi = kTRUE ;
-
-      AliDebug(2,Form("Leading:  charged pt %f, pi0 pt  %f",ptch,ptpi)) ;
-      
-      if (ptch > 0 || ptpi > 0)
-       {
-         if(insidech && (ptch > ptpi))
-           {
-             if (GetPrintInfo())
-               AliInfo("Leading found in CTS");
-             pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
-             AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
-             fhChargeRatio->Fill(ptg,ptch/pGamma->Pt());
-             fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich);
-             fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach);
-             return 1;
-           }
-         
-         else if((ptpi > ptch) && insidepi)
-           {
-             if (GetPrintInfo())
-               AliInfo("Leading found in EMCAL");
-             pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
-             AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
-             fhPi0Ratio     ->Fill(ptg,ptpi/ptg);
-             fhDeltaPhiPi0->Fill(ptg,phig-phipi);
-             fhDeltaEtaPi0->Fill(ptg,etag-etapi);
-             return 1;     
-           }
-
-         else{
-           AliDebug(3,"NO LEADING PARTICLE FOUND");}
-         return 0; 
-       }
-      else{
-       AliDebug(3,"NO LEADING PARTICLE FOUND");
-       return 0;
-      }
-    }
-
-  else
-    {
-      //No calorimeter present for Leading particle detection
-      GetLeadingCharge(plCTS, pGamma, pLeading) ;
-      Double_t ptch = pLeading->Pt(); 
-      Double_t phich = pLeading->Phi(); 
-      Double_t etach = pLeading->Eta(); 
-      if(ptch > 0){
-       fhChargeRatio->Fill(ptg,ptch/ptg);
-       fhDeltaPhiCharge->Fill(ptg,phig-phich);
-       fhDeltaEtaCharge->Fill(ptg,etag-etach);
-       AliDebug(3,Form("Leading found :  pt %f, phi %f, eta %f",ptch,phich,etach)) ;
-       return 1;
-      }
-      else
-       {
-         AliDebug(3,"NO LEADING PARTICLE FOUND");      
-         return 0;
-       }
-    }
-}
-
-  //____________________________________________________________________________
-void AliAnaGammaJet::Init(const Option_t * )
-{
-//   // Initialisation of branch container 
-  AliAnaGammaDirect::Init();
- 
-  //Initialize the parameters of the analysis.
-  //fCalorimeter="PHOS";
-  fAngleMaxParam.Set(4) ;
-  fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
-  fAngleMaxParam.AddAt(-0.25,1) ;
-  fAngleMaxParam.AddAt(0.025,2) ;
-  fAngleMaxParam.AddAt(-2e-4,3) ;
-  fSeveralConeAndPtCuts   = kFALSE ;
-  //fPrintInfo           = kTRUE;
-  fPbPb                = kFALSE ;
-  fInvMassMaxCut  = 0.16 ;
-  fInvMassMinCut  = 0.11 ;
-  //fJetsOnlyInCTS    = kTRUE ;
-  fEtaEMCALCut     = 0.7 ;
-  fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
-  fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
-  fPhiMaxCut      = 3.4 ;
-  fPhiMinCut      = 2.9 ;
-
-  //Jet selection parameters
-  //Fixed cut (old)
-  fRatioMaxCut    = 1.0 ;
-  fRatioMinCut    = 0.1 ; 
-  fJetRatioMaxCut = 1.2 ; 
-  fJetRatioMinCut = 0.3 ; 
-  fJetsOnlyInCTS = kFALSE ;
-  fJetCTSRatioMaxCut = 1.2 ;
-  fJetCTSRatioMinCut = 0.3 ;
-  fSelect         = 0  ;
-
-  //Cut depending on gamma energy
-
-  fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
-  //Reconstructed jet energy dependence parameters 
-  //e_jet = a1+e_gamma b2. 
-  //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
-  fJetE1[0] = -5.75; fJetE1[1] = -4.1;
-  fJetE2[0] = 1.005; fJetE2[1] = 1.05;
-
-  //Reconstructed sigma of jet energy dependence parameters 
-  //s_jet = a1+e_gamma b2. 
-  //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
-  fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
-  fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
-
-  //Background mean energy and RMS
-  //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
-  //Index 2-> (low pt jets)BKG > 0.5 GeV;
-  //Index > 2, same for CTS conf
-  fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
-  fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
-  fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0; 
-  fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2; 
-
-  //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
-  //limits for monoenergetic jets.
-  //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
-  //Index 2-> (low pt jets) BKG > 0.5 GeV;
-  //Index > 2, same for CTS conf
+  if(!fSeveralConeAndPtCuts){
+    
+    //Count
+    fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); 
+    fhNBkg->SetYTitle("counts");
+    fhNBkg->SetXTitle("N");
+    fOutputContainer->Add(fhNBkg) ; 
+    
+    fhNLeading  = new TH2F
+      ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); 
+    fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
+    fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
+    fOutputContainer->Add(fhNLeading) ; 
+    
+    fhNJet  = new TH1F("NJet","Accepted jets",240,0,120); 
+    fhNJet->SetYTitle("N");
+    fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
+    fOutputContainer->Add(fhNJet) ; 
+    
+    //Ratios and Pt dist of reconstructed (not selected) jets
+    //Jet
+    fhJetRatio  = new TH2F
+      ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
+       240,0,120,200,0,10);
+    fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+    fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhJetRatio) ; 
+    
+    fhJetPt  = new TH2F
+      ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
+    fhJetPt->SetYTitle("p_{T jet}");
+    fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhJetPt) ; 
+    
+    //Bkg
+    
+    fhBkgRatio  = new TH2F
+      ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
+       240,0,120,200,0,10);
+    fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
+    fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhBkgRatio) ;
+    
+    fhBkgPt  = new TH2F
+      ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
+    fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
+    fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhBkgPt) ;
+    
+    //Jet Distributions
+    
+    fhJetFragment  = 
+      new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
+              240,0.,120.,1000,0.,1.2); 
+    fhJetFragment->SetYTitle("x_{T}");
+    fhJetFragment->SetXTitle("p_{T #gamma}");
+    fOutputContainer->Add(fhJetFragment) ;
+    
+    fhBkgFragment  = new TH2F
+      ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
+       240,0.,120.,1000,0.,1.2);
+    fhBkgFragment->SetYTitle("x_{T}");
+    fhBkgFragment->SetXTitle("p_{T #gamma}");
+    fOutputContainer->Add(fhBkgFragment) ;
+    
+    fhJetPtDist  = 
+      new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
+    fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhJetPtDist) ;
+    
+    fhBkgPtDist  = new TH2F
+      ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
+    fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
+    fOutputContainer->Add(fhBkgPtDist) ;
 
-  fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; 
-  fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
-  fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; 
-  fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
-  fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ; 
-  fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
-  fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; 
-  fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
+  }
+  else{
+    //If we want to study the jet for different cones and pt
+    
+    for(Int_t icone = 0; icone<fNCone; icone++){
+      for(Int_t ipt = 0; ipt<fNPt;ipt++){ 
+       
+       //Jet
+       
+       fhJetRatios[icone][ipt]  = new TH2F
+         ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
+          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+          240,0,120,200,0,10);
+       fhJetRatios[icone][ipt]->
+         SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+       fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+       fOutputContainer->Add(fhJetRatios[icone][ipt]) ; 
+       
+       
+       fhJetPts[icone][ipt]  = new TH2F
+         ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
+          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+          240,0,120,400,0,200);
+       fhJetPts[icone][ipt]->
+         SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+       fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+       fOutputContainer->Add(fhJetPts[icone][ipt]) ; 
+       
+       //Bkg
+       fhBkgRatios[icone][ipt]  = new TH2F
+         ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
+          "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+          240,0,120,200,0,10);
+       fhBkgRatios[icone][ipt]->
+         SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
+       fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+       fOutputContainer->Add(fhBkgRatios[icone][ipt]) ; 
+       
+       fhBkgPts[icone][ipt]  = new TH2F
+         ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
+          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
+          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
+          240,0,120,400,0,200);
+       fhBkgPts[icone][ipt]->
+         SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
+       fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+       fOutputContainer->Add(fhBkgPts[icone][ipt]) ; 
+       
+       //Counts
+       fhNBkgs[icone][ipt]  = new TH1F
+         ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "bkg multiplicity cone ="+fNameCones[icone]+", pt>" 
+          +fNamePtThres[ipt]+" GeV/c",9000,0,9000); 
+       fhNBkgs[icone][ipt]->SetYTitle("counts");
+       fhNBkgs[icone][ipt]->SetXTitle("N");
+       fOutputContainer->Add(fhNBkgs[icone][ipt]) ; 
+       
+       fhNLeadings[icone][ipt]  = new TH2F
+         ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>" 
+          +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); 
+       fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
+       fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
+       fOutputContainer->Add(fhNLeadings[icone][ipt]) ; 
+       
+       fhNJets[icone][ipt]  = new TH1F
+         ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "Number of neutral jets, cone ="+fNameCones[icone]+", pt>" 
+          +fNamePtThres[ipt]+" GeV/c",120,0,120); 
+       fhNJets[icone][ipt]->SetYTitle("N");
+       fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
+       fOutputContainer->Add(fhNJets[icone][ipt]) ; 
+       
+       //Fragmentation Function
+       fhJetFragments[icone][ipt]  = new TH2F
+         ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" 
+          +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
+       fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
+       fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
+       fOutputContainer->Add(fhJetFragments[icone][ipt]) ; 
+       
+       fhBkgFragments[icone][ipt]  = new TH2F
+         ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" 
+          +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
+       fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
+       fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
+       fOutputContainer->Add(fhBkgFragments[icone][ipt]) ; 
+       
+       //Jet particle distribution
+       
+       fhJetPtDists[icone][ipt]  = new TH2F
+         ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
+          " GeV/c",120,0.,120.,120,0.,120.); 
+       fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+       fOutputContainer->Add(fhJetPtDists[icone][ipt]) ; 
+       
+       fhBkgPtDists[icone][ipt]  = new TH2F
+         ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
+          "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
+          " GeV/c",120,0.,120.,120,0.,120.); 
+       fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
+       fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ; 
+       
+      }//ipt
+    } //icone
+  }//If we want to study any cone or pt threshold
+}
 
 
-  //Different cones and pt thresholds to construct the jet
+//____________________________________________________________________________
+void AliAnaGammaJet::Exec(Option_t *) 
+{
+  
+  // Processing of one event
 
-  fCone        = 0.3  ;
-  fPtThreshold = 0.   ;
-  fNCone       = 3    ;
-  fNPt         = 3    ;
-  fCones[1]    = 0.2  ; fNameCones[1]   = "02" ;
-  fPtThres[0]  = 0.0  ; fNamePtThres[0] = "00" ;
-  fCones[0]    = 0.3  ; fNameCones[0]   = "03" ;
-  fPtThres[1]  = 0.5  ; fNamePtThres[1] = "05" ;
-  fCones[2]    = 0.4  ; fNameCones[2]   = "04" ;
-  fPtThres[2]  = 1.0  ; fNamePtThres[2] = "10" ;
-  //Fill particle lists when PID is ok
-  // fEMCALPID = kFALSE;
-  // fPHOSPID = kFALSE;
+  //Get ESDs
+  Long64_t entry = GetChain()->GetReadEntry() ;
+  
+  if (!GetESD()) {
+    AliError("fESD is not connected to the input!") ; 
+    return ; 
+  }
+  
+  if (GetPrintInfo()) 
+    AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ; 
+  
+  //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
+  
+  TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
+  TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
+  TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter
+  TClonesArray * plCalo     = new TClonesArray("TParticle",1000);  // All particles measured in Prompt Gamma calorimeter 
+  
+  
+  TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
+  TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
+  
+  Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
+  
+  //Fill lists with photons, neutral particles and charged particles
+  //look for the highest energy photon in the event inside fCalorimeter
+  AliDebug(2, "Fill particle lists, get prompt gamma");
+  
+  //Fill particle lists 
+  
+  if(GetCalorimeter() == "PHOS")
+    CreateParticleList(particleList, plCTS,plNe,plCalo); 
+  else if(GetCalorimeter() == "EMCAL")
+    CreateParticleList(particleList, plCTS,plCalo,plNe); 
+  else
+    AliError("No calorimeter selected");
+  
+  //Search highest energy prompt gamma in calorimeter
+  GetPromptGamma(plCalo,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
+  
+  AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
+  AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
+                 plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries()));
 
-  //Initialization of histograms 
+  //If there is any prompt photon in fCalorimeter, 
+  //search jet leading particle
+  
+  if(iIsInPHOSorEMCAL){
+    if (GetPrintInfo())
+      AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+ 
+    AliDebug(2, "Get Leading Particles, Make jet");
 
-  MakeHistos() ; 
+    //Search leading particles in CTS and EMCAL 
+    if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
+      if (GetPrintInfo())
+       AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
 
-}
+      //Search Jet
+      if(!fSeveralConeAndPtCuts)
+       MakeJet(particleList,pGamma,pLeading,"");
+      else{
+       for(Int_t icone = 0; icone<fNCone; icone++) {
+         for(Int_t ipt = 0; ipt<fNPt;ipt++) {  
+           TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt];
+           MakeJet(particleList, pGamma, pLeading,lastname);
+         }//icone
+       }//ipt
+      }//fSeveralConeAndPtCuts
+    }//Leading
+  }//Gamma in Calo
+     
+  AliDebug(2, "End of analysis, delete pointers");
 
-//__________________________________________________________________________-
-Bool_t AliAnaGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
-  //Check if the opening angle of the candidate pairs is inside 
-  //our selection windowd
+  particleList->Delete() ; 
+  plCTS->Delete() ;
+  plNe->Delete() ;
+  plCalo->Delete() ;
+  pLeading->Delete();
+  pGamma->Delete();
 
-  Bool_t result = kFALSE;
-  Double_t mpi0 = 0.1349766;
-  Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
-    +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
-  Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
-  Double_t min = 100. ;
-  if(arg>0.)
-    min = TMath::ACos(arg);
+  delete plNe ;
+  delete plCalo ;
+  delete plCTS ;
+  delete particleList ;
+  //  delete pLeading;
+  //  delete pGamma;
 
-  if((angle<max)&&(angle>=min))
-    result = kTRUE;
- 
-  return result;
-}
+  PostData(0, fOutputContainer);
+}    
 
-//__________________________________________________________________________-
-Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){
-  //Check if the energy of the reconstructed jet is within an energy window
+//____________________________________________________________________________
+void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
+{
+  //Fill histograms wth jet fragmentation 
+  //and number of selected jets and leading particles
+  //and the background multiplicity
+  TParticle * particle = 0 ;
+  Int_t ipr = 0;
+  Float_t  charge = 0;
 
-  Double_t par[6];
-  Double_t xmax[2];
-  Double_t xmin[2];
+  TIter next(pl) ; 
+  while ( (particle = dynamic_cast<TParticle*>(next())) ) {
+    ipr++ ;
+    Double_t pt = particle->Pt();
 
-  Int_t iCTS = 0;
-  if(fJetsOnlyInCTS)
-    iCTS = 3 ;
+    charge = TDatabasePDG::Instance()
+      ->GetParticle(particle->GetPdgCode())->Charge();
+    if(charge != 0){//Only jet Charged particles 
+      dynamic_cast<TH2F*>
+       (fOutputContainer->FindObject(type+"Fragment"+lastname))
+       ->Fill(ptg,pt/ptg);
+      dynamic_cast<TH2F*>
+       (fOutputContainer->FindObject(type+"PtDist"+lastname))
+       ->Fill(ptg,pt);
+    }//charged
 
-  if(!fPbPb){
-    //Phythia alone, jets with pt_th > 0.2, r = 0.3 
-    par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
-    //Energy of the jet peak
-    //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
-    par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
-    //Sigma  of the jet peak
-    //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
-    par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
-    //Parameters reserved for PbPb bkg.
-    xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
-    xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
-    //Factor that multiplies sigma to obtain the best limits, 
-    //by observation, of mono jet ratios (ptjet/ptg)
-    //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
-   
-  }
-  else{
-    if(ptg > fPtJetSelectionCut){
-      //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3 
-      par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
-      //Energy of the jet peak, same as in pp
-      //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
-      par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
-      //Sigma  of the jet peak, same as in pp
-      //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
-      par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
-      //Mean value and RMS of PbPb Bkg 
-      xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
-      xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
-      //Factor that multiplies sigma to obtain the best limits, 
-      //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
-      //pt_th > 2 GeV, r = 0.3
-      //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
-     
-    }
-    else{
-      //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
-      par[0] = fJetE1[1]; par[1] = fJetE2[1]; 
-      //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 
-      //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
-      par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
-      //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
-      //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
-      par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
-      //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
-      xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
-      xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
-      //Factor that multiplies sigma to obtain the best limits, 
-      //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
-      //pt_th > 2 GeV, r = 0.3
-      //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
-     
-    }//If low pt jet in bkg
-  }//if Bkg
+  }//while
 
- //Calculate minimum and maximum limits of the jet ratio.
-  Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
-  Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
+  if(type == "Bkg")
+    dynamic_cast<TH1F*>
+      (fOutputContainer->FindObject("NBkg"+lastname))
+      ->Fill(ipr);
+  else{
+    dynamic_cast<TH1F*>
+      (fOutputContainer->FindObject("NJet"+lastname))->
+      Fill(ptg);
+    dynamic_cast<TH2F*>
+      (fOutputContainer->FindObject("NLeading"+lastname))
+      ->Fill(ptg,ptl);
+  }
   
-  AliDebug(3,Form("Jet selection?  : Limits min %f, max %f,  pt_jet %f,  pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg));
+}
 
-  if(( min < ptj/ptg ) && ( max > ptj/ptg))
-    return kTRUE;
-  else
-    return kFALSE;
+//____________________________________________________________________________
+void  AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const 
+{  
+  //Search for the charged particle with highest with 
+  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
+  Double_t pt  = -100.;
+  Double_t phi = -100.;
+
+  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+
+    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+    Double_t ptl  = particle->Pt();
+    Double_t rat  = ptl/pGamma->Pt() ;
+    Double_t phil = particle->Phi() ;
+    Double_t phig = pGamma->Phi();
+
+    //Selection within angular and energy limits
+    if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
+        (rat > fRatioMinCut) && (rat < fRatioMaxCut)  && (ptl  > pt)) {
+       phi = phil ;
+       pt  = ptl ;
+       pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+       AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
+     }
+  }
+  
+  AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
 
 }
 
+
 //____________________________________________________________________________
-void AliAnaGammaJet::MakeHistos()
-{
-  // Create histograms to be saved in output file and 
-  // stores them in fOutputContainer
-  
-  fOutputContainer = new TObjArray(10000) ;
- 
-  TObjArray  * outputContainer =GetOutputContainer();
-  for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
-    fOutputContainer->Add(outputContainer->At(i)) ;
+void  AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)  
+{  
 
-  //
-  fhChargeRatio  = new TH2F
-    ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
-     120,0,120,120,0,1); 
-  fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
-  fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhChargeRatio) ;
-  
-  fhDeltaPhiCharge  = new TH2F
-    ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
-     200,0,120,200,0,6.4); 
-  fhDeltaPhiCharge->SetYTitle("#Delta #phi");
-  fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhDeltaPhiCharge) ; 
-  
-  fhDeltaEtaCharge  = new TH2F
-    ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
-     200,0,120,200,-2,2); 
-  fhDeltaEtaCharge->SetYTitle("#Delta #eta");
-  fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
-  fOutputContainer->Add(fhDeltaEtaCharge) ; 
+  //Search for the neutral pion with highest with 
+  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
+  Double_t pt  = -100.;
+  Double_t phi = -100.;
+  Double_t ptl = -100.;
+  Double_t rat = -100.; 
+  Double_t phil = -100. ;
+  Double_t ptg  = pGamma->Pt();
+  Double_t phig = pGamma->Phi();
+
+  TIter next(pl);
+  TParticle * particle = 0;
   
-  //
-  if(!fJetsOnlyInCTS){
-    fhPi0Ratio  = new TH2F
-      ("Pi0Ratio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
-       120,0,120,120,0,1); 
-    fhPi0Ratio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
-    fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhPi0Ratio) ; 
-    
-    fhDeltaPhiPi0  = new TH2F
-      ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
-       200,0,120,200,0,6.4); 
-    fhDeltaPhiPi0->SetYTitle("#Delta #phi");
-    fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhDeltaPhiPi0) ; 
-    
-    fhDeltaEtaPi0  = new TH2F
-      ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
-       200,0,120,200,-2,2); 
-    fhDeltaEtaPi0->SetYTitle("#Delta #eta");
-    fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhDeltaEtaPi0) ; 
- 
-    //
-    fhAnglePair  = new TH2F
-      ("AnglePair",
-       "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}",
-       200,0,50,200,0,0.2); 
-    fhAnglePair->SetYTitle("Angle (rad)");
-    fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePair) ; 
-    
-    fhAnglePairAccepted  = new TH2F
-      ("AnglePairAccepted",
-       "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
-       200,0,50,200,0,0.2); 
-    fhAnglePairAccepted->SetYTitle("Angle (rad)");
-    fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairAccepted) ; 
-    
-    fhAnglePairNoCut  = new TH2F
-      ("AnglePairNoCut",
-       "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2); 
-    fhAnglePairNoCut->SetYTitle("Angle (rad)");
-    fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairNoCut) ; 
-    
-    fhAnglePairLeadingCut  = new TH2F
-      ("AnglePairLeadingCut",
-       "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
-       200,0,50,200,0,0.2); 
-    fhAnglePairLeadingCut->SetYTitle("Angle (rad)");
-    fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairLeadingCut) ; 
-    
-    fhAnglePairAngleCut  = new TH2F
-      ("AnglePairAngleCut",
-       "Angle between all #gamma pair (angle + leading cut) vs p_{T  #pi^{0}}"
-       ,200,0,50,200,0,0.2); 
-    fhAnglePairAngleCut->SetYTitle("Angle (rad)");
-    fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairAngleCut) ;
-    
-    fhAnglePairAllCut  = new TH2F
-      ("AnglePairAllCut",
-       "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T  #pi^{0}}"
-       ,200,0,50,200,0,0.2); 
-    fhAnglePairAllCut->SetYTitle("Angle (rad)");
-    fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairAllCut) ; 
-    
-    fhAnglePairLeading  = new TH2F
-      ("AnglePairLeading",
-       "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
-       200,0,50,200,0,0.2); 
-    fhAnglePairLeading->SetYTitle("Angle (rad)");
-    fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
-    fOutputContainer->Add(fhAnglePairLeading) ; 
-    
-    //
-    fhInvMassPairNoCut  = new TH2F
-      ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairNoCut) ; 
-    
-    fhInvMassPairLeadingCut  = new TH2F
-      ("InvMassPairLeadingCut",
-       "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairLeadingCut) ; 
-    
-    fhInvMassPairAngleCut  = new TH2F
-      ("InvMassPairAngleCut",
-       "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairAngleCut) ; 
-    
-    fhInvMassPairAllCut  = new TH2F
-      ("InvMassPairAllCut",
-       "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairAllCut) ; 
+  Int_t iPrimary = -1;
+  TLorentzVector gammai,gammaj;
+  Double_t angle = 0., e = 0., invmass = 0.;
+  Double_t anglef = 0., ef = 0., invmassf = 0.;
+  Int_t ksPdg = 0;
+  Int_t jPrimary=-1;
+
+  while ( (particle = (TParticle*)next()) ) {
+    iPrimary++;          
     
-    fhInvMassPairLeading  = new TH2F
-      ("InvMassPairLeading",
-       "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
-       120,0,120,360,0,0.5); 
-    fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
-    fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhInvMassPairLeading) ; 
+    ksPdg = particle->GetPdgCode();
+    ptl  = particle->Pt();
+    if(ksPdg == 111){ //2 gamma overlapped, found with PID
+      rat = ptl/ptg ;
+      phil = particle->Phi() ;
+      //Selection within angular and energy limits
+      if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) && 
+        ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
+       phi = phil ;
+       pt  = ptl ;
+       pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+       AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
+      }// cuts
+    }// pdg = 111
+    else if(ksPdg == 22){//1 gamma
+      //Search the photon companion in case it comes from  a Pi0 decay
+      //Apply several cuts to select the good pair
+      particle->Momentum(gammai);
+      jPrimary=-1;
+      TIter next2(pl);
+      while ( (particle = (TParticle*)next2()) ) {
+       jPrimary++;
+       if(jPrimary>iPrimary){
+         ksPdg = particle->GetPdgCode();
+
+         if(ksPdg == 22){
+           particle->Momentum(gammaj);
+           //Info("GetLeadingPi0","Egammai %f, Egammaj %f", 
+           //gammai.Pt(),gammaj.Pt());
+           ptl  = (gammai+gammaj).Pt();
+           phil = (gammai+gammaj).Phi();
+           if(phil < 0)
+             phil+=TMath::TwoPi();
+           rat = ptl/ptg ;
+           invmass = (gammai+gammaj).M();
+           angle = gammaj.Angle(gammai.Vect());
+           //Info("GetLeadingPi0","Angle %f", angle);
+           e = (gammai+gammaj).E();
+           //Fill histograms with no cuts applied.
+           fhAnglePairNoCut->Fill(e,angle);
+           fhInvMassPairNoCut->Fill(ptg,invmass);
+           //First cut on the energy and azimuth of the pair
+           if((rat > fRatioMinCut) && (rat < fRatioMaxCut) && 
+              ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
+             
+             fhAnglePairLeadingCut->Fill(e,angle);
+             fhInvMassPairLeadingCut->Fill(ptg,invmass);
+             //Second cut on the aperture of the pair
+             if(IsAngleInWindow(angle,e)){
+               fhAnglePairAngleCut->Fill(e,angle);
+               fhInvMassPairAngleCut->Fill(ptg,invmass);
+               
+               //Info("GetLeadingPi0","InvMass %f", invmass);
+               //Third cut on the invariant mass of the pair
+               if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
+                 fhInvMassPairAllCut->Fill(ptg,invmass);
+                 fhAnglePairAllCut->Fill(e,angle);
+                 if(ptl > pt ){
+                   pt       = ptl;
+                   phi      = phil ;
+                   ef       = e ;
+                   anglef   = angle ;
+                   invmassf = invmass ;
+                   pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+                   AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
+                 }
+               }//cuts
+             }//(invmass>0.125) && (invmass<0.145)
+           }//gammaj.Angle(gammai.Vect())<0.04
+         }//if pdg = 22
+       }//iprimary<jprimary
+      }//while
+    }// if pdg = 22
+    //     }
+  }//while
+  
+  if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms
+    fhInvMassPairLeading->Fill(ptg,invmassf);
+    fhAnglePairLeading->Fill(ef,anglef);
   }
   
+  AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
+}
+
+//____________________________________________________________________________
+Bool_t  AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,  
+                                        TParticle * pGamma, TParticle * pLeading) 
+{
+  //Search Charged or Neutral leading particle, select the highest one.
   
-  if(!fSeveralConeAndPtCuts){
-    
-    //Count
-    fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); 
-    fhNBkg->SetYTitle("counts");
-    fhNBkg->SetXTitle("N");
-    fOutputContainer->Add(fhNBkg) ; 
-    
-    fhNLeading  = new TH2F
-      ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); 
-    fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
-    fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
-    fOutputContainer->Add(fhNLeading) ; 
-    
-    fhNJet  = new TH1F("NJet","Accepted jets",240,0,120); 
-    fhNJet->SetYTitle("N");
-    fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
-    fOutputContainer->Add(fhNJet) ; 
-    
-    //Ratios and Pt dist of reconstructed (not selected) jets
-    //Jet
-    fhJetRatio  = new TH2F
-      ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
-       240,0,120,200,0,10);
-    fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
-    fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhJetRatio) ; 
-    
-    fhJetPt  = new TH2F
-      ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
-    fhJetPt->SetYTitle("p_{T jet}");
-    fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhJetPt) ; 
-    
-    //Bkg
-    
-    fhBkgRatio  = new TH2F
-      ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
-       240,0,120,200,0,10);
-    fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
-    fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhBkgRatio) ;
-    
-    fhBkgPt  = new TH2F
-      ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
-    fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
-    fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhBkgPt) ;
-    
-    //Jet Distributions
-    
-    fhJetFragment  = 
-      new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
-              240,0.,120.,1000,0.,1.2); 
-    fhJetFragment->SetYTitle("x_{T}");
-    fhJetFragment->SetXTitle("p_{T #gamma}");
-    fOutputContainer->Add(fhJetFragment) ;
-    
-    fhBkgFragment  = new TH2F
-      ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
-       240,0.,120.,1000,0.,1.2);
-    fhBkgFragment->SetYTitle("x_{T}");
-    fhBkgFragment->SetXTitle("p_{T #gamma}");
-    fOutputContainer->Add(fhBkgFragment) ;
-    
-    fhJetPtDist  = 
-      new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
-    fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhJetPtDist) ;
-    
-    fhBkgPtDist  = new TH2F
-      ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
-    fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
-    fOutputContainer->Add(fhBkgPtDist) ;
+  TParticle * pLeadingCh = new TParticle();
+  TParticle * pLeadingPi0 = new TParticle();
+  
+  Double_t ptg  =  pGamma->Pt(); 
+  Double_t phig = pGamma->Phi(); 
+  Double_t etag = pGamma->Eta(); 
+  
+  if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS)
+    {
+      AliDebug(3, "GetLeadingPi0");
+      GetLeadingPi0   (plNe, pGamma, pLeadingPi0) ;
+      AliDebug(3, "GetLeadingCharge");
+      GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
+      
+      Double_t ptch = pLeadingCh->Pt(); 
+      Double_t phich = pLeadingCh->Phi(); 
+      Double_t etach = pLeadingCh->Eta(); 
+      Double_t ptpi = pLeadingPi0->Pt(); 
+      Double_t phipi = pLeadingPi0->Phi(); 
+      Double_t etapi = pLeadingPi0->Eta(); 
+
+      //Is leading cone inside EMCAL acceptance?
+      
+      Bool_t insidech = kFALSE ;
+      if((phich - fCone) >  fPhiEMCALCut[0] && 
+        (phich + fCone) <  fPhiEMCALCut[1] && 
+       (etach-fCone) < fEtaEMCALCut )
+       insidech = kTRUE ;
+      
+      Bool_t insidepi = kFALSE ;
+      if((phipi - fCone) >  fPhiEMCALCut[0] && 
+        (phipi + fCone) <  fPhiEMCALCut[1] &&
+       (etapi-fCone) < fEtaEMCALCut )
+       insidepi = kTRUE ;
+
+      AliDebug(2,Form("Leading:  charged pt %f, pi0 pt  %f",ptch,ptpi)) ;
+      
+      if (ptch > 0 || ptpi > 0)
+       {
+         if(insidech && (ptch > ptpi))
+           {
+             if (GetPrintInfo())
+               AliInfo("Leading found in CTS");
+             pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
+             AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
+             fhChargeRatio->Fill(ptg,ptch/pGamma->Pt());
+             fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich);
+             fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach);
+             return 1;
+           }
+         
+         else if((ptpi > ptch) && insidepi)
+           {
+             if (GetPrintInfo())
+               AliInfo("Leading found in EMCAL");
+             pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
+             AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
+             fhPi0Ratio     ->Fill(ptg,ptpi/ptg);
+             fhDeltaPhiPi0->Fill(ptg,phig-phipi);
+             fhDeltaEtaPi0->Fill(ptg,etag-etapi);
+             return 1;     
+           }
+
+         else{
+           AliDebug(3,"NO LEADING PARTICLE FOUND");}
+         return 0; 
+       }
+      else{
+       AliDebug(3,"NO LEADING PARTICLE FOUND");
+       return 0;
+      }
+    }
+
+  else
+    {
+      //No calorimeter present for Leading particle detection
+      GetLeadingCharge(plCTS, pGamma, pLeading) ;
+      Double_t ptch = pLeading->Pt(); 
+      Double_t phich = pLeading->Phi(); 
+      Double_t etach = pLeading->Eta(); 
+      if(ptch > 0){
+       fhChargeRatio->Fill(ptg,ptch/ptg);
+       fhDeltaPhiCharge->Fill(ptg,phig-phich);
+       fhDeltaEtaCharge->Fill(ptg,etag-etach);
+       AliDebug(3,Form("Leading found :  pt %f, phi %f, eta %f",ptch,phich,etach)) ;
+       return 1;
+      }
+      else
+       {
+         AliDebug(3,"NO LEADING PARTICLE FOUND");      
+         return 0;
+       }
+    }
+}
+
+  //____________________________________________________________________________
+void AliAnaGammaJet::InitParameters()
+{
+//   // Initialisation of branch container 
+  //AliAnaGammaDirect::InitParameters();
+ 
+  //Initialize the parameters of the analysis.
+  //fCalorimeter="PHOS";
+  fAngleMaxParam.Set(4) ;
+  fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
+  fAngleMaxParam.AddAt(-0.25,1) ;
+  fAngleMaxParam.AddAt(0.025,2) ;
+  fAngleMaxParam.AddAt(-2e-4,3) ;
+  fSeveralConeAndPtCuts   = kFALSE ;
+  //fPrintInfo           = kTRUE;
+  fPbPb                = kFALSE ;
+  fInvMassMaxCut  = 0.16 ;
+  fInvMassMinCut  = 0.11 ;
+  //fJetsOnlyInCTS    = kTRUE ;
+  fEtaEMCALCut     = 0.7 ;
+  fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
+  fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
+  fPhiMaxCut      = 3.4 ;
+  fPhiMinCut      = 2.9 ;
+
+  //Jet selection parameters
+  //Fixed cut (old)
+  fRatioMaxCut    = 1.0 ;
+  fRatioMinCut    = 0.1 ; 
+  fJetRatioMaxCut = 1.2 ; 
+  fJetRatioMinCut = 0.3 ; 
+  fJetsOnlyInCTS = kFALSE ;
+  fJetCTSRatioMaxCut = 1.2 ;
+  fJetCTSRatioMinCut = 0.3 ;
+  fSelect         = 0  ;
+
+  //Cut depending on gamma energy
+
+  fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
+  //Reconstructed jet energy dependence parameters 
+  //e_jet = a1+e_gamma b2. 
+  //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
+  fJetE1[0] = -5.75; fJetE1[1] = -4.1;
+  fJetE2[0] = 1.005; fJetE2[1] = 1.05;
 
+  //Reconstructed sigma of jet energy dependence parameters 
+  //s_jet = a1+e_gamma b2. 
+  //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
+  fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
+  fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
+
+  //Background mean energy and RMS
+  //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
+  //Index 2-> (low pt jets)BKG > 0.5 GeV;
+  //Index > 2, same for CTS conf
+  fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
+  fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
+  fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0; 
+  fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2; 
+
+  //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
+  //limits for monoenergetic jets.
+  //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
+  //Index 2-> (low pt jets) BKG > 0.5 GeV;
+  //Index > 2, same for CTS conf
+
+  fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; 
+  fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
+  fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; 
+  fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
+  fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ; 
+  fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
+  fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; 
+  fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
+
+
+  //Different cones and pt thresholds to construct the jet
+
+  fCone        = 0.3  ;
+  fPtThreshold = 0.   ;
+  fNCone       = 3    ;
+  fNPt         = 3    ;
+  fCones[1]    = 0.2  ; fNameCones[1]   = "02" ;
+  fPtThres[0]  = 0.0  ; fNamePtThres[0] = "00" ;
+  fCones[0]    = 0.3  ; fNameCones[0]   = "03" ;
+  fPtThres[1]  = 0.5  ; fNamePtThres[1] = "05" ;
+  fCones[2]    = 0.4  ; fNameCones[2]   = "04" ;
+  fPtThres[2]  = 1.0  ; fNamePtThres[2] = "10" ;
+  //Fill particle lists when PID is ok
+  // fEMCALPID = kFALSE;
+  // fPHOSPID = kFALSE;
+
+}
+
+//__________________________________________________________________________-
+Bool_t AliAnaGammaJet::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)
+    +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
+  Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
+  Double_t min = 100. ;
+  if(arg>0.)
+    min = TMath::ACos(arg);
+
+  if((angle<max)&&(angle>=min))
+    result = kTRUE;
+ 
+  return result;
+}
+
+//__________________________________________________________________________-
+Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){
+  //Check if the energy of the reconstructed jet is within an energy window
+
+  Double_t par[6];
+  Double_t xmax[2];
+  Double_t xmin[2];
+
+  Int_t iCTS = 0;
+  if(fJetsOnlyInCTS)
+    iCTS = 3 ;
+
+  if(!fPbPb){
+    //Phythia alone, jets with pt_th > 0.2, r = 0.3 
+    par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
+    //Energy of the jet peak
+    //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
+    par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
+    //Sigma  of the jet peak
+    //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
+    par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
+    //Parameters reserved for PbPb bkg.
+    xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
+    xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
+    //Factor that multiplies sigma to obtain the best limits, 
+    //by observation, of mono jet ratios (ptjet/ptg)
+    //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
+   
   }
   else{
-    //If we want to study the jet for different cones and pt
-    
-    for(Int_t icone = 0; icone<fNCone; icone++){
-      for(Int_t ipt = 0; ipt<fNPt;ipt++){ 
-       
-       //Jet
-       
-       fhJetRatios[icone][ipt]  = new TH2F
-         ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
-          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
-          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
-          240,0,120,200,0,10);
-       fhJetRatios[icone][ipt]->
-         SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
-       fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
-       fOutputContainer->Add(fhJetRatios[icone][ipt]) ; 
-       
-       
-       fhJetPts[icone][ipt]  = new TH2F
-         ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
-          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
-          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
-          240,0,120,400,0,200);
-       fhJetPts[icone][ipt]->
-         SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
-       fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
-       fOutputContainer->Add(fhJetPts[icone][ipt]) ; 
-       
-       //Bkg
-       fhBkgRatios[icone][ipt]  = new TH2F
-         ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
-          "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
-          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
-          240,0,120,200,0,10);
-       fhBkgRatios[icone][ipt]->
-         SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
-       fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
-       fOutputContainer->Add(fhBkgRatios[icone][ipt]) ; 
-       
-       fhBkgPts[icone][ipt]  = new TH2F
-         ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
-          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
-          +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
-          240,0,120,400,0,200);
-       fhBkgPts[icone][ipt]->
-         SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
-       fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
-       fOutputContainer->Add(fhBkgPts[icone][ipt]) ; 
-       
-       //Counts
-       fhNBkgs[icone][ipt]  = new TH1F
-         ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "bkg multiplicity cone ="+fNameCones[icone]+", pt>" 
-          +fNamePtThres[ipt]+" GeV/c",9000,0,9000); 
-       fhNBkgs[icone][ipt]->SetYTitle("counts");
-       fhNBkgs[icone][ipt]->SetXTitle("N");
-       fOutputContainer->Add(fhNBkgs[icone][ipt]) ; 
-       
-       fhNLeadings[icone][ipt]  = new TH2F
-         ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>" 
-          +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); 
-       fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
-       fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
-       fOutputContainer->Add(fhNLeadings[icone][ipt]) ; 
-       
-       fhNJets[icone][ipt]  = new TH1F
-         ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "Number of neutral jets, cone ="+fNameCones[icone]+", pt>" 
-          +fNamePtThres[ipt]+" GeV/c",120,0,120); 
-       fhNJets[icone][ipt]->SetYTitle("N");
-       fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
-       fOutputContainer->Add(fhNJets[icone][ipt]) ; 
-       
-       //Fragmentation Function
-       fhJetFragments[icone][ipt]  = new TH2F
-         ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" 
-          +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
-       fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
-       fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
-       fOutputContainer->Add(fhJetFragments[icone][ipt]) ; 
-       
-       fhBkgFragments[icone][ipt]  = new TH2F
-         ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" 
-          +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
-       fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
-       fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
-       fOutputContainer->Add(fhBkgFragments[icone][ipt]) ; 
-       
-       //Jet particle distribution
-       
-       fhJetPtDists[icone][ipt]  = new TH2F
-         ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
-          " GeV/c",120,0.,120.,120,0.,120.); 
-       fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
-       fOutputContainer->Add(fhJetPtDists[icone][ipt]) ; 
-       
-       fhBkgPtDists[icone][ipt]  = new TH2F
-         ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
-          "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
-          " GeV/c",120,0.,120.,120,0.,120.); 
-       fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
-       fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ; 
-       
-      }//ipt
-    } //icone
-  }//If we want to study any cone or pt threshold
+    if(ptg > fPtJetSelectionCut){
+      //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3 
+      par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
+      //Energy of the jet peak, same as in pp
+      //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
+      par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
+      //Sigma  of the jet peak, same as in pp
+      //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
+      par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
+      //Mean value and RMS of PbPb Bkg 
+      xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
+      xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
+      //Factor that multiplies sigma to obtain the best limits, 
+      //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
+      //pt_th > 2 GeV, r = 0.3
+      //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
+     
+    }
+    else{
+      //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
+      par[0] = fJetE1[1]; par[1] = fJetE2[1]; 
+      //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 
+      //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
+      par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
+      //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
+      //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
+      par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
+      //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
+      xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
+      xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
+      //Factor that multiplies sigma to obtain the best limits, 
+      //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
+      //pt_th > 2 GeV, r = 0.3
+      //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
+     
+    }//If low pt jet in bkg
+  }//if Bkg
+
+ //Calculate minimum and maximum limits of the jet ratio.
+  Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
+  Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
+  
+  AliDebug(3,Form("Jet selection?  : Limits min %f, max %f,  pt_jet %f,  pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg));
+
+  if(( min < ptj/ptg ) && ( max > ptj/ptg))
+    return kTRUE;
+  else
+    return kFALSE;
+
 }
 
 //____________________________________________________________________________

diff --git a/PWG4/AliAnaGammaJet.h b/PWG4/AliAnaGammaJet.h
index b6714c2b6695f671964c2de269b1027978fd4d7a..da86d63d622e0308355d118b203f737956a5ffe5 100644 (file)
--- a/PWG4/AliAnaGammaJet.h
+++ b/PWG4/AliAnaGammaJet.h
@@ -7,6 +7,9 @@
 /* History of cvs commits:
  *
  * $Log$
+ * Revision 1.1  2007/01/23 17:17:29  schutz
+ * New Gamma package
+ *
  *
  */
 
@@ -40,9 +43,11 @@ public:
   AliAnaGammaJet(const AliAnaGammaJet & gj) ; // cpy ctor
   virtual ~AliAnaGammaJet() ; //virtual dtor
   virtual void Exec(Option_t * opt = "") ;
-  virtual void Init(Option_t * opt = "");
+  virtual void ConnectInputData(Option_t *);
+  virtual void CreateOutputObjects();
   virtual void Terminate(Option_t * opt = "");
 
+  void InitParameters();
   Bool_t     AreJetOnlyInCTS() const {return fJetsOnlyInCTS ; } 
   Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
   Double_t GetEtaEMCALCut() const {return fEtaEMCALCut;}

diff --git a/PWG4/Makefile b/PWG4/Makefile
index fc2a10e7275326290c130f03842b2efeae50b08b..48fa25bca8c562e95599609f74a2fafccfae0cbb 100644 (file)
--- a/PWG4/Makefile
+++ b/PWG4/Makefile
@@ -1,5 +1,5 @@
 
-include $(ROOTSYS)/test/Makefile.arch
+include Makefile.arch
 
 default-target: libPWG4Gamma.so