setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetResponseMaker.cxx
index aafdd5e..23cc9ea 100644 (file)
 #include "AliJetResponseMaker.h"
 
 #include <TClonesArray.h>
-#include <TH1F.h>
 #include <TH2F.h>
-#include <TH3F.h>
+#include <THnSparse.h>
 #include <TLorentzVector.h>
 
+#include "AliAnalysisManager.h"
 #include "AliVCluster.h"
 #include "AliVTrack.h"
 #include "AliEmcalJet.h"
-#include "AliGenPythiaEventHeader.h"
 #include "AliLog.h"
-#include "AliMCEvent.h"
+#include "AliRhoParameter.h"
+#include "AliNamedArrayI.h"
+#include "AliJetContainer.h"
+#include "AliParticleContainer.h"
+#include "AliClusterContainer.h"
 
 ClassImp(AliJetResponseMaker)
 
 //________________________________________________________________________
 AliJetResponseMaker::AliJetResponseMaker() : 
   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
-  fMCTracksName("MCParticles"),
-  fMCJetsName("MCJets"),
-  fMaxDistance(0.25),
-  fDoWeighting(kFALSE),
-  fEventWeightHist(kFALSE),
-  fMCFiducial(kFALSE),
-  fMCminEta(0),
-  fMCmaxEta(0),
-  fMCminPhi(0),
-  fMCmaxPhi(0),
-  fSelectPtHardBin(-999),
-  fDoMatching(kTRUE),
-  fPythiaHeader(0),
-  fEventWeight(0),
-  fPtHardBin(0),
-  fNTrials(0),
-  fMCTracks(0),
-  fMCJets(0),
-  fHistNTrials(0),
-  fHistEvents(0),
-  fHistMCJetsPhiEta(0),
-  fHistMCJetsPtArea(0),
-  fHistMCJetsPhiEtaFiducial(0),
-  fHistMCJetsPtAreaFiducial(0),
-  fHistMCJetsNEFvsPt(0),
-  fHistMCJetsZvsPt(0),
-  fHistJetsPhiEta(0),
-  fHistJetsPtArea(0),
-  fHistJetsCorrPtArea(0),
-  fHistJetsNEFvsPt(0),
-  fHistJetsZvsPt(0),
-  fHistMatchingLevelMCPt(0),
-  fHistClosestDeltaEtaPhiMCPt(0),
-  fHistClosestDeltaPtMCPt(0),
-  fHistClosestDeltaCorrPtMCPt(0),
-  fHistNonMatchedMCJetsPtArea(0),
-  fHistNonMatchedJetsPtArea(0),
-  fHistNonMatchedJetsCorrPtArea(0),
-  fHistPartvsDetecPt(0),
-  fHistPartvsDetecCorrPt(0),
-  fHistMissedMCJetsPtArea(0)
+  fMatching(kNoMatching),
+  fMatchingPar1(0),
+  fMatchingPar2(0),
+  fUseCellsToMatch(kFALSE),
+  fMinJetMCPt(1),
+  fHistoType(0),
+  fDeltaPtAxis(0),
+  fDeltaEtaDeltaPhiAxis(0),
+  fNEFAxis(0),
+  fZAxis(0),
+  fIsJet1Rho(kFALSE),
+  fIsJet2Rho(kFALSE),
+  fHistJets1(0),
+  fHistJets2(0),
+  fHistMatching(0),
+  fHistJets1PhiEta(0),
+  fHistJets1PtArea(0),
+  fHistJets1CorrPtArea(0),
+  fHistJets1NEFvsPt(0),
+  fHistJets1CEFvsCEFPt(0),
+  fHistJets1ZvsPt(0),
+  fHistJets2PhiEta(0),
+  fHistJets2PtArea(0),
+  fHistJets2CorrPtArea(0),
+  fHistJets2NEFvsPt(0),
+  fHistJets2CEFvsCEFPt(0),
+  fHistJets2ZvsPt(0),
+  fHistCommonEnergy1vsJet1Pt(0),
+  fHistCommonEnergy2vsJet2Pt(0),
+  fHistDistancevsJet1Pt(0),
+  fHistDistancevsJet2Pt(0),
+  fHistDistancevsCommonEnergy1(0),
+  fHistDistancevsCommonEnergy2(0),
+  fHistCommonEnergy1vsCommonEnergy2(0),
+  fHistDeltaEtaDeltaPhi(0),
+  fHistJet2PtOverJet1PtvsJet2Pt(0),
+  fHistJet1PtOverJet2PtvsJet1Pt(0),
+  fHistDeltaPtvsJet1Pt(0),
+  fHistDeltaPtvsJet2Pt(0),
+  fHistDeltaPtOverJet1PtvsJet1Pt(0),
+  fHistDeltaPtOverJet2PtvsJet2Pt(0),
+  fHistDeltaPtvsDistance(0),
+  fHistDeltaPtvsCommonEnergy1(0),
+  fHistDeltaPtvsCommonEnergy2(0),
+  fHistDeltaPtvsArea1(0),
+  fHistDeltaPtvsArea2(0),
+  fHistDeltaPtvsDeltaArea(0),
+  fHistJet1PtvsJet2Pt(0),
+  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
+  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
+  fHistDeltaCorrPtvsJet1CorrPt(0),
+  fHistDeltaCorrPtvsJet2CorrPt(0),
+  fHistDeltaCorrPtvsDistance(0),
+  fHistDeltaCorrPtvsCommonEnergy1(0),
+  fHistDeltaCorrPtvsCommonEnergy2(0),
+  fHistDeltaCorrPtvsArea1(0),
+  fHistDeltaCorrPtvsArea2(0),
+  fHistDeltaCorrPtvsDeltaArea(0),
+  fHistJet1CorrPtvsJet2CorrPt(0),
+  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
+  fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
+  fHistDeltaMCPtvsJet1MCPt(0),
+  fHistDeltaMCPtvsJet2Pt(0),
+  fHistDeltaMCPtvsDistance(0),
+  fHistDeltaMCPtvsCommonEnergy1(0),
+  fHistDeltaMCPtvsCommonEnergy2(0),
+  fHistDeltaMCPtvsArea1(0),
+  fHistDeltaMCPtvsArea2(0),
+  fHistDeltaMCPtvsDeltaArea(0),
+  fHistJet1MCPtvsJet2Pt(0)
 {
   // Default constructor.
 
-  for (Int_t i = 0; i < 11; i++) {
-    fHistEventWeight[i] = 0;
-  }
-
   SetMakeGeneralHistograms(kTRUE);
 }
 
 //________________________________________________________________________
 AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
   AliAnalysisTaskEmcalJet(name, kTRUE),
