]> 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
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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();
15a3ff83 538 if(num.GetSize()>0) {
539 for(Int_t i = 0; i<nr; i++) {
540 Double_t r = i*wr + 0.5*wr;
541 fh3DeltaGRNumRPtMatch[fCentBin]->Fill(num[i],r,ptcorr);
542 fh3DeltaGRDenRPtMatch[fCentBin]->Fill(den[i],r,ptcorr);
543 fh2DeltaGRNumRPtMatch[fCentBin]->Fill(r,ptcorr,num[i]);
544 fh2DeltaGRDenRPtMatch[fCentBin]->Fill(r,ptcorr,den[i]);
545 fh2DeltaGRNumRPtTrueMatch[fCentBin]->Fill(r,jet2->Pt(),num[i]);
546 fh2DeltaGRDenRPtTrueMatch[fCentBin]->Fill(r,jet2->Pt(),den[i]);
547
548 Double_t dGR = 0.;
549 fh2PtTrueDeltaGR[fCentBin]->Fill(jet2->Pt(),dGR);
550 }
e1239602 551 }
552 }
15a3ff83 553
e1239602 554 if(fCreateTree) {
555 fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
556 if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
557 else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
558 fArea = (Float_t)jet1->Area();
559 fAreaPhi = (Float_t)jet1->AreaPhi();
560 fAreaEta = (Float_t)jet1->AreaEta();
561 fRho = (Float_t)jetCont->GetRhoVal();
562 fRhoM = (Float_t)jetCont->GetRhoMassVal();
563 fNConst = (Int_t)jet1->GetNumberOfTracks();
564 fTreeJetBkg->Fill();
565 }
566 } //jet1 loop
567 }//jetCont
568
569
570 return kTRUE;
571}
572
15a3ff83 573//________________________________________________________________________
574Bool_t AliAnalysisTaskJetShapeGR::FillTrueJets() {
575
576 AliEmcalJet* jet1 = NULL;
577 AliJetContainer *jetCont = GetJetContainer(fContainerTrue);
f831dc9f 578 if(!jetCont)
15a3ff83 579 return kFALSE;
15a3ff83 580
581 AliDebug(11,Form("NJets True: %d",jetCont->GetNJets()));
582
583 //create arrays
584 const Int_t nr = TMath::CeilNint(fMaxR/fDRStep);
585 TArrayF *fNum = new TArrayF(nr);
586 TArrayF *fDen = new TArrayF(nr);
587
588 if(jetCont) {
589 jetCont->ResetCurrentID();
590 while((jet1 = jetCont->GetNextAcceptJet())) {
591 fh1PtTrue[fCentBin]->Fill(jet1->Pt());
592
8d95aff1 593 //Double_t dev = CalcDeltaGR(jet1,fContainerTrue,fNum,fDen);//num,den);
15a3ff83 594 for(Int_t i = 0; i<nr; i++) {
595 Double_t r = i*fDRStep + 0.5*fDRStep;
596 fh3DeltaGRNumRPtTrue[fCentBin]->Fill(fNum->At(i),r,jet1->Pt());
597 fh3DeltaGRDenRPtTrue[fCentBin]->Fill(fDen->At(i),r,jet1->Pt());
598 fh2DeltaGRNumRPtTrue[fCentBin]->Fill(r,jet1->Pt(),fNum->At(i));
599 fh2DeltaGRDenRPtTrue[fCentBin]->Fill(r,jet1->Pt(),fDen->At(i));
600 }
601 }
602 }
603 if(fNum) delete fNum;
604 if(fDen) delete fDen;
605 return kTRUE;
606}
e1239602 607
608//________________________________________________________________________
15a3ff83 609Double_t AliAnalysisTaskJetShapeGR::CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArrayF *fNum, TArrayF *fDen) { //Double_t *num, Double_t *den) {
e1239602 610 //Calculate G(R)
15a3ff83 611
612 //First clear the arrays
613 const Int_t nr = TMath::CeilNint(fMaxR/fDRStep);
614 for(Int_t i = 0; i<nr; i++) {
615 fNum->SetAt(0.,i);
616 fDen->SetAt(0.,i);
617 }
618
e1239602 619 AliJetContainer *jetCont = GetJetContainer(ic);
620 AliVParticle *vp1 = 0x0;
621 AliVParticle *vp2 = 0x0;
15a3ff83 622 Double_t A = 0.; Double_t B = 0.;
623 if(jet->GetNumberOfTracks()<2) return 0.;
e1239602 624 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
625 vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
15a3ff83 626 if(!vp1) continue;
627 for(Int_t j=i+1; j<jet->GetNumberOfTracks(); j++) {
e1239602 628 vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
15a3ff83 629 if(!vp2) continue;
630 Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
631 Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
15a3ff83 632 if(dr2>0.) {
633 for(Int_t k = 0; k<nr; k++) {
634 Double_t r = k*fDRStep + 0.5*fDRStep;
635 // Double_t x = jetCont->GetJetRadius()-TMath::Sqrt(dr2);
636 Double_t dr = TMath::Sqrt(dr2);
637 Double_t x = r-dr;
638 //noisy function
639 Double_t noise = TMath::Exp(-x*x/(2*fDRStep*fDRStep))/(TMath::Sqrt(2.*TMath::Pi())*fDRStep);
640 //error function
641 Double_t erf = 0.5*(1.+TMath::Erf(x/(TMath::Sqrt(2.)*fDRStep)));
642
643 A = vp1->Pt()*vp2->Pt()*dr2*noise;
644 B = vp1->Pt()*vp2->Pt()*dr2*erf;
645 fNum->AddAt(fNum->At(k)+A,k);
646 fDen->AddAt(fDen->At(k)+B,k);
647 }
648 }
e1239602 649 }
650 }
15a3ff83 651
652 Double_t deltaGR = 0.;
653 if(B>0.) deltaGR = A/B; //useless
654 return deltaGR;
e1239602 655}
656
657//________________________________________________________________________
15a3ff83 658Double_t AliAnalysisTaskJetShapeGR::CalcGR(AliEmcalJet *jet, Int_t ic) {
e1239602 659 //Calculate G(R)
660 AliJetContainer *jetCont = GetJetContainer(ic);
661 AliVParticle *vp1 = 0x0;
662 AliVParticle *vp2 = 0x0;
663 Double_t gR = 0.;
664 Double_t wr = 0.04;
665 const Int_t nr = TMath::CeilNint(jetCont->GetJetRadius()/wr);
666 Double_t grArr[nr];
f831dc9f 667 for(Int_t i = 0; i<nr; i++)
e1239602 668 grArr[i] = 0.;
f831dc9f 669
e1239602 670 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
671 vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
672 for(Int_t j=i; j<jet->GetNumberOfTracks(); j++) {
673 vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
15a3ff83 674 Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
675 Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
e1239602 676 Int_t bin = TMath::FloorNint(TMath::Sqrt(dr2)/wr);
677 Double_t gr = vp1->Pt()*vp2->Pt()*dr2;
678 if(bin<nr) grArr[bin]+=gr;
679
680 if(TMath::Sqrt(dr2)<jetCont->GetJetRadius())
681 gR += gr;
682 }
683 }
e1239602 684 return gR;
685}
686
e1239602 687//________________________________________________________________________
688AliVParticle* AliAnalysisTaskJetShapeGR::GetEmbeddedConstituent(AliEmcalJet *jet) {
689
690 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
691 AliVParticle *vp = 0x0;
692 AliVParticle *vpe = 0x0; //embedded particle
693 Int_t nc = 0;
694 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
695 vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
696 if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
697 if(!vpe) vpe = vp;
698 else if(vp->Pt()>vpe->Pt()) vpe = vp;
699 nc++;
700 }
701 AliDebug(11,Form("Found %d embedded particles",nc));
702 return vpe;
703}
704
15a3ff83 705//________________________________________________________________________
706Double_t AliAnalysisTaskJetShapeGR::GetDeltaPhi(Double_t phi1,Double_t phi2) {
707 //
708 // Calculate azimuthal angle difference
709 //
710
711 Double_t dPhi = phi1-phi2;
712
713 if(dPhi <-1.*TMath::Pi()) dPhi += TMath::TwoPi();
714 if(dPhi > 1.*TMath::Pi()) dPhi -= TMath::TwoPi();
715
716 return dPhi;
717}
718
e1239602 719
720//________________________________________________________________________
721Bool_t AliAnalysisTaskJetShapeGR::RetrieveEventObjects() {
722 //
723 // retrieve event objects
724 //
725
726 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
727 return kFALSE;
728
729 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
730 jetCont->LoadRhoMass(InputEvent());
731
732 return kTRUE;
733}
734
735//_______________________________________________________________________
736void AliAnalysisTaskJetShapeGR::Terminate(Option_t *)
737{
738 // Called once at the end of the analysis.
739}
740