]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeDeriv.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
change binning delta histos jet shape
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetShapeDeriv.cxx
CommitLineData
f4b02da3 1//
2// Do subtraction for jet shapes using derivatives arXiv:1211:2811
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 "AliAnalysisTaskJetShapeDeriv.h"
36
37ClassImp(AliAnalysisTaskJetShapeDeriv)
38
39//________________________________________________________________________
40AliAnalysisTaskJetShapeDeriv::AliAnalysisTaskJetShapeDeriv() :
41 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeDeriv", kTRUE),
42 fContainerBase(0),
a5a1c51a 43 fContainerNoEmb(1),
f4b02da3 44 fMinFractionShared(0),
45 fSingleTrackEmb(kFALSE),
46 fCreateTree(kFALSE),
d3660d18 47 fJetMassVarType(kMass),
f4b02da3 48 fTreeJetBkg(),
49 fJet1Vec(new TLorentzVector()),
50 fJet2Vec(new TLorentzVector()),
51 fArea(0),
52 fAreaPhi(0),
53 fAreaEta(0),
54 fRho(0),
55 fRhoM(0),
56 fNConst(0),
57 fM1st(0),
58 fM2nd(0),
59 fDeriv1st(0),
60 fDeriv2nd(0),
61 fMatch(0),
62 fh2MSubMatch(0x0),
63 fh2MSubPtRawAll(0x0),
b09cf300 64 fh3MSubPtRawDRMatch(0x0),
65 fh3MSubPtTrueDR(0x0),
66 fh3MTruePtTrueDR(0x0),
67 fh3PtTrueDeltaMDR(0x0),
68 fh3PtTrueDeltaMRelDR(0x0),
f4b02da3 69 fhnMassResponse(0x0),
70 fh2PtTrueSubFacV1(0x0),
71 fh2PtRawSubFacV1(0x0),
72 fh2PtCorrSubFacV1(0x0),
73 fh2NConstSubFacV1(0x0),
74 fh2PtTrueSubFacV2(0x0),
75 fh2PtRawSubFacV2(0x0),
76 fh2PtCorrSubFacV2(0x0),
77 fh2NConstSubFacV2(0x0)
78{
79 // Default constructor.
80
b09cf300 81 fh2MSubMatch = new TH2F*[fNcentBins];
82 fh2MSubPtRawAll = new TH2F*[fNcentBins];
83 fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
84 fh3MSubPtTrueDR = new TH3F*[fNcentBins];
85 fh3MTruePtTrueDR = new TH3F*[fNcentBins];
86 fh3PtTrueDeltaMDR = new TH3F*[fNcentBins];
87 fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
88 fhnMassResponse = new THnSparse*[fNcentBins];
89 fh2PtTrueSubFacV1 = new TH2F*[fNcentBins];
90 fh2PtRawSubFacV1 = new TH2F*[fNcentBins];
91 fh2PtCorrSubFacV1 = new TH2F*[fNcentBins];
92 fh2NConstSubFacV1 = new TH2F*[fNcentBins];
93 fh2PtTrueSubFacV2 = new TH2F*[fNcentBins];
94 fh2PtRawSubFacV2 = new TH2F*[fNcentBins];
95 fh2PtCorrSubFacV2 = new TH2F*[fNcentBins];
96 fh2NConstSubFacV2 = new TH2F*[fNcentBins];
f4b02da3 97
98 for (Int_t i = 0; i < fNcentBins; i++) {
b09cf300 99 fh2MSubMatch[i] = 0;
100 fh2MSubPtRawAll[i] = 0;
101 fh3MSubPtRawDRMatch[i] = 0;
102 fh3MSubPtTrueDR[i] = 0;
103 fh3MTruePtTrueDR[i] = 0;
104 fh3PtTrueDeltaMDR[i] = 0;
105 fh3PtTrueDeltaMRelDR[i] = 0;
106 fhnMassResponse[i] = 0;
107 fh2PtTrueSubFacV1[i] = 0;
108 fh2PtRawSubFacV1[i] = 0;
109 fh2PtCorrSubFacV1[i] = 0;
110 fh2NConstSubFacV1[i] = 0;
111 fh2PtTrueSubFacV2[i] = 0;
112 fh2PtRawSubFacV2[i] = 0;
113 fh2PtCorrSubFacV2[i] = 0;
114 fh2NConstSubFacV2[i] = 0;
f4b02da3 115 }
116
117 SetMakeGeneralHistograms(kTRUE);
713c5683 118 if(fCreateTree) DefineOutput(2, TTree::Class());
f4b02da3 119}
120
121//________________________________________________________________________
122AliAnalysisTaskJetShapeDeriv::AliAnalysisTaskJetShapeDeriv(const char *name) :
123 AliAnalysisTaskEmcalJet(name, kTRUE),
124 fContainerBase(0),
a5a1c51a 125 fContainerNoEmb(1),
f4b02da3 126 fMinFractionShared(0),
127 fSingleTrackEmb(kFALSE),
128 fCreateTree(kFALSE),
d3660d18 129 fJetMassVarType(kMass),
f4b02da3 130 fTreeJetBkg(0),
131 fJet1Vec(new TLorentzVector()),
132 fJet2Vec(new TLorentzVector()),
133 fArea(0),
134 fAreaPhi(0),
135 fAreaEta(0),
136 fRho(0),
137 fRhoM(0),
138 fNConst(0),
139 fM1st(0),
140 fM2nd(0),
141 fDeriv1st(0),
142 fDeriv2nd(0),
143 fMatch(0),
144 fh2MSubMatch(0x0),
145 fh2MSubPtRawAll(0x0),
b09cf300 146 fh3MSubPtRawDRMatch(0x0),
147 fh3MSubPtTrueDR(0x0),
148 fh3MTruePtTrueDR(0x0),
149 fh3PtTrueDeltaMDR(0x0),
150 fh3PtTrueDeltaMRelDR(0x0),
f4b02da3 151 fhnMassResponse(0x0),
152 fh2PtTrueSubFacV1(0x0),
153 fh2PtRawSubFacV1(0x0),
154 fh2PtCorrSubFacV1(0x0),
155 fh2NConstSubFacV1(0x0),
156 fh2PtTrueSubFacV2(0x0),
157 fh2PtRawSubFacV2(0x0),
158 fh2PtCorrSubFacV2(0x0),
159 fh2NConstSubFacV2(0x0)
160{
161 // Standard constructor.
162
b09cf300 163 fh2MSubMatch = new TH2F*[fNcentBins];
164 fh2MSubPtRawAll = new TH2F*[fNcentBins];
165 fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
166 fh3MSubPtTrueDR = new TH3F*[fNcentBins];
167 fh3MTruePtTrueDR = new TH3F*[fNcentBins];
168 fh3PtTrueDeltaMDR = new TH3F*[fNcentBins];
169 fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
170 fhnMassResponse = new THnSparse*[fNcentBins];
171 fh2PtTrueSubFacV1 = new TH2F*[fNcentBins];
172 fh2PtRawSubFacV1 = new TH2F*[fNcentBins];
173 fh2PtCorrSubFacV1 = new TH2F*[fNcentBins];
174 fh2NConstSubFacV1 = new TH2F*[fNcentBins];
175 fh2PtTrueSubFacV2 = new TH2F*[fNcentBins];
176 fh2PtRawSubFacV2 = new TH2F*[fNcentBins];
177 fh2PtCorrSubFacV2 = new TH2F*[fNcentBins];
178 fh2NConstSubFacV2 = new TH2F*[fNcentBins];
f4b02da3 179
180 for (Int_t i = 0; i < fNcentBins; i++) {
b09cf300 181 fh2MSubMatch[i] = 0;
182 fh2MSubPtRawAll[i] = 0;
183 fh3MSubPtRawDRMatch[i] = 0;
184 fh3MSubPtTrueDR[i] = 0;
185 fh3MTruePtTrueDR[i] = 0;
186 fh3PtTrueDeltaMDR[i] = 0;
187 fh3PtTrueDeltaMRelDR[i] = 0;
188 fhnMassResponse[i] = 0;
189 fh2PtTrueSubFacV1[i] = 0;
190 fh2PtRawSubFacV1[i] = 0;
191 fh2PtCorrSubFacV1[i] = 0;
192 fh2NConstSubFacV1[i] = 0;
193 fh2PtTrueSubFacV2[i] = 0;
194 fh2PtRawSubFacV2[i] = 0;
195 fh2PtCorrSubFacV2[i] = 0;
196 fh2NConstSubFacV2[i] = 0;
f4b02da3 197 }
198
199 SetMakeGeneralHistograms(kTRUE);
713c5683 200 if(fCreateTree) DefineOutput(2, TTree::Class());
f4b02da3 201}
202
203//________________________________________________________________________
204AliAnalysisTaskJetShapeDeriv::~AliAnalysisTaskJetShapeDeriv()
205{
206 // Destructor.
207}
208
209//________________________________________________________________________
210void AliAnalysisTaskJetShapeDeriv::UserCreateOutputObjects()
211{
212 // Create user output.
213
214 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
215
216 Bool_t oldStatus = TH1::AddDirectoryStatus();
217 TH1::AddDirectory(kFALSE);
218
219 const Int_t nBinsPt = 200;
220 const Double_t minPt = -50.;
221 const Double_t maxPt = 150.;
222
d3660d18 223 Int_t nBinsM = 100;
224 Double_t minM = -20.;
225 Double_t maxM = 80.;
226 if(fJetMassVarType==kRatMPt) {
227 nBinsM = 100;
228 minM = -0.2;
229 maxM = 0.8;
230 }
f4b02da3 231
d3660d18 232 Int_t nBinsDM = 100;
a3c67475 233 Double_t minDM = -25.;
234 Double_t maxDM = 25.;
d3660d18 235 if(fJetMassVarType==kRatMPt) {
236 nBinsDM = 100;
237 minDM = -0.5;
238 maxDM = 0.5;
239 }
a5a1c51a 240
241 const Int_t nBinsDRToLJ = 20; //distance to leading jet in Pb-Pb only event
242 const Double_t minDRToLJ = 0.;
243 const Double_t maxDRToLJ = 1.;
8d95aff1 244
f4b02da3 245 const Int_t nBinsV1 = 60;
246 const Double_t minV1 = -60.;
247 const Double_t maxV1 = 0.;
248
249 const Int_t nBinsV2 = 60;
250 const Double_t minV2 = -30.;
251 const Double_t maxV2 = 0.;
252
253 //Binning for THnSparse
8d95aff1 254 const Int_t nBinsSparse0 = 5;
a5a1c51a 255 const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsDRToLJ};
256 const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt, minDRToLJ};
257 const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt, maxDRToLJ};
f4b02da3 258
259 TString histName = "";
260 TString histTitle = "";
d3660d18 261 TString varName = "#it{M}_{jet}";
262 if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";
f4b02da3 263
264 for (Int_t i = 0; i < fNcentBins; i++) {
265 histName = Form("fh2MSubMatch_%d",i);
d3660d18 266 histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
f4b02da3 267 fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
268 fOutput->Add(fh2MSubMatch[i]);
269
270 histName = Form("fh2MSubPtRawAll_%d",i);
d3660d18 271 histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
f4b02da3 272 fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
273 fOutput->Add(fh2MSubPtRawAll[i]);
274
b09cf300 275 histName = Form("fh3MSubPtRawDRMatch_%d",i);
d3660d18 276 histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300 277 fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
278 fOutput->Add(fh3MSubPtRawDRMatch[i]);
f4b02da3 279
b09cf300 280 histName = Form("fh3MSubPtTrueDR_%d",i);
d3660d18 281 histTitle = Form("fh3MSubPtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300 282 fh3MSubPtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
283 fOutput->Add(fh3MSubPtTrueDR[i]);
f4b02da3 284
b09cf300 285 histName = Form("fh3MTruePtTrueDR_%d",i);
d3660d18 286 histTitle = Form("fh3MTruePtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300 287 fh3MTruePtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
288 fOutput->Add(fh3MTruePtTrueDR[i]);
f4b02da3 289
b09cf300 290 histName = Form("fh3PtTrueDeltaMDR_%d",i);
d3660d18 291 histTitle = Form("fh3PtTrueDeltaMDR_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
292 fh3PtTrueDeltaMDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
b09cf300 293 fOutput->Add(fh3PtTrueDeltaMDR[i]);
f4b02da3 294
b09cf300 295 histName = Form("fh3PtTrueDeltaMRelDR_%d",i);
d3660d18 296 histTitle = Form("fh3PtTrueDeltaMRelDR_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
a3c67475 297 fh3PtTrueDeltaMRelDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,400,-1.,3.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
b09cf300 298 fOutput->Add(fh3PtTrueDeltaMRelDR[i]);
0947b103 299
f4b02da3 300 histName = Form("fhnMassResponse_%d",i);
d3660d18 301 histTitle = Form("fhnMassResponse_%d;%s sub;%s true;#it{p}_{T,sub};#it{p}_{T,true};#Delta R_{LJ}",i,varName.Data(),varName.Data());
f4b02da3 302 fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
303 fOutput->Add(fhnMassResponse[i]);
304
305 //derivative histograms
306 histName = Form("fh2PtTrueSubFacV1_%d",i);
307 histTitle = Form("fh2PtTrueSubFacV1_%d;#it{p}_{T,true};-(#rho+#rho_{m})V_{1}",i);
308 fh2PtTrueSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV1,minV1,maxV1);
309
310 histName = Form("fh2PtRawSubFacV1_%d",i);
311 histTitle = Form("fh2PtRawSubFacV1_%d;#it{p}_{T,raw};-(#rho+#rho_{m})V_{1}",i);
312 fh2PtRawSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV1,minV1,maxV1);
313
314 histName = Form("fh2PtCorrSubFacV1_%d",i);
315 histTitle = Form("fh2PtCorrSubFacV1_%d;#it{p}_{T,corr};-(#rho+#rho_{m})V_{1}",i);
316 fh2PtCorrSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV1,minV1,maxV1);
317
318 histName = Form("fh2NConstSubFacV1_%d",i);
319 histTitle = Form("fh2NConstSubFacV1_%d;#it{N}_{const};-(#rho+#rho_{m})V_{1}",i);
320 fh2NConstSubFacV1[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,200.);
321
322 histName = Form("fh2PtTrueSubFacV2_%d",i);
323 histTitle = Form("fh2PtTrueSubFacV2_%d;#it{p}_{T,true};0.5(#rho+#rho_{m})^{2}V_{2}",i);
324 fh2PtTrueSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV2,minV2,maxV2);
325
326 histName = Form("fh2PtRawSubFacV2_%d",i);
327 histTitle = Form("fh2PtRawSubFacV2_%d;#it{p}_{T,raw};0.5(#rho+#rho_{m})^{2}V_{2}",i);
328 fh2PtRawSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV2,minV2,maxV2);
329
330 histName = Form("fh2PtCorrSubFacV2_%d",i);
331 histTitle = Form("fh2PtCorrSubFacV2_%d;#it{p}_{T,corr};0.5(#rho+#rho_{m})^{2}V_{2}",i);
332 fh2PtCorrSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsV2,minV2,maxV2);
333
334 histName = Form("fh2NConstSubFacV2_%d",i);
335 histTitle = Form("fh2NConstSubFacV2_%d;#it{N}_{const};0.5(#rho+#rho_{m})^{2}V_{2}",i);
336 fh2NConstSubFacV2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,100,0.,200.);
337
338 }
339
340 // =========== Switch on Sumw2 for all histos ===========
341 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
342 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
343 if (h1){
344 h1->Sumw2();
345 continue;
346 }
347 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
348 if(hn)hn->Sumw2();
349 }
350
351 TH1::AddDirectory(oldStatus);
352
353 // Create a tree.
354 if(fCreateTree) {
355 fTreeJetBkg = new TTree("fTreeJetBkg", "fTreeJetBkg");
356 fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
357 fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
358 fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
359 fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
360 fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
361 fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
362 fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
363 fTreeJetBkg->Branch("fNConst",&fNConst,"fNConst/I");
364 fTreeJetBkg->Branch("fM1st",&fM1st,"fM1st/F");
365 fTreeJetBkg->Branch("fM2nd",&fM2nd,"fM2nd/F");
366 fTreeJetBkg->Branch("fDeriv1st",&fDeriv1st,"fDeriv1st/F");
367 fTreeJetBkg->Branch("fDeriv2nd",&fDeriv2nd,"fDeriv2nd/F");
368 fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
369 }
370
371 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
372 if(fCreateTree) PostData(2, fTreeJetBkg);
373}
374
375//________________________________________________________________________
376Bool_t AliAnalysisTaskJetShapeDeriv::Run()
377{
378 // Run analysis code here, if needed. It will be executed before FillHistograms().
379
380 return kTRUE;
381}
382
383//________________________________________________________________________
384Bool_t AliAnalysisTaskJetShapeDeriv::FillHistograms()
385{
386 // Fill histograms.
387
8d95aff1 388 AliEmcalJet* jet1 = NULL; //AA jet
389 AliEmcalJet *jet2 = NULL; //Embedded Pythia jet
a5a1c51a 390 // AliEmcalJet *jet1T = NULL; //tagged AA jet
8d95aff1 391 // AliEmcalJet *jet2T = NULL; //tagged Pythia jet
f4b02da3 392 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
393 fRho = (Float_t)jetCont->GetRhoVal();
394 fRhoM = (Float_t)jetCont->GetRhoMassVal();
a5a1c51a 395
396 //Get leading jet in Pb-Pb event without embedded objects
397 AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
398 AliEmcalJet *jetL = NULL;
399 if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
400
f4b02da3 401 if(jetCont) {
402 jetCont->ResetCurrentID();
403 while((jet1 = jetCont->GetNextAcceptJet())) {
8d95aff1 404 jet2 = NULL;
d3660d18 405
406 Double_t mjet1 = jet1->GetSecondOrderSubtracted();
407 Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
408 Double_t var = mjet1;
409 if(fJetMassVarType==kRatMPt) {
410 if(ptjet1>0. || ptjet1<0.) var = mjet1/ptjet1;
411 else var = -999.;
412 }
413
f4b02da3 414 //Fill histograms for all AA jets
d3660d18 415 fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);
f4b02da3 416 fh2PtRawSubFacV1[fCentBin]->Fill(jet1->Pt(),-1.*(fRho+fRhoM)*jet1->GetFirstDerivative());
417 fh2PtCorrSubFacV1[fCentBin]->Fill(jet1->Pt()-fRho*jet1->Area(),-1.*(fRho+fRhoM)*jet1->GetFirstDerivative());
418 fh2NConstSubFacV1[fCentBin]->Fill(jet1->GetNumberOfTracks(),-1.*(fRho+fRhoM)*jet1->GetFirstDerivative());
419 fh2PtRawSubFacV2[fCentBin]->Fill(jet1->Pt(),0.5*(fRho+fRhoM)*(fRho+fRhoM)*jet1->GetSecondDerivative());
420 fh2PtCorrSubFacV2[fCentBin]->Fill(jet1->Pt()-fRho*jet1->Area(),0.5*(fRho+fRhoM)*(fRho+fRhoM)*jet1->GetSecondDerivative());
421 fh2NConstSubFacV2[fCentBin]->Fill(jet1->GetNumberOfTracks(),0.5*(fRho+fRhoM)*(fRho+fRhoM)*jet1->GetSecondDerivative());
422
423 Double_t fraction = 0.;
424 fMatch = 0;
425 fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
426 if(fSingleTrackEmb) {
427 AliVParticle *vp = GetEmbeddedConstituent(jet1);
428 if(vp) {
429 fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
430 fMatch = 1;
431 }
432 } else {
433 jet2 = jet1->ClosestJet();
434 fraction = jetCont->GetFractionSharedPt(jet1);
435 fMatch = 1;
436 if(fMinFractionShared>0.) {
437 if(fraction>fMinFractionShared) {
438 fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
439 fMatch = 1;
440 } else
441 fMatch = 0;
442 }
443 }
444
8d95aff1 445 // if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
446
f4b02da3 447 //Fill histograms for matched jets
d3660d18 448 fh2MSubMatch[fCentBin]->Fill(var,fMatch);
f4b02da3 449 if(fMatch==1) {
b09cf300 450 Double_t drToLJ = -1.;
451 if(jetL) drToLJ = jet1->DeltaR(jetL);
d3660d18 452 fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
f4b02da3 453 if(jet2) {
d3660d18 454 Double_t var2 = -999.;
455 if(jet2->Pt()>0. || jet2->Pt()<0.) var2 = jet2->M()/jet2->Pt();
456 fh3MSubPtTrueDR[fCentBin]->Fill(var,jet2->Pt(),drToLJ);
457 fh3MTruePtTrueDR[fCentBin]->Fill(var2,jet2->Pt(),drToLJ);
458 fh3PtTrueDeltaMDR[fCentBin]->Fill(jet2->Pt(),var-var2,drToLJ);
459 if(jet2->M()>0.) fh3PtTrueDeltaMRelDR[fCentBin]->Fill(jet2->Pt(),(var-var2)/var2,drToLJ);
460 Double_t varsp[5] = {var,var2,ptjet1,jet2->Pt(),drToLJ};
461 fhnMassResponse[fCentBin]->Fill(varsp);
f4b02da3 462 }
463 }
464
f4b02da3 465 if(fCreateTree) {
466 fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
467 fArea = (Float_t)jet1->Area();
468 fAreaPhi = (Float_t)jet1->AreaPhi();
469 fAreaEta = (Float_t)jet1->AreaEta();
470 fNConst = (Int_t)jet1->GetNumberOfTracks();
471 fM1st = (Float_t)jet1->GetFirstOrderSubtracted();
472 fM2nd = (Float_t)jet1->GetSecondOrderSubtracted();
473 fDeriv1st = (Float_t)jet1->GetFirstDerivative();
474 fDeriv2nd = (Float_t)jet1->GetSecondDerivative();
475 fTreeJetBkg->Fill();
476 }
477 }
478 }
479 return kTRUE;
480}
481
482//________________________________________________________________________
483AliVParticle* AliAnalysisTaskJetShapeDeriv::GetEmbeddedConstituent(AliEmcalJet *jet) {
484
485 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
486 AliVParticle *vp = 0x0;
487 AliVParticle *vpe = 0x0; //embedded particle
488 Int_t nc = 0;
489 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
490 vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
491 if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
492 if(!vpe) vpe = vp;
493 else if(vp->Pt()>vpe->Pt()) vpe = vp;
494 nc++;
495 }
496 AliDebug(11,Form("Found %d embedded particles",nc));
497 return vpe;
498}
499
500
501//________________________________________________________________________
502Bool_t AliAnalysisTaskJetShapeDeriv::RetrieveEventObjects() {
503 //
504 // retrieve event objects
505 //
506
507 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
508 return kFALSE;
509
510 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
511 jetCont->LoadRhoMass(InputEvent());
512
513 return kTRUE;
514}
515
516//_______________________________________________________________________
517void AliAnalysisTaskJetShapeDeriv::Terminate(Option_t *)
518{
519 // Called once at the end of the analysis.
520}
521