update to master versions
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskJetShapeConst.cxx
CommitLineData
f4b02da3 1//
2// Do subtraction for jet shapes using constituents arXiv:1403.3108
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 "AliAnalysisTaskJetShapeConst.h"
36
37ClassImp(AliAnalysisTaskJetShapeConst)
38
39//________________________________________________________________________
40AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst() :
41 AliAnalysisTaskEmcalJet("AliAnalysisTaskJetShapeConst", kTRUE),
42 fContainerBase(0),
43 fContainerSub(1),
a5a1c51a 44 fContainerNoEmb(2),
f4b02da3 45 fMinFractionShared(0),
46 fSingleTrackEmb(kFALSE),
47 fCreateTree(kFALSE),
d3660d18 48 fJetMassVarType(kMass),
f4b02da3 49 fTreeJetBkg(),
50 fJet1Vec(new TLorentzVector()),
51 fJet2Vec(new TLorentzVector()),
52 fJetSubVec(new TLorentzVector()),
53 fArea(0),
54 fAreaPhi(0),
55 fAreaEta(0),
56 fRho(0),
57 fRhoM(0),
58 fNConst(0),
59 fMatch(0),
60 fh2MSubMatch(0x0),
61 fh2MSubPtRawAll(0x0),
b09cf300 62 fh3MSubPtRawDRMatch(0x0),
63 fh3MSubPtTrueDR(0x0),
64 fh3MTruePtTrueDR(0x0),
65 fh3PtTrueDeltaMDR(0x0),
66 fh3PtTrueDeltaMRelDR(0x0),
f4b02da3 67 fhnMassResponse(0x0)
68{
69 // Default constructor.
70
b09cf300 71 fh2MSubMatch = new TH2F*[fNcentBins];
72 fh2MSubPtRawAll = new TH2F*[fNcentBins];
73 fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
74 fh3MSubPtTrueDR = new TH3F*[fNcentBins];
75 fh3MTruePtTrueDR = new TH3F*[fNcentBins];
76 fh3PtTrueDeltaMDR = new TH3F*[fNcentBins];
77 fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
78 fhnMassResponse = new THnSparse*[fNcentBins];
f4b02da3 79
80 for (Int_t i = 0; i < fNcentBins; i++) {
b09cf300 81 fh2MSubMatch[i] = 0;
82 fh2MSubPtRawAll[i] = 0;
83 fh3MSubPtRawDRMatch[i] = 0;
84 fh3MSubPtTrueDR[i] = 0;
85 fh3MTruePtTrueDR[i] = 0;
86 fh3PtTrueDeltaMDR[i] = 0;
87 fh3PtTrueDeltaMRelDR[i] = 0;
88 fhnMassResponse[i] = 0;
f4b02da3 89 }
90
91 SetMakeGeneralHistograms(kTRUE);
713c5683 92 if(fCreateTree) DefineOutput(2, TTree::Class());
f4b02da3 93}
94
95//________________________________________________________________________
96AliAnalysisTaskJetShapeConst::AliAnalysisTaskJetShapeConst(const char *name) :
97 AliAnalysisTaskEmcalJet(name, kTRUE),
98 fContainerBase(0),
99 fContainerSub(1),
a5a1c51a 100 fContainerNoEmb(2),
f4b02da3 101 fMinFractionShared(0),
102 fSingleTrackEmb(kFALSE),
103 fCreateTree(kFALSE),
d3660d18 104 fJetMassVarType(kMass),
f4b02da3 105 fTreeJetBkg(0),
106 fJet1Vec(new TLorentzVector()),
107 fJet2Vec(new TLorentzVector()),
108 fJetSubVec(new TLorentzVector()),
109 fArea(0),
110 fAreaPhi(0),
111 fAreaEta(0),
112 fRho(0),
113 fRhoM(0),
114 fNConst(0),
115 fMatch(0),
116 fh2MSubMatch(0x0),
117 fh2MSubPtRawAll(0x0),
b09cf300 118 fh3MSubPtRawDRMatch(0x0),
119 fh3MSubPtTrueDR(0x0),
120 fh3MTruePtTrueDR(0x0),
121 fh3PtTrueDeltaMDR(0x0),
122 fh3PtTrueDeltaMRelDR(0x0),
f4b02da3 123 fhnMassResponse(0x0)
124{
125 // Standard constructor.
126
b09cf300 127 fh2MSubMatch = new TH2F*[fNcentBins];
128 fh2MSubPtRawAll = new TH2F*[fNcentBins];
129 fh3MSubPtRawDRMatch = new TH3F*[fNcentBins];
130 fh3MSubPtTrueDR = new TH3F*[fNcentBins];
131 fh3MTruePtTrueDR = new TH3F*[fNcentBins];
132 fh3PtTrueDeltaMDR = new TH3F*[fNcentBins];
133 fh3PtTrueDeltaMRelDR = new TH3F*[fNcentBins];
134 fhnMassResponse = new THnSparse*[fNcentBins];
f4b02da3 135
136 for (Int_t i = 0; i < fNcentBins; i++) {
b09cf300 137 fh2MSubMatch[i] = 0;
138 fh2MSubPtRawAll[i] = 0;
139 fh3MSubPtRawDRMatch[i] = 0;
140 fh3MSubPtTrueDR[i] = 0;
141 fh3MTruePtTrueDR[i] = 0;
142 fh3PtTrueDeltaMDR[i] = 0;
143 fh3PtTrueDeltaMRelDR[i] = 0;
144 fhnMassResponse[i] = 0;
f4b02da3 145 }
146
147 SetMakeGeneralHistograms(kTRUE);
713c5683 148 if(fCreateTree) DefineOutput(2, TTree::Class());
f4b02da3 149}
150
151//________________________________________________________________________
152AliAnalysisTaskJetShapeConst::~AliAnalysisTaskJetShapeConst()
153{
154 // Destructor.
155}
156
157//________________________________________________________________________
158void AliAnalysisTaskJetShapeConst::UserCreateOutputObjects()
159{
160 // Create user output.
161
162 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
163
164 Bool_t oldStatus = TH1::AddDirectoryStatus();
165 TH1::AddDirectory(kFALSE);
166
167 const Int_t nBinsPt = 200;
168 const Double_t minPt = -50.;
169 const Double_t maxPt = 150.;
170
d3660d18 171 Int_t nBinsM = 100;
172 Double_t minM = -20.;
173 Double_t maxM = 80.;
174 if(fJetMassVarType==kRatMPt) {
175 nBinsM = 100;
176 minM = -0.2;
177 maxM = 0.8;
178 }
f4b02da3 179
d3660d18 180 Int_t nBinsDM = 100;
5ca97880 181 Double_t minDM = -25.;
182 Double_t maxDM = 25.;
d3660d18 183 if(fJetMassVarType==kRatMPt) {
184 nBinsDM = 100;
185 minDM = -0.5;
186 maxDM = 0.5;
187 }
a5a1c51a 188
189 const Int_t nBinsDRToLJ = 20; //distance to leading jet in Pb-Pb only event
190 const Double_t minDRToLJ = 0.;
191 const Double_t maxDRToLJ = 1.;
8d95aff1 192
f4b02da3 193 //Binning for THnSparse
8d95aff1 194 const Int_t nBinsSparse0 = 5;
a5a1c51a 195 const Int_t nBins0[nBinsSparse0] = {nBinsM,nBinsM,nBinsPt,nBinsPt,nBinsDRToLJ};
196 const Double_t xmin0[nBinsSparse0] = { minM, minM, minPt, minPt, minDRToLJ};
197 const Double_t xmax0[nBinsSparse0] = { maxM, maxM, maxPt, maxPt, maxDRToLJ};
f4b02da3 198
199 TString histName = "";
200 TString histTitle = "";
d3660d18 201 TString varName = "#it{M}_{jet}";
202 if(fJetMassVarType==kRatMPt) varName = "#it{M}_{jet}/#it{p}_{T,jet}";
203
f4b02da3 204 for (Int_t i = 0; i < fNcentBins; i++) {
205 histName = Form("fh2MSubMatch_%d",i);
d3660d18 206 histTitle = Form("fh2MSubMatch_%d;%s;match",i,varName.Data());
f4b02da3 207 fh2MSubMatch[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,2,-0.5,1.5);
208 fOutput->Add(fh2MSubMatch[i]);
209
210 histName = Form("fh2MSubPtRawAll_%d",i);
d3660d18 211 histTitle = Form("fh2MSubPtRawAll_%d;%s;#it{p}_{T}",i,varName.Data());
f4b02da3 212 fh2MSubPtRawAll[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt);
213 fOutput->Add(fh2MSubPtRawAll[i]);
214
b09cf300 215 histName = Form("fh3MSubPtRawDRMatch_%d",i);
d3660d18 216 histTitle = Form("fh3MSubPtRawDRMatch_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300 217 fh3MSubPtRawDRMatch[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
218 fOutput->Add(fh3MSubPtRawDRMatch[i]);
f4b02da3 219
b09cf300 220 histName = Form("fh3MSubPtTrueDR_%d",i);
d3660d18 221 histTitle = Form("fh3MSubPtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300 222 fh3MSubPtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
223 fOutput->Add(fh3MSubPtTrueDR[i]);
f4b02da3 224
b09cf300 225 histName = Form("fh3MTruePtTrueDR_%d",i);
d3660d18 226 histTitle = Form("fh3MTruePtTrueDR_%d;%s;#it{p}_{T}",i,varName.Data());
b09cf300 227 fh3MTruePtTrueDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsM,minM,maxM,nBinsPt,minPt,maxPt,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
228 fOutput->Add(fh3MTruePtTrueDR[i]);
f4b02da3 229
b09cf300 230 histName = Form("fh3PtTrueDeltaMDR_%d",i);
d3660d18 231 histTitle = Form("fh3PtTrueDeltaMDR_%d;#it{p}_{T,true};#Delta %s",i,varName.Data());
232 fh3PtTrueDeltaMDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDM,minDM,maxDM,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
b09cf300 233 fOutput->Add(fh3PtTrueDeltaMDR[i]);
f4b02da3 234
b09cf300 235 histName = Form("fh3PtTrueDeltaMRelDR_%d",i);
d3660d18 236 histTitle = Form("fh3PtTrueDeltaMRelDR_%d;#it{p}_{T,true};Rel #Delta %s",i,varName.Data());
5ca97880 237 fh3PtTrueDeltaMRelDR[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,400,-1.,3.,nBinsDRToLJ,minDRToLJ,maxDRToLJ);
b09cf300 238 fOutput->Add(fh3PtTrueDeltaMRelDR[i]);
0947b103 239
f4b02da3 240 histName = Form("fhnMassResponse_%d",i);
d3660d18 241 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 242 fhnMassResponse[i] = new THnSparseF(histName.Data(),histTitle.Data(),nBinsSparse0,nBins0,xmin0,xmax0);
243 fOutput->Add(fhnMassResponse[i]);
244 }
245
246 // =========== Switch on Sumw2 for all histos ===========
247 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
248 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
249 if (h1){
250 h1->Sumw2();
251 continue;
252 }
253 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
254 if(hn)hn->Sumw2();
255 }
256
257 TH1::AddDirectory(oldStatus);
258
259 // Create a tree.
260 if(fCreateTree) {
261 fTreeJetBkg = new TTree("fTreeJetSubConst", "fTreeJetSubConst");
262 fTreeJetBkg->Branch("fJet1Vec","TLorentzVector",&fJet1Vec);
263 fTreeJetBkg->Branch("fJet2Vec","TLorentzVector",&fJet2Vec);
264 fTreeJetBkg->Branch("fJetSubVec","TLorentzVector",&fJetSubVec);
265 fTreeJetBkg->Branch("fArea",&fArea,"fArea/F");
266 fTreeJetBkg->Branch("fAreaPhi",&fAreaPhi,"fAreaPhi/F");
267 fTreeJetBkg->Branch("fAreaEta",&fAreaEta,"fAreaEta/F");
268 fTreeJetBkg->Branch("fRho",&fRho,"fRho/F");
269 fTreeJetBkg->Branch("fRhoM",&fRhoM,"fRhoM/F");
270 fTreeJetBkg->Branch("fMatch",&fMatch,"fMatch/I");
271 }
272
273 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
274 if(fCreateTree) PostData(2, fTreeJetBkg);
275}
276
277//________________________________________________________________________
278Bool_t AliAnalysisTaskJetShapeConst::Run()
279{
280 // Run analysis code here, if needed. It will be executed before FillHistograms().
281
282 return kTRUE;
283}
284
285//________________________________________________________________________
286Bool_t AliAnalysisTaskJetShapeConst::FillHistograms()
287{
288 // Fill histograms.
289
8d95aff1 290 AliEmcalJet* jet1 = NULL; //AA jet
291 AliEmcalJet *jet2 = NULL; //Embedded Pythia jet
a5a1c51a 292 // AliEmcalJet *jet1T = NULL; //tagged AA jet
8d95aff1 293 // AliEmcalJet *jet2T = NULL; //tagged Pythia jet
294 AliEmcalJet *jetS = NULL; //subtracted jet
f4b02da3 295 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
296 AliJetContainer *jetContS = GetJetContainer(fContainerSub);
a5a1c51a 297
298 //Get leading jet in Pb-Pb event without embedded objects
299 AliJetContainer *jetContNoEmb = GetJetContainer(fContainerNoEmb);
300 AliEmcalJet *jetL = NULL;
301 if(jetContNoEmb) jetL = jetContNoEmb->GetLeadingJet("rho");
302
f4b02da3 303 AliDebug(11,Form("NJets Incl: %d Csub: %d",jetCont->GetNJets(),jetContS->GetNJets()));
304 if(jetCont) {
305 jetCont->ResetCurrentID();
306 while((jet1 = jetCont->GetNextAcceptJet())) {
8d95aff1 307 jet2 = NULL;
8d95aff1 308 jetS = NULL;
d3660d18 309
f4b02da3 310 //Get constituent subtacted version of jet
311 Int_t ifound = 0;
312 Int_t ilab = -1;
313 for(Int_t i = 0; i<jetContS->GetNJets(); i++) {
314 //if(ifound==1) continue;
315 jetS = jetContS->GetJet(i);
316 if(jetS->GetLabel()==jet1->GetLabel()) { // && jetS->Pt()>0.) {
317 ifound++;
318 if(ifound==1) ilab = i;
319 }
320 }
321 if(ifound>1) AliDebug(2,Form("Found %d partners",ifound));
322 if(ifound==0) jetS = 0x0;
323 else jetS = jetContS->GetJet(ilab);
324
d3660d18 325 Double_t mjetS = jetS->M();
326 Double_t ptjet1 = jet1->Pt()-jetCont->GetRhoVal()*jet1->Area();
327 Double_t var = mjetS;
328 if(fJetMassVarType==kRatMPt) {
329 if(ptjet1>0. || ptjet1<0.) var = mjetS/ptjet1;
330 else var = -999.;
331 }
332
f4b02da3 333 //Fill histograms for all AA jets
d3660d18 334 fh2MSubPtRawAll[fCentBin]->Fill(var,ptjet1);
f4b02da3 335
336 Double_t fraction = 0.;
337 fMatch = 0;
338 fJet2Vec->SetPtEtaPhiM(0.,0.,0.,0.);
339 if(fSingleTrackEmb) {
340 AliVParticle *vp = GetEmbeddedConstituent(jet1);
341 if(vp) {
342 fJet2Vec->SetPxPyPzE(vp->Px(),vp->Py(),vp->Pz(),vp->E());
343 fMatch = 1;
344 }
345 } else {
346 jet2 = jet1->ClosestJet();
347 fraction = jetCont->GetFractionSharedPt(jet1);
348 fMatch = 1;
349 if(fMinFractionShared>0.) {
350 if(fraction>fMinFractionShared) {
351 fJet2Vec->SetPxPyPzE(jet2->Px(),jet2->Py(),jet2->Pz(),jet2->E());
352 fMatch = 1;
353 } else
354 fMatch = 0;
355 }
356 }
357
8d95aff1 358 // if(fMatch==1 && jet2->GetTagStatus()>0) jet2T = jet2->GetTaggedJet();
359
f4b02da3 360 //Fill histograms for matched jets
d3660d18 361 fh2MSubMatch[fCentBin]->Fill(var,fMatch);
f4b02da3 362 if(fMatch==1) {
b09cf300 363 Double_t drToLJ = -1.;
364 if(jetL) drToLJ = jet1->DeltaR(jetL);
d3660d18 365 fh3MSubPtRawDRMatch[fCentBin]->Fill(var,ptjet1,drToLJ);
f4b02da3 366 if(jet2) {
d3660d18 367 Double_t var2 = -999.;
368 if(jet2->Pt()>0. || jet2->Pt()<0.) var2 = jet2->M()/jet2->Pt();
369 fh3MSubPtTrueDR[fCentBin]->Fill(var,jet2->Pt(),drToLJ);
370 fh3MTruePtTrueDR[fCentBin]->Fill(var2,jet2->Pt(),drToLJ);
371 fh3PtTrueDeltaMDR[fCentBin]->Fill(jet2->Pt(),var-var2,drToLJ);
372 if(jet2->M()>0.) fh3PtTrueDeltaMRelDR[fCentBin]->Fill(jet2->Pt(),(var-var2)/var2,drToLJ);
373 Double_t varsp[5] = {var,var2,ptjet1,jet2->Pt(),drToLJ};
374 fhnMassResponse[fCentBin]->Fill(varsp);
f4b02da3 375 }
376 }
377
378 if(fCreateTree) {
379 fJet1Vec->SetPxPyPzE(jet1->Px(),jet1->Py(),jet1->Pz(),jet1->E());
380 if(jetS && jetS->Pt()>0.) fJetSubVec->SetPtEtaPhiM(jetS->Pt(),jetS->Eta(),jetS->Phi(),jetS->M());
381 else fJetSubVec->SetPtEtaPhiM(0.,0.,0.,0.);
382 fArea = (Float_t)jet1->Area();
383 fAreaPhi = (Float_t)jet1->AreaPhi();
384 fAreaEta = (Float_t)jet1->AreaEta();
385 fRho = (Float_t)jetCont->GetRhoVal();
386 fRhoM = (Float_t)jetCont->GetRhoMassVal();
387 fNConst = (Int_t)jet1->GetNumberOfTracks();
388 fTreeJetBkg->Fill();
389 }
390 } //jet1 loop
391 }//jetCont
392
393
394 return kTRUE;
395}
396
397//________________________________________________________________________
398AliVParticle* AliAnalysisTaskJetShapeConst::GetEmbeddedConstituent(AliEmcalJet *jet) {
399
400 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
401 AliVParticle *vp = 0x0;
402 AliVParticle *vpe = 0x0; //embedded particle
403 Int_t nc = 0;
404 for(Int_t i=0; i<jet->GetNumberOfTracks(); i++) {
405 vp = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray())); //check if fTracks is the correct track branch
406 if (vp->TestBits(TObject::kBitMask) != (Int_t)(TObject::kBitMask) ) continue;
407 if(!vpe) vpe = vp;
408 else if(vp->Pt()>vpe->Pt()) vpe = vp;
409 nc++;
410 }
411 AliDebug(11,Form("Found %d embedded particles",nc));
412 return vpe;
413}
414
415
416//________________________________________________________________________
417Bool_t AliAnalysisTaskJetShapeConst::RetrieveEventObjects() {
418 //
419 // retrieve event objects
420 //
421
422 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
423 return kFALSE;
424
425 AliJetContainer *jetCont = GetJetContainer(fContainerBase);
426 jetCont->LoadRhoMass(InputEvent());
427
428 return kTRUE;
429}
430
431//_______________________________________________________________________
432void AliAnalysisTaskJetShapeConst::Terminate(Option_t *)
433{
434 // Called once at the end of the analysis.
435}
436