]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaGeneratorKine.cxx
Dalitz: New histogram for photon effiVsRadius
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaGeneratorKine.cxx
index 28dd79e406116244d4a9bd54d7796f72a0d4efa4..112309e51e16134e361bcd3484f7fc31d52a1f72 100755 (executable)
@@ -37,10 +37,14 @@ ClassImp(AliAnaGeneratorKine)
 //__________________________________________
 AliAnaGeneratorKine::AliAnaGeneratorKine() : 
 AliAnaCaloTrackCorrBaseClass(), 
+fTriggerDetector(),  fTriggerDetectorString(),
+fFidCutTrigger(0),
+fMinChargedPt(0),    fMinNeutralPt(0),
 fStack(0),
 fParton2(0),         fParton3(0), 
-fParton6(0),         fParton7(0),   
+fParton6(0),         fParton7(0),
 fJet6(),             fJet7(),
+fTrigger(),          fLVTmp(),
 fPtHard(0),
 fhPtHard(0),         fhPtParton(0),    fhPtJet(0),
 fhPtPartonPtHard(0), fhPtJetPtHard(0), fhPtJetPtParton(0),
@@ -53,38 +57,50 @@ fhPtPhoton(0),       fhPtPi0(0)
   
   for(Int_t i = 0; i < 4; i++)
   {
-    fhPtPhotonLeading[i]              = fhPtPi0Leading[i]              = 0; 
-    fhPtPhotonLeadingIsolated[i]      = fhPtPi0LeadingIsolated[i]      = 0; 
-    fhZHardPhotonLeading[i]           = fhZHardPi0Leading[i]           = 0;            
-    fhZHardPhotonLeadingIsolated[i]   = fhZHardPi0LeadingIsolated[i]   = 0; 
-    fhZPartonPhotonLeading[i]         = fhZPartonPi0Leading[i]         = 0;            
-    fhZPartonPhotonLeadingIsolated[i] = fhZPartonPi0LeadingIsolated[i] = 0; 
-    fhZJetPhotonLeading[i]            = fhZJetPi0Leading[i]            = 0;            
-    fhZJetPhotonLeadingIsolated[i]    = fhZJetPi0LeadingIsolated[i]    = 0; 
-    fhXEPhotonLeading[i]              = fhXEPi0Leading[i]              = 0;            
-    fhXEPhotonLeadingIsolated[i]      = fhXEPi0LeadingIsolated[i]      = 0; 
-    fhXEUEPhotonLeading[i]            = fhXEUEPi0Leading[i]            = 0;            
-    fhXEUEPhotonLeadingIsolated[i]    = fhXEUEPi0LeadingIsolated[i]    = 0; 
+    fhPtPhotonLeading[i]          = fhPtPi0Leading[i]          = 0;
+    fhPtPhotonLeadingSumPt[i]     = fhPtPi0LeadingSumPt[i]     = 0;
+    fhPtPhotonLeadingIsolated[i]  = fhPtPi0LeadingIsolated[i]  = 0;
+    for(Int_t j = 0; j < 2; j++)
+    {
+    fhZHardPhoton[j][i]           = fhZHardPi0[j][i]           = 0;            
+    fhZHardPhotonIsolated[j][i]   = fhZHardPi0Isolated[j][i]   = 0; 
+    fhZPartonPhoton[j][i]         = fhZPartonPi0[j][i]         = 0;            
+    fhZPartonPhotonIsolated[j][i] = fhZPartonPi0Isolated[j][i] = 0; 
+    fhZJetPhoton[j][i]            = fhZJetPi0[j][i]            = 0;            
+    fhZJetPhotonIsolated[j][i]    = fhZJetPi0Isolated[j][i]    = 0; 
+    fhXEPhoton[j][i]              = fhXEPi0[j][i]              = 0;            
+    fhXEPhotonIsolated[j][i]      = fhXEPi0Isolated[j][i]      = 0; 
+    fhXEUEPhoton[j][i]            = fhXEUEPi0[j][i]            = 0;            
+    fhXEUEPhotonIsolated[j][i]    = fhXEUEPi0Isolated[j][i]    = 0; 
+
+    fhPtPartonTypeNearPhoton[j][i]         = fhPtPartonTypeNearPi0[j][i]         = 0;            
+    fhPtPartonTypeNearPhotonIsolated[j][i] = fhPtPartonTypeNearPi0Isolated[j][i] = 0; 
+
+    fhPtPartonTypeAwayPhoton[j][i]         = fhPtPartonTypeAwayPi0[j][i]         = 0;            
+    fhPtPartonTypeAwayPhotonIsolated[j][i] = fhPtPartonTypeAwayPi0Isolated[j][i] = 0; 
+    }
   }
   
 }
 
-//_______________________________________________________________________________
-void  AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger, 
-                                                    const Int_t   indexTrig,                     
-                                                    const Int_t   pdgTrig, 
-                                                    const Bool_t  leading[4],
-                                                    const Bool_t  isolated[4], 
-                                                    Int_t & iparton )  
+//___________________________________________________________________________
+Bool_t  AliAnaGeneratorKine::CorrelateWithPartonOrJet(Int_t   indexTrig,
+                                                      Int_t   pdgTrig,
+                                                      Bool_t  leading[4],
+                                                      Bool_t  isolated[4],
+                                                      Int_t & iparton )  
 {
   //Correlate trigger with partons or jets, get z
   
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::CorrelateWithPartonOrJet() - Start \n");
+  
   //Get the index of the mother
   iparton =  (fStack->Particle(indexTrig))->GetFirstMother();
   TParticle * mother = fStack->Particle(iparton);
   while (iparton > 7) 
   {
     iparton   = mother->GetFirstMother();
+    if(iparton < 0) { printf("AliAnaGeneratorKine::CorrelateWithPartonOrJet() - Negative index, skip event\n"); return kFALSE; }
     mother = fStack->Particle(iparton);
   }
   
@@ -93,10 +109,10 @@ void  AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger
   if(iparton < 6)
   {
     //printf("This particle is not from hard process - pdg %d, parton index %d\n",pdgTrig, iparton);
-    return; 
+    return kFALSE
   }
   
-  Float_t ptTrig   = trigger.Pt(); 
+  Float_t ptTrig   = fTrigger.Pt(); 
   Float_t partonPt = fParton6->Pt();
   Float_t jetPt    = fJet6.Pt();
   if(iparton==7)
@@ -105,10 +121,63 @@ void  AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger
     jetPt    = fJet6.Pt();
   }
   
+  //Get id of parton in near and away side
   
-  fhPtPartonPtHard->Fill(fPtHard, partonPt/fPtHard);
-  fhPtJetPtHard   ->Fill(fPtHard, jetPt/fPtHard);
-  fhPtJetPtParton ->Fill(fPtHard, jetPt/partonPt);
+  Int_t away = -1;
+  Int_t near = -1;
+  Int_t nearPDG = -1;
+  Int_t awayPDG = -1;
+  
+  //printf("parton 6 pdg = %d, parton 7 pdg = %d\n",fParton6->GetPdgCode(),fParton7->GetPdgCode());
+  
+  if(iparton==6)
+  {
+    nearPDG = fParton6->GetPdgCode();
+    awayPDG = fParton7->GetPdgCode();
+  }
+  else 
+  {
+    nearPDG = fParton7->GetPdgCode();
+    awayPDG = fParton6->GetPdgCode();
+  }
+
+  if     (nearPDG == 22) near = 0;
+  else if(nearPDG == 21) near = 1;
+  else                   near = 2;
+  
+  if     (awayPDG == 22) away = 0;
+  else if(awayPDG == 21) away = 1;
+  else                   away = 2;
+  
+  for( Int_t i = 0; i < 4; i++ )
+  {
+    if(pdgTrig==111)
+    {
+      
+      fhPtPartonTypeNearPi0[leading[i]][i]->Fill(ptTrig,near);
+      fhPtPartonTypeAwayPi0[leading[i]][i]->Fill(ptTrig,away);
+      if(isolated[i])
+      {
+        fhPtPartonTypeNearPi0Isolated[leading[i]][i]->Fill(ptTrig,near);
+        fhPtPartonTypeAwayPi0Isolated[leading[i]][i]->Fill(ptTrig,away);
+      }
+      
+    }// pi0
+    else if(pdgTrig==22)
+    {
+      fhPtPartonTypeNearPhoton[leading[i]][i]->Fill(ptTrig,near);
+      fhPtPartonTypeAwayPhoton[leading[i]][i]->Fill(ptTrig,away);
+      if(isolated[i])
+      {
+        fhPtPartonTypeNearPhotonIsolated[leading[i]][i]->Fill(ptTrig,near);
+        fhPtPartonTypeAwayPhotonIsolated[leading[i]][i]->Fill(ptTrig,away);
+      }
+      
+    } // photon
+  } // conditions loop
+  
+  
+  // RATIOS
 
   Float_t zHard = ptTrig / fPtHard;
   Float_t zPart = ptTrig / partonPt;
