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