Charged jets (pPb): Improved trackcut analysis
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskDeltaPtJEmb.cxx
CommitLineData
3a02e8fb 1// $Id$
2//
3// Jet embedding deltaPt task.
4//
5// Author: M. Verweij
6
7#include <TClonesArray.h>
8#include <TH1F.h>
9#include <TH2F.h>
10#include <TH3F.h>
11#include <TList.h>
12#include <TLorentzVector.h>
13#include <TRandom3.h>
14
15#include "AliVCluster.h"
16#include "AliVParticle.h"
17#include "AliVTrack.h"
18#include "AliEmcalJet.h"
19#include "AliRhoParameter.h"
20#include "AliLog.h"
21#include "AliJetContainer.h"
22#include "AliParticleContainer.h"
23#include "AliClusterContainer.h"
24
25#include "AliAnalysisTaskDeltaPtJEmb.h"
26
27ClassImp(AliAnalysisTaskDeltaPtJEmb)
28
29//________________________________________________________________________
30AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb() :
31 AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPtJEmb", kTRUE),
32 fJetsCont(0),
33 fTracksCont(0),
34 fCaloClustersCont(0),
35 fMinFractionShared(0.5),
36 fHistEmbJetsPtArea(0),
37 fHistEmbJetsCorrPtArea(0),
38 fHistEmbPartPtvsJetPt(0),
39 fHistEmbPartPtvsJetCorrPt(0),
40 fHistJetPtvsJetCorrPt(0),
41 fHistDistLeadPart2JetAxis(0),
42 fHistEmbBkgArea(0),
43 fHistRhoVSEmbBkg(0),
44 fHistDeltaPtEmbArea(0),
45 fHistDeltaPtEmbvsEP(0),
46 fHistDeltaPtEmbPtProbe(0),
47 fHistEmbJetsPhiEta(0),
48 fHistLeadPartPhiEta(0)
49{
50 // Default constructor.
51
52 fHistEmbJetsPtArea = new TH3*[fNcentBins];
53 fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
54 fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
55 fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
56 fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
57 fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
58 fHistEmbBkgArea = new TH2*[fNcentBins];
59 fHistRhoVSEmbBkg = new TH2*[fNcentBins];
60 fHistDeltaPtEmbArea = new TH2*[fNcentBins];
61 fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
62 fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];
63
64 for (Int_t i = 0; i < fNcentBins; i++) {
65 fHistEmbJetsPtArea[i] = 0;
66 fHistEmbJetsCorrPtArea[i] = 0;
67 fHistEmbPartPtvsJetPt[i] = 0;
68 fHistEmbPartPtvsJetCorrPt[i] = 0;
69 fHistJetPtvsJetCorrPt[i] = 0;
70 fHistDistLeadPart2JetAxis[i] = 0;
71 fHistEmbBkgArea[i] = 0;
72 fHistRhoVSEmbBkg[i] = 0;
73 fHistDeltaPtEmbArea[i] = 0;
74 fHistDeltaPtEmbvsEP[i] = 0;
75 fHistDeltaPtEmbPtProbe[i] = 0;
76 }
77
78 SetMakeGeneralHistograms(kTRUE);
79}
80
81//________________________________________________________________________
82AliAnalysisTaskDeltaPtJEmb::AliAnalysisTaskDeltaPtJEmb(const char *name) :
83 AliAnalysisTaskEmcalJet(name, kTRUE),
84 fJetsCont(0),
85 fTracksCont(0),
86 fCaloClustersCont(0),
87 fMinFractionShared(0.5),
88 fHistEmbJetsPtArea(0),
89 fHistEmbJetsCorrPtArea(0),
90 fHistEmbPartPtvsJetPt(0),
91 fHistEmbPartPtvsJetCorrPt(0),
92 fHistJetPtvsJetCorrPt(0),
93 fHistDistLeadPart2JetAxis(0),
94 fHistEmbBkgArea(0),
95 fHistRhoVSEmbBkg(0),
96 fHistDeltaPtEmbArea(0),
97 fHistDeltaPtEmbvsEP(0),
98 fHistDeltaPtEmbPtProbe(0),
99 fHistEmbJetsPhiEta(0),
100 fHistLeadPartPhiEta(0)
101{
102 // Standard constructor.
103
104 fHistEmbJetsPtArea = new TH3*[fNcentBins];
105 fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
106 fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
107 fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
108 fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
109 fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
110 fHistEmbBkgArea = new TH2*[fNcentBins];
111 fHistRhoVSEmbBkg = new TH2*[fNcentBins];
112 fHistDeltaPtEmbArea = new TH2*[fNcentBins];
113 fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
114 fHistDeltaPtEmbPtProbe = new TH2*[fNcentBins];
115
116 for (Int_t i = 0; i < fNcentBins; i++) {
117 fHistEmbJetsPtArea[i] = 0;
118 fHistEmbJetsCorrPtArea[i] = 0;
119 fHistEmbPartPtvsJetPt[i] = 0;
120 fHistEmbPartPtvsJetCorrPt[i] = 0;
121 fHistJetPtvsJetCorrPt[i] = 0;
122 fHistDistLeadPart2JetAxis[i] = 0;
123 fHistEmbBkgArea[i] = 0;
124 fHistRhoVSEmbBkg[i] = 0;
125 fHistDeltaPtEmbArea[i] = 0;
126 fHistDeltaPtEmbvsEP[i] = 0;
127 fHistDeltaPtEmbPtProbe[i] = 0;
128 }
129
130 SetMakeGeneralHistograms(kTRUE);
131}
132
133//________________________________________________________________________
134void AliAnalysisTaskDeltaPtJEmb::UserCreateOutputObjects()
135{
136 // Create user output.
137
138 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
139
140 fJetsCont = GetJetContainer("Jets");
141 fTracksCont = GetParticleContainer("Tracks");
142 fCaloClustersCont = GetClusterContainer("CaloClusters");
143
144 fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
145 fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
146 fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
147 fOutput->Add(fHistEmbJetsPhiEta);
148
149 fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
150 fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
151 fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
152 fOutput->Add(fHistLeadPartPhiEta);
153
154 TString histname;
155
156 const Int_t nbinsZ = 12;
157 Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
158
159 Float_t *binsPt = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
160 Float_t *binsCorrPt = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
161 Float_t *binsArea = GenerateFixedBinArray(50, 0, 2);
162
163 for (Int_t i = 0; i < fNcentBins; i++) {
164 histname = "fHistEmbJetsPtArea_";
165 histname += i;
166 fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
167 fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
168 fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
169 fOutput->Add(fHistEmbJetsPtArea[i]);
170
171 histname = "fHistEmbJetsCorrPtArea_";
172 histname += i;
173 fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
174 fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
175 fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
176 fOutput->Add(fHistEmbJetsCorrPtArea[i]);
177
178 histname = "fHistEmbPartPtvsJetPt_";
179 histname += i;
180 fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
181 fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
182 fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
183 fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
184 fOutput->Add(fHistEmbPartPtvsJetPt[i]);
185
186 histname = "fHistEmbPartPtvsJetCorrPt_";
187 histname += i;
188 fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
189 fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
190 fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
191 fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
192 fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
193 fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
194
195 histname = "fHistJetPtvsJetCorrPt_";
196 histname += i;
197 fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
198 fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
199 fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
200 fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
201 fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
202 fOutput->Add(fHistJetPtvsJetCorrPt[i]);
203
204 histname = "fHistDistLeadPart2JetAxis_";
205 histname += i;
206 fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
207 fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
208 fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
209 fOutput->Add(fHistDistLeadPart2JetAxis[i]);
210
211 histname = "fHistEmbBkgArea_";
212 histname += i;
213 fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
214 fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
215 fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
216 fOutput->Add(fHistEmbBkgArea[i]);
217
218 histname = "fHistRhoVSEmbBkg_";
219 histname += i;
220 fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
221 fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
222 fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
223 fOutput->Add(fHistRhoVSEmbBkg[i]);
224
225 histname = "fHistDeltaPtEmbArea_";
226 histname += i;
227 fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
228 50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
229 fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
230 fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
231 fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
232 fOutput->Add(fHistDeltaPtEmbArea[i]);
233
234 histname = "fHistDeltaPtEmbvsEP_";
235 histname += i;
236 fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(),
237 402, -TMath::Pi()*1.01, TMath::Pi()*3.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
238 fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
239 fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
240 fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
241 fOutput->Add(fHistDeltaPtEmbvsEP[i]);
242
243 histname = "fHistDeltaPtEmbPtProbe_";
244 histname += i;
245 fHistDeltaPtEmbPtProbe[i] = new TH2F(histname.Data(), histname.Data(),
246 fNbins, fMinBinPt, fMaxBinPt, fNbins * 2, -fMaxBinPt, fMaxBinPt);
247 fHistDeltaPtEmbPtProbe[i]->GetXaxis()->SetTitle("#it{p}_{T,probe} (GeV/#it{c})");
248 fHistDeltaPtEmbPtProbe[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
249 fHistDeltaPtEmbPtProbe[i]->GetZaxis()->SetTitle("counts");
250 fOutput->Add(fHistDeltaPtEmbPtProbe[i]);
251 }
252
253 delete[] binsPt;
254 delete[] binsCorrPt;
255 delete[] binsArea;
256
257 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
258}
259
260//________________________________________________________________________
261Bool_t AliAnalysisTaskDeltaPtJEmb::FillHistograms()
262{
263 // Fill histograms.
264
265 // ************
266 // Embedding
267 // Jets should already be matched, GetClosestJet gives matched partner. Cut on shared fraction to be applied
268 // _________________________________
269
270 if (fJetsCont) {
271 AliEmcalJet *jet1 = NULL;
272 AliEmcalJet *jet2 = NULL;
273 fJetsCont->ResetCurrentID();
274 while((jet1 = fJetsCont->GetNextAcceptJet())) {
275
276 jet2 = jet1->ClosestJet();
277 if(!jet2) continue;
278
279 Double_t fraction = fJetsCont->GetFractionSharedPt(jet1);
280 if(fMinFractionShared>0. && fraction<fMinFractionShared) continue;
281
282 TLorentzVector mom;
283 fJetsCont->GetLeadingHadronMomentum(mom,jet1);
284
285 Double_t distLeading2Jet = TMath::Sqrt((jet1->Eta() - mom.Eta()) * (jet1->Eta() - mom.Eta()) + (jet1->Phi() - mom.Phi()+TMath::Pi()) * (jet1->Phi() - mom.Phi()+TMath::Pi()));
286
287 fHistEmbPartPtvsJetPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt());
288 fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal);
289 fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi()+TMath::Pi());
290 fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
291
292 fHistEmbJetsPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt(), mom.Pt());
293 fHistEmbJetsCorrPtArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - fRhoVal * jet1->Area(), mom.Pt());
294 fHistEmbJetsPhiEta->Fill(jet1->Eta(), jet1->Phi());
295 fHistJetPtvsJetCorrPt[fCentBin]->Fill(jet1->Pt(), jet1->Pt() - fRhoVal * jet1->Area());
296
297 fHistEmbBkgArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet2->Pt());
298 fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * jet1->Area(), jet1->Pt() - jet2->Pt());
299 fHistDeltaPtEmbArea[fCentBin]->Fill(jet1->Area(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
300 fHistDeltaPtEmbvsEP[fCentBin]->Fill(jet1->Phi() - fEPV0, jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
301 fHistDeltaPtEmbPtProbe[fCentBin]->Fill(jet2->Pt(), jet1->Pt() - jet1->Area() * fRhoVal - jet2->Pt());
302 }
303 }
304
305 return kTRUE;
306}
307
308//________________________________________________________________________
309void AliAnalysisTaskDeltaPtJEmb::ExecOnce()
310{
311 // Initialize the analysis.
312
313 AliAnalysisTaskEmcalJet::ExecOnce();
314
315 if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
316 if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
317 if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
318}