]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/PartCorrDep/AliAnaParticleJetLeadingConeCorrelation.cxx
change summary plots name accoring to convention proposed by Anton (Markus F)
[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaParticleJetLeadingConeCorrelation.cxx
index 6ebf5e429e1197fc275effb4e407dd2fba7133fd..903a399f7922c7a494daaf998aa0d27df0770d29 100755 (executable)
@@ -17,7 +17,7 @@
 //_________________________________________________________________________
 // Class that contains the algorithm for the reconstruction of jet, cone around leading particle
 // The seed is a backward particle (direct photon)
-// 1)Take the trigger particle stored in AliAODPWG4ParticleCorrelation,
+// 1) Take the trigger particle stored in AliAODPWG4ParticleCorrelation,
 // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
 // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
 //
 
 //---- Analysis system ----
 #include "AliAODTrack.h"
-#include "AliAODCaloCluster.h"
+#include "AliVCluster.h"
 #include "AliCaloTrackReader.h"
 #include "AliNeutralMesonSelection.h"
 #include "AliAnaParticleJetLeadingConeCorrelation.h"  
 #include "AliCaloPID.h"
 #include "AliAODPWG4ParticleCorrelation.h"
-#include "AliFidutialCut.h"
+#include "AliFiducialCut.h"
 #include "AliAODTrack.h"
-#include "AliAODCaloCluster.h"
 
 ClassImp(AliAnaParticleJetLeadingConeCorrelation)
 
@@ -70,8 +69,8 @@ ClassImp(AliAnaParticleJetLeadingConeCorrelation)
     fhNeutralLeadingPt(0),fhNeutralLeadingPhi(0),fhNeutralLeadingEta(0),
     fhNeutralLeadingDeltaPt(0),fhNeutralLeadingDeltaPhi(0),fhNeutralLeadingDeltaEta(0),
     fhNeutralLeadingRatioPt(0),fhChargedLeadingXi(0), fhNeutralLeadingXi(0),
-       fhChargedLeadingDeltaPhiRatioPt30(0), fhNeutralLeadingDeltaPhiRatioPt30(0),
-       fhChargedLeadingDeltaPhiRatioPt50(0), fhNeutralLeadingDeltaPhiRatioPt50(0),
+    fhChargedLeadingDeltaPhiRatioPt30(0), fhNeutralLeadingDeltaPhiRatioPt30(0),
+    fhChargedLeadingDeltaPhiRatioPt50(0), fhNeutralLeadingDeltaPhiRatioPt50(0),
     //Jet
     fhJetPt(0),fhJetRatioPt(0),fhJetDeltaPhi(0), fhJetDeltaEta(0),
     fhJetLeadingRatioPt(0),fhJetLeadingDeltaPhi(0),fhJetLeadingDeltaEta(0),
@@ -141,7 +140,7 @@ ClassImp(AliAnaParticleJetLeadingConeCorrelation)
   InitParameters();
 
 }
-
+/*
 //____________________________________________________________________________
 AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation(const AliAnaParticleJetLeadingConeCorrelation & jetlc) :   
   AliAnaPartCorrBaseClass(jetlc), fJetsOnlyInCTS(jetlc.fJetsOnlyInCTS), fPbPb(jetlc.fPbPb), 
@@ -350,25 +349,11 @@ AliAnaParticleJetLeadingConeCorrelation & AliAnaParticleJetLeadingConeCorrelatio
   return *this;
 
 }
-
+*/
 //____________________________________________________________________________
 AliAnaParticleJetLeadingConeCorrelation::~AliAnaParticleJetLeadingConeCorrelation() 
 {
    // Remove all pointers except analysis output pointers.
-  delete [] fJetE1;  
-  delete [] fJetE2;    
-  delete [] fJetSigma1;
-  delete [] fJetSigma2;
-  delete [] fBkgMean; 
-  delete [] fBkgRMS;  
-  delete [] fJetXMin1;
-  delete [] fJetXMin2;
-  delete [] fJetXMax1;
-  delete [] fJetXMax2; 
-  delete [] fJetCones;         
-  delete [] fJetNameCones;   
-  delete [] fJetPtThres;       
-  delete [] fJetNamePtThres;  
 }
 
 //____________________________________________________________________________
@@ -403,32 +388,41 @@ void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCo
   Double_t etaJet = jet.Eta();
   Double_t etaLead = leading.Eta();
   
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"Pt"+lastname))->
-    Fill(ptTrig,ptJet);
+  TH2F *h1 = 0x0;
+  h1 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h1)h1->Fill(ptTrig,ptJet);
   
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"RatioPt"+lastname))->
-    Fill(ptTrig,ptJet/ptTrig);
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"LeadingRatioPt"+lastname))->
-    Fill(ptTrig,ptLead/ptJet);
-  //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"Phi"+lastname))->
-  //     Fill(ptTrig,phiJet);
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"DeltaPhi"+lastname))->
-    Fill(ptTrig,phiJet-phiTrig);
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"LeadingDeltaPhi"+lastname))->
-    Fill(ptTrig,phiJet-phiLead);
+  TH2F *h2 = 0x0;
+  h2 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h2) h2->Fill(ptTrig,ptJet/ptTrig);
+  
+  TH2F *h3 = 0x0;
+  h3 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h3)h3->Fill(ptTrig,ptLead/ptJet);
   
