Removed depricated analysis framework, ie renamed
[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 <TH2F.h>
11 #include <THnSparse.h>
12 #include <TLorentzVector.h>
13
14 #include "AliAnalysisManager.h"
15 #include "AliVCluster.h"
16 #include "AliVTrack.h"
17 #include "AliEmcalJet.h"
18 #include "AliLog.h"
19 #include "AliRhoParameter.h"
20 #include "AliNamedArrayI.h"
21 #include "AliJetContainer.h"
22 #include "AliParticleContainer.h"
23 #include "AliClusterContainer.h"
24
25 ClassImp(AliJetResponseMaker)
26
27 //________________________________________________________________________
28 AliJetResponseMaker::AliJetResponseMaker() : 
29   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
30   fMatching(kNoMatching),
31   fMatchingPar1(0),
32   fMatchingPar2(0),
33   fUseCellsToMatch(kFALSE),
34   fMinJetMCPt(1),
35   fHistoType(0),
36   fDeltaPtAxis(0),
37   fDeltaEtaDeltaPhiAxis(0),
38   fNEFAxis(0),
39   fZAxis(0),
40   fDoJet2Histogram(0),
41   fIsJet1Rho(kFALSE),
42   fIsJet2Rho(kFALSE),
43   fHistLeadingJets1PtArea(0),
44   fHistLeadingJets1CorrPtArea(0),
45   fHistLeadingJets2PtArea(0),
46   fHistLeadingJets2CorrPtArea(0),
47   fHistLeadingJets2PtAreaAcceptance(0),
48   fHistLeadingJets2CorrPtAreaAcceptance(0),
49   fHistJets1(0),
50   fHistJets2(0),
51   fHistJets2Acceptance(0),
52   fHistMatching(0),
53   fHistJets1PhiEta(0),
54   fHistJets1PtArea(0),
55   fHistJets1CorrPtArea(0),
56   fHistJets1NEFvsPt(0),
57   fHistJets1CEFvsCEFPt(0),
58   fHistJets1ZvsPt(0),
59   fHistJets2PhiEta(0),
60   fHistJets2PtArea(0),
61   fHistJets2CorrPtArea(0),
62   fHistJets2PhiEtaAcceptance(0),
63   fHistJets2PtAreaAcceptance(0),
64   fHistJets2CorrPtAreaAcceptance(0),
65   fHistJets2NEFvsPt(0),
66   fHistJets2CEFvsCEFPt(0),
67   fHistJets2ZvsPt(0),
68   fHistCommonEnergy1vsJet1Pt(0),
69   fHistCommonEnergy2vsJet2Pt(0),
70   fHistDistancevsJet1Pt(0),
71   fHistDistancevsJet2Pt(0),
72   fHistDistancevsCommonEnergy1(0),
73   fHistDistancevsCommonEnergy2(0),
74   fHistCommonEnergy1vsCommonEnergy2(0),
75   fHistDeltaEtaDeltaPhi(0),
76   fHistJet2PtOverJet1PtvsJet2Pt(0),
77   fHistJet1PtOverJet2PtvsJet1Pt(0),
78   fHistDeltaPtvsJet1Pt(0),
79   fHistDeltaPtvsJet2Pt(0),
80   fHistDeltaPtOverJet1PtvsJet1Pt(0),
81   fHistDeltaPtOverJet2PtvsJet2Pt(0),
82   fHistDeltaPtvsDistance(0),
83   fHistDeltaPtvsCommonEnergy1(0),
84   fHistDeltaPtvsCommonEnergy2(0),
85   fHistDeltaPtvsArea1(0),
86   fHistDeltaPtvsArea2(0),
87   fHistDeltaPtvsDeltaArea(0),
88   fHistJet1PtvsJet2Pt(0),
89   fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
90   fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
91   fHistDeltaCorrPtvsJet1CorrPt(0),
92   fHistDeltaCorrPtvsJet2CorrPt(0),
93   fHistDeltaCorrPtvsDistance(0),
94   fHistDeltaCorrPtvsCommonEnergy1(0),
95   fHistDeltaCorrPtvsCommonEnergy2(0),
96   fHistDeltaCorrPtvsArea1(0),
97   fHistDeltaCorrPtvsArea2(0),
98   fHistDeltaCorrPtvsDeltaArea(0),
99   fHistJet1CorrPtvsJet2CorrPt(0),
100   fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
101   fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
102   fHistDeltaMCPtvsJet1MCPt(0),
103   fHistDeltaMCPtvsJet2Pt(0),
104   fHistDeltaMCPtvsDistance(0),
105   fHistDeltaMCPtvsCommonEnergy1(0),
106   fHistDeltaMCPtvsCommonEnergy2(0),
107   fHistDeltaMCPtvsArea1(0),
108   fHistDeltaMCPtvsArea2(0),
109   fHistDeltaMCPtvsDeltaArea(0),
110   fHistJet1MCPtvsJet2Pt(0)
111 {
112   // Default constructor.
113
114   SetMakeGeneralHistograms(kTRUE);
115 }
116
117 //________________________________________________________________________
118 AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
119   AliAnalysisTaskEmcalJet(name, kTRUE),
120   fMatching(kNoMatching),
121   fMatchingPar1(0),
122   fMatchingPar2(0),
123   fUseCellsToMatch(kFALSE),
124   fMinJetMCPt(1),
125   fHistoType(0),
126   fDeltaPtAxis(0),
127   fDeltaEtaDeltaPhiAxis(0),
128   fNEFAxis(0),
129   fZAxis(0),
130   fDoJet2Histogram(0),
131   fIsJet1Rho(kFALSE),
132   fIsJet2Rho(kFALSE),
133   fHistLeadingJets1PtArea(0),
134   fHistLeadingJets1CorrPtArea(0),
135   fHistLeadingJets2PtArea(0),
136   fHistLeadingJets2CorrPtArea(0),
137   fHistLeadingJets2PtAreaAcceptance(0),
138   fHistLeadingJets2CorrPtAreaAcceptance(0),
139   fHistJets1(0),
140   fHistJets2(0),
141   fHistJets2Acceptance(0),
142   fHistMatching(0),
143   fHistJets1PhiEta(0),
144   fHistJets1PtArea(0),
145   fHistJets1CorrPtArea(0),
146   fHistJets1NEFvsPt(0),
147   fHistJets1CEFvsCEFPt(0),
148   fHistJets1ZvsPt(0),
149   fHistJets2PhiEta(0),
150   fHistJets2PtArea(0),
151   fHistJets2CorrPtArea(0),
152   fHistJets2PhiEtaAcceptance(0),
153   fHistJets2PtAreaAcceptance(0),
154   fHistJets2CorrPtAreaAcceptance(0),
155   fHistJets2NEFvsPt(0),
156   fHistJets2CEFvsCEFPt(0),
157   fHistJets2ZvsPt(0),
158   fHistCommonEnergy1vsJet1Pt(0),
159   fHistCommonEnergy2vsJet2Pt(0),
160   fHistDistancevsJet1Pt(0),
161   fHistDistancevsJet2Pt(0),
162   fHistDistancevsCommonEnergy1(0),
163   fHistDistancevsCommonEnergy2(0),
164   fHistCommonEnergy1vsCommonEnergy2(0),
165   fHistDeltaEtaDeltaPhi(0),
166   fHistJet2PtOverJet1PtvsJet2Pt(0),
167   fHistJet1PtOverJet2PtvsJet1Pt(0),
168   fHistDeltaPtvsJet1Pt(0),
169   fHistDeltaPtvsJet2Pt(0),
170   fHistDeltaPtOverJet1PtvsJet1Pt(0),
171   fHistDeltaPtOverJet2PtvsJet2Pt(0),
172   fHistDeltaPtvsDistance(0),
173   fHistDeltaPtvsCommonEnergy1(0),
174   fHistDeltaPtvsCommonEnergy2(0),
175   fHistDeltaPtvsArea1(0),
176   fHistDeltaPtvsArea2(0),
177   fHistDeltaPtvsDeltaArea(0),
178   fHistJet1PtvsJet2Pt(0),
179   fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt(0),
180   fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt(0),
181   fHistDeltaCorrPtvsJet1CorrPt(0),
182   fHistDeltaCorrPtvsJet2CorrPt(0),
183   fHistDeltaCorrPtvsDistance(0),
184   fHistDeltaCorrPtvsCommonEnergy1(0),
185   fHistDeltaCorrPtvsCommonEnergy2(0),
186   fHistDeltaCorrPtvsArea1(0),
187   fHistDeltaCorrPtvsArea2(0),
188   fHistDeltaCorrPtvsDeltaArea(0),
189   fHistJet1CorrPtvsJet2CorrPt(0),
190   fHistDeltaMCPtOverJet1MCPtvsJet1MCPt(0),
191   fHistDeltaMCPtOverJet2PtvsJet2Pt(0),
192   fHistDeltaMCPtvsJet1MCPt(0),
193   fHistDeltaMCPtvsJet2Pt(0),
194   fHistDeltaMCPtvsDistance(0),
195   fHistDeltaMCPtvsCommonEnergy1(0),
196   fHistDeltaMCPtvsCommonEnergy2(0),
197   fHistDeltaMCPtvsArea1(0),
198   fHistDeltaMCPtvsArea2(0),
199   fHistDeltaMCPtvsDeltaArea(0),
200   fHistJet1MCPtvsJet2Pt(0)
201 {
202   // Standard constructor.
203
204   SetMakeGeneralHistograms(kTRUE);
205 }
206
207 //________________________________________________________________________
208 AliJetResponseMaker::~AliJetResponseMaker()
209 {
210   // Destructor
211 }
212
213
214 //________________________________________________________________________
215 void AliJetResponseMaker::AllocateTH2()
216 {
217   // Allocate TH2 histograms.
218
219   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
220   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
221   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
222   fOutput->Add(fHistJets1PhiEta);
223   
224   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
225   fHistJets1PtArea->GetXaxis()->SetTitle("area");
226   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
227   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
228   fOutput->Add(fHistJets1PtArea);
229  
230   if (fIsJet1Rho) {
231     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 
232                                     fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
233     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
234     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
235     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
236     fOutput->Add(fHistJets1CorrPtArea);
237   }
238
239   fHistJets1ZvsPt = new TH2F("fHistJets1ZvsPt", "fHistJets1ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
240   fHistJets1ZvsPt->GetXaxis()->SetTitle("Z");
241   fHistJets1ZvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
242   fHistJets1ZvsPt->GetZaxis()->SetTitle("counts");
243   fOutput->Add(fHistJets1ZvsPt);
244   
245   fHistJets1NEFvsPt = new TH2F("fHistJets1NEFvsPt", "fHistJets1NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
246   fHistJets1NEFvsPt->GetXaxis()->SetTitle("NEF");
247   fHistJets1NEFvsPt->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
248   fHistJets1NEFvsPt->GetZaxis()->SetTitle("counts");
249   fOutput->Add(fHistJets1NEFvsPt);
250   
251   fHistJets1CEFvsCEFPt = new TH2F("fHistJets1CEFvsCEFPt", "fHistJets1CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
252   fHistJets1CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
253   fHistJets1CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,1} (GeV/c)");
254   fHistJets1CEFvsCEFPt->GetZaxis()->SetTitle("counts");
255   fOutput->Add(fHistJets1CEFvsCEFPt);
256
257   // Jets 2 histograms
258
259   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
260   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
261   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
262   fOutput->Add(fHistJets2PhiEta);
263
264   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
265   fHistJets2PtArea->GetXaxis()->SetTitle("area");
266   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
267   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
268   fOutput->Add(fHistJets2PtArea);
269
270   if (fIsJet2Rho) {
271     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 
272                                     fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
273     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
274     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
275     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
276     fOutput->Add(fHistJets2CorrPtArea);
277   }
278
279   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
280   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
281   fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
282   fOutput->Add(fHistJets2PhiEtaAcceptance);
283
284   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 
285                                         fNbins/2, 0, 1.5, fNbins, fMinBinPt, fMaxBinPt);
286   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
287   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
288   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
289   fOutput->Add(fHistJets2PtAreaAcceptance);
290
291   if (fIsJet2Rho) {
292     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 
293                                               fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
294     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
295     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
296     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
297     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
298   }
299
300   fHistJets2ZvsPt = new TH2F("fHistJets2ZvsPt", "fHistJets2ZvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
301   fHistJets2ZvsPt->GetXaxis()->SetTitle("Z");
302   fHistJets2ZvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
303   fHistJets2ZvsPt->GetZaxis()->SetTitle("counts");
304   fOutput->Add(fHistJets2ZvsPt);
305   
306   fHistJets2NEFvsPt = new TH2F("fHistJets2NEFvsPt", "fHistJets2NEFvsPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
307   fHistJets2NEFvsPt->GetXaxis()->SetTitle("NEF");
308   fHistJets2NEFvsPt->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
309   fHistJets2NEFvsPt->GetZaxis()->SetTitle("counts");
310   fOutput->Add(fHistJets2NEFvsPt);
311   
312   fHistJets2CEFvsCEFPt = new TH2F("fHistJets2CEFvsCEFPt", "fHistJets2CEFvsCEFPt", 120, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
313   fHistJets2CEFvsCEFPt->GetXaxis()->SetTitle("1-NEF");
314   fHistJets2CEFvsCEFPt->GetYaxis()->SetTitle("(1-NEF)*p_{T,2} (GeV/c)");
315   fHistJets2CEFvsCEFPt->GetZaxis()->SetTitle("counts");
316   fOutput->Add(fHistJets2CEFvsCEFPt);
317
318   // Matching histograms
319
320   fHistCommonEnergy1vsJet1Pt = new TH2F("fHistCommonEnergy1vsJet1Pt", "fHistCommonEnergy1vsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
321   fHistCommonEnergy1vsJet1Pt->GetXaxis()->SetTitle("Common energy 1 (%)");
322   fHistCommonEnergy1vsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");  
323   fHistCommonEnergy1vsJet1Pt->GetZaxis()->SetTitle("counts");
324   fOutput->Add(fHistCommonEnergy1vsJet1Pt);
325
326   fHistCommonEnergy2vsJet2Pt = new TH2F("fHistCommonEnergy2vsJet2Pt", "fHistCommonEnergy2vsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
327   fHistCommonEnergy2vsJet2Pt->GetXaxis()->SetTitle("Common energy 2 (%)");
328   fHistCommonEnergy2vsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");  
329   fHistCommonEnergy2vsJet2Pt->GetZaxis()->SetTitle("counts");
330   fOutput->Add(fHistCommonEnergy2vsJet2Pt);
331
332   fHistDistancevsJet1Pt = new TH2F("fHistDistancevsJet1Pt", "fHistDistancevsJet1Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
333   fHistDistancevsJet1Pt->GetXaxis()->SetTitle("Distance");
334   fHistDistancevsJet1Pt->GetYaxis()->SetTitle("p_{T,1}");  
335   fHistDistancevsJet1Pt->GetZaxis()->SetTitle("counts");
336   fOutput->Add(fHistDistancevsJet1Pt);
337
338   fHistDistancevsJet2Pt = new TH2F("fHistDistancevsJet2Pt", "fHistDistancevsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
339   fHistDistancevsJet2Pt->GetXaxis()->SetTitle("Distance");
340   fHistDistancevsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");  
341   fHistDistancevsJet2Pt->GetZaxis()->SetTitle("counts");
342   fOutput->Add(fHistDistancevsJet2Pt);
343
344   fHistDistancevsCommonEnergy1 = new TH2F("fHistDistancevsCommonEnergy1", "fHistDistancevsCommonEnergy1", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
345   fHistDistancevsCommonEnergy1->GetXaxis()->SetTitle("Distance");
346   fHistDistancevsCommonEnergy1->GetYaxis()->SetTitle("Common energy 1 (%)");  
347   fHistDistancevsCommonEnergy1->GetZaxis()->SetTitle("counts");
348   fOutput->Add(fHistDistancevsCommonEnergy1);
349
350   fHistDistancevsCommonEnergy2 = new TH2F("fHistDistancevsCommonEnergy2", "fHistDistancevsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
351   fHistDistancevsCommonEnergy2->GetXaxis()->SetTitle("Distance");
352   fHistDistancevsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");  
353   fHistDistancevsCommonEnergy2->GetZaxis()->SetTitle("counts");
354   fOutput->Add(fHistDistancevsCommonEnergy2);
355
356   fHistCommonEnergy1vsCommonEnergy2 = new TH2F("fHistCommonEnergy1vsCommonEnergy2", "fHistCommonEnergy1vsCommonEnergy2", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
357   fHistCommonEnergy1vsCommonEnergy2->GetXaxis()->SetTitle("Common energy 1 (%)");
358   fHistCommonEnergy1vsCommonEnergy2->GetYaxis()->SetTitle("Common energy 2 (%)");  
359   fHistCommonEnergy1vsCommonEnergy2->GetZaxis()->SetTitle("counts");
360   fOutput->Add(fHistCommonEnergy1vsCommonEnergy2);
361
362   fHistDeltaEtaDeltaPhi = new TH2F("fHistDeltaEtaDeltaPhi", "fHistDeltaEtaDeltaPhi", fNbins/4, -1, 1, fNbins/4, -TMath::Pi()/2, TMath::Pi()*3/2);
363   fHistDeltaEtaDeltaPhi->GetXaxis()->SetTitle("Common energy 1 (%)");
364   fHistDeltaEtaDeltaPhi->GetYaxis()->SetTitle("Common energy 2 (%)");  
365   fHistDeltaEtaDeltaPhi->GetZaxis()->SetTitle("counts");
366   fOutput->Add(fHistDeltaEtaDeltaPhi);
367
368   fHistJet2PtOverJet1PtvsJet2Pt = new TH2F("fHistJet2PtOverJet1PtvsJet2Pt", "fHistJet2PtOverJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
369   fHistJet2PtOverJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
370   fHistJet2PtOverJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2} / p_{T,1}");
371   fHistJet2PtOverJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
372   fOutput->Add(fHistJet2PtOverJet1PtvsJet2Pt);
373
374   fHistJet1PtOverJet2PtvsJet1Pt = new TH2F("fHistJet1PtOverJet2PtvsJet1Pt", "fHistJet1PtOverJet2PtvsJet1Pt", fNbins, fMinBinPt, fMaxBinPt, 300, 0, 1.5);
375   fHistJet1PtOverJet2PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
376   fHistJet1PtOverJet2PtvsJet1Pt->GetYaxis()->SetTitle("p_{T,1} / p_{T,2}");
377   fHistJet1PtOverJet2PtvsJet1Pt->GetZaxis()->SetTitle("counts");
378   fOutput->Add(fHistJet1PtOverJet2PtvsJet1Pt);
379
380   fHistDeltaPtvsJet1Pt = new TH2F("fHistDeltaPtvsJet1Pt", "fHistDeltaPtvsJet1Pt", 
381                                   fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
382   fHistDeltaPtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
383   fHistDeltaPtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
384   fHistDeltaPtvsJet1Pt->GetZaxis()->SetTitle("counts");
385   fOutput->Add(fHistDeltaPtvsJet1Pt);
386
387   fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", 
388                                   fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
389   fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
390   fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
391   fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
392   fOutput->Add(fHistDeltaPtvsJet2Pt);
393
394   fHistDeltaPtOverJet1PtvsJet1Pt = new TH2F("fHistDeltaPtOverJet1PtvsJet1Pt", "fHistDeltaPtOverJet1PtvsJet1Pt", 
395                                             fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
396   fHistDeltaPtOverJet1PtvsJet1Pt->GetXaxis()->SetTitle("p_{T,1}");  
397   fHistDeltaPtOverJet1PtvsJet1Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}");
398   fHistDeltaPtOverJet1PtvsJet1Pt->GetZaxis()->SetTitle("counts");
399   fOutput->Add(fHistDeltaPtOverJet1PtvsJet1Pt);
400
401   fHistDeltaPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaPtOverJet2PtvsJet2Pt", "fHistDeltaPtOverJet2PtvsJet2Pt", 
402                                             fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
403   fHistDeltaPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
404   fHistDeltaPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
405   fHistDeltaPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
406   fOutput->Add(fHistDeltaPtOverJet2PtvsJet2Pt);
407
408   fHistDeltaPtvsDistance = new TH2F("fHistDeltaPtvsDistance", "fHistDeltaPtvsDistance", 
409                                     fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
410   fHistDeltaPtvsDistance->GetXaxis()->SetTitle("Distance");  
411   fHistDeltaPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
412   fHistDeltaPtvsDistance->GetZaxis()->SetTitle("counts");
413   fOutput->Add(fHistDeltaPtvsDistance);
414
415   fHistDeltaPtvsCommonEnergy1 = new TH2F("fHistDeltaPtvsCommonEnergy1", "fHistDeltaPtvsCommonEnergy1",
416                                          fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
417   fHistDeltaPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");  
418   fHistDeltaPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
419   fHistDeltaPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
420   fOutput->Add(fHistDeltaPtvsCommonEnergy1);
421
422   fHistDeltaPtvsCommonEnergy2 = new TH2F("fHistDeltaPtvsCommonEnergy2", "fHistDeltaPtvsCommonEnergy2", 
423                                          fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
424   fHistDeltaPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");  
425   fHistDeltaPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
426   fHistDeltaPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
427   fOutput->Add(fHistDeltaPtvsCommonEnergy2);
428
429   fHistDeltaPtvsArea1 = new TH2F("fHistDeltaPtvsArea1", "fHistDeltaPtvsArea1", 
430                                  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
431   fHistDeltaPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
432   fHistDeltaPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
433   fHistDeltaPtvsArea1->GetZaxis()->SetTitle("counts");
434   fOutput->Add(fHistDeltaPtvsArea1);
435
436   fHistDeltaPtvsArea2 = new TH2F("fHistDeltaPtvsArea2", "fHistDeltaPtvsArea2", 
437                                  fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
438   fHistDeltaPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
439   fHistDeltaPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
440   fHistDeltaPtvsArea2->GetZaxis()->SetTitle("counts");
441   fOutput->Add(fHistDeltaPtvsArea2);
442
443   fHistDeltaPtvsDeltaArea = new TH2F("fHistDeltaPtvsDeltaArea", "fHistDeltaPtvsDeltaArea", 
444                                      fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
445   fHistDeltaPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
446   fHistDeltaPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
447   fHistDeltaPtvsDeltaArea->GetZaxis()->SetTitle("counts");
448   fOutput->Add(fHistDeltaPtvsDeltaArea);
449
450   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
451   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
452   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
453   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
454   fOutput->Add(fHistJet1PtvsJet2Pt);
455
456   if (fIsJet1Rho || fIsJet2Rho) {
457     if (!fIsJet1Rho) 
458       fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", 
459                                               fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
460     else
461       fHistDeltaCorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtvsJet1CorrPt", "fHistDeltaCorrPtvsJet1CorrPt", 
462                                               2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
463       
464     fHistDeltaCorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");  
465     fHistDeltaCorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
466     fHistDeltaCorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
467     fOutput->Add(fHistDeltaCorrPtvsJet1CorrPt);
468
469     if (!fIsJet2Rho) 
470       fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", 
471                                               fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);                                         
472     else
473       fHistDeltaCorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtvsJet2CorrPt", "fHistDeltaCorrPtvsJet2CorrPt", 
474                                               2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
475
476     fHistDeltaCorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");  
477     fHistDeltaCorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
478     fHistDeltaCorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
479     fOutput->Add(fHistDeltaCorrPtvsJet2CorrPt);
480
481     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt = new TH2F("fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", "fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt", 
482                                                           2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
483     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");  
484     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,1}^{corr}");
485     fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->GetZaxis()->SetTitle("counts");
486     fOutput->Add(fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt);
487
488     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt = new TH2F("fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", "fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt", 
489                                                           2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, -5, 5);
490     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,2}^{corr}");  
491     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("#deltap_{T}^{corr} / p_{T,2}^{corr}");
492     fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
493     fOutput->Add(fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt);
494
495     fHistDeltaCorrPtvsDistance = new TH2F("fHistDeltaCorrPtvsDistance", "fHistDeltaCorrPtvsDistance", 
496                                           fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
497     fHistDeltaCorrPtvsDistance->GetXaxis()->SetTitle("Distance");  
498     fHistDeltaCorrPtvsDistance->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
499     fHistDeltaCorrPtvsDistance->GetZaxis()->SetTitle("counts");
500     fOutput->Add(fHistDeltaCorrPtvsDistance);
501
502     fHistDeltaCorrPtvsCommonEnergy1 = new TH2F("fHistDeltaCorrPtvsCommonEnergy1", "fHistDeltaCorrPtvsCommonEnergy1", 
503                                                fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
504     fHistDeltaCorrPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");  
505     fHistDeltaCorrPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
506     fHistDeltaCorrPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
507     fOutput->Add(fHistDeltaCorrPtvsCommonEnergy1);
508
509     fHistDeltaCorrPtvsCommonEnergy2 = new TH2F("fHistDeltaCorrPtvsCommonEnergy2", "fHistDeltaCorrPtvsCommonEnergy2", 
510                                                fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
511     fHistDeltaCorrPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");  
512     fHistDeltaCorrPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
513     fHistDeltaCorrPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
514     fOutput->Add(fHistDeltaCorrPtvsCommonEnergy2);
515
516     fHistDeltaCorrPtvsArea1 = new TH2F("fHistDeltaCorrPtvsArea1", "fHistDeltaCorrPtvsArea1", 
517                                        fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
518     fHistDeltaCorrPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
519     fHistDeltaCorrPtvsArea1->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
520     fHistDeltaCorrPtvsArea1->GetZaxis()->SetTitle("counts");
521     fOutput->Add(fHistDeltaCorrPtvsArea1);
522     
523     fHistDeltaCorrPtvsArea2 = new TH2F("fHistDeltaCorrPtvsArea2", "fHistDeltaCorrPtvsArea2", 
524                                        fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
525     fHistDeltaCorrPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
526     fHistDeltaCorrPtvsArea2->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
527     fHistDeltaCorrPtvsArea2->GetZaxis()->SetTitle("counts");
528     fOutput->Add(fHistDeltaCorrPtvsArea2);
529     
530     fHistDeltaCorrPtvsDeltaArea = new TH2F("fHistDeltaCorrPtvsDeltaArea", "fHistDeltaCorrPtvsDeltaArea", 
531                                            fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
532     fHistDeltaCorrPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
533     fHistDeltaCorrPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T}^{corr} (GeV/c)");
534     fHistDeltaCorrPtvsDeltaArea->GetZaxis()->SetTitle("counts");
535     fOutput->Add(fHistDeltaCorrPtvsDeltaArea);
536     
537     if (!fIsJet1Rho) 
538       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 
539                                              fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
540     else if (!fIsJet2Rho) 
541       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 
542                                              2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
543     else
544       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 
545                                              2*fNbins, -fMaxBinPt, fMaxBinPt, 
546                                              2*fNbins, -fMaxBinPt, fMaxBinPt);
547     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
548     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
549     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
550     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
551   }
552
553   if (fIsEmbedded) {
554     fHistDeltaMCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtvsJet1MCPt", "fHistDeltaMCPtvsJet1MCPt", 
555                                         fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
556     fHistDeltaMCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");  
557     fHistDeltaMCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
558     fHistDeltaMCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
559     fOutput->Add(fHistDeltaMCPtvsJet1MCPt);
560
561     fHistDeltaMCPtvsJet2Pt = new TH2F("fHistDeltaMCPtvsJet2Pt", "fHistDeltaMCPtvsJet2Pt", 
562                                       fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
563     fHistDeltaMCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
564     fHistDeltaMCPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
565     fHistDeltaMCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
566     fOutput->Add(fHistDeltaMCPtvsJet2Pt);
567
568     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt = new TH2F("fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", "fHistDeltaMCPtOverJet1MCPtvsJet1MCPt", 
569                                                     fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
570     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetXaxis()->SetTitle("p_{T,1}^{MC}");  
571     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,1}^{MC}");
572     fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->GetZaxis()->SetTitle("counts");
573     fOutput->Add(fHistDeltaMCPtOverJet1MCPtvsJet1MCPt);
574
575     fHistDeltaMCPtOverJet2PtvsJet2Pt = new TH2F("fHistDeltaMCPtOverJet2PtvsJet2Pt", "fHistDeltaMCPtOverJet2PtvsJet2Pt", 
576                                                 fNbins, fMinBinPt, fMaxBinPt, fNbins, -5, 5);
577     fHistDeltaMCPtOverJet2PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
578     fHistDeltaMCPtOverJet2PtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} / p_{T,2}");
579     fHistDeltaMCPtOverJet2PtvsJet2Pt->GetZaxis()->SetTitle("counts");
580     fOutput->Add(fHistDeltaMCPtOverJet2PtvsJet2Pt);
581
582     fHistDeltaMCPtvsDistance = new TH2F("fHistDeltaMCPtvsDistance", "fHistDeltaMCPtvsDistance", 
583                                         fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
584     fHistDeltaMCPtvsDistance->GetXaxis()->SetTitle("Distance");  
585     fHistDeltaMCPtvsDistance->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
586     fHistDeltaMCPtvsDistance->GetZaxis()->SetTitle("counts");
587     fOutput->Add(fHistDeltaMCPtvsDistance);
588
589     fHistDeltaMCPtvsCommonEnergy1 = new TH2F("fHistDeltaMCPtvsCommonEnergy1", "fHistDeltaMCPtvsCommonEnergy1", 
590                                              fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
591     fHistDeltaMCPtvsCommonEnergy1->GetXaxis()->SetTitle("Common energy 1 (%)");  
592     fHistDeltaMCPtvsCommonEnergy1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
593     fHistDeltaMCPtvsCommonEnergy1->GetZaxis()->SetTitle("counts");
594     fOutput->Add(fHistDeltaMCPtvsCommonEnergy1);
595
596     fHistDeltaMCPtvsCommonEnergy2 = new TH2F("fHistDeltaMCPtvsCommonEnergy2", "fHistDeltaMCPtvsCommonEnergy2", 
597                                              fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
598     fHistDeltaMCPtvsCommonEnergy2->GetXaxis()->SetTitle("Common energy 2 (%)");  
599     fHistDeltaMCPtvsCommonEnergy2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
600     fHistDeltaMCPtvsCommonEnergy2->GetZaxis()->SetTitle("counts");
601     fOutput->Add(fHistDeltaMCPtvsCommonEnergy2);
602
603     fHistDeltaMCPtvsArea1 = new TH2F("fHistDeltaMCPtvsArea1", "fHistDeltaMCPtvsArea1", 
604                                      fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
605     fHistDeltaMCPtvsArea1->GetXaxis()->SetTitle("A_{jet,1}");
606     fHistDeltaMCPtvsArea1->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
607     fHistDeltaMCPtvsArea1->GetZaxis()->SetTitle("counts");
608     fOutput->Add(fHistDeltaMCPtvsArea1);
609
610     fHistDeltaMCPtvsArea2 = new TH2F("fHistDeltaMCPtvsArea2", "fHistDeltaMCPtvsArea2", 
611                                      fNbins/2, 0, 1.5, 2*fNbins, -fMaxBinPt, fMaxBinPt);
612     fHistDeltaMCPtvsArea2->GetXaxis()->SetTitle("A_{jet,2}");
613     fHistDeltaMCPtvsArea2->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
614     fHistDeltaMCPtvsArea2->GetZaxis()->SetTitle("counts");
615     fOutput->Add(fHistDeltaMCPtvsArea2);
616
617     fHistDeltaMCPtvsDeltaArea = new TH2F("fHistDeltaMCPtvsDeltaArea", "fHistDeltaMCPtvsDeltaArea", 
618                                          fNbins, -1.+1./fNbins, 1.+1./fNbins, 2*fNbins, -fMaxBinPt, fMaxBinPt);
619     fHistDeltaMCPtvsDeltaArea->GetXaxis()->SetTitle("#deltaA_{jet}");
620     fHistDeltaMCPtvsDeltaArea->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
621     fHistDeltaMCPtvsDeltaArea->GetZaxis()->SetTitle("counts");
622     fOutput->Add(fHistDeltaMCPtvsDeltaArea);
623
624     fHistJet1MCPtvsJet2Pt = new TH2F("fHistJet1MCPtvsJet2Pt", "fHistJet1MCPtvsJet2Pt", 
625                                      fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
626     fHistJet1MCPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}^{MC}");
627     fHistJet1MCPtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
628     fHistJet1MCPtvsJet2Pt->GetZaxis()->SetTitle("counts");
629     fOutput->Add(fHistJet1MCPtvsJet2Pt);
630   }
631 }
632
633 //________________________________________________________________________
634 void AliJetResponseMaker::AllocateTHnSparse()
635 {
636   // Allocate THnSparse histograms.
637
638   TString title[20]= {""};
639   Int_t nbins[20]  = {0};
640   Double_t min[20] = {0.};
641   Double_t max[20] = {0.};
642   Int_t dim = 0;
643
644   title[dim] = "#phi";
645   nbins[dim] = fNbins/4;
646   min[dim] = 0;
647   max[dim] = 2*TMath::Pi()*(1 + 1./(nbins[dim]-1));
648   dim++;
649
650   title[dim] = "#eta";
651   nbins[dim] = fNbins/4;
652   min[dim] = -1;
653   max[dim] = 1;
654   dim++;
655
656   title[dim] = "p_{T}";
657   nbins[dim] = fNbins;
658   min[dim] = 0;
659   max[dim] = 250;
660   dim++;
661
662   title[dim] = "A_{jet}";
663   nbins[dim] = fNbins/4;
664   min[dim] = 0;
665   max[dim] = 1.5;
666   dim++;
667
668   title[dim] = "NEF";
669   nbins[dim] = fNbins/4;
670   min[dim] = 0;
671   max[dim] = 1.2;
672   dim++;
673
674   title[dim] = "Z";
675   nbins[dim] = fNbins/4;
676   min[dim] = 0;
677   max[dim] = 1.2;
678   dim++;
679
680   Int_t dim1 = dim, dim2 = dim;
681
682   if (fIsJet1Rho) {
683     title[dim1] = "p_{T}^{corr}";
684     nbins[dim1] = fNbins*2;
685     min[dim1] = -250;
686     max[dim1] = 250;
687     dim1++;
688   }
689
690   if (fIsEmbedded) {
691     title[dim1] = "p_{T}^{MC}";
692     nbins[dim1] = fNbins;
693     min[dim1] = 0;
694     max[dim1] = 250;
695     dim1++;
696   }
697
698   fHistJets1 = new THnSparseD("fHistJets1","fHistJets1",dim1,nbins,min,max);
699   for (Int_t i = 0; i < dim1; i++) 
700     fHistJets1->GetAxis(i)->SetTitle(title[i]);
701   fOutput->Add(fHistJets1);
702
703   if (fIsJet2Rho) {
704     title[dim2] = "p_{T}^{corr}";
705     nbins[dim2] = fNbins*2;
706     min[dim2] = -250;
707     max[dim2] = 250;
708     dim2++;
709   }
710
711   if (fDoJet2Histogram) {
712     fHistJets2 = new THnSparseD("fHistJets2","fHistJets2",dim2,nbins,min,max);
713     for (Int_t i = 0; i < dim2; i++) 
714       fHistJets2->GetAxis(i)->SetTitle(title[i]);
715     fOutput->Add(fHistJets2);
716   }
717
718   fHistJets2Acceptance = new THnSparseD("fHistJets2Acceptance","fHistJets2Acceptance",dim2,nbins,min,max);
719   for (Int_t i = 0; i < dim2; i++) 
720     fHistJets2Acceptance->GetAxis(i)->SetTitle(title[i]);
721   fOutput->Add(fHistJets2Acceptance);
722   
723   // Matching
724
725   dim = 0;
726
727   title[dim] = "p_{T,1}";
728   nbins[dim] = fNbins;
729   min[dim] = 0;
730   max[dim] = 250;
731   dim++;
732
733   title[dim] = "p_{T,2}";
734   nbins[dim] = fNbins;
735   min[dim] = 0;
736   max[dim] = 250;
737   dim++;
738
739   title[dim] = "A_{jet,1}";
740   nbins[dim] = fNbins/4;
741   min[dim] = 0;
742   max[dim] = 1.5;
743   dim++;
744
745   title[dim] = "A_{jet,2}";
746   nbins[dim] = fNbins/4;
747   min[dim] = 0;
748   max[dim] = 1.5;
749   dim++;
750
751   title[dim] = "distance";
752   nbins[dim] = fNbins/2;
753   min[dim] = 0;
754   max[dim] = 1.2;
755   dim++;
756
757   title[dim] = "CE1";
758   nbins[dim] = fNbins/2;
759   min[dim] = 0;
760   max[dim] = 1.2;
761   dim++;
762
763   title[dim] = "CE2";
764   nbins[dim] = fNbins/2;
765   min[dim] = 0;
766   max[dim] = 1.2;
767   dim++;
768
769   if (fDeltaPtAxis) {
770     title[dim] = "#deltaA_{jet}";
771     nbins[dim] = fNbins/2;
772     min[dim] = -1;
773     max[dim] = 1;
774     dim++;
775
776     title[dim] = "#deltap_{T}";
777     nbins[dim] = fNbins*2;
778     min[dim] = -250;
779     max[dim] = 250;
780     dim++;
781   }
782   if (fIsJet1Rho) {
783     title[dim] = "p_{T,1}^{corr}";
784     nbins[dim] = fNbins*2;
785     min[dim] = -250;
786     max[dim] = 250;
787     dim++;
788   }
789   if (fIsJet2Rho) {
790     title[dim] = "p_{T,2}^{corr}";
791     nbins[dim] = fNbins*2;
792     min[dim] = -250;
793     max[dim] = 250;
794     dim++;
795   }
796   if (fDeltaPtAxis && (fIsJet1Rho || fIsJet2Rho)) {
797     title[dim] = "#deltap_{T}^{corr}";
798     nbins[dim] = fNbins*2;
799     min[dim] = -250;
800     max[dim] = 250;
801     dim++;
802   }
803   if (fDeltaEtaDeltaPhiAxis) {
804     title[dim] = "#delta#eta";
805     nbins[dim] = fNbins/2;
806     min[dim] = -1;
807     max[dim] = 1;
808     dim++;
809
810     title[dim] = "#delta#phi";
811     nbins[dim] = fNbins/2;
812     min[dim] = -TMath::Pi()/2;
813     max[dim] = TMath::Pi()*3/2;
814     dim++;
815   }
816   if (fIsEmbedded) {
817     title[dim] = "p_{T,1}^{MC}";
818     nbins[dim] = fNbins;
819     min[dim] = 0;
820     max[dim] = 250;
821     dim++;
822
823     if (fDeltaPtAxis) {
824       title[dim] = "#deltap_{T}^{MC}";
825       nbins[dim] = fNbins*2;
826       min[dim] = -250;
827       max[dim] = 250;
828       dim++;
829     }
830   }
831
832   if (fNEFAxis) {
833     title[dim] = "NEF_{1}";
834     nbins[dim] = fNbins/4;
835     min[dim] = 0;
836     max[dim] = 1.2;
837     dim++;
838
839     title[dim] = "NEF_{2}";
840     nbins[dim] = fNbins/4;
841     min[dim] = 0;
842     max[dim] = 1.2;
843     dim++;
844   }
845
846   if (fZAxis) {
847     title[dim] = "Z_{1}";
848     nbins[dim] = fNbins/4;
849     min[dim] = 0;
850     max[dim] = 1.2;
851     dim++;
852
853     title[dim] = "Z_{2}";
854     nbins[dim] = fNbins/4;
855     min[dim] = 0;
856     max[dim] = 1.2;
857     dim++;
858   }
859       
860   fHistMatching = new THnSparseD("fHistMatching","fHistMatching",dim,nbins,min,max);
861     
862   for (Int_t i = 0; i < dim; i++)
863     fHistMatching->GetAxis(i)->SetTitle(title[i]);
864    
865   fOutput->Add(fHistMatching);
866 }
867
868 //________________________________________________________________________
869 void AliJetResponseMaker::UserCreateOutputObjects()
870 {
871   // Create user objects.
872
873   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
874
875   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
876   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
877
878   if (!jets1 || !jets2) return;
879
880   if (jets1->GetRhoName().IsNull()) fIsJet1Rho = kFALSE;
881   else fIsJet1Rho = kTRUE;
882   if (jets2->GetRhoName().IsNull()) fIsJet2Rho = kFALSE;
883   else fIsJet2Rho = kTRUE;
884
885   if (fHistoType==0)
886     AllocateTH2();
887   else 
888     AllocateTHnSparse();
889
890   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
891 }
892
893 //________________________________________________________________________
894 void AliJetResponseMaker::FillJetHisto(Double_t Phi, Double_t Eta, Double_t Pt, Double_t A, Double_t NEF, Double_t Z, Double_t CorrPt, Double_t MCPt, Int_t Set)
895 {
896   if (fHistoType==1) {
897     THnSparse *histo = 0;
898     if (Set==1)
899       histo = fHistJets1;
900     else if (Set==2)
901       histo = fHistJets2;
902     else if (Set==3)
903       histo = fHistJets2Acceptance;
904
905     if (!histo)
906       return;
907
908     Double_t contents[20]={0};
909
910     for (Int_t i = 0; i < histo->GetNdimensions(); i++) {
911       TString title(histo->GetAxis(i)->GetTitle());
912       if (title=="#phi")
913         contents[i] = Phi;
914       else if (title=="#eta")
915         contents[i] = Eta;
916       else if (title=="p_{T}")
917         contents[i] = Pt;
918       else if (title=="A_{jet}")
919         contents[i] = A;
920       else if (title=="NEF")
921         contents[i] = NEF;
922       else if (title=="Z")
923         contents[i] = Z;
924       else if (title=="p_{T}^{corr}")
925         contents[i] = CorrPt;
926       else if (title=="p_{T}^{MC}")
927         contents[i] = MCPt;
928       else 
929         AliWarning(Form("Unable to fill dimension %s!",title.Data()));
930     }
931
932     histo->Fill(contents);
933   }
934   else {
935     if (Set == 1) {
936       fHistJets1PtArea->Fill(A, Pt);
937       fHistJets1PhiEta->Fill(Eta, Phi);
938       fHistJets1ZvsPt->Fill(Z, Pt);
939       fHistJets1NEFvsPt->Fill(NEF, Pt);
940       fHistJets1CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
941       if (fIsJet1Rho)
942         fHistJets1CorrPtArea->Fill(A, CorrPt);
943     }
944     else if (Set == 2) {
945       fHistJets2PtArea->Fill(A, Pt);
946       fHistJets2PhiEta->Fill(Eta, Phi);
947       if (fIsJet2Rho)
948         fHistJets2CorrPtArea->Fill(A, CorrPt);
949     }
950     else if (Set == 3) {
951       fHistJets2PtAreaAcceptance->Fill(A, Pt);
952       fHistJets2PhiEtaAcceptance->Fill(Eta, Phi);
953       fHistJets2ZvsPt->Fill(Z, Pt);
954       fHistJets2NEFvsPt->Fill(NEF, Pt);
955       fHistJets2CEFvsCEFPt->Fill(1-NEF, (1-NEF)*Pt);
956       if (fIsJet2Rho)
957         fHistJets2CorrPtAreaAcceptance->Fill(A, CorrPt);
958     }
959   }
960 }
961
962 //________________________________________________________________________
963 void AliJetResponseMaker::FillMatchingHistos(Double_t Pt1, Double_t Pt2, Double_t Eta1, Double_t Eta2, Double_t Phi1, Double_t Phi2, 
964                                              Double_t A1, Double_t A2, Double_t d, Double_t CE1, Double_t CE2, Double_t CorrPt1, Double_t CorrPt2, 
965                                              Double_t MCPt1, Double_t NEF1, Double_t NEF2, Double_t Z1, Double_t Z2)
966 {
967   if (fHistoType==1) {
968     Double_t contents[20]={0};
969
970     for (Int_t i = 0; i < fHistMatching->GetNdimensions(); i++) {
971       TString title(fHistMatching->GetAxis(i)->GetTitle());
972       if (title=="p_{T,1}")
973         contents[i] = Pt1;
974       else if (title=="p_{T,2}")
975         contents[i] = Pt2;
976       else if (title=="A_{jet,1}")
977         contents[i] = A1;
978       else if (title=="A_{jet,2}")
979         contents[i] = A2;
980       else if (title=="distance")
981         contents[i] = d;
982       else if (title=="CE1")
983         contents[i] = CE1;
984       else if (title=="CE2")
985         contents[i] = CE2;
986       else if (title=="#deltaA_{jet}")
987         contents[i] = A1-A2;
988       else if (title=="#deltap_{T}")
989         contents[i] = Pt1-Pt2;
990       else if (title=="#delta#eta")
991         contents[i] = Eta1-Eta2;
992       else if (title=="#delta#phi")
993         contents[i] = Phi1-Phi2;
994       else if (title=="p_{T,1}^{corr}")
995         contents[i] = CorrPt1;
996       else if (title=="p_{T,2}^{corr}")
997         contents[i] = CorrPt2;
998       else if (title=="#deltap_{T}^{corr}")
999         contents[i] = CorrPt1-CorrPt2;
1000       else if (title=="p_{T,1}^{MC}")
1001         contents[i] = MCPt1;
1002       else if (title=="#deltap_{T}^{MC}")
1003         contents[i] = MCPt1-Pt2;
1004       else if (title=="NEF_{1}")
1005         contents[i] = NEF1;
1006       else if (title=="NEF_{2}")
1007         contents[i] = NEF2;
1008       else if (title=="Z_{1}")
1009         contents[i] = Z1;
1010       else if (title=="Z_{2}")
1011         contents[i] = Z2;
1012       else 
1013         AliWarning(Form("Unable to fill dimension %s!",title.Data()));
1014     }
1015
1016     fHistMatching->Fill(contents);
1017   }
1018   else {
1019     fHistCommonEnergy1vsJet1Pt->Fill(CE1, Pt1);
1020     fHistCommonEnergy2vsJet2Pt->Fill(CE2, Pt2);
1021
1022     fHistDistancevsJet1Pt->Fill(d, Pt1);
1023     fHistDistancevsJet2Pt->Fill(d, Pt2);
1024
1025     fHistDistancevsCommonEnergy1->Fill(d, CE1);
1026     fHistDistancevsCommonEnergy2->Fill(d, CE2);
1027     fHistCommonEnergy1vsCommonEnergy2->Fill(CE1, CE2);
1028
1029     fHistDeltaEtaDeltaPhi->Fill(Eta1-Eta2,Phi1-Phi2);
1030
1031     fHistJet2PtOverJet1PtvsJet2Pt->Fill(Pt2, Pt2 / Pt1);
1032     fHistJet1PtOverJet2PtvsJet1Pt->Fill(Pt1, Pt1 / Pt2);
1033
1034     Double_t dpt = Pt1 - Pt2;
1035     fHistDeltaPtvsJet1Pt->Fill(Pt1, dpt);
1036     fHistDeltaPtvsJet2Pt->Fill(Pt2, dpt);
1037     fHistDeltaPtOverJet1PtvsJet1Pt->Fill(Pt1, dpt/Pt1);
1038     fHistDeltaPtOverJet2PtvsJet2Pt->Fill(Pt2, dpt/Pt2);
1039
1040     fHistDeltaPtvsDistance->Fill(d, dpt);
1041     fHistDeltaPtvsCommonEnergy1->Fill(CE1, dpt);
1042     fHistDeltaPtvsCommonEnergy2->Fill(CE2, dpt);
1043
1044     fHistDeltaPtvsArea1->Fill(A1, dpt);
1045     fHistDeltaPtvsArea2->Fill(A2, dpt);
1046
1047     Double_t darea = A1 - A2;
1048     fHistDeltaPtvsDeltaArea->Fill(darea, dpt);
1049
1050     fHistJet1PtvsJet2Pt->Fill(Pt1, Pt2);
1051
1052     if (fIsJet1Rho || fIsJet2Rho) {
1053       Double_t dcorrpt = CorrPt1 - CorrPt2;
1054       fHistDeltaCorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt);
1055       fHistDeltaCorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt);
1056       fHistDeltaCorrPtOverJet1CorrPtvsJet1CorrPt->Fill(CorrPt1, dcorrpt/CorrPt1);
1057       fHistDeltaCorrPtOverJet2CorrPtvsJet2CorrPt->Fill(CorrPt2, dcorrpt/CorrPt2);
1058       fHistDeltaCorrPtvsDistance->Fill(d, dcorrpt);
1059       fHistDeltaCorrPtvsCommonEnergy1->Fill(CE1, dcorrpt);
1060       fHistDeltaCorrPtvsCommonEnergy2->Fill(CE2, dcorrpt);
1061       fHistDeltaCorrPtvsArea1->Fill(A1, dcorrpt);
1062       fHistDeltaCorrPtvsArea2->Fill(A2, dcorrpt);
1063       fHistDeltaCorrPtvsDeltaArea->Fill(darea, dcorrpt);
1064       fHistJet1CorrPtvsJet2CorrPt->Fill(CorrPt1, CorrPt2);
1065     }
1066
1067     if (fIsEmbedded) {
1068       Double_t dmcpt = MCPt1 - Pt2;
1069       fHistDeltaMCPtvsJet1MCPt->Fill(MCPt1, dmcpt);
1070       fHistDeltaMCPtvsJet2Pt->Fill(Pt2, dmcpt);
1071       fHistDeltaMCPtOverJet1MCPtvsJet1MCPt->Fill(MCPt1, dmcpt/MCPt1);
1072       fHistDeltaMCPtOverJet2PtvsJet2Pt->Fill(Pt2, dmcpt/Pt2);
1073       fHistDeltaMCPtvsDistance->Fill(d, dmcpt);
1074       fHistDeltaMCPtvsCommonEnergy1->Fill(CE1, dmcpt);
1075       fHistDeltaMCPtvsCommonEnergy2->Fill(CE2, dmcpt);
1076       fHistDeltaMCPtvsArea1->Fill(A1, dmcpt);
1077       fHistDeltaMCPtvsArea2->Fill(A2, dmcpt);
1078       fHistDeltaMCPtvsDeltaArea->Fill(darea, dmcpt);
1079       fHistJet1MCPtvsJet2Pt->Fill(MCPt1, Pt2);
1080     }
1081   }
1082 }
1083
1084 //________________________________________________________________________
1085 void AliJetResponseMaker::ExecOnce()
1086 {
1087   // Execute once.
1088
1089   AliAnalysisTaskEmcalJet::ExecOnce();
1090
1091   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1092   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1093
1094   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1095
1096   if (fMatching == kMCLabel && (!jets2->GetIsParticleLevel() || jets1->GetIsParticleLevel())) {
1097     if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel()) {
1098       AliWarning("Changing matching type from MC label to same collection...");
1099       fMatching = kSameCollections;
1100     }
1101     else {
1102       AliWarning("Changing matching type from MC label to geometrical...");
1103       fMatching = kGeometrical;
1104     }
1105   }
1106   else if (fMatching == kSameCollections && jets1->GetIsParticleLevel() != jets2->GetIsParticleLevel()) {
1107     if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) {
1108       AliWarning("Changing matching type from same collection to MC label...");
1109       fMatching = kMCLabel;
1110     }
1111     else {
1112       AliWarning("Changing matching type from same collection to geometrical...");
1113       fMatching = kGeometrical;
1114     }
1115   }
1116 }
1117
1118 //________________________________________________________________________
1119 Bool_t AliJetResponseMaker::Run()
1120 {
1121   // Find the closest jets
1122
1123   if (fMatching == kNoMatching) 
1124     return kTRUE;
1125   else
1126     return DoJetMatching();
1127 }
1128
1129 //________________________________________________________________________
1130 Bool_t AliJetResponseMaker::DoJetMatching()
1131 {
1132   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1133   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1134
1135   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1136
1137   DoJetLoop();
1138
1139   AliEmcalJet* jet1 = 0;
1140
1141   jets1->ResetCurrentID();
1142   while ((jet1 = jets1->GetNextJet())) {
1143
1144     AliEmcalJet *jet2 = jet1->ClosestJet();
1145
1146     if (!jet2) continue;
1147     if (jet2->ClosestJet() != jet1) continue;
1148     if (jet1->ClosestJetDistance() > fMatchingPar1 || jet2->ClosestJetDistance() > fMatchingPar2) continue;
1149
1150     // Matched jet found
1151     jet1->SetMatchedToClosest(fMatching);
1152     jet2->SetMatchedToClosest(fMatching);
1153     AliDebug(2,Form("Found matching: jet1 pt = %f, eta = %f, phi = %f, jet2 pt = %f, eta = %f, phi = %f", 
1154                     jet1->Pt(), jet1->Eta(), jet1->Phi(), 
1155                     jet2->Pt(), jet2->Eta(), jet2->Phi()));
1156   }
1157
1158   return kTRUE;
1159 }
1160
1161 //________________________________________________________________________
1162 void AliJetResponseMaker::DoJetLoop()
1163 {
1164   // Do the jet loop.
1165
1166   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1167   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1168
1169   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1170
1171   AliEmcalJet* jet1 = 0;
1172   AliEmcalJet* jet2 = 0;
1173
1174   jets2->ResetCurrentID();
1175   while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
1176     
1177   jets1->ResetCurrentID();
1178   while ((jet1 = jets1->GetNextJet())) {
1179     jet1->ResetMatching();
1180     
1181     if (jet1->MCPt() < fMinJetMCPt) continue;
1182
1183     jets2->ResetCurrentID();
1184     while ((jet2 = jets2->GetNextJet())) {
1185       SetMatchingLevel(jet1, jet2, fMatching);
1186     } // jet2 loop
1187   } // jet1 loop
1188 }
1189
1190 //________________________________________________________________________
1191 void AliJetResponseMaker::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
1192 {
1193   Double_t deta = jet2->Eta() - jet1->Eta();
1194   Double_t dphi = jet2->Phi() - jet1->Phi();
1195   d = TMath::Sqrt(deta * deta + dphi * dphi);
1196 }
1197
1198 //________________________________________________________________________
1199 void AliJetResponseMaker::GetMCLabelMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1200
1201   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1202   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1203
1204   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1205
1206   AliParticleContainer *tracks1   = jets1->GetParticleContainer();
1207   AliClusterContainer  *clusters1 = jets1->GetClusterContainer();
1208   AliParticleContainer *tracks2   = jets2->GetParticleContainer();
1209
1210   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1211   d1 = jet1->Pt();
1212   d2 = jet2->Pt();
1213   Double_t totalPt1 = d1; // the total pt of the reconstructed jet will be cleaned from the background
1214
1215   // remove completely tracks that are not MC particles (label == 0)
1216   if (tracks1 && tracks1->GetArray()) {
1217     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1218       AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1219       if (!track) {
1220         AliWarning(Form("Could not find track %d!", iTrack));
1221         continue;
1222       }
1223
1224       Int_t MClabel = TMath::Abs(track->GetLabel());
1225       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1226       if (MClabel != 0) continue;
1227
1228       // this is not a MC particle; remove it completely
1229       AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1230       totalPt1 -= track->Pt();
1231       d1 -= track->Pt();
1232     }
1233   }
1234
1235   // remove completely clusters that are not MC particles (label == 0)
1236   if (fUseCellsToMatch && fCaloCells) { 
1237     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1238       AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1239       if (!clus) {
1240         AliWarning(Form("Could not find cluster %d!", iClus));
1241         continue;
1242       }
1243       TLorentzVector part;
1244       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1245           
1246       for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1247         Int_t cellId = clus->GetCellAbsId(iCell);
1248         Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1249
1250         Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1251         if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1252         if (MClabel != 0) continue;
1253
1254         // this is not a MC particle; remove it completely
1255         AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1256         totalPt1 -= part.Pt() * cellFrac;
1257         d1 -= part.Pt() * cellFrac;
1258       }
1259     }
1260   }
1261   else {
1262     for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1263       AliVCluster *clus = jet1->ClusterAt(iClus,clusters1->GetArray());
1264       if (!clus) {
1265         AliWarning(Form("Could not find cluster %d!", iClus));
1266         continue;
1267       }
1268       TLorentzVector part;
1269       clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1270           
1271       Int_t MClabel = TMath::Abs(clus->GetLabel());
1272       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1273       if (MClabel != 0) continue;
1274
1275       // this is not a MC particle; remove it completely
1276       AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1277       totalPt1 -= part.Pt();
1278       d1 -= part.Pt();
1279     }
1280   }
1281
1282   for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1283     Bool_t track2Found = kFALSE;
1284     Int_t index2 = jet2->TrackAt(iTrack2);
1285
1286     // now look for common particles in the track array
1287     for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
1288       AliVParticle *track = jet1->TrackAt(iTrack,tracks1->GetArray());
1289       if (!track) {
1290         AliWarning(Form("Could not find track %d!", iTrack));
1291         continue;
1292       }
1293       Int_t MClabel = TMath::Abs(track->GetLabel());
1294       if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;      
1295       if (MClabel <= 0) continue;
1296
1297       Int_t index = -1;
1298       index = tracks2->GetIndexFromLabel(MClabel);
1299       if (index < 0) {
1300         AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
1301         continue;
1302       }
1303         
1304       if (index2 != index) continue;
1305       
1306       // found common particle
1307       d1 -= track->Pt();
1308
1309       if (!track2Found) {
1310         AliVParticle *MCpart = tracks2->GetParticle(index2);
1311         AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1312                         iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1313         d2 -= MCpart->Pt();
1314       }
1315
1316       track2Found = kTRUE;
1317     }
1318
1319     // now look for common particles in the cluster array
1320     if (fUseCellsToMatch && fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
1321       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1322         AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1323         if (!clus) {
1324           AliWarning(Form("Could not find cluster %d!", iClus));
1325           continue;
1326         }
1327         TLorentzVector part;
1328         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1329       
1330         for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
1331           Int_t cellId = clus->GetCellAbsId(iCell);
1332           Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
1333         
1334           Int_t MClabel = TMath::Abs(fCaloCells->GetCellMCLabel(cellId));
1335           if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;
1336           if (MClabel <= 0) continue;
1337           
1338           Int_t index1 = -1;
1339           index1 = tracks2->GetIndexFromLabel(MClabel);
1340           if (index1 < 0) {
1341             AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
1342             continue;
1343           }
1344         
1345           if (index2 != index1) continue;
1346           
1347           // found common particle
1348           d1 -= part.Pt() * cellFrac;
1349         
1350           if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1351             AliVParticle *MCpart = tracks2->GetParticle(index2);
1352             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)!",
1353                             iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));               
1354             d2 -= MCpart->Pt() * cellFrac;
1355           }
1356
1357           track2Found = kTRUE;
1358         }
1359       }
1360     }
1361     else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
1362       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
1363         AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
1364         if (!clus) {
1365           AliWarning(Form("Could not find cluster %d!", iClus));
1366           continue;
1367         }
1368         TLorentzVector part;
1369         clus->GetMomentum(part, const_cast<Double_t*>(fVertex));
1370       
1371         Int_t MClabel = TMath::Abs(clus->GetLabel());
1372         if (MClabel > fMCLabelShift) MClabel -= fMCLabelShift;      
1373         if (MClabel <= 0) continue;
1374
1375         Int_t index = -1;
1376         index = tracks2->GetIndexFromLabel(MClabel);
1377       
1378         if (index < 0) {
1379           AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
1380           continue;
1381         }
1382
1383         if (index2 != index) continue;
1384       
1385         // found common particle
1386         d1 -= part.Pt();
1387       
1388         if (!track2Found) { // only if it is not already found among charged tracks (charged particles are most likely already found)
1389           AliVParticle *MCpart = tracks2->GetParticle(index2);
1390           AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
1391                           iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
1392           
1393           d2 -= MCpart->Pt();
1394         }
1395
1396         track2Found = kTRUE;
1397       }
1398     }
1399   }
1400
1401   if (d1 < 0)
1402     d1 = 0;
1403
1404   if (d2 < 0)
1405     d2 = 0;
1406
1407   if (totalPt1 < 1)
1408     d1 = -1;
1409   else
1410     d1 /= totalPt1;
1411
1412   if (jet2->Pt() < 1)
1413     d2 = -1;
1414   else
1415     d2 /= jet2->Pt();
1416 }
1417
1418 //________________________________________________________________________
1419 void AliJetResponseMaker::GetSameCollectionsMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d1, Double_t &d2) const
1420
1421   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1422   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1423
1424   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
1425
1426   AliParticleContainer *tracks1   = jets1->GetParticleContainer();
1427   AliClusterContainer  *clusters1 = jets1->GetClusterContainer();
1428   AliParticleContainer *tracks2   = jets2->GetParticleContainer();
1429   AliClusterContainer  *clusters2 = jets2->GetClusterContainer();
1430
1431   // d1 and d2 represent the matching level: 0 = maximum level of matching, 1 = the two jets are completely unrelated
1432   d1 = jet1->Pt();
1433   d2 = jet2->Pt();
1434
1435   if (tracks1 && tracks2) {
1436
1437     for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
1438       Int_t index2 = jet2->TrackAt(iTrack2);
1439       for (Int_t iTrack1 = 0; iTrack1 < jet1->GetNumberOfTracks(); iTrack1++) {
1440         Int_t index1 = jet1->TrackAt(iTrack1);
1441         if (index2 == index1) { // found common particle
1442           AliVParticle *part1 = tracks1->GetParticle(index1);
1443           if (!part1) {
1444             AliWarning(Form("Could not find track %d!", index1));
1445             continue;
1446           }
1447           AliVParticle *part2 = tracks2->GetParticle(index2);
1448           if (!part2) {
1449             AliWarning(Form("Could not find track %d!", index2));
1450             continue;
1451           }
1452
1453           d1 -= part1->Pt();
1454           d2 -= part2->Pt();
1455           break;
1456         }
1457       }
1458     }
1459
1460   }
1461
1462   if (clusters1 && clusters2) {
1463
1464     if (fUseCellsToMatch) {
1465       const Int_t nClus1 = jet1->GetNumberOfClusters();
1466
1467       Int_t ncells1[nClus1];
1468       UShort_t *cellsId1[nClus1];
1469       Double_t *cellsFrac1[nClus1];
1470       Int_t *sortedIndexes1[nClus1];
1471       Double_t ptClus1[nClus1];
1472       for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1473         Int_t index1 = jet1->ClusterAt(iClus1);
1474         AliVCluster *clus1 = clusters1->GetCluster(index1);
1475         if (!clus1) {
1476           AliWarning(Form("Could not find cluster %d!", index1));
1477           ncells1[iClus1] = 0;
1478           cellsId1[iClus1] = 0;
1479           cellsFrac1[iClus1] = 0;
1480           sortedIndexes1[iClus1] = 0;
1481           ptClus1[iClus1] = 0;
1482           continue;
1483         }
1484         TLorentzVector part1;
1485         clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1486
1487         ncells1[iClus1] = clus1->GetNCells();
1488         cellsId1[iClus1] = clus1->GetCellsAbsId();
1489         cellsFrac1[iClus1] = clus1->GetCellsAmplitudeFraction();
1490         sortedIndexes1[iClus1] = new Int_t[ncells1[iClus1]];
1491         ptClus1[iClus1] = part1.Pt();
1492
1493         TMath::Sort(ncells1[iClus1], cellsId1[iClus1], sortedIndexes1[iClus1], kFALSE);
1494       }
1495       
1496       const Int_t nClus2 = jet2->GetNumberOfClusters();
1497
1498       const Int_t maxNcells2 = 11520;
1499       Int_t sortedIndexes2[maxNcells2];
1500       for (Int_t iClus2 = 0; iClus2 < nClus2; iClus2++) {
1501         Int_t index2 = jet2->ClusterAt(iClus2);
1502         AliVCluster *clus2 =  clusters2->GetCluster(index2);
1503         if (!clus2) {
1504           AliWarning(Form("Could not find cluster %d!", index2));
1505           continue;
1506         }
1507         Int_t ncells2 = clus2->GetNCells();
1508         if (ncells2 >= maxNcells2) {
1509           AliError(Form("Number of cells in the cluster %d >= %d",ncells2,maxNcells2));
1510           continue;
1511         }
1512         UShort_t *cellsId2 = clus2->GetCellsAbsId();
1513         Double_t *cellsFrac2 = clus2->GetCellsAmplitudeFraction();
1514
1515         TLorentzVector part2;
1516         clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1517         Double_t ptClus2 = part2.Pt();
1518
1519         TMath::Sort(ncells2, cellsId2, sortedIndexes2, kFALSE);
1520
1521         for (Int_t iClus1 = 0; iClus1 < nClus1; iClus1++) {
1522           if (sortedIndexes1[iClus1] == 0)
1523             continue;
1524           Int_t iCell1 = 0, iCell2 = 0;
1525           while (iCell1 < ncells1[iClus1] && iCell2 < ncells2) {
1526             if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] == cellsId2[sortedIndexes2[iCell2]]) { // found a common cell
1527               d1 -= cellsFrac1[iClus1][sortedIndexes1[iClus1][iCell1]] * ptClus1[iClus1];
1528               d2 -= cellsFrac2[sortedIndexes2[iCell2]] * ptClus2;
1529               iCell1++;
1530               iCell2++;
1531             }
1532             else if (cellsId1[iClus1][sortedIndexes1[iClus1][iCell1]] > cellsId2[sortedIndexes2[iCell2]]) { 
1533               iCell2++;
1534             }
1535             else {
1536               iCell1++;
1537             }
1538           }
1539         }
1540       }
1541     }
1542     else {
1543       for (Int_t iClus2 = 0; iClus2 < jet2->GetNumberOfClusters(); iClus2++) {
1544         Int_t index2 = jet2->ClusterAt(iClus2);
1545         for (Int_t iClus1 = 0; iClus1 < jet1->GetNumberOfClusters(); iClus1++) {
1546           Int_t index1 = jet1->ClusterAt(iClus1);
1547           if (index2 == index1) { // found common particle
1548             AliVCluster *clus1 = clusters1->GetCluster(index1);
1549             if (!clus1) {
1550               AliWarning(Form("Could not find cluster %d!", index1));
1551               continue;
1552             }
1553             AliVCluster *clus2 =  clusters2->GetCluster(index2);
1554             if (!clus2) {
1555               AliWarning(Form("Could not find cluster %d!", index2));
1556               continue;
1557             }
1558             TLorentzVector part1, part2;
1559             clus1->GetMomentum(part1, const_cast<Double_t*>(fVertex));
1560             clus2->GetMomentum(part2, const_cast<Double_t*>(fVertex));
1561             
1562             d1 -= part1.Pt();
1563             d2 -= part2.Pt();
1564             break;
1565           }
1566         }
1567       }
1568     }
1569   }
1570
1571   if (d1 < 0)
1572     d1 = 0;
1573
1574   if (d2 < 0)
1575     d2 = 0;
1576
1577   if (jet1->Pt() > 0)
1578     d1 /= jet1->Pt();
1579   else
1580     d1 = -1;
1581
1582   if (jet2->Pt() > 0)
1583     d2 /= jet2->Pt();
1584   else
1585     d2 = -1;
1586 }
1587
1588 //________________________________________________________________________
1589 void AliJetResponseMaker::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
1590 {
1591   Double_t d1 = -1;
1592   Double_t d2 = -1;
1593
1594   switch (matching) {
1595   case kGeometrical:
1596     GetGeometricalMatchingLevel(jet1,jet2,d1);
1597     d2 = d1;
1598     break;
1599   case kMCLabel: // jet1 = detector level and jet2 = particle level!
1600     GetMCLabelMatchingLevel(jet1,jet2,d1,d2);
1601     break;
1602   case kSameCollections:
1603     GetSameCollectionsMatchingLevel(jet1,jet2,d1,d2);
1604     break;
1605   default:
1606     ;
1607   }
1608
1609   if (d1 >= 0) {
1610
1611     if (d1 < jet1->ClosestJetDistance()) {
1612       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
1613       jet1->SetClosestJet(jet2, d1);
1614     }
1615     else if (d1 < jet1->SecondClosestJetDistance()) {
1616       jet1->SetSecondClosestJet(jet2, d1);
1617     }
1618   }
1619   
1620   if (d2 >= 0) {
1621     
1622     if (d2 < jet2->ClosestJetDistance()) {
1623       jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
1624       jet2->SetClosestJet(jet1, d2);
1625     }
1626     else if (d2 < jet2->SecondClosestJetDistance()) {
1627       jet2->SetSecondClosestJet(jet1, d2);
1628     }
1629   }
1630 }
1631
1632 //________________________________________________________________________
1633 Bool_t AliJetResponseMaker::FillHistograms()
1634 {
1635   // Fill histograms.
1636
1637   AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1638   AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1639
1640   if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
1641
1642   AliEmcalJet* jet1 = 0;  
1643   AliEmcalJet* jet2 = 0;
1644   
1645   jets2->ResetCurrentID();
1646   while ((jet2 = jets2->GetNextJet())) {
1647     AliDebug(2,Form("Processing jet (2) %d", jets2->GetCurrentID()));
1648
1649     if (fDoJet2Histogram)
1650       FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), 
1651                    jet2->Pt() - jets2->GetRhoVal() * jet2->Area(), jet2->MCPt(), 2);
1652
1653     // Verify accept jet
1654     if (jets2->AcceptJet(jet2))
1655       FillJetHisto(jet2->Phi(), jet2->Eta(), jet2->Pt(), jet2->Area(), jet2->NEF(), jet2->MaxPartPt()/jet2->Pt(), 
1656                    jet2->Pt() - jets2->GetRhoVal() * jet2->Area(), jet2->MCPt(), 3);
1657
1658     jet1 = jet2->MatchedJet();
1659
1660     if (!jet1) continue;
1661     if (!jets1->AcceptJet(jet1)) continue;
1662     if (jet1->MCPt() < fMinJetMCPt) continue;
1663
1664     Double_t d=-1, ce1=-1, ce2=-1;
1665     if (jet2->GetMatchingType() == kGeometrical) {
1666       if (jets2->GetIsParticleLevel() && !jets1->GetIsParticleLevel()) // the other way around is not supported
1667         GetMCLabelMatchingLevel(jet1, jet2, ce1, ce2);
1668       else if (jets1->GetIsParticleLevel() == jets2->GetIsParticleLevel())
1669         GetSameCollectionsMatchingLevel(jet1, jet2, ce1, ce2);
1670       
1671       d = jet2->ClosestJetDistance();
1672     }
1673     else if (jet2->GetMatchingType() == kMCLabel || jet2->GetMatchingType() == kSameCollections) {
1674       GetGeometricalMatchingLevel(jet1, jet2, d);
1675
1676       ce1 = jet1->ClosestJetDistance();
1677       ce2 = jet2->ClosestJetDistance();
1678     }
1679
1680     Double_t corrpt1 = jet1->Pt() - jets1->GetRhoVal() * jet1->Area();
1681     Double_t corrpt2 = jet2->Pt() - jets2->GetRhoVal() * jet2->Area();
1682
1683     FillMatchingHistos(jet1->Pt(), jet2->Pt(), jet1->Eta(), jet2->Eta(), jet1->Phi(), jet2->Phi(), 
1684                        jet1->Area(), jet2->Area(), d, ce1, ce2, corrpt1, corrpt2, jet1->MCPt(), 
1685                        jet1->NEF(), jet2->NEF(), jet1->MaxPartPt()/jet1->Pt(), jet2->MaxPartPt()/jet2->Pt());
1686   }
1687
1688   jets1->ResetCurrentID();
1689   while ((jet1 = jets1->GetNextAcceptJet())) {
1690     if (jet1->MCPt() < fMinJetMCPt) continue;
1691     AliDebug(2,Form("Processing jet (1) %d", jets1->GetCurrentID()));
1692
1693     FillJetHisto(jet1->Phi(), jet1->Eta(), jet1->Pt(), jet1->Area(), jet1->NEF(), jet1->MaxPartPt()/jet1->Pt(), 
1694                  jet1->Pt() - jets1->GetRhoVal() * jet1->Area(), jet1->MCPt(), 1);
1695   }
1696
1697   return kTRUE;
1698 }