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