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