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