-  fMCTracksName("MCParticles"),
-  fMCJetsName("MCJets"),
-  fMaxDistance(0.25),
-  fDoWeighting(kFALSE),
-  fEventWeightHist(kFALSE),
-  fMCFiducial(kFALSE),
-  fMCminEta(0),
-  fMCmaxEta(0),
-  fMCminPhi(0),
-  fMCmaxPhi(0),
-  fSelectPtHardBin(-999),
-  fDoMatching(kTRUE),
-  fPythiaHeader(0),
-  fEventWeight(0),
-  fPtHardBin(0),
-  fNTrials(0),
-  fMCTracks(0),
-  fMCJets(0),
-  fHistNTrials(0),
-  fHistEvents(0),
-  fHistMCJetsPhiEta(0),
-  fHistMCJetsPtArea(0),
-  fHistMCJetsPhiEtaFiducial(0),
-  fHistMCJetsPtAreaFiducial(0),
-  fHistMCJetsNEFvsPt(0),
-  fHistMCJetsZvsPt(0),
-  fHistJetsPhiEta(0),
-  fHistJetsPtArea(0),
-  fHistJetsCorrPtArea(0),
-  fHistJetsNEFvsPt(0),
-  fHistJetsZvsPt(0),
-  fHistMatchingLevelMCPt(0),
-  fHistClosestDeltaEtaPhiMCPt(0),
-  fHistClosestDeltaPtMCPt(0),
-  fHistClosestDeltaCorrPtMCPt(0),
-  fHistNonMatchedMCJetsPtArea(0),
-  fHistNonMatchedJetsPtArea(0),
-  fHistNonMatchedJetsCorrPtArea(0),
-  fHistPartvsDetecPt(0),
-  fHistPartvsDetecCorrPt(0),
-  fHistMissedMCJetsPtArea(0)
+  fMatching(kNoMatching),
+  fMatchingPar1(0),
+  fMatchingPar2(0),
+  fUseCellsToMatch(kFALSE),
+  fMinJetMCPt(1),
+  fHistoType(0),
+  fDeltaPtAxis(0),
+  fDeltaEtaDeltaPhiAxis(0),
+  fNEFAxis(0),
+  fZAxis(0),
+  fIsJet1Rho(kFALSE),
+  fIsJet2Rho(kFALSE),
+  fHistJets1(0),
+  fHistJets2(0),
+  fHistMatching(0),
+  fHistJets1PhiEta(0),
+  fHistJets1PtArea(0),
+  fHistJets1CorrPtArea(0),
+  fHistJets1NEFvsPt(0),
+  fHistJets1CEFvsCEFPt(0),
+  fHistJets1ZvsPt(0),
+  fHistJets2PhiEta(0),
+  fHistJets2PtArea(0),
+  fHistJets2CorrPtArea(0),
+  fHistJets2NEFvsPt(0),
+  fHistJets2CEFvsCEFPt(0),
+  fHistJets2ZvsPt(0),
+  fHistCommonEnergy1vsJet1Pt(0),
+  fHistCommonEnergy2vsJet2Pt(0),
+  fHistDistancevsJet1Pt(0),
+  fHistDistancevsJet2Pt(0),
+  fHistDistancevsCommonEnergy1(0),
+  fHistDistancevsCommonEnergy2(0),
+  fHistCommonEnergy1vsCommonEnergy2(0),
+  fHistDeltaEtaDeltaPhi(0),
+  fHistJet2PtOverJet1PtvsJet2Pt(0),
+  fHistJet1PtOverJet2PtvsJet1Pt(0),
+  fHistDeltaPtvsJet1Pt(0),
+  fHistDeltaPtvsJet2Pt(0),
+  fHistDeltaPtOverJet1PtvsJet1Pt(0),
+  fHistDeltaPtOverJet2PtvsJet2Pt(0),
+  fHistDeltaPtvsDistance(0),
+  fHistDeltaPtvsCommonEnergy1(0),
+  fHistDeltaPtvsCommonEnergy2(0),
+  fHistDeltaPtvsArea1(0),
+  fHistDeltaPtvsArea2(0),
+  fHistDeltaPtvsDeltaArea(0),
+  fHistJet1PtvsJet2Pt(0),
+  fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
+  fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
+  fHistDeltaCorrPtvsJet1CorrPt(0),
+  fHistDeltaCorrPtvsJet2CorrPt(0),
+  fHistDeltaCorrPtvsDistance(0),
+  fHistDeltaCorrPtvsCommonEnergy1(0),
+  fHistDeltaCorrPtvsCommonEnergy2(0),
+  fHistDeltaCorrPtvsArea1(0),
+  fHistDeltaCorrPtvsArea2(0),
+  fHistDeltaCorrPtvsDeltaArea(0),
+  fHistJet1CorrPtvsJet2CorrPt(0),
+  fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
+  fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
+  fHistDeltaMCPtvsJet1MCPt(0),
+  fHistDeltaMCPtvsJet2Pt(0),
+  fHistDeltaMCPtvsDistance(0),
+  fHistDeltaMCPtvsCommonEnergy1(0),
+  fHistDeltaMCPtvsCommonEnergy2(0),
+  fHistDeltaMCPtvsArea1(0),
+  fHistDeltaMCPtvsArea2(0),
+  fHistDeltaMCPtvsDeltaArea(0),
+  fHistJet1MCPtvsJet2Pt(0)
 {
   // Standard constructor.
 
-  for (Int_t i = 0; i < 11; i++) {
-    fHistEventWeight[i] = 0;
-  }
-
   SetMakeGeneralHistograms(kTRUE);
 }
 
@@ -135,263 +188,901 @@ AliJetResponseMaker::~AliJetResponseMaker()
   // Destructor
 }
 
+
 //________________________________________________________________________
-void AliJetResponseMaker::UserCreateOutputObjects()
+void AliJetResponseMaker::AllocateTH2()
 {
-  // Create user objects.
+  // Allocate TH2 histograms.
 
-  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+  fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
+  fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
+  fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistJets1PhiEta);
+  
+  fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets1PtArea->GetXaxis()->SetTitle("area");
+  fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
+  fHistJets1PtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets1PtArea);
+  if (fIsJet1Rho) {
+    fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 
+                                   fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
+    fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
+    fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistJets1CorrPtArea);
+  }
 
-  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
-  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
+  fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
+  fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
+  fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets1ZvsPt);
+  
+  fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
+  fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
+  fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets1NEFvsPt);
+  
+  fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
+  fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
+  fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets1CEFvsCEFPt);
+
+  // Jets 2 histograms
+
+  fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
+  fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
+  fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
+  fOutput->Add(fHistJets2PhiEta);
+
+  fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets2PtArea->GetXaxis()->SetTitle("area");
+  fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
+  fHistJets2PtArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets2PtArea);
+
+  if (fIsJet2Rho) {
+    fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 
+                                   fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
+    fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
+    fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistJets2CorrPtArea);
+  }
 
