embedding (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
15 #include "AliVCluster.h"
16 #include "AliVTrack.h"
17 #include "AliEmcalJet.h"
18 #include "AliGenPythiaEventHeader.h"
19 #include "AliLog.h"
20 #include "AliMCEvent.h"
21 #include "AliRhoParameter.h"
22
23 ClassImp(AliJetResponseMaker)
24
25 //________________________________________________________________________
26 AliJetResponseMaker::AliJetResponseMaker() : 
27   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
28   fTracks2Name(""),
29   fCalo2Name(""),
30   fJets2Name(""),
31   fRho2Name(""),
32   fPtBiasJet2Track(0),
33   fPtBiasJet2Clus(0),
34   fAreCollections1MC(kFALSE),  
35   fAreCollections2MC(kTRUE),
36   fMatching(kNoMatching),
37   fMatchingPar(0),
38   fJet2MinEta(-999),
39   fJet2MaxEta(-999),
40   fJet2MinPhi(-999),
41   fJet2MaxPhi(-999),
42   fSelectPtHardBin(-999),
43   fPythiaHeader(0),
44   fPtHardBin(0),
45   fNTrials(0),
46   fTracks2(0),
47   fCaloClusters2(0),
48   fJets2(0),
49   fRho2(0),
50   fRho2Val(0),
51   fTracks2Map(0),
52   fHistNTrials(0),
53   fHistEvents(0),
54   fHistJets1PhiEta(0),
55   fHistJets1PtArea(0),
56   fHistJets1CorrPtArea(0),
57   fHistJets2PhiEta(0),
58   fHistJets2PtArea(0),
59   fHistJets2CorrPtArea(0),
60   fHistJets2PhiEtaAcceptance(0),
61   fHistJets2PtAreaAcceptance(0),
62   fHistJets2CorrPtAreaAcceptance(0),
63   fHistMatchingLevelvsJet2Pt(0),
64   fHistDistancevsCommonEnergy(0),
65   fHistDeltaEtaPhivsJet2Pt(0),
66   fHistDeltaPtvsJet2Pt(0),
67   fHistDeltaPtvsMatchingLevel(0),
68   fHistDeltaCorrPtvsJet2Pt(0),
69   fHistDeltaCorrPtvsMatchingLevel(0),
70   fHistNonMatchedJets1PtArea(0),
71   fHistNonMatchedJets2PtArea(0),
72   fHistNonMatchedJets1CorrPtArea(0),
73   fHistNonMatchedJets2CorrPtArea(0),
74   fHistJet1PtvsJet2Pt(0),
75   fHistJet1CorrPtvsJet2CorrPt(0),
76   fHistMissedJets2PtArea(0)
77 {
78   // Default constructor.
79
80   SetMakeGeneralHistograms(kTRUE);
81 }
82
83 //________________________________________________________________________
84 AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
85   AliAnalysisTaskEmcalJet(name, kTRUE),
86   fTracks2Name("MCParticles"),
87   fCalo2Name(""),
88   fJets2Name("MCJets"),
89   fRho2Name(""),
90   fPtBiasJet2Track(0),
91   fPtBiasJet2Clus(0),
92   fAreCollections1MC(kFALSE),  
93   fAreCollections2MC(kTRUE),
94   fMatching(kNoMatching),
95   fMatchingPar(0.25),
96   fJet2MinEta(-999),
97   fJet2MaxEta(-999),
98   fJet2MinPhi(-999),
99   fJet2MaxPhi(-999),
100   fSelectPtHardBin(-999),
101   fPythiaHeader(0),
102   fPtHardBin(0),
103   fNTrials(0),
104   fTracks2(0),
105   fCaloClusters2(0),
106   fJets2(0),
107   fRho2(0),
108   fRho2Val(0),
109   fTracks2Map(0),
110   fHistNTrials(0),
111   fHistEvents(0),
112   fHistJets1PhiEta(0),
113   fHistJets1PtArea(0),
114   fHistJets1CorrPtArea(0),
115   fHistJets2PhiEta(0),
116   fHistJets2PtArea(0),
117   fHistJets2CorrPtArea(0),
118   fHistJets2PhiEtaAcceptance(0),
119   fHistJets2PtAreaAcceptance(0),
120   fHistJets2CorrPtAreaAcceptance(0),
121   fHistMatchingLevelvsJet2Pt(0),
122   fHistDistancevsCommonEnergy(0),
123   fHistDeltaEtaPhivsJet2Pt(0),
124   fHistDeltaPtvsJet2Pt(0),
125   fHistDeltaPtvsMatchingLevel(0),
126   fHistDeltaCorrPtvsJet2Pt(0),
127   fHistDeltaCorrPtvsMatchingLevel(0),
128   fHistNonMatchedJets1PtArea(0),
129   fHistNonMatchedJets2PtArea(0),
130   fHistNonMatchedJets1CorrPtArea(0),
131   fHistNonMatchedJets2CorrPtArea(0),
132   fHistJet1PtvsJet2Pt(0),
133   fHistJet1CorrPtvsJet2CorrPt(0),
134   fHistMissedJets2PtArea(0)
135 {
136   // Standard constructor.
137
138   SetMakeGeneralHistograms(kTRUE);
139 }
140
141 //________________________________________________________________________
142 AliJetResponseMaker::~AliJetResponseMaker()
143 {
144   // Destructor
145 }
146
147 //________________________________________________________________________
148 void AliJetResponseMaker::UserCreateOutputObjects()
149 {
150   // Create user objects.
151
152   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
153
154   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
155   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
156
157   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
158   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
159   fHistNTrials->GetYaxis()->SetTitle("trials");
160   fOutput->Add(fHistNTrials);
161
162   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
163   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
164   fHistEvents->GetYaxis()->SetTitle("total events");
165   fOutput->Add(fHistEvents);
166
167   for (Int_t i = 1; i < 12; i++) {
168     fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
169     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
170   }
171
172   fHistJets1PhiEta = new TH2F("fHistJets1PhiEta", "fHistJets1PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
173   fHistJets1PhiEta->GetXaxis()->SetTitle("#eta");
174   fHistJets1PhiEta->GetYaxis()->SetTitle("#phi");
175   fOutput->Add(fHistJets1PhiEta);
176   
177   fHistJets1PtArea = new TH2F("fHistJets1PtArea", "fHistJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
178   fHistJets1PtArea->GetXaxis()->SetTitle("area");
179   fHistJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
180   fHistJets1PtArea->GetZaxis()->SetTitle("counts");
181   fOutput->Add(fHistJets1PtArea);
182
183   if (!fRhoName.IsNull()) {
184     fHistJets1CorrPtArea = new TH2F("fHistJets1CorrPtArea", "fHistJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
185     fHistJets1CorrPtArea->GetXaxis()->SetTitle("area");
186     fHistJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1}^{corr} (GeV/c)");
187     fHistJets1CorrPtArea->GetZaxis()->SetTitle("counts");
188     fOutput->Add(fHistJets1CorrPtArea);
189   }
190
191   fHistJets2PhiEta = new TH2F("fHistJets2PhiEta", "fHistJets2PhiEta", 40, -1, 1, 40, 0, TMath::Pi()*2);
192   fHistJets2PhiEta->GetXaxis()->SetTitle("#eta");
193   fHistJets2PhiEta->GetYaxis()->SetTitle("#phi");
194   fOutput->Add(fHistJets2PhiEta);
195
196   fHistJets2PhiEtaAcceptance = new TH2F("fHistJets2PhiEtaAcceptance", "fHistJets2PhiEtaAcceptance", 40, -1, 1, 40, 0, TMath::Pi()*2);
197   fHistJets2PhiEtaAcceptance->GetXaxis()->SetTitle("#eta");
198   fHistJets2PhiEtaAcceptance->GetYaxis()->SetTitle("#phi");
199   fOutput->Add(fHistJets2PhiEtaAcceptance);
200   
201   fHistJets2PtArea = new TH2F("fHistJets2PtArea", "fHistJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
202   fHistJets2PtArea->GetXaxis()->SetTitle("area");
203   fHistJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
204   fHistJets2PtArea->GetZaxis()->SetTitle("counts");
205   fOutput->Add(fHistJets2PtArea);
206
207   fHistJets2PtAreaAcceptance = new TH2F("fHistJets2PtAreaAcceptance", "fHistJets2PtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
208   fHistJets2PtAreaAcceptance->GetXaxis()->SetTitle("area");
209   fHistJets2PtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
210   fHistJets2PtAreaAcceptance->GetZaxis()->SetTitle("counts");
211   fOutput->Add(fHistJets2PtAreaAcceptance);
212
213   if (!fRho2Name.IsNull()) {
214     fHistJets2CorrPtArea = new TH2F("fHistJets2CorrPtArea", "fHistJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
215     fHistJets2CorrPtArea->GetXaxis()->SetTitle("area");
216     fHistJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
217     fHistJets2CorrPtArea->GetZaxis()->SetTitle("counts");
218     fOutput->Add(fHistJets2CorrPtArea);
219
220     fHistJets2CorrPtAreaAcceptance = new TH2F("fHistJets2CorrPtAreaAcceptance", "fHistJets2CorrPtAreaAcceptance", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
221     fHistJets2CorrPtAreaAcceptance->GetXaxis()->SetTitle("area");
222     fHistJets2CorrPtAreaAcceptance->GetYaxis()->SetTitle("p_{T,2}^{corr} (GeV/c)");
223     fHistJets2CorrPtAreaAcceptance->GetZaxis()->SetTitle("counts");
224     fOutput->Add(fHistJets2CorrPtAreaAcceptance);
225   }
226
227   fHistMatchingLevelvsJet2Pt = new TH2F("fHistMatchingLevelvsJet2Pt", "fHistMatchingLevelvsJet2Pt", fNbins/2, 0, 1.2, fNbins/2, fMinBinPt, fMaxBinPt);
228   fHistMatchingLevelvsJet2Pt->GetXaxis()->SetTitle("Matching level");
229   fHistMatchingLevelvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");  
230   fHistMatchingLevelvsJet2Pt->GetZaxis()->SetTitle("counts");
231   fOutput->Add(fHistMatchingLevelvsJet2Pt);
232
233   fHistDistancevsCommonEnergy = new TH2F("fHistDistancevsCommonEnergy", "fHistDistancevsCommonEnergy", fNbins/2, 0, 1.2, fNbins/2, 0, 1.2);
234   fHistDistancevsCommonEnergy->GetXaxis()->SetTitle("Distance");
235   fHistDistancevsCommonEnergy->GetYaxis()->SetTitle("Common energy (%)");  
236   fHistDistancevsCommonEnergy->GetZaxis()->SetTitle("counts");
237   fOutput->Add(fHistDistancevsCommonEnergy);
238
239   fHistDeltaEtaPhivsJet2Pt = new TH3F("fHistDeltaEtaPhivsJet2Pt", "fHistDeltaEtaPhivsJet2Pt", 40, -1, 1, 128, -1.6, 4.8, fNbins/2, fMinBinPt, fMaxBinPt);
240   fHistDeltaEtaPhivsJet2Pt->GetXaxis()->SetTitle("#Delta#eta");
241   fHistDeltaEtaPhivsJet2Pt->GetYaxis()->SetTitle("#Delta#phi");
242   fHistDeltaEtaPhivsJet2Pt->GetZaxis()->SetTitle("p_{T,2}");
243   fOutput->Add(fHistDeltaEtaPhivsJet2Pt);
244
245   fHistDeltaPtvsJet2Pt = new TH2F("fHistDeltaPtvsJet2Pt", "fHistDeltaPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
246   fHistDeltaPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
247   fHistDeltaPtvsJet2Pt->GetYaxis()->SetTitle("#deltap_{T} (GeV/c)");
248   fHistDeltaPtvsJet2Pt->GetZaxis()->SetTitle("counts");
249   fOutput->Add(fHistDeltaPtvsJet2Pt);
250
251   fHistDeltaPtvsMatchingLevel = new TH2F("fHistDeltaPtvsMatchingLevel", "fHistDeltaPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
252   fHistDeltaPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");  
253   fHistDeltaPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T} (GeV/c)");
254   fHistDeltaPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
255   fOutput->Add(fHistDeltaPtvsMatchingLevel);
256
257   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {  
258     fHistDeltaCorrPtvsJet2Pt = new TH2F("fHistDeltaCorrPtvsJet2Pt", "fHistDeltaCorrPtvsJet2Pt", fNbins/2, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
259     fHistDeltaCorrPtvsJet2Pt->GetXaxis()->SetTitle("p_{T,2}");  
260     fHistDeltaCorrPtvsJet2Pt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
261     fHistDeltaCorrPtvsJet2Pt->GetZaxis()->SetTitle("counts");
262     fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
263
264     fHistDeltaCorrPtvsMatchingLevel = new TH2F("fHistDeltaCorrPtvsMatchingLevel", "fHistDeltaCorrPtvsMatchingLevel", fNbins/2, 0, 1.2, 2*fNbins, -fMaxBinPt, fMaxBinPt);
265     fHistDeltaCorrPtvsMatchingLevel->GetXaxis()->SetTitle("Matching level");  
266     fHistDeltaCorrPtvsMatchingLevel->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)");
267     fHistDeltaCorrPtvsMatchingLevel->GetZaxis()->SetTitle("counts");
268     fOutput->Add(fHistDeltaCorrPtvsJet2Pt);
269   }
270
271   fHistNonMatchedJets1PtArea = new TH2F("fHistNonMatchedJets1PtArea", "fHistNonMatchedJets1PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
272   fHistNonMatchedJets1PtArea->GetXaxis()->SetTitle("area");
273   fHistNonMatchedJets1PtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
274   fHistNonMatchedJets1PtArea->GetZaxis()->SetTitle("counts");
275   fOutput->Add(fHistNonMatchedJets1PtArea);
276
277   fHistNonMatchedJets2PtArea = new TH2F("fHistNonMatchedJets2PtArea", "fHistNonMatchedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
278   fHistNonMatchedJets2PtArea->GetXaxis()->SetTitle("area");
279   fHistNonMatchedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
280   fHistNonMatchedJets2PtArea->GetZaxis()->SetTitle("counts");
281   fOutput->Add(fHistNonMatchedJets2PtArea);
282
283   if (!fRhoName.IsNull()) {  
284     fHistNonMatchedJets1CorrPtArea = new TH2F("fHistNonMatchedJets1CorrPtArea", "fHistNonMatchedJets1CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
285     fHistNonMatchedJets1CorrPtArea->GetXaxis()->SetTitle("area");
286     fHistNonMatchedJets1CorrPtArea->GetYaxis()->SetTitle("p_{T,1} (GeV/c)");
287     fHistNonMatchedJets1CorrPtArea->GetZaxis()->SetTitle("counts");
288     fOutput->Add(fHistNonMatchedJets1CorrPtArea);
289   }
290
291   if (!fRho2Name.IsNull()) {  
292     fHistNonMatchedJets2CorrPtArea = new TH2F("fHistNonMatchedJets2CorrPtArea", "fHistNonMatchedJets2CorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, 2*fNbins, -fMaxBinPt, fMaxBinPt);
293     fHistNonMatchedJets2CorrPtArea->GetXaxis()->SetTitle("area");
294     fHistNonMatchedJets2CorrPtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
295     fHistNonMatchedJets2CorrPtArea->GetZaxis()->SetTitle("counts");
296     fOutput->Add(fHistNonMatchedJets2CorrPtArea);
297   }
298
299   fHistJet1PtvsJet2Pt = new TH2F("fHistJet1PtvsJet2Pt", "fHistJet1PtvsJet2Pt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
300   fHistJet1PtvsJet2Pt->GetXaxis()->SetTitle("p_{T,1}");
301   fHistJet1PtvsJet2Pt->GetYaxis()->SetTitle("p_{T,2}");
302   fHistJet1PtvsJet2Pt->GetZaxis()->SetTitle("counts");
303   fOutput->Add(fHistJet1PtvsJet2Pt);
304
305   if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
306     if (fRhoName.IsNull()) 
307       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", fNbins, fMinBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
308     else if (fRho2Name.IsNull()) 
309       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
310     else
311       fHistJet1CorrPtvsJet2CorrPt = new TH2F("fHistJet1CorrPtvsJet2CorrPt", "fHistJet1CorrPtvsJet2CorrPt", 2*fNbins, -fMaxBinPt, fMaxBinPt, 2*fNbins, -fMaxBinPt, fMaxBinPt);
312     fHistJet1CorrPtvsJet2CorrPt->GetXaxis()->SetTitle("p_{T,1}^{corr}");
313     fHistJet1CorrPtvsJet2CorrPt->GetYaxis()->SetTitle("p_{T,2}^{corr}");
314     fHistJet1CorrPtvsJet2CorrPt->GetZaxis()->SetTitle("counts");
315     fOutput->Add(fHistJet1CorrPtvsJet2CorrPt);
316   }
317
318   fHistMissedJets2PtArea = new TH2F("fHistMissedJets2PtArea", "fHistMissedJets2PtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
319   fHistMissedJets2PtArea->GetXaxis()->SetTitle("area");  
320   fHistMissedJets2PtArea->GetYaxis()->SetTitle("p_{T,2} (GeV/c)");
321   fHistMissedJets2PtArea->GetZaxis()->SetTitle("counts");
322   fOutput->Add(fHistMissedJets2PtArea);
323
324   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
325 }
326
327 //________________________________________________________________________
328 Bool_t AliJetResponseMaker::AcceptJet(AliEmcalJet *jet) const
329 {   
330   // Return true if jet is accepted.
331
332   if (jet->Pt() <= fJetPtCut)
333     return kFALSE;
334   if (jet->Area() <= fJetAreaCut)
335     return kFALSE;
336
337   return kTRUE;
338 }
339
340 //________________________________________________________________________
341 Bool_t AliJetResponseMaker::AcceptBiasJet2(AliEmcalJet *jet) const
342
343   // Accept jet with a bias.
344
345   if (fLeadingHadronType == 0) {
346     if (jet->MaxTrackPt() < fPtBiasJet2Track) return kFALSE;
347   }
348   else if (fLeadingHadronType == 1) {
349     if (jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
350   }
351   else {
352     if (jet->MaxTrackPt() < fPtBiasJet2Track && jet->MaxClusterPt() < fPtBiasJet2Clus) return kFALSE;
353   }
354
355   return kTRUE;
356 }
357
358 //________________________________________________________________________
359 void AliJetResponseMaker::ExecOnce()
360 {
361   // Execute once.
362
363   if (!fJets2Name.IsNull() && !fJets2) {
364     fJets2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJets2Name));
365     if (!fJets2) {
366       AliError(Form("%s: Could not retrieve jets2 %s!", GetName(), fJets2Name.Data()));
367       return;
368     }
369     else if (!fJets2->GetClass()->GetBaseClass("AliEmcalJet")) {
370       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJets2Name.Data())); 
371       fJets2 = 0;
372       return;
373     }
374   }
375
376   if (!fTracks2Name.IsNull() && !fTracks2) {
377     fTracks2 = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracks2Name));
378     if (!fTracks2) {
379       AliError(Form("%s: Could not retrieve tracks2 %s!", GetName(), fTracks2Name.Data())); 
380       return;
381     }
382     else {
383       TClass *cl = fTracks2->GetClass();
384       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
385         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fTracks2Name.Data())); 
386         fTracks2 = 0;
387         return;
388       }
389     }
390
391     if (fAreCollections2MC) {
392       fTracks2Map = dynamic_cast<TH1*>(InputEvent()->FindListObject(fTracks2Name + "_Map"));
393       if (!fTracks2Map) {
394         AliError(Form("%s: Could not retrieve map for tracks2 %s!", GetName(), fTracks2Name.Data())); 
395         return;
396       }
397     }
398   }
399
400   if (!fCalo2Name.IsNull() && !fCaloClusters2) {
401     fCaloClusters2 =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCalo2Name));
402     if (!fCaloClusters2) {
403       AliError(Form("%s: Could not retrieve clusters %s!", GetName(), fCalo2Name.Data())); 
404       return;
405     } else {
406       TClass *cl = fCaloClusters2->GetClass();
407       if (!cl->GetBaseClass("AliVCluster") && !cl->GetBaseClass("AliEmcalParticle")) {
408         AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fCalo2Name.Data())); 
409         fCaloClusters2 = 0;
410         return;
411       }
412     }
413   }
414
415   if (!fRho2Name.IsNull() && !fRho2) {
416     fRho2 = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRho2Name));
417     if (!fRho2) {
418       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRho2Name.Data()));
419       fInitialized = kFALSE;
420       return;
421     }
422   }
423
424   if (fJet2MinEta == -999)
425     fJet2MinEta = fJetMinEta - fJetRadius;
426   if (fJet2MaxEta == -999)
427     fJet2MaxEta = fJetMaxEta + fJetRadius;
428   if (fJet2MinPhi == -999)
429     fJet2MinPhi = fJetMinPhi - fJetRadius;
430   if (fJet2MaxPhi == -999)
431     fJet2MaxPhi = fJetMaxPhi + fJetRadius;
432
433   AliAnalysisTaskEmcalJet::ExecOnce();
434 }
435
436 //________________________________________________________________________
437 Bool_t AliJetResponseMaker::IsEventSelected()
438 {
439   // Check if event is selected
440
441   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) 
442     return kFALSE;
443
444   return AliAnalysisTaskEmcalJet::IsEventSelected();
445 }
446
447 //________________________________________________________________________
448 Bool_t AliJetResponseMaker::RetrieveEventObjects()
449 {
450   // Retrieve event objects.
451
452   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
453     return kFALSE;
454
455   if (fRho2)
456     fRho2Val = fRho2->GetVal();
457   
458   fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
459
460   if (!fPythiaHeader)
461     return kFALSE;
462
463   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
464   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
465
466   Double_t pthard = fPythiaHeader->GetPtHard();
467
468   for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
469     if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
470       break;
471   }
472
473   fNTrials = fPythiaHeader->Trials();
474
475   return kTRUE;
476 }
477
478 //________________________________________________________________________
479 Bool_t AliJetResponseMaker::Run()
480 {
481   // Find the closest jets
482
483   switch (fMatching) {
484   case kGeometrical:
485     return GeometricalMatching();
486   case kMCLabel: 
487     return MCLabelMatching();
488   default:
489     return kTRUE;
490   }
491 }
492
493 //________________________________________________________________________
494 Bool_t AliJetResponseMaker::GeometricalMatching()
495 {
496   DoJetLoop(kFALSE);
497   DoJetLoop(kTRUE);
498
499   const Int_t nJets2 = fJets2->GetEntriesFast();
500
501   for (Int_t i = 0; i < nJets2; i++) {
502
503     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
504
505     if (!jet2) {
506       AliError(Form("Could not receive jet %d", i));
507       continue;
508     }  
509
510     if (!AcceptJet(jet2))
511       continue;
512
513     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
514       continue;
515
516     if (jet2->Pt() > fMaxBinPt)
517       continue;
518
519     if (jet2->ClosestJet() && jet2->ClosestJet()->ClosestJet() == jet2 && 
520         jet2->ClosestJetDistance() < fMatchingPar) {    // Matched jet found
521       jet2->SetMatchedToClosest(fMatching);
522       jet2->ClosestJet()->SetMatchedToClosest(fMatching);
523     }
524   }
525
526   return kTRUE;
527 }
528
529 //________________________________________________________________________
530 Bool_t AliJetResponseMaker::MCLabelMatching()
531 {
532   DoJetLoop(kFALSE);
533
534   const Int_t nJets1 = fJets->GetEntriesFast();
535
536   for (Int_t i = 0; i < nJets1; i++) {
537
538     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
539
540     if (!jet1) {
541       AliError(Form("Could not receive jet %d", i));
542       continue;
543     }  
544
545     if (!AcceptJet(jet1))
546       continue;
547
548     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
549       continue;
550
551     if (jet1->Pt() > fMaxBinPt)
552       continue;
553
554     if (jet1->ClosestJet() && jet1->ClosestJetDistance() < fMatchingPar) {    // Matched jet found
555       jet1->SetMatchedToClosest(fMatching);
556       jet1->ClosestJet()->SetClosestJet(jet1, jet1->ClosestJetDistance());
557       jet1->ClosestJet()->SetMatchedToClosest(fMatching);
558     }
559   }
560
561   return kTRUE;
562 }
563
564 //________________________________________________________________________
565 void AliJetResponseMaker::DoJetLoop(Bool_t order)
566 {
567   // Do the jet loop.
568
569   TClonesArray *jets1 = 0;
570   TClonesArray *jets2 = 0;
571
572   if (order) {
573     jets1 = fJets2;
574     jets2 = fJets;
575   }
576   else {
577     jets1 = fJets;
578     jets2 = fJets2;
579   }
580
581   Int_t nJets1 = jets1->GetEntriesFast();
582   Int_t nJets2 = jets2->GetEntriesFast();
583
584   for (Int_t i = 0; i < nJets1; i++) {
585
586     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
587
588     if (!jet1) {
589       AliError(Form("Could not receive jet %d", i));
590       continue;
591     }  
592
593     if (!AcceptJet(jet1))
594       continue;
595
596     if (order) {
597      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
598         continue;
599     }
600     else {
601       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
602         continue;
603     }
604
605     for (Int_t j = 0; j < nJets2; j++) {
606       
607       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
608       
609       if (!jet2) {
610         AliError(Form("Could not receive jet %d", j));
611         continue;
612       }  
613       
614       if (!AcceptJet(jet2))
615         continue;
616
617       if (order) {
618         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
619           continue;
620       }
621       else {
622         if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
623           continue;
624       }
625
626       Double_t d = GetMatchingLevel(jet1, jet2, fMatching);
627
628       if (d < 0)
629         continue;
630
631       if (d < jet1->ClosestJetDistance()) {
632         jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
633         jet1->SetClosestJet(jet2, d);
634       }
635       else if (d < jet1->SecondClosestJetDistance()) {
636         jet1->SetSecondClosestJet(jet2, d);
637       }
638     }
639   }
640 }
641
642 //________________________________________________________________________
643 Double_t AliJetResponseMaker::GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
644 {
645   Double_t r = -1;
646
647   switch (matching) {
648   case kGeometrical:
649     {
650       Double_t deta = jet2->Eta() - jet1->Eta();
651       Double_t dphi = jet2->Phi() - jet1->Phi();
652       r = TMath::Sqrt(deta * deta + dphi * dphi);
653     }
654     break;
655   case kMCLabel: // jet1 should be detector level and jet2 particle level!
656     { 
657       if (!fTracks2Map) {
658         fTracks2Map = new TH1I("tracksMap","tracksMap",1000,0,1);
659         for (Int_t i = 0; i < 1000; i++) {
660           fTracks2Map->SetBinContent(i,i);
661         }
662       }
663       r = jet1->Pt();
664       for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
665         AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
666         if (!track) {
667           AliWarning(Form("Could not find track %d!", iTrack));
668           continue;
669         }
670         Int_t MClabel = track->GetLabel();
671         Int_t index = fTracks2Map->GetBinContent(MClabel);
672         if (index < 0) {
673           AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
674           continue;
675         }
676         for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
677           Int_t index2 = jet2->TrackAt(iTrack2);
678           if (index2 == index) { // found common particle
679             r -= track->Pt();
680             break;
681           }
682         }
683       }
684       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
685         AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
686         if (!clus) {
687           AliWarning(Form("Could not find cluster %d!", iClus));
688           continue;
689         }
690         TLorentzVector part;
691         clus->GetMomentum(part, fVertex);
692         Int_t *MClabels = clus->GetLabels();
693         UInt_t nMClabels = clus->GetNLabels();
694         for (UInt_t j = 0; j < nMClabels; j++) {
695           Int_t index = fTracks2Map->GetBinContent(MClabels[j]);
696           if (index < 0) {
697             AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabels[j]));
698             continue;
699           }
700           for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
701             Int_t index2 = jet2->TrackAt(iTrack2);
702             if (index2 == index) { // found common particle
703               r -= part.Pt();
704               j = nMClabels; // stop looking for other MC particles
705               break;
706             }
707           }
708         }
709       }
710       if (r < 0)
711         r = 0;
712       r /= jet1->Pt();
713     }
714     break;
715   default:
716     ;
717   }
718   
719   return r;
720 }
721
722 //________________________________________________________________________
723 Bool_t AliJetResponseMaker::FillHistograms()
724 {
725   // Fill histograms.
726
727   fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
728   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
729
730   const Int_t nJets2 = fJets2->GetEntriesFast();
731
732   for (Int_t i = 0; i < nJets2; i++) {
733
734     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
735
736     if (!jet2) {
737       AliError(Form("Could not receive jet2 %d", i));
738       continue;
739     }
740
741     if (jet2->Pt() > fMaxBinPt)
742       continue;
743
744     if (!AcceptJet(jet2))
745       continue;
746
747     if (AcceptBiasJet(jet2) &&
748         (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
749       
750       fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
751       fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
752       
753       if (!fRho2Name.IsNull())
754         fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
755     }
756
757     if (!AcceptBiasJet2(jet2))
758       continue;
759
760     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
761       continue;
762
763     fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
764     fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
765
766     if (!fRho2Name.IsNull())
767       fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
768
769     if (jet2->MatchedJet()) {
770
771       if (!AcceptBiasJet(jet2->MatchedJet()) || 
772           jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
773           jet2->MatchedJet()->Pt() > fMaxBinPt) {
774         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
775       }
776       else {
777         if (jet2->GetMatchingType() == kGeometrical)
778           fHistDistancevsCommonEnergy->Fill(jet2->ClosestJetDistance(), GetMatchingLevel(jet2->MatchedJet(), jet2, kMCLabel));
779         else if (jet2->GetMatchingType() == kMCLabel)
780           fHistDistancevsCommonEnergy->Fill(GetMatchingLevel(jet2->MatchedJet(), jet2, kGeometrical), jet2->ClosestJetDistance());
781         else
782           fHistDistancevsCommonEnergy->Fill(GetMatchingLevel(jet2->MatchedJet(), jet2, kGeometrical), GetMatchingLevel(jet2->MatchedJet(), jet2, kMCLabel));
783           
784         fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
785
786         Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
787         Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
788         fHistDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
789
790         Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
791         fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
792         fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
793
794         fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
795         
796         if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
797           dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
798           fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
799           fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
800           fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
801         }
802       }
803     }
804     else {
805       fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
806       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
807
808       if (!fRho2Name.IsNull())
809         fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
810     }
811   }
812
813   const Int_t nJets1 = fJets->GetEntriesFast();
814
815   for (Int_t i = 0; i < nJets1; i++) {
816
817     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
818
819     if (!jet1) {
820       AliError(Form("Could not receive jet %d", i));
821       continue;
822     }  
823     
824     if (!AcceptJet(jet1))
825       continue;
826
827     if (!AcceptBiasJet(jet1))
828       continue;
829
830     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
831       continue;
832
833     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
834       continue;
835
836     if (!jet1->MatchedJet()) {
837       fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
838       if (!fRhoName.IsNull())
839         fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
840     }
841
842     fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
843     fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
844
845     if (!fRhoName.IsNull())
846       fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
847   }
848
849   return kTRUE;
850 }