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