-  fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
-  fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
-  fHistNTrials->GetYaxis()->SetTitle("trials");
-  fOutput->Add(fHistNTrials);
+  fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
+  fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
+  fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets2ZvsPt);
+  
+  fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
+  fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
+  fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets2NEFvsPt);
+  
+  fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
+  fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
+  fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJets2CEFvsCEFPt);
+
+  // Matching histograms
+
+  fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
+  fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
+  fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");  
+  fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistCommonEnergy1vsJet1Pt);
+
+  fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
+  fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
+  fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");  
+  fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistCommonEnergy2vsJet2Pt);
+
+  fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
+  fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
+  fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");  
+  fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDistancevsJet1Pt);
+
+  fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
+  fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
+  fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");  
+  fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDistancevsJet2Pt);
+
+  fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
+  fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
+  fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");  
+  fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDistancevsCommonEnergy1);
+
+  fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
+  fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
+  fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");  
+  fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDistancevsCommonEnergy2);
+
+  fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
+  fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
+  fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");  
+  fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
+
+  fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
+  fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
+  fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");  
+  fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaEtaDeltaPhi);
+
+  fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
+  fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
+  fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
+  fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
+
+  fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
+  fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
+  fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
+  fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
+
+  fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", 
+                                 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
+  fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsJet1Pt);
+
+  fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", 
+                                 fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
+  fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsJet2Pt);
+
+  fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt", 
+                                           fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
+  fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
+  fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}");
+  fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt);
+
+  fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt", 
+                                           fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
+  fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
+  fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
+  fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt);
+
+  fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance", 
+                                   fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");  
+  fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsDistance);
+
+  fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1",
+                                        fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");  
+  fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsCommonEnergy1);
+
+  fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2", 
+                                        fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");  
+  fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsCommonEnergy2);
+
+  fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1", 
+                                fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
+  fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsArea1);
+
+  fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2", 
+                                fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
+  fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsArea2);
+
+  fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea", 
+                                    fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+  fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
+  fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+  fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistDeltaPtvsDeltaArea);
+
+  fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+  fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
+  fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
+  fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
+  fOutput->Add(fHistJet1PtvsJet2Pt);
+
+  if (fIsJet1Rho || fIsJet2Rho) {
+    if (!fIsJet1Rho) 
+      fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", 
+                                             fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    else
+      fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", 
+                                             2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+      
+    fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");  
+    fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt);
+
+    if (!fIsJet2Rho) 
+      fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", 
+                                             fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);                                         
+    else
+      fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", 
+                                             2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+
+    fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");  
+    fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt);
+
+    fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", 
+                                                         2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
+    fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");  
+    fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}");
+    fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt);
+
+    fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", 
+                                                         2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
+    fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");  
+    fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}");
+    fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt);
+
+    fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance", 
+                                         fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");  
+    fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsDistance);
+
+    fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1", 
+                                              fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");  
+    fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1);
+
+    fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2", 
+                                              fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");  
+    fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
+
+    fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1", 
+                                      fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
+    fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsArea1);
+    
+    fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2", 
+                                      fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
+    fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsArea2);
+    
+    fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea", 
+                                          fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
+    fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
+    fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
+    
+    if (!fIsJet1Rho) 
+      fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 
+                                            fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    else if (!fIsJet2Rho) 
+      fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 
+                                            2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+    else
+      fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 
+                                            2*fNbins, -fMaxBinPt, fMaxBinPt, 
+                                            2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
+    fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
+    fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
+  }
 
-  fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
-  fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
-  fHistEvents->GetYaxis()->SetTitle("total events");
-  fOutput->Add(fHistEvents);
+  if (fIsEmbedded) {
+    fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt", 
+                                       fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");  
+    fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsJet1MCPt);
+
+    fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt", 
+                                     fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
+    fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsJet2Pt);
+
+    fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", 
+                                                   fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
+    fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");  
+    fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}");
+    fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt);
+
+    fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt", 
+                                               fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
+    fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
+    fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
+    fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt);
+
+    fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance", 
+                                       fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");  
+    fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsDistance);
+
+    fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1", 
+                                            fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");  
+    fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsCommonEnergy1);
+
+    fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2", 
+                                            fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");  
+    fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
+
+    fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1", 
+                                    fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
+    fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsArea1);
+
+    fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2", 
+                                    fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
+    fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsArea2);
+
+    fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea", 
+                                        fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
+    fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
+    fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
+    fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistDeltaMCPtvsDeltaArea);
+
+    fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt", 
+                                    fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+    fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
+    fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
+    fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
+    fOutput->Add(fHistJet1MCPtvsJet2Pt);
+  }
+}
 
-  if (fEventWeightHist) {
-    for (Int_t i = 0; i < 11; i++) {
-      TString name(Form("fHistEventWeight_%d", i+1));
-      fHistEventWeight[i] = new TH1F(name, name, 10, 0, 10);
-      fOutput->Add(fHistEventWeight[i]);
-    }
+//________________________________________________________________________
+void AliJetResponseMaker::AllocateTHnSparse()
+{
+  // Allocate THnSparse histograms.
+
+  TString title[25]= {""};
+  Int_t nbins[25]  = {0};
+  Double_t min[25] = {0.};
+  Double_t max[25] = {0.};
+  Int_t dim = 0;
+
+  title[dim] = "#phi";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
+  dim++;
+
+  title[dim] = "#eta";
+  nbins[dim] = fNbins/4;
+  min[dim] = -1;
+  max[dim] = 1;
+  dim++;
+
+  title[dim] = "p_{T}";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "A_{jet}";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1.5;
+  dim++;
+
+  if (fNEFAxis) {
+    title[dim] = "NEF";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
   }
 
-  for (Int_t i = 1; i < 12; i++) {
-    fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
-    fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
+  if (fZAxis) {
+    title[dim] = "Z";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
   }
 
-  fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 20, -2, 2, 32, 0, 6.4);
-  fHistJetsPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistJetsPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistJetsPhiEta);
-  
-  fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
-  fHistJetsPtArea->GetXaxis()->SetTitle("area");
-  fHistJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
-  fHistJetsPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistJetsPtArea);
-
-  fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt);
-  fHistJetsCorrPtArea->GetXaxis()->SetTitle("area");
-  fHistJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
-  fHistJetsCorrPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistJetsCorrPtArea);
-  
-  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
-  fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
-  fOutput->Add(fHistJetsZvsPt);
-  
-  fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
-  fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
-  fOutput->Add(fHistJetsNEFvsPt);
-
-  fHistMCJetsPhiEta = new TH2F("fHistMCJetsPhiEta", "fHistMCJetsPhiEta", 20, -2, 2, 32, 0, 6.4);
-  fHistMCJetsPhiEta->GetXaxis()->SetTitle("#eta");
-  fHistMCJetsPhiEta->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistMCJetsPhiEta);
-  
-  fHistMCJetsPtArea = new TH2F("fHistMCJetsPtArea", "fHistMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMCJetsPtArea->GetXaxis()->SetTitle("area");
-  fHistMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
-  fHistMCJetsPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistMCJetsPtArea);
-
-  fHistMCJetsPhiEtaFiducial = new TH2F("fHistMCJetsPhiEtaFiducial", "fHistMCJetsPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4);
-  fHistMCJetsPhiEtaFiducial->GetXaxis()->SetTitle("#eta");
-  fHistMCJetsPhiEtaFiducial->GetYaxis()->SetTitle("#phi");
-  fOutput->Add(fHistMCJetsPhiEtaFiducial);
-  
-  fHistMCJetsPtAreaFiducial = new TH2F("fHistMCJetsPtAreaFiducial", "fHistMCJetsPtAreaFiducial", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMCJetsPtAreaFiducial->GetXaxis()->SetTitle("area");
-  fHistMCJetsPtAreaFiducial->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
-  fHistMCJetsPtAreaFiducial->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistMCJetsPtAreaFiducial);
+  title[dim] = "p_{T,particle}^{leading} (GeV/c)";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 120;
+  dim++;
+
+  Int_t dim1 = dim, dim2 = dim;
+
+  if (fIsJet1Rho) {
+    title[dim1] = "p_{T}^{corr}";
+    nbins[dim1] = fNbins*2;
+    min[dim1] = -250;
+    max[dim1] = 250;
+    dim1++;
+  }
+
+  if (fIsEmbedded) {
+    title[dim1] = "p_{T}^{MC}";
+    nbins[dim1] = fNbins;
+    min[dim1] = 0;
+    max[dim1] = 250;
+    dim1++;
+  }
+
+  fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
+  for (Int_t i = 0; i < dim1; i++) 
+    fHistJets1->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistJets1);
+
+  if (fIsJet2Rho) {
+    title[dim2] = "p_{T}^{corr}";
+    nbins[dim2] = fNbins*2;
+    min[dim2] = -250;
+    max[dim2] = 250;
+    dim2++;
+  }
+
+  fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
+  for (Int_t i = 0; i < dim2; i++) 
+    fHistJets2->GetAxis(i)->SetTitle(title[i]);
+  fOutput->Add(fHistJets2);
   
