update to master versions
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetShapeDeriv.cxx
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
37 ClassImp(AliAnalysisTaskJetShapeDeriv)
38
39 //________________________________________________________________________
40 AliAnalysisTaskJetShapeDeriv::AliAnalysisTaskJetShapeDeriv() : 
41   AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeDeriv", kTRUE),
42   fContainerBase(0),
43   fContainerNoEmb(1),
44   fMinFractionShared(0),
45   fSingleTrackEmb(kFALSE),
46   fCreateTree(kFALSE),
47   fJetMassVarType(kMass),
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),
64   fh3MSubPtRawDRMatch(0x0),
65   fh3MSubPtTrueDR(0x0),
66   fh3MTruePtTrueDR(0x0),
67   fh3PtTrueDeltaMDR(0x0),
68   fh3PtTrueDeltaMRelDR(0x0),
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
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];
97
98   for (Int_t i = 0; i < fNcentBins; i++) {
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;
115   }
116
117   SetMakeGeneralHistograms(kTRUE);
118   if(fCreateTree) DefineOutput(2, TTree::Class());
119 }
120
121 //________________________________________________________________________
122 AliAnalysisTaskJetShapeDeriv::AliAnalysisTaskJetShapeDeriv(const char *name) : 
123   AliAnalysisTaskEmcalJet(name, kTRUE),  
124   fContainerBase(0),
125   fContainerNoEmb(1),
126   fMinFractionShared(0),
127   fSingleTrackEmb(kFALSE),
128   fCreateTree(kFALSE),
129   fJetMassVarType(kMass),
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),
146   fh3MSubPtRawDRMatch(0x0),
147   fh3MSubPtTrueDR(0x0),
148   fh3MTruePtTrueDR(0x0),
149   fh3PtTrueDeltaMDR(0x0),
150   fh3PtTrueDeltaMRelDR(0x0),
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
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];
179
180   for (Int_t i = 0; i < fNcentBins; i++) {
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;
197   }
198
199   SetMakeGeneralHistograms(kTRUE);
200   if(fCreateTree) DefineOutput(2, TTree::Class());
201 }
202
203 //________________________________________________________________________
204 AliAnalysisTaskJetShapeDeriv::~AliAnalysisTaskJetShapeDeriv()
205 {
206   // Destructor.
207 }
208
209 //________________________________________________________________________
210 void 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
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   }
231
232   Int_t nBinsDM  = 100;
233   Double_t minDM = -25.;
234   Double_t maxDM = 25.;
235   if(fJetMassVarType==kRatMPt) {
236     nBinsDM = 100;
237     minDM   = -0.5;
238     maxDM   = 0.5;
239   }
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.;
244
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
254   const Int_t nBinsSparse0 = 5;
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};
258
259   TString histName = "";
260   TString histTitle = "";
261   TString varName = "#it{M}_{jet}";
262   if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";
263
264   for (Int_t i = 0; i < fNcentBins; i++) {
265     histName = Form("fh2MSubMatch_%d",i);
266     histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
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);
271     histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
272     fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
273     fOutput->Add(fh2MSubPtRawAll[i]);
274
275     histName = Form("fh3MSubPtRawDRMatch_%d",i);
276     histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
277     fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
278     fOutput->Add(fh3MSubPtRawDRMatch[i]);
279
280     histName = Form("fh3MSubPtTrueDR_%d",i);
281     histTitle = Form("fh3MSubPtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
282     fh3MSubPtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
283     fOutput->Add(fh3MSubPtTrueDR[i]);
284
285     histName = Form("fh3MTruePtTrueDR_%d",i);
286     histTitle = Form("fh3MTruePtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
287     fh3MTruePtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
288     fOutput->Add(fh3MTruePtTrueDR[i]);
289
290     histName = Form("fh3PtTrueDeltaMDR_%d",i);
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);
293     fOutput->Add(fh3PtTrueDeltaMDR[i]);
294
295     histName = Form("fh3PtTrueDeltaMRelDR_%d",i);
296     histTitle = Form("fh3PtTrueDeltaMRelDR_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
297     fh3PtTrueDeltaMRelDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,400,-1.,3.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
298     fOutput->Add(fh3PtTrueDeltaMRelDR[i]);
299
300     histName = Form("fhnMassResponse_%d",i);
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());
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 //________________________________________________________________________
376 Bool_t AliAnalysisTaskJetShapeDeriv::Run()
377 {
378   // Run analysis code here, if needed. It will be executed before FillHistograms().
379
380   return kTRUE;
381 }
382
383 //________________________________________________________________________
384 Bool_t AliAnalysisTaskJetShapeDeriv::FillHistograms()
385 {
386   // Fill histograms.
387
388   AliEmcalJet* jet1  = NULL; //AA jet
389   AliEmcalJet *jet2  = NULL; //Embedded Pythia jet
390   //  AliEmcalJet *jet1T = NULL; //tagged AA jet
391   //  AliEmcalJet *jet2T = NULL; //tagged Pythia jet
392   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
393   fRho  = (Float_t)jetCont->GetRhoVal();
394   fRhoM = (Float_t)jetCont->GetRhoMassVal();
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
401   if(jetCont) {
402     jetCont->ResetCurrentID();
403     while((jet1 = jetCont->GetNextAcceptJet())) {
404       jet2 = NULL;
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
414       //Fill histograms for all AA jets
415       fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);
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
445       //      if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
446
447       //Fill histograms for matched jets
448       fh2MSubMatch[fCentBin]->Fill(var,fMatch);
449       if(fMatch==1) {
450         Double_t drToLJ = -1.;
451         if(jetL) drToLJ = jet1->DeltaR(jetL);
452         fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
453         if(jet2) {
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);
462         }
463       }
464
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 //________________________________________________________________________
483 AliVParticle* 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 //________________________________________________________________________
502 Bool_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 //_______________________________________________________________________
517 void AliAnalysisTaskJetShapeDeriv::Terminate(Option_t *) 
518 {
519   // Called once at the end of the analysis.
520 }
521