PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx

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