-  fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
-  fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
-  fOutput->Add(fHistMCJetsZvsPt);
-
-  fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
-  fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
-  fOutput->Add(fHistMCJetsNEFvsPt);
-
-  fHistMatchingLevelMCPt = new TH2F("fHistMatchingLevelMCPt", "fHistMatchingLevelMCPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMatchingLevelMCPt->GetXaxis()->SetTitle("Matching level");
-  fHistMatchingLevelMCPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
-  fOutput->Add(fHistMatchingLevelMCPt);
-
-  fHistClosestDeltaEtaPhiMCPt = new TH3F("fHistClosestDeltaEtaPhiMCPt", "fHistClosestDeltaEtaPhiMCPt", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2, 64, -1.6, 4.8, fNbins, fMinBinPt, fMaxBinPt);
-  fHistClosestDeltaEtaPhiMCPt->GetXaxis()->SetTitle("#Delta#eta");
-  fHistClosestDeltaEtaPhiMCPt->GetYaxis()->SetTitle("#Delta#phi");
-  fHistClosestDeltaEtaPhiMCPt->GetZaxis()->SetTitle("p_{T}^{gen}");
-  fOutput->Add(fHistClosestDeltaEtaPhiMCPt);
-
-  fHistClosestDeltaPtMCPt = new TH2F("fHistClosestDeltaPtMCPt", "fHistClosestDeltaPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt);
-  fHistClosestDeltaPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}");  
-  fHistClosestDeltaPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{rec} (GeV/c)");
-  fHistClosestDeltaPtMCPt->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistClosestDeltaPtMCPt);
-
-  fHistClosestDeltaCorrPtMCPt = new TH2F("fHistClosestDeltaCorrPtMCPt", "fHistClosestDeltaCorrPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt);
-  fHistClosestDeltaCorrPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}");  
-  fHistClosestDeltaCorrPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
-  fHistClosestDeltaCorrPtMCPt->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistClosestDeltaCorrPtMCPt);
-
-  fHistNonMatchedMCJetsPtArea = new TH2F("fHistNonMatchedMCJetsPtArea", "fHistNonMatchedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
-  fHistNonMatchedMCJetsPtArea->GetXaxis()->SetTitle("area");
-  fHistNonMatchedMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
-  fHistNonMatchedMCJetsPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedMCJetsPtArea);
-
-  fHistNonMatchedJetsPtArea = new TH2F("fHistNonMatchedJetsPtArea", "fHistNonMatchedJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
-  fHistNonMatchedJetsPtArea->GetXaxis()->SetTitle("area");
-  fHistNonMatchedJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)");
-  fHistNonMatchedJetsPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedJetsPtArea);
-
-  fHistNonMatchedJetsCorrPtArea = new TH2F("fHistNonMatchedJetsCorrPtArea", "fHistNonMatchedJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt);
-  fHistNonMatchedJetsCorrPtArea->GetXaxis()->SetTitle("area");
-  fHistNonMatchedJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
-  fHistNonMatchedJetsCorrPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistNonMatchedJetsCorrPtArea);
-
-  fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
-  fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
-  fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
-  fOutput->Add(fHistPartvsDetecPt);
-
-  fHistPartvsDetecCorrPt = new TH2F("fHistPartvsDetecCorrPt", "fHistPartvsDetecCorrPt", (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
-  fHistPartvsDetecCorrPt->GetXaxis()->SetTitle("p_{T}^{corr}");
-  fHistPartvsDetecCorrPt->GetYaxis()->SetTitle("p_{T}^{gen}");
-  fOutput->Add(fHistPartvsDetecCorrPt);
-
-  fHistMissedMCJetsPtArea = new TH2F("fHistMissedMCJetsPtArea", "fHistMissedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
-  fHistMissedMCJetsPtArea->GetXaxis()->SetTitle("area");  
-  fHistMissedMCJetsPtArea->GetYaxis()->SetTitle("p_{T} (GeV/c)");
-  fHistMissedMCJetsPtArea->GetZaxis()->SetTitle("counts");
-  fOutput->Add(fHistMissedMCJetsPtArea);
+  // Matching
+
+  dim = 0;
+
+  title[dim] = "p_{T,1}";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "p_{T,2}";
+  nbins[dim] = fNbins;
+  min[dim] = 0;
+  max[dim] = 250;
+  dim++;
+
+  title[dim] = "A_{jet,1}";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1.5;
+  dim++;
+
+  title[dim] = "A_{jet,2}";
+  nbins[dim] = fNbins/4;
+  min[dim] = 0;
+  max[dim] = 1.5;
+  dim++;
+
+  title[dim] = "distance";
+  nbins[dim] = fNbins/2;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "CE1";
+  nbins[dim] = fNbins/2;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "CE2";
+  nbins[dim] = fNbins/2;
+  min[dim] = 0;
+  max[dim] = 1.2;
+  dim++;
+
+  title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 120;
+  dim++;
+
+  title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
+  nbins[dim] = 120;
+  min[dim] = 0;
+  max[dim] = 120;
+  dim++;
+
+  if (fDeltaPtAxis) {
+    title[dim] = "#deltaA_{jet}";
+    nbins[dim] = fNbins/2;
+    min[dim] = -1;
+    max[dim] = 1;
+    dim++;
+
+    title[dim] = "#deltap_{T}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (fIsJet1Rho) {
+    title[dim] = "p_{T,1}^{corr}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (fIsJet2Rho) {
+    title[dim] = "p_{T,2}^{corr}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
+    title[dim] = "#deltap_{T}^{corr}";
+    nbins[dim] = fNbins*2;
+    min[dim] = -250;
+    max[dim] = 250;
+    dim++;
+  }
+  if (fDeltaEtaDeltaPhiAxis) {
+    title[dim] = "#delta#eta";
+    nbins[dim] = fNbins/2;
+    min[dim] = -1;
+    max[dim] = 1;
+    dim++;
+
+    title[dim] = "#delta#phi";
+    nbins[dim] = fNbins/2;
+    min[dim] = -TMath::Pi()/2;
+    max[dim] = TMath::Pi()*3/2;
+    dim++;
+  }
+  if (fIsEmbedded) {
+    title[dim] = "p_{T,1}^{MC}";
+    nbins[dim] = fNbins;
+    min[dim] = 0;
+    max[dim] = 250;
+    dim++;
+
+    if (fDeltaPtAxis) {
+      title[dim] = "#deltap_{T}^{MC}";
+      nbins[dim] = fNbins*2;
+      min[dim] = -250;
+      max[dim] = 250;
+      dim++;
+    }
+  }
 
-  PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
+  if (fNEFAxis) {
+    title[dim] = "NEF_{1}";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
+
+    title[dim] = "NEF_{2}";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
+  }
+
+  if (fZAxis) {
+    title[dim] = "Z_{1}";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
+
+    title[dim] = "Z_{2}";
+    nbins[dim] = fNbins/4;
+    min[dim] = 0;
+    max[dim] = 1.2;
+    dim++;
+  }
+      
+  fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
+    
+  for (Int_t i = 0; i < dim; i++)
+    fHistMatching->GetAxis(i)->SetTitle(title[i]);
+   
+  fOutput->Add(fHistMatching);
 }
 
 //________________________________________________________________________
-Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
-{   
-  // Return true if jet is accepted.
+void AliJetResponseMaker::UserCreateOutputObjects()
+{
+  // Create user objects.
 
-  if (jet->Pt() <= fJetPtCut)
-    return kFALSE;
-  if (jet->Area() <= fJetAreaCut)
-    return kFALSE;
+  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 
-  return kTRUE;
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
+
+  if (!jets1 || !jets2) return;
+
+  if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
+  else fIsJet1Rho = kTRUE;
+  if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
+  else fIsJet2Rho = kTRUE;
+
+  if (fHistoType==0)
+    AllocateTH2();
+  else 
+    AllocateTHnSparse();
+
+  PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 }
 
 //________________________________________________________________________
-void AliJetResponseMaker::ExecOnce()
+void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t LeadingPt, Double_t CorrPt, Double_t MCPt, Int_t Set)
 {
-  // Execute once.
-
-  if (!fMCJetsName.IsNull() && !fMCJets) {
-    fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
-    if (!fMCJets) {
-      AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
-      return;
-    }
-    else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
-      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data())); 
-      fMCJets = 0;
+  if (fHistoType==1) {
+    THnSparse *histo = 0;
+    if (Set==1)
+      histo = fHistJets1;
+    else if (Set==2)
+      histo = fHistJets2;
+
+    if (!histo)
       return;
-    }
-  }
 
-  if (!fMCTracksName.IsNull() && !fMCTracks) {
-    fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
-    if (!fMCTracks) {
-      AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data())); 
-      return;
-    }
-    else {
-      TClass *cl = fMCTracks->GetClass();
-      if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
-       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data())); 
-       fMCTracks = 0;
-       return;
-      }
+    Double_t contents[20]={0};
+
+    for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
+      TString title(histo->GetAxis(i)->GetTitle());
+      if (title=="#phi")
+       contents[i] = Phi;
+      else if (title=="#eta")
+       contents[i] = Eta;
+      else if (title=="p_{T}")
+       contents[i] = Pt;
+      else if (title=="A_{jet}")
+       contents[i] = A;
+      else if (title=="NEF")
+       contents[i] = NEF;
+      else if (title=="Z")
+       contents[i] = LeadingPt/Pt;
+      else if (title=="p_{T}^{corr}")
+       contents[i] = CorrPt;
+      else if (title=="p_{T}^{MC}")
+       contents[i] = MCPt;
+      else if (title=="p_{T,particle}^{leading} (GeV/c)")
+       contents[i] = LeadingPt;
+      else 
+       AliWarning(Form("Unable to fill dimension %s!",title.Data()));
     }
