]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetAna.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalDiJetAna.cxx
CommitLineData
6e8b6371 1//
2// Dijet analysis task.
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 <TList.h>
12#include <TLorentzVector.h>
13#include <TProfile.h>
14#include <TChain.h>
15#include <TSystem.h>
16#include <TFile.h>
17#include <TKey.h>
18
19#include "AliVCluster.h"
20#include "AliVTrack.h"
21#include "AliEmcalJet.h"
22#include "AliRhoParameter.h"
23#include "AliLog.h"
6e8b6371 24#include "AliEmcalParticle.h"
25#include "AliMCEvent.h"
26#include "AliGenPythiaEventHeader.h"
27#include "AliAODMCHeader.h"
28#include "AliMCEvent.h"
29#include "AliAnalysisManager.h"
30#include "AliJetContainer.h"
eca9052c 31#include "AliCentrality.h"
6e8b6371 32
33#include "AliAnalysisTaskEmcalDiJetAna.h"
34
35ClassImp(AliAnalysisTaskEmcalDiJetAna)
36
37//________________________________________________________________________
38AliAnalysisTaskEmcalDiJetAna::AliAnalysisTaskEmcalDiJetAna() :
39 AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetAna"),
40 fDoMatchFullCharged(kTRUE),
677e90b9 41 fNKtBins(30),
a29bc208
CL		 42 fNDiJetEtaBins(1),
43 fNAjBins(1),
6e8b6371 44 fh2CentRhoCh(0),
45 fh2CentRhoScaled(0),
46 fh3PtEtaPhiJetFull(0),
47 fh3PtEtaPhiJetCharged(0),
48 fhnDiJetVarsFull(0),
49 fhnDiJetVarsCh(0),
50 fhnDiJetVarsFullCharged(0),
51 fhnMatchingFullCharged(0),
8e49a788 52 fh3PtTrigKt1Kt2Ch(0),
53 fh3PtTrigKt1Kt2FuCh(0),
54 fh3PtTrigDPhi1DPhi2Ch(0),
55 fh3PtTrigDPhi1DPhi2FuCh(0)
6e8b6371 56{
57 // Default constructor.
bb84b374 58
eca9052c 59 for(Int_t i=0; i<4; i++) {
8e49a788 60 fh3DiJetKtNEFPtAssoc[i] = 0;
61 fCentCorrPtAssocCh[i] = 0;
62 fCentCorrPtAssocFuCh[i] = 0;
63 fAjPtAssocCentCh[i] = 0;
64 fAjPtAssocCentFuCh[i] = 0;
65 fh3PtAssoc1PtAssoc2DPhi23Ch[i] = 0;
66 fh3PtAssoc1PtAssoc2DPhi23FuCh[i] = 0;
eca9052c 67 }
68
6e8b6371 69 SetMakeGeneralHistograms(kTRUE);
70}
71
72//________________________________________________________________________
73AliAnalysisTaskEmcalDiJetAna::AliAnalysisTaskEmcalDiJetAna(const char *name) :
74 AliAnalysisTaskEmcalDiJetBase(name),
75 fDoMatchFullCharged(kTRUE),
677e90b9 76 fNKtBins(30),
a29bc208
CL		 77 fNDiJetEtaBins(1),
78 fNAjBins(1),
6e8b6371 79 fh2CentRhoCh(0),
80 fh2CentRhoScaled(0),
81 fh3PtEtaPhiJetFull(0),
82 fh3PtEtaPhiJetCharged(0),
83 fhnDiJetVarsFull(0),
84 fhnDiJetVarsCh(0),
85 fhnDiJetVarsFullCharged(0),
86 fhnMatchingFullCharged(0),
8e49a788 87 fh3PtTrigKt1Kt2Ch(0),
88 fh3PtTrigKt1Kt2FuCh(0),
89 fh3PtTrigDPhi1DPhi2Ch(0),
90 fh3PtTrigDPhi1DPhi2FuCh(0)
6e8b6371 91{
92 // Standard constructor.
93
eca9052c 94 for(Int_t i=0; i<4; i++) {
8e49a788 95 fh3DiJetKtNEFPtAssoc[i] = 0;
96 fCentCorrPtAssocCh[i] = 0;
97 fCentCorrPtAssocFuCh[i] = 0;
98 fAjPtAssocCentCh[i] = 0;
99 fAjPtAssocCentFuCh[i] = 0;
100 fh3PtAssoc1PtAssoc2DPhi23Ch[i] = 0;
101 fh3PtAssoc1PtAssoc2DPhi23FuCh[i] = 0;
eca9052c 102 }
bb84b374 103
6e8b6371 104 SetMakeGeneralHistograms(kTRUE);
105}
106
107//________________________________________________________________________
108AliAnalysisTaskEmcalDiJetAna::~AliAnalysisTaskEmcalDiJetAna()
109{
110 // Destructor.
111}
112
113//________________________________________________________________________
114Bool_t AliAnalysisTaskEmcalDiJetAna::RetrieveEventObjects() {
115 //
116 // retrieve event objects
117 //
118
119 if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
120 return kFALSE;
121
122 return kTRUE;
123}
124
125//________________________________________________________________________
126void AliAnalysisTaskEmcalDiJetAna::UserCreateOutputObjects()
127{
128 // Create user output.
129
6e8b6371 130 AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
131
132 Bool_t oldStatus = TH1::AddDirectoryStatus();
133 TH1::AddDirectory(kFALSE);
134
0f4704b4 135 const Int_t nBinsCent = 100;
6e8b6371 136 Double_t minCent = 0.;
137 Double_t maxCent = 100.;
138
139 const Int_t nBinsRho = 200;
140 Double_t minRho = 0.;
141 Double_t maxRho = 20.;
142 fh2CentRhoCh = new TH2F("fh2CentRhoCh","fh2CentRhoCh;centrality;#rho_{ch}",nBinsCent,minCent,maxCent,nBinsRho,minRho,maxRho);
143 fOutput->Add(fh2CentRhoCh);
144 fh2CentRhoScaled = new TH2F("fh2CentRhoScaled","fh2CentRhoScaled;centrality;s_{EMC}#rho_{ch}",nBinsCent,minCent,maxCent,nBinsRho,minRho,maxRho);
145 fOutput->Add(fh2CentRhoScaled);
146
a822e422 147 const Int_t nBinsPt = 150;
6e8b6371 148 Double_t minPt = -20.;
a822e422 149 Double_t maxPt = 130.;
6e8b6371 150 const Int_t nBinsEta = 40;
151 Double_t minEta = -1.;
152 Double_t maxEta = 1.;
153 const Int_t nBinsPhi = 18*6;
154 Double_t minPhi = 0.;
155 Double_t maxPhi = TMath::TwoPi();
156
157 fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
158 fOutput->Add(fh3PtEtaPhiJetFull);
159
160 fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
161 fOutput->Add(fh3PtEtaPhiJetCharged);
162
e50460de 163 const Int_t nBinsSparse0 = 7;
eca9052c 164 const Int_t nBinsPtW = 30;
6e8b6371 165 const Int_t nBinsDPhi = 72;
677e90b9 166 const Int_t nBinsKt = fNKtBins;
a29bc208 167 const Int_t nBinsDiJetEta = fNDiJetEtaBins;//40;
8daeee93 168 const Int_t nBinsCentr = fNcentBins;
a29bc208 169 const Int_t nBinsAj = fNAjBins;//20
e50460de 170 const Int_t nBins0[nBinsSparse0] = {nBinsPtW,nBinsPtW,nBinsDPhi,nBinsKt,nBinsDiJetEta,nBinsCentr,nBinsAj};
6e8b6371 171 //pT1, pT2, deltaPhi, kT
e50460de 172 const Double_t xmin0[nBinsSparse0] = { minPt, minPt, -0.5*TMath::Pi(), 0.,-1.,0. , 0.};
6c90d6c8 173 const Double_t xmax0[nBinsSparse0] = { maxPt, maxPt, 1.5*TMath::Pi(), 120., 1.,100., 1.};
8daeee93 174 const Double_t centArrayBins[8] = {0.,2.,5.,10.,20.,40.,60.,100.};
6e8b6371 175
176 if(fDoChargedCharged) {
177 fhnDiJetVarsCh = new THnSparseF("fhnDiJetVarsCh",
e50460de 178 "fhnDiJetVarsCh;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality;#it{A}_{j}",
6e8b6371 179 nBinsSparse0,nBins0,xmin0,xmax0);
8daeee93 180 if(fNcentBins==7) fhnDiJetVarsCh->SetBinEdges(5,centArrayBins);
6e8b6371 181 fOutput->Add(fhnDiJetVarsCh);
182 }
183
184 if(fDoFullCharged) {
185 fhnDiJetVarsFullCharged = new THnSparseF("fhnDiJetVarsFullCharged",
e50460de 186 "fhnDiJetVarsFullCharged;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality;#it{A}_{j}",
6e8b6371 187 nBinsSparse0,nBins0,xmin0,xmax0);
8daeee93 188 if(fNcentBins==7) fhnDiJetVarsFullCharged->SetBinEdges(5,centArrayBins);
6e8b6371 189 fOutput->Add(fhnDiJetVarsFullCharged);
190 }
191
bb84b374 192 if(fDoFullFull) {
193 fhnDiJetVarsFull = new THnSparseF("fhnDiJetVarsFull",
e50460de 194 "fhnDiJetVarsFull;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality;#it{A}_{j}",
bb84b374 195 nBinsSparse0,nBins0,xmin0,xmax0);
196 fOutput->Add(fhnDiJetVarsFull);
197 }
198
6c90d6c8 199 fh3PtTrigKt1Kt2Ch = new TH3F("fh3PtTrigKt1Kt2Ch","fh3PtTrigKt1Kt2Ch;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,nBinsKt,0.,120.,nBinsKt,0.,120.);
8e49a788 200 fOutput->Add(fh3PtTrigKt1Kt2Ch);
201
6c90d6c8 202 fh3PtTrigKt1Kt2FuCh = new TH3F("fh3PtTrigKt1Kt2FuCh","fh3PtTrigKt1Kt2FuCh;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,nBinsKt,0.,120.,nBinsKt,0.,120.);
8e49a788 203 fOutput->Add(fh3PtTrigKt1Kt2FuCh);
204
6c90d6c8 205 fh3PtTrigDPhi1DPhi2Ch = new TH3F("fh3PtTrigDPhi1DPhi2Ch","fh3PtTrigDPhi1DPhi2Ch;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,36,-0.5*TMath::Pi(),1.5*TMath::Pi(),36,-0.5*TMath::Pi(),1.5*TMath::Pi());
8e49a788 206 fOutput->Add(fh3PtTrigDPhi1DPhi2Ch);
207
6c90d6c8 208 fh3PtTrigDPhi1DPhi2FuCh = new TH3F("fh3PtTrigDPhi1DPhi2FuCh","fh3PtTrigDPhi1DPhi2FuCh;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,36,-0.5*TMath::Pi(),1.5*TMath::Pi(),36,-0.5*TMath::Pi(),1.5*TMath::Pi());
8e49a788 209 fOutput->Add(fh3PtTrigDPhi1DPhi2FuCh);
210
211
bb84b374 212 for(Int_t i=0; i<4; i++) {
213 TString histoName = Form("fh3DiJetKtNEFPtAssoc_TrigBin%d",i);
6c90d6c8 214 fh3DiJetKtNEFPtAssoc[i] = new TH3F(histoName.Data(),histoName.Data(),nBinsKt,0.,120.,50,0.,1.,nBinsPt,minPt,maxPt);
bb84b374 215 fOutput->Add(fh3DiJetKtNEFPtAssoc[i]);
eca9052c 216
217 histoName = Form("fCentCorrPtAssocCh_TrigBin%d",i);
6c90d6c8 218 fCentCorrPtAssocCh[i] = new TH3F(histoName.Data(),histoName.Data(),10,0.,100.,10,0.,100.,nBinsPt,minPt,maxPt);
eca9052c 219 fOutput->Add(fCentCorrPtAssocCh[i]);
220
221 histoName = Form("fCentCorrPtAssocFuCh_TrigBin%d",i);
6c90d6c8 222 fCentCorrPtAssocFuCh[i] = new TH3F(histoName.Data(),histoName.Data(),10,0.,100.,10,0.,100.,nBinsPt,minPt,maxPt);
eca9052c 223 fOutput->Add(fCentCorrPtAssocFuCh[i]);
224
0327741d 225 histoName = Form("fAjPtAssocCentCh_TrigBin%d",i);
6c90d6c8 226 fAjPtAssocCentCh[i] = new TH3F(histoName.Data(),histoName.Data(),50,0.,1.,nBinsPt,minPt,maxPt,20,0.,100.);
0327741d 227 fOutput->Add(fAjPtAssocCentCh[i]);
228
229 histoName = Form("fAjPtAssocCentFuCh_TrigBin%d",i);
6c90d6c8 230 fAjPtAssocCentFuCh[i] = new TH3F(histoName.Data(),histoName.Data(),50,0.,1.,nBinsPt,minPt,maxPt,20,0.,100.);
0327741d 231 fOutput->Add(fAjPtAssocCentFuCh[i]);
232
8e49a788 233 histoName = Form("fh3PtAssoc1PtAssoc2DPhi23Ch_TrigBin%d",i);
234 fh3PtAssoc1PtAssoc2DPhi23Ch[i] = new TH3F(histoName.Data(),histoName.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
235 fOutput->Add(fh3PtAssoc1PtAssoc2DPhi23Ch[i]);
0327741d 236
8e49a788 237 histoName = Form("fh3PtAssoc1PtAssoc2DPhi23FuCh_TrigBin%d",i);
238 fh3PtAssoc1PtAssoc2DPhi23FuCh[i] = new TH3F(histoName.Data(),histoName.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
239 fOutput->Add(fh3PtAssoc1PtAssoc2DPhi23FuCh[i]);
bb84b374 240 }
241
6e8b6371 242 const Int_t nBinsSparseMatch = 7;
243 const Int_t nBinsDPhiMatch = 80;
244 const Int_t nBinsDEtaMatch = 80;
245 const Int_t nBinsDR = 20;
246 const Int_t nBinsFraction = 21;
247 const Int_t nBinsType = 3;
248 const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsFraction,nBinsType};
249 //pTfull, pTch, deltaPhi, deltaEta, deltaR, fraction, jet type (leading,subleading,other)
6c90d6c8 250 const Double_t xminMatch[nBinsSparseMatch] = { minPt, minPt, -0.5,-0.5, 0. ,0. ,0};
6e8b6371 251 const Double_t xmaxMatch[nBinsSparseMatch] = { maxPt, maxPt, 0.5, 0.5, 0.5,1.05,3};
8daeee93 252 if(fDoMatchFullCharged) {
253 fhnMatchingFullCharged = new THnSparseF("fhnMatchingFullCharged","fhnMatchingFullCharged;#it{p}_{T,full} (GeV/#it{c});#it{p}_{T,ch} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;f_{ch};type",
6e8b6371 254 nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
8daeee93 255 fOutput->Add(fhnMatchingFullCharged);
256 }
6e8b6371 257
258 // =========== Switch on Sumw2 for all histos ===========
259 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
260 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
261 if (h1){
262 h1->Sumw2();
263 continue;
264 }
265 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
266 if(hn) {
267 hn->Sumw2();
268 continue;
269 }
270 }
271
272 TH1::AddDirectory(oldStatus);
273
274 PostData(1, fOutput); // Post data for ALL output slots > 0 here.
275}
276
277//________________________________________________________________________
278Bool_t AliAnalysisTaskEmcalDiJetAna::FillHistograms()
279{
280 // Fill histograms.
281
282 fh2CentRhoCh->Fill(fCent,fRhoChVal);
283 fh2CentRhoScaled->Fill(fCent,fRhoFullVal);
284
6e8b6371 285 Int_t nJetsFull = GetNJets(fContainerFull);
286 for(Int_t ij=0; ij<nJetsFull; ij++) {
287 AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, fContainerFull));
288 if(!jet) continue; //jet not selected
289
290 Double_t jetPt = GetJetPt(jet,0);
291 fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
6e8b6371 292 }
293
294 Int_t nJetsCh = GetNJets(fContainerCharged);
295 for(Int_t ij=0; ij<nJetsCh; ij++) {
296 AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, fContainerCharged));
297 if(!jet) continue; //jet not selected
298
299 Double_t jetPt = GetJetPt(jet,1);
300 fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
6e8b6371 301 }
6e8b6371 302
303 return kTRUE;
304}
305
306//________________________________________________________________________
307Bool_t AliAnalysisTaskEmcalDiJetAna::Run()
308{
309 // Run analysis code here, if needed. It will be executed before FillHistograms().
310
311 //Check if event is selected (vertex & pile-up)
312 if(!SelectEvent())
313 return kFALSE;
314
6e8b6371 315 if(fRhoType==0) {
316 fRhoFullVal = 0.;
317 fRhoChVal = 0.;
318 }
319 if(fRhoType==1) {
320 fRhoFullVal = GetRhoVal(fContainerFull);
321 fRhoChVal = GetRhoVal(fContainerCharged);
322 }
323
4bd08270 324 if(fDoFullFull)
325 CorrelateJets(0);
bb84b374 326
6e8b6371 327 // MatchFullAndChargedJets();
328 if(fDoChargedCharged) CorrelateJets(1);
329
330 if(fDoMatchFullCharged) FillMatchFullChargedHistos(fContainerFull,fContainerCharged);
331
332 if(fDoFullCharged) {
333 SetChargedFractionIndex();
334 CorrelateJets(2);
335 }
336
337 return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
338}
339
c66364a6 340//________________________________________________________________________
341void AliAnalysisTaskEmcalDiJetAna::CorrelateJets(const Int_t type) {
342 //
343 // Correlate jets and fill histos
344 //
345
346 if( fJetCorrelationType==kCorrelateAll )
347 CorrelateAllJets(type);
348 else if( fJetCorrelationType==kCorrelateTwo )
349 CorrelateTwoJets(type);
350 else if( fJetCorrelationType==kCorrelateLS )
351 CorrelateLeadingSubleadingJets(type);
352
353 return;
354
355}
356
4bd08270 357//________________________________________________________________________
358void AliAnalysisTaskEmcalDiJetAna::CorrelateTwoJets(const Int_t type) {
359 //
360 // Correlate jets and fill histos
361 //
362
363 Int_t typet = 0;
364 Int_t typea = 0;
365 if(type==0) { //full-full
366 typet = fContainerFull;
367 typea = fContainerFull;
368 }
369 else if(type==1) { //charged-charged
370 typet = fContainerCharged;
371 typea = fContainerCharged;
372 }
373 else if(type==2) { //full-charged
374 typet = fContainerFull;
375 typea = fContainerCharged;
376 }
377 else {
378 AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
379 return;
380 }
381
6c90d6c8 382 Int_t nJetsTrig = GetNJets(typet);
4bd08270 383 for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
384
385 AliEmcalJet *jetTrig = NULL;
386 if(type==0) {
387 jetTrig = static_cast<AliEmcalJet*>(GetJetFromArray(ijt, typet));
388 if(TMath::Abs(jetTrig->Eta())>0.5)
389 jetTrig = NULL;
390 }
391 else
392 jetTrig = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typet));
393
4bd08270 394 if(!jetTrig)
395 continue; //jet not selected
396
397 Double_t jetTrigPt = GetJetPt(jetTrig,typet);
398
399 if(jetTrigPt<fPtMinTriggerJet)
400 continue;
401
c66364a6 402 AliEmcalJet *jetAssoc = GetLeadingAssociatedJet(typea,jetTrig);
4bd08270 403 if(!jetAssoc)
404 continue;
405
f9de9b73 406 if(fDoPtBias) {
407 if(type==0 || type==1) {
408 if(GetJetPt(jetAssoc,typea)>jetTrigPt)
409 continue;
410 }
3c709670 411 }
412
4bd08270 413 FillDiJetHistos(jetTrig,jetAssoc, type);
8e49a788 414
415 //Look for second jet on away side - 3-jet events
416 AliEmcalJet *jetAssoc2 = GetSecondLeadingAssociatedJet(typea,jetTrig);
417 if(jetAssoc2)
418 FillThreeJetHistos(jetTrig,jetAssoc,jetAssoc2,type);
4bd08270 419
4bd08270 420 }
4bd08270 421}
422
423//________________________________________________________________________
c66364a6 424AliEmcalJet* AliAnalysisTaskEmcalDiJetAna::GetLeadingJet(const Int_t type) {
425
426 //Get associated jet which is the leading jet in the opposite hemisphere
427
428 Int_t cont = 0;
429 if(type==0) //full-full
430 cont = fContainerFull;
431 else if(type==1) //charged-charged
432 cont = fContainerCharged;
433 else if(type==2) //full-charged
434 cont = fContainerFull;
435
436 Int_t nJets = GetNJets(cont);
437 Double_t ptLead = -999;
438 Int_t iJetLead = -1;
439 for(Int_t ij=0; ij<nJets; ij++) {
c66364a6 440 AliEmcalJet *jet = NULL;
441 if(type==0) {
442 jet = static_cast<AliEmcalJet*>(GetJetFromArray(ij, cont));
443 if(TMath::Abs(jet->Eta())>0.5)
444 jet = NULL;
445 }
446 else
447 jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, cont));
448
449 if(!jet)
450 continue;
451
452 Double_t jetPt = GetJetPt(jet,cont);
453
454 if(jetPt>ptLead) {
455 ptLead = jetPt;
456 iJetLead = ij;
457 }
c66364a6 458 }
459
460 AliEmcalJet *jetLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, cont));
461
462 return jetLead;
c66364a6 463}
464
465//________________________________________________________________________
466AliEmcalJet* AliAnalysisTaskEmcalDiJetAna::GetLeadingAssociatedJet(const Int_t type, AliEmcalJet *jetTrig) {
4bd08270 467
468 //Get associated jet which is the leading jet in the opposite hemisphere
469
470 Int_t typea = 0;
471 if(type==0) //full-full
472 typea = fContainerFull;
473 else if(type==1) //charged-charged
474 typea = fContainerCharged;
475 else if(type==2) //full-charged
476 typea = fContainerCharged;
4bd08270 477
6ab30d5f 478 AliEmcalJet *jetAssocLead = GetLeadingJetOppositeHemisphere(type, typea, jetTrig);
4bd08270 479
4bd08270 480 return jetAssocLead;
4bd08270 481}
482
8e49a788 483//________________________________________________________________________
484AliEmcalJet* AliAnalysisTaskEmcalDiJetAna::GetSecondLeadingAssociatedJet(const Int_t type, AliEmcalJet *jetTrig) {
485
486 //Get associated jet which is the leading jet in the opposite hemisphere
487
488 Int_t typea = 0;
489 if(type==0) //full-full
490 typea = fContainerFull;
491 else if(type==1) //charged-charged
492 typea = fContainerCharged;
493 else if(type==2) //full-charged
494 typea = fContainerCharged;
495
496 AliEmcalJet *jetAssocLead2 = GetSecondLeadingJetOppositeHemisphere(type, typea, jetTrig);
497
498 return jetAssocLead2;
8e49a788 499}
500
6e8b6371 501//________________________________________________________________________
c66364a6 502void AliAnalysisTaskEmcalDiJetAna::CorrelateAllJets(const Int_t type) {
6e8b6371 503 //
504 // Correlate jets and fill histos
505 //
506
507 Int_t typet = 0;
508 Int_t typea = 0;
509 if(type==0) { //full-full
510 typet = fContainerFull;
511 typea = fContainerFull;
512 }
513 else if(type==1) { //charged-charged
514 typet = fContainerCharged;
515 typea = fContainerCharged;
516 }
517 else if(type==2) { //full-charged
518 typet = fContainerFull;
519 typea = fContainerCharged;
520 }
521 else {
522 AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
523 return;
524 }
525
6c90d6c8 526 Int_t nJetsTrig = GetNJets(typet);
527 Int_t nJetsAssoc = GetNJets(typea);
6e8b6371 528
529 for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
bb84b374 530 AliEmcalJet *jetTrig = NULL;
531 if(type==0) {
532 jetTrig = static_cast<AliEmcalJet*>(GetJetFromArray(ijt, typet));
533 if(TMath::Abs(jetTrig->Eta())>0.5)
534 jetTrig = NULL;
535 }
536 else
537 jetTrig = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typet));
538
6e8b6371 539 if(!jetTrig)
540 continue; //jet not selected
541
542 Double_t jetTrigPt = GetJetPt(jetTrig,typet);
543
544 if(jetTrigPt<fPtMinTriggerJet)
545 continue;
546
547 for(Int_t ija=0; ija<nJetsAssoc; ija++) {
548 if(IsSameJet(ijt,ija,type)) continue;
549
bb84b374 550 AliEmcalJet *jetAssoc = NULL;
551 if(type==0) {
552 jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
553 if(TMath::Abs(jetAssoc->Eta())>0.5)
554 jetAssoc = NULL;
555 }
556 else
557 jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
558
6e8b6371 559 if(!jetAssoc)
560 continue;
561
562 Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
563
564 if(jetTrigPt>jetAssocPt)
565 FillDiJetHistos(jetTrig,jetAssoc, type);
566
567 } // associate jet loop
568 }//trigger jet loop
569
570}
571
572//________________________________________________________________________
c66364a6 573void AliAnalysisTaskEmcalDiJetAna::CorrelateLeadingSubleadingJets(const Int_t type) {
574 //
575 // Correlate leading jet in event with leading jet in opposite hemisphere
576 //
577
578 Int_t typet = 0;
579 Int_t typea = 0;
580 if(type==0) { //full-full
581 typet = fContainerFull;
582 typea = fContainerFull;
583 }
584 else if(type==1) { //charged-charged
585 typet = fContainerCharged;
586 typea = fContainerCharged;
587 }
588 else if(type==2) { //full-charged
589 typet = fContainerFull;
590 typea = fContainerCharged;
591 }
592 else {
593 AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
594 return;
595 }
596
6c90d6c8 597 AliEmcalJet *jetTrig = GetLeadingJet(typet);
c66364a6 598 if(!jetTrig)
599 return;
600
601 Double_t jetTrigPt = GetJetPt(jetTrig,typet);
602
603 if(jetTrigPt<fPtMinTriggerJet)
604 return;
605
606 AliEmcalJet *jetAssoc = GetLeadingAssociatedJet(typea,jetTrig);
607 if(!jetAssoc)
608 return;
609
c66364a6 610 FillDiJetHistos(jetTrig,jetAssoc, type);
c66364a6 611}
612
613//________________________________________________________________________
6e8b6371 614void AliAnalysisTaskEmcalDiJetAna::FillDiJetHistos(const AliEmcalJet *jet1, const AliEmcalJet *jet2, const Int_t mode) {
615 //
616 // Fill histos
bb84b374 617 // mode: full vs full = 0
618 // charged vs charged = 1
619 // full vs charged = 2
6e8b6371 620 //
621
622 Int_t typet = 0;
623 Int_t typea = 0;
6c90d6c8 624 if(mode==0) { //full-full
625 typet = fContainerFull;
626 typea = fContainerFull;
6e8b6371 627 }
6c90d6c8 628 else if(mode==1) { //charged-charged
629 typet = fContainerCharged;
630 typea = fContainerCharged;
6e8b6371 631 }
6c90d6c8 632 else if(mode==2) { //full-charged
633 typet = fContainerFull;
634 typea = fContainerCharged;
6e8b6371 635 }
636 else {
637 AliWarning(Form("%s: mode %d of dijet correlation not defined!",GetName(),mode));
638 return;
639 }
640
641 Double_t jetTrigPt = GetJetPt(jet1,typet);
642 Double_t jetAssocPt = GetJetPt(jet2,typea);
643
644 Double_t deltaPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
6e8b6371 645
e50460de 646 Double_t kT = TMath::Abs(jetTrigPt*TMath::Sin(deltaPhi));
6e8b6371 647
648 Double_t dijetEta = (jet1->Eta()+jet2->Eta())/2.;
649
6c90d6c8 650 Double_t aj = 0.;
651 if((jetTrigPt+jetAssocPt)>0.) aj = (jetTrigPt-jetAssocPt)/(jetTrigPt+jetAssocPt);
6e8b6371 652
e50460de 653 Double_t diJetVars[7] = {jetTrigPt,jetAssocPt,deltaPhi,kT,dijetEta,fCent,aj};
6e8b6371 654
655 if(mode==0)
656 fhnDiJetVarsFull->Fill(diJetVars);
657 else if(mode==1)
658 fhnDiJetVarsCh->Fill(diJetVars);
659 else if(mode==2)
660 fhnDiJetVarsFullCharged->Fill(diJetVars);
661
8e49a788 662 Double_t dPhiMin = TMath::Pi() - 1./3.*TMath::Pi();
663 Double_t dPhiMax = TMath::Pi() + 1./3.*TMath::Pi();
eca9052c 664 Int_t trigBin = GetPtTriggerBin(jetTrigPt);
bb84b374 665 if(mode==2) {
bb84b374 666 if(trigBin>-1 && trigBin<4) {
bb84b374 667 if(deltaPhi>dPhiMin && deltaPhi<dPhiMax)
668 fh3DiJetKtNEFPtAssoc[trigBin]->Fill(kT, jet1->NEF(), jetAssocPt);
669 }
670 }
671
eca9052c 672 //Fill centrality correlation histos in case a dijet is present in acceptance
eca9052c 673 Double_t centZNA = -1.;
674 AliCentrality *aliCent = InputEvent()->GetCentrality();
675 if (aliCent) {
676 centZNA = aliCent->GetCentralityPercentile("ZNA");
eca9052c 677 if(trigBin>-1 && trigBin<4) {
0327741d 678 if(deltaPhi>dPhiMin && deltaPhi<dPhiMax) {
679 if(mode==1) {
680 fCentCorrPtAssocCh[trigBin]->Fill(fCent,centZNA,jetAssocPt);
681 fAjPtAssocCentCh[trigBin]->Fill(aj,jetAssocPt,fCent);
682 }
683 else if(mode==2) {
684 fCentCorrPtAssocFuCh[trigBin]->Fill(fCent,centZNA,jetAssocPt);
685 fAjPtAssocCentFuCh[trigBin]->Fill(aj,jetAssocPt,fCent);
686 }
687 }
eca9052c 688 }
689 }
690
691}
8e49a788 692
693//________________________________________________________________________
694void AliAnalysisTaskEmcalDiJetAna::FillThreeJetHistos(const AliEmcalJet *jet1, const AliEmcalJet *jet2, const AliEmcalJet *jet3, const Int_t mode) {
695 //
696 // Fill histos
697 // mode: full vs full = 0
698 // charged vs charged = 1
699 // full vs charged = 2
700 //
701
702 Int_t typet = 0;
703 Int_t typea = 0;
6c90d6c8 704 if(mode==0) { //full-full
705 typet = fContainerFull;
706 typea = fContainerFull;
8e49a788 707 }
6c90d6c8 708 else if(mode==1) { //charged-charged
709 typet = fContainerCharged;
710 typea = fContainerCharged;
8e49a788 711 }
6c90d6c8 712 else if(mode==2) { //full-charged
713 typet = fContainerFull;
714 typea = fContainerCharged;
8e49a788 715 }
716 else {
717 AliWarning(Form("%s: mode %d of dijet correlation not defined!",GetName(),mode));
718 return;
719 }
720
721 Double_t jetTrigPt = GetJetPt(jet1,typet);
722 Double_t jetAssoc2Pt = GetJetPt(jet2,typea);
723 Double_t jetAssoc3Pt = GetJetPt(jet3,typea);
724
725 Double_t deltaPhi12 = GetDeltaPhi(jet1->Phi(),jet2->Phi());
726 Double_t deltaPhi13 = GetDeltaPhi(jet1->Phi(),jet3->Phi());
727 Double_t deltaPhi23 = GetDeltaPhi(jet2->Phi(),jet3->Phi());
728
729 Double_t kT12 = TMath::Abs(jetTrigPt*TMath::Sin(deltaPhi12));
730 Double_t kT13 = TMath::Abs(jetTrigPt*TMath::Sin(deltaPhi13));
731
732 Double_t dPhiMin = TMath::Pi() - 1./3.*TMath::Pi();
733 Double_t dPhiMax = TMath::Pi() + 1./3.*TMath::Pi();
734
735 Int_t trigBin = GetPtTriggerBin(jetTrigPt);
736
737 if(jetAssoc2Pt>20. && jetAssoc3Pt>20.) {
738 if(mode==1) {
739 fh3PtTrigDPhi1DPhi2Ch->Fill(jetTrigPt,deltaPhi12,deltaPhi13);
740 fh3PtAssoc1PtAssoc2DPhi23Ch[trigBin]->Fill(jetAssoc2Pt,jetAssoc3Pt,deltaPhi23);
741 }
742 else if(mode==1) {
743 fh3PtTrigDPhi1DPhi2FuCh->Fill(jetTrigPt,deltaPhi12,deltaPhi13);
744 fh3PtAssoc1PtAssoc2DPhi23FuCh[trigBin]->Fill(jetAssoc2Pt,jetAssoc3Pt,deltaPhi23);
745 }
746
747 if(deltaPhi12>dPhiMin && deltaPhi12<dPhiMax) {
748 if(mode==1)
749 fh3PtTrigKt1Kt2Ch->Fill(jetTrigPt,kT12,kT13);
750 else if(mode==2)
751 fh3PtTrigKt1Kt2FuCh->Fill(jetTrigPt,kT12,kT13);
752 }
8e49a788 753 }
754
755}
756
bb84b374 757//________________________________________________________________________
758Int_t AliAnalysisTaskEmcalDiJetAna::GetPtTriggerBin(Double_t pt) {
759
760 Int_t binTrig = -1;
761 if(pt>=20 && pt<40)
762 binTrig = 0;
763 else if(pt>=40 && pt<60)
764 binTrig = 1;
765 else if(pt>=60 && pt<80)
766 binTrig = 2;
767 else if(pt>=80 && pt<100)
768 binTrig = 3;
769
770 return binTrig;
6e8b6371 771}
772
773//________________________________________________________________________
774void AliAnalysisTaskEmcalDiJetAna::FillMatchFullChargedHistos(Int_t cFull,Int_t cCharged) {
775 //
776 // Match full to charged jets and fill histo
777 //
778
779 Int_t match = MatchFullAndChargedJets(cFull,cCharged);
780 if(match==0) {
0560ee0f 781 AliDebug(11,Form("%s: matching failed",GetName()));
6e8b6371 782 return;
783 }
784
785 for(int ig = 0;ig < GetNJets(cFull);++ig){
786 AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ig, cFull));
787 if(!jetFull) continue;
788
789 AliEmcalJet *jetCh = jetFull->ClosestJet();
790 if(!jetCh) continue;
791
792 Double_t shFraction = GetFractionSharedPt(jetFull,jetCh);
6e8b6371 793 Double_t matchVars[7] = {
794 jetFull->Pt(),
795 jetCh->Pt(),
796 GetDeltaPhi(jetFull->Phi(),jetCh->Phi()),
797 jetFull->Eta()-jetCh->Eta(),
798 GetDeltaR(jetFull,jetCh),
799 shFraction,TMath::Min((Float_t)ig+0.5,2.5)
800 };
801 fhnMatchingFullCharged->Fill(matchVars);
802
803 }//loop over full jets
804
805}
806
807
808//________________________________________________________________________
809Int_t AliAnalysisTaskEmcalDiJetAna::MatchFullAndChargedJets(Int_t cFull, Int_t cCharged) {
810 //
811 // Match charged jets to full jets
812 //
813
814 if(GetNJets(cFull)<1) {
0560ee0f 815 AliDebug(2,Form("%s: no full jets: %d", GetName(),GetNJets(cFull)));
6e8b6371 816 return 0;
817 }
818
819 if(GetNJets(cCharged)<1) {
0560ee0f 820 AliDebug(2,Form("%s: no charged jets: %d", GetName(),GetNJets(cCharged)));
6e8b6371 821 return 0;
822 }
823
824 TClonesArray *cJetsFull = GetJetArray(cFull);
825 TClonesArray *cJetsCharged = GetJetArray(cCharged);
826
827 if(!cJetsFull) {
0560ee0f 828 AliDebug(2,Form("%s: no full jet array",GetName()));
6e8b6371 829 return 0;
830 }
831
832 if(!cJetsCharged) {
0560ee0f 833 AliDebug(2,Form("%s: no charged jet array",GetName()));
6e8b6371 834 return 0;
835 }
836
6e8b6371 837 if(!fMatchingDone) {
0560ee0f 838 MatchJetsGeo(cFull, cCharged, 0);
6e8b6371 839 return 1;
840 } else {
0560ee0f 841 AliDebug(11,Form("%s: Matching already done before",GetName()));
6e8b6371 842 return 1;
843 }
844
845}
846
847//_______________________________________________________________________
848void AliAnalysisTaskEmcalDiJetAna::Terminate(Option_t *)
849{
850 // Called once at the end of the analysis.
851}
852
853
854