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 <TChain.h>
  10 #include <TClonesArray.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TList.h>
  14 #include <TLorentzVector.h>
  15
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliFJWrapper.h"
  19 #include "AliVCluster.h"
  20 #include "AliVTrack.h"
  21 #include "AliEmcalJet.h"
  22 #include "AliGenPythiaEventHeader.h"
  23 #include "AliMCEvent.h"
  24
  25 ClassImp(AliJetResponseMaker)
  26
  27 //________________________________________________________________________
  28 AliJetResponseMaker::AliJetResponseMaker() :
  29   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  30   fMCTracksName("MCParticles"),
  31   fMCJetsName("MCJets"),
  32   fMaxDistance(0.25),
  33   fDoWeighting(kFALSE),
  34   fEventWeightHist(kFALSE),
  35   fMCFiducial(kFALSE),
  36   fMCminEta(0),
  37   fMCmaxEta(0),
  38   fMCminPhi(0),
  39   fMCmaxPhi(0),
  40   fSelectPtHardBin(-999),
  41   fDoMatching(kTRUE),
  42   fPythiaHeader(0),
  43   fEventWeight(0),
  44   fPtHardBin(0),
  45   fNTrials(0),
  46   fMCTracks(0),
  47   fMCJets(0),
  48   fHistNTrials(0),
  49   fHistEvents(0),
  50   fHistMCJetPhiEta(0),
  51   fHistMCJetsPt(0),
  52   fHistMCJetPhiEtaFiducial(0),
  53   fHistMCJetsPtFiducial(0),
  54   fHistMCJetsNEFvsPt(0),
  55   fHistMCJetsZvsPt(0),
  56   fHistJetPhiEta(0),
  57   fHistJetsPt(0),
  58   fHistJetsNEFvsPt(0),
  59   fHistJetsZvsPt(0),
  60   fHistClosestDistance(0),
  61   fHistClosestDeltaPhi(0),
  62   fHistClosestDeltaEta(0),
  63   fHistClosestDeltaPt(0),
  64   fHistNonMatchedMCJetPt(0),
  65   fHistNonMatchedJetPt(0),
  66   fHistPartvsDetecPt(0),
  67   fHistMissedMCJets(0)
  68 {
  69   // Default constructor.
  70
  71   for (Int_t i = 0; i < 11; i++) {
  72     fHistEventWeight[i] = 0;
  73   }
  74
  75   SetMakeGeneralHistograms(kTRUE);
  76 }
  77
  78 //________________________________________________________________________
  79 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
  80   AliAnalysisTaskEmcalJet(name, kTRUE),
  81   fMCTracksName("MCParticles"),
  82   fMCJetsName("MCJets"),
  83   fMaxDistance(0.25),
  84   fDoWeighting(kFALSE),
  85   fEventWeightHist(kFALSE),
  86   fMCFiducial(kFALSE),
  87   fMCminEta(0),
  88   fMCmaxEta(0),
  89   fMCminPhi(0),
  90   fMCmaxPhi(0),
  91   fSelectPtHardBin(-999),
  92   fDoMatching(kTRUE),
  93   fPythiaHeader(0),
  94   fEventWeight(0),
  95   fPtHardBin(0),
  96   fNTrials(0),
  97   fMCTracks(0),
  98   fMCJets(0),
  99   fHistNTrials(0),
 100   fHistEvents(0),
 101   fHistMCJetPhiEta(0),
 102   fHistMCJetsPt(0),
 103   fHistMCJetPhiEtaFiducial(0),
 104   fHistMCJetsPtFiducial(0),
 105   fHistMCJetsNEFvsPt(0),
 106   fHistMCJetsZvsPt(0),
 107   fHistJetPhiEta(0),
 108   fHistJetsPt(0),
 109   fHistJetsNEFvsPt(0),
 110   fHistJetsZvsPt(0),
 111   fHistClosestDistance(0),
 112   fHistClosestDeltaPhi(0),
 113   fHistClosestDeltaEta(0),
 114   fHistClosestDeltaPt(0),
 115   fHistNonMatchedMCJetPt(0),
 116   fHistNonMatchedJetPt(0),
 117   fHistPartvsDetecPt(0),
 118   fHistMissedMCJets(0)
 119 {
 120   // Standard constructor.
 121
 122   for (Int_t i = 0; i < 11; i++) {
 123     fHistEventWeight[i] = 0;
 124   }
 125
 126   SetMakeGeneralHistograms(kTRUE);
 127 }
 128
 129 //________________________________________________________________________
 130 AliJetResponseMaker::~AliJetResponseMaker()
 131 {
 132   // Destructor
 133 }
 134
 135 //________________________________________________________________________
 136 void AliJetResponseMaker::UserCreateOutputObjects()
 137 {
 138   // Create user objects.
 139
 140   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 141
 142   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 143   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 144
 145   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
 146   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 147   fHistNTrials->GetYaxis()->SetTitle("trials");
 148   fOutput->Add(fHistNTrials);
 149
 150   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
 151   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 152   fHistEvents->GetYaxis()->SetTitle("total events");
 153   fOutput->Add(fHistEvents);
 154
 155   if (fEventWeightHist) {
 156     for (Int_t i = 0; i < 11; i++) {
 157       TString name(Form("fHistEventWeight_%d", i+1));
 158       fHistEventWeight[i] = new TH1F(name, name, 10, 0, 10);
 159       fOutput->Add(fHistEventWeight[i]);
 160     }
 161   }
 162
 163   for (Int_t i = 1; i < 12; i++) {
 164     fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 165     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 166   }
 167
 168   fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
 169   fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
 170   fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
 171   fOutput->Add(fHistJetPhiEta);
 172
 173   fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
 174   fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 175   fHistJetsPt->GetYaxis()->SetTitle("counts");
 176   fOutput->Add(fHistJetsPt);
 177
 178   fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 179   fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
 180   fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 181   fOutput->Add(fHistJetsZvsPt);
 182
 183   fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 184   fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
 185   fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 186   fOutput->Add(fHistJetsNEFvsPt);
 187
 188   fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
 189   fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
 190   fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
 191   fOutput->Add(fHistMCJetPhiEta);
 192
 193   fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
 194   fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 195   fHistMCJetsPt->GetYaxis()->SetTitle("counts");
 196   fOutput->Add(fHistMCJetsPt);
 197
 198   fHistMCJetPhiEtaFiducial = new TH2F("fHistMCJetPhiEtaFiducial", "fHistMCJetPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4);
 199   fHistMCJetPhiEtaFiducial->GetXaxis()->SetTitle("#eta");
 200   fHistMCJetPhiEtaFiducial->GetYaxis()->SetTitle("#phi");
 201   fOutput->Add(fHistMCJetPhiEtaFiducial);
 202
 203   fHistMCJetsPtFiducial = new TH1F("fHistMCJetsPtFiducial", "fHistMCJetsPtFiducial", fNbins, fMinBinPt, fMaxBinPt);
 204   fHistMCJetsPtFiducial->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 205   fHistMCJetsPtFiducial->GetYaxis()->SetTitle("counts");
 206   fOutput->Add(fHistMCJetsPtFiducial);
 207
 208   fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 209   fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
 210   fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 211   fOutput->Add(fHistMCJetsZvsPt);
 212
 213   fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 214   fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
 215   fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 216   fOutput->Add(fHistMCJetsNEFvsPt);
 217
 218   fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
 219   fHistClosestDistance->GetXaxis()->SetTitle("d");
 220   fHistClosestDistance->GetYaxis()->SetTitle("counts");
 221   fOutput->Add(fHistClosestDistance);
 222
 223   fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
 224   fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
 225   fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
 226   fOutput->Add(fHistClosestDeltaPhi);
 227
 228   fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2);
 229   fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
 230   fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
 231   fOutput->Add(fHistClosestDeltaEta);
 232
 233   fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
 234   fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
 235   fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
 236   fOutput->Add(fHistClosestDeltaPt);
 237
 238   fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
 239   fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 240   fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
 241   fOutput->Add(fHistNonMatchedMCJetPt);
 242
 243   fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
 244   fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 245   fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
 246   fOutput->Add(fHistNonMatchedJetPt);
 247
 248   fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 249   fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
 250   fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
 251   fOutput->Add(fHistPartvsDetecPt);
 252
 253   fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
 254   fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 255   fHistMissedMCJets->GetYaxis()->SetTitle("counts");
 256   fOutput->Add(fHistMissedMCJets);
 257
 258   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 259 }
 260
 261 //________________________________________________________________________
 262 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
 263 {
 264   // Return true if jet is accepted.
 265
 266   if (jet->Pt() <= fJetPtCut)
 267     return kFALSE;
 268   if (jet->Area() <= fJetAreaCut)
 269     return kFALSE;
 270
 271   return kTRUE;
 272 }
 273
 274 //________________________________________________________________________
 275 void AliJetResponseMaker::ExecOnce()
 276 {
 277   // Execute once.
 278
 279   if (!fMCJetsName.IsNull() && !fMCJets) {
 280     fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
 281     if (!fMCJets) {
 282       AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
 283       return;
 284     }
 285     else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 286       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data()));
 287       fMCJets = 0;
 288       return;
 289     }
 290   }
 291
 292   if (!fMCTracksName.IsNull() && !fMCTracks) {
 293     fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
 294     if (!fMCTracks) {
 295       AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data()));
 296       return;
 297     }
 298     else {
 299       TClass *cl = fMCTracks->GetClass();
 300       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
 301         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data()));
 302         fMCTracks = 0;
 303         return;
 304       }
 305     }
 306   }
 307
 308   AliAnalysisTaskEmcalJet::ExecOnce();
 309
 310   if (fMCFiducial) {
 311     fMCminEta = fJetMinEta;
 312     fMCmaxEta = fJetMaxEta;
 313     fMCminPhi = fJetMinPhi;
 314     fMCmaxPhi = fJetMaxPhi;
 315   }
 316   else {
 317     fMCminEta = fJetMinEta - fJetRadius;
 318     fMCmaxEta = fJetMaxEta + fJetRadius;
 319     fMCminPhi = fJetMinPhi - fJetRadius;
 320     fMCmaxPhi = fJetMaxPhi + fJetRadius;
 321   }
 322 }
 323
 324 //________________________________________________________________________
 325 Bool_t AliJetResponseMaker::IsEventSelected()
 326 {
 327   // Check if event is selected
 328
 329   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)
 330     return kFALSE;
 331
 332   return AliAnalysisTaskEmcalJet::IsEventSelected();
 333 }
 334
 335 //________________________________________________________________________
 336 Bool_t AliJetResponseMaker::RetrieveEventObjects()
 337 {
 338   // Retrieve event objects.
 339
 340   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 341     return kFALSE;
 342
 343   fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 344
 345   if (!fPythiaHeader)
 346     return kFALSE;
 347
 348   if (fDoWeighting)
 349     fEventWeight = fPythiaHeader->EventWeight();
 350   else
 351     fEventWeight = 1;
 352
 353   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 354   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 355
 356   Double_t pthard = fPythiaHeader->GetPtHard();
 357
 358   for (fPtHardBin = 0; fPtHardBin <= 11; fPtHardBin++) {
 359     if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
 360       break;
 361   }
 362
 363   fNTrials = fPythiaHeader->Trials();
 364
 365   return kTRUE;
 366 }
 367
 368 //________________________________________________________________________
 369 Bool_t AliJetResponseMaker::Run()
 370 {
 371   // Find the closest jets
 372
 373   if (!fDoMatching)
 374     return kTRUE;
 375
 376   DoJetLoop(fJets, fMCJets, kFALSE);
 377   DoJetLoop(fMCJets, fJets, kTRUE);
 378
 379   const Int_t nMCJets = fMCJets->GetEntriesFast();
 380
 381   for (Int_t i = 0; i < nMCJets; i++) {
 382
 383     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
 384
 385     if (!jet) {
 386       AliError(Form("Could not receive jet %d", i));
 387       continue;
 388     }
 389
 390     if (!AcceptJet(jet))
 391       continue;
 392
 393     if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi)
 394       continue;
 395
 396     if (jet->Pt() > fMaxBinPt)
 397       continue;
 398
 399     if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
 400         jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
 401       jet->SetMatchedToClosest();
 402       jet->ClosestJet()->SetMatchedToClosest();
 403     }
 404   }
 405
 406   return kTRUE;
 407 }
 408
 409 //________________________________________________________________________
 410 void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc)
 411 {
 412   // Do the jet loop.
 413
 414   Int_t nJets1 = jets1->GetEntriesFast();
 415   Int_t nJets2 = jets2->GetEntriesFast();
 416
 417   for (Int_t i = 0; i < nJets1; i++) {
 418
 419     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
 420
 421     if (!jet1) {
 422       AliError(Form("Could not receive jet %d", i));
 423       continue;
 424     }
 425
 426     if (!AcceptJet(jet1))
 427       continue;
 428
 429     if (!mc) {
 430       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
 431         continue;
 432     }
 433     else {
 434       if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi)
 435         continue;
 436     }
 437
 438     for (Int_t j = 0; j < nJets2; j++) {
 439
 440       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
 441
 442       if (!jet2) {
 443         AliError(Form("Could not receive jet %d", j));
 444         continue;
 445       }
 446
 447       if (!AcceptJet(jet2))
 448         continue;
 449
 450       if (mc) {
 451         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
 452           continue;
 453       }
 454       else {
 455         if (jet1->Eta() < fMCminEta || jet1->Eta() > fMCmaxEta || jet1->Phi() < fMCminPhi || jet1->Phi() > fMCmaxPhi)
 456           continue;
 457       }
 458
 459       Double_t deta = jet2->Eta() - jet1->Eta();
 460       Double_t dphi = jet2->Phi() - jet1->Phi();
 461       Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
 462
 463       if (d < jet1->ClosestJetDistance()) {
 464         jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
 465         jet1->SetClosestJet(jet2, d);
 466       }
 467       else if (d < jet1->SecondClosestJetDistance()) {
 468         jet1->SetSecondClosestJet(jet2, d);
 469       }
 470     }
 471   }
 472 }
 473
 474 //________________________________________________________________________
 475 Bool_t AliJetResponseMaker::FillHistograms()
 476 {
 477   // Fill histograms.
 478
 479   fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
 480   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
 481   if (fEventWeightHist)
 482     fHistEventWeight[fPtHardBin]->Fill(fPythiaHeader->EventWeight());
 483
 484   const Int_t nMCJets = fMCJets->GetEntriesFast();
 485
 486   for (Int_t i = 0; i < nMCJets; i++) {
 487
 488     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
 489
 490     if (!jet) {
 491       AliError(Form("Could not receive jet %d", i));
 492       continue;
 493     }
 494
 495     if (!AcceptJet(jet))
 496       continue;
 497
 498     if (jet->Eta() < fMCminEta || jet->Eta() > fMCmaxEta || jet->Phi() < fMCminPhi || jet->Phi() > fMCmaxPhi)
 499       continue;
 500
 501     if (jet->Pt() > fMaxBinPt)
 502       continue;
 503
 504     if (jet->MatchedJet()) {
 505
 506       if (!AcceptBiasJet(jet->MatchedJet()) ||
 507           jet->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
 508           jet->MatchedJet()->Pt() > fMaxBinPt) {
 509         fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
 510       }
 511       else {
 512         fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
 513
 514         Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
 515         fHistClosestDeltaEta->Fill(deta, fEventWeight);
 516
 517         Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
 518         fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
 519
 520         Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
 521         fHistClosestDeltaPt->Fill(dpt, fEventWeight);
 522
 523         fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
 524       }
 525     }
 526     else {
 527       fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
 528       fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
 529     }
 530
 531     fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
 532     fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 533
 534     fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
 535
 536     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 537       AliVParticle *track = jet->TrackAt(it, fMCTracks);
 538       if (track)
 539         fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
 540     }
 541
 542     if (!AcceptBiasJet(jet))
 543       continue;
 544     if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
 545       continue;
 546
 547     fHistMCJetsPtFiducial->Fill(jet->Pt(), fEventWeight);
 548     fHistMCJetPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 549   }
 550
 551   const Int_t nJets = fJets->GetEntriesFast();
 552
 553   for (Int_t i = 0; i < nJets; i++) {
 554
 555     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
 556
 557     if (!jet) {
 558       AliError(Form("Could not receive mc jet %d", i));
 559       continue;
 560     }
 561
 562     if (!AcceptJet(jet))
 563       continue;
 564     if (!AcceptBiasJet(jet))
 565       continue;
 566     if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
 567       continue;
 568     if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi)
 569       continue;
 570
 571     if (!jet->MatchedJet()) {
 572       fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
 573     }
 574
 575     fHistJetsPt->Fill(jet->Pt(), fEventWeight);
 576
 577     fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 578
 579     fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
 580
 581     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 582       AliVParticle *track = jet->TrackAt(it, fTracks);
 583       if (track)
 584         fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
 585     }
 586
 587     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 588       AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 589       if (cluster) {
 590         TLorentzVector nP;
 591         cluster->GetMomentum(nP, fVertex);
 592         fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight);
 593       }
 594     }
 595   }
 596
 597   return kTRUE;
 598 }