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