AliPhysics  ced2227 (ced2227)
AliAnalysisTaskJetShapeGR.cxx
Go to the documentation of this file.
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 
36 
38 
39 //________________________________________________________________________
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),
59  fDRStep(0.04),
60  fMaxR(2.),
61  fh2PtTrueDeltaGR(0x0),
62  fh2PtTrueDeltaGRRel(0x0),
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)
86 {
87  // Default constructor.
88 
89  fh2PtTrueDeltaGR = new TH2F*[fNcentBins];
90  fh2PtTrueDeltaGRRel = new TH2F*[fNcentBins];
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];
114 
115  for (Int_t i = 0; i < fNcentBins; i++) {
116  fh2PtTrueDeltaGR[i] = 0;
117  fh2PtTrueDeltaGRRel[i] = 0;
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;
141  }
142 
143  SetMakeGeneralHistograms(kTRUE);
144  DefineOutput(2, TTree::Class());
145 }
146 
147 //________________________________________________________________________
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),
167  fDRStep(0.04),
168  fMaxR(2.),
169  fh2PtTrueDeltaGR(0x0),
170  fh2PtTrueDeltaGRRel(0x0),
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)
194 {
195  // Standard constructor.
196 
199  fhnGRResponse = new THnSparse*[fNcentBins];
200  fh1PtTrue = new TH1F*[fNcentBins];
205  fh1PtRaw = new TH1F*[fNcentBins];
210  fh1PtRawMatch = new TH1F*[fNcentBins];
215  fh1PtMatch = new TH1F*[fNcentBins];
222 
223  for (Int_t i = 0; i < fNcentBins; i++) {
224  fh2PtTrueDeltaGR[i] = 0;
225  fh2PtTrueDeltaGRRel[i] = 0;
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;
242  fh1PtMatch[i] = 0;
243  fh3DeltaGRNumRPtMatch[i] = 0;
244  fh3DeltaGRDenRPtMatch[i] = 0;
245  fh2DeltaGRNumRPtMatch[i] = 0;
246  fh2DeltaGRDenRPtMatch[i] = 0;
249  }
250 
252  DefineOutput(2, TTree::Class());
253 }
254 
255 //________________________________________________________________________
257 {
258  // Destructor.
259 }
260 
261 //________________________________________________________________________
263 {
264  // Create user output.
265 
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 
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 
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++) {
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]);
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 
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 //________________________________________________________________________
460 {
461  // Run analysis code here, if needed. It will be executed before FillHistograms().
462 
463  return kTRUE;
464 }
465 
466 //________________________________________________________________________
468 {
469  // Fill histograms.
470  FillTrueJets();
471 
472  AliEmcalJet* jet1 = NULL;
473  AliEmcalJet *jet2 = NULL;
474  AliEmcalJet *jetS = NULL;
477  AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
478  if(jetCont && jetContS) {
479  jetCont->ResetCurrentID();
480  Double_t rmax = jetCont->GetJetRadius()+0.2;
481  Double_t wr = 0.04;
482  Int_t nr = TMath::CeilNint(rmax/wr);
483  while((jet1 = jetCont->GetNextAcceptJet())) {
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();
495  if(!jet2) continue;
496 
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 
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->GetShapeProperties()->GetGRNumerator();
514  TArrayF denRaw = jet1->GetShapeProperties()->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 
531  //Fill histograms for matched jets
532  if(fMatch==1) {
533  fh1PtMatch[fCentBin]->Fill(ptcorr);
534 
535  //now get second derivative vs R and do final calculation
536  TArrayF num = jet1->GetShapeProperties()->GetGRNumeratorSub();
537  TArrayF den = jet1->GetShapeProperties()->GetGRDenominatorSub();
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  }
551  }
552  }
553 
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 
573 //________________________________________________________________________
575 
576  AliEmcalJet* jet1 = NULL;
578  if(!jetCont)
579  return kFALSE;
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 
593  //Double_t dev = CalcDeltaGR(jet1,fContainerTrue,fNum,fDen);//num,den);
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 }
607 
608 //________________________________________________________________________
609 Double_t AliAnalysisTaskJetShapeGR::CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArrayF *fNum, TArrayF *fDen) { //Double_t *num, Double_t *den) {
610  //Calculate G(R)
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 
619  AliJetContainer *jetCont = GetJetContainer(ic);
620  AliVParticle *vp1 = 0x0;
621  AliVParticle *vp2 = 0x0;
622  Double_t A = 0.; Double_t B = 0.;
623  if(jet->GetNumberOfTracks()<2) return 0.;
624  for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
625  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
626  if(!vp1) continue;
627  for(Int_t j=i+1; j<jet->GetNumberOfTracks(); j++) {
628  vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
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;
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  }
649  }
650  }
651 
652  Double_t deltaGR = 0.;
653  if(B>0.) deltaGR = A/B; //useless
654  return deltaGR;
655 }
656 
657 //________________________________________________________________________
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[999] = {0.};
667 
668  for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
669  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
670  for(Int_t j=i; j<jet->GetNumberOfTracks(); j++) {
671  vp2 = static_cast<AliVParticle*>(jet->TrackAt(j, jetCont->GetParticleContainer()->GetArray()));
672  Double_t dphi = GetDeltaPhi(vp1->Phi(),vp2->Phi());
673  Double_t dr2 = (vp1->Eta()-vp2->Eta())*(vp1->Eta()-vp2->Eta()) + dphi*dphi;
674  Int_t bin = TMath::FloorNint(TMath::Sqrt(dr2)/wr);
675  Double_t gr = vp1->Pt()*vp2->Pt()*dr2;
676  if(bin<nr) grArr[bin]+=gr;
677 
678  if(TMath::Sqrt(dr2)<jetCont->GetJetRadius())
679  gR += gr;
680  }
681  }
682  return gR;
683 }
684 
685 //________________________________________________________________________
687 
689  AliVParticle *vp = 0x0;
690  AliVParticle *vpe = 0x0; //embedded particle
691  Int_t nc = 0;
692  for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
693  vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
694  if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
695  if(!vpe) vpe = vp;
696  else if(vp->Pt()>vpe->Pt()) vpe = vp;
697  nc++;
698  }
699  AliDebug(11,Form("Found %d embedded particles",nc));
700  return vpe;
701 }
702 
703 //________________________________________________________________________
705  //
706  // Calculate azimuthal angle difference
707  //
708 
709  Double_t dPhi = phi1-phi2;
710 
711  if(dPhi <-1.*TMath::Pi()) dPhi += TMath::TwoPi();
712  if(dPhi > 1.*TMath::Pi()) dPhi -= TMath::TwoPi();
713 
714  return dPhi;
715 }
716 
717 
718 //________________________________________________________________________
720  //
721  // retrieve event objects
722  //
723 
725  return kFALSE;
726 
728  jetCont->LoadRhoMass(InputEvent());
729 
730  return kTRUE;
731 }
732 
733 //_______________________________________________________________________
735 {
736  // Called once at the end of the analysis.
737 }
738 
Double_t Area() const
Definition: AliEmcalJet.h:130
Double_t GetRhoVal() const
TH2F ** fh2DeltaGRDenRPtRaw
Numerator of DeltaGR vs R vs Pt : filled with weights of sum.
double Double_t
Definition: External.C:58
TH2F ** fh2DeltaGRDenRPtMatch
Numerator of DeltaGR vs R vs Pt : filled with weights of sum.
TH3F ** fh3DeltaGRNumRPtRaw
bookkeep number of jets vs Pt
Definition: External.C:260
Definition: External.C:236
TH3F ** fh3DeltaGRDenRPtTrue
Numerator of DeltaGR vs R vs Pt.
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
TH2F ** fh2DeltaGRNumRPtTrueMatch
Denomerator of DeltaGR vs R vs Pt : filled with weights of sum.
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Double_t GetDeltaPhi(Double_t phi1, Double_t phi2)
TH2F ** fh2DeltaGRNumRPtRawMatch
Denomerator of DeltaGR vs R vs Pt.
Int_t fCentBin
!event centrality bin
TH3F ** fh3DeltaGRNumRPtTrue
bookkeep number of jets vs Pt
void LoadRhoMass(const AliVEvent *event)
Double_t E() const
Definition: AliEmcalJet.h:119
THnSparse ** fhnGRResponse
true jet pT vs (Msub - Mtrue)/Mtrue
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
AliVParticle * GetEmbeddedConstituent(AliEmcalJet *jet)
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t Px() const
Definition: AliEmcalJet.h:106
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
TH2F ** fh2DeltaGRNumRPtMatch
Denomerator of DeltaGR vs R vs Pt.
AliParticleContainer * GetParticleContainer() const
Double_t CalcDeltaGR(AliEmcalJet *jet, Int_t ic, TArrayF *fNum, TArrayF *fDen)
Double_t AreaEta() const
Definition: AliEmcalJet.h:132
TH2F ** fh2DeltaGRNumRPtRaw
Denomerator of DeltaGR vs R vs Pt.
int Int_t
Definition: External.C:63
TH2F ** fh2PtTrueDeltaGRRel
true jet pT vs (Msub - Mtrue)
float Float_t
Definition: External.C:68
TH3F ** fh3DeltaGRDenRPtRaw
Numerator of DeltaGR vs R vs Pt.
TH2F ** fh2DeltaGRDenRPtTrueMatch
Numerator of DeltaGR vs R vs Pt : filled with weights of sum.
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
TH3F ** fh3DeltaGRDenRPtMatch
Numerator of DeltaGR vs R vs Pt.
TH3F ** fh3DeltaGRNumRPtMatch
bookkeep number of jets vs Pt
Int_t GetNJets() const
TH1F ** fh1PtMatch
Denomerator of DeltaGR vs R vs Pt : filled with weights of sum.
Int_t fNcentBins
how many centrality bins
TLorentzVector * fJet1Vec
tree with jet and bkg variables
Double_t GetRhoMassVal() const
TH1F ** fh1PtRawMatch
Denomerator of DeltaGR vs R vs Pt : filled with weights of sum.
TH2F ** fh2DeltaGRDenRPtRawMatch
Numerator of DeltaGR vs R vs Pt : filled with weights of sum.
TH1F ** fh1PtTrue
Msub vs Mtrue vs PtCorr vs PtTrue.
AliEmcalJet * GetNextAcceptJet()
TH2F ** fh2DeltaGRNumRPtTrue
Denomerator of DeltaGR vs R vs Pt.
Double_t AreaPhi() const
Definition: AliEmcalJet.h:133
Double_t Pt() const
Definition: AliEmcalJet.h:109
Bool_t FillHistograms()
Function filling histograms.
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Double_t Pz() const
Definition: AliEmcalJet.h:108
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
TH1F ** fh1PtRaw
Denomerator of DeltaGR vs R vs Pt : filled with weights of sum.
TH3F ** fh3DeltaGRDenRPtRawMatch
Numerator of DeltaGR vs R vs Pt.
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
TH3F ** fh3DeltaGRNumRPtRawMatch
bookkeep number of jets vs Pt
Double_t CalcGR(AliEmcalJet *jet, Int_t ic)
TH2F ** fh2DeltaGRDenRPtTrue
Numerator of DeltaGR vs R vs Pt : filled with weights of sum.
Double_t M() const
Definition: AliEmcalJet.h:120
Container for jet within the EMCAL jet framework.
Definition: External.C:196