-  }
-
-  AliAnalysisTaskEmcalJet::ExecOnce();
 
-  if (fMCFiducial) {
-    fMCminEta = fJetMinEta;
-    fMCmaxEta = fJetMaxEta;
-    fMCminPhi = fJetMinPhi;
-    fMCmaxPhi = fJetMaxPhi;
+    histo->Fill(contents);
   }
   else {
-    fMCminEta = fJetMinEta - fJetRadius;
-    fMCmaxEta = fJetMaxEta + fJetRadius;
-    fMCminPhi = fJetMinPhi - fJetRadius;
-    fMCmaxPhi = fJetMaxPhi + fJetRadius;
+    if (Set == 1) {
+      fHistJets1PtArea->Fill(A, Pt);
+      fHistJets1PhiEta->Fill(Eta, Phi);
+      fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt);
+      fHistJets1NEFvsPt->Fill(NEF, Pt);
+      fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
+      if (fIsJet1Rho)
+       fHistJets1CorrPtArea->Fill(A, CorrPt);
+    }
+    else if (Set == 2) {
+      fHistJets2PtArea->Fill(A, Pt);
+      fHistJets2PhiEta->Fill(Eta, Phi);
+      fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt);
+      fHistJets2NEFvsPt->Fill(NEF, Pt);
+      fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
+      if (fIsJet2Rho)
+       fHistJets2CorrPtArea->Fill(A, CorrPt);
+    }
   }
 }
 
 //________________________________________________________________________
-Bool_t AliJetResponseMaker::IsEventSelected()
+void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, 
+                                            Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, 
+                                            Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
 {
-  // Check if event is selected
+  if (fHistoType==1) {
+    Double_t contents[20]={0};
+
+    for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
+      TString title(fHistMatching->GetAxis(i)->GetTitle());
+      if (title=="p_{T,1}")
+       contents[i] = Pt1;
+      else if (title=="p_{T,2}")
+       contents[i] = Pt2;
+      else if (title=="A_{jet,1}")
+       contents[i] = A1;
+      else if (title=="A_{jet,2}")
+       contents[i] = A2;
+      else if (title=="distance")
+       contents[i] = d;
+      else if (title=="CE1")
+       contents[i] = CE1;
+      else if (title=="CE2")
+       contents[i] = CE2;
+      else if (title=="#deltaA_{jet}")
+       contents[i] = A1-A2;
+      else if (title=="#deltap_{T}")
+       contents[i] = Pt1-Pt2;
+      else if (title=="#delta#eta")
+       contents[i] = Eta1-Eta2;
+      else if (title=="#delta#phi")
+       contents[i] = Phi1-Phi2;
+      else if (title=="p_{T,1}^{corr}")
+       contents[i] = CorrPt1;
+      else if (title=="p_{T,2}^{corr}")
+       contents[i] = CorrPt2;
+      else if (title=="#deltap_{T}^{corr}")
+       contents[i] = CorrPt1-CorrPt2;
+      else if (title=="p_{T,1}^{MC}")
+       contents[i] = MCPt1;
+      else if (title=="#deltap_{T}^{MC}")
+       contents[i] = MCPt1-Pt2;
+      else if (title=="NEF_{1}")
+       contents[i] = NEF1;
+      else if (title=="NEF_{2}")
+       contents[i] = NEF2;
+      else if (title=="Z_{1}")
+       contents[i] = LeadingPt1/Pt1;
+      else if (title=="Z_{2}")
+       contents[i] = LeadingPt2/Pt2;
+      else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
+       contents[i] = LeadingPt1;
+      else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
+       contents[i] = LeadingPt2;
+      else 
+       AliWarning(Form("Unable to fill dimension %s!",title.Data()));
+    }
 
-  if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) 
-    return kFALSE;
+    fHistMatching->Fill(contents);
+  }
+  else {
+    fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
+    fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
+
+    fHistDistancevsJet1Pt->Fill(d, Pt1);
+    fHistDistancevsJet2Pt->Fill(d, Pt2);
+
+    fHistDistancevsCommonEnergy1->Fill(d, CE1);
+    fHistDistancevsCommonEnergy2->Fill(d, CE2);
+    fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
+
+    fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
+
+    fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
+    fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
+
+    Double_t dpt = Pt1 - Pt2;
+    fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
+    fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
+    fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
+    fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
+
+    fHistDeltaPtvsDistance->Fill(d, dpt);
+    fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
+    fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
+
+    fHistDeltaPtvsArea1->Fill(A1, dpt);
+    fHistDeltaPtvsArea2->Fill(A2, dpt);
+
+    Double_t darea = A1 - A2;
+    fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
+
+    fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
+
+    if (fIsJet1Rho || fIsJet2Rho) {
+      Double_t dcorrpt = CorrPt1 - CorrPt2;
+      fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
+      fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
+      fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
+      fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
+      fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
+      fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
+      fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
+      fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
+      fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
+      fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
+      fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
+    }
 
-  return AliAnalysisTaskEmcalJet::IsEventSelected();
+    if (fIsEmbedded) {
+      Double_t dmcpt = MCPt1 - Pt2;
+      fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
+      fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
+      fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
+      fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
+      fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
+      fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
+      fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
+      fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
+      fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
+      fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
+      fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
+    }
+  }
 }
 
 //________________________________________________________________________
