PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx

   1 // $Id$
   2 //
   3 // Emcal jet response matrix maker task.
   4 //
   5 // Author: S. Aiola
   6
   7 #include "AliJetResponseMaker.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TH2F.h>
  11 #include <THnSparse.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliAnalysisManager.h"
  15 #include "AliVCluster.h"
  16 #include "AliVTrack.h"
  17 #include "AliEmcalJet.h"
  18 #include "AliLog.h"
  19 #include "AliRhoParameter.h"
  20 #include "AliNamedArrayI.h"
  21 #include "AliJetContainer.h"
  22 #include "AliParticleContainer.h"
  23 #include "AliClusterContainer.h"
  24
  25 ClassImp(AliJetResponseMaker)
  26
  27 //________________________________________________________________________
  28 AliJetResponseMaker::AliJetResponseMaker() :
  29   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  30   fMatching(kNoMatching),
  31   fMatchingPar1(0),
  32   fMatchingPar2(0),
  33   fUseCellsToMatch(kFALSE),
  34   fMinJetMCPt(1),
  35   fHistoType(0),
  36   fDeltaPtAxis(0),
  37   fDeltaEtaDeltaPhiAxis(0),
  38   fNEFAxis(0),
  39   fZAxis(0),
  40   fIsJet1Rho(kFALSE),
  41   fIsJet2Rho(kFALSE),
  42   fHistLeadingJets1PtArea(0),
  43   fHistLeadingJets1CorrPtArea(0),
  44   fHistLeadingJets2PtArea(0),
  45   fHistLeadingJets2CorrPtArea(0),
  46   fHistLeadingJets2PtAreaAcceptance(0),
  47   fHistLeadingJets2CorrPtAreaAcceptance(0),
  48   fHistJets1(0),
  49   fHistJets2(0),
  50   fHistMatching(0),
  51   fHistJets1PhiEta(0),
  52   fHistJets1PtArea(0),
  53   fHistJets1CorrPtArea(0),
  54   fHistJets1NEFvsPt(0),
  55   fHistJets1CEFvsCEFPt(0),
  56   fHistJets1ZvsPt(0),
  57   fHistJets2PhiEta(0),
  58   fHistJets2PtArea(0),
  59   fHistJets2CorrPtArea(0),
  60   fHistJets2PhiEtaAcceptance(0),
  61   fHistJets2PtAreaAcceptance(0),
  62   fHistJets2CorrPtAreaAcceptance(0),
  63   fHistJets2NEFvsPt(0),
  64   fHistJets2CEFvsCEFPt(0),
  65   fHistJets2ZvsPt(0),
  66   fHistCommonEnergy1vsJet1Pt(0),
  67   fHistCommonEnergy2vsJet2Pt(0),
  68   fHistDistancevsJet1Pt(0),
  69   fHistDistancevsJet2Pt(0),
  70   fHistDistancevsCommonEnergy1(0),
  71   fHistDistancevsCommonEnergy2(0),
  72   fHistCommonEnergy1vsCommonEnergy2(0),
  73   fHistDeltaEtaDeltaPhi(0),
  74   fHistJet2PtOverJet1PtvsJet2Pt(0),
  75   fHistJet1PtOverJet2PtvsJet1Pt(0),
  76   fHistDeltaPtvsJet1Pt(0),
  77   fHistDeltaPtvsJet2Pt(0),
  78   fHistDeltaPtOverJet1PtvsJet1Pt(0),
  79   fHistDeltaPtOverJet2PtvsJet2Pt(0),
  80   fHistDeltaPtvsDistance(0),
  81   fHistDeltaPtvsCommonEnergy1(0),
  82   fHistDeltaPtvsCommonEnergy2(0),
  83   fHistDeltaPtvsArea1(0),
  84   fHistDeltaPtvsArea2(0),
  85   fHistDeltaPtvsDeltaArea(0),
  86   fHistJet1PtvsJet2Pt(0),
  87   fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
  88   fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
  89   fHistDeltaCorrPtvsJet1CorrPt(0),
  90   fHistDeltaCorrPtvsJet2CorrPt(0),
  91   fHistDeltaCorrPtvsDistance(0),
  92   fHistDeltaCorrPtvsCommonEnergy1(0),
  93   fHistDeltaCorrPtvsCommonEnergy2(0),
  94   fHistDeltaCorrPtvsArea1(0),
  95   fHistDeltaCorrPtvsArea2(0),
  96   fHistDeltaCorrPtvsDeltaArea(0),
  97   fHistJet1CorrPtvsJet2CorrPt(0),
  98   fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
  99   fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
 100   fHistDeltaMCPtvsJet1MCPt(0),
 101   fHistDeltaMCPtvsJet2Pt(0),
 102   fHistDeltaMCPtvsDistance(0),
 103   fHistDeltaMCPtvsCommonEnergy1(0),
 104   fHistDeltaMCPtvsCommonEnergy2(0),
 105   fHistDeltaMCPtvsArea1(0),
 106   fHistDeltaMCPtvsArea2(0),
 107   fHistDeltaMCPtvsDeltaArea(0),
 108   fHistJet1MCPtvsJet2Pt(0)
 109 {
 110   // Default constructor.
 111
 112   SetMakeGeneralHistograms(kTRUE);
 113 }
 114
 115 //________________________________________________________________________
 116 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
 117   AliAnalysisTaskEmcalJet(name, kTRUE),
 118   fMatching(kNoMatching),
 119   fMatchingPar1(0),
 120   fMatchingPar2(0),
 121   fUseCellsToMatch(kFALSE),
 122   fMinJetMCPt(1),
 123   fHistoType(0),
 124   fDeltaPtAxis(0),
 125   fDeltaEtaDeltaPhiAxis(0),
 126   fNEFAxis(0),
 127   fZAxis(0),
 128   fIsJet1Rho(kFALSE),
 129   fIsJet2Rho(kFALSE),
 130   fHistLeadingJets1PtArea(0),
 131   fHistLeadingJets1CorrPtArea(0),
 132   fHistLeadingJets2PtArea(0),
 133   fHistLeadingJets2CorrPtArea(0),
 134   fHistLeadingJets2PtAreaAcceptance(0),
 135   fHistLeadingJets2CorrPtAreaAcceptance(0),
 136   fHistJets1(0),
 137   fHistJets2(0),
 138   fHistMatching(0),
 139   fHistJets1PhiEta(0),
 140   fHistJets1PtArea(0),
 141   fHistJets1CorrPtArea(0),
 142   fHistJets1NEFvsPt(0),
 143   fHistJets1CEFvsCEFPt(0),
 144   fHistJets1ZvsPt(0),
 145   fHistJets2PhiEta(0),
 146   fHistJets2PtArea(0),
 147   fHistJets2CorrPtArea(0),
 148   fHistJets2PhiEtaAcceptance(0),
 149   fHistJets2PtAreaAcceptance(0),
 150   fHistJets2CorrPtAreaAcceptance(0),
 151   fHistJets2NEFvsPt(0),
 152   fHistJets2CEFvsCEFPt(0),
 153   fHistJets2ZvsPt(0),
 154   fHistCommonEnergy1vsJet1Pt(0),
 155   fHistCommonEnergy2vsJet2Pt(0),
 156   fHistDistancevsJet1Pt(0),
 157   fHistDistancevsJet2Pt(0),
 158   fHistDistancevsCommonEnergy1(0),
 159   fHistDistancevsCommonEnergy2(0),
 160   fHistCommonEnergy1vsCommonEnergy2(0),
 161   fHistDeltaEtaDeltaPhi(0),
 162   fHistJet2PtOverJet1PtvsJet2Pt(0),
 163   fHistJet1PtOverJet2PtvsJet1Pt(0),
 164   fHistDeltaPtvsJet1Pt(0),
 165   fHistDeltaPtvsJet2Pt(0),
 166   fHistDeltaPtOverJet1PtvsJet1Pt(0),
 167   fHistDeltaPtOverJet2PtvsJet2Pt(0),
 168   fHistDeltaPtvsDistance(0),
 169   fHistDeltaPtvsCommonEnergy1(0),
 170   fHistDeltaPtvsCommonEnergy2(0),
 171   fHistDeltaPtvsArea1(0),
 172   fHistDeltaPtvsArea2(0),
 173   fHistDeltaPtvsDeltaArea(0),
 174   fHistJet1PtvsJet2Pt(0),
 175   fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
 176   fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
 177   fHistDeltaCorrPtvsJet1CorrPt(0),
 178   fHistDeltaCorrPtvsJet2CorrPt(0),
 179   fHistDeltaCorrPtvsDistance(0),
 180   fHistDeltaCorrPtvsCommonEnergy1(0),
 181   fHistDeltaCorrPtvsCommonEnergy2(0),
 182   fHistDeltaCorrPtvsArea1(0),
 183   fHistDeltaCorrPtvsArea2(0),
 184   fHistDeltaCorrPtvsDeltaArea(0),
 185   fHistJet1CorrPtvsJet2CorrPt(0),
 186   fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
 187   fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
 188   fHistDeltaMCPtvsJet1MCPt(0),
 189   fHistDeltaMCPtvsJet2Pt(0),
 190   fHistDeltaMCPtvsDistance(0),
 191   fHistDeltaMCPtvsCommonEnergy1(0),
 192   fHistDeltaMCPtvsCommonEnergy2(0),
 193   fHistDeltaMCPtvsArea1(0),
 194   fHistDeltaMCPtvsArea2(0),
 195   fHistDeltaMCPtvsDeltaArea(0),
 196   fHistJet1MCPtvsJet2Pt(0)
 197 {
 198   // Standard constructor.
 199
 200   SetMakeGeneralHistograms(kTRUE);
 201 }
 202
 203 //________________________________________________________________________
 204 AliJetResponseMaker::~AliJetResponseMaker()
 205 {
 206   // Destructor
 207 }
 208
 209
 210 //________________________________________________________________________
 211 void AliJetResponseMaker::AllocateTH2()
 212 {
 213   // Allocate TH2 histograms.
 214
 215   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
 216   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
 217   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
 218   fOutput->Add(fHistJets1PhiEta);
 219
 220   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
 221   fHistJets1PtArea->GetXaxis()->SetTitle("area");
 222   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 223   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
 224   fOutput->Add(fHistJets1PtArea);
 225
 226   if (fIsJet1Rho) {
 227     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea",
 228                                     fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 229     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
 230     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
 231     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 232     fOutput->Add(fHistJets1CorrPtArea);
 233   }
 234
 235   fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 236   fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
 237   fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 238   fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
 239   fOutput->Add(fHistJets1ZvsPt);
 240
 241   fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 242   fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
 243   fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 244   fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
 245   fOutput->Add(fHistJets1NEFvsPt);
 246
 247   fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 248   fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
 249   fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
 250   fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
 251   fOutput->Add(fHistJets1CEFvsCEFPt);
 252
 253   // Jets 2 histograms
 254
 255   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
 256   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
 257   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
 258   fOutput->Add(fHistJets2PhiEta);
 259
 260   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
 261   fHistJets2PtArea->GetXaxis()->SetTitle("area");
 262   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 263   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
 264   fOutput->Add(fHistJets2PtArea);
 265
 266   if (fIsJet2Rho) {
 267     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea",
 268                                     fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 269     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
 270     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 271     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 272     fOutput->Add(fHistJets2CorrPtArea);
 273   }
 274
 275   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
 276   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
 277   fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
 278   fOutput->Add(fHistJets2PhiEtaAcceptance);
 279
 280   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance",
 281                                         fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
 282   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
 283   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 284   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
 285   fOutput->Add(fHistJets2PtAreaAcceptance);
 286
 287   if (fIsJet2Rho) {
 288     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance",
 289                                               fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 290     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
 291     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 292     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
 293     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
 294   }
 295
 296   fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 297   fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
 298   fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 299   fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
 300   fOutput->Add(fHistJets2ZvsPt);
 301
 302   fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 303   fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
 304   fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 305   fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
 306   fOutput->Add(fHistJets2NEFvsPt);
 307
 308   fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 309   fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
 310   fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
 311   fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
 312   fOutput->Add(fHistJets2CEFvsCEFPt);
 313
 314   // Matching histograms
 315
 316   fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 317   fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
 318   fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
 319   fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
 320   fOutput->Add(fHistCommonEnergy1vsJet1Pt);
 321
 322   fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 323   fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
 324   fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 325   fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
 326   fOutput->Add(fHistCommonEnergy2vsJet2Pt);
 327
 328   fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 329   fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
 330   fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");
 331   fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
 332   fOutput->Add(fHistDistancevsJet1Pt);
 333
 334   fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 335   fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
 336   fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 337   fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
 338   fOutput->Add(fHistDistancevsJet2Pt);
 339
 340   fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
 341   fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
 342   fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");
 343   fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
 344   fOutput->Add(fHistDistancevsCommonEnergy1);
 345
 346   fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
 347   fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
 348   fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
 349   fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
 350   fOutput->Add(fHistDistancevsCommonEnergy2);
 351
 352   fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
 353   fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
 354   fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");
 355   fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
 356   fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
 357
 358   fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
 359   fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
 360   fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");
 361   fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
 362   fOutput->Add(fHistDeltaEtaDeltaPhi);
 363
 364   fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
 365   fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 366   fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
 367   fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 368   fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
 369
 370   fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
 371   fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
 372   fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
 373   fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
 374   fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
 375
 376   fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt",
 377                                   fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 378   fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
 379   fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 380   fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
 381   fOutput->Add(fHistDeltaPtvsJet1Pt);
 382
 383   fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt",
 384                                   fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 385   fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 386   fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 387   fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 388   fOutput->Add(fHistDeltaPtvsJet2Pt);
 389
 390   fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt",
 391                                             fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
 392   fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");
 393   fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}");
 394   fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts");
 395   fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt);
 396
 397   fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt",
 398                                             fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
 399   fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 400   fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
 401   fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 402   fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt);
 403
 404   fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance",
 405                                     fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 406   fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");
 407   fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 408   fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
 409   fOutput->Add(fHistDeltaPtvsDistance);
 410
 411   fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1",
 412                                          fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 413   fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
 414   fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 415   fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
 416   fOutput->Add(fHistDeltaPtvsCommonEnergy1);
 417
 418   fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2",
 419                                          fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 420   fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
 421   fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 422   fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
 423   fOutput->Add(fHistDeltaPtvsCommonEnergy2);
 424
 425   fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1",
 426                                  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 427   fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
 428   fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 429   fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
 430   fOutput->Add(fHistDeltaPtvsArea1);
 431
 432   fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2",
 433                                  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 434   fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
 435   fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 436   fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
 437   fOutput->Add(fHistDeltaPtvsArea2);
 438
 439   fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea",
 440                                      fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 441   fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
 442   fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 443   fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
 444   fOutput->Add(fHistDeltaPtvsDeltaArea);
 445
 446   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 447   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
 448   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 449   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 450   fOutput->Add(fHistJet1PtvsJet2Pt);
 451
 452   if (fIsJet1Rho || fIsJet2Rho) {
 453     if (!fIsJet1Rho)
 454       fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
 455                                               fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 456     else
 457       fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt",
 458                                               2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 459
 460     fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
 461     fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 462     fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
 463     fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt);
 464
 465     if (!fIsJet2Rho)
 466       fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
 467                                               fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 468     else
 469       fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt",
 470                                               2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 471
 472     fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
 473     fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 474     fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
 475     fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt);
 476
 477     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt",
 478                                                           2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
 479     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
 480     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}");
 481     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
 482     fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt);
 483
 484     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt",
 485                                                           2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
 486     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");
 487     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}");
 488     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
 489     fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt);
 490
 491     fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance",
 492                                           fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 493     fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");
 494     fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 495     fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
 496     fOutput->Add(fHistDeltaCorrPtvsDistance);
 497
 498     fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1",
 499                                                fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 500     fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
 501     fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 502     fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
 503     fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1);
 504
 505     fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2",
 506                                                fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 507     fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
 508     fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 509     fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
 510     fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
 511
 512     fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1",
 513                                        fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 514     fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
 515     fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 516     fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
 517     fOutput->Add(fHistDeltaCorrPtvsArea1);
 518
 519     fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2",
 520                                        fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 521     fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
 522     fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 523     fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
 524     fOutput->Add(fHistDeltaCorrPtvsArea2);
 525
 526     fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea",
 527                                            fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 528     fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
 529     fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
 530     fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
 531     fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
 532
 533     if (!fIsJet1Rho)
 534       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
 535                                              fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 536     else if (!fIsJet2Rho)
 537       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
 538                                              2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 539     else
 540       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt",
 541                                              2*fNbins, -fMaxBinPt, fMaxBinPt,
 542                                              2*fNbins, -fMaxBinPt, fMaxBinPt);
 543     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
 544     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
 545     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
 546     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
 547   }
 548
 549   if (fIsEmbedded) {
 550     fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt",
 551                                         fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 552     fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
 553     fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 554     fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
 555     fOutput->Add(fHistDeltaMCPtvsJet1MCPt);
 556
 557     fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt",
 558                                       fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 559     fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 560     fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 561     fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 562     fOutput->Add(fHistDeltaMCPtvsJet2Pt);
 563
 564     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt",
 565                                                     fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
 566     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
 567     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}");
 568     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
 569     fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt);
 570
 571     fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt",
 572                                                 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
 573     fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 574     fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
 575     fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 576     fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt);
 577
 578     fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance",
 579                                         fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 580     fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");
 581     fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 582     fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
 583     fOutput->Add(fHistDeltaMCPtvsDistance);
 584
 585     fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1",
 586                                              fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 587     fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");
 588     fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 589     fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
 590     fOutput->Add(fHistDeltaMCPtvsCommonEnergy1);
 591
 592     fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2",
 593                                              fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 594     fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");
 595     fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 596     fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
 597     fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
 598
 599     fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1",
 600                                      fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 601     fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
 602     fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 603     fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
 604     fOutput->Add(fHistDeltaMCPtvsArea1);
 605
 606     fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2",
 607                                      fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 608     fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
 609     fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 610     fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
 611     fOutput->Add(fHistDeltaMCPtvsArea2);
 612
 613     fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea",
 614                                          fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 615     fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
 616     fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 617     fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
 618     fOutput->Add(fHistDeltaMCPtvsDeltaArea);
 619
 620     fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt",
 621                                      fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 622     fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
 623     fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 624     fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 625     fOutput->Add(fHistJet1MCPtvsJet2Pt);
 626   }
 627 }
 628
 629 //________________________________________________________________________
 630 void AliJetResponseMaker::AllocateTHnSparse()
 631 {
 632   // Allocate THnSparse histograms.
 633
 634   TString title[25]= {""};
 635   Int_t nbins[25]  = {0};
 636   Double_t min[25] = {0.};
 637   Double_t max[25] = {0.};
 638   Int_t dim = 0;
 639
 640   title[dim] = "#phi";
 641   nbins[dim] = fNbins/4;
 642   min[dim] = 0;
 643   max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
 644   dim++;
 645
 646   title[dim] = "#eta";
 647   nbins[dim] = fNbins/4;
 648   min[dim] = -1;
 649   max[dim] = 1;
 650   dim++;
 651
 652   title[dim] = "p_{T}";
 653   nbins[dim] = fNbins;
 654   min[dim] = 0;
 655   max[dim] = 250;
 656   dim++;
 657
 658   title[dim] = "A_{jet}";
 659   nbins[dim] = fNbins/4;
 660   min[dim] = 0;
 661   max[dim] = 1.5;
 662   dim++;
 663
 664   if (fNEFAxis) {
 665     title[dim] = "NEF";
 666     nbins[dim] = fNbins/4;
 667     min[dim] = 0;
 668     max[dim] = 1.2;
 669     dim++;
 670   }
 671
 672   if (fZAxis) {
 673     title[dim] = "Z";
 674     nbins[dim] = fNbins/4;
 675     min[dim] = 0;
 676     max[dim] = 1.2;
 677     dim++;
 678   }
 679
 680   title[dim] = "p_{T,particle}^{leading} (GeV/c)";
 681   nbins[dim] = 120;
 682   min[dim] = 0;
 683   max[dim] = 120;
 684   dim++;
 685
 686   Int_t dim1 = dim, dim2 = dim;
 687
 688   if (fIsJet1Rho) {
 689     title[dim1] = "p_{T}^{corr}";
 690     nbins[dim1] = fNbins*2;
 691     min[dim1] = -250;
 692     max[dim1] = 250;
 693     dim1++;
 694   }
 695
 696   if (fIsEmbedded) {
 697     title[dim1] = "p_{T}^{MC}";
 698     nbins[dim1] = fNbins;
 699     min[dim1] = 0;
 700     max[dim1] = 250;
 701     dim1++;
 702   }
 703
 704   fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
 705   for (Int_t i = 0; i < dim1; i++)
 706     fHistJets1->GetAxis(i)->SetTitle(title[i]);
 707   fOutput->Add(fHistJets1);
 708
 709   if (fIsJet2Rho) {
 710     title[dim2] = "p_{T}^{corr}";
 711     nbins[dim2] = fNbins*2;
 712     min[dim2] = -250;
 713     max[dim2] = 250;
 714     dim2++;
 715   }
 716
 717   fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
 718   for (Int_t i = 0; i < dim2; i++)
 719     fHistJets2->GetAxis(i)->SetTitle(title[i]);
 720   fOutput->Add(fHistJets2);
 721
 722   // Matching
 723
 724   dim = 0;
 725
 726   title[dim] = "p_{T,1}";
 727   nbins[dim] = fNbins;
 728   min[dim] = 0;
 729   max[dim] = 250;
 730   dim++;
 731
 732   title[dim] = "p_{T,2}";
 733   nbins[dim] = fNbins;
 734   min[dim] = 0;
 735   max[dim] = 250;
 736   dim++;
 737
 738   title[dim] = "A_{jet,1}";
 739   nbins[dim] = fNbins/4;
 740   min[dim] = 0;
 741   max[dim] = 1.5;
 742   dim++;
 743
 744   title[dim] = "A_{jet,2}";
 745   nbins[dim] = fNbins/4;
 746   min[dim] = 0;
 747   max[dim] = 1.5;
 748   dim++;
 749
 750   title[dim] = "distance";
 751   nbins[dim] = fNbins/2;
 752   min[dim] = 0;
 753   max[dim] = 1.2;
 754   dim++;
 755
 756   title[dim] = "CE1";
 757   nbins[dim] = fNbins/2;
 758   min[dim] = 0;
 759   max[dim] = 1.2;
 760   dim++;
 761
 762   title[dim] = "CE2";
 763   nbins[dim] = fNbins/2;
 764   min[dim] = 0;
 765   max[dim] = 1.2;
 766   dim++;
 767
 768   title[dim] = "p_{T,particle,1}^{leading} (GeV/c)";
 769   nbins[dim] = 120;
 770   min[dim] = 0;
 771   max[dim] = 120;
 772   dim++;
 773
 774   title[dim] = "p_{T,particle,2}^{leading} (GeV/c)";
 775   nbins[dim] = 120;
 776   min[dim] = 0;
 777   max[dim] = 120;
 778   dim++;
 779
 780   if (fDeltaPtAxis) {
 781     title[dim] = "#deltaA_{jet}";
 782     nbins[dim] = fNbins/2;
 783     min[dim] = -1;
 784     max[dim] = 1;
 785     dim++;
 786
 787     title[dim] = "#deltap_{T}";
 788     nbins[dim] = fNbins*2;
 789     min[dim] = -250;
 790     max[dim] = 250;
 791     dim++;
 792   }
 793   if (fIsJet1Rho) {
 794     title[dim] = "p_{T,1}^{corr}";
 795     nbins[dim] = fNbins*2;
 796     min[dim] = -250;
 797     max[dim] = 250;
 798     dim++;
 799   }
 800   if (fIsJet2Rho) {
 801     title[dim] = "p_{T,2}^{corr}";
 802     nbins[dim] = fNbins*2;
 803     min[dim] = -250;
 804     max[dim] = 250;
 805     dim++;
 806   }
 807   if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
 808     title[dim] = "#deltap_{T}^{corr}";
 809     nbins[dim] = fNbins*2;
 810     min[dim] = -250;
 811     max[dim] = 250;
 812     dim++;
 813   }
 814   if (fDeltaEtaDeltaPhiAxis) {
 815     title[dim] = "#delta#eta";
 816     nbins[dim] = fNbins/2;
 817     min[dim] = -1;
 818     max[dim] = 1;
 819     dim++;
 820
 821     title[dim] = "#delta#phi";
 822     nbins[dim] = fNbins/2;
 823     min[dim] = -TMath::Pi()/2;
 824     max[dim] = TMath::Pi()*3/2;
 825     dim++;
 826   }
 827   if (fIsEmbedded) {
 828     title[dim] = "p_{T,1}^{MC}";
 829     nbins[dim] = fNbins;
 830     min[dim] = 0;
 831     max[dim] = 250;
 832     dim++;
 833
 834     if (fDeltaPtAxis) {
 835       title[dim] = "#deltap_{T}^{MC}";
 836       nbins[dim] = fNbins*2;
 837       min[dim] = -250;
 838       max[dim] = 250;
 839       dim++;
 840     }
 841   }
 842
 843   if (fNEFAxis) {
 844     title[dim] = "NEF_{1}";
 845     nbins[dim] = fNbins/4;
 846     min[dim] = 0;
 847     max[dim] = 1.2;
 848     dim++;
 849
 850     title[dim] = "NEF_{2}";
 851     nbins[dim] = fNbins/4;
 852     min[dim] = 0;
 853     max[dim] = 1.2;
 854     dim++;
 855   }
 856
 857   if (fZAxis) {
 858     title[dim] = "Z_{1}";
 859     nbins[dim] = fNbins/4;
 860     min[dim] = 0;
 861     max[dim] = 1.2;
 862     dim++;
 863
 864     title[dim] = "Z_{2}";
 865     nbins[dim] = fNbins/4;
 866     min[dim] = 0;
 867     max[dim] = 1.2;
 868     dim++;
 869   }
 870
 871   fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
 872
 873   for (Int_t i = 0; i < dim; i++)
 874     fHistMatching->GetAxis(i)->SetTitle(title[i]);
 875
 876   fOutput->Add(fHistMatching);
 877 }
 878
 879 //________________________________________________________________________
 880 void AliJetResponseMaker::UserCreateOutputObjects()
 881 {
 882   // Create user objects.
 883
 884   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 885
 886   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
 887   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
 888
 889   if (!jets1 || !jets2) return;
 890
 891   if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
 892   else fIsJet1Rho = kTRUE;
 893   if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
 894   else fIsJet2Rho = kTRUE;
 895
 896   if (fHistoType==0)
 897     AllocateTH2();
 898   else
 899     AllocateTHnSparse();
 900
 901   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 902 }
 903
 904 //________________________________________________________________________
 905 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)
 906 {
 907   if (fHistoType==1) {
 908     THnSparse *histo = 0;
 909     if (Set==1)
 910       histo = fHistJets1;
 911     else if (Set==2)
 912       histo = fHistJets2;
 913
 914     if (!histo)
 915       return;
 916
 917     Double_t contents[20]={0};
 918
 919     for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
 920       TString title(histo->GetAxis(i)->GetTitle());
 921       if (title=="#phi")
 922         contents[i] = Phi;
 923       else if (title=="#eta")
 924         contents[i] = Eta;
 925       else if (title=="p_{T}")
 926         contents[i] = Pt;
 927       else if (title=="A_{jet}")
 928         contents[i] = A;
 929       else if (title=="NEF")
 930         contents[i] = NEF;
 931       else if (title=="Z")
 932         contents[i] = LeadingPt/Pt;
 933       else if (title=="p_{T}^{corr}")
 934         contents[i] = CorrPt;
 935       else if (title=="p_{T}^{MC}")
 936         contents[i] = MCPt;
 937       else if (title=="p_{T,particle}^{leading} (GeV/c)")
 938         contents[i] = LeadingPt;
 939       else
 940         AliWarning(Form("Unable to fill dimension %s!",title.Data()));
 941     }
 942
 943     histo->Fill(contents);
 944   }
 945   else {
 946     if (Set == 1) {
 947       fHistJets1PtArea->Fill(A, Pt);
 948       fHistJets1PhiEta->Fill(Eta, Phi);
 949       fHistJets1ZvsPt->Fill(LeadingPt/Pt, Pt);
 950       fHistJets1NEFvsPt->Fill(NEF, Pt);
 951       fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
 952       if (fIsJet1Rho)
 953         fHistJets1CorrPtArea->Fill(A, CorrPt);
 954     }
 955     else if (Set == 2) {
 956       fHistJets2PtArea->Fill(A, Pt);
 957       fHistJets2PhiEta->Fill(Eta, Phi);
 958       if (fIsJet2Rho)
 959         fHistJets2CorrPtArea->Fill(A, CorrPt);
 960     }
 961     else if (Set == 3) {
 962       fHistJets2PtAreaAcceptance->Fill(A, Pt);
 963       fHistJets2PhiEtaAcceptance->Fill(Eta, Phi);
 964       fHistJets2ZvsPt->Fill(LeadingPt/Pt, Pt);
 965       fHistJets2NEFvsPt->Fill(NEF, Pt);
 966       fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
 967       if (fIsJet2Rho)
 968         fHistJets2CorrPtAreaAcceptance->Fill(A, CorrPt);
 969     }
 970   }
 971 }
 972
 973 //________________________________________________________________________
 974 void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2,
 975                                              Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2,
 976                                              Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t LeadingPt1, Double_t LeadingPt2)
 977 {
 978   if (fHistoType==1) {
 979     Double_t contents[20]={0};
 980
 981     for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
 982       TString title(fHistMatching->GetAxis(i)->GetTitle());
 983       if (title=="p_{T,1}")
 984         contents[i] = Pt1;
 985       else if (title=="p_{T,2}")
 986         contents[i] = Pt2;
 987       else if (title=="A_{jet,1}")
 988         contents[i] = A1;
 989       else if (title=="A_{jet,2}")
 990         contents[i] = A2;
 991       else if (title=="distance")
 992         contents[i] = d;
 993       else if (title=="CE1")
 994         contents[i] = CE1;
 995       else if (title=="CE2")
 996         contents[i] = CE2;
 997       else if (title=="#deltaA_{jet}")
 998         contents[i] = A1-A2;
 999       else if (title=="#deltap_{T}")
1000         contents[i] = Pt1-Pt2;
1001       else if (title=="#delta#eta")
1002         contents[i] = Eta1-Eta2;
1003       else if (title=="#delta#phi")
1004         contents[i] = Phi1-Phi2;
1005       else if (title=="p_{T,1}^{corr}")
1006         contents[i] = CorrPt1;
1007       else if (title=="p_{T,2}^{corr}")
1008         contents[i] = CorrPt2;
1009       else if (title=="#deltap_{T}^{corr}")
1010         contents[i] = CorrPt1-CorrPt2;
1011       else if (title=="p_{T,1}^{MC}")
1012         contents[i] = MCPt1;
1013       else if (title=="#deltap_{T}^{MC}")
1014         contents[i] = MCPt1-Pt2;
1015       else if (title=="NEF_{1}")
1016         contents[i] = NEF1;
1017       else if (title=="NEF_{2}")
1018         contents[i] = NEF2;
1019       else if (title=="Z_{1}")
1020         contents[i] = LeadingPt1/Pt1;
1021       else if (title=="Z_{2}")
1022         contents[i] = LeadingPt2/Pt2;
1023       else if (title=="p_{T,particle,1}^{leading} (GeV/c)")
1024         contents[i] = LeadingPt1;
1025       else if (title=="p_{T,particle,2}^{leading} (GeV/c)")
1026         contents[i] = LeadingPt2;
1027       else
1028         AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1029     }
1030
1031     fHistMatching->Fill(contents);
1032   }
1033   else {
1034     fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
1035     fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
1036
1037     fHistDistancevsJet1Pt->Fill(d, Pt1);
1038     fHistDistancevsJet2Pt->Fill(d, Pt2);
1039
1040     fHistDistancevsCommonEnergy1->Fill(d, CE1);
1041     fHistDistancevsCommonEnergy2->Fill(d, CE2);
1042     fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
1043
1044     fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
1045
1046     fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
1047     fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
1048
1049     Double_t dpt = Pt1 - Pt2;
1050     fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
1051     fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
1052     fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
1053     fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
1054
1055     fHistDeltaPtvsDistance->Fill(d, dpt);
1056     fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
1057     fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
1058
1059     fHistDeltaPtvsArea1->Fill(A1, dpt);
1060     fHistDeltaPtvsArea2->Fill(A2, dpt);
1061
1062     Double_t darea = A1 - A2;
1063     fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
1064
1065     fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
1066
1067     if (fIsJet1Rho || fIsJet2Rho) {
1068       Double_t dcorrpt = CorrPt1 - CorrPt2;
1069       fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
1070       fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
1071       fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
1072       fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
1073       fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1074       fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
1075       fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
1076       fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
1077       fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
1078       fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1079       fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
1080     }
1081
1082     if (fIsEmbedded) {
1083       Double_t dmcpt = MCPt1 - Pt2;
1084       fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
1085       fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
1086       fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
1087       fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
1088       fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1089       fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
1090       fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
1091       fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
1092       fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
1093       fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1094       fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
1095     }
1096   }
1097 }
1098
1099 //________________________________________________________________________
1100 void AliJetResponseMaker::ExecOnce()
1101 {
1102   // Execute once.
1103
1104   AliAnalysisTaskEmcalJet::ExecOnce();
1105
1106   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1107   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1108
1109   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1110
1111   if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
1112     if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
1113       AliWarning("Changing matching type from MC label to same collection...");
1114       fMatching = kSameCollections;
1115     }
1116     else {
1117       AliWarning("Changing matching type from MC label to geometrical...");
1118       fMatching = kGeometrical;
1119     }
1120   }
1121   else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
1122     if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
1123       AliWarning("Changing matching type from same collection to MC label...");
1124       fMatching = kMCLabel;
1125     }
1126     else {
1127       AliWarning("Changing matching type from same collection to geometrical...");
1128       fMatching = kGeometrical;
1129     }
1130   }
1131 }
1132
1133 //________________________________________________________________________
1134 Bool_t AliJetResponseMaker::Run()
1135 {
1136   // Find the closest jets
1137
1138   if (fMatching == kNoMatching)
1139     return kTRUE;
1140   else
1141     return DoJetMatching();
1142 }
1143
1144 //________________________________________________________________________
1145 Bool_t AliJetResponseMaker::DoJetMatching()
1146 {
1147   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1148   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1149
1150   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1151
1152   DoJetLoop();
1153
1154   AliEmcalJet* jet1 = 0;
1155
1156   jets1->ResetCurrentID();
1157   while ((jet1 = jets1->GetNextJet())) {
1158
1159     AliEmcalJet *jet2 = jet1->ClosestJet();
1160
1161     if (!jet2) continue;
1162     if (jet2->ClosestJet() != jet1) continue;
1163     if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
1164
1165     // Matched jet found
1166     jet1->SetMatchedToClosest(fMatching);
1167     jet2->SetMatchedToClosest(fMatching);
1168     AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f",
1169                     jet1->Pt(), jet1->Eta(), jet1->Phi(),
1170                     jet2->Pt(), jet2->Eta(), jet2->Phi()));
1171   }
1172
1173   return kTRUE;
1174 }
1175
1176 //________________________________________________________________________
1177 void AliJetResponseMaker::DoJetLoop()
1178 {
1179   // Do the jet loop.
1180
1181   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1182   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1183
1184   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1185
1186   AliEmcalJet* jet1 = 0;
1187   AliEmcalJet* jet2 = 0;
1188
1189   jets2->ResetCurrentID();
1190   while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
1191
1192   jets1->ResetCurrentID();
1193   while ((jet1 = jets1->GetNextJet())) {
1194     jet1->ResetMatching();
1195
1196     if (jet1->MCPt() < fMinJetMCPt) continue;
1197
1198     jets2->ResetCurrentID();
1199     while ((jet2 = jets2->GetNextJet())) {
1200       SetMatchingLevel(jet1, jet2, fMatching);
1201     } // jet2 loop
1202   } // jet1 loop
1203 }
1204
1205 //________________________________________________________________________
1206 void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
1207 {
1208   Double_t deta = jet2->Eta() - jet1->Eta();
1209   Double_t dphi = jet2->Phi() - jet1->Phi();
1210   d = TMath::Sqrt(deta * deta + dphi * dphi);
1211 }
1212
1213 //________________________________________________________________________
1214 void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1215 {
1216   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1217   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1218
1219   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1220
1221   AliParticleContainer *tracks1   = jets1->GetParticleContainer();
1222   AliClusterContainer  *clusters1 = jets1->GetClusterContainer();
1223   AliParticleContainer *tracks2   = jets2->GetParticleContainer();
1224
1225   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1226   d1 = jet1->Pt();
1227   d2 = jet2->Pt();
1228   Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1229
1230   // remove completely tracks that are not MC particles (label == 0)
1231   if (tracks1 && tracks1->GetArray()) {
1232     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1233       AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1234       if (!track) {
1235         AliWarning(Form("Could not find track %d!", iTrack));
1236         continue;
1237       }
1238
1239       Int_t MClabel = TMath::Abs(track->GetLabel());
1240       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1241       if (MClabel != 0) continue;
1242
1243       // this is not a MC particle; remove it completely
1244       AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1245       totalPt1 -= track->Pt();
1246       d1 -= track->Pt();
1247     }
1248   }
1249
1250   // remove completely clusters that are not MC particles (label == 0)
1251   if (fUseCellsToMatch && fCaloCells) {
1252     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1253       AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1254       if (!clus) {
1255         AliWarning(Form("Could not find cluster %d!", iClus));
1256         continue;
1257       }
1258       TLorentzVector part;
1259       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1260
1261       for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1262         Int_t cellId = clus->GetCellAbsId(iCell);
1263         Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1264
1265         Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1266         if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1267         if (MClabel != 0) continue;
1268
1269         // this is not a MC particle; remove it completely
1270         AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1271         totalPt1 -= part.Pt() * cellFrac;
1272         d1 -= part.Pt() * cellFrac;
1273       }
1274     }
1275   }
1276   else {
1277     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1278       AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1279       if (!clus) {
1280         AliWarning(Form("Could not find cluster %d!", iClus));
1281         continue;
1282       }
1283       TLorentzVector part;
1284       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1285
1286       Int_t MClabel = TMath::Abs(clus->GetLabel());
1287       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1288       if (MClabel != 0) continue;
1289
1290       // this is not a MC particle; remove it completely
1291       AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1292       totalPt1 -= part.Pt();
1293       d1 -= part.Pt();
1294     }
1295   }
1296
1297   for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1298     Bool_t track2Found = kFALSE;
1299     Int_t index2 = jet2->TrackAt(iTrack2);
1300
1301     // now look for common particles in the track array
1302     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1303       AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1304       if (!track) {
1305         AliWarning(Form("Could not find track %d!", iTrack));
1306         continue;
1307       }
1308       Int_t MClabel = TMath::Abs(track->GetLabel());
1309       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1310       if (MClabel <= 0) continue;
1311
1312       Int_t index = -1;
1313       index = tracks2->GetIndexFromLabel(MClabel);
1314       if (index < 0) {
1315         AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1316         continue;
1317       }
1318
1319       if (index2 != index) continue;
1320
1321       // found common particle
1322       d1 -= track->Pt();
1323
1324       if (!track2Found) {
1325         AliVParticle *MCpart = tracks2->GetParticle(index2);
1326         AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1327                         iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1328         d2 -= MCpart->Pt();
1329       }
1330
1331       track2Found = kTRUE;
1332     }
1333
1334     // now look for common particles in the cluster array
1335     if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1336       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1337         AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1338         if (!clus) {
1339           AliWarning(Form("Could not find cluster %d!", iClus));
1340           continue;
1341         }
1342         TLorentzVector part;
1343         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1344
1345         for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1346           Int_t cellId = clus->GetCellAbsId(iCell);
1347           Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1348
1349           Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1350           if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1351           if (MClabel <= 0) continue;
1352
1353           Int_t index1 = -1;
1354           index1 = tracks2->GetIndexFromLabel(MClabel);
1355           if (index1 < 0) {
1356             AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1357             continue;
1358           }
1359
1360           if (index2 != index1) continue;
1361
1362           // found common particle
1363           d1 -= part.Pt() * cellFrac;
1364
1365           if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1366             AliVParticle *MCpart = tracks2->GetParticle(index2);
1367             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)!",
1368                             iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1369             d2 -= MCpart->Pt() * cellFrac;
1370           }
1371
1372           track2Found = kTRUE;
1373         }
1374       }
1375     }
1376     else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
1377       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1378         AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1379         if (!clus) {
1380           AliWarning(Form("Could not find cluster %d!", iClus));
1381           continue;
1382         }
1383         TLorentzVector part;
1384         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1385
1386         Int_t MClabel = TMath::Abs(clus->GetLabel());
1387         if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1388         if (MClabel <= 0) continue;
1389
1390         Int_t index = -1;
1391         index = tracks2->GetIndexFromLabel(MClabel);
1392
1393         if (index < 0) {
1394           AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1395           continue;
1396         }
1397
1398         if (index2 != index) continue;
1399
1400         // found common particle
1401         d1 -= part.Pt();
1402
1403         if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1404           AliVParticle *MCpart = tracks2->GetParticle(index2);
1405           AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1406                           iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1407
1408           d2 -= MCpart->Pt();
1409         }
1410
1411         track2Found = kTRUE;
1412       }
1413     }
1414   }
1415
1416   if (d1 < 0)
1417     d1 = 0;
1418
1419   if (d2 < 0)
1420     d2 = 0;
1421
1422   if (totalPt1 < 1)
1423     d1 = -1;
1424   else
1425     d1 /= totalPt1;
1426
1427   if (jet2->Pt() < 1)
1428     d2 = -1;
1429   else
1430     d2 /= jet2->Pt();
1431 }
1432
1433 //________________________________________________________________________
1434 void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1435 {
1436   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1437   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1438
1439   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1440
1441   AliParticleContainer *tracks1   = jets1->GetParticleContainer();
1442   AliClusterContainer  *clusters1 = jets1->GetClusterContainer();
1443   AliParticleContainer *tracks2   = jets2->GetParticleContainer();
1444   AliClusterContainer  *clusters2 = jets2->GetClusterContainer();
1445
1446   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1447   d1 = jet1->Pt();
1448   d2 = jet2->Pt();
1449
1450   if (tracks1 && tracks2) {
1451
1452     for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1453       Int_t index2 = jet2->TrackAt(iTrack2);
1454       for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
1455         Int_t index1 = jet1->TrackAt(iTrack1);
1456         if (index2 == index1) { // found common particle
1457           AliVParticle *part1 = tracks1->GetParticle(index1);
1458           if (!part1) {
1459             AliWarning(Form("Could not find track %d!", index1));
1460             continue;
1461           }
1462           AliVParticle *part2 = tracks2->GetParticle(index2);
1463           if (!part2) {
1464             AliWarning(Form("Could not find track %d!", index2));
1465             continue;
1466           }
1467
1468           d1 -= part1->Pt();
1469           d2 -= part2->Pt();
1470           break;
1471         }
1472       }
1473     }
1474
1475   }
1476
1477   if (clusters1 && clusters2) {
1478
1479     if (fUseCellsToMatch) {
1480       const Int_t nClus1 = jet1->GetNumberOfClusters();
1481
1482       Int_t ncells1[nClus1];
1483       UShort_t *cellsId1[nClus1];
1484       Double_t *cellsFrac1[nClus1];
1485       Int_t *sortedIndexes1[nClus1];
1486       Double_t ptClus1[nClus1];
1487       for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1488         Int_t index1 = jet1->ClusterAt(iClus1);
1489         AliVCluster *clus1 = clusters1->GetCluster(index1);
1490         if (!clus1) {
1491           AliWarning(Form("Could not find cluster %d!", index1));
1492           ncells1[iClus1] = 0;
1493           cellsId1[iClus1] = 0;
1494           cellsFrac1[iClus1] = 0;
1495           sortedIndexes1[iClus1] = 0;
1496           ptClus1[iClus1] = 0;
1497           continue;
1498         }
1499         TLorentzVector part1;
1500         clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1501
1502         ncells1[iClus1] = clus1->GetNCells();
1503         cellsId1[iClus1] = clus1->GetCellsAbsId();
1504         cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
1505         sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
1506         ptClus1[iClus1] = part1.Pt();
1507
1508         TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
1509       }
1510
1511       const Int_t nClus2 = jet2->GetNumberOfClusters();
1512
1513       const Int_t maxNcells2 = 11520;
1514       Int_t sortedIndexes2[maxNcells2];
1515       for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
1516         Int_t index2 = jet2->ClusterAt(iClus2);
1517         AliVCluster *clus2 =  clusters2->GetCluster(index2);
1518         if (!clus2) {
1519           AliWarning(Form("Could not find cluster %d!", index2));
1520           continue;
1521         }
1522         Int_t ncells2 = clus2->GetNCells();
1523         if (ncells2 >= maxNcells2) {
1524           AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
1525           continue;
1526         }
1527         UShort_t *cellsId2 = clus2->GetCellsAbsId();
1528         Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
1529
1530         TLorentzVector part2;
1531         clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1532         Double_t ptClus2 = part2.Pt();
1533
1534         TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
1535
1536         for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1537           if (sortedIndexes1[iClus1] == 0)
1538             continue;
1539           Int_t iCell1 = 0, iCell2 = 0;
1540           while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
1541             if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
1542               d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
1543               d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2;
1544               iCell1++;
1545               iCell2++;
1546             }
1547             else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) {
1548               iCell2++;
1549             }
1550             else {
1551               iCell1++;
1552             }
1553           }
1554         }
1555       }
1556     }
1557     else {
1558       for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1559         Int_t index2 = jet2->ClusterAt(iClus2);
1560         for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
1561           Int_t index1 = jet1->ClusterAt(iClus1);
1562           if (index2 == index1) { // found common particle
1563             AliVCluster *clus1 = clusters1->GetCluster(index1);
1564             if (!clus1) {
1565               AliWarning(Form("Could not find cluster %d!", index1));
1566               continue;
1567             }
1568             AliVCluster *clus2 =  clusters2->GetCluster(index2);
1569             if (!clus2) {
1570               AliWarning(Form("Could not find cluster %d!", index2));
1571               continue;
1572             }
1573             TLorentzVector part1, part2;
1574             clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1575             clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1576
1577             d1 -= part1.Pt();
1578             d2 -= part2.Pt();
1579             break;
1580           }
1581         }
1582       }
1583     }
1584   }
1585
1586   if (d1 < 0)
1587     d1 = 0;
1588
1589   if (d2 < 0)
1590     d2 = 0;
1591
1592   if (jet1->Pt() > 0)
1593     d1 /= jet1->Pt();
1594   else
1595     d1 = -1;
1596
1597   if (jet2->Pt() > 0)
1598     d2 /= jet2->Pt();
1599   else
1600     d2 = -1;
1601 }
1602
1603 //________________________________________________________________________
1604 void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
1605 {
1606   Double_t d1 = -1;
1607   Double_t d2 = -1;
1608
1609   switch (matching) {
1610   case kGeometrical:
1611     GetGeometricalMatchingLevel(jet1,jet2,d1);
1612     d2 = d1;
1613     break;
1614   case kMCLabel: // jet1 = detector level and jet2 = particle level!
1615     GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1616     break;
1617   case kSameCollections:
1618     GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
1619     break;
1620   default:
1621     ;
1622   }
1623
1624   if (d1 >= 0) {
1625
1626     if (d1 < jet1->ClosestJetDistance()) {
1627       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1628       jet1->SetClosestJet(jet2, d1);
1629     }
1630     else if (d1 < jet1->SecondClosestJetDistance()) {
1631       jet1->SetSecondClosestJet(jet2, d1);
1632     }
1633   }
1634
1635   if (d2 >= 0) {
1636
1637     if (d2 < jet2->ClosestJetDistance()) {
1638       jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1639       jet2->SetClosestJet(jet1, d2);
1640     }
1641     else if (d2 < jet2->SecondClosestJetDistance()) {
1642       jet2->SetSecondClosestJet(jet1, d2);
1643     }
1644   }
1645 }
1646
1647 //________________________________________________________________________
1648 Bool_t AliJetResponseMaker::FillHistograms()
1649 {
1650   // Fill histograms.
1651
1652   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1653   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1654
1655   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1656
1657   AliEmcalJet* jet1 = 0;
1658   AliEmcalJet* jet2 = 0;
1659
1660   jets2->ResetCurrentID();
1661   while ((jet2 = jets2->GetNextJet())) {
1662
1663     AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
1664
1665     Double_t ptLeading2 = jets2->GetLeadingHadronPt(jet2);
1666     Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
1667
1668     if (jets2->AcceptJet(jet2))
1669       FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), ptLeading2,
1670                    corrpt2, jet2->MCPt(), 2);
1671
1672     jet1 = jet2->MatchedJet();
1673
1674     if (!jet1) continue;
1675     if (!jets1->AcceptJet(jet1)) continue;
1676     if (jet1->MCPt() < fMinJetMCPt) continue;
1677
1678     Double_t d=-1, ce1=-1, ce2=-1;
1679     if (jet2->GetMatchingType() == kGeometrical) {
1680       if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
1681         GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
1682       else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
1683         GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
1684
1685       d = jet2->ClosestJetDistance();
1686     }
1687     else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1688       GetGeometricalMatchingLevel(jet1, jet2, d);
1689
1690       ce1 = jet1->ClosestJetDistance();
1691       ce2 = jet2->ClosestJetDistance();
1692     }
1693
1694     Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
1695     Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1696
1697     FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(),
1698                        jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(),
1699                        jet1->NEF(), jet2->NEF(), ptLeading1, ptLeading2);
1700   }
1701
1702   jets1->ResetCurrentID();
1703   while ((jet1 = jets1->GetNextAcceptJet())) {
1704     if (jet1->MCPt() < fMinJetMCPt) continue;
1705     AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
1706
1707     Double_t ptLeading1 = jets1->GetLeadingHadronPt(jet1);
1708     Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1709
1710     FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), ptLeading1,
1711                  corrpt1, jet1->MCPt(), 1);
1712   }
1713   return kTRUE;
1714 }