up 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 <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 "AliMCEvent.h"
20 #include "AliLog.h"
21 #include "AliRhoParameter.h"
22
23 ClassImp(AliJetResponseMaker)
24
25 //________________________________________________________________________
26 AliJetResponseMaker::AliJetResponseMaker() : 
27   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
28   fTracks2Name(""),
29   fCalo2Name(""),
30   fJets2Name(""),
31   fRho2Name(""),
32   fPtBiasJet2Track(0),
33   fPtBiasJet2Clus(0),
34   fAreCollections1MC(kFALSE),  
35   fAreCollections2MC(kTRUE),
36   fMatching(kNoMatching),
37   fMatchingPar1(0),
38   fMatchingPar2(0),
39   fJet2MinEta(-999),
40   fJet2MaxEta(-999),
41   fJet2MinPhi(-999),
42   fJet2MaxPhi(-999),
43   fSelectPtHardBin(-999),
44   fPythiaHeader(0),
45   fPtHardBin(0),
46   fNTrials(0),
47   fTracks2(0),
48   fCaloClusters2(0),
49   fJets2(0),
50   fRho2(0),
51   fRho2Val(0),
52   fTracks2Map(0),
53   fHistNTrials(0),
54   fHistEvents(0),
55   fHistJets1PhiEta(0),
56   fHistJets1PtArea(0),
57   fHistJets1CorrPtArea(0),
58   fHistJets2PhiEta(0),
59   fHistJets2PtArea(0),
60   fHistJets2CorrPtArea(0),
61   fHistJets2PhiEtaAcceptance(0),
62   fHistJets2PtAreaAcceptance(0),
63   fHistJets2CorrPtAreaAcceptance(0),
64   fHistMatchingLevelvsJet2Pt(0),
65   fHistDistancevsCommonEnergy(0),
66   fHistDeltaEtaPhivsJet2Pt(0),
67   fHistDeltaPtvsJet2Pt(0),
68   fHistDeltaPtvsMatchingLevel(0),
69   fHistDeltaCorrPtvsJet2Pt(0),
70   fHistDeltaCorrPtvsMatchingLevel(0),
71   fHistNonMatchedJets1PtArea(0),
72   fHistNonMatchedJets2PtArea(0),
73   fHistNonMatchedJets1CorrPtArea(0),
74   fHistNonMatchedJets2CorrPtArea(0),
75   fHistJet1PtvsJet2Pt(0),
76   fHistJet1CorrPtvsJet2CorrPt(0),
77   fHistMissedJets2PtArea(0)
78 {
79   // Default constructor.
80
81   SetMakeGeneralHistograms(kTRUE);
82 }
83
84 //________________________________________________________________________
85 AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
86   AliAnalysisTaskEmcalJet(name, kTRUE),
87   fTracks2Name("MCParticles"),
88   fCalo2Name(""),
89   fJets2Name("MCJets"),
90   fRho2Name(""),
91   fPtBiasJet2Track(0),
92   fPtBiasJet2Clus(0),
93   fAreCollections1MC(kFALSE),  
94   fAreCollections2MC(kTRUE),
95   fMatching(kNoMatching),
96   fMatchingPar1(0),
97   fMatchingPar2(0),
98   fJet2MinEta(-999),
99   fJet2MaxEta(-999),
100   fJet2MinPhi(-999),
101   fJet2MaxPhi(-999),
102   fSelectPtHardBin(-999),
103   fPythiaHeader(0),
104   fPtHardBin(0),
105   fNTrials(0),
106   fTracks2(0),
107   fCaloClusters2(0),
108   fJets2(0),
109   fRho2(0),
110   fRho2Val(0),
111   fTracks2Map(0),
112   fHistNTrials(0),
113   fHistEvents(0),
114   fHistJets1PhiEta(0),
115   fHistJets1PtArea(0),
116   fHistJets1CorrPtArea(0),
117   fHistJets2PhiEta(0),
118   fHistJets2PtArea(0),
119   fHistJets2CorrPtArea(0),
120   fHistJets2PhiEtaAcceptance(0),
121   fHistJets2PtAreaAcceptance(0),
122   fHistJets2CorrPtAreaAcceptance(0),
123   fHistMatchingLevelvsJet2Pt(0),
124   fHistDistancevsCommonEnergy(0),
125   fHistDeltaEtaPhivsJet2Pt(0),
126   fHistDeltaPtvsJet2Pt(0),
127   fHistDeltaPtvsMatchingLevel(0),
128   fHistDeltaCorrPtvsJet2Pt(0),
129   fHistDeltaCorrPtvsMatchingLevel(0),
130   fHistNonMatchedJets1PtArea(0),
131   fHistNonMatchedJets2PtArea(0),
132   fHistNonMatchedJets1CorrPtArea(0),
133   fHistNonMatchedJets2CorrPtArea(0),
134   fHistJet1PtvsJet2Pt(0),
135   fHistJet1CorrPtvsJet2CorrPt(0),
136   fHistMissedJets2PtArea(0)
137 {
138   // Standard constructor.
139
140   SetMakeGeneralHistograms(kTRUE);
141 }
142
143 //________________________________________________________________________
144 AliJetResponseMaker::~AliJetResponseMaker()
145 {
146   // Destructor
147 }
148
149 //________________________________________________________________________
150 void AliJetResponseMaker::UserCreateOutputObjects()
151 {
152   // Create user objects.
153
154   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
155
156   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
157   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
158
159   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
160   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
161   fHistNTrials->GetYaxis()->SetTitle("trials");
162   fOutput->Add(fHistNTrials);
163
164   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
165   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
166   fHistEvents->GetYaxis()->SetTitle("total events");
167   fOutput->Add(fHistEvents);
168
169   for (Int_t i = 1; i < 12; i++) {
170     fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
171     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
172   }
173
174   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
175   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
176   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
177   fOutput->Add(fHistJets1PhiEta);
178   
179   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
180   fHistJets1PtArea->GetXaxis()->SetTitle("area");
181   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
182   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
183   fOutput->Add(fHistJets1PtArea);
184
185   if (!fRhoName.IsNull()) {
186     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
187     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
188     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
189     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
190     fOutput->Add(fHistJets1CorrPtArea);
191   }
192
193   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
194   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
195   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
196   fOutput->Add(fHistJets2PhiEta);
197
198   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
199   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
200   fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
201   fOutput->Add(fHistJets2PhiEtaAcceptance);
202   
203   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
204   fHistJets2PtArea->GetXaxis()->SetTitle("area");
205   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
206   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
207   fOutput->Add(fHistJets2PtArea);
208
209   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
210   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
211   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
212   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
213   fOutput->Add(fHistJets2PtAreaAcceptance);
214
215   if (!fRho2Name.IsNull()) {
216     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
217     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
218     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
219     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
220     fOutput->Add(fHistJets2CorrPtArea);
221
222     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
223     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
224     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
225     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
226     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
227   }
228
229   fHistMatchingLevelvsJet2Pt = new TH2F("fHistMatchingLevelvsJet2Pt", "fHistMatchingLevelvsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
230   fHistMatchingLevelvsJet2Pt->GetXaxis()->SetTitle("Matching level");
231   fHistMatchingLevelvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");  
232   fHistMatchingLevelvsJet2Pt->GetZaxis()->SetTitle("counts");
233   fOutput->Add(fHistMatchingLevelvsJet2Pt);
234
235   fHistDistancevsCommonEnergy = new TH2F("fHistDistancevsCommonEnergy", "fHistDistancevsCommonEnergy", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
236   fHistDistancevsCommonEnergy->GetXaxis()->SetTitle("Distance");
237   fHistDistancevsCommonEnergy->GetYaxis()->SetTitle("Common energy (%)");  
238   fHistDistancevsCommonEnergy->GetZaxis()->SetTitle("counts");
239   fOutput->Add(fHistDistancevsCommonEnergy);
240
241   fHistDeltaEtaPhivsJet2Pt = new TH3F("fHistDeltaEtaPhivsJet2Pt", "fHistDeltaEtaPhivsJet2Pt", 40, -1, 1, 128, -1.6, 4.8, fNbins/2, fMinBinPt, fMaxBinPt);
242   fHistDeltaEtaPhivsJet2Pt->GetXaxis()->SetTitle("#Delta#eta");
243   fHistDeltaEtaPhivsJet2Pt->GetYaxis()->SetTitle("#Delta#phi");
244   fHistDeltaEtaPhivsJet2Pt->GetZaxis()->SetTitle("p_{T,2}");
245   fOutput->Add(fHistDeltaEtaPhivsJet2Pt);
246
247   fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
248   fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
249   fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
250   fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
251   fOutput->Add(fHistDeltaPtvsJet2Pt);
252
253   fHistDeltaPtvsMatchingLevel = new TH2F("fHistDeltaPtvsMatchingLevel", "fHistDeltaPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
254   fHistDeltaPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");  
255   fHistDeltaPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T} (GeV/c)");
256   fHistDeltaPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
257   fOutput->Add(fHistDeltaPtvsMatchingLevel);
258
259   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {  
260     fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
261     fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
262     fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
263     fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
264     fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
265
266     fHistDeltaCorrPtvsMatchingLevel = new TH2F("fHistDeltaCorrPtvsMatchingLevel", "fHistDeltaCorrPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
267     fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");  
268     fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
269     fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
270     fOutput->Add(fHistDeltaCorrPtvsMatchingLevel);
271   }
272
273   fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
274   fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
275   fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
276   fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
277   fOutput->Add(fHistNonMatchedJets1PtArea);
278
279   fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
280   fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
281   fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
282   fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
283   fOutput->Add(fHistNonMatchedJets2PtArea);
284
285   if (!fRhoName.IsNull()) {  
286     fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
287     fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
288     fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
289     fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
290     fOutput->Add(fHistNonMatchedJets1CorrPtArea);
291   }
292
293   if (!fRho2Name.IsNull()) {  
294     fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
295     fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
296     fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
297     fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
298     fOutput->Add(fHistNonMatchedJets2CorrPtArea);
299   }
300
301   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
302   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
303   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
304   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
305   fOutput->Add(fHistJet1PtvsJet2Pt);
306
307   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
308     if (fRhoName.IsNull()) 
309       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
310     else if (fRho2Name.IsNull()) 
311       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
312     else
313       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
314     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
315     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
316     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
317     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
318   }
319
320   fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
321   fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");  
322   fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
323   fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
324   fOutput->Add(fHistMissedJets2PtArea);
325
326   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
327 }
328
329 //________________________________________________________________________
330 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
331 {   
332   // Return true if jet is accepted.
333
334   if (jet->Pt() <= fJetPtCut)
335     return kFALSE;
336   if (jet->Area() <= fJetAreaCut)
337     return kFALSE;
338
339   return kTRUE;
340 }
341
342 //________________________________________________________________________
343 Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
344
345   // Accept jet with a bias.
346
347   if (fLeadingHadronType == 0) {
348     if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
349   }
350   else if (fLeadingHadronType == 1) {
351     if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
352   }
353   else {
354     if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
355   }
356
357   return kTRUE;
358 }
359
360 //________________________________________________________________________
361 void AliJetResponseMaker::ExecOnce()
362 {
363   // Execute once.
364
365   if (!fJets2Name.IsNull() && !fJets2) {
366     fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
367     if (!fJets2) {
368       AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
369       return;
370     }
371     else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
372       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data())); 
373       fJets2 = 0;
374       return;
375     }
376   }
377
378   if (!fTracks2Name.IsNull() && !fTracks2) {
379     fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
380     if (!fTracks2) {
381       AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data())); 
382       return;
383     }
384     else {
385       TClass *cl = fTracks2->GetClass();
386       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
387         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data())); 
388         fTracks2 = 0;
389         return;
390       }
391     }
392
393     if (fAreCollections2MC) {
394       fTracks2Map = dynamic_cast<TH1*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
395       // this is needed to map the MC labels with the indexes of the MC particle collection
396       // if teh map is not given, the MC labels are assumed to be consistent with the indexes (which is not the case if AliEmcalMCTrackSelector is used)
397       if (!fTracks2Map) {
398         AliWarning(Form("%s: Could not retrieve map for tracks2 %s! Will assume MC labels consistent with indexes...", GetName(), fTracks2Name.Data())); 
399         fTracks2Map = new TH1I("tracksMap","tracksMap",9999,0,1);
400         for (Int_t i = 0; i < 9999; i++) {
401           fTracks2Map->SetBinContent(i,i);
402         }
403       }
404     }
405   }
406
407   if (!fCalo2Name.IsNull() && !fCaloClusters2) {
408     fCaloClusters2 =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
409     if (!fCaloClusters2) {
410       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data())); 
411       return;
412     } else {
413       TClass *cl = fCaloClusters2->GetClass();
414       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
415         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data())); 
416         fCaloClusters2 = 0;
417         return;
418       }
419     }
420   }
421
422   if (!fRho2Name.IsNull() && !fRho2) {
423     fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
424     if (!fRho2) {
425       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
426       fInitialized = kFALSE;
427       return;
428     }
429   }
430
431   if (fJet2MinEta == -999)
432     fJet2MinEta = fJetMinEta - fJetRadius;
433   if (fJet2MaxEta == -999)
434     fJet2MaxEta = fJetMaxEta + fJetRadius;
435   if (fJet2MinPhi == -999)
436     fJet2MinPhi = fJetMinPhi - fJetRadius;
437   if (fJet2MaxPhi == -999)
438     fJet2MaxPhi = fJetMaxPhi + fJetRadius;
439
440   AliAnalysisTaskEmcalJet::ExecOnce();
441 }
442
443 //________________________________________________________________________
444 Bool_t AliJetResponseMaker::IsEventSelected()
445 {
446   // Check if event is selected
447
448   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) 
449     return kFALSE;
450
451   return AliAnalysisTaskEmcalJet::IsEventSelected();
452 }
453
454 //________________________________________________________________________
455 Bool_t AliJetResponseMaker::RetrieveEventObjects()
456 {
457   // Retrieve event objects.
458
459   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
460     return kFALSE;
461
462   if (fRho2)
463     fRho2Val = fRho2->GetVal();
464
465   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
466   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
467   
468   if (MCEvent())
469     fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
470
471   if (fPythiaHeader) {
472     Double_t pthard = fPythiaHeader->GetPtHard();
473     
474     for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
475       if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
476         break;
477     }
478     
479     fNTrials = fPythiaHeader->Trials();
480   }
481
482   return kTRUE;
483 }
484
485 //________________________________________________________________________
486 Bool_t AliJetResponseMaker::Run()
487 {
488   // Find the closest jets
489
490   if (fMatching == kNoMatching) 
491     return kTRUE;
492   else
493     return DoJetMatching();
494 }
495
496 //________________________________________________________________________
497 Bool_t AliJetResponseMaker::DoJetMatching()
498 {
499   DoJetLoop(kFALSE);
500
501   const Int_t nJets = fJets->GetEntriesFast();
502
503   for (Int_t i = 0; i < nJets; i++) {
504
505     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
506
507     if (!jet1) {
508       AliError(Form("Could not receive jet %d", i));
509       continue;
510     }  
511
512     if (!AcceptJet(jet1))
513       continue;
514
515     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
516       continue;
517
518     if (jet1->Pt() > fMaxBinPt)
519       continue;
520
521     if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 && 
522         jet1->ClosestJetDistance() < fMatchingPar1 && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar2) {    // Matched jet found
523       jet1->SetMatchedToClosest(fMatching);
524       jet1->ClosestJet()->SetMatchedToClosest(fMatching);
525     }
526   }
527
528   return kTRUE;
529 }
530
531 //________________________________________________________________________
532 void AliJetResponseMaker::DoJetLoop(Bool_t order)
533 {
534   // Do the jet loop.
535
536   TClonesArray *jets1 = 0;
537   TClonesArray *jets2 = 0;
538
539   if (order) {
540     jets1 = fJets2;
541     jets2 = fJets;
542   }
543   else {
544     jets1 = fJets;
545     jets2 = fJets2;
546   }
547
548   Int_t nJets1 = jets1->GetEntriesFast();
549   Int_t nJets2 = jets2->GetEntriesFast();
550
551   for (Int_t i = 0; i < nJets1; i++) {
552
553     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
554
555     if (!jet1) {
556       AliError(Form("Could not receive jet %d", i));
557       continue;
558     }  
559
560     if (!AcceptJet(jet1))
561       continue;
562
563     if (order) {
564      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
565         continue;
566     }
567     else {
568       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
569         continue;
570     }
571
572     for (Int_t j = 0; j < nJets2; j++) {
573       
574       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
575       
576       if (!jet2) {
577         AliError(Form("Could not receive jet %d", j));
578         continue;
579       }  
580       
581       if (!AcceptJet(jet2))
582         continue;
583
584       if (order) {
585         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
586           continue;
587       }
588       else {
589         if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
590           continue;
591       }
592
593       SetMatchingLevel(jet1, jet2, fMatching);
594     }
595   }
596 }
597
598 //________________________________________________________________________
599 void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
600 {
601   Double_t deta = jet2->Eta() - jet1->Eta();
602   Double_t dphi = jet2->Phi() - jet1->Phi();
603   d = TMath::Sqrt(deta * deta + dphi * dphi);
604 }
605
606 //________________________________________________________________________
607 void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
608
609   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
610   d1 = jet1->Pt();
611   d2 = jet2->Pt();
612   Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
613
614   for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
615     Bool_t track2Found = kFALSE;
616     Int_t index2 = jet2->TrackAt(iTrack2);
617     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
618       AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
619       if (!track) {
620         AliWarning(Form("Could not find track %d!", iTrack));
621         continue;
622       }
623       Int_t MClabel = TMath::Abs(track->GetLabel());
624       Int_t index = -1;
625           
626       if (MClabel == 0) {// this is not a MC particle; remove it completely
627         AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
628         totalPt1 -= track->Pt();
629         d1 -= track->Pt();
630         continue;
631       }
632       else if (MClabel < fTracks2Map->GetNbinsX()-2) {
633         index = fTracks2Map->GetBinContent(MClabel);
634       }
635           
636       if (index < 0) {
637         AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
638         continue;
639       }
640
641       if (index2 == index) { // found common particle
642         track2Found = kTRUE;
643         d1 -= track->Pt();
644         AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
645         AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
646                         iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
647         d2 -= MCpart->Pt();
648         break;
649       }
650     }
651     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
652       AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
653       if (!clus) {
654         AliWarning(Form("Could not find cluster %d!", iClus));
655         continue;
656       }
657       TLorentzVector part;
658       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
659           
660       if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
661         for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
662           Int_t cellId = clus->GetCellAbsId(iCell);
663           Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
664
665           Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
666           Int_t index = -1;
667           
668           if (MClabel == 0) {// this is not a MC particle; remove it completely
669             AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
670             totalPt1 -= part.Pt() * cellFrac;
671             d1 -= part.Pt() * cellFrac;
672             continue;
673           }
674           else if (MClabel < fTracks2Map->GetNbinsX()-2) {
675             index = fTracks2Map->GetBinContent(MClabel);
676           }
677
678           if (index < 0) {
679             AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
680             continue;
681           }
682           if (index2 == index) { // found common particle
683             d1 -= part.Pt() * cellFrac;
684                 
685             if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
686               AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
687               AliDebug(3,Form("Cell %d belonging to cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
688                               iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));             
689               d2 -= MCpart->Pt() * cellFrac;
690             }
691             break;
692           }
693         }
694       }
695       else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
696         Int_t MClabel = TMath::Abs(clus->GetLabel());
697         Int_t index = -1;
698             
699         if (MClabel == 0) {// this is not a MC particle; remove it completely
700           AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
701           totalPt1 -= part.Pt();
702           d1 -= part.Pt();
703           continue;
704         }
705         else if (MClabel < fTracks2Map->GetNbinsX()-2) {
706           index = fTracks2Map->GetBinContent(MClabel);
707         }
708          
709         if (index < 0) {
710           AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
711           continue;
712         }
713         if (index2 == index) { // found common particle
714           d1 -= part.Pt();
715
716           if (!track2Found) {// only if it is not already found among charged tracks (charged particles are most likely already found)
717             AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
718             AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
719                             iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
720                 
721             d2 -= MCpart->Pt();
722           }
723           break;
724         }
725       }
726     }
727   }
728   if (d1 <= 0 || totalPt1 < 1)
729     d1 = 0;
730   else
731     d1 /= totalPt1;
732
733   if (jet2->Pt() > 0 && d2 > 0)
734     d2 /= jet2->Pt();
735   else
736     d2 = 0;
737 }
738
739 //________________________________________________________________________
740 void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
741
742   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
743   d1 = jet1->Pt();
744   d2 = jet2->Pt();
745
746   if (fTracks && fTracks2) {
747
748     for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
749       Int_t index2 = jet2->TrackAt(iTrack2);
750       for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
751         Int_t index = jet1->TrackAt(iTrack);
752         if (index2 == index) { // found common particle
753           AliVParticle *part = static_cast<AliVParticle*>(fTracks->At(index));
754           if (!part) {
755             AliWarning(Form("Could not find track %d!", index));
756             continue;
757           }
758           AliVParticle *part2 = static_cast<AliVParticle*>(fTracks2->At(index2));
759           if (!part2) {
760             AliWarning(Form("Could not find track %d!", index2));
761             continue;
762           }
763
764           d1 -= part->Pt();
765           d2 -= part2->Pt();
766           break;
767         }
768       }
769     }
770
771   }
772
773   if (fCaloClusters && fCaloClusters2) {
774
775     for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
776       Int_t index2 = jet2->ClusterAt(iClus2);
777       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
778         Int_t index = jet1->ClusterAt(iClus);
779         AliVCluster *clus =  static_cast<AliVCluster*>(fCaloClusters->At(index));
780         if (!clus) {
781           AliWarning(Form("Could not find cluster %d!", index));
782           continue;
783         }
784         AliVCluster *clus2 =  static_cast<AliVCluster*>(fCaloClusters2->At(index2));
785         if (!clus2) {
786           AliWarning(Form("Could not find cluster %d!", index2));
787           continue;
788         }
789         TLorentzVector part, part2;
790         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
791         clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
792
793         d1 -= part.Pt();
794         d2 -= part2.Pt();
795         break;
796       }
797     }
798
799   }
800
801   if (jet1->Pt() > 0 && d1 > 0)
802     d1 /= jet1->Pt();
803   else
804     d1 = 0;
805
806   if (jet2->Pt() > 0 && d2 > 0)
807     d2 /= jet2->Pt();
808   else
809     d2 = 0;
810 }
811
812 //________________________________________________________________________
813 void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
814 {
815   Double_t d1 = -1;
816   Double_t d2 = -1;
817
818   switch (matching) {
819   case kGeometrical:
820     GetGeometricalMatchingLevel(jet1,jet2,d1);
821     d2 = d1;
822     break;
823   case kMCLabel: // jet1 = detector level and jet2 = particle level!
824     GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
825     break;
826   case kSameCollections:
827     GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
828     break;
829   default:
830     ;
831   }
832
833   if (d1 > 0) {
834
835     if (d1 < jet1->ClosestJetDistance()) {
836       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
837       jet1->SetClosestJet(jet2, d1);
838     }
839     else if (d1 < jet1->SecondClosestJetDistance()) {
840       jet1->SetSecondClosestJet(jet2, d1);
841     }
842   }
843   
844   if (d2 > 0) {
845     
846     if (d2 < jet2->ClosestJetDistance()) {
847       jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
848       jet2->SetClosestJet(jet1, d2);
849     }
850     else if (d2 < jet2->SecondClosestJetDistance()) {
851       jet2->SetSecondClosestJet(jet1, d2);
852     }
853   }
854 }
855
856 //________________________________________________________________________
857 Bool_t AliJetResponseMaker::FillHistograms()
858 {
859   // Fill histograms.
860
861   fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
862   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
863
864   const Int_t nJets2 = fJets2->GetEntriesFast();
865
866   for (Int_t i = 0; i < nJets2; i++) {
867
868     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
869
870     if (!jet2) {
871       AliError(Form("Could not receive jet2 %d", i));
872       continue;
873     }
874
875     if (jet2->Pt() > fMaxBinPt)
876       continue;
877
878     if (!AcceptJet(jet2))
879       continue;
880
881     if (AcceptBiasJet(jet2) &&
882         (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
883       
884       fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
885       fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
886       
887       if (!fRho2Name.IsNull())
888         fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
889     }
890
891     if (!AcceptBiasJet2(jet2))
892       continue;
893
894     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
895       continue;
896
897     fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
898     fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
899
900     if (!fRho2Name.IsNull())
901       fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
902
903     if (jet2->MatchedJet()) {
904
905       if (!AcceptBiasJet(jet2->MatchedJet()) || 
906           jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
907           jet2->MatchedJet()->Pt() > fMaxBinPt) {
908         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
909       }
910       else {
911         Double_t d1=-1, d2=-1;
912         if (jet2->GetMatchingType() == kGeometrical) {
913
914           if (fAreCollections2MC && !fAreCollections1MC)
915             GetMCLabelMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
916           else if ((fAreCollections1MC && fAreCollections2MC) || (!fAreCollections1MC && !fAreCollections2MC))
917             GetSameCollectionsMatchingLevel(jet2->MatchedJet(), jet2, d1, d2);
918
919           fHistDistancevsCommonEnergy->Fill(jet2->ClosestJetDistance(), d2);
920         }
921         else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
922           GetGeometricalMatchingLevel(jet2->MatchedJet(), jet2, d1);
923           fHistDistancevsCommonEnergy->Fill(d1, jet2->ClosestJetDistance());
924         }
925           
926         fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
927
928         Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
929         Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
930         fHistDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
931
932         Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
933         fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
934         fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
935
936         fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
937         
938         if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
939           dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
940           fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
941           fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
942           fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
943         }
944       }
945     }
946     else {
947       fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
948       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
949
950       if (!fRho2Name.IsNull())
951         fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
952     }
953   }
954
955   const Int_t nJets1 = fJets->GetEntriesFast();
956
957   for (Int_t i = 0; i < nJets1; i++) {
958
959     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
960
961     if (!jet1) {
962       AliError(Form("Could not receive jet %d", i));
963       continue;
964     }  
965     
966     if (!AcceptJet(jet1))
967       continue;
968
969     if (!AcceptBiasJet(jet1))
970       continue;
971
972     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
973       continue;
974
975     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
976       continue;
977
978     if (!jet1->MatchedJet()) {
979       fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
980       if (!fRhoName.IsNull())
981         fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
982     }
983
984     fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
985     fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
986
987     if (!fRhoName.IsNull())
988       fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
989   }
990
991   return kTRUE;
992 }