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