-Bool_t AliJetResponseMaker::RetrieveEventObjects()
+void AliJetResponseMaker::ExecOnce()
 {
-  // Retrieve event objects.
-
-  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
-    return kFALSE;
-  
-  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
-
-  if (!fPythiaHeader)
-    return kFALSE;
+  // Execute once.
 
-  if (fDoWeighting)
-    fEventWeight = fPythiaHeader->EventWeight();
-  else
-    fEventWeight = 1;
+  AliAnalysisTaskEmcalJet::ExecOnce();
 
-  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
-  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
 
-  Double_t pthard = fPythiaHeader->GetPtHard();
+  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
 
-  for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
-    if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
-      break;
+  if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
+    if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
+      AliWarning("Changing matching type from MC label to same collection...");
+      fMatching = kSameCollections;
+    }
+    else {
+      AliWarning("Changing matching type from MC label to geometrical...");
+      fMatching = kGeometrical;
+    }
+  }
+  else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
+    if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
+      AliWarning("Changing matching type from same collection to MC label...");
+      fMatching = kMCLabel;
+    }
+    else {
+      AliWarning("Changing matching type from same collection to geometrical...");
+      fMatching = kGeometrical;
+    }
   }
-
-  fNTrials = fPythiaHeader->Trials();
-
-  return kTRUE;
 }
 
 //________________________________________________________________________
@@ -399,231 +1090,580 @@ Bool_t AliJetResponseMaker::Run()
 {
   // Find the closest jets
 
-  if (!fDoMatching)
+  if (fMatching == kNoMatching) 
     return kTRUE;
+  else
+    return DoJetMatching();
+}
 
-  DoJetLoop(fJets, fMCJets, kFALSE);
-  DoJetLoop(fMCJets, fJets, kTRUE);
-
-  const Int_t nMCJets = fMCJets->GetEntriesFast();
+//________________________________________________________________________
+Bool_t AliJetResponseMaker::DoJetMatching()
+{
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
 
-  for (Int_t i = 0; i < nMCJets; i++) {
+  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
 
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
+  DoJetLoop();
 
-    if (!jet) {
-      AliError(Form("Could not receive jet %d", i));
-      continue;
-    }  
+  AliEmcalJet* jet1 = 0;
 
-    if (!AcceptJet(jet))
-      continue;
+  jets1->ResetCurrentID();
+  while ((jet1 = jets1->GetNextJet())) {
 
-    if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi)
-      continue;
+    AliEmcalJet *jet2 = jet1->ClosestJet();
 
-    if (jet->Pt() > fMaxBinPt)
-      continue;
+    if (!jet2) continue;
+    if (jet2->ClosestJet() != jet1) continue;
+    if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
 
-    if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && 
-        jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
-      jet->SetMatchedToClosest();
-      jet->ClosestJet()->SetMatchedToClosest();
-    }
+    // Matched jet found
+    jet1->SetMatchedToClosest(fMatching);
+    jet2->SetMatchedToClosest(fMatching);
+    AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f", 
+                   jet1->Pt(), jet1->Eta(), jet1->Phi(), 
+                   jet2->Pt(), jet2->Eta(), jet2->Phi()));
   }
 
   return kTRUE;
 }
 
 //________________________________________________________________________
-void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc)
+void AliJetResponseMaker::DoJetLoop()
 {
   // Do the jet loop.
 
-  Int_t nJets1 = jets1->GetEntriesFast();
-  Int_t nJets2 = jets2->GetEntriesFast();
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
+
+  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
 
-  for (Int_t i = 0; i < nJets1; i++) {
+  AliEmcalJet* jet1 = 0;
+  AliEmcalJet* jet2 = 0;
+
+  jets2->ResetCurrentID();
+  while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
+    
+  jets1->ResetCurrentID();
+  while ((jet1 = jets1->GetNextJet())) {
+    jet1->ResetMatching();
+    
+    if (jet1->MCPt() < fMinJetMCPt) continue;
+
+    jets2->ResetCurrentID();
+    while ((jet2 = jets2->GetNextJet())) {
+      SetMatchingLevel(jet1, jet2, fMatching);
+    } // jet2 loop
+  } // jet1 loop
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
+{
+  Double_t deta = jet2->Eta() - jet1->Eta();
+  Double_t dphi = jet2->Phi() - jet1->Phi();
+  d = TMath::Sqrt(deta * deta + dphi * dphi);
+}
 
-    AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
+//________________________________________________________________________
+void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
+{ 
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
+
+  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
+
+  AliParticleContainer *tracks1   = jets1->GetParticleContainer();
+  AliClusterContainer  *clusters1 = jets1->GetClusterContainer();
+  AliParticleContainer *tracks2   = jets2->GetParticleContainer();
+
+  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
+  d1 = jet1->Pt();
+  d2 = jet2->Pt();
+  Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
+
+  // remove completely tracks that are not MC particles (label == 0)
+  if (tracks1 && tracks1->GetArray()) {
+    for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
+      AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
+      if (!track) {
+       AliWarning(Form("Could not find track %d!", iTrack));
+       continue;
+      }
 
-    if (!jet1) {
-      AliError(Form("Could not receive jet %d", i));
-      continue;
-    }  
+      Int_t MClabel = TMath::Abs(track->GetLabel());
+      MClabel -= fMCLabelShift;
+      if (MClabel != 0) continue;
 
-    if (!AcceptJet(jet1))
-      continue;
+      // this is not a MC particle; remove it completely
+      AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
+      totalPt1 -= track->Pt();
+      d1 -= track->Pt();
+    }
+  }
 
-    if (!mc) {
-      if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
+  // remove completely clusters that are not MC particles (label == 0)
+  if (fUseCellsToMatch && fCaloCells) { 
+    for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+      AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
+      if (!clus) {
+       AliWarning(Form("Could not find cluster %d!", iClus));
        continue;
+      }
+      TLorentzVector part;
+      clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+         
+      for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
+       Int_t cellId = clus->GetCellAbsId(iCell);
+       Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
+
+       Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
+       MClabel -= fMCLabelShift;
+       if (MClabel != 0) continue;
+
+       // this is not a MC particle; remove it completely
+       AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
+       totalPt1 -= part.Pt() * cellFrac;
+       d1 -= part.Pt() * cellFrac;
+      }
     }
-    else {
-      if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi)
+  }
+  else {
+    for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+      AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
+      if (!clus) {
+       AliWarning(Form("Could not find cluster %d!", iClus));
        continue;
+      }
+      TLorentzVector part;
+      clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+         
+      Int_t MClabel = TMath::Abs(clus->GetLabel());
+      MClabel -= fMCLabelShift;
+      if (MClabel != 0) continue;
+
+      // this is not a MC particle; remove it completely
+      AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
+      totalPt1 -= part.Pt();
+      d1 -= part.Pt();
     }
+  }
 
