move EMCALJetTasks from PWGGA to PWGJE
[u/mrichter/AliRoot.git] / 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 }