]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeGR.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update GR analysis task
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetShapeGR.cxx
CommitLineData
e1239602 1//
2// Analysis task for angular jet shape G(R) arXiv:1201.2688
3//
4// Author: M.Verweij
5
6#include <TClonesArray.h>
7#include <TH1F.h>
8#include <TH2F.h>
9#include <TH3F.h>
10#include <THnSparse.h>
11#include <TF1.h>
12#include <TList.h>
13#include <TLorentzVector.h>
14#include <TProfile.h>
15#include <TChain.h>
16#include <TSystem.h>
17#include <TFile.h>
18#include <TKey.h>
19#include <TTree.h>
20
21#include "AliVCluster.h"
22#include "AliVTrack.h"
23#include "AliEmcalJet.h"
24#include "AliRhoParameter.h"
25#include "AliLog.h"
26#include "AliEmcalParticle.h"
27#include "AliMCEvent.h"
28#include "AliAODEvent.h"
29#include "AliGenPythiaEventHeader.h"
30#include "AliAODMCHeader.h"
31#include "AliAnalysisManager.h"
32#include "AliJetContainer.h"
33#include "AliParticleContainer.h"
34
35#include "AliAnalysisTaskJetShapeGR.h"
36
37ClassImp(AliAnalysisTaskJetShapeGR)
38
39//________________________________________________________________________
40AliAnalysisTaskJetShapeGR::AliAnalysisTaskJetShapeGR() :
41 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeGR", kTRUE),
42 fContainerBase(0),
43 fContainerSub(1),
44 fContainerTrue(2),
45 fMinFractionShared(0),
46 fSingleTrackEmb(kFALSE),
47 fCreateTree(kFALSE),
48 fTreeJetBkg(),
49 fJet1Vec(new TLorentzVector()),
50 fJet2Vec(new TLorentzVector()),
51 fJetSubVec(new TLorentzVector()),
52 fArea(0),
53 fAreaPhi(0),
54 fAreaEta(0),
55 fRho(0),
56 fRhoM(0),
57 fNConst(0),
58 fMatch(0),
15a3ff83 59 fDRStep(0.04),
60 fMaxR(2.),
e1239602 61 fh2PtTrueDeltaGR(0x0),
62 fh2PtTrueDeltaGRRel(0x0),
15a3ff83 63 fhnGRResponse(0x0),
64 fh1PtTrue(0x0),
65 fh3DeltaGRNumRPtTrue(0x0),
66 fh3DeltaGRDenRPtTrue(0x0),
67 fh2DeltaGRNumRPtTrue(0x0),
68 fh2DeltaGRDenRPtTrue(0x0),
69 fh1PtRaw(0x0),
70 fh3DeltaGRNumRPtRaw(0x0),
71 fh3DeltaGRDenRPtRaw(0x0),
72 fh2DeltaGRNumRPtRaw(0x0),
73 fh2DeltaGRDenRPtRaw(0x0),
74 fh1PtRawMatch(0x0),
75 fh3DeltaGRNumRPtRawMatch(0x0),
76 fh3DeltaGRDenRPtRawMatch(0x0),
77 fh2DeltaGRNumRPtRawMatch(0x0),
78 fh2DeltaGRDenRPtRawMatch(0x0),
79 fh1PtMatch(0x0),
80 fh3DeltaGRNumRPtMatch(0x0),
81 fh3DeltaGRDenRPtMatch(0x0),
82 fh2DeltaGRNumRPtMatch(0x0),
83 fh2DeltaGRDenRPtMatch(0x0),
84 fh2DeltaGRNumRPtTrueMatch(0x0),
85 fh2DeltaGRDenRPtTrueMatch(0x0)
e1239602 86{
87 // Default constructor.
88
e1239602 89 fh2PtTrueDeltaGR = new TH2F*[fNcentBins];
90 fh2PtTrueDeltaGRRel = new TH2F*[fNcentBins];
15a3ff83 91 fhnGRResponse = new THnSparse*[fNcentBins];
92 fh1PtTrue = new TH1F*[fNcentBins];
93 fh3DeltaGRNumRPtTrue = new TH3F*[fNcentBins];
94 fh3DeltaGRDenRPtTrue = new TH3F*[fNcentBins];
95 fh2DeltaGRNumRPtTrue = new TH2F*[fNcentBins];
96 fh2DeltaGRDenRPtTrue = new TH2F*[fNcentBins];
97 fh1PtRaw = new TH1F*[fNcentBins];
98 fh3DeltaGRNumRPtRaw = new TH3F*[fNcentBins];
99 fh3DeltaGRDenRPtRaw = new TH3F*[fNcentBins];
100 fh2DeltaGRNumRPtRaw = new TH2F*[fNcentBins];
101 fh2DeltaGRDenRPtRaw = new TH2F*[fNcentBins];
102 fh1PtRawMatch = new TH1F*[fNcentBins];
103 fh3DeltaGRNumRPtRawMatch = new TH3F*[fNcentBins];
104 fh3DeltaGRDenRPtRawMatch = new TH3F*[fNcentBins];
105 fh2DeltaGRNumRPtRawMatch = new TH2F*[fNcentBins];
106 fh2DeltaGRDenRPtRawMatch = new TH2F*[fNcentBins];
107 fh1PtMatch = new TH1F*[fNcentBins];
108 fh3DeltaGRNumRPtMatch = new TH3F*[fNcentBins];
109 fh3DeltaGRDenRPtMatch = new TH3F*[fNcentBins];
110 fh2DeltaGRNumRPtMatch = new TH2F*[fNcentBins];
111 fh2DeltaGRDenRPtMatch = new TH2F*[fNcentBins];
112 fh2DeltaGRNumRPtTrueMatch = new TH2F*[fNcentBins];
113 fh2DeltaGRDenRPtTrueMatch = new TH2F*[fNcentBins];
e1239602 114
115 for (Int_t i = 0; i < fNcentBins; i++) {
e1239602 116 fh2PtTrueDeltaGR[i] = 0;
117 fh2PtTrueDeltaGRRel[i] = 0;
15a3ff83 118 fhnGRResponse[i] = 0;
119 fh1PtTrue[i] = 0;
120 fh3DeltaGRNumRPtTrue[i] = 0;
121 fh3DeltaGRDenRPtTrue[i] = 0;
122 fh2DeltaGRNumRPtTrue[i] = 0;
123 fh2DeltaGRDenRPtTrue[i] = 0;
124 fh1PtRaw[i] = 0;
125 fh3DeltaGRNumRPtRaw[i] = 0;
126 fh3DeltaGRDenRPtRaw[i] = 0;
127 fh2DeltaGRNumRPtRaw[i] = 0;
128 fh2DeltaGRDenRPtRaw[i] = 0;
129 fh1PtRawMatch[i] = 0;
130 fh3DeltaGRNumRPtRawMatch[i] = 0;
131 fh3DeltaGRDenRPtRawMatch[i] = 0;
132 fh2DeltaGRNumRPtRawMatch[i] = 0;
133 fh2DeltaGRDenRPtRawMatch[i] = 0;
134 fh1PtMatch[i] = 0;
135 fh3DeltaGRNumRPtMatch[i] = 0;
136 fh3DeltaGRDenRPtMatch[i] = 0;
137 fh2DeltaGRNumRPtMatch[i] = 0;
138 fh2DeltaGRDenRPtMatch[i] = 0;
139 fh2DeltaGRNumRPtTrueMatch[i] = 0;
140 fh2DeltaGRDenRPtTrueMatch[i] = 0;
e1239602 141 }
142
143 SetMakeGeneralHistograms(kTRUE);
144 DefineOutput(2, TTree::Class());
145}
146
147//________________________________________________________________________
148AliAnalysisTaskJetShapeGR::AliAnalysisTaskJetShapeGR(const char *name) :
149 AliAnalysisTaskEmcalJet(name, kTRUE),
150 fContainerBase(0),
151 fContainerSub(1),
152 fContainerTrue(2),
153 fMinFractionShared(0),
154 fSingleTrackEmb(kFALSE),
155 fCreateTree(kFALSE),
156 fTreeJetBkg(0),
157 fJet1Vec(new TLorentzVector()),
158 fJet2Vec(new TLorentzVector()),
159 fJetSubVec(new TLorentzVector()),
160 fArea(0),
161 fAreaPhi(0),
162 fAreaEta(0),
163 fRho(0),
164 fRhoM(0),
165 fNConst(0),
166 fMatch(0),
15a3ff83 167 fDRStep(0.04),
168 fMaxR(2.),
e1239602 169 fh2PtTrueDeltaGR(0x0),
170 fh2PtTrueDeltaGRRel(0x0),
15a3ff83 171 fhnGRResponse(0x0),
172 fh1PtTrue(0x0),
173 fh3DeltaGRNumRPtTrue(0x0),
174 fh3DeltaGRDenRPtTrue(0x0),
175 fh2DeltaGRNumRPtTrue(0x0),
176 fh2DeltaGRDenRPtTrue(0x0),
177 fh1PtRaw(0x0),
178 fh3DeltaGRNumRPtRaw(0x0),
179 fh3DeltaGRDenRPtRaw(0x0),
180 fh2DeltaGRNumRPtRaw(0x0),
181 fh2DeltaGRDenRPtRaw(0x0),
182 fh1PtRawMatch(0x0),
183 fh3DeltaGRNumRPtRawMatch(0x0),
184 fh3DeltaGRDenRPtRawMatch(0x0),
185 fh2DeltaGRNumRPtRawMatch(0x0),
186 fh2DeltaGRDenRPtRawMatch(0x0),
187 fh1PtMatch(0x0),
188 fh3DeltaGRNumRPtMatch(0x0),
189 fh3DeltaGRDenRPtMatch(0x0),
190 fh2DeltaGRNumRPtMatch(0x0),
191 fh2DeltaGRDenRPtMatch(0x0),
192 fh2DeltaGRNumRPtTrueMatch(0x0),
193 fh2DeltaGRDenRPtTrueMatch(0x0)
e1239602 194{
195 // Standard constructor.
196
e1239602 197 fh2PtTrueDeltaGR = new TH2F*[fNcentBins];
198 fh2PtTrueDeltaGRRel = new TH2F*[fNcentBins];
15a3ff83 199 fhnGRResponse = new THnSparse*[fNcentBins];
200 fh1PtTrue = new TH1F*[fNcentBins];
201 fh3DeltaGRNumRPtTrue = new TH3F*[fNcentBins];
202 fh3DeltaGRDenRPtTrue = new TH3F*[fNcentBins];
203 fh2DeltaGRNumRPtTrue = new TH2F*[fNcentBins];
204 fh2DeltaGRDenRPtTrue = new TH2F*[fNcentBins];
205 fh1PtRaw = new TH1F*[fNcentBins];
206 fh3DeltaGRNumRPtRaw = new TH3F*[fNcentBins];
207 fh3DeltaGRDenRPtRaw = new TH3F*[fNcentBins];
208 fh2DeltaGRNumRPtRaw = new TH2F*[fNcentBins];
209 fh2DeltaGRDenRPtRaw = new TH2F*[fNcentBins];
210 fh1PtRawMatch = new TH1F*[fNcentBins];
211 fh3DeltaGRNumRPtRawMatch = new TH3F*[fNcentBins];
212 fh3DeltaGRDenRPtRawMatch = new TH3F*[fNcentBins];
213 fh2DeltaGRNumRPtRawMatch = new TH2F*[fNcentBins];
214 fh2DeltaGRDenRPtRawMatch = new TH2F*[fNcentBins];
215 fh1PtMatch = new TH1F*[fNcentBins];
216 fh3DeltaGRNumRPtMatch = new TH3F*[fNcentBins];
217 fh3DeltaGRDenRPtMatch = new TH3F*[fNcentBins];
218 fh2DeltaGRNumRPtMatch = new TH2F*[fNcentBins];
219 fh2DeltaGRDenRPtMatch = new TH2F*[fNcentBins];
220 fh2DeltaGRNumRPtTrueMatch = new TH2F*[fNcentBins];
221 fh2DeltaGRDenRPtTrueMatch = new TH2F*[fNcentBins];
e1239602 222
223 for (Int_t i = 0; i < fNcentBins; i++) {
e1239602 224 fh2PtTrueDeltaGR[i] = 0;
225 fh2PtTrueDeltaGRRel[i] = 0;
15a3ff83 226 fhnGRResponse[i] = 0;
227 fh1PtTrue[i] = 0;
228 fh3DeltaGRNumRPtTrue[i] = 0;
229 fh3DeltaGRDenRPtTrue[i] = 0;
230 fh2DeltaGRNumRPtTrue[i] = 0;
231 fh2DeltaGRDenRPtTrue[i] = 0;
232 fh1PtRaw[i] = 0;
233 fh3DeltaGRNumRPtRaw[i] = 0;
234 fh3DeltaGRDenRPtRaw[i] = 0;
235 fh2DeltaGRNumRPtRaw[i] = 0;
236 fh2DeltaGRDenRPtRaw[i] = 0;
237 fh1PtRawMatch[i] = 0;
238 fh3DeltaGRNumRPtRawMatch[i] = 0;
239 fh3DeltaGRDenRPtRawMatch[i] = 0;
240 fh2DeltaGRNumRPtRawMatch[i] = 0;
241 fh2DeltaGRDenRPtRawMatch[i] = 0;
242 fh1PtMatch[i] = 0;
243 fh3DeltaGRNumRPtMatch[i] = 0;
244 fh3DeltaGRDenRPtMatch[i] = 0;
245 fh2DeltaGRNumRPtMatch[i] = 0;
246 fh2DeltaGRDenRPtMatch[i] = 0;
247 fh2DeltaGRNumRPtTrueMatch[i] = 0;
248 fh2DeltaGRDenRPtTrueMatch[i] = 0;
e1239602 249 }
250
251 SetMakeGeneralHistograms(kTRUE);
252 DefineOutput(2, TTree::Class());
253}
254
255//________________________________________________________________________
256AliAnalysisTaskJetShapeGR::~AliAnalysisTaskJetShapeGR()
257{
258 // Destructor.
259}
260
261//________________________________________________________________________
262void AliAnalysisTaskJetShapeGR::UserCreateOutputObjects()
263{
264 // Create user output.
265
266 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
267
268 Bool_t oldStatus = TH1::AddDirectoryStatus();
269 TH1::AddDirectory(kFALSE);
270
271 const Int_t nBinsPt = 200;
272 const Double_t minPt = -50.;
273 const Double_t maxPt = 150.;
274
275 const Int_t nBinsM = 150;
276 const Double_t minM = -50.;
277 const Double_t maxM = 100.;
278
15a3ff83 279 const Int_t nBinsR = 50;
280 const Double_t minR = 0.;
281 const Double_t maxR = 2.;
282
283 const Int_t nBinsDGR = 100;
284 const Double_t minDGR = 0.;
285 const Double_t maxDGR = 10.;
286
e1239602 287 //Binning for THnSparse
288 const Int_t nBinsSparse0 = 4;
289 const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt};
290 const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt};
291 const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt};
292
293 TString histName = "";
294 TString histTitle = "";
295 for (Int_t i = 0; i < fNcentBins; i++) {
e1239602 296 histName = Form("fh2PtTrueDeltaGR_%d",i);
297 histTitle = Form("fh2PtTrueDeltaGR_%d;#it{p}_{T,true};#it{G(R)}_{sub}-#it{G(R)}_{true}",i);
298 fh2PtTrueDeltaGR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,-50.,50.);
299 fOutput->Add(fh2PtTrueDeltaGR[i]);
300
301 histName = Form("fh2PtTrueDeltaGRRel_%d",i);
302 histTitle = Form("fh2PtTrueDeltaGRRel_%d;#it{p}_{T,true};(#it{G(R)}_{sub}-#it{G(R)}_{true})/#it{G(R)}_{true}",i);
303 fh2PtTrueDeltaGRRel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,200,-1.,1.);
304 fOutput->Add(fh2PtTrueDeltaGRRel[i]);
305
306 histName = Form("fhnGRResponse_%d",i);
307 histTitle = Form("fhnGRResponse_%d;#it{G(R)}_{sub};#it{G(R)}_{true};#it{p}_{T,sub};#it{p}_{T,true}",i);
308 fhnGRResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
309 fOutput->Add(fhnGRResponse[i]);
15a3ff83 310
311 //Histos for true jets
312 histName = Form("fh1PtTrue_%d",i);
313 histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
314 fh1PtTrue[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
315 fOutput->Add(fh1PtTrue[i]);
316
317 histName = Form("fh3DeltaGRNumRPtTrue_%d",i);
318 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
319 fh3DeltaGRNumRPtTrue[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
320 fOutput->Add(fh3DeltaGRNumRPtTrue[i]);
321
322 histName = Form("fh3DeltaGRDenRPtTrue_%d",i);
323 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
324 fh3DeltaGRDenRPtTrue[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
325 fOutput->Add(fh3DeltaGRDenRPtTrue[i]);
326
327 histName = Form("fh2DeltaGRNumRPtTrue_%d",i);
328 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
329 fh2DeltaGRNumRPtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
330 fOutput->Add(fh2DeltaGRNumRPtTrue[i]);
331
332 histName = Form("fh2DeltaGRDenRPtTrue_%d",i);
333 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
334 fh2DeltaGRDenRPtTrue[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
335 fOutput->Add(fh2DeltaGRDenRPtTrue[i]);
336
337 //Histos for raw AA jets
338 histName = Form("fh1PtRaw_%d",i);
339 histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
340 fh1PtRaw[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
341 fOutput->Add(fh1PtRaw[i]);
342
343 histName = Form("fh3DeltaGRNumRPtRaw_%d",i);
344 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
345 fh3DeltaGRNumRPtRaw[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
346 fOutput->Add(fh3DeltaGRNumRPtRaw[i]);
347
348 histName = Form("fh3DeltaGRDenRPtRaw_%d",i);
349 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
350 fh3DeltaGRDenRPtRaw[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
351 fOutput->Add(fh3DeltaGRDenRPtRaw[i]);
352
353 histName = Form("fh2DeltaGRNumRPtRaw_%d",i);
354 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
355 fh2DeltaGRNumRPtRaw[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
356 fOutput->Add(fh2DeltaGRNumRPtRaw[i]);
357
358 histName = Form("fh2DeltaGRDenRPtRaw_%d",i);
359 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
360 fh2DeltaGRDenRPtRaw[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
361 fOutput->Add(fh2DeltaGRDenRPtRaw[i]);
362
363 //Histos for raw AA jets matched to MC jet
364 histName = Form("fh1PtRawMatch_%d",i);
365 histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
366 fh1PtRawMatch[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
367 fOutput->Add(fh1PtRawMatch[i]);
368
369 histName = Form("fh3DeltaGRNumRPtRawMatch_%d",i);
370 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
371 fh3DeltaGRNumRPtRawMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
372 fOutput->Add(fh3DeltaGRNumRPtRawMatch[i]);
373
374 histName = Form("fh3DeltaGRDenRPtRawMatch_%d",i);
375 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
376 fh3DeltaGRDenRPtRawMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
377 fOutput->Add(fh3DeltaGRDenRPtRawMatch[i]);
378
379 histName = Form("fh2DeltaGRNumRPtRawMatch_%d",i);
380 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
381 fh2DeltaGRNumRPtRawMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
382 fOutput->Add(fh2DeltaGRNumRPtRawMatch[i]);
383
384 histName = Form("fh2DeltaGRDenRPtRawMatch_%d",i);
385 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
386 fh2DeltaGRDenRPtRawMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
387 fOutput->Add(fh2DeltaGRDenRPtRawMatch[i]);
388
389 //Histos for matched jets
390 histName = Form("fh1PtMatch_%d",i);
391 histTitle = Form("%s;#it{p}_{T};#it{N}",histName.Data());
392 fh1PtMatch[i] = new TH1F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
393 fOutput->Add(fh1PtMatch[i]);
394
395 histName = Form("fh3DeltaGRNumRPtMatch_%d",i);
396 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#delta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
397 fh3DeltaGRNumRPtMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
398 fOutput->Add(fh3DeltaGRNumRPtMatch[i]);
399
400 histName = Form("fh3DeltaGRDenRPtMatch_%d",i);
401 histTitle = Form("%s;#Sigma#it{p}_{Tk,i}#it{p}_{Tk,j}#Delta#it{R}_{ij}^{2}#Theta_{dR}(#it{r}-#it{R}_{ij});#it{r};#it{p}_{T,true}",histName.Data());
402 fh3DeltaGRDenRPtMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsDGR,minDGR,maxDGR,nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
403 fOutput->Add(fh3DeltaGRDenRPtMatch[i]);
404
405 histName = Form("fh2DeltaGRNumRPtMatch_%d",i);
406 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
407 fh2DeltaGRNumRPtMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
408 fOutput->Add(fh2DeltaGRNumRPtMatch[i]);
409
410 histName = Form("fh2DeltaGRDenRPtMatch_%d",i);
411 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
412 fh2DeltaGRDenRPtMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
413 fOutput->Add(fh2DeltaGRDenRPtMatch[i]);
414
415 histName = Form("fh2DeltaGRNumRPtTrueMatch_%d",i);
416 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
417 fh2DeltaGRNumRPtTrueMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
418 fOutput->Add(fh2DeltaGRNumRPtTrueMatch[i]);
419
420 histName = Form("fh2DeltaGRDenRPtTrueMatch_%d",i);
421 histTitle = Form("%s;#it{r};#it{p}_{T,true}",histName.Data());
422 fh2DeltaGRDenRPtTrueMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsR,minR,maxR,nBinsPt,minPt,maxPt);
423 fOutput->Add(fh2DeltaGRDenRPtTrueMatch[i]);
424
e1239602 425 }
426
427 // =========== Switch on Sumw2 for all histos ===========
428 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
429 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
430 if (h1){
431 h1->Sumw2();
432 continue;
433 }
434 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
435 if(hn)hn->Sumw2();
436 }
437
438 TH1::AddDirectory(oldStatus);
439
440 // Create a tree.
441 if(fCreateTree) {
442 fTreeJetBkg = new TTree("fTreeJetSubGR", "fTreeJetSubGR");
443 fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
444 fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
445 fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
446 fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
447 fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
448 fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
449 fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
450 fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
451 fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
452 }
453
454 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
455 if(fCreateTree) PostData(2, fTreeJetBkg);
456}
457
458//________________________________________________________________________
459Bool_t AliAnalysisTaskJetShapeGR::Run()
460{
461 // Run analysis code here, if needed. It will be executed before FillHistograms().
462
463 return kTRUE;
464}
465
466//________________________________________________________________________
467Bool_t AliAnalysisTaskJetShapeGR::FillHistograms()
468{
469 // Fill histograms.
15a3ff83 470 FillTrueJets();
e1239602 471
472 AliEmcalJet* jet1 = NULL;
473 AliEmcalJet *jet2 = NULL;
474 AliEmcalJet *jetS = NULL;
475 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
476 AliJetContainer *jetContS = GetJetContainer(fContainerSub);
477 AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
15a3ff83 478 if(jetCont && jetContS) {
e1239602 479 jetCont->ResetCurrentID();
15a3ff83 480 Double_t rmax = jetCont->GetJetRadius()+0.2;
481 Double_t wr = 0.04;
482 Int_t nr = TMath::CeilNint(rmax/wr);
e1239602 483 while((jet1 = jetCont->GetNextAcceptJet())) {
e1239602 484 Double_t fraction = 0.;
485 fMatch = 0;
486 fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
487 if(fSingleTrackEmb) {
488 AliVParticle *vp = GetEmbeddedConstituent(jet1);
489 if(vp) {
490 fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
491 fMatch = 1;
492 }
493 } else {
494 jet2 = jet1->ClosestJet();
15a3ff83 495 if(!jet2) continue;
496
e1239602 497 fraction = jetCont->GetFractionSharedPt(jet1);
498 fMatch = 1;
499 if(fMinFractionShared>0.) {
500 if(fraction>fMinFractionShared) {
501 fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
502 fMatch = 1;
503 } else
504 fMatch = 0;
505 }
506 }
507
15a3ff83 508 //Fill histograms for all AA jets
509 Double_t ptcorr = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
510 fh1PtRaw[fCentBin]->Fill(ptcorr);
511 if(fMatch==1) fh1PtRawMatch[fCentBin]->Fill(ptcorr);
512
513 TArrayF numRaw = jet1->GetGRNumerator();
514 TArrayF denRaw = jet1->GetGRDenominator();
515 if(numRaw.GetSize()>0) {
516 for(Int_t i = 0; i<nr; i++) {
517 Double_t r = i*wr + 0.5*wr;
518 fh3DeltaGRNumRPtRaw[fCentBin]->Fill(numRaw[i],r,ptcorr);
519 fh3DeltaGRDenRPtRaw[fCentBin]->Fill(denRaw[i],r,ptcorr);
520 fh2DeltaGRNumRPtRaw[fCentBin]->Fill(r,ptcorr,numRaw[i]);
521 fh2DeltaGRDenRPtRaw[fCentBin]->Fill(r,ptcorr,denRaw[i]);
522 if(fMatch==1) {
523 fh3DeltaGRNumRPtRawMatch[fCentBin]->Fill(numRaw[i],r,ptcorr);
524 fh3DeltaGRDenRPtRawMatch[fCentBin]->Fill(denRaw[i],r,ptcorr);
525 fh2DeltaGRNumRPtRawMatch[fCentBin]->Fill(r,ptcorr,numRaw[i]);
526 fh2DeltaGRDenRPtRawMatch[fCentBin]->Fill(r,ptcorr,denRaw[i]);
527 }
528 }
529 }
530
e1239602 531 //Fill histograms for matched jets
e1239602 532 if(fMatch==1) {
15a3ff83 533 fh1PtMatch[fCentBin]->Fill(ptcorr);
534
535 //now get second derivative vs R and do final calculation
536 TArrayF num = jet1->GetGRNumeratorSub();
537 TArrayF den = jet1->GetGRDenominatorSub();
538 Printf("Got numerator size: nr: %d",nr);
539 cout << "size: " << num.GetSize() << endl;
540 if(num.GetSize()>0) {
541 for(Int_t i = 0; i<nr; i++) {
542 Double_t r = i*wr + 0.5*wr;
543 fh3DeltaGRNumRPtMatch[fCentBin]->Fill(num[i],r,ptcorr);
544 fh3DeltaGRDenRPtMatch[fCentBin]->Fill(den[i],r,ptcorr);
545 fh2DeltaGRNumRPtMatch[fCentBin]->Fill(r,ptcorr,num[i]);
546 fh2DeltaGRDenRPtMatch[fCentBin]->Fill(r,ptcorr,den[i]);
547 fh2DeltaGRNumRPtTrueMatch[fCentBin]->Fill(r,jet2->Pt(),num[i]);
548 fh2DeltaGRDenRPtTrueMatch[fCentBin]->Fill(r,jet2->Pt(),den[i]);
549
550 Double_t dGR = 0.;
551 fh2PtTrueDeltaGR[fCentBin]->Fill(jet2->Pt(),dGR);
552 }
e1239602 553 }
554 }
15a3ff83 555
e1239602 556 if(fCreateTree) {
557 fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
558 if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
559 else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
560 fArea = (Float_t)jet1->Area();
561 fAreaPhi = (Float_t)jet1->AreaPhi();
562 fAreaEta = (Float_t)jet1->AreaEta();
563 fRho = (Float_t)jetCont->GetRhoVal();
564 fRhoM = (Float_t)jetCont->GetRhoMassVal();
565 fNConst = (Int_t)jet1->GetNumberOfTracks();
566 fTreeJetBkg->Fill();
567 }
568 } //jet1 loop
569 }//jetCont
570
571
572 return kTRUE;
573}
574
15a3ff83 575//________________________________________________________________________
576Bool_t AliAnalysisTaskJetShapeGR::FillTrueJets() {
577
578 AliEmcalJet* jet1 = NULL;
579 AliJetContainer *jetCont = GetJetContainer(fContainerTrue);
580 if(!jetCont) {
581 Printf("cannot find container");
582 return kFALSE;
583 }
584
585 AliDebug(11,Form("NJets True: %d",jetCont->GetNJets()));
586
587 //create arrays
588 const Int_t nr = TMath::CeilNint(fMaxR/fDRStep);
589 TArrayF *fNum = new TArrayF(nr);
590 TArrayF *fDen = new TArrayF(nr);
591
592 if(jetCont) {
593 jetCont->ResetCurrentID();
594 while((jet1 = jetCont->GetNextAcceptJet())) {
595 fh1PtTrue[fCentBin]->Fill(jet1->Pt());
596
597 Double_t dev = CalcDeltaGR(jet1,fContainerTrue,fNum,fDen);//num,den);
598 for(Int_t i = 0; i<nr; i++) {
599 Double_t r = i*fDRStep + 0.5*fDRStep;
600 fh3DeltaGRNumRPtTrue[fCentBin]->Fill(fNum->At(i),r,jet1->Pt());
601 fh3DeltaGRDenRPtTrue[fCentBin]->Fill(fDen->At(i),r,jet1->Pt());
602 fh2DeltaGRNumRPtTrue[fCentBin]->Fill(r,jet1->Pt(),fNum->At(i));
603 fh2DeltaGRDenRPtTrue[fCentBin]->Fill(r,jet1->Pt(),fDen->At(i));
604 }
605 }
606 }
607 if(fNum) delete fNum;
608 if(fDen) delete fDen;
609 return kTRUE;
610}
e1239602 611
612//________________________________________________________________________
15a3ff83 613Double_t AliAnalysisTaskJetShapeGR::CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArrayF *fNum, TArrayF *fDen) { //Double_t *num, Double_t *den) {
e1239602 614 //Calculate G(R)
15a3ff83 615
616 //First clear the arrays
617 const Int_t nr = TMath::CeilNint(fMaxR/fDRStep);
618 for(Int_t i = 0; i<nr; i++) {
619 fNum->SetAt(0.,i);
620 fDen->SetAt(0.,i);
621 }
622
e1239602 623 AliJetContainer *jetCont = GetJetContainer(ic);
624 AliVParticle *vp1 = 0x0;
625 AliVParticle *vp2 = 0x0;
15a3ff83 626 Double_t A = 0.; Double_t B = 0.;
627 if(jet->GetNumberOfTracks()<2) return 0.;
628 // Printf("jet pt: %f nconst: %d",jet->Pt(),jet->GetNumberOfTracks());
e1239602 629 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
630 vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
15a3ff83 631 if(!vp1) continue;
632 for(Int_t j=i+1; j<jet->GetNumberOfTracks(); j++) {
e1239602 633 vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
15a3ff83 634 if(!vp2) continue;
635 Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
636 Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
637 // Printf("dr2: %f dEta: %f dphi: %f",dr2,vp1->Eta()-vp2->Eta(),dphi);
638 if(dr2>0.) {
639 for(Int_t k = 0; k<nr; k++) {
640 Double_t r = k*fDRStep + 0.5*fDRStep;
641 // Double_t x = jetCont->GetJetRadius()-TMath::Sqrt(dr2);
642 Double_t dr = TMath::Sqrt(dr2);
643 Double_t x = r-dr;
644 //noisy function
645 Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
646 //error function
647 Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
648
649 A = vp1->Pt()*vp2->Pt()*dr2*noise;
650 B = vp1->Pt()*vp2->Pt()*dr2*erf;
651 fNum->AddAt(fNum->At(k)+A,k);
652 fDen->AddAt(fDen->At(k)+B,k);
653 }
654 }
e1239602 655 }
656 }
15a3ff83 657
658 Double_t deltaGR = 0.;
659 if(B>0.) deltaGR = A/B; //useless
660 return deltaGR;
e1239602 661}
662
663//________________________________________________________________________
15a3ff83 664Double_t AliAnalysisTaskJetShapeGR::CalcGR(AliEmcalJet *jet, Int_t ic) {
e1239602 665 //Calculate G(R)
666 AliJetContainer *jetCont = GetJetContainer(ic);
667 AliVParticle *vp1 = 0x0;
668 AliVParticle *vp2 = 0x0;
669 Double_t gR = 0.;
670 Double_t wr = 0.04;
671 const Int_t nr = TMath::CeilNint(jetCont->GetJetRadius()/wr);
672 Double_t grArr[nr];
673 for(Int_t i = 0; i<nr; i++) {
674 grArr[i] = 0.;
15a3ff83 675 //Printf("bin up edge %d=%f",i,wr+i*wr);
e1239602 676 }
15a3ff83 677 // Printf("jet pt: %f nconst: %d",jet->Pt(),jet->GetNumberOfTracks());
e1239602 678 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
679 vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
680 for(Int_t j=i; j<jet->GetNumberOfTracks(); j++) {
681 vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
15a3ff83 682 Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
683 Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
e1239602 684 Int_t bin = TMath::FloorNint(TMath::Sqrt(dr2)/wr);
685 Double_t gr = vp1->Pt()*vp2->Pt()*dr2;
686 if(bin<nr) grArr[bin]+=gr;
687
688 if(TMath::Sqrt(dr2)<jetCont->GetJetRadius())
689 gR += gr;
690 }
691 }
e1239602 692 return gR;
693}
694
e1239602 695//________________________________________________________________________
696AliVParticle* AliAnalysisTaskJetShapeGR::GetEmbeddedConstituent(AliEmcalJet *jet) {
697
698 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
699 AliVParticle *vp = 0x0;
700 AliVParticle *vpe = 0x0; //embedded particle
701 Int_t nc = 0;
702 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
703 vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
704 if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
705 if(!vpe) vpe = vp;
706 else if(vp->Pt()>vpe->Pt()) vpe = vp;
707 nc++;
708 }
709 AliDebug(11,Form("Found %d embedded particles",nc));
710 return vpe;
711}
712
15a3ff83 713//________________________________________________________________________
714Double_t AliAnalysisTaskJetShapeGR::GetDeltaPhi(Double_t phi1,Double_t phi2) {
715 //
716 // Calculate azimuthal angle difference
717 //
718
719 Double_t dPhi = phi1-phi2;
720
721 if(dPhi <-1.*TMath::Pi()) dPhi += TMath::TwoPi();
722 if(dPhi > 1.*TMath::Pi()) dPhi -= TMath::TwoPi();
723
724 return dPhi;
725}
726
e1239602 727
728//________________________________________________________________________
729Bool_t AliAnalysisTaskJetShapeGR::RetrieveEventObjects() {
730 //
731 // retrieve event objects
732 //
733
734 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
735 return kFALSE;
736
737 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
738 jetCont->LoadRhoMass(InputEvent());
739
740 return kTRUE;
741}
742
743//_______________________________________________________________________
744void AliAnalysisTaskJetShapeGR::Terminate(Option_t *)
745{
746 // Called once at the end of the analysis.
747}
748