-  //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"Eta"+lastname))->
+  //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
+  //     Fill(ptTrig,phiJet);
+  TH2F *h4 = 0x0;
+  h4 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h4) h4->Fill(ptTrig,phiJet-phiTrig);
+  TH2F *h5 = 0x0;
+  h5 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h5) h5->Fill(ptTrig,phiJet-phiLead);
+  
+  //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
   //     Fill(ptTrig,etaJet);
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"DeltaEta"+lastname))->
-    Fill(ptTrig,etaJet-etaTrig);
-  dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"LeadingDeltaEta"+lastname))->
-    Fill(ptTrig,etaJet-etaLead);
+  TH2F *h6 = 0x0;
+  h6 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h6) h6->Fill(ptTrig,etaJet-etaTrig);
+  TH2F *h7 = 0x0;
+  h7 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+  if(h7) h7->Fill(ptTrig,etaJet-etaLead);
   
   //Construct fragmentation function
-  TRefArray * pl = new TRefArray;
+  TObjArray * pl = new TObjArray;
   
-  if(type == "Jet") pl = particle->GetRefArray(GetAODRefArrayName()+"Tracks");
-  else if(type == "Bkg") particle->GetRefArray(GetAODRefArrayName()+"TracksBkg");
+  if(type == "Jet") pl = particle->GetObjArray(Form("%sTracks",GetAODObjArrayName().Data()));
+  else if(type == "Bkg") particle->GetObjArray(Form("%sTracksBkg",GetAODObjArrayName().Data()));
   
   if(!pl) return ;
   
@@ -442,7 +436,8 @@ void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCo
   
   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
     AliAODTrack* track = dynamic_cast<AliAODTrack *>(pl->At(ipr)) ;
-    p3.SetXYZ(track->Px(),track->Py(),track->Pz());
+    if(track)p3.SetXYZ(track->Px(),track->Py(),track->Pz());
+    else printf("AliAnaParticleJetLeadingConeCorrelation::FillJetHistos() - Track not available\n");
     
     //Recheck if particle is in jet cone
     if(fReMakeJet || fSeveralConeAndPtCuts)
@@ -450,17 +445,19 @@ void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCo
     
     nTracksInCone++; 
     
-    dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"FFz"+lastname))
-      ->Fill(ptTrig,p3.Pt()/ptTrig);
-    dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"FFxi"+lastname))
-      ->Fill(ptTrig,TMath::Log(ptTrig/p3.Pt()));
-    dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"FFpt"+lastname))
-      ->Fill(ptTrig,p3.Pt());
+    TH2F *ha =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFz%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+    if(ha) ha->Fill(ptTrig,p3.Pt()/ptTrig);
+    TH2F *hb  =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFxi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+    if(hb) hb->Fill(ptTrig,TMath::Log(ptTrig/p3.Pt()));
+    TH2F *hc =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFpt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+    if(hc) hc->Fill(ptTrig,p3.Pt());
     
   }//track loop
   
-  if(nTracksInCone > 0) dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(GetAddedHistogramsStringToName()+type+"NTracksInCone"+lastname))
-    ->Fill(ptTrig, nTracksInCone);
+  if(nTracksInCone > 0) {
+    TH2F *hd = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sNTracksInCone%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
+    hd->Fill(ptTrig, nTracksInCone);
+  }
   
 }
 
@@ -475,9 +472,9 @@ TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
   fOutCont = new TList() ; 
   fOutCont->SetName("ParticleJetLeadingInConeCorrelationHistograms") ; 
   
-  Int_t nptbins  = GetHistoNPtBins();
-  Int_t nphibins = GetHistoNPhiBins();
-  Int_t netabins = GetHistoNEtaBins();
+  Int_t nptbins  = GetHistoPtBins();
+  Int_t nphibins = GetHistoPhiBins();
+  Int_t netabins = GetHistoEtaBins();
   Float_t ptmax  = GetHistoPtMax();
   Float_t phimax = GetHistoPhiMax();
   Float_t etamax = GetHistoEtaMax();
@@ -717,47 +714,47 @@ TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
        TString lastnametitle =", cone ="+fJetNameCones[icone]+", pt > " +fJetNamePtThres[ipt]+" GeV/c";
        
        //Jet Distributions
