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