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