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