New Gamma package
authorschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Jan 2007 17:17:29 +0000 (17:17 +0000)
committerschutz <schutz@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 23 Jan 2007 17:17:29 +0000 (17:17 +0000)
ESDCheck/PROOF-INF.AnalysisCheck/SETUP.C
PWG4/AliAnaGammaDirect.cxx [new file with mode: 0644]
PWG4/AliAnaGammaDirect.h [new file with mode: 0644]
PWG4/AliAnaGammaHadron.cxx [new file with mode: 0644]
PWG4/AliAnaGammaHadron.h [new file with mode: 0644]
PWG4/AliAnaGammaJet.cxx [new file with mode: 0644]
PWG4/AliAnaGammaJet.h [new file with mode: 0644]
PWG4/GammaLinkDef.h [new file with mode: 0644]
PWG4/Makefile [new file with mode: 0644]
PWG4/libGamma.pkg [new file with mode: 0644]

index d63773c..028db33 100644 (file)
@@ -2,12 +2,12 @@ void SETUP()
 {
 
    // Load the ESD library
-   gSystem->Load("libAnalysisCheck");
+   gSystem->Load("libGamma");
 
    // Set the Include paths
-   gSystem->SetIncludePath("-I$ROOTSYS/include -IAnalysisCheck");
-   gROOT->ProcessLine(".include AnalysisCheck");
+   gSystem->SetIncludePath("-I$ROOTSYS/include -IPWG4");
+   gROOT->ProcessLine(".include PWG4");
 
    // Set our location, so that other packages can find us
-   gSystem->Setenv("AnalysisCheck_INCLUDE", "AnalysisCheck");
+   gSystem->Setenv("Gamma_INCLUDE", "PWG4");
 }