-    for (Int_t j = 0; j < nJets2; j++) {
-      
-      AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
-      
-      if (!jet2) {
-       AliError(Form("Could not receive jet %d", j));
+  for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
+    Bool_t track2Found = kFALSE;
+    Int_t index2 = jet2->TrackAt(iTrack2);
+
+    // now look for common particles in the track array
+    for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
+      AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
+      if (!track) {
+       AliWarning(Form("Could not find track %d!", iTrack));
        continue;
-      }  
-      
-      if (!AcceptJet(jet2))
+      }
+      Int_t MClabel = TMath::Abs(track->GetLabel());
+      MClabel -= fMCLabelShift;          
+      if (MClabel <= 0) continue;
+
+      Int_t index = -1;
+      index = tracks2->GetIndexFromLabel(MClabel);
+      if (index < 0) {
+       AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
        continue;
+      }
+       
+      if (index2 != index) continue;
+      
+      // found common particle
+      d1 -= track->Pt();
+
+      if (!track2Found) {
+       AliVParticle *MCpart = tracks2->GetParticle(index2);
+       AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+                       iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
+       d2 -= MCpart->Pt();
+      }
 
-      if (mc) {
-       if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
+      track2Found = kTRUE;
+    }
+
+    // now look for common particles in the cluster array
+    if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
+      for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+       AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
+       if (!clus) {
+         AliWarning(Form("Could not find cluster %d!", iClus));
          continue;
+       }
+       TLorentzVector part;
+       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
+      
+       for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
+         Int_t cellId = clus->GetCellAbsId(iCell);
+         Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
+       
+         Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
+         MClabel -= fMCLabelShift;
+         if (MClabel <= 0) continue;
+         
+         Int_t index1 = -1;
+         index1 = tracks2->GetIndexFromLabel(MClabel);
+         if (index1 < 0) {
+           AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
+           continue;
+         }
+       
+         if (index2 != index1) continue;
+         
+         // found common particle
+         d1 -= part.Pt() * cellFrac;
+       
+         if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
+           AliVParticle *MCpart = tracks2->GetParticle(index2);
+           AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+                           iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));               
+           d2 -= MCpart->Pt() * cellFrac;
+         }
+
+         track2Found = kTRUE;
+       }
       }
-      else {
-       if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi)
+    }
+    else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
+      for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
+       AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
+       if (!clus) {
+         AliWarning(Form("Could not find cluster %d!", iClus));
          continue;
-      }
+       }
+       TLorentzVector part;
+       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
       
-      Double_t deta = jet2->Eta() - jet1->Eta();
-      Double_t dphi = jet2->Phi() - jet1->Phi();
-      Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
+       Int_t MClabel = TMath::Abs(clus->GetLabel());
+       MClabel -= fMCLabelShift;           
+       if (MClabel <= 0) continue;
 
-      if (d < jet1->ClosestJetDistance()) {
-       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
-       jet1->SetClosestJet(jet2, d);
-      }
-      else if (d < jet1->SecondClosestJetDistance()) {
-       jet1->SetSecondClosestJet(jet2, d);
+       Int_t index = -1;
+       index = tracks2->GetIndexFromLabel(MClabel);
+      
+       if (index < 0) {
+         AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
+         continue;
+       }
+
+       if (index2 != index) continue;
+      
+       // found common particle
+       d1 -= part.Pt();
+      
+       if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
+         AliVParticle *MCpart = tracks2->GetParticle(index2);
+         AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
+                         iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
+         
+         d2 -= MCpart->Pt();
+       }
+
+       track2Found = kTRUE;
       }
     }
   }
-}
-
-//________________________________________________________________________
-Bool_t AliJetResponseMaker::FillHistograms()
-{
-  // Fill histograms.
-
-  fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
-  fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
-  if (fEventWeightHist)
-    fHistEventWeight[fPtHardBin]->Fill(fPythiaHeader->EventWeight());
 
-  const Int_t nMCJets = fMCJets->GetEntriesFast();
+  if (d1 < 0)
+    d1 = 0;
 
-  for (Int_t i = 0; i < nMCJets; i++) {
+  if (d2 < 0)
+    d2 = 0;
 
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
+  if (totalPt1 < 1)
+    d1 = -1;
+  else
+    d1 /= totalPt1;
 
-    if (!jet) {
-      AliError(Form("Could not receive jet %d", i));
-      continue;
-    }  
+  if (jet2->Pt() < 1)
+    d2 = -1;
+  else
+    d2 /= jet2->Pt();
+}
 
-    if (!AcceptJet(jet))
-      continue;
+//________________________________________________________________________
+void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
+{ 
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
+
+  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
+
+  AliParticleContainer *tracks1   = jets1->GetParticleContainer();
+  AliClusterContainer  *clusters1 = jets1->GetClusterContainer();
+  AliParticleContainer *tracks2   = jets2->GetParticleContainer();
+  AliClusterContainer  *clusters2 = jets2->GetClusterContainer();
+
+  // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
+  d1 = jet1->Pt();
+  d2 = jet2->Pt();
+
+  if (tracks1 && tracks2) {
+
+    for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
+      Int_t index2 = jet2->TrackAt(iTrack2);
+      for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
+       Int_t index1 = jet1->TrackAt(iTrack1);
+       if (index2 == index1) { // found common particle
+         AliVParticle *part1 = tracks1->GetParticle(index1);
+         if (!part1) {
+           AliWarning(Form("Could not find track %d!", index1));
+           continue;
+         }
+         AliVParticle *part2 = tracks2->GetParticle(index2);
+         if (!part2) {
+           AliWarning(Form("Could not find track %d!", index2));
+           continue;
+         }
+
+         d1 -= part1->Pt();
+         d2 -= part2->Pt();
+         break;
+       }
+      }
+    }
 
-    if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi)
-      continue;
+  }
 
