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 <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TH3F.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliVCluster.h"
  16 #include "AliVTrack.h"
  17 #include "AliEmcalJet.h"
  18 #include "AliGenPythiaEventHeader.h"
  19 #include "AliLog.h"
  20 #include "AliMCEvent.h"
  21 #include "AliRhoParameter.h"
  22
  23 ClassImp(AliJetResponseMaker)
  24
  25 //________________________________________________________________________
  26 AliJetResponseMaker::AliJetResponseMaker() :
  27   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  28   fTracks2Name(""),
  29   fCalo2Name(""),
  30   fJets2Name(""),
  31   fRho2Name(""),
  32   fPtBiasJet2Track(0),
  33   fPtBiasJet2Clus(0),
  34   fAreCollections1MC(kFALSE),
  35   fAreCollections2MC(kTRUE),
  36   fMatching(kNoMatching),
  37   fMatchingPar1(0),
  38   fMatchingPar2(0),
  39   fJet2MinEta(-999),
  40   fJet2MaxEta(-999),
  41   fJet2MinPhi(-999),
  42   fJet2MaxPhi(-999),
  43   fSelectPtHardBin(-999),
  44   fPythiaHeader(0),
  45   fPtHardBin(0),
  46   fNTrials(0),
  47   fTracks2(0),
  48   fCaloClusters2(0),
  49   fJets2(0),
  50   fRho2(0),
  51   fRho2Val(0),
  52   fTracks2Map(0),
  53   fHistNTrials(0),
  54   fHistEvents(0),
  55   fHistJets1PhiEta(0),
  56   fHistJets1PtArea(0),
  57   fHistJets1CorrPtArea(0),
  58   fHistJets2PhiEta(0),
  59   fHistJets2PtArea(0),
  60   fHistJets2CorrPtArea(0),
  61   fHistJets2PhiEtaAcceptance(0),
  62   fHistJets2PtAreaAcceptance(0),
  63   fHistJets2CorrPtAreaAcceptance(0),
  64   fHistMatchingLevelvsJet2Pt(0),
  65   fHistDistancevsCommonEnergy(0),
  66   fHistDeltaEtaPhivsJet2Pt(0),
  67   fHistDeltaPtvsJet2Pt(0),
  68   fHistDeltaPtvsMatchingLevel(0),
  69   fHistDeltaCorrPtvsJet2Pt(0),
  70   fHistDeltaCorrPtvsMatchingLevel(0),
  71   fHistNonMatchedJets1PtArea(0),
  72   fHistNonMatchedJets2PtArea(0),
  73   fHistNonMatchedJets1CorrPtArea(0),
  74   fHistNonMatchedJets2CorrPtArea(0),
  75   fHistJet1PtvsJet2Pt(0),
  76   fHistJet1CorrPtvsJet2CorrPt(0),
  77   fHistMissedJets2PtArea(0)
  78 {
  79   // Default constructor.
  80
  81   SetMakeGeneralHistograms(kTRUE);
  82 }
  83
  84 //________________________________________________________________________
  85 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
  86   AliAnalysisTaskEmcalJet(name, kTRUE),
  87   fTracks2Name("MCParticles"),
  88   fCalo2Name(""),
  89   fJets2Name("MCJets"),
  90   fRho2Name(""),
  91   fPtBiasJet2Track(0),
  92   fPtBiasJet2Clus(0),
  93   fAreCollections1MC(kFALSE),
  94   fAreCollections2MC(kTRUE),
  95   fMatching(kNoMatching),
  96   fMatchingPar1(0),
  97   fMatchingPar2(0),
  98   fJet2MinEta(-999),
  99   fJet2MaxEta(-999),
 100   fJet2MinPhi(-999),
 101   fJet2MaxPhi(-999),
 102   fSelectPtHardBin(-999),
 103   fPythiaHeader(0),
 104   fPtHardBin(0),
 105   fNTrials(0),
 106   fTracks2(0),
 107   fCaloClusters2(0),
 108   fJets2(0),
 109   fRho2(0),
 110   fRho2Val(0),
 111   fTracks2Map(0),
 112   fHistNTrials(0),
 113   fHistEvents(0),
 114   fHistJets1PhiEta(0),
 115   fHistJets1PtArea(0),
 116   fHistJets1CorrPtArea(0),
 117   fHistJets2PhiEta(0),
 118   fHistJets2PtArea(0),
 119   fHistJets2CorrPtArea(0),
 120   fHistJets2PhiEtaAcceptance(0),
 121   fHistJets2PtAreaAcceptance(0),
 122   fHistJets2CorrPtAreaAcceptance(0),
 123   fHistMatchingLevelvsJet2Pt(0),
 124   fHistDistancevsCommonEnergy(0),
 125   fHistDeltaEtaPhivsJet2Pt(0),
 126   fHistDeltaPtvsJet2Pt(0),
 127   fHistDeltaPtvsMatchingLevel(0),
 128   fHistDeltaCorrPtvsJet2Pt(0),
 129   fHistDeltaCorrPtvsMatchingLevel(0),
 130   fHistNonMatchedJets1PtArea(0),
 131   fHistNonMatchedJets2PtArea(0),
 132   fHistNonMatchedJets1CorrPtArea(0),
 133   fHistNonMatchedJets2CorrPtArea(0),
 134   fHistJet1PtvsJet2Pt(0),
 135   fHistJet1CorrPtvsJet2CorrPt(0),
 136   fHistMissedJets2PtArea(0)
 137 {
 138   // Standard constructor.
 139
 140   SetMakeGeneralHistograms(kTRUE);
 141 }
 142
 143 //________________________________________________________________________
 144 AliJetResponseMaker::~AliJetResponseMaker()
 145 {
 146   // Destructor
 147 }
 148
 149 //________________________________________________________________________
 150 void AliJetResponseMaker::UserCreateOutputObjects()
 151 {
 152   // Create user objects.
 153
 154   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 155
 156   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 157   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 158
 159   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
 160   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 161   fHistNTrials->GetYaxis()->SetTitle("trials");
 162   fOutput->Add(fHistNTrials);
 163
 164   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
 165   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 166   fHistEvents->GetYaxis()->SetTitle("total events");
 167   fOutput->Add(fHistEvents);
 168
 169   for (Int_t i = 1; i < 12; i++) {
 170     fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 171     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 172   }
 173
 174   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
 175   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
 176   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
 177   fOutput->Add(fHistJets1PhiEta);
 178
 179   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 180   fHistJets1PtArea->GetXaxis()->SetTitle("area");
 181   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 182   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
 183   fOutput->Add(fHistJets1PtArea);
 184
 185   if (!fRhoName.IsNull()) {
 186     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 187     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
 188     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
 189     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 190     fOutput->Add(fHistJets1CorrPtArea);
 191   }
 192
 193   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
 194   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
 195   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
 196   fOutput->Add(fHistJets2PhiEta);
 197
 198   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
 199   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
 200   fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
 201   fOutput->Add(fHistJets2PhiEtaAcceptance);
 202
 203   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 204   fHistJets2PtArea->GetXaxis()->SetTitle("area");
 205   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 206   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
 207   fOutput->Add(fHistJets2PtArea);
 208
 209   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 210   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
 211   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 212   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
 213   fOutput->Add(fHistJets2PtAreaAcceptance);
 214
 215   if (!fRho2Name.IsNull()) {
 216     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 217     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
 218     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 219     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 220     fOutput->Add(fHistJets2CorrPtArea);
 221
 222     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 223     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
 224     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
 225     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
 226     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
 227   }
 228
 229   fHistMatchingLevelvsJet2Pt = new TH2F("fHistMatchingLevelvsJet2Pt", "fHistMatchingLevelvsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
 230   fHistMatchingLevelvsJet2Pt->GetXaxis()->SetTitle("Matching level");
 231   fHistMatchingLevelvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 232   fHistMatchingLevelvsJet2Pt->GetZaxis()->SetTitle("counts");
 233   fOutput->Add(fHistMatchingLevelvsJet2Pt);
 234
 235   fHistDistancevsCommonEnergy = new TH2F("fHistDistancevsCommonEnergy", "fHistDistancevsCommonEnergy", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
 236   fHistDistancevsCommonEnergy->GetXaxis()->SetTitle("Distance");
 237   fHistDistancevsCommonEnergy->GetYaxis()->SetTitle("Common energy (%)");
 238   fHistDistancevsCommonEnergy->GetZaxis()->SetTitle("counts");
 239   fOutput->Add(fHistDistancevsCommonEnergy);
 240
 241   fHistDeltaEtaPhivsJet2Pt = new TH3F("fHistDeltaEtaPhivsJet2Pt", "fHistDeltaEtaPhivsJet2Pt", 40, -1, 1, 128, -1.6, 4.8, fNbins/2, fMinBinPt, fMaxBinPt);
 242   fHistDeltaEtaPhivsJet2Pt->GetXaxis()->SetTitle("#Delta#eta");
 243   fHistDeltaEtaPhivsJet2Pt->GetYaxis()->SetTitle("#Delta#phi");
 244   fHistDeltaEtaPhivsJet2Pt->GetZaxis()->SetTitle("p_{T,2}");
 245   fOutput->Add(fHistDeltaEtaPhivsJet2Pt);
 246
 247   fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 248   fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 249   fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
 250   fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 251   fOutput->Add(fHistDeltaPtvsJet2Pt);
 252
 253   fHistDeltaPtvsMatchingLevel = new TH2F("fHistDeltaPtvsMatchingLevel", "fHistDeltaPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 254   fHistDeltaPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
 255   fHistDeltaPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T} (GeV/c)");
 256   fHistDeltaPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
 257   fOutput->Add(fHistDeltaPtvsMatchingLevel);
 258
 259   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 260     fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 261     fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");
 262     fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
 263     fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
 264     fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
 265
 266     fHistDeltaCorrPtvsMatchingLevel = new TH2F("fHistDeltaCorrPtvsMatchingLevel", "fHistDeltaCorrPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 267     fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");
 268     fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
 269     fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
 270     fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
 271   }
 272
 273   fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 274   fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
 275   fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 276   fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
 277   fOutput->Add(fHistNonMatchedJets1PtArea);
 278
 279   fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 280   fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
 281   fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 282   fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
 283   fOutput->Add(fHistNonMatchedJets2PtArea);
 284
 285   if (!fRhoName.IsNull()) {
 286     fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 287     fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
 288     fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
 289     fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
 290     fOutput->Add(fHistNonMatchedJets1CorrPtArea);
 291   }
 292
 293   if (!fRho2Name.IsNull()) {
 294     fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 295     fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
 296     fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 297     fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
 298     fOutput->Add(fHistNonMatchedJets2CorrPtArea);
 299   }
 300
 301   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 302   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
 303   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
 304   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
 305   fOutput->Add(fHistJet1PtvsJet2Pt);
 306
 307   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 308     if (fRhoName.IsNull())
 309       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 310     else if (fRho2Name.IsNull())
 311       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 312     else
 313       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
 314     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
 315     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
 316     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
 317     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
 318   }
 319
 320   fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 321   fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");
 322   fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
 323   fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
 324   fOutput->Add(fHistMissedJets2PtArea);
 325
 326   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 327 }
 328
 329 //________________________________________________________________________
 330 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
 331 {
 332   // Return true if jet is accepted.
 333
 334   if (jet->Pt() <= fJetPtCut)
 335     return kFALSE;
 336   if (jet->Area() <= fJetAreaCut)
 337     return kFALSE;
 338
 339   return kTRUE;
 340 }
 341
 342 //________________________________________________________________________
 343 Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
 344 {
 345   // Accept jet with a bias.
 346
 347   if (fLeadingHadronType == 0) {
 348     if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
 349   }
 350   else if (fLeadingHadronType == 1) {
 351     if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
 352   }
 353   else {
 354     if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
 355   }
 356
 357   return kTRUE;
 358 }
 359
 360 //________________________________________________________________________
 361 void AliJetResponseMaker::ExecOnce()
 362 {
 363   // Execute once.
 364
 365   if (!fJets2Name.IsNull() && !fJets2) {
 366     fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
 367     if (!fJets2) {
 368       AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
 369       return;
 370     }
 371     else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
 372       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data()));
 373       fJets2 = 0;
 374       return;
 375     }
 376   }
 377
 378   if (!fTracks2Name.IsNull() && !fTracks2) {
 379     fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
 380     if (!fTracks2) {
 381       AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data()));
 382       return;
 383     }
 384     else {
 385       TClass *cl = fTracks2->GetClass();
 386       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
 387         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data()));
 388         fTracks2 = 0;
 389         return;
 390       }
 391     }
 392
 393     if (fAreCollections2MC) {
 394       fTracks2Map = dynamic_cast<TH1*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
 395       // this is needed to map the MC labels with the indexes of the MC particle collection
 396       // if teh map is not given, the MC labels are assumed to be consistent with the indexes (which is not the case if AliEmcalMCTrackSelector is used)
 397       if (!fTracks2Map) {
 398         AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data()));
 399         fTracks2Map = new TH1I("tracksMap","tracksMap",9999,0,1);
 400         for (Int_t i = 0; i < 9999; i++) {
 401           fTracks2Map->SetBinContent(i,i);
 402         }
 403       }
 404     }
 405   }
 406
 407   if (!fCalo2Name.IsNull() && !fCaloClusters2) {
 408     fCaloClusters2 =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
 409     if (!fCaloClusters2) {
 410       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data()));
 411       return;
 412     } else {
 413       TClass *cl = fCaloClusters2->GetClass();
 414       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
 415         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data()));
 416         fCaloClusters2 = 0;
 417         return;
 418       }
 419     }
 420   }
 421
 422   if (!fRho2Name.IsNull() && !fRho2) {
 423     fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
 424     if (!fRho2) {
 425       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
 426       fInitialized = kFALSE;
 427       return;
 428     }
 429   }
 430
 431   if (fJet2MinEta == -999)
 432     fJet2MinEta = fJetMinEta - fJetRadius;
 433   if (fJet2MaxEta == -999)
 434     fJet2MaxEta = fJetMaxEta + fJetRadius;
 435   if (fJet2MinPhi == -999)
 436     fJet2MinPhi = fJetMinPhi - fJetRadius;
 437   if (fJet2MaxPhi == -999)
 438     fJet2MaxPhi = fJetMaxPhi + fJetRadius;
 439
 440   AliAnalysisTaskEmcalJet::ExecOnce();
 441 }
 442
 443 //________________________________________________________________________
 444 Bool_t AliJetResponseMaker::IsEventSelected()
 445 {
 446   // Check if event is selected
 447
 448   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
 449     return kFALSE;
 450
 451   return AliAnalysisTaskEmcalJet::IsEventSelected();
 452 }
 453
 454 //________________________________________________________________________
 455 Bool_t AliJetResponseMaker::RetrieveEventObjects()
 456 {
 457   // Retrieve event objects.
 458
 459   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 460     return kFALSE;
 461
 462   if (fRho2)
 463     fRho2Val = fRho2->GetVal();
 464
 465   fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 466
 467   if (!fPythiaHeader)
 468     return kFALSE;
 469
 470   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 471   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 472
 473   Double_t pthard = fPythiaHeader->GetPtHard();
 474
 475   for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 476     if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
 477       break;
 478   }
 479
 480   fNTrials = fPythiaHeader->Trials();
 481
 482   return kTRUE;
 483 }
 484
 485 //________________________________________________________________________
 486 Bool_t AliJetResponseMaker::Run()
 487 {
 488   // Find the closest jets
 489
 490   if (fMatching == kNoMatching)
 491     return kTRUE;
 492   else
 493     return DoJetMatching();
 494 }
 495
 496 //________________________________________________________________________
 497 Bool_t AliJetResponseMaker::DoJetMatching()
 498 {
 499   DoJetLoop(kFALSE);
 500
 501   const Int_t nJets = fJets->GetEntriesFast();
 502
 503   for (Int_t i = 0; i < nJets; i++) {
 504
 505     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
 506
 507     if (!jet1) {
 508       AliError(Form("Could not receive jet %d", i));
 509       continue;
 510     }
 511
 512     if (!AcceptJet(jet1))
 513       continue;
 514
 515     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 516       continue;
 517
 518     if (jet1->Pt() > fMaxBinPt)
 519       continue;
 520
 521     if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 &&
 522         jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) {    // Matched jet found
 523       jet1->SetMatchedToClosest(fMatching);
 524       jet1->ClosestJet()->SetMatchedToClosest(fMatching);
 525     }
 526   }
 527
 528   return kTRUE;
 529 }
 530
 531 //________________________________________________________________________
 532 void AliJetResponseMaker::DoJetLoop(Bool_t order)
 533 {
 534   // Do the jet loop.
 535
 536   TClonesArray *jets1 = 0;
 537   TClonesArray *jets2 = 0;
 538
 539   if (order) {
 540     jets1 = fJets2;
 541     jets2 = fJets;
 542   }
 543   else {
 544     jets1 = fJets;
 545     jets2 = fJets2;
 546   }
 547
 548   Int_t nJets1 = jets1->GetEntriesFast();
 549   Int_t nJets2 = jets2->GetEntriesFast();
 550
 551   for (Int_t i = 0; i < nJets1; i++) {
 552
 553     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
 554
 555     if (!jet1) {
 556       AliError(Form("Could not receive jet %d", i));
 557       continue;
 558     }
 559
 560     if (!AcceptJet(jet1))
 561       continue;
 562
 563     if (order) {
 564      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
 565         continue;
 566     }
 567     else {
 568       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 569         continue;
 570     }
 571
 572     for (Int_t j = 0; j < nJets2; j++) {
 573
 574       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
 575
 576       if (!jet2) {
 577         AliError(Form("Could not receive jet %d", j));
 578         continue;
 579       }
 580
 581       if (!AcceptJet(jet2))
 582         continue;
 583
 584       if (order) {
 585         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
 586           continue;
 587       }
 588       else {
 589         if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
 590           continue;
 591       }
 592
 593       SetMatchingLevel(jet1, jet2, fMatching);
 594     }
 595   }
 596 }
 597
 598 //________________________________________________________________________
 599 void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
 600 {
 601   Double_t deta = jet2->Eta() - jet1->Eta();
 602   Double_t dphi = jet2->Phi() - jet1->Phi();
 603   d = TMath::Sqrt(deta * deta + dphi * dphi);
 604 }
 605
 606 //________________________________________________________________________
 607 void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
 608 {
 609   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
 610   d1 = jet1->Pt();
 611   d2 = jet2->Pt();
 612   Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
 613
 614   for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
 615     Bool_t track2Found = kFALSE;
 616     Int_t index2 = jet2->TrackAt(iTrack2);
 617     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
 618       AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
 619       if (!track) {
 620         AliWarning(Form("Could not find track %d!", iTrack));
 621         continue;
 622       }
 623       Int_t MClabel = track->GetLabel();
 624       Int_t index = -1;
 625
 626       if (MClabel <= 0) {// this is not a MC particle; remove it completely
 627         AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
 628         totalPt1 -= track->Pt();
 629         d1 -= track->Pt();
 630         continue;
 631       }
 632       else if (MClabel < fTracks2Map->GetNbinsX()-2) {
 633         index = fTracks2Map->GetBinContent(MClabel);
 634       }
 635
 636       if (index < 0) {
 637         AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
 638         continue;
 639       }
 640
 641       if (index2 == index) { // found common particle
 642         track2Found = kTRUE;
 643         d1 -= track->Pt();
 644         AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
 645         AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
 646                         iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
 647         d2 -= MCpart->Pt();
 648         break;
 649       }
 650     }
 651     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
 652       AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
 653       if (!clus) {
 654         AliWarning(Form("Could not find cluster %d!", iClus));
 655         continue;
 656       }
 657       TLorentzVector part;
 658       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
 659
 660       if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
 661         for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
 662           Int_t cellId = clus->GetCellAbsId(iCell);
 663           Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
 664
 665           Int_t MClabel = fCaloCells->GetCellMCLabel(cellId);
 666           Int_t index = -1;
 667
 668           if (MClabel <= 0) {// this is not a MC particle; remove it completely
 669             AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
 670             totalPt1 -= part.Pt() * cellFrac;
 671             d1 -= part.Pt() * cellFrac;
 672             continue;
 673           }
 674           else if (MClabel < fTracks2Map->GetNbinsX()-2) {
 675             index = fTracks2Map->GetBinContent(MClabel);
 676           }
 677
 678           if (index < 0) {
 679             AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
 680             continue;
 681           }
 682           if (index2 == index) { // found common particle
 683             d1 -= part.Pt() * cellFrac;
 684
 685             if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
 686               AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
 687               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)!",
 688                               iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
 689               d2 -= MCpart->Pt() * cellFrac;
 690             }
 691             break;
 692           }
 693         }
 694       }
 695       else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
 696         Int_t MClabel = clus->GetLabel();
 697         Int_t index = -1;
 698
 699         if (MClabel <= 0) {// this is not a MC particle; remove it completely
 700           AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
 701           totalPt1 -= part.Pt();
 702           d1 -= part.Pt();
 703           continue;
 704         }
 705         else if (MClabel < fTracks2Map->GetNbinsX()-2) {
 706           index = fTracks2Map->GetBinContent(MClabel);
 707         }
 708
 709         if (index < 0) {
 710           AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
 711           continue;
 712         }
 713         if (index2 == index) { // found common particle
 714           d1 -= part.Pt();
 715
 716           if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
 717             AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
 718             AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
 719                             iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
 720
 721             d2 -= MCpart->Pt();
 722           }
 723           break;
 724         }
 725       }
 726     }
 727   }
 728   if (d1 <= 0 || totalPt1 < 1)
 729     d1 = 0;
 730   else
 731     d1 /= totalPt1;
 732
 733   if (jet2->Pt() > 0 && d2 > 0)
 734     d2 /= jet2->Pt();
 735   else
 736     d2 = 0;
 737 }
 738
 739 //________________________________________________________________________
 740 void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
 741 {
 742   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
 743   d1 = jet1->Pt();
 744   d2 = jet2->Pt();
 745
 746   if (fTracks && fTracks2) {
 747
 748     for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
 749       Int_t index2 = jet2->TrackAt(iTrack2);
 750       for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
 751         Int_t index = jet1->TrackAt(iTrack);
 752         if (index2 == index) { // found common particle
 753           AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
 754           if (!part) {
 755             AliWarning(Form("Could not find track %d!", index));
 756             continue;
 757           }
 758           AliVParticle *part2 = static_cast<AliVParticle*>(fTracks2->At(index2));
 759           if (!part2) {
 760             AliWarning(Form("Could not find track %d!", index2));
 761             continue;
 762           }
 763
 764           d1 -= part->Pt();
 765           d2 -= part2->Pt();
 766           break;
 767         }
 768       }
 769     }
 770
 771   }
 772
 773   if (fCaloClusters && fCaloClusters2) {
 774
 775     for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
 776       Int_t index2 = jet2->ClusterAt(iClus2);
 777       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
 778         Int_t index = jet1->ClusterAt(iClus);
 779         AliVCluster *clus =  static_cast<AliVCluster*>(fCaloClusters->At(index));
 780         if (!clus) {
 781           AliWarning(Form("Could not find cluster %d!", index));
 782           continue;
 783         }
 784         AliVCluster *clus2 =  static_cast<AliVCluster*>(fCaloClusters2->At(index2));
 785         if (!clus2) {
 786           AliWarning(Form("Could not find cluster %d!", index2));
 787           continue;
 788         }
 789         TLorentzVector part, part2;
 790         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
 791         clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
 792
 793         d1 -= part.Pt();
 794         d2 -= part2.Pt();
 795         break;
 796       }
 797     }
 798
 799   }
 800
 801   if (jet1->Pt() > 0 && d1 > 0)
 802     d1 /= jet1->Pt();
 803   else
 804     d1 = 0;
 805
 806   if (jet2->Pt() > 0 && d2 > 0)
 807     d2 /= jet2->Pt();
 808   else
 809     d2 = 0;
 810 }
 811
 812 //________________________________________________________________________
 813 void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching)
 814 {
 815   Double_t d1 = -1;
 816   Double_t d2 = -1;
 817
 818   switch (matching) {
 819   case kGeometrical:
 820     GetGeometricalMatchingLevel(jet1,jet2,d1);
 821     d2 = d1;
 822     break;
 823   case kMCLabel: // jet1 = detector level and jet2 = particle level!
 824     GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
 825     break;
 826   case kSameCollections:
 827     GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
 828     break;
 829   default:
 830     ;
 831   }
 832
 833   if (d1 > 0) {
 834
 835     if (d1 < jet1->ClosestJetDistance()) {
 836       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
 837       jet1->SetClosestJet(jet2, d1);
 838     }
 839     else if (d1 < jet1->SecondClosestJetDistance()) {
 840       jet1->SetSecondClosestJet(jet2, d1);
 841     }
 842   }
 843
 844   if (d2 > 0) {
 845
 846     if (d2 < jet2->ClosestJetDistance()) {
 847       jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
 848       jet2->SetClosestJet(jet1, d2);
 849     }
 850     else if (d2 < jet2->SecondClosestJetDistance()) {
 851       jet2->SetSecondClosestJet(jet1, d2);
 852     }
 853   }
 854 }
 855
 856 //________________________________________________________________________
 857 Bool_t AliJetResponseMaker::FillHistograms()
 858 {
 859   // Fill histograms.
 860
 861   fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
 862   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
 863
 864   const Int_t nJets2 = fJets2->GetEntriesFast();
 865
 866   for (Int_t i = 0; i < nJets2; i++) {
 867
 868     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
 869
 870     if (!jet2) {
 871       AliError(Form("Could not receive jet2 %d", i));
 872       continue;
 873     }
 874
 875     if (jet2->Pt() > fMaxBinPt)
 876       continue;
 877
 878     if (!AcceptJet(jet2))
 879       continue;
 880
 881     if (AcceptBiasJet(jet2) &&
 882         (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
 883
 884       fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
 885       fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
 886
 887       if (!fRho2Name.IsNull())
 888         fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
 889     }
 890
 891     if (!AcceptBiasJet2(jet2))
 892       continue;
 893
 894     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
 895       continue;
 896
 897     fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 898     fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
 899
 900     if (!fRho2Name.IsNull())
 901       fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
 902
 903     if (jet2->MatchedJet()) {
 904
 905       if (!AcceptBiasJet(jet2->MatchedJet()) ||
 906           jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
 907           jet2->MatchedJet()->Pt() > fMaxBinPt) {
 908         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 909       }
 910       else {
 911         Double_t d1=-1, d2=-1;
 912         if (jet2->GetMatchingType() == kGeometrical) {
 913
 914           if (fAreCollections2MC && !fAreCollections1MC)
 915             GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
 916           else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
 917             GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
 918
 919           fHistDistancevsCommonEnergy->Fill(jet2->ClosestJetDistance(), d2);
 920         }
 921         else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
 922           GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d1);
 923           fHistDistancevsCommonEnergy->Fill(d1, jet2->ClosestJetDistance());
 924         }
 925
 926         fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
 927
 928         Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
 929         Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
 930         fHistDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
 931
 932         Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
 933         fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
 934         fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
 935
 936         fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
 937
 938         if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
 939           dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
 940           fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
 941           fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
 942           fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
 943         }
 944       }
 945     }
 946     else {
 947       fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 948       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
 949
 950       if (!fRho2Name.IsNull())
 951         fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
 952     }
 953   }
 954
 955   const Int_t nJets1 = fJets->GetEntriesFast();
 956
 957   for (Int_t i = 0; i < nJets1; i++) {
 958
 959     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
 960
 961     if (!jet1) {
 962       AliError(Form("Could not receive jet %d", i));
 963       continue;
 964     }
 965
 966     if (!AcceptJet(jet1))
 967       continue;
 968
 969     if (!AcceptBiasJet(jet1))
 970       continue;
 971
 972     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
 973       continue;
 974
 975     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 976       continue;
 977
 978     if (!jet1->MatchedJet()) {
 979       fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
 980       if (!fRhoName.IsNull())
 981         fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
 982     }
 983
 984     fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
 985     fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
 986
 987     if (!fRhoName.IsNull())
 988       fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
 989   }
 990
 991   return kTRUE;
 992 }