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