@@ -116,44 +185,45 @@ void  AliAnaGeneratorKine::CorrelateWithPartonOrJet(const TLorentzVector trigger
 
   //if(zHard > 1 ) printf("*** Particle energy larger than pT hard z=%f\n",zHard); 
   
-  //printf("Z : hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
+  //printf("Z: hard %2.2f, parton %2.2f, jet %2.2f\n",zHard,zPart,zJet);
   
   for( Int_t i = 0; i < 4; i++ )
   {
     if(pdgTrig==111)
     {
-      if(leading[i])
-      { 
-        fhZHardPi0Leading[i]  ->Fill(ptTrig,zHard);
-        fhZPartonPi0Leading[i]->Fill(ptTrig,zPart);
-        fhZJetPi0Leading[i]   ->Fill(ptTrig,zJet );
-        
-        if(isolated[i])
-        {
-          fhZHardPi0LeadingIsolated[i]  ->Fill(ptTrig,zHard);
-          fhZPartonPi0LeadingIsolated[i]->Fill(ptTrig,zPart);
-          fhZJetPi0LeadingIsolated[i]   ->Fill(ptTrig,zJet);
-        }
+      fhZHardPi0[leading[i]][i]  ->Fill(ptTrig,zHard);
+      fhZPartonPi0[leading[i]][i]->Fill(ptTrig,zPart);
+      fhZJetPi0[leading[i]][i]   ->Fill(ptTrig,zJet );
+      
+      if(isolated[i])
+      {
+        fhZHardPi0Isolated[leading[i]][i]  ->Fill(ptTrig,zHard);
+        fhZPartonPi0Isolated[leading[i]][i]->Fill(ptTrig,zPart);
+        fhZJetPi0Isolated[leading[i]][i]   ->Fill(ptTrig,zJet);
       }
+      
     }// pi0
     else if(pdgTrig==22)
     {
-      if(leading[i])
-      { 
-        fhZHardPhotonLeading[i]  ->Fill(ptTrig,zHard);
-        fhZPartonPhotonLeading[i]->Fill(ptTrig,zPart);
-        fhZJetPhotonLeading[i]   ->Fill(ptTrig,zJet );
-        
-        if(isolated[i])
-        {
-          fhZHardPhotonLeadingIsolated[i]  ->Fill(ptTrig,zHard);
-          fhZPartonPhotonLeadingIsolated[i]->Fill(ptTrig,zPart);
-          fhZJetPhotonLeadingIsolated[i]   ->Fill(ptTrig,zJet);
-        }
-      }      
+      
+      fhZHardPhoton[leading[i]][i]  ->Fill(ptTrig,zHard);
+      fhZPartonPhoton[leading[i]][i]->Fill(ptTrig,zPart);
+      fhZJetPhoton[leading[i]][i]   ->Fill(ptTrig,zJet );
+      
+      if(isolated[i])
+      {
+        fhZHardPhotonIsolated[leading[i]][i]  ->Fill(ptTrig,zHard);
+        fhZPartonPhotonIsolated[leading[i]][i]->Fill(ptTrig,zPart);
+        fhZJetPhotonIsolated[leading[i]][i]   ->Fill(ptTrig,zJet);
+      }
+      
     } // photon
-  } // conditions loop  
-}  
+  } // conditions loop
+  
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::CorrelateWithPartonOrJet() - End TRUE \n");
+  
+  return kTRUE;
+}
 
 
 //____________________________________________________
@@ -164,233 +234,339 @@ TList *  AliAnaGeneratorKine::GetCreateOutputObjects()
   TList * outputContainer = new TList() ; 
   outputContainer->SetName("GenKineHistos") ; 
   
-  Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();  
-  Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  
-  Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin(); 
+  Int_t   nptbins    = GetHistogramRanges()->GetHistoPtBins();
+  Float_t ptmax      = GetHistogramRanges()->GetHistoPtMax();
+  Float_t ptmin      = GetHistogramRanges()->GetHistoPtMin();
+
+  Int_t   nptsumbins = GetHistogramRanges()->GetHistoNPtSumBins();
+  Float_t ptsummax   = GetHistogramRanges()->GetHistoPtSumMax();
+  Float_t ptsummin   = GetHistogramRanges()->GetHistoPtSumMin();
+
   
   fhPtHard  = new TH1F("hPtHard"," pt hard for selected triggers",nptbins,ptmin,ptmax); 
-  fhPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
+  fhPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
   outputContainer->Add(fhPtHard);
   
   fhPtParton  = new TH1F("hPtParton"," pt parton for selected triggers",nptbins,ptmin,ptmax); 
-  fhPtParton->SetXTitle("p_{T}^{parton} (GeV/c)");
+  fhPtParton->SetXTitle("#it{p}_{T}^{parton} (GeV/#it{c})");
   outputContainer->Add(fhPtParton);
   
   fhPtJet  = new TH1F("hPtJet"," pt jet for selected triggers",nptbins,ptmin,ptmax); 
