Updates Salvatore Aiola
[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   if (fMatching == kNoMatching) 
484     return kTRUE;
485   else
486     return DoJetMatching();
487 }
488
489 //________________________________________________________________________
490 Bool_t AliJetResponseMaker::DoJetMatching()
491 {
492   DoJetLoop(kFALSE);
493
494   const Int_t nJets = fJets->GetEntriesFast();
495
496   for (Int_t i = 0; i < nJets; i++) {
497
498     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
499
500     if (!jet1) {
501       AliError(Form("Could not receive jet %d", i));
502       continue;
503     }  
504
505     if (!AcceptJet(jet1))
506       continue;
507
508     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
509       continue;
510
511     if (jet1->Pt() > fMaxBinPt)
512       continue;
513
514     if (jet1->ClosestJet() && jet1->ClosestJet()->ClosestJet() == jet1 && 
515         jet1->ClosestJetDistance() < fMatchingPar && jet1->ClosestJet()->ClosestJetDistance() < fMatchingPar) {    // Matched jet found
516       jet1->SetMatchedToClosest(fMatching);
517       jet1->ClosestJet()->SetMatchedToClosest(fMatching);
518     }
519   }
520
521   return kTRUE;
522 }
523
524 //________________________________________________________________________
525 void AliJetResponseMaker::DoJetLoop(Bool_t order)
526 {
527   // Do the jet loop.
528
529   TClonesArray *jets1 = 0;
530   TClonesArray *jets2 = 0;
531
532   if (order) {
533     jets1 = fJets2;
534     jets2 = fJets;
535   }
536   else {
537     jets1 = fJets;
538     jets2 = fJets2;
539   }
540
541   Int_t nJets1 = jets1->GetEntriesFast();
542   Int_t nJets2 = jets2->GetEntriesFast();
543
544   for (Int_t i = 0; i < nJets1; i++) {
545
546     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
547
548     if (!jet1) {
549       AliError(Form("Could not receive jet %d", i));
550       continue;
551     }  
552
553     if (!AcceptJet(jet1))
554       continue;
555
556     if (order) {
557      if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
558         continue;
559     }
560     else {
561       if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
562         continue;
563     }
564
565     for (Int_t j = 0; j < nJets2; j++) {
566       
567       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
568       
569       if (!jet2) {
570         AliError(Form("Could not receive jet %d", j));
571         continue;
572       }  
573       
574       if (!AcceptJet(jet2))
575         continue;
576
577       if (order) {
578         if (jet2->Eta() < fJetMinEta || jet2->Eta() > fJetMaxEta || jet2->Phi() < fJetMinPhi || jet2->Phi() > fJetMaxPhi)
579           continue;
580       }
581       else {
582         if (jet1->Eta() < fJet2MinEta || jet1->Eta() > fJet2MaxEta || jet1->Phi() < fJet2MinPhi || jet1->Phi() > fJet2MaxPhi)
583           continue;
584       }
585
586       GetMatchingLevel(jet1, jet2, fMatching);
587     }
588   }
589 }
590
591 //________________________________________________________________________
592 Double_t AliJetResponseMaker::GetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, MatchingType matching) 
593 {
594   Double_t d1 = -1;
595   Double_t d2 = -1;
596
597   switch (matching) {
598   case kGeometrical:
599     {
600       Double_t deta = jet2->Eta() - jet1->Eta();
601       Double_t dphi = jet2->Phi() - jet1->Phi();
602       d1 = TMath::Sqrt(deta * deta + dphi * dphi);
603       d2 = d1;
604     }
605     break;
606   case kMCLabel: // jet1 = detector level and jet2 = particle level!
607     { 
608       if (!fTracks2Map) {
609         fTracks2Map = new TH1I("tracksMap","tracksMap",1000,0,1);
610         for (Int_t i = 0; i < 1000; i++) {
611           fTracks2Map->SetBinContent(i,i);
612         }
613       }
614       d1 = jet1->Pt();
615       d2 = jet2->Pt();
616       Double_t totalPt1 = d1;
617       for (Int_t iTrack = 0; iTrack < jet1->GetNumberOfTracks(); iTrack++) {
618         AliVParticle *track = jet1->TrackAt(iTrack,fTracks);
619         if (!track) {
620           AliWarning(Form("Could not find track %d!", iTrack));
621           continue;
622         }
623         Int_t MClabel = track->GetLabel();
624         if (MClabel < 0) {// this is not a MC particle; remove it completely
625           AliDebug(3,Form("Track %d (pT = %f) is not a MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
626           totalPt1 -= track->Pt();
627           d1 -= track->Pt();
628           continue;
629         }
630         Int_t index = fTracks2Map->GetBinContent(MClabel);
631         if (index < 0) {
632           AliDebug(2,Form("Track %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iTrack,track->Pt(),MClabel));
633           continue;
634         }
635         for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
636           Int_t index2 = jet2->TrackAt(iTrack2);
637           if (index2 == index) { // found common particle
638             d1 -= track->Pt();
639             AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
640             AliDebug(3,Form("Track %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
641                             iTrack,track->Pt(),track->Eta(),track->Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
642             d2 -= MCpart->Pt();
643             break;
644           }
645         }
646       }
647       for (Int_t iClus = 0; iClus < jet1->GetNumberOfClusters(); iClus++) {
648         AliVCluster *clus = jet1->ClusterAt(iClus,fCaloClusters);
649         if (!clus) {
650           AliWarning(Form("Could not find cluster %d!", iClus));
651           continue;
652         }
653         TLorentzVector part;
654         clus->GetMomentum(part, fVertex);
655
656         if (fCaloCells) { // if the cell colection is available, look for cells with a matched MC particle
657           for (Int_t iCell = 0; iCell < clus->GetNCells(); iCell++) {
658             Int_t cellId = clus->GetCellAbsId(iCell);
659             Double_t cellFrac = clus->GetCellAmplitudeFraction(iCell);
660             Int_t MClabel = fCaloCells->GetCellMCLabel(cellId);
661
662             if (MClabel < 0) {// this is not a MC particle; remove it competely
663               AliDebug(3,Form("Cell %d (frac = %f) is not a MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
664               totalPt1 -= part.Pt() * cellFrac;
665               d1 -= part.Pt() * cellFrac;
666               continue;
667             }
668
669             Int_t index = fTracks2Map->GetBinContent(MClabel);
670             if (index < 0) {
671               AliDebug(3,Form("Cell %d (frac = %f) does not have an associated MC particle (MClabel = %d)!",iCell,cellFrac,MClabel));
672               continue;
673             }
674             for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
675               Int_t index2 = jet2->TrackAt(iTrack2);
676               if (index2 == index) { // found common particle
677                 d1 -= part.Pt() * cellFrac;
678                 AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
679                 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)!",
680                                 iCell,iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
681                 if (MCpart->Charge() != 0) // only if it is a neutral particle (charged particles are most likely already removed, to be fixed)
682                   d2 -= MCpart->Pt() * cellFrac;
683                 break;
684               }
685             }
686           }
687         }
688         else { //otherwise look for the first contributor to the cluster, and if matched to a MC label remove it
689           Int_t MClabel = clus->GetLabel();
690           
691           if (MClabel < 0) {// this is not a MC particle; remove it competely
692             AliDebug(3,Form("Cluster %d (pT = %f) is not a MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
693             totalPt1 -= part.Pt();
694             d1 -= part.Pt();
695             continue;
696           }
697           
698           Int_t index = fTracks2Map->GetBinContent(MClabel);
699           if (index < 0) {
700             AliDebug(3,Form("Cluster %d (pT = %f) does not have an associated MC particle (MClabel = %d)!",iClus,part.Pt(),MClabel));
701             continue;
702           }
703           for (Int_t iTrack2 = 0; iTrack2 < jet2->GetNumberOfTracks(); iTrack2++) {
704             Int_t index2 = jet2->TrackAt(iTrack2);
705             if (index2 == index) { // found common particle
706               d1 -= part.Pt();
707               AliVParticle *MCpart = static_cast<AliVParticle*>(fTracks2->At(index2));
708               AliDebug(3,Form("Cluster %d (pT = %f, eta = %f, phi = %f) is associated with the MC particle %d (pT = %f, eta = %f, phi = %f)!",
709                               iClus,part.Pt(),part.Eta(),part.Phi(),MClabel,MCpart->Pt(),MCpart->Eta(),MCpart->Phi()));
710               if (MCpart->Charge() != 0) // only if it is a neutral particle (charged particles are most likely already removed, to be fixed)
711                 d2 -= MCpart->Pt();
712               break;
713             }
714           }
715         }
716       }
717       if (d1 <= 0 || totalPt1 < 1)
718         d1 = 0;
719       else
720         d1 /= totalPt1;
721       if (jet2->Pt() > 0)
722         d2 /= jet2->Pt();
723       else
724         d2 = 0;
725     }
726     break;
727   default:
728     ;
729   }
730
731   if (d1 > 0) {
732
733     if (d1 < jet1->ClosestJetDistance()) {
734       jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
735       jet1->SetClosestJet(jet2, d1);
736     }
737     else if (d1 < jet1->SecondClosestJetDistance()) {
738       jet1->SetSecondClosestJet(jet2, d1);
739     }
740   }
741   
742   if (d2 > 0) {
743     
744     if (d2 < jet2->ClosestJetDistance()) {
745       jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
746       jet2->SetClosestJet(jet1, d2);
747     }
748     else if (d2 < jet2->SecondClosestJetDistance()) {
749       jet2->SetSecondClosestJet(jet1, d2);
750     }
751   }
752   
753   return d1;
754 }
755
756 //________________________________________________________________________
757 Bool_t AliJetResponseMaker::FillHistograms()
758 {
759   // Fill histograms.
760
761   fHistEvents->SetBinContent(fPtHardBin + 1, fHistEvents->GetBinContent(fPtHardBin + 1) + 1);
762   fHistNTrials->SetBinContent(fPtHardBin + 1, fHistNTrials->GetBinContent(fPtHardBin + 1) + fNTrials);
763
764   const Int_t nJets2 = fJets2->GetEntriesFast();
765
766   for (Int_t i = 0; i < nJets2; i++) {
767
768     AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(fJets2->At(i));
769
770     if (!jet2) {
771       AliError(Form("Could not receive jet2 %d", i));
772       continue;
773     }
774
775     if (jet2->Pt() > fMaxBinPt)
776       continue;
777
778     if (!AcceptJet(jet2))
779       continue;
780
781     if (AcceptBiasJet(jet2) &&
782         (jet2->Eta() > fJetMinEta && jet2->Eta() < fJetMaxEta && jet2->Phi() > fJetMinPhi && jet2->Phi() < fJetMaxPhi)) {
783       
784       fHistJets2PtAreaAcceptance->Fill(jet2->Area(), jet2->Pt());
785       fHistJets2PhiEtaAcceptance->Fill(jet2->Eta(), jet2->Phi());
786       
787       if (!fRho2Name.IsNull())
788         fHistJets2CorrPtAreaAcceptance->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
789     }
790
791     if (!AcceptBiasJet2(jet2))
792       continue;
793
794     if (jet2->Eta() < fJet2MinEta || jet2->Eta() > fJet2MaxEta || jet2->Phi() < fJet2MinPhi || jet2->Phi() > fJet2MaxPhi)
795       continue;
796
797     fHistJets2PtArea->Fill(jet2->Area(), jet2->Pt());
798     fHistJets2PhiEta->Fill(jet2->Eta(), jet2->Phi());
799
800     if (!fRho2Name.IsNull())
801       fHistJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRho2Val * jet2->Area());
802
803     if (jet2->MatchedJet()) {
804
805       if (!AcceptBiasJet(jet2->MatchedJet()) || 
806           jet2->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet2->MatchedJet()->MaxClusterPt() > fMaxClusterPt ||
807           jet2->MatchedJet()->Pt() > fMaxBinPt) {
808         fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
809       }
810       else {
811         if (jet2->GetMatchingType() == kGeometrical)
812           fHistDistancevsCommonEnergy->Fill(jet2->ClosestJetDistance(), GetMatchingLevel(jet2->MatchedJet(), jet2, kMCLabel));
813         else if (jet2->GetMatchingType() == kMCLabel)
814           fHistDistancevsCommonEnergy->Fill(GetMatchingLevel(jet2->MatchedJet(), jet2, kGeometrical), jet2->ClosestJetDistance());
815         else
816           fHistDistancevsCommonEnergy->Fill(GetMatchingLevel(jet2->MatchedJet(), jet2, kGeometrical), GetMatchingLevel(jet2->MatchedJet(), jet2, kMCLabel));
817           
818         fHistMatchingLevelvsJet2Pt->Fill(jet2->ClosestJetDistance(), jet2->Pt());
819
820         Double_t deta = jet2->MatchedJet()->Eta() - jet2->Eta();
821         Double_t dphi = jet2->MatchedJet()->Phi() - jet2->Phi();
822         fHistDeltaEtaPhivsJet2Pt->Fill(deta, dphi, jet2->Pt());
823
824         Double_t dpt = jet2->MatchedJet()->Pt() - jet2->Pt();
825         fHistDeltaPtvsJet2Pt->Fill(jet2->Pt(), dpt);
826         fHistDeltaPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
827
828         fHistJet1PtvsJet2Pt->Fill(jet2->MatchedJet()->Pt(), jet2->Pt());
829         
830         if (!fRhoName.IsNull() || !fRho2Name.IsNull()) {
831           dpt -= fRhoVal * jet2->MatchedJet()->Area() - fRho2Val * jet2->Area();
832           fHistDeltaCorrPtvsJet2Pt->Fill(jet2->Pt(), dpt);
833           fHistDeltaCorrPtvsMatchingLevel->Fill(jet2->ClosestJetDistance(), dpt);
834           fHistJet1CorrPtvsJet2CorrPt->Fill(jet2->MatchedJet()->Pt() - fRhoVal * jet2->MatchedJet()->Area(), jet2->Pt() - fRho2Val * jet2->Area());
835         }
836       }
837     }
838     else {
839       fHistNonMatchedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
840       fHistMissedJets2PtArea->Fill(jet2->Area(), jet2->Pt());
841
842       if (!fRho2Name.IsNull())
843         fHistNonMatchedJets2CorrPtArea->Fill(jet2->Area(), jet2->Pt() - fRhoVal * jet2->Area());
844     }
845   }
846
847   const Int_t nJets1 = fJets->GetEntriesFast();
848
849   for (Int_t i = 0; i < nJets1; i++) {
850
851     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(fJets->At(i));
852
853     if (!jet1) {
854       AliError(Form("Could not receive jet %d", i));
855       continue;
856     }  
857     
858     if (!AcceptJet(jet1))
859       continue;
860
861     if (!AcceptBiasJet(jet1))
862       continue;
863
864     if (jet1->MaxTrackPt() > fMaxTrackPt || jet1->MaxClusterPt() > fMaxClusterPt)
865       continue;
866
867     if (jet1->Eta() < fJetMinEta || jet1->Eta() > fJetMaxEta || jet1->Phi() < fJetMinPhi || jet1->Phi() > fJetMaxPhi)
868       continue;
869
870     if (!jet1->MatchedJet()) {
871       fHistNonMatchedJets1PtArea->Fill(jet1->Area(), jet1->Pt());
872       if (!fRhoName.IsNull())
873         fHistNonMatchedJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
874     }
875
876     fHistJets1PtArea->Fill(jet1->Area(), jet1->Pt());
877     fHistJets1PhiEta->Fill(jet1->Eta(), jet1->Phi());
878
879     if (!fRhoName.IsNull())
880       fHistJets1CorrPtArea->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area());
881   }
882
883   return kTRUE;
884 }