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