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