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