-    if (jet->Pt() > fMaxBinPt)
-      continue;
+  if (clusters1 && clusters2) {
+
+    if (fUseCellsToMatch) {
+      const Int_t nClus1 = jet1->GetNumberOfClusters();
+
+      Int_t ncells1[nClus1];
+      UShort_t *cellsId1[nClus1];
+      Double_t *cellsFrac1[nClus1];
+      Int_t *sortedIndexes1[nClus1];
+      Double_t ptClus1[nClus1];
+      for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
+       Int_t index1 = jet1->ClusterAt(iClus1);
+       AliVCluster *clus1 = clusters1->GetCluster(index1);
+       if (!clus1) {
+         AliWarning(Form("Could not find cluster %d!", index1));
+         ncells1[iClus1] = 0;
+         cellsId1[iClus1] = 0;
+         cellsFrac1[iClus1] = 0;
+         sortedIndexes1[iClus1] = 0;
+         ptClus1[iClus1] = 0;
+         continue;
+       }
+       TLorentzVector part1;
+       clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
 
-    if (jet->MatchedJet()) {
+       ncells1[iClus1] = clus1->GetNCells();
+       cellsId1[iClus1] = clus1->GetCellsAbsId();
+       cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
+       sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
+       ptClus1[iClus1] = part1.Pt();
 
-      if (!AcceptBiasJet(jet->MatchedJet()) || 
-         jet->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
-         jet->MatchedJet()->Pt() > fMaxBinPt) {
-       fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
+       TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
       }
-      else {
-       fHistMatchingLevelMCPt->Fill(jet->ClosestJetDistance(), jet->Pt(), fEventWeight);
-
-       Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
-       Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
-       fHistClosestDeltaEtaPhiMCPt->Fill(deta, dphi, jet->Pt(), fEventWeight);
-
-       Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
-       fHistClosestDeltaPtMCPt->Fill(jet->Pt(), dpt, fEventWeight);
-       fHistClosestDeltaCorrPtMCPt->Fill(jet->Pt(), dpt - fRhoVal * jet->MatchedJet()->Area(), fEventWeight);
-
-       fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
-       fHistPartvsDetecCorrPt->Fill(jet->MatchedJet()->Pt() - fRhoVal * jet->MatchedJet()->Area(), jet->Pt(), fEventWeight);
+      
+      const Int_t nClus2 = jet2->GetNumberOfClusters();
+
+      const Int_t maxNcells2 = 11520;
+      Int_t sortedIndexes2[maxNcells2];
+      for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
+       Int_t index2 = jet2->ClusterAt(iClus2);
+       AliVCluster *clus2 =  clusters2->GetCluster(index2);
+       if (!clus2) {
+         AliWarning(Form("Could not find cluster %d!", index2));
+         continue;
+       }
+       Int_t ncells2 = clus2->GetNCells();
+       if (ncells2 >= maxNcells2) {
+         AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
+         continue;
+       }
+       UShort_t *cellsId2 = clus2->GetCellsAbsId();
+       Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
+
+       TLorentzVector part2;
+       clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
+       Double_t ptClus2 = part2.Pt();
+
+       TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
+
+       for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
+         if (sortedIndexes1[iClus1] == 0)
+           continue;
+         Int_t iCell1 = 0, iCell2 = 0;
+         while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
+           if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
+             d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
+             d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2;
+             iCell1++;
+             iCell2++;
+           }
+           else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) { 
+             iCell2++;
+           }
+           else {
+             iCell1++;
+           }
+         }
+       }
       }
     }
     else {
-      fHistNonMatchedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
-      fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
+      for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
+       Int_t index2 = jet2->ClusterAt(iClus2);
+       for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
+         Int_t index1 = jet1->ClusterAt(iClus1);
+         if (index2 == index1) { // found common particle
+           AliVCluster *clus1 = clusters1->GetCluster(index1);
+           if (!clus1) {
+             AliWarning(Form("Could not find cluster %d!", index1));
+             continue;
+           }
+           AliVCluster *clus2 =  clusters2->GetCluster(index2);
+           if (!clus2) {
+             AliWarning(Form("Could not find cluster %d!", index2));
+             continue;
+           }
+           TLorentzVector part1, part2;
+           clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
+           clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
+           
+           d1 -= part1.Pt();
+           d2 -= part2.Pt();
+           break;
+         }
+       }
+      }
     }
+  }
 
-    fHistMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
-    fHistMCJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
+  if (d1 < 0)
+    d1 = 0;
 
-    fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
-      
-    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
-      AliVParticle *track = jet->TrackAt(it, fMCTracks);
-      if (track)
-       fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
-    }
+  if (d2 < 0)
+    d2 = 0;
 
-    if (!AcceptBiasJet(jet))
-      continue;
-    if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
-      continue;
-    
-    fHistMCJetsPtAreaFiducial->Fill(jet->Area(), jet->Pt(), fEventWeight);
-    fHistMCJetsPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight);
-  }
+  if (jet1->Pt() > 0)
+    d1 /= jet1->Pt();
+  else
+    d1 = -1;
 
-  const Int_t nJets = fJets->GetEntriesFast();
+  if (jet2->Pt() > 0)
+    d2 /= jet2->Pt();
+  else
+    d2 = -1;
+}
 
-  for (Int_t i = 0; i < nJets; i++) {
+//________________________________________________________________________
+void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
+{
+  Double_t d1 = -1;
+  Double_t d2 = -1;
+
+  switch (matching) {
+  case kGeometrical:
+    GetGeometricalMatchingLevel(jet1,jet2,d1);
+    d2 = d1;
+    break;
+  case kMCLabel: // jet1 = detector level and jet2 = particle level!
+    GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
+    break;
+  case kSameCollections:
+    GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
+    break;
+  default:
+    ;
+  }
 
-    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
+  if (d1 >= 0) {
 
-    if (!jet) {
-      AliError(Form("Could not receive mc jet %d", i));
-      continue;
-    }  
+    if (d1 < jet1->ClosestJetDistance()) {
+      jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
+      jet1->SetClosestJet(jet2, d1);
+    }
+    else if (d1 < jet1->SecondClosestJetDistance()) {
+      jet1->SetSecondClosestJet(jet2, d1);
+    }
+  }
+  
+  if (d2 >= 0) {
     
-    if (!AcceptJet(jet))
-      continue;
-    if (!AcceptBiasJet(jet))
-      continue;
-    if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
-      continue;
-    if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
-      continue;
-
-    if (!jet->MatchedJet()) {
-      fHistNonMatchedJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
-      fHistNonMatchedJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight);
+    if (d2 < jet2->ClosestJetDistance()) {
+      jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
+      jet2->SetClosestJet(jet1, d2);
     }
+    else if (d2 < jet2->SecondClosestJetDistance()) {
+      jet2->SetSecondClosestJet(jet1, d2);
+    }
+  }
+}
 
-    fHistJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight);
-    fHistJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight);
+//________________________________________________________________________
+Bool_t AliJetResponseMaker::FillHistograms()
+{
+  // Fill histograms.
+
+  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
+  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
 
-    fHistJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
+  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
 
-    fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
+  AliEmcalJet* jet1 = 0;  
+  AliEmcalJet* jet2 = 0;
+  
+  jets2->ResetCurrentID();
+  while ((jet2 = jets2->GetNextJet())) {
 
-    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
-      AliVParticle *track = jet->TrackAt(it, fTracks);
-      if (track)
-       fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
+    AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
+
+    Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
+    Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
+
+    if (jets2->AcceptJet(jet2))
+      FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2, 
+                  corrpt2, jet2->MCPt(), 2);
+
+    jet1 = jet2->MatchedJet();
+
+    if (!jet1) continue;
+    if (!jets1->AcceptJet(jet1)) continue;
+    if (jet1->MCPt() < fMinJetMCPt) continue;
+
+    Double_t d=-1, ce1=-1, ce2=-1;
+    if (jet2->GetMatchingType() == kGeometrical) {
+      if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
+       GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
+      else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
+       GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
+      
+      d = jet2->ClosestJetDistance();
     }
+    else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
+      GetGeometricalMatchingLevel(jet1, jet2, d);
 
-    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
-      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
-      if (cluster) {
-       TLorentzVector nP;
-       cluster->GetMomentum(nP, fVertex);
-       fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight);
-      }
+      ce1 = jet1->ClosestJetDistance();
+      ce2 = jet2->ClosestJetDistance();
     }
+
+    Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
+    Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
+
+    FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(), 
+                      jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(), 
+                      jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2);
   }
 
+  jets1->ResetCurrentID();
+  while ((jet1 = jets1->GetNextAcceptJet())) {
+    if (jet1->MCPt() < fMinJetMCPt) continue;
+    AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
+
+    Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
+    Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
+
+    FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1, 
+                corrpt1, jet1->MCPt(), 1);
+  }
   return kTRUE;
 }