-       fhJetPts[icone][ipt] = new TH2F("JetPt"+lastnamehist,"p_{T  jet} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
+       fhJetPts[icone][ipt] = new TH2F(Form("JetPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
        fhJetPts[icone][ipt]->SetYTitle("p_{T  jet}");
        fhJetPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetRatioPts[icone][ipt] = new TH2F("JetRatioPt"+lastnamehist,"p_{T  jet}/p_{T trigger} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,2); 
+       fhJetRatioPts[icone][ipt] = new TH2F(Form("JetRatioPt%s",lastnamehist.Data()),Form("p_{T  jet}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
        fhJetRatioPts[icone][ipt]->SetYTitle("p_{T  jet}/p_{T trigger}");
        fhJetRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetDeltaPhis[icone][ipt] = new TH2F("JetDeltaPhi"+lastnamehist,"#phi_{jet} - #phi_{trigger} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
+       fhJetDeltaPhis[icone][ipt] = new TH2F(Form("JetDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
        fhJetDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
        fhJetDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetDeltaEtas[icone][ipt] = new TH2F("JetDeltaEta"+lastnamehist,"#eta_{jet} - #eta_{trigger} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,-2,2); 
+       fhJetDeltaEtas[icone][ipt] = new TH2F(Form("JetDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
        fhJetDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
        fhJetDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetLeadingRatioPts[icone][ipt] = new TH2F("JetLeadingRatioPt"+lastnamehist,"p_{T  jet} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,2); 
+       fhJetLeadingRatioPts[icone][ipt] = new TH2F(Form("JetLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
        fhJetLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T jet}");
        fhJetLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetLeadingDeltaPhis[icone][ipt] = new TH2F("JetLeadingDeltaPhi"+lastnamehist,"#phi_{jet} - #phi_{leading} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
+       fhJetLeadingDeltaPhis[icone][ipt] = new TH2F(Form("JetLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
        fhJetLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
        fhJetLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetLeadingDeltaEtas[icone][ipt] = new TH2F("JetLeadingDeltaEta"+lastnamehist,"#eta_{jet} - #eta_{leading} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,-2,2); 
+       fhJetLeadingDeltaEtas[icone][ipt] = new TH2F(Form("JetLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
        fhJetLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
        fhJetLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhJetFFzs[icone][ipt] = new TH2F("JetFFz"+lastnamehist,"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2); 
+       fhJetFFzs[icone][ipt] = new TH2F(Form("JetFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2); 
        fhJetFFzs[icone][ipt]->SetYTitle("z");
        fhJetFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
        
-       fhJetFFxis[icone][ipt] = new TH2F("JetFFxi"+lastnamehist,"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
+       fhJetFFxis[icone][ipt] = new TH2F(Form("JetFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
        fhJetFFxis[icone][ipt]->SetYTitle("#xi");
        fhJetFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
        
-       fhJetFFpts[icone][ipt] = new TH2F("JetFFpt"+lastnamehist,"p_{T charged hadron } in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
+       fhJetFFpts[icone][ipt] = new TH2F(Form("JetFFpt%s",lastnamehist.Data()),"p_{T charged hadron } in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
        fhJetFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
        fhJetFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
        
-       fhJetNTracksInCones[icone][ipt] = new TH2F("JetNTracksInCone"+lastnamehist,"N particles in cone vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,5000,0, 5000); 
+       fhJetNTracksInCones[icone][ipt] = new TH2F(Form("JetNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000); 
        fhJetNTracksInCones[icone][ipt]->SetYTitle("N tracks in jet cone");
        fhJetNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
@@ -774,47 +771,47 @@ TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
        fOutCont->Add(fhJetNTracksInCones[icone][ipt]) ;
        
        //Bkg Distributions
-       fhBkgPts[icone][ipt] = new TH2F("BkgPt"+lastnamehist,"p_{T  bkg} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
+       fhBkgPts[icone][ipt] = new TH2F(Form("BkgPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
        fhBkgPts[icone][ipt]->SetYTitle("p_{T  bkg}");
        fhBkgPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgRatioPts[icone][ipt] = new TH2F("BkgRatioPt"+lastnamehist,"p_{T  bkg}/p_{T trigger} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,2); 
+       fhBkgRatioPts[icone][ipt] = new TH2F(Form("BkgRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
        fhBkgRatioPts[icone][ipt]->SetYTitle("p_{T  bkg}/p_{T trigger}");
        fhBkgRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgDeltaPhis[icone][ipt] = new TH2F("BkgDeltaPhi"+lastnamehist,"#phi_{bkg} - #phi_{trigger} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
+       fhBkgDeltaPhis[icone][ipt] = new TH2F(Form("BkgDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
        fhBkgDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
        fhBkgDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgDeltaEtas[icone][ipt] = new TH2F("BkgDeltaEta"+lastnamehist,"#eta_{bkg} - #eta_{trigger} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,-2,2); 
+       fhBkgDeltaEtas[icone][ipt] = new TH2F(Form("BkgDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
        fhBkgDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
        fhBkgDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgLeadingRatioPts[icone][ipt] = new TH2F("BkgLeadingRatioPt"+lastnamehist,"p_{T  bkg} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,2); 
+       fhBkgLeadingRatioPts[icone][ipt] = new TH2F(Form("BkgLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2); 
        fhBkgLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T bkg}");
        fhBkgLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgLeadingDeltaPhis[icone][ipt] = new TH2F("BkgLeadingDeltaPhi"+lastnamehist,"#phi_{bkg} - #phi_{leading} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
+       fhBkgLeadingDeltaPhis[icone][ipt] = new TH2F(Form("BkgLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi()); 
        fhBkgLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
        fhBkgLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgLeadingDeltaEtas[icone][ipt] = new TH2F("BkgLeadingDeltaEta"+lastnamehist,"#eta_{bkg} - #eta_{leading} vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,120,-2,2); 
+       fhBkgLeadingDeltaEtas[icone][ipt] = new TH2F(Form("BkgLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2); 
        fhBkgLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
        fhBkgLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
-       fhBkgFFzs[icone][ipt] = new TH2F("BkgFFz"+lastnamehist,"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2); 
+       fhBkgFFzs[icone][ipt] = new TH2F(Form("BkgFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2); 
        fhBkgFFzs[icone][ipt]->SetYTitle("z");
        fhBkgFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
        
-       fhBkgFFxis[icone][ipt] = new TH2F("BkgFFxi"+lastnamehist,"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
+       fhBkgFFxis[icone][ipt] = new TH2F(Form("BkgFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.); 
        fhBkgFFxis[icone][ipt]->SetYTitle("#xi");
        fhBkgFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
        
-       fhBkgFFpts[icone][ipt] = new TH2F("BkgFFpt"+lastnamehist,"p_{T charged hadron} in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
+       fhBkgFFpts[icone][ipt] = new TH2F(Form("BkgFFpt%s",lastnamehist.Data()),"p_{T charged hadron} in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.); 
        fhBkgFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
        fhBkgFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
        
-       fhBkgNTracksInCones[icone][ipt] = new TH2F("BkgNTracksInCone"+lastnamehist,"N particles in cone vs p_{T trigger}"+lastnametitle,nptbins,ptmin,ptmax,5000,0, 5000); 
+       fhBkgNTracksInCones[icone][ipt] = new TH2F(Form("BkgNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000); 
        fhBkgNTracksInCones[icone][ipt]->SetYTitle("N tracks in bkg cone");
        fhBkgNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
        
@@ -840,6 +837,7 @@ TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) fOutCont->Add(nmsHistos->At(i)) ;
+       delete nmsHistos;
   }
   
   
@@ -856,7 +854,7 @@ TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
 
 //____________________________________________________________________________
 Bool_t  AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle(AliAODPWG4ParticleCorrelation *particle, TLorentzVector & pLeading) 
-  const {
+{
   //Search Charged or Neutral leading particle, select the highest one and fill AOD
   
   TLorentzVector pLeadingCh(0,0,0,0) ;
@@ -898,7 +896,7 @@ Bool_t  AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle(AliAODPWG4Pa
 }
 
 //____________________________________________________________________________
-void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4ParticleCorrelation * particle, TLorentzVector & pLeading) const
+void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4ParticleCorrelation* const particle, TLorentzVector & pLeading) const
 {  
   //Search for the charged particle with highest pt and with 
   //Phi=Phi_trigger-Pi and pT=0.1E_gamma 
@@ -938,7 +936,7 @@ void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4Partic
 }
 
 //____________________________________________________________________________
-void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleCorrelation * particle, TLorentzVector & pLeading) const
+void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleCorrelation* const particle, TLorentzVector & pLeading)
 {  
   //Search for the neutral pion with highest pt and with 
   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
@@ -954,85 +952,105 @@ void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleC
     
     TLorentzVector gammai;
     TLorentzVector gammaj;
-    
-    Double_t vertex[] = {0,0,0};
-    if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
-    
+         
+    //Get vertex for photon momentum calculation
+    Double_t vertex [] = {0,0,0} ; //vertex 
+    //Double_t vertex2[] = {0,0,0} ; //vertex of second input AOD 
+    if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) 
+    {
+      GetVertex(vertex);
+      //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
+    }
+         
     //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
     for(Int_t iclus = 0;iclus < GetAODEMCAL()->GetEntriesFast() ; iclus ++ ){
-      AliAODCaloCluster * calo = (AliAODCaloCluster *)(GetAODEMCAL()->At(iclus)) ;
+      AliVCluster * calo = (AliVCluster *)(GetAODEMCAL()->At(iclus)) ;
       
-      //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
+      //Input from second AOD?
+      //Int_t inputi = 0;
+      //         if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) inputi = 1 ;
+      //         else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) inputi = 1;
+      
+      //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
       Int_t pdgi=0;
-      if(!SelectCluster(calo,vertex, gammai, pdgi)) continue ;
+      //if     (inputi == 0 && !SelectCluster(calo, vertex,  gammai, pdgi))  continue ;
+      //else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdgi))  continue ;        
+      if(!SelectCluster(calo, vertex,  gammai, pdgi))  continue ;
       
       if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster: pt %2.3f, phi %2.3f \n", 
-                               gammai.Pt(),gammai.Phi());
+                                gammai.Pt(),gammai.Phi());
       
       //2 gamma overlapped, found with PID
       if(pdgi == AliCaloPID::kPi0){ 
-       
-       if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster ID as Pi0 \n");
-       
-       pt  = gammai.Pt();
-       rat = pt/ptTrig;
-       phi = gammai.Phi();
-       if(phi < 0) phi+=TMath::TwoPi();
-       
-       //Selection within angular and energy limits
-       Float_t deltaphi = TMath::Abs(phiTrig-phi);
-       if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  && 
-          deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut )
-         {
-           phil = phi ;
-           ptl  = pt ;
-           pLeading.SetPxPyPzE(gammai.Px(),gammai.Py(),gammai.Pz(),gammai.E());
-         }// cuts
+        
+        if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster ID as Pi0 \n");
+        
+        pt  = gammai.Pt();
+        rat = pt/ptTrig;
+        phi = gammai.Phi();
+        if(phi < 0) phi+=TMath::TwoPi();
+        
+        //Selection within angular and energy limits
+        Float_t deltaphi = TMath::Abs(phiTrig-phi);
+        if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  && 
+           deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut )
+        {
+          phil = phi ;
+          ptl  = pt ;
+          pLeading.SetPxPyPzE(gammai.Px(),gammai.Py(),gammai.Pz(),gammai.E());
+        }// cuts
       }// pdg = AliCaloPID::kPi0
       //Make invariant mass analysis
       else if(pdgi == AliCaloPID::kPhoton){    
-       //Search the photon companion in case it comes from  a Pi0 decay
-       //Apply several cuts to select the good pair
-       for(Int_t jclus = iclus+1; jclus < GetAODEMCAL()->GetEntriesFast() ; jclus ++ ){
-         AliAODCaloCluster * calo2 = (AliAODCaloCluster *) (GetAODEMCAL()->At(jclus)) ;
-         
-         //Cluster selection, not charged with photon or pi0 id and in fidutial cut
-         Int_t pdgj=0;
-         if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
-         
-         if(pdgj == AliCaloPID::kPhoton ){
-           
-           pt  = (gammai+gammaj).Pt();
-           phi = (gammai+gammaj).Phi();
-               if(phi < 0) phi+=TMath::TwoPi();  
-               rat = pt/ptTrig;
-               
-               //Selection within angular and energy limits
-               Float_t deltaphi = TMath::Abs(phiTrig-phi);  
-               if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, |phiTrig-phi| %2.3f, pt/ptTrig %2.3f, M %2.3f\n",
-                                          pt,phi,(gammai+gammaj).Eta(), deltaphi, rat, (gammai+gammaj).M());
-               
-               if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  && 
-                  deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut ){
-                 //Select good pair (aperture and invariant mass)
-                 if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
-                   phil = phi ;
-                   ptl  = pt ;
-                   pLeading=(gammai+gammaj);   
-                   
-                   if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
-                                              ptl,phil,(gammai+gammaj).Eta(), (gammai+gammaj).M());
-                 }//pi0 selection
-                 
-                 
-               }//Pair selected as leading
-         }//if pair of gammas
-       }//2nd loop
+        //Search the photon companion in case it comes from  a Pi0 decay
+        //Apply several cuts to select the good pair
+        for(Int_t jclus = iclus+1; jclus < GetAODEMCAL()->GetEntriesFast() ; jclus ++ ){
+          AliVCluster * calo2 = (AliVCluster *) (GetAODEMCAL()->At(jclus)) ;
+          
+          //Input from second AOD?
+          //Int_t inputj = 0;
+          //     if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= jclus) inputj = 1;
+          //     else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= jclus) inputj = 1;
+          
+          //Cluster selection, not charged with photon or pi0 id and in fiducial cut
+          Int_t pdgj=0;
+          //if     (inputj == 0 && !SelectCluster(calo2, vertex,  gammaj, pdgj))  continue ;
+          //else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))  continue ;
+          if     (!SelectCluster(calo2, vertex,  gammaj, pdgj))  continue ;
+
+          if(pdgj == AliCaloPID::kPhoton ){
+            
+            pt  = (gammai+gammaj).Pt();
+            phi = (gammai+gammaj).Phi();
+            if(phi < 0) phi+=TMath::TwoPi();  
+            rat = pt/ptTrig;
+            
+            //Selection within angular and energy limits
+            Float_t deltaphi = TMath::Abs(phiTrig-phi);  
+            if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, |phiTrig-phi| %2.3f, pt/ptTrig %2.3f, M %2.3f\n",
+                                       pt,phi,(gammai+gammaj).Eta(), deltaphi, rat, (gammai+gammaj).M());
+            
+            if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  && 
+               deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut ){
+              //Select good pair (aperture and invariant mass)
+              if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
+                phil = phi ;
+                ptl  = pt ;
+                pLeading=(gammai+gammaj);      
+                
+                if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
+                                           ptl,phil,(gammai+gammaj).Eta(), (gammai+gammaj).M());
+              }//pi0 selection
+              
+              
+            }//Pair selected as leading
+          }//if pair of gammas
+        }//2nd loop
       }// if pdg = 22
     }// 1st Loop
     
     if(GetDebug() > 2 && pLeading.Pt() > 0 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Leading EMCAL: pt %2.3f eta %2.3f phi %2.3f pt/Eg %2.3f \n",  
-                                                   pLeading.Pt(), pLeading.Eta(), phil,  pLeading.Pt()/ptTrig) ;
+                                                    pLeading.Pt(), pLeading.Eta(), phil,  pLeading.Pt()/ptTrig) ;
     
   }//EMCAL list exists
 }
@@ -1042,7 +1060,7 @@ void AliAnaParticleJetLeadingConeCorrelation::InitParameters()
 {
   //Initialize the parameters of the analysis.
   SetInputAODName("PWG4Particle");
-  SetAODRefArrayName("JetLeadingCone");    
+  SetAODObjArrayName("JetLeadingCone");    
   AddToHistogramsName("AnaJetLCCorr_");
   
   fJetsOnlyInCTS      = kFALSE ;
@@ -1425,8 +1443,6 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
   Double_t phil     = pLeading.Phi();
   if(phil<0) phil+=TMath::TwoPi();
   Double_t etal     = pLeading.Eta();
-  Bool_t   first    = kTRUE;
-  Bool_t   firstbkg = kTRUE;
 
   //Different pt cut for jet particles in different collisions systems
   Float_t ptcut = fJetPtThreshold;
@@ -1441,8 +1457,8 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
   //Fill jet with tracks
   TVector3 p3;
   //Initialize reference arrays that will contain jet and background tracks
-  TRefArray * reftracks  = new TRefArray;
-  TRefArray * reftracksbkg  = new TRefArray;
+  TObjArray * reftracks  = new TObjArray;
+  TObjArray * reftracksbkg  = new TObjArray;
   
   for(Int_t ipr = 0;ipr < (GetAODCTS())->GetEntriesFast() ; ipr ++ ){
     AliAODTrack* track = (AliAODTrack *)((GetAODCTS())->At(ipr)) ;
@@ -1450,11 +1466,6 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
     
     //Particles in jet 
     if(IsParticleInJetCone(p3.Eta(), p3.Phi(), etal, phil)){
-
-      if(first) {
-       new (reftracks) TRefArray(TProcessID::GetProcessWithUID(track)); 
-       first = kFALSE;
-      }
       
       reftracks->Add(track); 
       
@@ -1467,11 +1478,6 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
     //Background around (phi_gamma-pi, eta_leading)
     else if(IsParticleInJetCone(p3.Eta(),p3.Phi(),etal, phiTrig)) { 
       
-      if(firstbkg) {
-       new (reftracksbkg) TRefArray(TProcessID::GetProcessWithUID(track)); 
-       firstbkg = kFALSE;
-      }
-      
       reftracksbkg->Add(track); 
 
       if(p3.Pt() > ptcut ){
@@ -1483,43 +1489,49 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
   
   //Add referenced tracks to AOD
   if(reftracks->GetEntriesFast() > 0 ){
-    reftracks->SetName(GetAODRefArrayName()+"Tracks");
-    particle->AddRefArray(reftracks);
+    reftracks->SetName(Form("%sTracks",GetAODObjArrayName().Data()));
+    particle->AddObjArray(reftracks);
   }
   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No tracks in jet cone\n");
   if(reftracksbkg->GetEntriesFast() > 0 ){
-    reftracksbkg->SetName(GetAODRefArrayName()+"TracksBkg");
-    particle->AddRefArray(reftracksbkg);
+    reftracksbkg->SetName(Form("%sTracksBkg",GetAODObjArrayName().Data()));
+    particle->AddObjArray(reftracksbkg);
   }
   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background tracks in jet cone\n");
   
   //Add neutral particles to jet
   //Initialize reference arrays that will contain jet and background tracks
-  TRefArray * refclusters  = new TRefArray;
-  TRefArray * refclustersbkg  = new TRefArray;
+  TObjArray * refclusters  = new TObjArray;
+  TObjArray * refclustersbkg  = new TObjArray;
   if(!fJetsOnlyInCTS && GetAODEMCAL()){
     
-    Double_t vertex[] = {0,0,0};
-    if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
-    
-    first = kTRUE;
-    firstbkg = kTRUE;
-    
+       //Get vertex for photon momentum calculation
+       Double_t vertex[]  = {0,0,0} ; //vertex 
+       //Double_t vertex2[] = {0,0,0} ; //vertex of second input aod
+       if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
+       {
+               GetReader()->GetVertex(vertex);
+               //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
+       }
+         
     for(Int_t iclus = 0;iclus < (GetAODEMCAL())->GetEntriesFast() ; iclus ++ ){
-      AliAODCaloCluster * calo = (AliAODCaloCluster *) (GetAODEMCAL()->At(iclus)) ;
+      AliVCluster * calo = (AliVCluster *) (GetAODEMCAL()->At(iclus)) ;
       
       //Cluster selection, not charged
-      if(calo->GetNTracksMatched() > 0) continue ;
+      if(IsTrackMatched(calo)) continue ;
       
-      calo->GetMomentum(lv,vertex);
+         //Input from second AOD?
+         Int_t input = 0;
+//       if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1 ;
+//       else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) input = 1;
+//             
+         //Get Momentum vector, 
+         if     (input == 0) calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
+         //else if(input == 1) calo->GetMomentum(lv,vertex2);//Assume that come from vertex in straight line  
+               
       //Particles in jet 
       if(IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)){
-       
-       if(first) {
-         new (refclusters) TRefArray(TProcessID::GetProcessWithUID(calo)); 
-         first = kFALSE;
-       }
-       
+
        refclusters->Add(calo); 
        
        if(lv.Pt() > ptcut )  jet+=lv;
@@ -1527,10 +1539,6 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
       //Background around (phi_gamma-pi, eta_leading)
       else if(IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)){
        
-       if(firstbkg) {
-         new (refclustersbkg) TRefArray(TProcessID::GetProcessWithUID(calo)); 
-         firstbkg = kFALSE;
-       }
        
        refclustersbkg->Add(calo); 
        
@@ -1541,13 +1549,13 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorre
   
   //Add referenced clusters to AOD
   if(refclusters->GetEntriesFast() > 0 ){
-    refclusters->SetName(GetAODRefArrayName()+"Clusters");
-    particle->AddRefArray(refclusters);
+    refclusters->SetName(Form("%sClusters",GetAODObjArrayName().Data()));
+    particle->AddObjArray(refclusters);
   }
   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No clusters in jet cone\n");
   if(refclustersbkg->GetEntriesFast() > 0 ){
-    refclustersbkg->SetName(GetAODRefArrayName()+"ClustersBkg");
-    particle->AddRefArray(refclustersbkg);
+    refclustersbkg->SetName(Form("%sClustersBkg",GetAODObjArrayName().Data()));
+    particle->AddObjArray(refclustersbkg);
   }
   else if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background clusters in jet cone\n");
   
@@ -1576,10 +1584,10 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleC
   if(phil < 0) phil+=TMath::TwoPi();
   Double_t etal = pLeading.Eta();
   
-  TRefArray * refclusters = particle->GetRefArray(GetAODRefArrayName()+"Clusters");
-  TRefArray * reftracks   = particle->GetRefArray(GetAODRefArrayName()+"Tracks");
-  TRefArray * refclustersbkg = particle->GetRefArray(GetAODRefArrayName()+"ClustersBkg");
-  TRefArray * reftracksbkg   = particle->GetRefArray(GetAODRefArrayName()+"TracksBkg");
+  TObjArray * refclusters    = particle->GetObjArray(Form("Clusters%s"   ,GetAODObjArrayName().Data()));
+  TObjArray * reftracks      = particle->GetObjArray(Form("Tracks%s"     ,GetAODObjArrayName().Data()));
+  TObjArray * refclustersbkg = particle->GetObjArray(Form("ClustersBkg%s",GetAODObjArrayName().Data()));
+  TObjArray * reftracksbkg   = particle->GetObjArray(Form("TracksBkg%s"  ,GetAODObjArrayName().Data()));
   
   //Different pt cut for jet particles in different collisions systems
   Float_t ptcut = fJetPtThreshold;
@@ -1615,24 +1623,48 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleC
   //Add neutral particles to jet
   if(!fJetsOnlyInCTS && refclusters){
     
-    Double_t vertex[] = {0,0,0};
-    if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
-    
+       //Get vertex for photon momentum calculation
+       Double_t vertex[]  = {0,0,0} ; //vertex 
+       //Double_t vertex2[] = {0,0,0} ; //vertex of second input aod
+       if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) 
+       {
+               GetReader()->GetVertex(vertex);
+               //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
+       }
+         
     //Loop on jet particles
     if(refclusters){
       for(Int_t iclus = 0;iclus < refclusters->GetEntriesFast() ; iclus ++ ){
-       AliAODCaloCluster * calo = (AliAODCaloCluster *) refclusters->At(iclus) ;
-       calo->GetMomentum(lv,vertex); 
-       if(lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)) jet+=lv;   
+                 AliVCluster * calo = (AliVCluster *) refclusters->At(iclus) ;
+                 
+                 //Input from second AOD?
+                 Int_t input = 0;
+       //        if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1 ;
+//               else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) input = 1;
+                 
+                 //Get Momentum vector, 
+                 if     (input == 0) calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
+                 //else if(input == 1) calo->GetMomentum(lv,vertex2);//Assume that come from vertex in straight line  
+                 
+                 if(lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)) jet+=lv;   
       }//jet cluster loop
     }
     
     //Loop on background particles
     if(refclustersbkg){
       for(Int_t iclus = 0;iclus < refclustersbkg->GetEntriesFast() ; iclus ++ ){
-       AliAODCaloCluster * calo = (AliAODCaloCluster *) refclustersbkg->At(iclus) ;
-       calo->GetMomentum(lv,vertex);
-       if( lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)) bkg+=lv;
+                 AliVCluster * calo = (AliVCluster *) refclustersbkg->At(iclus) ;
+
+                 //Input from second AOD?
+                 Int_t input = 0;
+//               if     (particle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) input = 1 ;
+//               else if(particle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) input = 1;
+                 
+                 //Get Momentum vector, 
+                 if     (input == 0) calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
+                 //else if(input == 1) calo->GetMomentum(lv,vertex2);//Assume that come from vertex in straight line  
+                 
+                 if( lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)) bkg+=lv;
       }//background cluster loop 
     }
   }//clusters in jet
@@ -1649,40 +1681,40 @@ void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleC
 }
 
 //____________________________________________________________________________
-Bool_t  AliAnaParticleJetLeadingConeCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) const {
-  //Select cluster depending on its pid and acceptance selections
-  
-  //Skip matched clusters with tracks
-  if(calo->GetNTracksMatched() > 0) return kFALSE;
-  
-  //Check PID
-  calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
-  pdg = AliCaloPID::kPhoton;   
-  if(IsCaloPIDOn()){
-    //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
-    //or redo PID, recommended option for EMCal.
-    if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
-      pdg = GetCaloPID()->GetPdg("EMCAL",calo->PID(),mom.E());//PID with weights
-    else
-      pdg = GetCaloPID()->GetPdg("EMCAL",mom,calo);//PID recalculated
-    
-  //  if(GetDebug() > 3) printf("AliAnaParticleJetLeadingConeCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
-    //If it does not pass pid, skip
-    if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) 
-      return kFALSE ;
-  }//CaloPID
-  
-   //Check acceptance selection
-  if(IsFidutialCutOn()){
-    Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"EMCAL") ;
-    if(! in ) return kFALSE ;
-  }
-  
-  //if(GetDebug() > 3) printf("AliAnaParticleJetLeadingConeCorrelation::SelectCluster() - Cluster selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
-  
-  return kTRUE; 
-  
-}
+//Bool_t  AliAnaParticleJetLeadingConeCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {
+//  //Select cluster depending on its pid and acceptance selections
+//  
+//  //Skip matched clusters with tracks
+//  if(IsTrackMatched(calo)) return kFALSE;
+//  
+//  //Check PID
+//  calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
+//  pdg = AliCaloPID::kPhoton;   
+//  if(IsCaloPIDOn()){
+//    //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
+//    //or redo PID, recommended option for EMCal.
+//    if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
+//      pdg = GetCaloPID()->GetPdg("EMCAL",calo->GetPID(),mom.E());//PID with weights
+//    else
+//      pdg = GetCaloPID()->GetPdg("EMCAL",mom,calo);//PID recalculated
+//    
+//  //  if(GetDebug() > 3) printf("AliAnaParticleJetLeadingConeCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
+//    //If it does not pass pid, skip
+//    if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) 
+//      return kFALSE ;
+//  }//CaloPID
+//  
+//   //Check acceptance selection
+//  if(IsFiducialCutOn()){
+//    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"EMCAL") ;
+//    if(! in ) return kFALSE ;
+//  }
+//  
+//  //if(GetDebug() > 3) printf("AliAnaParticleJetLeadingConeCorrelation::SelectCluster() - Cluster selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
+//  
+//  return kTRUE; 
+//  
+//}
 
 //__________________________________________________________________
 void AliAnaParticleJetLeadingConeCorrelation::Print(const Option_t * opt) const