-  fhPtJet->SetXTitle("p_{T}^{jet} (GeV/c)");
+  fhPtJet->SetXTitle("#it{p}_{T}^{jet} (GeV/#it{c})");
   outputContainer->Add(fhPtJet);
   
   fhPtPartonPtHard  = new TH2F("hPtPartonPtHard","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2); 
-  fhPtPartonPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
-  fhPtPartonPtHard->SetYTitle("p_{T}^{parton}/p_{T}^{hard}");
+  fhPtPartonPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
+  fhPtPartonPtHard->SetYTitle("#it{p}_{T}^{parton}/#it{p}_{T}^{hard}");
   outputContainer->Add(fhPtPartonPtHard);
   
   fhPtJetPtHard  = new TH2F("hPtJetPtHard","jet pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2); 
-  fhPtJetPtHard->SetXTitle("p_{T}^{hard} (GeV/c)");
-  fhPtJetPtHard->SetYTitle("p_{T}^{jet}/p_{T}^{hard}");
+  fhPtJetPtHard->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
+  fhPtJetPtHard->SetYTitle("#it{p}_{T}^{jet}/#it{p}_{T}^{hard}");
   outputContainer->Add(fhPtJetPtHard);
   
   fhPtJetPtParton  = new TH2F("hPtJetPtParton","parton pt / pt hard for selected triggers",nptbins,ptmin,ptmax,200,0,2); 
-  fhPtJetPtParton->SetXTitle("p_{T}^{hard} (GeV/c)");
-  fhPtJetPtParton->SetYTitle("p_{T}^{jet}/p_{T}^{parton}");
+  fhPtJetPtParton->SetXTitle("#it{p}_{T}^{hard} (GeV/#it{c})");
+  fhPtJetPtParton->SetYTitle("#it{p}_{T}^{jet}/#it{p}_{T}^{parton}");
   outputContainer->Add(fhPtJetPtParton);
   
-  
-  fhPtPhoton  = new TH1F("hPtPhoton","Input Photon",nptbins,ptmin,ptmax); 
-  fhPtPhoton->SetXTitle("p_{T} (GeV/c)");
+  fhPtPhoton  = new TH1F("hPtPhoton","Input #gamma",nptbins,ptmin,ptmax);
+  fhPtPhoton->SetXTitle("#it{p}_{T} (GeV/#it{c})");
   outputContainer->Add(fhPtPhoton);
 
-  fhPtPi0  = new TH1F("hPtPi0","Input Pi0",nptbins,ptmin,ptmax); 
-  fhPtPi0->SetXTitle("p_{T} (GeV/c)");
+  fhPtPi0  = new TH1F("hPtPi0","Input #pi^{0}",nptbins,ptmin,ptmax);
+  fhPtPi0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
   outputContainer->Add(fhPtPi0);
   
-  TString name[]  = {"","_EMC","_Photon","_EMC_Photon"};
-  TString title[] = {"",", neutral in EMCal",", neutral only photon like",", neutral in EMCal and only photon like"};
-
+  TString name   [] = {"","_EMC","_Photon","_EMC_Photon"};
+  TString title  [] = {"",", neutral in EMCal",", neutral only #gamma-like",", neutral in EMCal and only #gamma-like"};
+  TString leading[] = {"NotLeading","Leading"};
+  
   for(Int_t i = 0; i < 4; i++)
   {
     
     // Pt
     
     fhPtPhotonLeading[i]  = new TH1F(Form("hPtPhotonLeading%s",name[i].Data()),
-                                  Form("Photon : Leading of all particles%s",title[i].Data()),
-                                  nptbins,ptmin,ptmax); 
-    fhPtPhotonLeading[i]->SetXTitle("p_{T} (GeV/c)");
+                                     Form("#gamma: Leading of all particles%s",title[i].Data()),
+                                     nptbins,ptmin,ptmax);
+    fhPtPhotonLeading[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhPtPhotonLeading[i]);
     
     fhPtPi0Leading[i]  = new TH1F(Form("hPtPi0Leading%s",name[i].Data()),
-                               Form("Pi0 : Leading of all particles%s",title[i].Data()),
-                               nptbins,ptmin,ptmax); 
-    fhPtPi0Leading[i]->SetXTitle("p_{T} (GeV/c)");
-    outputContainer->Add(fhPtPi0Leading[i]);  
+                                  Form("#pi^{0}: Leading of all particles%s",title[i].Data()),
+                                  nptbins,ptmin,ptmax);
+    fhPtPi0Leading[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtPi0Leading[i]);
     
     fhPtPhotonLeadingIsolated[i]  = new TH1F(Form("hPtPhotonLeadingIsolated%s",name[i].Data()),
-                                     Form("Photon : Leading of all particles%s, isolated",title[i].Data()),
-                                     nptbins,ptmin,ptmax); 
-    fhPtPhotonLeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
+                                             Form("#gamma: Leading of all particles%s, isolated",title[i].Data()),
+                                             nptbins,ptmin,ptmax);
+    fhPtPhotonLeadingIsolated[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhPtPhotonLeadingIsolated[i]);
     
     fhPtPi0LeadingIsolated[i]  = new TH1F(Form("hPtPi0LeadingIsolated%s",name[i].Data()),
-                                  Form("Pi0 : Leading of all particles%s, isolated",title[i].Data()),
-                                  nptbins,ptmin,ptmax); 
-    fhPtPi0LeadingIsolated[i]->SetXTitle("p_{T} (GeV/c)");
-    outputContainer->Add(fhPtPi0LeadingIsolated[i]);  
-    
-    
-    
-    // zHard
-    
-    fhZHardPhotonLeading[i]  = new TH2F(Form("hZHardPhotonLeading%s",name[i].Data()),
-                                     Form("Z-Hard of Photon : Leading of all particles%s",title[i].Data()),
-                                     nptbins,ptmin,ptmax,200,0,2); 
-    fhZHardPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZHardPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZHardPhotonLeading[i]);
-    
-    fhZHardPi0Leading[i]  = new TH2F(Form("hZHardPi0Leading%s",name[i].Data()),
-                                  Form("Z-Hard of Pi0 : Leading of all particles%s",title[i].Data()),
-                                  nptbins,ptmin,ptmax,200,0,2); 
-    fhZHardPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZHardPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZHardPi0Leading[i]);  
-    
-    fhZHardPhotonLeadingIsolated[i]  = new TH2F(Form("hZHardPhotonLeadingIsolated%s",name[i].Data()),
-                                             Form("Z-Hard of Photon : Leading of all particles%s, isolated",title[i].Data()),
-                                             nptbins,ptmin,ptmax,200,0,2); 
-    fhZHardPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZHardPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZHardPhotonLeadingIsolated[i]);
-    
-    fhZHardPi0LeadingIsolated[i]  = new TH2F(Form("hZHardPi0LeadingIsolated%s",name[i].Data()),
-                                          Form("Z-Hard of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
-                                          nptbins,ptmin,ptmax,200,0,2); 
-    fhZHardPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZHardPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZHardPi0LeadingIsolated[i]);  
-    
-    // zHard
-    
-    fhZPartonPhotonLeading[i]  = new TH2F(Form("hZPartonPhotonLeading%s",name[i].Data()),
-                                        Form("Z-Parton of Photon : Leading of all particles%s",title[i].Data()),
-                                        nptbins,ptmin,ptmax,200,0,2); 
-    fhZPartonPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZPartonPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZPartonPhotonLeading[i]);
-    
-    fhZPartonPi0Leading[i]  = new TH2F(Form("hZPartonPi0Leading%s",name[i].Data()),
-                                     Form("Z-Parton of Pi0 : Leading of all particles%s",title[i].Data()),
-                                     nptbins,ptmin,ptmax,200,0,2); 
-    fhZPartonPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZPartonPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZPartonPi0Leading[i]);  
-    
-    fhZPartonPhotonLeadingIsolated[i]  = new TH2F(Form("hZPartonPhotonLeadingIsolated%s",name[i].Data()),
-                                                Form("Z-Parton of Photon : Leading of all particles%s, isolated",title[i].Data()),
-                                                nptbins,ptmin,ptmax,200,0,2); 
-    fhZPartonPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZPartonPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZPartonPhotonLeadingIsolated[i]);
-    
-    fhZPartonPi0LeadingIsolated[i]  = new TH2F(Form("hZPartonPi0LeadingIsolated%s",name[i].Data()),
-                                             Form("Z-Parton of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
-                                             nptbins,ptmin,ptmax,200,0,2); 
-    fhZPartonPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZPartonPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZPartonPi0LeadingIsolated[i]);  
-    
-    
-    // zJet
-    
-    fhZJetPhotonLeading[i]  = new TH2F(Form("hZJetPhotonLeading%s",name[i].Data()),
-                                        Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
-                                        nptbins,ptmin,ptmax,200,0,2); 
-    fhZJetPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZJetPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZJetPhotonLeading[i]);
-    
-    fhZJetPi0Leading[i]  = new TH2F(Form("hZJetPi0Leading%s",name[i].Data()),
-                                     Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
-                                     nptbins,ptmin,ptmax,200,0,2); 
-    fhZJetPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZJetPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZJetPi0Leading[i]);  
-    
-    fhZJetPhotonLeadingIsolated[i]  = new TH2F(Form("hZJetPhotonLeadingIsolated%s",name[i].Data()),
-                                                Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
-                                                nptbins,ptmin,ptmax,200,0,2); 
-    fhZJetPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZJetPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZJetPhotonLeadingIsolated[i]);
-    
-    fhZJetPi0LeadingIsolated[i]  = new TH2F(Form("hZJetPi0LeadingIsolated%s",name[i].Data()),
-                                             Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
-                                             nptbins,ptmin,ptmax,200,0,2); 
-    fhZJetPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhZJetPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhZJetPi0LeadingIsolated[i]);    
-    
-    
-    // XE
-    
-    fhXEPhotonLeading[i]  = new TH2F(Form("hXEPhotonLeading%s",name[i].Data()),
-                                       Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
-                                       nptbins,ptmin,ptmax,200,0,2); 
-    fhXEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEPhotonLeading[i]);
-    
-    fhXEPi0Leading[i]  = new TH2F(Form("hXEPi0Leading%s",name[i].Data()),
-                                    Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
-                                    nptbins,ptmin,ptmax,200,0,2); 
-    fhXEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEPi0Leading[i]);  
-    
-    fhXEPhotonLeadingIsolated[i]  = new TH2F(Form("hXEPhotonLeadingIsolated%s",name[i].Data()),
-                                               Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
-                                               nptbins,ptmin,ptmax,200,0,2); 
-    fhXEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEPhotonLeadingIsolated[i]);
-    
-    fhXEPi0LeadingIsolated[i]  = new TH2F(Form("hXEPi0LeadingIsolated%s",name[i].Data()),
-                                            Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
-                                            nptbins,ptmin,ptmax,200,0,2); 
-    fhXEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEPi0LeadingIsolated[i]);          
-    
-    
-    // XE from UE
-    
-    fhXEUEPhotonLeading[i]  = new TH2F(Form("hXEUEPhotonLeading%s",name[i].Data()),
-                                       Form("Z-Jet of Photon : Leading of all particles%s",title[i].Data()),
-                                       nptbins,ptmin,ptmax,200,0,2); 
-    fhXEUEPhotonLeading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEUEPhotonLeading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEUEPhotonLeading[i]);
-    
-    fhXEUEPi0Leading[i]  = new TH2F(Form("hXEUEPi0Leading%s",name[i].Data()),
-                                    Form("Z-Jet of Pi0 : Leading of all particles%s",title[i].Data()),
-                                    nptbins,ptmin,ptmax,200,0,2); 
-    fhXEUEPi0Leading[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEUEPi0Leading[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEUEPi0Leading[i]);  
-    
-    fhXEUEPhotonLeadingIsolated[i]  = new TH2F(Form("hXEUEPhotonLeadingIsolated%s",name[i].Data()),
-                                               Form("Z-Jet of Photon : Leading of all particles%s, isolated",title[i].Data()),
-                                               nptbins,ptmin,ptmax,200,0,2); 
-    fhXEUEPhotonLeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEUEPhotonLeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEUEPhotonLeadingIsolated[i]);
-    
-    fhXEUEPi0LeadingIsolated[i]  = new TH2F(Form("hXEUEPi0LeadingIsolated%s",name[i].Data()),
-                                            Form("Z-Jet of Pi0 : Leading of all particles%s, isolated",title[i].Data()),
-                                            nptbins,ptmin,ptmax,200,0,2); 
-    fhXEUEPi0LeadingIsolated[i]->SetYTitle("p_{T}^{particle}/p_{T}^{hard}");
-    fhXEUEPi0LeadingIsolated[i]->SetXTitle("p_{T}^{particle} (GeV/c)");
-    outputContainer->Add(fhXEUEPi0LeadingIsolated[i]);          
+                                          Form("#pi^{0}: Leading of all particles%s, isolated",title[i].Data()),
+                                          nptbins,ptmin,ptmax);
+    fhPtPi0LeadingIsolated[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtPi0LeadingIsolated[i]);
+    
+    fhPtPhotonLeadingSumPt[i]  = new TH2F(Form("hPtPhotonLeadingSumPt%s",name[i].Data()),
+                                     Form("#gamma: Leading of all particles%s",title[i].Data()),
+                                     nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+    fhPtPhotonLeadingSumPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhPtPhotonLeadingSumPt[i]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtPhotonLeadingSumPt[i]);
+    
+    fhPtPi0LeadingSumPt[i]  = new TH2F(Form("hPtPi0LeadingSumPt%s",name[i].Data()),
+                                  Form("#pi^{0}: Leading of all particles%s",title[i].Data()),
+                                  nptbins,ptmin,ptmax,nptsumbins,ptsummin,ptsummax);
+    fhPtPi0LeadingSumPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhPtPi0LeadingSumPt[i]->SetYTitle("#Sigma #it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtPi0LeadingSumPt[i]);
+
     
+    // Leading or not loop
+    for(Int_t j = 0; j < 2; j++)
+    {
+      // Near side parton
+      
+      fhPtPartonTypeNearPhoton[j][i]  = new TH2F(Form("hPtPartonTypeNearPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                                 Form("#gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                                 nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeNearPhoton[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeNearPhoton[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeNearPhoton[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeNearPhoton[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeNearPhoton[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeNearPhoton[j][i]);
+      
+      fhPtPartonTypeNearPi0[j][i]  = new TH2F(Form("hPtPartonTypeNearPi0%s%s",leading[j].Data(),name[i].Data()),
+                                              Form("#pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                              nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeNearPi0[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeNearPi0[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeNearPi0[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeNearPi0[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeNearPi0[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeNearPi0[j][i]);
+      
+      fhPtPartonTypeNearPhotonIsolated[j][i]  = new TH2F(Form("hPtPartonTypeNearPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                                         Form("#gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                                         nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeNearPhotonIsolated[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeNearPhotonIsolated[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeNearPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeNearPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeNearPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeNearPhotonIsolated[j][i]);
+      
+      fhPtPartonTypeNearPi0Isolated[j][i]  = new TH2F(Form("hPtPartonTypeNearPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                                      Form("#pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                                      nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeNearPi0Isolated[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeNearPi0Isolated[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeNearPi0Isolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeNearPi0Isolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeNearPi0Isolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeNearPi0Isolated[j][i]);
+      
+      
+      // Away side parton
+      
+      fhPtPartonTypeAwayPhoton[j][i]  = new TH2F(Form("hPtPartonTypeAwayPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                                 Form("#gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                                 nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeAwayPhoton[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeAwayPhoton[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeAwayPhoton[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeAwayPhoton[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeAwayPhoton[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeAwayPhoton[j][i]);
+      
+      fhPtPartonTypeAwayPi0[j][i]  = new TH2F(Form("hPtPartonTypeAwayPi0%s%s",leading[j].Data(),name[i].Data()),
+                                              Form("#pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                              nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeAwayPi0[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeAwayPi0[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeAwayPi0[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeAwayPi0[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeAwayPi0[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeAwayPi0[j][i]);
+      
+      fhPtPartonTypeAwayPhotonIsolated[j][i]  = new TH2F(Form("hPtPartonTypeAwayPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                                         Form("#gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                                         nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeAwayPhotonIsolated[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeAwayPhotonIsolated[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeAwayPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeAwayPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeAwayPhotonIsolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeAwayPhotonIsolated[j][i]);
+      
+      fhPtPartonTypeAwayPi0Isolated[j][i]  = new TH2F(Form("hPtPartonTypeAwayPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                                      Form("#pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                                      nptbins,ptmin,ptmax,3,0,3);
+      fhPtPartonTypeAwayPi0Isolated[j][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtPartonTypeAwayPi0Isolated[j][i]->SetYTitle("Parton type");
+      fhPtPartonTypeAwayPi0Isolated[j][i]->GetYaxis()->SetBinLabel(1,"#gamma");
+      fhPtPartonTypeAwayPi0Isolated[j][i]->GetYaxis()->SetBinLabel(2,"g");
+      fhPtPartonTypeAwayPi0Isolated[j][i]->GetYaxis()->SetBinLabel(3,"q");
+      outputContainer->Add(fhPtPartonTypeAwayPi0Isolated[j][i]);
+      
+      // zHard
+      
+      fhZHardPhoton[j][i]  = new TH2F(Form("hZHardPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                      Form("#it{z}_{Hard} of #gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                      nptbins,ptmin,ptmax,200,0,2);
+      fhZHardPhoton[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZHardPhoton[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZHardPhoton[j][i]);
+      
+      fhZHardPi0[j][i]  = new TH2F(Form("hZHardPi0%s%s",leading[j].Data(),name[i].Data()),
+                                   Form("#it{z}_{Hard} of #pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                   nptbins,ptmin,ptmax,200,0,2);
+      fhZHardPi0[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZHardPi0[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZHardPi0[j][i]);
+      
+      fhZHardPhotonIsolated[j][i]  = new TH2F(Form("hZHardPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                              Form("#it{z}_{Hard} of #gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                              nptbins,ptmin,ptmax,200,0,2);
+      fhZHardPhotonIsolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZHardPhotonIsolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZHardPhotonIsolated[j][i]);
+      
+      fhZHardPi0Isolated[j][i]  = new TH2F(Form("hZHardPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                           Form("#it{z}_{Hard} of #pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                           nptbins,ptmin,ptmax,200,0,2);
+      fhZHardPi0Isolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZHardPi0Isolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZHardPi0Isolated[j][i]);
+      
+      // zHard
+      
+      fhZPartonPhoton[j][i]  = new TH2F(Form("hZPartonPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                        Form("#it{z}_{Parton} of #gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                        nptbins,ptmin,ptmax,200,0,2);
+      fhZPartonPhoton[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZPartonPhoton[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZPartonPhoton[j][i]);
+      
+      fhZPartonPi0[j][i]  = new TH2F(Form("hZPartonPi0%s%s",leading[j].Data(),name[i].Data()),
+                                     Form("#it{z}_{Parton} of #pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                     nptbins,ptmin,ptmax,200,0,2);
+      fhZPartonPi0[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZPartonPi0[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZPartonPi0[j][i]);
+      
+      fhZPartonPhotonIsolated[j][i]  = new TH2F(Form("hZPartonPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                                Form("#it{z}_{Parton} of #gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                                nptbins,ptmin,ptmax,200,0,2);
+      fhZPartonPhotonIsolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZPartonPhotonIsolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZPartonPhotonIsolated[j][i]);
+      
+      fhZPartonPi0Isolated[j][i]  = new TH2F(Form("hZPartonPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                             Form("#it{z}_{Parton} of #pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                             nptbins,ptmin,ptmax,200,0,2);
+      fhZPartonPi0Isolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZPartonPi0Isolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZPartonPi0Isolated[j][i]);
+      
+      
+      // zJet
+      
+      fhZJetPhoton[j][i]  = new TH2F(Form("hZJetPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                     Form("#it{z}_{Jet} of #gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                     nptbins,ptmin,ptmax,200,0,2);
+      fhZJetPhoton[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZJetPhoton[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZJetPhoton[j][i]);
+      
+      fhZJetPi0[j][i]  = new TH2F(Form("hZJetPi0%s%s",leading[j].Data(),name[i].Data()),
+                                  Form("#it{z}_{Jet} of #pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                  nptbins,ptmin,ptmax,200,0,2);
+      fhZJetPi0[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZJetPi0[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZJetPi0[j][i]);
+      
+      fhZJetPhotonIsolated[j][i]  = new TH2F(Form("hZJetPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                             Form("#it{z}_{Jet} of #gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                             nptbins,ptmin,ptmax,200,0,2);
+      fhZJetPhotonIsolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZJetPhotonIsolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZJetPhotonIsolated[j][i]);
+      
+      fhZJetPi0Isolated[j][i]  = new TH2F(Form("hZJetPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                          Form("#it{z}_{Jet} of #pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                          nptbins,ptmin,ptmax,200,0,2);
+      fhZJetPi0Isolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhZJetPi0Isolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhZJetPi0Isolated[j][i]);
+      
+      
+      // XE
+      
+      fhXEPhoton[j][i]  = new TH2F(Form("hXEPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                   Form("#it{z}_{Jet} of #gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                   nptbins,ptmin,ptmax,200,0,2);
+      fhXEPhoton[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEPhoton[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEPhoton[j][i]);
+      
+      fhXEPi0[j][i]  = new TH2F(Form("hXEPi0%s%s",leading[j].Data(),name[i].Data()),
+                                Form("#it{z}_{Jet} of #pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                nptbins,ptmin,ptmax,200,0,2);
+      fhXEPi0[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEPi0[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEPi0[j][i]);
+      
+      fhXEPhotonIsolated[j][i]  = new TH2F(Form("hXEPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                           Form("#it{z}_{Jet} of #gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                           nptbins,ptmin,ptmax,200,0,2);
+      fhXEPhotonIsolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEPhotonIsolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEPhotonIsolated[j][i]);
+      
+      fhXEPi0Isolated[j][i]  = new TH2F(Form("hXEPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                        Form("#it{z}_{Jet} of #pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                        nptbins,ptmin,ptmax,200,0,2);
+      fhXEPi0Isolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEPi0Isolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEPi0Isolated[j][i]);
+      
+      
+      // XE from UE
+      
+      fhXEUEPhoton[j][i]  = new TH2F(Form("hXEUEPhoton%s%s",leading[j].Data(),name[i].Data()),
+                                     Form("#it{z}_{Jet} of #gamma: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                     nptbins,ptmin,ptmax,200,0,2);
+      fhXEUEPhoton[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEUEPhoton[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEUEPhoton[j][i]);
+      
+      fhXEUEPi0[j][i]  = new TH2F(Form("hXEUEPi0%s%s",leading[j].Data(),name[i].Data()),
+                                  Form("#it{z}_{Jet} of #pi^{0}: %s of all particles%s",leading[j].Data(),title[i].Data()),
+                                  nptbins,ptmin,ptmax,200,0,2);
+      fhXEUEPi0[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEUEPi0[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEUEPi0[j][i]);
+      
+      fhXEUEPhotonIsolated[j][i]  = new TH2F(Form("hXEUEPhoton%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                             Form("#it{z}_{Jet} of #gamma: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                             nptbins,ptmin,ptmax,200,0,2);
+      fhXEUEPhotonIsolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEUEPhotonIsolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEUEPhotonIsolated[j][i]);
+      
+      fhXEUEPi0Isolated[j][i]  = new TH2F(Form("hXEUEPi0%sIsolated%s",leading[j].Data(),name[i].Data()),
+                                          Form("#it{z}_{Jet} of #pi^{0}: %s of all particles%s, isolated",leading[j].Data(),title[i].Data()),
+                                          nptbins,ptmin,ptmax,200,0,2); 
+      fhXEUEPi0Isolated[j][i]->SetYTitle("#it{p}_{T}^{particle}/#it{p}_{T}^{hard}");
+      fhXEUEPi0Isolated[j][i]->SetXTitle("#it{p}_{T}^{particle} (GeV/#it{c})");
+      outputContainer->Add(fhXEUEPi0Isolated[j][i]);          
+    }
   }
   
   return outputContainer;
@@ -402,6 +578,8 @@ void  AliAnaGeneratorKine::GetPartonsAndJets()
 {
   // Fill data members with partons,jets and generated pt hard 
   
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::GetPartonsAndJets() - Start \n");
+
   fStack =  GetMCStack() ;
   
   if(!fStack) 
@@ -444,25 +622,25 @@ void  AliAnaGeneratorKine::GetPartonsAndJets()
     {
       pygeh->TriggerJet(ijet, tmpjet);
       
-      TLorentzVector jet(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
-      Float_t jphi = jet.Phi();
+      fLVTmp.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
+      Float_t jphi = fLVTmp.Phi();
       if(jphi < 0) jphi +=TMath::TwoPi();
       
-      Double_t radius6 = GetIsolationCut()->Radius(fParton6->Eta(), p6phi, jet.Eta() , jphi) ;
-      Double_t radius7 = GetIsolationCut()->Radius(fParton7->Eta(), p7phi, jet.Eta() , jphi) ;
+      Double_t radius6 = GetIsolationCut()->Radius(fParton6->Eta(), p6phi, fLVTmp.Eta() , jphi) ;
+      Double_t radius7 = GetIsolationCut()->Radius(fParton7->Eta(), p7phi, fLVTmp.Eta() , jphi) ;
       
       //printf("jet %d: pt %2.2f, eta %2.2f, phi %2.2f, r6 %2.2f, r7 %2.2f\n",ijet,jet.Pt(),jet.Eta(),jphi,radius6, radius7);
       
       if (radius6 < jet6R)
       {
         jet6R = radius6;
-        fJet6 = jet;
+        fJet6 = fLVTmp;
         
       }
       if (radius7 < jet7R) 
       {
         jet7R = radius7;
-        fJet7 = jet;
+        fJet7 = fLVTmp;
       }
             
     } // jet loop
@@ -478,30 +656,36 @@ void  AliAnaGeneratorKine::GetPartonsAndJets()
   fhPtParton ->Fill(fParton6->Pt());
   fhPtParton ->Fill(fParton7->Pt());
 
+  fhPtPartonPtHard->Fill(fPtHard, fParton6->Pt()/fPtHard);
+  fhPtPartonPtHard->Fill(fPtHard, fParton7->Pt()/fPtHard);
+  fhPtJetPtHard   ->Fill(fPtHard, fJet6.Pt()/fPtHard);
+  fhPtJetPtHard   ->Fill(fPtHard, fJet7.Pt()/fPtHard);
+  fhPtJetPtParton ->Fill(fPtHard, fJet6.Pt()/fParton6->Pt());
+  fhPtJetPtParton ->Fill(fPtHard, fJet7.Pt()/fParton7->Pt());
+  
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::GetPartonsAndJets() - End \n");
+
 }
 
-//___________________________________________________________
-void AliAnaGeneratorKine::GetXE(const TLorentzVector trigger,  
-                                const Int_t   indexTrig,                     
-                                const Int_t   pdgTrig, 
-                                const Bool_t  leading[4], 
-                                const Bool_t  isolated[4], 
-                                const Int_t   iparton)     
+//_____________________________________________________
+void AliAnaGeneratorKine::GetXE(Int_t   indexTrig,
+                                Int_t   pdgTrig,
+                                Bool_t  leading[4],
+                                Bool_t  isolated[4],
+                                Int_t   iparton)
 {
 
   // Calculate the real XE and the UE XE
 
-  Float_t ptThresTrack = 0.2;
-
-  Float_t ptTrig  = trigger.Pt();
-  Float_t etaTrig = trigger.Eta();
-  Float_t phiTrig = trigger.Phi();
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::GetXE() - Start \n");
+  
+  Float_t ptTrig  = fTrigger.Pt();
+  Float_t phiTrig = fTrigger.Phi();
   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
   
   //Loop on primaries, start from position 8, no partons
   for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
   {
-    
     TParticle * particle = fStack->Particle(ipr) ;
     
     if(ipr==indexTrig) continue;
@@ -510,105 +694,100 @@ void AliAnaGeneratorKine::GetXE(const TLorentzVector trigger,
     Int_t   status = particle->GetStatusCode();
         
     // Compare trigger with final state particles
-    if( status != 1) continue ;
+    if( status != 1 ) continue ;
     
     Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
     
-    if(charge==0) continue; // construct xe only with charged        
+    if(charge == 0 ) continue; // construct xe only with charged
     
     Float_t pt     = particle->Pt();
-    Float_t eta    = particle->Eta();
     Float_t phi    = particle->Phi();
     if(phi < 0 ) phi += TMath::TwoPi();
     
-    if( pt < ptThresTrack)    continue ;
+    if( pt < fMinChargedPt)    continue ;
     
-    if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
+    particle->Momentum(fLVTmp);
+    Bool_t inTPC = GetFiducialCut()->IsInFiducialCut(fLVTmp.Eta(),fLVTmp.Phi(),kCTS) ;
     
-    //Isolation
-    Double_t radius = GetIsolationCut()->Radius(etaTrig, phiTrig, eta , phi) ;
+    if(!inTPC) continue;
     
     Float_t xe = -pt/ptTrig*TMath::Cos(phi-phiTrig);
     
-    //Get the index of the mother
+    // ---------------------------------------------------
+    // Get the index of the mother, get from what parton
     Int_t ipartonAway =  particle->GetFirstMother();
+    if(ipartonAway < 0) return;
+    
     TParticle * mother = fStack->Particle(ipartonAway);
     while (ipartonAway > 7) 
     {
       ipartonAway   = mother->GetFirstMother();
+      if(ipartonAway < 0) break;
       mother = fStack->Particle(ipartonAway);
     }
     
+    //-----------------------------------------
+    // Get XE of particles belonging to the jet
+    // on the opposite side of the trigger
     if((ipartonAway==6 || ipartonAway==7) && iparton!=ipartonAway) 
     {
-      //printf("xE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
-      if(radius > 1 ) continue; // avoid particles too far from trigger
-      
       for( Int_t i = 0; i < 4; i++ )
       {
         if(pdgTrig==111)
         {
-          if(leading[i])
-          { 
-            fhXEPi0Leading[i]  ->Fill(ptTrig,xe);
-            
-            if(isolated[i])
-            {
-              fhXEPi0LeadingIsolated[i]  ->Fill(ptTrig,xe);
-            }
+          fhXEPi0[leading[i]][i]  ->Fill(ptTrig,xe);
+          
+          if(isolated[i])
+          {
+            fhXEPi0Isolated[leading[i]][i]  ->Fill(ptTrig,xe);
           }
+          
         }// pi0
         else if(pdgTrig==22)
         {
-          if(leading[i])
-          { 
-            fhXEPhotonLeading[i]  ->Fill(ptTrig,xe);
-            
-            if(isolated[i])
-            {
-              fhXEPhotonLeadingIsolated[i]  ->Fill(ptTrig,xe);
-            }
-          }      
+          fhXEPhoton[leading[i]][i]  ->Fill(ptTrig,xe);
+          
+          if(isolated[i])
+          {
+            fhXEPhotonIsolated[leading[i]][i]  ->Fill(ptTrig,xe);
+          }
+          
         } // photon
-      } // conditions loop  
+      } // conditions loop
     } // Away side
 
-    if(ipartonAway!=6 && ipartonAway!=7) 
+    //----------------------------------------------------------
+    // Get the XE from particles not attached to any of the jets
+    if(ipartonAway!=6 && ipartonAway!=7)
     {
-      
-      //printf("xE UE : iparton %d, ipartonAway %d\n",iparton,ipartonAway);
-
       for( Int_t i = 0; i < 4; i++ )
       {
         if(pdgTrig==111)
         {
-          if(leading[i])
-          { 
-            fhXEUEPi0Leading[i]  ->Fill(ptTrig,xe);
-            
-            if(isolated[i])
-            {
-              fhXEUEPi0LeadingIsolated[i]  ->Fill(ptTrig,xe);
-            }
+          fhXEUEPi0[leading[i]][i]  ->Fill(ptTrig,xe);
+          
+          if(isolated[i])
+          {
+            fhXEUEPi0Isolated[leading[i]][i]  ->Fill(ptTrig,xe);
           }
+          
         }// pi0
         else if(pdgTrig==22)
         {
-          if(leading[i])
-          { 
-            fhXEUEPhotonLeading[i]  ->Fill(ptTrig,xe);
-            
-            if(isolated[i])
-            {
-              fhXEUEPhotonLeadingIsolated[i]  ->Fill(ptTrig,xe);
-            }
-          }      
+          fhXEUEPhoton[leading[i]][i]  ->Fill(ptTrig,xe);
+          
+          if(isolated[i])
+          {
+            fhXEUEPhotonIsolated[leading[i]][i]  ->Fill(ptTrig,xe);
+          }
+          
         } // photon
       } // conditions loop  
     } // Away side    
     
   } // primary loop
 
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::GetPartonsAndJets() - End \n");
 
 }
 
@@ -619,18 +798,24 @@ void AliAnaGeneratorKine::InitParameters()
   //Initialize the parameters of the analysis.
   AddToHistogramsName("AnaGenKine_");
   
+  fTriggerDetector = kEMCAL;
+  
+  fMinChargedPt    = 0.2;
+  fMinNeutralPt    = 0.3;
+  
 }
 
-//___________________________________________________________________________
-void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger, 
-                                                const Int_t indexTrig, 
-                                                const Int_t pdgTrig, 
+//_____________________________________________________________________
+void  AliAnaGeneratorKine::IsLeadingAndIsolated(Int_t indexTrig,
+                                                Int_t pdgTrig,
                                                 Bool_t leading[4],
                                                 Bool_t isolated[4]) 
 {
-  // Check if the trigger is the leading particle
+  // Check if the trigger is the leading particle and if it is isolated
   // In case of neutral particles check all neutral or neutral in EMCAL acceptance
   
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::GetIsLeadingAndIsolated() - Start \n");
+
   Float_t ptMaxCharged       = 0; // all charged
   Float_t ptMaxNeutral       = 0; // all neutral
   Float_t ptMaxNeutEMCAL     = 0; // for neutral, select them in EMCAL acceptance
@@ -647,24 +832,33 @@ void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
   isolated[2] = 0;
   isolated[3] = 0;
   
-  Float_t ptTrig  = trigger.Pt();
-  Float_t etaTrig = trigger.Eta();
-  Float_t phiTrig = trigger.Phi();
+  Float_t ptTrig  = fTrigger.Pt();
+  Float_t etaTrig = fTrigger.Eta();
+  Float_t phiTrig = fTrigger.Phi();
   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
 
   // Minimum track or cluster energy
-  Float_t ptThresTrack = 0.2;
-  Float_t ptThresCalo  = 0.3;
 
   //Isolation cuts
-  Float_t ptThresIC    = 0.5;
-  Float_t rThresIC     = 0.4;
+  Float_t ptThresIC    = GetIsolationCut()->GetPtThreshold();
+  Float_t sumThresIC   = GetIsolationCut()->GetPtThreshold();
+  Float_t rThresIC     = GetIsolationCut()->GetConeSize();
+  Float_t isoMethod    = GetIsolationCut()->GetICMethod();
+  
+  //Counters
   Int_t   nICTrack     = 0;
   Int_t   nICNeutral   = 0;
   Int_t   nICNeutEMCAL = 0;
   Int_t   nICNeutPhot  = 0;
   Int_t   nICNeutEMCALPhot = 0;
   
+  // Sum of pT
+  Float_t sumNePt         = 0;
+  Float_t sumNePtPhot     = 0;
+  Float_t sumNePtEMC      = 0;
+  Float_t sumNePtEMCPhot  = 0;
+  Float_t sumChPt         = 0;
+  
   //Loop on primaries, start from position 8, no partons
   for(Int_t ipr = 8; ipr < fStack->GetNprimary(); ipr ++ )
   {
@@ -681,14 +875,16 @@ void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
      
     // Compare trigger with final state particles
     if( status != 1) continue ;
+
+    // Select all particles in at least the TPC acceptance
+    Bool_t inTPC = GetFiducialCut()->IsInFiducialCut(fTrigger.Eta(),fTrigger.Phi(),kCTS) ;
+    if(!inTPC) continue;
     
     Float_t pt     = particle->Pt();
     Float_t eta    = particle->Eta();
     Float_t phi    = particle->Phi();
     if(phi < 0 ) phi += TMath::TwoPi();
     
-    if(TMath::Abs(eta) > 0.8) continue ; // TPC acceptance cut
-
     Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
     
     //Isolation
@@ -696,50 +892,64 @@ void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
     
     if(charge==0)
     {
-      if(pt < ptThresCalo)  continue ;
+      if(pt < fMinNeutralPt)  continue ;
+      
+      if( ptMaxNeutral < pt ) ptMaxNeutral = pt;
+      
+      if( radius < rThresIC )
+      {
+        if( pt > ptThresIC ) nICNeutral++ ;
+        sumNePt+= pt;
+      }
       
-      if( ptMaxNeutral < pt )                   ptMaxNeutral = pt;
-      if( pt > ptThresIC && radius < rThresIC ) nICNeutral++ ;
-
       Bool_t phPDG = kFALSE;
       if(pdg==22 || pdg==111) phPDG = kTRUE;
     
       //if(pt > ptTrig) printf(" --- pdg %d, phPDG %d pT %2.2f, pTtrig %2.2f, eta %2.2f, phi %2.2f\n",pdg,phPDG,pt,ptTrig,particle->Eta(), particle->Phi()*TMath::RadToDeg());
       if(phPDG)
       {
-        if( ptMaxNeutPhot < pt)                   ptMaxNeutPhot = pt;
-        if( pt > ptThresIC && radius < rThresIC ) nICNeutPhot++ ;
+        if( ptMaxNeutPhot < pt) ptMaxNeutPhot = pt;
+        
+        if( radius < rThresIC )
+        {
+          if(pt > ptThresIC) nICNeutPhot++ ;
+          sumNePtPhot += pt;
+        }
       }
       
-      //EMCAL acceptance
-      Bool_t inEMCAL = kTRUE;
-      if(TMath::Abs(particle->Eta()) > 0.7
-         || particle->Phi() > TMath::DegToRad()*180
-         || particle->Phi() < TMath::DegToRad()*80 ) inEMCAL = kFALSE ;
+      //Calorimeter acceptance
+      Bool_t inCalo = GetFiducialCut()->IsInFiducialCut(fTrigger.Eta(),fTrigger.Phi(),GetCalorimeter()) ;
+      if(!inCalo) continue;
       
-      if(inEMCAL)
+      if( ptMaxNeutEMCAL < pt ) ptMaxNeutEMCAL = pt;
+      if( radius < rThresIC )
       {
-        if( ptMaxNeutEMCAL < pt )                 ptMaxNeutEMCAL = pt;
-        if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCAL++ ;
-
-        if(phPDG)
+        if( pt > ptThresIC ) nICNeutEMCAL++ ;
+        sumNePtEMC += pt;
+      }
+      
+      if(phPDG)
+      {
+        if( ptMaxNeutEMCALPhot < pt )             ptMaxNeutEMCALPhot = pt;
+        if(  radius < rThresIC )
         {
-          if( ptMaxNeutEMCALPhot < pt )             ptMaxNeutEMCALPhot = pt;
-          if( pt > ptThresIC && radius < rThresIC ) nICNeutEMCALPhot++ ;
+          if (pt > ptThresIC) nICNeutEMCALPhot++ ;
+          sumNePtEMCPhot += pt;
         }
       }
     }
-    else  
+    else
     {
-      if( pt < ptThresTrack)  continue ;
-
+      if( pt < fMinChargedPt)  continue ;
+      
       if( ptMaxCharged < pt )   ptMaxCharged   = pt;
       
-      if( pt > ptThresIC && radius < rThresIC )  
+      if( radius < rThresIC )
       {
         //printf("UE track? pTtrig %2.2f, pt %2.2f, etaTrig %2.2f,  eta %2.2f, phiTrig %2.2f,  phi %2.2f, radius %2.2f\n",
         //       ptTrig, pt,etaTrig, eta, phiTrig*TMath::RadToDeg(), phi*TMath::RadToDeg(),radius);
-        nICTrack++ ;
+        if( pt > ptThresIC ) nICTrack++ ;
+        sumChPt += pt;
       }
     }
 
@@ -756,18 +966,31 @@ void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
     if(ptTrig > ptMaxNeutEMCALPhot) leading[3] = kTRUE ;
   }
   
-  //printf("N in cone over threshold : tracks  %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n", 
+  //printf("N in cone over threshold: tracks  %d, neutral %d, neutral emcal %d, photon %d, photon emcal %d\n", 
   //       nICTrack, nICNeutral ,nICNeutEMCAL,nICNeutPhot, nICNeutEMCALPhot);
   
+  //------------------
   // Isolation decision
-  if(nICTrack == 0)
+  if( isoMethod == AliIsolationCut::kPtThresIC )
   {
-    if(nICNeutral       == 0 ) isolated[0] = kTRUE ;
-    if(nICNeutEMCAL     == 0 ) isolated[1] = kTRUE ;
-    if(nICNeutPhot      == 0 ) isolated[2] = kTRUE ;
-    if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
+    if( nICTrack == 0 )
+    {
+      if(nICNeutral       == 0 ) isolated[0] = kTRUE ;
+      if(nICNeutEMCAL     == 0 ) isolated[1] = kTRUE ;
+      if(nICNeutPhot      == 0 ) isolated[2] = kTRUE ;
+      if(nICNeutEMCALPhot == 0 ) isolated[3] = kTRUE ;
+    }
   }
+  else if( isoMethod == AliIsolationCut::kSumPtIC )
+  {
+    if(sumChPt + sumNePt        < sumThresIC ) isolated[0] = kTRUE ;
+    if(sumChPt + sumNePtEMC     < sumThresIC ) isolated[1] = kTRUE ;
+    if(sumChPt + sumNePtPhot    < sumThresIC ) isolated[2] = kTRUE ;
+    if(sumChPt + sumNePtEMCPhot < sumThresIC ) isolated[3] = kTRUE ;
+  }
+  
   
+  //----------------------------------------------------
   // Fill histograms if conditions apply for all 4 cases
   for( Int_t i = 0; i < 4; i++ )
   {
@@ -776,6 +999,12 @@ void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
       if(leading[i])
       { 
         fhPtPi0Leading[i]->Fill(ptTrig);
+        
+        if     (i == 0) fhPtPi0LeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePt);
+        else if(i == 1) fhPtPi0LeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePtEMC);
+        else if(i == 2) fhPtPi0LeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePtPhot);
+        else if(i == 3) fhPtPi0LeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePtEMCPhot);
+
         if(isolated[i]) fhPtPi0LeadingIsolated[i]->Fill(ptTrig);
       }
     }// pi0
@@ -784,11 +1013,19 @@ void  AliAnaGeneratorKine::IsLeadingAndIsolated(const TLorentzVector trigger,
       if(leading[i])
       { 
         fhPtPhotonLeading[i]->Fill(ptTrig);
+        
+        if     (i == 0) fhPtPhotonLeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePt);
+        else if(i == 1) fhPtPhotonLeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePtEMC);
+        else if(i == 2) fhPtPhotonLeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePtPhot);
+        else if(i == 3) fhPtPhotonLeadingSumPt[i]->Fill(ptTrig, sumChPt + sumNePtEMCPhot);
+        
         if(isolated[i]) fhPtPhotonLeadingIsolated[i]->Fill(ptTrig);
       }      
     } // photon
   } // conditions loop
-    
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::IsLeadingAndIsolated() - End \n");
+  
 }
   
 //_____________________________________________________
@@ -796,7 +1033,7 @@ void  AliAnaGeneratorKine::MakeAnalysisFillHistograms()
 {
   //Particle-Parton Correlation Analysis, fill histograms
   
-  TLorentzVector trigger;
+  if(GetDebug() > 1) printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - Start \n");
   
   GetPartonsAndJets();
   
@@ -820,21 +1057,21 @@ void  AliAnaGeneratorKine::MakeAnalysisFillHistograms()
     else if(pdgTrig != 111) continue;
     
     // Acceptance and kinematical cuts
-    if( ptTrig < 8 )    continue ;
+    if( ptTrig < GetMinPt() ) continue ;
+    
+    // Recover the kinematics:
+    particle->Momentum(fTrigger);
     
-    //EMCAL acceptance, a bit less
-    if(TMath::Abs(particle->Eta()) > 0.6) continue ;
-    if(particle->Phi() > TMath::DegToRad()*176) continue ;
-    if(particle->Phi() < TMath::DegToRad()*74 ) continue ;
+    Bool_t in = GetFiducialCutForTrigger()->IsInFiducialCut(fTrigger.Eta(),fTrigger.Phi(),fTriggerDetector) ;
+    if(! in ) continue ;
 
-//    printf("Particle %d : pdg %d status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f \n",
-//           ipr, pdgTrig, statusTrig, imother, ptTrig, particle->Eta(), particle->Phi()*TMath::RadToDeg());
+    if( GetDebug() > 2) printf("Select trigger particle %d: pdg %d status %d, mother index %d, pT %2.2f, eta %2.2f, phi %2.2f \n",
+                               ipr, pdgTrig, statusTrig, imother, ptTrig, particle->Eta(), particle->Phi()*TMath::RadToDeg());
     
 //    if(pdgTrig==111)
 //    {
 //      printf("\t pi0 daughters %d, %d\n", particle->GetDaughter(0), particle->GetDaughter(1));
 //    }
-    
 
     if     (pdgTrig==22 ) fhPtPhoton->Fill(ptTrig);
     else if(pdgTrig==111) fhPtPi0   ->Fill(ptTrig);
@@ -843,17 +1080,49 @@ void  AliAnaGeneratorKine::MakeAnalysisFillHistograms()
     Bool_t leading[4] ;
     Bool_t isolated[4] ;
 
-    particle->Momentum(trigger);
-    
-    IsLeadingAndIsolated(trigger, ipr, pdgTrig, leading, isolated);
+    IsLeadingAndIsolated(ipr, pdgTrig, leading, isolated);
     
     Int_t iparton = -1;
-    CorrelateWithPartonOrJet(trigger, ipr, pdgTrig, leading, isolated, iparton); 
+    Int_t ok = CorrelateWithPartonOrJet(ipr, pdgTrig, leading, isolated, iparton);
+    if(!ok) continue;
     
-    GetXE(trigger,ipr,pdgTrig,leading,isolated,iparton) ;    
+    GetXE(ipr,pdgTrig,leading,isolated,iparton) ;    
     
   }
   
   if(GetDebug() > 1) printf("AliAnaGeneratorKine::MakeAnalysisFillHistograms() - End fill histograms \n");
   
-} 
+}
+
+//_________________________________________________________
+void AliAnaGeneratorKine::SetTriggerDetector(TString & det)
+{
+  // Set the detrimeter for the analysis
+  
+  fTriggerDetectorString = det;
+  
+  if     (det=="EMCAL") fTriggerDetector = kEMCAL;
+  else if(det=="PHOS" ) fTriggerDetector = kPHOS;
+  else if(det=="CTS")   fTriggerDetector = kCTS;
+  else if(det=="DCAL")  fTriggerDetector = kDCAL;
+  else if(det.Contains("DCAL") && det.Contains("PHOS")) fTriggerDetector = kDCALPHOS;
+  else AliFatal(Form("Detector < %s > not known!", det.Data()));
+  
+}
+
+//_____________________________________________________
+void AliAnaGeneratorKine::SetTriggerDetector(Int_t det)
+{
+  // Set the detrimeter for the analysis
+  
+  fTriggerDetector = det;
+  
+  if     (det==kEMCAL)    fTriggerDetectorString = "EMCAL";
+  else if(det==kPHOS )    fTriggerDetectorString = "PHOS";
+  else if(det==kCTS)      fTriggerDetectorString = "CTS";
+  else if(det==kDCAL)     fTriggerDetectorString = "DCAL";
+  else if(det==kDCALPHOS) fTriggerDetectorString = "DCAL_PHOS";
+  else AliFatal(Form("Detector < %d > not known!", det));
+  
+}
+