diff --git a/PWG4/AliAnaGammaDirect.cxx b/PWG4/AliAnaGammaDirect.cxx
new file mode 100644 (file)
index 0000000..7c6ee3b
--- /dev/null
@@ -0,0 +1,642 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma correlations (gamma-jet, 
+// gamma-hadron and isolation cut.
+// This class contains 3 main methods: one to fill lists of particles (ESDs) comming 
+//  from the CTS (ITS+TPC) and the calorimeters;  the other one tags a candidate 
+//  cluster as isolated;  the last one search in the 
+//  corresponing calorimeter for the highest energy cluster, identified it as 
+//  prompt photon;
+//
+//  Class created from old AliPHOSGammaJet 
+//  (see AliRoot versions previous Release 4-09)
+//
+//*-- Author: Gustavo Conesa (LNF-INFN) 
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include <TFile.h>
+#include <TParticle.h>
+#include <TH2.h>
+
+#include "AliAnaGammaDirect.h" 
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliESDCaloCluster.h"
+#include "Riostream.h"
+#include "AliLog.h"
+
+ClassImp(AliAnaGammaDirect)
+
+//____________________________________________________________________________
+  AliAnaGammaDirect::AliAnaGammaDirect(const char *name) : 
+    AliAnalysisTask(name,""), fChain(0), fESD(0),
+    fOutputContainer(new TObjArray(100)), 
+    fPrintInfo(0), fMinGammaPt(0.),
+    fCalorimeter(""), fEMCALPID(0),fPHOSPID(0),
+    fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
+    fNCones(0),fNPtThres(0),fMakeICMethod(0)
+{
+  //Ctor        
+  TList * list = gDirectory->GetListOfKeys() ; 
+  TIter next(list) ; 
+  TH2F * h = 0 ;
+  Int_t index ; 
+  for (index = 0 ; index < list->GetSize()-1 ; index++) { 
+    //-1 to avoid GammaJet Task
+    h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; 
+    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
+  DefineOutput(0,  TObjArray::Class()) ; 
+  
+}
+
+
+//____________________________________________________________________________
+AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) : 
+  AliAnalysisTask(g), fChain(g.fChain), fESD(g.fESD),
+  fOutputContainer(g. fOutputContainer),  fPrintInfo(g.fPrintInfo),
+  fMinGammaPt(g.fMinGammaPt), fCalorimeter(g.fCalorimeter),
+  fEMCALPID(g.fEMCALPID),fPHOSPID(g.fPHOSPID),
+  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];
+  }
+}
+
+//____________________________________________________________________________
+AliAnaGammaDirect::~AliAnaGammaDirect() 
+{
+  // Remove all pointers
+  fOutputContainer->Clear() ; 
+  delete fOutputContainer ;
+  delete fhNGamma    ;  
+  delete fhPhiGamma  ; 
+  delete fhEtaGamma   ;  
+  delete fhPtCandidate ;
+  delete [] fhPtThresIsolated ;
+  delete [] fhPtSumIsolated ;
+
+}
+
+//____________________________________________________________________________
+void AliAnaGammaDirect::CreateParticleList(TClonesArray * pl, 
+                                       TClonesArray * plCTS, 
+                                       TClonesArray * plEMCAL,  
+                                       TClonesArray * plPHOS){
+  
+  //Create a list of particles from the ESD. These particles have been measured 
+  //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 
+  
+  Int_t index = pl->GetEntries() ; 
+  Int_t npar  = 0 ;
+  Float_t *pid = new Float_t[AliPID::kSPECIESN];  
+  AliDebug(3,"Fill particle lists");
+  if(fPrintInfo)
+    AliInfo(Form("fCalorimeter %s",fCalorimeter.Data()));
+
+  Double_t v[3] ; //vertex ;
+  fESD->GetVertex()->GetXYZ(v) ; 
+
+  //########### PHOS ##############
+  if(fCalorimeter == "PHOS"){
+
+    Int_t begphos = fESD->GetFirstPHOSCluster();  
+    Int_t endphos = fESD->GetFirstPHOSCluster() + 
+      fESD->GetNumberOfPHOSClusters() ;  
+    Int_t indexNePHOS = plPHOS->GetEntries() ;
+    AliDebug(3,Form("First PHOS particle %d, last particle %d", begphos,endphos));
+
+    if(fCalorimeter == "PHOS"){
+      for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
+       AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
+
+       //Create a TParticle to fill the particle list
+       
+       Float_t en = clus->GetClusterEnergy() ;
+       Float_t *p = new Float_t();
+       clus->GetGlobalPosition(p) ;
+       TVector3 pos(p[0],p[1],p[2]) ; 
+       Double_t phi  = pos.Phi();
+       Double_t theta= pos.Theta();
+       Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
+       Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
+       Double_t pz = en*TMath::Cos(theta);
+
+       TParticle * particle = new TParticle() ;
+       particle->SetMomentum(px,py,pz,en) ;
+       AliDebug(4,Form("PHOS clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+       
+       //Select only photons
+       
+       pid=clus->GetPid();
+       //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
+       if( !fPHOSPID)
+         new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
+       else if( pid[AliPID::kPhoton] > 0.75)
+         new((*plPHOS)[indexNePHOS++])   TParticle(*particle) ;
+      }
+    }
+  }
+
+  //########### CTS (TPC+ITS) #####################
+  Int_t begtpc   = 0 ;  
+  Int_t endtpc   = fESD->GetNumberOfTracks() ;
+  Int_t indexCh  = plCTS->GetEntries() ;
+  AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
+  
+  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
+    AliESDtrack * track = fESD->GetTrack(npar) ; // retrieve track from esd
+    //We want tracks fitted in the detectors:
+    ULong_t status=AliESDtrack::kTPCrefit;
+    status|=AliESDtrack::kITSrefit;
+   
+    //We want tracks whose PID bit is set:
+    //     ULong_t status =AliESDtrack::kITSpid;
+    //     status|=AliESDtrack::kTPCpid;
+  
+    if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
+      // Do something with the tracks which were successfully
+      // re-fitted 
+      Double_t en = 0; //track ->GetTPCsignal() ;
+      Double_t mom[3];
+      track->GetPxPyPz(mom) ;
+      Double_t px = mom[0];
+      Double_t py = mom[1];
+      Double_t pz = mom[2]; //Check with TPC people if this is correct.
+      Int_t pdg = 11; //Give any charged PDG code, in this case electron.
+      //I just want to tag the particle as charged
+       TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                                px, py, pz, en, v[0], v[1], v[2], 0);
+  
+      //TParticle * particle = new TParticle() ;
+      //particle->SetMomentum(px,py,pz,en) ;
+      new((*plCTS)[indexCh++])       TParticle(*particle) ;    
+      new((*pl)[index++])           TParticle(*particle) ;
+    }
+  }
+  
+  //################ EMCAL ##############
+  
+  Int_t begem = fESD->GetFirstEMCALCluster();  
+  Int_t endem = fESD->GetFirstEMCALCluster() + 
+    fESD->GetNumberOfEMCALClusters() ;  
+  Int_t indexNe  = plEMCAL->GetEntries() ; 
+  
+  AliDebug(3,Form("First EMCAL particle %d, last particle %d",begem,endem));
+    
+    for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
+      AliESDCaloCluster * clus = fESD->GetCaloCluster(npar) ; // retrieve track from esd
+      Int_t clustertype= clus->GetClusterType();
+      if(clustertype == AliESDCaloCluster::kClusterv1){
+       Float_t en = clus->GetClusterEnergy() ;
+       Float_t *p = new Float_t();
+       clus->GetGlobalPosition(p) ;
+       TVector3 pos(p[0],p[1],p[2]) ;
+       Double_t phi  = pos.Phi();
+       Double_t theta= pos.Theta();
+       Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
+       Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
+       Double_t pz = en*TMath::Cos(theta);
+
+       pid=clus->GetPid();
+       if(fCalorimeter == "EMCAL")
+         {
+           TParticle * particle = new TParticle() ;
+           particle->SetMomentum(px,py,pz,en) ;
+           AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+           if(!fEMCALPID) //Only identified particles
+             new((*plEMCAL)[indexNe++])       TParticle(*particle) ;
+           else if(pid[AliPID::kPhoton] > 0.75)
+             new((*plEMCAL)[indexNe++])       TParticle(*particle) ;       
+         }
+       else
+         {
+           Int_t pdg = 0;
+           if(fEMCALPID) 
+             {
+               if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
+                 pdg = 22;
+               else if( pid[AliPID::kPi0] > 0.75)
+                 pdg = 111;
+             }
+           else
+             pdg = 22; //No PID, assume all photons
+           
+           TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
+                                                px, py, pz, en, v[0], v[1], v[2], 0);
+           AliDebug(4,Form("EMCAL clusters: pt %f, phi %f, eta %f", particle->Pt(),particle->Phi(),particle->Eta()));
+           
+           new((*plEMCAL)[indexNe++])       TParticle(*particle) ; 
+           new((*pl)[index++])           TParticle(*particle) ;
+         }
+      }
+    }
+
+    AliDebug(3,"Particle lists filled");
+    
+}
+
+
+
+//____________________________________________________________________________
+void AliAnaGammaDirect::Exec(Option_t *) 
+{
+  
+  // Processing of one event
+    
+  //Get ESDs
+  Long64_t entry = fChain->GetReadEntry() ;
+  
+  if (!fESD) {
+    AliError("fESD is not connected to the input!") ; 
+    return ; 
+  }
+  
+  if (fPrintInfo) 
+     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(fChain))->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 (EMCAL)
+  TClonesArray * plPHOS     = new TClonesArray("TParticle",1000);  // All particles measured in PHOS as Gamma calorimeter
+  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){
+    
+    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
+  
+  AliDebug(2, "End of analysis, delete pointers");
+  
+  particleList->Delete() ; 
+  plCTS->Delete() ;
+  plNe->Delete() ;
+  plPHOS->Delete() ;
+  pLeading->Delete();
+  pGamma->Delete();
+
+  delete plNe ;
+  delete plCTS ;
+  delete particleList ;
+  //  delete pLeading;
+  //  delete pGamma;
+
+  PostData(0, fOutputContainer);
+}    
+
+
+//____________________________________________________________________________
+void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &Is) const 
+{
+  //Search for the prompt photon in Calorimeter with pt > fMinGammaPt
+
+  Double_t pt = 0;
+  Int_t index = -1; 
+  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
+      index = ipr ;
+      pt = particle->Pt();
+      pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
+      AliDebug(4,Form("Cluster in calo: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+      Is  = kTRUE;
+    }
+  }
+
+  //Do Isolation?
+  if(fMakeICMethod && Is)
+    {
+      Float_t coneptsum = 0 ;
+      Bool_t  icPtThres   = kFALSE;
+      Bool_t  icPtSum     = kFALSE;
+      MakeIsolationCut(plCTS,pl, pGamma, index, 
+                      icPtThres, icPtSum,coneptsum);
+      if(fMakeICMethod == 1) //Pt thres method
+       Is = icPtThres ;
+      if(fMakeICMethod == 2) //Pt cone sum method
+       Is = icPtSum ;
+    }
+  
+  if(Is){
+    AliDebug(3,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
+    AliDebug(3,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
+    //Fill prompt gamma histograms
+    fhNGamma  ->Fill(pGamma->Pt());
+    fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
+    fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
+  }
+  else
+    AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
+}
+
+  //____________________________________________________________________________
+void AliAnaGammaDirect::Init(const Option_t * )
+{
+  // 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.
+  fCalorimeter="PHOS";
+  fPrintInfo           = kTRUE;
+  fMinGammaPt  = 10. ;
+
+  //Fill particle lists when PID is ok
+  fEMCALPID = kFALSE;
+  fPHOSPID = kFALSE;
+
+  fConeSize             = 0.2 ; 
+  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() ;
+}
+
+//__________________________________________________________________
+void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
+                                          TClonesArray * plNe, 
+                                          TParticle * pCandidate, 
+                                          Int_t index, 
+                                          Bool_t  &icmpt,  Bool_t  &icms, 
+                                          Float_t &coneptsum) const 
+{  
+  //Search in cone around a candidate particle if it is isolated 
+  Float_t phiC  = pCandidate->Phi() ;
+  Float_t etaC = pCandidate->Eta() ;
+  Float_t pt     = -100. ;
+  Float_t eta   = -100.  ;
+  Float_t phi    = -100.  ;
+  Float_t rad   = -100 ;
+  Int_t    n        = 0 ;
+  TParticle * particle  = new TParticle();
+
+  coneptsum = 0; 
+  icmpt = kFALSE;
+  icms  = kFALSE;
+
+  //Check charged particles in cone.
+  for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
+    particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
+    pt    = particle->Pt();
+    eta  = particle->Eta();
+    phi  = particle->Phi() ;
+    
+    //Check if there is any particle inside cone with pt larger than  fPtThreshold
+    rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
+                     TMath::Power((phi-phiC),2));
+    if(rad<fConeSize){
+      AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
+      coneptsum+=pt;
+      if(pt > fPtThreshold ) n++;
+    }
+  }// charged particle loop
+  
+  //Check neutral particles in cone.
+  for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
+    if(ipr != index){//Do not count the candidate
+      particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
+      pt    = particle->Pt();
+      eta  = particle->Eta();
+      phi  = particle->Phi() ;
+      
+      //Check if there is any particle inside cone with pt larger than  fPtThreshold
+      rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
+                       TMath::Power((phi-phiC),2));
+      if(rad<fConeSize){
+       AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
+       coneptsum+=pt;
+       if(pt > fPtThreshold ) n++;
+      }
+    }
+  }// neutral particle loop
+  
+  if(n == 0) 
+    icmpt =  kTRUE ;
+  if(coneptsum < fPtSumThreshold)
+    icms  =  kTRUE ;
+
+}
+
+//___________________________________________________________________
+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
+{
+
+  //Print some relevant parameters set for the analysis
+  if(! opt)
+    return;
+
+  Info("Print", "%s %s", GetName(), GetTitle() ) ;
+  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("Calorimeter            =     %s\n", fCalorimeter.Data()) ; 
+} 
+
+void AliAnaGammaDirect::Terminate(Option_t *)
+{
+   // The Terminate() function is the last function to be called during
+   // a query. It always runs on the client, it can be used to present
+   // the results graphically or save the results to file.
+    
+
+}
diff --git a/PWG4/AliAnaGammaDirect.h b/PWG4/AliAnaGammaDirect.h
new file mode 100644 (file)
index 0000000..65a983d
--- /dev/null
@@ -0,0 +1,132 @@
+#ifndef ALIANAGAMMADIRECT_H
+#define ALIANAGAMMADIRECT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+
+// 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 
+//  corresponing calorimeter for the highest energy cluster, identify it as 
+//  prompt photon(Shower Shape (done) and Isolation Cut (TODO in new class?)).
+//
+//  Class created from old AliPHOSGammaJet
+//  (see AliRoot versions previous Release 4-09)
+
+//*-- 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 "AliAnalysisTask.h" 
+
+class AliESD ; 
+class AliAnaGammaDirect : public AliAnalysisTask {
+
+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 Terminate(Option_t * opt = "");
+  
+  TTree *     GetChain()                const {return fChain ; }
+  AliESD *    GetESD()                  const {return fESD ; }
+  TObjArray * GetOutputContainer()      const {return fOutputContainer ; }
+  Double_t  GetMinGammaPt()    const {return fMinGammaPt ; }
+  TString    GetCalorimeter()       const {return fCalorimeter ; }
+  Bool_t      GetPrintInfo()           const {return fPrintInfo ; }
+  Float_t     GetConeSize()          const {return fConeSize ; }
+  Float_t     GetPtThreshold()      const {return fPtThreshold ; }
+  Float_t     GetPtSumThres()     const {return fPtSumThreshold ; }
+  Int_t        GetICMethod()          const {return fMakeICMethod ; }
+
+  Bool_t   IsEMCALPIDOn() const {return fEMCALPID ; }
+  Bool_t   IsPHOSPIDOn() const {return fPHOSPID ; }
+
+  void Print(const Option_t * opt)const;
+
+  void SetMinGammaPt(Double_t ptcut){fMinGammaPt =ptcut;}
+  void SetCalorimeter(TString calo){ fCalorimeter= calo ; }
+  void SetPrintInfo(Bool_t print){ fPrintInfo = print ; }
+  void SetConeSize(Float_t r)              {fConeSize = r ; }
+  void SetPtThreshold(Float_t pt)        {fPtThreshold = pt; };
+  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 ; }
+
+  void CreateParticleList(TClonesArray * particleList, 
+                         TClonesArray * plCh, TClonesArray * plNe, 
+                         TClonesArray * plNePHOS);
+  
+  
+  void GetPromptGamma(TClonesArray * plNe, TClonesArray * plCTS, TParticle * pGamma, Bool_t &Is)  const;
+  
+  void MakeIsolationCut(TClonesArray * plCTS, TClonesArray * plNe, 
+                       TParticle *pCandidate, Int_t index, 
+                       Bool_t &imcpt, Bool_t &icms, Float_t &ptsum) const ;
+  
+  void MakeHistos() ;
+  
+  
+ private:
+
+  TTree       *fChain ;   //!pointer to the analyzed TTree or TChain
+  AliESD       *fESD ;     //! Declaration of leave types
+  TObjArray  *fOutputContainer ; //! output data container
+  Bool_t        fPrintInfo ;      //Print most interesting information on screen
+  Double_t    fMinGammaPt ;  // Min pt in Calorimeter
+  TString      fCalorimeter ; //PHOS or EMCAL detects Gamma
+  Bool_t       fEMCALPID ;//Fill EMCAL particle lists with particles with corresponding pid
+  Bool_t       fPHOSPID;  //Fill PHOS particle lists with particles with corresponding pid
+  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)
+} ;
+
+#endif //ALIANAGAMMADIRECT_H
+
+
+
diff --git a/PWG4/AliAnaGammaHadron.cxx b/PWG4/AliAnaGammaHadron.cxx
new file mode 100644 (file)
index 0000000..8b0049d
--- /dev/null
@@ -0,0 +1,600 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations 
+//  Basically it seaches for a prompt photon in the Calorimeters acceptance, 
+//  if so we construct a jet around the highest pt particle in the opposite 
+//  side in azimuth, inside the Central Tracking System (ITS+TPC) and 
+//  EMCAL acceptances (only when PHOS detects the gamma). First the leading 
+//  particle and then the jet have to fullfill several conditions 
+//  (energy, direction ..) to be accepted. Then the fragmentation function 
+//  of this jet is constructed   
+//  Class created from old AliPHOSGammaPion 
+//  (see AliRoot versions previous Release 4-09)
+//
+//*-- Author: Gustavo Conesa (LNF-INFN) 
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include <TFile.h>
+#include <TParticle.h>
+#include <TH2.h>
+
+#include "AliAnaGammaHadron.h" 
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliESDCaloCluster.h"
+#include "Riostream.h"
+#include "AliLog.h"
+
+ClassImp(AliAnaGammaHadron)
+
+//____________________________________________________________________________
+AliAnaGammaHadron::AliAnaGammaHadron(const char *name) : 
+  AliAnaGammaDirect(name), 
+  fPhiMaxCut(0.), fPhiMinCut(0.), 
+  fInvMassMaxCut(0.), fInvMassMinCut(0.),
+  fMinPtPion(0),
+  fAngleMaxParam()
+{
+
+  // ctor
+  fAngleMaxParam.Set(4) ;
+  fAngleMaxParam.Reset(0.);
+        
+  TList * list = gDirectory->GetListOfKeys() ; 
+  TIter next(list) ; 
+  TH2F * h = 0 ;
+  Int_t index ; 
+  for (index = 0 ; index < list->GetSize()-1 ; index++) { 
+    //-1 to avoid GammaPion Task
+    h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; 
+    fOutputContainer->Add(h) ; 
+  }
+
+  // Input slot #0 works with an Ntuple
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  DefineOutput(0,  TObjArray::Class()) ; 
+  
+}
+
+
+//____________________________________________________________________________
+AliAnaGammaHadron::AliAnaGammaHadron(const AliAnaGammaHadron & gj) : 
+  AliAnaGammaDirect(gj), 
+  fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut), 
+  fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
+  fMinPtPion(gj.fMinPtPion),
+  fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam)
+
+{
+  // cpy ctor
+  SetName (gj.GetName()) ; 
+  SetTitle(gj.GetTitle()) ; 
+
+}
+
+//____________________________________________________________________________
+AliAnaGammaHadron::~AliAnaGammaHadron() 
+{
+  // Remove all pointers
+  fOutputContainer->Clear() ; 
+  delete fOutputContainer ;
+  
+  delete fhPhiCharged  ;  
+  delete fhPhiNeutral   ; 
+  delete fhEtaCharged  ; 
+  delete fhEtaNeutral  ; 
+  delete fhDeltaPhiGammaCharged  ;  
+  delete fhDeltaPhiGammaNeutral   ; 
+  delete fhDeltaEtaGammaCharged  ; 
+  delete fhDeltaEtaGammaNeutral  ; 
+
+  delete fhCorrelationGammaNeutral  ; 
+  delete fhCorrelationGammaCharged  ; 
+
+  delete fhAnglePairNoCut  ; 
+  delete fhAnglePairAzimuthCut  ; 
+  delete fhAnglePairOpeningAngleCut   ; 
+  delete fhAnglePairAllCut   ;  
+  delete fhInvMassPairNoCut    ; 
+  delete fhInvMassPairAzimuthCut  ; 
+  delete fhInvMassPairOpeningAngleCut  ; 
+  delete fhInvMassPairAllCut   ;    
+
+}
+
+
+
+//____________________________________________________________________________
+void AliAnaGammaHadron::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
+
+  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
+  //Fill particle lists 
+  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, "Make correlation");
+    
+    //Search correlation 
+    MakeGammaChargedCorrelation(plCTS, pGamma);
+    MakeGammaNeutralCorrelation(plNe, pGamma);
+
+  }//Gamma in Calo
+     
+  AliDebug(2, "End of analysis, delete pointers");
+
+  particleList->Delete() ; 
+  plCTS->Delete() ;
+  plNe->Delete() ;
+  plCalo->Delete() ;
+  pGamma->Delete();
+
+  delete plNe ;
+  delete plCalo ;
+  delete plCTS ;
+  delete particleList ;
+  //  delete pGamma;
+
+  PostData(0, fOutputContainer);
+}    
+
+
+//____________________________________________________________________________
+void  AliAnaGammaHadron::MakeGammaChargedCorrelation(TClonesArray * pl, TParticle * pGamma) const 
+{  
+  //Search for the charged particle with highest with 
+  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
+  Double_t ptg  = pGamma->Pt();
+  Double_t phig = pGamma->Phi();
+  Double_t pt    = -100.;
+  Double_t rat   = -100.; 
+  Double_t phi   = -100. ;
+
+  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
+    
+    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
+
+    pt    = particle->Pt();
+    rat   = pt/ptg ;
+    phi   = particle->Phi() ;
+    
+    AliDebug(3,Form("pt %f, phi %f, phi gamma %f. Cuts:  delta phi min %f,  max%f, pT min %f",pt,phi,phig,fPhiMinCut,fPhiMaxCut,fMinPtPion));
+    
+    fhEtaCharged->Fill(ptg,particle->Eta());
+    fhPhiCharged->Fill(ptg,phi);
+    fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
+    fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
+    //Selection within angular and energy limits
+    if(((phig-phi)> fPhiMinCut) && ((phig-phi)<fPhiMaxCut) && pt > fMinPtPion){
+      AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
+      fhCorrelationGammaCharged->Fill(ptg,rat);
+    } 
+  }//particle loop
+}
+
+//____________________________________________________________________________
+void  AliAnaGammaHadron::MakeGammaNeutralCorrelation(TClonesArray * pl, TParticle * pGamma)  
+{  
+
+  //Search for the neutral pion with highest with 
+  //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
+  Double_t pt = -100.;
+  Double_t rat = -100.; 
+  Double_t phi = -100. ;
+  Double_t ptg  = pGamma->Pt();
+  Double_t phig = pGamma->Phi();
+
+  TIter next(pl);
+  TParticle * particle = 0;
+  
+  Int_t iPrimary = -1;
+  TLorentzVector gammai,gammaj;
+  Double_t angle = 0., e = 0., invmass = 0.;
+  Int_t ksPdg = 0;
+  Int_t jPrimary=-1;
+
+  while ( (particle = (TParticle*)next()) ) {
+    iPrimary++;          
+    ksPdg = particle->GetPdgCode();
+
+    //2 gamma overlapped, found with PID
+    if(ksPdg == 111){ 
+      pt  = particle->Pt();
+      rat = pt/ptg ;
+      phi = particle->Phi() ;
+      fhEtaCharged->Fill(ptg,particle->Eta());
+      fhPhiCharged->Fill(ptg,phi);
+      fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-particle->Eta());
+      fhDeltaPhiGammaCharged->Fill(ptg,phig-phi);
+      //AliDebug(3,Form("pt %f, phi %f",pt,phi));
+      if (GetPrintInfo())
+       AliInfo(Form("pt %f, phi %f",pt,phi));
+      //Selection within angular and energy limits
+      if((pt> ptg)&& ((phig-phi)>fPhiMinCut)&&((phig-phi)<fPhiMaxCut)){
+       fhCorrelationGammaNeutral ->Fill(ptg,rat);
+       //AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
+       if (GetPrintInfo())
+         AliInfo(Form("Selected: pt %f, phi %f",pt,phi));
+      }// cuts
+    }// pdg = 111
+
+    //Make invariant mass analysis
+    else if(ksPdg == 22){
+      //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());
+           pt  = (gammai+gammaj).Pt();
+           phi = (gammai+gammaj).Phi();
+           if(phi < 0)
+             phi+=TMath::TwoPi();
+           rat          = pt/ptg ;
+           invmass = (gammai+gammaj).M();
+           angle      = gammaj.Angle(gammai.Vect());
+           e             = (gammai+gammaj).E();
+           fhEtaNeutral->Fill(ptg,(gammai+gammaj).Eta());
+           fhPhiNeutral->Fill(ptg,phi);
+           fhDeltaEtaGammaNeutral->Fill(ptg,pGamma->Eta()-(gammai+gammaj).Eta());
+           fhDeltaPhiGammaNeutral->Fill(ptg,phig-phi);
+           // AliDebug(3,Form("pt %f, phi %f",pt,phi));
+           if (GetPrintInfo())
+             AliInfo(Form("pt %f, phi %f",pt,phi));
+
+           //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((phig-phi) > fPhiMinCut && (phig-phi) < fPhiMaxCut 
+              && pt > fMinPtPion){
+             fhAnglePairAzimuthCut     ->Fill(e,angle);
+             fhInvMassPairAzimuthCut->Fill(ptg,invmass);
+             AliDebug(3,Form("1st cut: pt %f, phi %f",pt,phi));
+
+             //Second cut on the aperture of the pair
+             if(IsAngleInWindow(angle,e)){
+               fhAnglePairOpeningAngleCut     ->Fill(e,angle);
+               fhInvMassPairOpeningAngleCut->Fill(ptg,invmass);
+                AliDebug(3,Form("2nd cut: pt %f, phi %f",pt,phi));
+
+               //Third cut on the invariant mass of the pair
+               if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
+                 fhInvMassPairAllCut  ->Fill(ptg,invmass);
+                 fhAnglePairAllCut       ->Fill(e,angle);
+                 //Fill correlation histogram
+                 fhCorrelationGammaNeutral ->Fill(ptg,rat);
+                 //AliDebug(2,Form("Selected: pt %f, phi %f",pt,phi));
+                 if (GetPrintInfo())
+                   AliInfo(Form("Selected: pt %f, phi %f",pt,phi));
+               }//(invmass>0.125) && (invmass<0.145)
+             }//Opening angle cut
+           }//Azimuth and pt cut.
+         }//if pdg = 22
+       }//iprimary<jprimary
+      }//while
+    }// if pdg = 22
+  }//while
+}
+
+  //____________________________________________________________________________
+void AliAnaGammaHadron::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) ;
+
+  //fPrintInfo           = kTRUE;
+  fInvMassMaxCut  = 0.16 ;
+  fInvMassMinCut  = 0.11 ;
+  fPhiMaxCut      = 4.5;
+  fPhiMinCut      = 1.5 ;
+
+  fMinPtPion = 0.   ;
+
+  //Fill particle lists when PID is ok
+  // fEMCALPID = kFALSE;
+  // fPHOSPID = kFALSE;
+
+  //Initialization of histograms 
+
+  MakeHistos() ; 
+
+}
+
+//__________________________________________________________________________-
+Bool_t AliAnaGammaHadron::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;
+}
+
+//____________________________________________________________________________
+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
+{
+
+  //Print some relevant parameters set for the analysis
+  if(! opt)
+    return;
+
+  Info("Print", "%s %s", GetName(), GetTitle() ) ;
+  printf("pT Pion       >    %f\n", fMinPtPion) ; 
+  printf("Phi Pion      <     %f\n", fPhiMaxCut) ; 
+  printf("Phi Pion      >     %f\n", fPhiMinCut) ;
+  printf("M_pair        <     %f\n", fInvMassMaxCut) ; 
+  printf("M_pair        >     %f\n", fInvMassMinCut) ; 
+} 
+
+//__________________________________________
+void AliAnaGammaHadron::Terminate(Option_t *)
+{
+   // The Terminate() function is the last function to be called during
+   // a query. It always runs on the client, it can be used to present
+   // the results graphically or save the results to file.
+    
+
+}
diff --git a/PWG4/AliAnaGammaHadron.h b/PWG4/AliAnaGammaHadron.h
new file mode 100644 (file)
index 0000000..cf33471
--- /dev/null
@@ -0,0 +1,120 @@
+#ifndef ALIANAGAMMAHADRON_H
+#define ALIANAGAMMAHADRON_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+//  Class for the analysis of gamma-jet correlations.     
+//  Basically it seaches for a prompt photon in the Calorimeters acceptance, 
+//  if so we construct a jet around the highest pt particle in the opposite 
+//  side in azimuth. This jet has to fullfill several conditions to be 
+//  accepted. Then the fragmentation function of this jet is constructed 
+//  Class created from old AliPHOSGammaPion
+
+//*-- 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 "AliAnaGammaDirect.h" 
+
+class AliESD ; 
+class AliAnaGammaHadron : public AliAnaGammaDirect {
+
+public: 
+
+  AliAnaGammaHadron(const char *name) ; // default ctor
+  AliAnaGammaHadron(const AliAnaGammaHadron & gj) ; // cpy ctor
+  virtual ~AliAnaGammaHadron() ; //virtual dtor
+  virtual void Exec(Option_t * opt = "") ;
+  virtual void Init(Option_t * opt = "");
+  virtual void Terminate(Option_t * opt = "");
+  Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
+  Double_t GetInvMassMaxCut() const {return fInvMassMaxCut ; }
+  Double_t GetInvMassMinCut() const {return fInvMassMinCut ; }
+  Double_t GetPhiMaxCut() const {return fPhiMaxCut ; }
+  Double_t GetPhiMinCut() const {return fPhiMinCut ; }
+  Float_t    GetMinPtPion() const {return fMinPtPion ; }
+
+  void Print(const Option_t * opt)const;
+  
+  void SetAngleMaxParam(Int_t i, Double_t par)
+  {fAngleMaxParam.AddAt(par,i) ; }
+  void SetMinPtPion(Float_t pt){fMinPtPion = pt; };
+  void SetInvMassCutRange(Double_t invmassmin, Double_t invmassmax)
+    {fInvMassMaxCut =invmassmax;  fInvMassMinCut =invmassmin;} 
+  void SetPhiCutRange(Double_t phimin, Double_t phimax)
+  {fPhiMaxCut =phimax;  fPhiMinCut =phimin;}
+
+  private:
+
+  Bool_t IsAngleInWindow(const Float_t angle, const Float_t e);
+  void MakeGammaChargedCorrelation(TClonesArray * pl, TParticle *pGamma) const ;
+  void MakeGammaNeutralCorrelation(TClonesArray * pl, TParticle *pGamma)  ;
+  void MakeHistos() ;
+
+ private:
+
+  //  TTree       *fChain ;   //!pointer to the analyzed TTree or TChain
+  //AliESD       *fESD ;     //! Declaration of leave types
+
+  Double_t   fPhiMaxCut ;      // 
+  Double_t   fPhiMinCut ;      // 
+  // Double_t   fGammaPtCut ;  // Min pt in Calorimeter
+  Double_t   fInvMassMaxCut ;  // Invariant Mass cut maximum
+  Double_t   fInvMassMinCut ;  // Invariant Masscut minimun
+  Double_t   fMinPtPion;       // Minimum pt of pion
+  TObjArray  *fOutputContainer ; //! output data container
+  TString     fNamePtThres[10];   // String name of pt th to append to histos
+  TArrayD    fAngleMaxParam ; //Max opening angle selection parameters
+  //TString     fCalorimeter ; //PHOS or EMCAL detects Gamma
+  //Bool_t      fEMCALPID ; //Fill EMCAL particle lists with particles with corresponding pid
+  //Bool_t      fPHOSPID;   //Fill PHOS particle lists with particles with corresponding pid
+
+  //Histograms
+
+  TH2F * fhPhiCharged  ; 
+  TH2F * fhPhiNeutral   ; 
+  TH2F * fhEtaCharged  ; 
+  TH2F * fhEtaNeutral   ; 
+  TH2F * fhDeltaPhiGammaCharged  ;  
+  TH2F * fhDeltaPhiGammaNeutral   ; 
+  TH2F * fhDeltaEtaGammaCharged  ; 
+  TH2F * fhDeltaEtaGammaNeutral  ; 
+
+  TH2F * fhCorrelationGammaNeutral  ; 
+  TH2F * fhCorrelationGammaCharged  ; 
+
+  TH2F * fhAnglePairAccepted  ; 
+  TH2F * fhAnglePairNoCut  ; 
+  TH2F * fhAnglePairAzimuthCut  ; 
+  TH2F * fhAnglePairOpeningAngleCut   ; 
+  TH2F * fhAnglePairAllCut   ; 
+  TH2F * fhInvMassPairNoCut    ; 
+  TH2F * fhInvMassPairAzimuthCut  ; 
+  TH2F * fhInvMassPairOpeningAngleCut  ; 
+  TH2F * fhInvMassPairAllCut   ; 
+  
+  ClassDef(AliAnaGammaHadron,0)
+} ;
+
+#endif //ALIANAGAMMAHADRON_H
+
+
+
diff --git a/PWG4/AliAnaGammaJet.cxx b/PWG4/AliAnaGammaJet.cxx
new file mode 100644 (file)
index 0000000..66a4c3a
--- /dev/null
@@ -0,0 +1,1354 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+// Class for the analysis of gamma-jet correlations 
+//  Basically it seaches for a prompt photon in the Calorimeters acceptance, 
+//  if so we construct a jet around the highest pt particle in the opposite 
+//  side in azimuth, inside the Central Tracking System (ITS+TPC) and 
+//  EMCAL acceptances (only when PHOS detects the gamma). First the leading 
+//  particle and then the jet have to fullfill several conditions 
+//  (energy, direction ..) to be accepted. Then the fragmentation function 
+//  of this jet is constructed   
+//  Class created from old AliPHOSGammaJet 
+//  (see AliRoot versions previous Release 4-09)
+//
+//*-- Author: Gustavo Conesa (LNF-INFN) 
+//////////////////////////////////////////////////////////////////////////////
+
+
+// --- ROOT system ---
+
+#include <TFile.h>
+#include <TParticle.h>
+#include <TH2.h>
+
+#include "AliAnaGammaJet.h" 
+#include "AliESD.h"
+#include "AliESDtrack.h"
+#include "AliESDCaloCluster.h"
+#include "Riostream.h"
+#include "AliLog.h"
+
+ClassImp(AliAnaGammaJet)
+
+//____________________________________________________________________________
+AliAnaGammaJet::AliAnaGammaJet(const char *name) : 
+  AliAnaGammaDirect(name), 
+  fSeveralConeAndPtCuts(0),
+  fPbPb(kFALSE), 
+  fJetsOnlyInCTS(0),
+  fEtaEMCALCut(0.),fPhiMaxCut(0.),
+  fPhiMinCut(0.), 
+  fInvMassMaxCut(0.), fInvMassMinCut(0.),
+  fJetCTSRatioMaxCut(0.),
+  fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
+  fJetRatioMinCut(0.), fNCone(0),
+  fNPt(0), fCone(0), fPtThreshold(0),
+  fPtJetSelectionCut(0.0),
+  fAngleMaxParam(), fSelect(0)
+{
+
+  // ctor
+  fAngleMaxParam.Set(4) ;
+  fAngleMaxParam.Reset(0.);
+
+  for(Int_t i = 0; i<10; i++){
+    fCones[i]         = 0.0 ;
+    fNameCones[i]     = ""  ;
+    fPtThres[i]      = 0.0 ;
+    fNamePtThres[i]  = ""  ;
+    if( i < 6 ){
+      fJetXMin1[i]     = 0.0 ;
+      fJetXMin2[i]     = 0.0 ;
+      fJetXMax1[i]     = 0.0 ;
+      fJetXMax2[i]     = 0.0 ;
+      fBkgMean[i]      = 0.0 ;
+      fBkgRMS[i]       = 0.0 ;
+      if( i < 2 ){
+       fJetE1[i]        = 0.0 ;
+       fJetE2[i]        = 0.0 ;
+       fJetSigma1[i]    = 0.0 ;
+       fJetSigma2[i]    = 0.0 ;
+       fPhiEMCALCut[i]  = 0.0 ;
+      }
+    }
+  }
+        
+  TList * list = gDirectory->GetListOfKeys() ; 
+  TIter next(list) ; 
+  TH2F * h = 0 ;
+  Int_t index ; 
+  for (index = 0 ; index < list->GetSize()-1 ; index++) { 
+    //-1 to avoid GammaJet Task
+    h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; 
+    fOutputContainer->Add(h) ; 
+  }
+
+  // Input slot #0 works with an Ntuple
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  DefineOutput(0,  TObjArray::Class()) ; 
+  
+}
+
+
+//____________________________________________________________________________
+AliAnaGammaJet::AliAnaGammaJet(const AliAnaGammaJet & gj) : 
+  AliAnaGammaDirect(gj), 
+  fSeveralConeAndPtCuts(gj.fSeveralConeAndPtCuts), 
+  fPbPb(gj.fPbPb), fJetsOnlyInCTS(gj.fJetsOnlyInCTS),
+  fEtaEMCALCut(gj.fEtaEMCALCut),
+  fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut), 
+  fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
+  fRatioMinCut(gj.fRatioMinCut), 
+  fJetCTSRatioMaxCut(gj.fJetCTSRatioMaxCut),
+  fJetCTSRatioMinCut(gj.fJetCTSRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
+  fJetRatioMinCut(gj.fJetRatioMinCut),  fNCone(gj.fNCone),
+  fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
+  fPtJetSelectionCut(gj.fPtJetSelectionCut),
+  fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam), 
+  fSelect(gj.fSelect)
+{
+  // cpy ctor
+  SetName (gj.GetName()) ; 
+  SetTitle(gj.GetTitle()) ; 
+
+  for(Int_t i = 0; i<10; i++){
+    fCones[i]        = gj.fCones[i] ;
+    fNameCones[i]    = gj.fNameCones[i] ;
+    fPtThres[i]      = gj.fPtThres[i] ;
+    fNamePtThres[i]  = gj.fNamePtThres[i] ;
+    if( i < 6 ){
+      fJetXMin1[i]       = gj.fJetXMin1[i] ;
+      fJetXMin2[i]       = gj.fJetXMin2[i] ;
+      fJetXMax1[i]       = gj.fJetXMax1[i] ;
+      fJetXMax2[i]       = gj.fJetXMax2[i] ;
+      fBkgMean[i]        = gj.fBkgMean[i] ;
+      fBkgRMS[i]         = gj.fBkgRMS[i] ;
+      if( i < 2 ){
+       fJetE1[i]        = gj.fJetE1[i] ;
+       fJetE2[i]        = gj.fJetE2[i] ;
+       fJetSigma1[i]    = gj.fJetSigma1[i] ;
+       fJetSigma2[i]    = gj.fJetSigma2[i] ;
+       fPhiEMCALCut[i]  = gj.fPhiEMCALCut[i] ;
+      }
+    }          
+  } 
+}
+
+//____________________________________________________________________________
+AliAnaGammaJet::~AliAnaGammaJet() 
+{
+  // Remove all pointers
+  fOutputContainer->Clear() ; 
+  delete fOutputContainer ;
+  
+  delete fhChargeRatio  ; 
+  delete fhPi0Ratio   ; 
+  delete fhDeltaPhiCharge  ;  
+  delete fhDeltaPhiPi0   ; 
+  delete fhDeltaEtaCharge  ; 
+  delete fhDeltaEtaPi0  ; 
+  delete fhAnglePair  ; 
+  delete fhAnglePairAccepted  ; 
+  delete fhAnglePairNoCut  ; 
+  delete fhAnglePairLeadingCut  ; 
+  delete fhAnglePairAngleCut   ; 
+  delete fhAnglePairAllCut   ; 
+  delete fhAnglePairLeading  ; 
+  delete fhInvMassPairNoCut    ; 
+  delete fhInvMassPairLeadingCut  ; 
+  delete fhInvMassPairAngleCut  ; 
+  delete fhInvMassPairAllCut   ; 
+  delete fhInvMassPairLeading  ; 
+  delete fhNBkg   ; 
+  delete fhNLeading  ; 
+  delete fhNJet  ; 
+  delete fhJetRatio  ; 
+  delete fhJetPt   ; 
+  delete fhBkgRatio   ; 
+  delete fhBkgPt  ; 
+  delete fhJetFragment  ; 
+  delete fhBkgFragment  ; 
+  delete fhJetPtDist  ; 
+  delete fhBkgPtDist  ; 
+  
+  delete [] fhJetRatios;  
+  delete [] fhJetPts;  
+  delete [] fhBkgRatios;
+  delete [] fhBkgPts;  
+
+  delete [] fhNLeadings;
+  delete [] fhNJets;  
+  delete [] fhNBkgs;
+  
+  delete [] fhJetFragments;
+  delete [] fhBkgFragments;
+  delete [] fhJetPtDists;
+  delete [] fhBkgPtDists;
+  
+}
+
+//____________________________________________________________________________
+Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg, 
+                                                const Double_t *par, 
+                                                const Double_t *x) {
+
+  //Parametrized cut for the energy of the jet.
+
+  Double_t epp = par[0] + par[1] * ptg ;
+  Double_t spp = par[2] + par[3] * ptg ;
+  Double_t f   = x[0]   + x[1]   * ptg ;
+  Double_t epb = epp + par[4] ;
+  Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
+  Double_t rat = (epb - spb * f) / ptg ;
+  //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
+  //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
+  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)
+{
+  //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();
+
+    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 
+{  
+  //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::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();
+
+  TIter next(pl);
+  TParticle * particle = 0;
+  
+  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++;          
+    
+    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.
+  
+  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
+
+  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;
+
+  //Initialization of histograms 
+
+  MakeHistos() ; 
+
+}
+
+//__________________________________________________________________________-
+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(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;
+
+}
+
+//____________________________________________________________________________
+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)) ;
+
+  //
+  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) ; 
+  
+  //
+  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) ; 
+    
+    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) ; 
+  }
+  
+  
+  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) ;
+
+  }
+  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
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::MakeJet(TClonesArray * pl, TParticle * pGamma, TParticle* pLeading,TString lastname)
+{
+  //Fill the jet with the particles around the leading particle with 
+  //R=fCone and pt_th = fPtThres. Caclulate the energy of the jet and 
+  //check if we select it. Fill jet histograms
+  
+  TClonesArray * jetList = new TClonesArray("TParticle",1000);
+  TClonesArray * bkgList = new TClonesArray("TParticle",1000);
+
+  TLorentzVector jet   (0,0,0,0);  
+  TLorentzVector bkg(0,0,0,0);
+  TLorentzVector lv (0,0,0,0);
+
+  Double_t ptjet = 0.0;
+  Double_t ptbkg = 0.0;
+  Int_t n0 = 0;
+  Int_t n1 = 0;  
+  Bool_t b1 = kFALSE;
+  Bool_t b0 = kFALSE;
+  
+  Double_t ptg  = pGamma->Pt();
+  Double_t phig = pGamma->Phi();
+  Double_t ptl  = pLeading->Pt();
+  Double_t phil = pLeading->Phi();
+  Double_t etal = pLeading->Eta();
+
+  Float_t ptcut = 0. ;
+  if(fPbPb){
+    if(ptg > fPtJetSelectionCut)  ptcut = 2. ;
+    else                          ptcut = 0.5;
+  }
+
+  TIter next(pl) ; 
+  TParticle * particle = 0 ; 
+  while ( (particle = dynamic_cast<TParticle*>(next())) ) {
+    
+    b0 = kFALSE;
+    b1 = kFALSE;
+
+    //Particles in jet 
+    SetJet(particle, b0, fCone, etal, phil) ;  
+
+    if(b0){
+      new((*jetList)[n0++]) TParticle(*particle) ;
+      particle->Momentum(lv);
+      if(particle->Pt() > ptcut ){
+       jet+=lv;
+       ptjet+=particle->Pt();
+      }
+    }
+
+    //Background around (phi_gamma-pi, eta_leading)
+    SetJet(particle, b1, fCone,etal, phig) ;
+
+    if(b1) { 
+      new((*bkgList)[n1++]) TParticle(*particle) ;
+      particle->Momentum(lv);
+      if(particle->Pt() > ptcut ){
+       bkg+=lv;
+       ptbkg+=particle->Pt();    
+      }  
+    }
+  }
+  
+  ptjet = jet.Pt();
+  ptbkg = bkg.Pt();
+
+  if(ptjet > 0.) {
+
+    AliDebug(2,Form("Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
+    
+    //Fill histograms
+    
+    Double_t ratjet   = ptjet/ptg ;
+    Double_t ratbkg  = ptbkg/ptg ;
+    
+    dynamic_cast<TH2F*>
+      (fOutputContainer->FindObject("JetRatio"+lastname))
+      ->Fill(ptg,ratjet);       
+    dynamic_cast<TH2F*>
+      (fOutputContainer->FindObject("JetPt"+lastname))
+      ->Fill(ptg,ptjet);
+    
+    dynamic_cast<TH2F*>
+      (fOutputContainer->FindObject("BkgRatio"+lastname))
+      ->Fill(ptg,ratbkg);
+    
+    dynamic_cast<TH2F*>
+      (fOutputContainer->FindObject("BkgPt"+lastname))
+      ->Fill(ptg,ptbkg);
+
+
+    //Jet selection
+    Bool_t kSelect = kFALSE;
+    if(fSelect == 0)
+      kSelect = kTRUE; //Accept all jets, no restriction
+    else if(fSelect == 1){
+      //Selection with parametrized cuts
+      if(IsJetSelected(ptg,ptjet))   kSelect = kTRUE;
+    }
+    else if(fSelect == 2){
+      //Simple selection
+      if(!fJetsOnlyInCTS){
+       if((ratjet <  fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
+      }
+      else{
+       if((ratjet <  fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
+      }
+    }
+    else
+      AliError("Jet selection option larger than 2, DONT SELECT JETS");
+    
+    
+    if(kSelect){
+      if (GetPrintInfo())
+       AliInfo(Form("Jet Selected: pt %f ", ptjet)) ;
+      
+      FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
+      FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
+    }
+  } //ptjet > 0
+  
+  jetList ->Delete();
+  bkgList ->Delete();
+  
+}
+
+//____________________________________________________________________________
+void AliAnaGammaJet::Print(const Option_t * opt) const
+{
+
+  //Print some relevant parameters set for the analysis
+  if(! opt)
+    return;
+
+  Info("Print", "%s %s", GetName(), GetTitle() ) ;
+  if(!fJetsOnlyInCTS)
+    printf("EMCAL Acceptance cut  : phi min %f, phi max %f, eta %f \n", fPhiEMCALCut[0],  fPhiEMCALCut[1], fEtaEMCALCut) ;
+  
+  printf("Phi_Leading                                <     %f\n", fPhiMaxCut) ; 
+  printf("Phi_Leading                                >     %f\n", fPhiMinCut) ;
+  printf("pT Leading / pT Gamma             <     %f\n", fRatioMaxCut) ; 
+  printf("pT Leading / pT Gamma             >     %f\n", fRatioMinCut) ;
+  printf("M_pair                                        <     %f\n", fInvMassMaxCut) ; 
+  printf("M_pair                                        >     %f\n", fInvMassMinCut) ; 
+  if(fSelect == 2){
+    printf("pT Jet / pT Gamma                     <    %f\n", fJetRatioMaxCut) ; 
+    printf("pT Jet / pT Gamma                     >    %f\n", fJetRatioMinCut) ;
+    printf("pT Jet (Only CTS)/ pT Gamma   <    %f\n", fJetCTSRatioMaxCut) ; 
+    printf("pT Jet (Only CTS)/ pT Gamma   >    %f\n", fJetCTSRatioMinCut) ;
+  }
+
+} 
+
+//___________________________________________________________________
+void AliAnaGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone, 
+                            Double_t eta, Double_t phi)
+{
+
+  //Check if the particle is inside the cone defined by the leading particle
+  b = kFALSE;
+  
+  if(phi > TMath::TwoPi())
+    phi-=TMath::TwoPi();
+  if(phi < 0.)
+    phi+=TMath::TwoPi();
+  
+  Double_t  rad = 10000 + cone;
+  
+  if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
+    rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
+                     TMath::Power(part->Phi()-phi,2));
+  else{
+    if(part->Phi()-phi > TMath::TwoPi() - cone)
+      rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
+                       TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
+    if(part->Phi()-phi < -(TMath::TwoPi() - cone))
+      rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
+                       TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
+  }
+
+  if(rad < cone )
+    b = kTRUE;
+  
+}
+
+
+void AliAnaGammaJet::Terminate(Option_t *)
+{
+   // The Terminate() function is the last function to be called during
+   // a query. It always runs on the client, it can be used to present
+   // the results graphically or save the results to file.
+    
+
+}
diff --git a/PWG4/AliAnaGammaJet.h b/PWG4/AliAnaGammaJet.h
new file mode 100644 (file)
index 0000000..b6714c2
--- /dev/null
@@ -0,0 +1,233 @@
+#ifndef ALIANAGAMMAJET_H
+#define ALIANAGAMMAJET_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+/* $Id$ */
+
+/* History of cvs commits:
+ *
+ * $Log$
+ *
+ */
+
+//_________________________________________________________________________
+//  Class for the analysis of gamma-jet correlations.     
+//  Basically it seaches for a prompt photon in the Calorimeters acceptance, 
+//  if so we construct a jet around the highest pt particle in the opposite 
+//  side in azimuth. This jet has to fullfill several conditions to be 
+//  accepted. Then the fragmentation function of this jet is constructed 
+//  Class created from old AliPHOSGammaJet
+
+//*-- 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 "AliAnaGammaDirect.h" 
+
+class AliESD ; 
+class AliAnaGammaJet : public AliAnaGammaDirect {
+
+public: 
+
+  AliAnaGammaJet(const char *name) ; // default ctor
+  AliAnaGammaJet(const AliAnaGammaJet & gj) ; // cpy ctor
+  virtual ~AliAnaGammaJet() ; //virtual dtor
+  virtual void Exec(Option_t * opt = "") ;
+  virtual void Init(Option_t * opt = "");
+  virtual void Terminate(Option_t * opt = "");
+
+  Bool_t     AreJetOnlyInCTS() const {return fJetsOnlyInCTS ; } 
+  Double_t GetAngleMaxParam(Int_t i) const {return fAngleMaxParam.At(i) ; }
+  Double_t GetEtaEMCALCut() const {return fEtaEMCALCut;}
+  Double_t GetPhiEMCALCut(Int_t i) const {return fPhiEMCALCut[i];}
+  Double_t GetInvMassMaxCut() const {return fInvMassMaxCut ; }
+  Double_t GetInvMassMinCut() const {return fInvMassMinCut ; }
+  Double_t GetPhiMaxCut() const {return fPhiMaxCut ; }
+  Double_t GetPhiMinCut() const {return fPhiMinCut ; }
+  Double_t GetPtJetSelectionCut() const {return fPtJetSelectionCut ; }
+  Double_t GetJetRatioMaxCut() const {return fJetRatioMaxCut ; }
+  Double_t GetJetRatioMinCut() const {return fJetRatioMinCut ; }
+  Double_t GetRatioMaxCut() const {return fRatioMaxCut ; }
+  Double_t GetRatioMinCut() const {return fRatioMinCut ; }
+  Int_t       GetNCones() const {return fNCone ; }
+  Int_t       GetNPtThres() const {return fNPt ; }
+  Float_t    GetCone() const {return fCone ; }
+  Float_t    GetPtThreshold() const {return fPtThreshold ; }
+  Float_t    GetCones(Int_t i) const {return fCones[i] ; }
+  Float_t    GetPtThreshold(Int_t i) const {return fPtThres[i] ; }
+  TString   GetConeName(Int_t i) const {return fNameCones[i] ; }
+  TString   GetPtThresName(Int_t i) const {return fNamePtThres[i] ; }
+  
+  Bool_t   AreSeveralConeAndPtCuts() const {return fSeveralConeAndPtCuts ; }
+  Bool_t   IsPbPb() const {return fPbPb ; }
+
+
+  void Print(const Option_t * opt)const;
+  
+  void SetAngleMaxParam(Int_t i, Double_t par)
+  {fAngleMaxParam.AddAt(par,i) ; }
+  void SetSeveralConeAndPtCuts(Bool_t several){fSeveralConeAndPtCuts = several ;}
+  void SetEtaEMCALCut(Double_t etacut) {fEtaEMCALCut = etacut ; }
+  void SetPhiEMCALCut(Double_t phi, Int_t i){fPhiEMCALCut[i] = phi; }
+  void SetPbPb(Bool_t opt){fPbPb = opt; }
+  void SetNCones(Int_t n){fNCone = n ; }
+  void SetNPtThresholds(Int_t n){fNPt = n ; }
+  void SetCones(Int_t i, Float_t cone, TString sc)
+    {fCones[i] = cone ; fNameCones[i] = sc; };
+  void SetCone(Float_t cone)
+    {fCone = cone; }
+  void SetPtThreshold(Float_t pt){fPtThreshold = pt; };
+  void SetPtThresholds(Int_t i,Float_t pt, TString spt){fPtThres[i] = pt ; 
+  fNamePtThres[i] = spt; };
+  void SetInvMassCutRange(Double_t invmassmin, Double_t invmassmax)
+    {fInvMassMaxCut =invmassmax;  fInvMassMinCut =invmassmin;} 
+  void SetJetRatioCutRange(Double_t ratiomin, Double_t ratiomax)
+    {fJetRatioMaxCut =ratiomax;  fJetRatioMinCut = ratiomin ; }
+  void SetJetsOnlyInCTS(Bool_t opt){fJetsOnlyInCTS = opt; }
+  void SetJetCTSRatioCutRange(Double_t ratiomin, Double_t ratiomax)
+    {fJetCTSRatioMaxCut =ratiomax;  fJetCTSRatioMinCut = ratiomin ; }
+  void SetPhiCutRange(Double_t phimin, Double_t phimax)
+  {fPhiMaxCut =phimax;  fPhiMinCut =phimin;}
+  void SetPtJetSelectionCut(Double_t cut){fPtJetSelectionCut = cut; }
+  void SetJetSelection(UInt_t select){ fSelect= select ; }
+  void SetRatioCutRange(Double_t ratiomin, Double_t ratiomax)
+  {fRatioMaxCut = ratiomax;  fRatioMinCut = ratiomin;}
+
+
+  private:
+
+  Double_t CalculateJetRatioLimit(const Double_t ptg, const Double_t *param, 
+                                 const Double_t *x);
+  
+  void FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl,TString type, TString lastname);
+
+  Bool_t IsAngleInWindow(const Float_t angle, const Float_t e);
+  Bool_t IsJetSelected(const Double_t ptg, const Double_t ptjet);
+
+  void MakeJet(TClonesArray * particleList, TParticle * pGamma, TParticle* pLeading, TString lastname); 
+
+  void GetLeadingCharge(TClonesArray * pl, TParticle *pGamma, TParticle * pLeading) const ;
+  void GetLeadingPi0   (TClonesArray * pl, TParticle *pGamma, TParticle * pLeading)  ;
+  Bool_t GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe, TParticle *pGamma, TParticle * pLeading) ;
+  void MakeHistos() ;
+  void SetJet(TParticle * part, Bool_t & b, Float_t cone, Double_t eta, 
+             Double_t phi);
+
+ private:
+
+  //  TTree       *fChain ;   //!pointer to the analyzed TTree or TChain
+  //AliESD       *fESD ;     //! Declaration of leave types
+  Bool_t       fSeveralConeAndPtCuts;     //  To play with the jet cone size and pt th.
+  //Bool_t       fPrintInfo ;       //Print most interesting information on screen
+                              // and give interesting information
+
+  Bool_t       fPbPb;          // PbPb event
+  Bool_t       fJetsOnlyInCTS ;    // Jets measured only in TPC+ITS.
+  Double_t   fEtaEMCALCut ;         // Eta EMCAL acceptance
+  Double_t   fPhiEMCALCut[2] ; // Phi cut maximum
+  Double_t   fPhiMaxCut ;      // Phi EMCAL maximum acceptance
+  Double_t   fPhiMinCut ;      // Phi EMCAL minimum acceptance
+  // Double_t   fGammaPtCut ;  // Min pt in Calorimeter
+  Double_t   fInvMassMaxCut ;  // Invariant Mass cut maximum
+  Double_t   fInvMassMinCut ;  // Invariant Masscut minimun
+  Double_t   fRatioMaxCut ;    // Leading particle/gamma Ratio cut maximum
+  Double_t   fRatioMinCut ;    // Leading particle/gamma Ratio cut minimum
+
+  //Jet selection parameters
+  //Fixed cuts (old)
+  Double_t   fJetCTSRatioMaxCut ; // Leading particle/gamma Ratio cut maximum
+  Double_t   fJetCTSRatioMinCut ; // Leading particle/gamma Ratio cut minimum
+  Double_t   fJetRatioMaxCut ; // Jet/gamma Ratio cut maximum
+  Double_t   fJetRatioMinCut ; // Jet/gamma Ratio cut minimum
+
+  //Cuts depending on jet pt
+  Double_t fJetE1[2];    //Rec. jet energy parameters
+  Double_t fJetE2[2];    //Rec. jet energy parameters
+  Double_t fJetSigma1[2];//Rec. sigma of jet energy  parameters
+  Double_t fJetSigma2[2];//Rec. sigma of jet energy  parameters
+  Double_t fBkgMean[6];  //Background mean energy 
+  Double_t fBkgRMS[6];   //Background RMS
+  Double_t fJetXMin1[6]; //X Factor to set jet min limit for pp
+  Double_t fJetXMin2[6]; //X Factor to set jet min limit for PbPb
+  Double_t fJetXMax1[6]; //X Factor to set jet max limit for pp
+  Double_t fJetXMax2[6]; //X Factor to set jet max limit for PbPb
+
+  Int_t         fNCone ;            // Number of jet cones sizes
+  Int_t         fNPt   ;            // Number of jet particle pT threshold
+  Double_t   fCone  ;            // Jet cone sizes under study (!fSeveralConeAndPtCuts)
+  Double_t   fCones[10];         // Jet cone sizes under study (fSeveralConeAndPtCuts)
+  TString     fNameCones[10];     // String name of cone to append to histos
+  Double_t   fPtThreshold;       // Jet pT threshold under study(!fSeveralConeAndPtCuts)
+  Double_t   fPtThres[10];       // Jet pT threshold under study(fSeveralConeAndPtCuts)
+  Double_t   fPtJetSelectionCut; // Jet pt to change to low pt jets analysis
+  TObjArray  *fOutputContainer ; //! output data container
+  TString     fNamePtThres[10];   // String name of pt th to append to histos
+  TArrayD    fAngleMaxParam ; //Max opening angle selection parameters
+  UInt_t       fSelect  ;   //kTRUE: Selects all jets, no limits.
+  //TString     fCalorimeter ; //PHOS or EMCAL detects Gamma
+  //Bool_t      fEMCALPID ; //Fill EMCAL particle lists with particles with corresponding pid
+  //Bool_t      fPHOSPID;   //Fill PHOS particle lists with particles with corresponding pid
+
+  //Histograms
+
+  TH2F * fhChargeRatio  ; 
+  TH2F * fhPi0Ratio   ; 
+  TH2F * fhDeltaPhiCharge  ;  
+  TH2F * fhDeltaPhiPi0   ; 
+  TH2F * fhDeltaEtaCharge  ; 
+  TH2F * fhDeltaEtaPi0  ; 
+  TH2F * fhAnglePair  ; 
+  TH2F * fhAnglePairAccepted  ; 
+  TH2F * fhAnglePairNoCut  ; 
+  TH2F * fhAnglePairLeadingCut  ; 
+  TH2F * fhAnglePairAngleCut   ; 
+  TH2F * fhAnglePairAllCut   ; 
+  TH2F * fhAnglePairLeading  ; 
+  TH2F * fhInvMassPairNoCut    ; 
+  TH2F * fhInvMassPairLeadingCut  ; 
+  TH2F * fhInvMassPairAngleCut  ; 
+  TH2F * fhInvMassPairAllCut   ; 
+  TH2F * fhInvMassPairLeading  ; 
+  TH1F * fhNBkg   ; 
+  TH2F * fhNLeading  ; 
+  TH1F * fhNJet  ; 
+  TH2F * fhJetRatio  ; 
+  TH2F * fhJetPt   ; 
+  TH2F * fhBkgRatio   ; 
+  TH2F * fhBkgPt  ; 
+  TH2F * fhJetFragment  ; 
+  TH2F * fhBkgFragment  ; 
+  TH2F * fhJetPtDist  ; 
+  TH2F * fhBkgPtDist  ; 
+
+  TH2F * fhJetRatios[5][5];
+  TH2F * fhJetPts[5][5];
+  TH2F * fhBkgRatios[5][5];  
+  TH2F * fhBkgPts[5][5];
+  
+  TH2F * fhNLeadings[5][5];
+  TH1F * fhNJets[5][5];
+  TH1F * fhNBkgs[5][5];
+  
+  TH2F * fhJetFragments[5][5];
+  TH2F * fhBkgFragments[5][5];
+  TH2F * fhJetPtDists[5][5];
+  TH2F * fhBkgPtDists[5][5];
+  
+  ClassDef(AliAnaGammaJet,0)
+} ;
+
+#endif //ALIANAGAMMAJET_H
+
+
+
diff --git a/PWG4/GammaLinkDef.h b/PWG4/GammaLinkDef.h
new file mode 100644 (file)
index 0000000..70fa9ca
--- /dev/null
@@ -0,0 +1,11 @@
+#ifdef __CINT__
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliAnaGammaDirect+;
+#pragma link C++ class AliAnaGammaHadron+;
+#pragma link C++ class AliAnaGammaJet+;
+
+#endif
diff --git a/PWG4/Makefile b/PWG4/Makefile
new file mode 100644 (file)
index 0000000..d896156
--- /dev/null
@@ -0,0 +1,89 @@
+
+include $(ROOTSYS)/test/Makefile.arch
+
+default-target: libGamma.so
+
+ALICEINC      = -I.
+
+### define include dir for local case and par case
+ifneq ($(ANALYSIS_NEW_INCLUDE),)
+  ALICEINC += -I../$(ESD_INCLUDE) -I../$(ANALYSIS_NEW_INCLUDE)
+else
+  ifneq ($(ALICE_ROOT),)
+    ALICEINC += -I$(ALICE_ROOT)/include
+  endif
+endif
+
+# for building of Gamma.par
+ifneq ($(GAMMA_INCLUDE),)
+  ALICEINC += -I../$(GAMMA_INCLUDE)
+endif
+
+CXXFLAGS     += $(ALICEINC) -g
+
+PACKAGE = Gamma
+include lib$(PACKAGE).pkg
+
+DHDR_Gamma := $(DHDR)
+HDRS_Gamma := $(HDRS)
+SRCS_Gamma := $(SRCS) G__$(PACKAGE).cxx
+OBJS_Gamma := $(SRCS_Gamma:.cxx=.o)
+
+PARFILE       = $(PACKAGE).par
+
+
+lib$(PACKAGE).so: $(OBJS_Gamma)
+       @echo "Linking" $@ ...
+       @/bin/rm -f $@
+ifeq ($(PLATFORM),macosx)
+       @$(LD) -bundle -undefined $(UNDEFOPT) $(LDFLAGS) $^ -o $@
+else
+       @$(LD) $(SOFLAGS) $(LDFLAGS) $^ -o $@
+endif
+       @chmod a+x $@
+       @echo "done"
+
+%.o:    %.cxx %.h
+       $(CXX) $(CXXFLAGS) -c $< -o $@
+
+clean:
+       @rm -f $(OBJS_Gamma) *.so G__$(PACKAGE).* $(PARFILE)
+
+G__$(PACKAGE).cxx G__$(PACKAGE).h: $(HDRS) $(DHDR)
+       @echo "Generating dictionary ..."
+       rootcint -f $@ -c $(ALICEINC) $^
+
+### CREATE PAR FILE
+
+$(PARFILE): $(patsubst %,$(PACKAGE)/%,$(filter-out G__%, $(HDRS_Gamma) $(SRCS_Gamma) $(DHDR_Gamma) Makefile Makefile.arch lib$(PACKAGE).pkg PROOF-INF))
+       @echo "Creating archive" $@ ...
+       @tar cfzh $@ $(PACKAGE)
+       @rm -rf $(PACKAGE)
+       @echo "done"
+
+$(PACKAGE)/Makefile: Makefile #.$(PACKAGE)
+       @echo Copying $< to $@ with transformations
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @sed 's/include \$$(ROOTSYS)\/test\/Makefile.arch/include Makefile.arch/' < $^ > $@
+
+$(PACKAGE)/Makefile.arch: $(ROOTSYS)/test/Makefile.arch
+       @echo Copying $< to $@
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @cp -a $^ $@
+
+$(PACKAGE)/PROOF-INF: PROOF-INF.$(PACKAGE)
+       @echo Copying $< to $@
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @cp -a -r $^ $@
+
+$(PACKAGE)/%: %
+       @echo Copying $< to $@
+       @[ -d $(dir $@) ] || mkdir -p $(dir $@)
+       @cp -a $< $@
+
+test-%.par: %.par
+       @echo "INFO: The file $< is now tested, in case of an error check in par-tmp."
+       @mkdir -p par-tmp
+       @cd par-tmp; tar xfz ../$<;     cd $(subst .par,,$<); PROOF-INF/BUILD.sh
+       @rm -rf par-tmp
+       @echo "INFO: Testing succeeded (already cleaned up)"
diff --git a/PWG4/libGamma.pkg b/PWG4/libGamma.pkg
new file mode 100644 (file)
index 0000000..95caba9
--- /dev/null
@@ -0,0 +1,9 @@
+SRCS = AliAnaGammaDirect.cxx AliAnaGammaHadron.cxx AliAnaGammaJet.cxx
+
+HDRS:= $(SRCS:.cxx=.h) 
+
+DHDR:= GammaLinkDef.h
+
+EXPORT:=$(SRCS:.cxx=.h)
+
+