]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHMEC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
cmake: Fix installation for AliRoot-config.cmake
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetHMEC.cxx
CommitLineData
7f71d351 1// $Id$
5b38fc6d 2
1af24178 3//////////
4//Measure Jet-hadron correlations
5//Does event Mixing using AliEventPoolManager
6/////////
7
5b38fc6d 8#include "AliAnalysisTaskEmcalJetHMEC.h"
9
1101e468 10#include "TChain.h"
11#include "TTree.h"
12#include "TList.h"
13#include "TH1F.h"
14#include "TH2F.h"
15#include "THnSparse.h"
16#include "TCanvas.h"
17#include <TClonesArray.h>
18#include <TParticle.h>
19#include "AliVTrack.h"
20#include "TParameter.h"
21
22#include "AliAODEvent.h"
23#include "AliAnalysisTask.h"
24#include "AliAnalysisManager.h"
25
26#include "AliESDEvent.h"
27#include "AliESDInputHandler.h"
28#include "AliESDCaloCluster.h"
29#include "AliESDVertex.h"
30#include "AliCentrality.h"
31#include "AliAODJet.h"
32#include "AliEmcalJet.h"
33#include "AliESDtrackCuts.h"
1101e468 34
5b38fc6d 35#include "TVector3.h"
36#include "AliPicoTrack.h"
55f64c2f 37#include "AliEventPoolManager.h"
5b38fc6d 38
1101e468 39ClassImp(AliAnalysisTaskEmcalJetHMEC)
40
41//________________________________________________________________________
42AliAnalysisTaskEmcalJetHMEC::AliAnalysisTaskEmcalJetHMEC() :
f569a5d2 43 AliAnalysisTaskEmcalJet("HMEC",kFALSE),
44 fTracksName(""),
45 fJetsName(""),
1101e468 46 fPhimin(-10),
47 fPhimax(10),
48 fEtamin(-0.9),
49 fEtamax(0.9),
50 fAreacut(0.0),
965c985f 51 fTrkBias(5),
52 fClusBias(5),
22d3f9c2 53 fTrkEta(0.9),
fab37389 54 fDoEventMixing(0),
b1121dd6 55 fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
56 fTriggerEventType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
57 fDoLessSparseAxes(0), fDoWiderTrackBin(0),
58 fCentBinSize(1),
f569a5d2 59 fESD(0),
60 fAOD(0),
55f64c2f 61 fPoolMgr(0x0),
1101e468 62 fHistTrackPt(0),
63 fHistCentrality(0),
64 fHistJetEtaPhi(0),
fab37389 65 fHistJetHEtaPhi(0),
22d3f9c2 66 fhnMixedEvents(0x0),
67 fhnJH(0x0)
1101e468 68{
69 // Default Constructor
70
22d3f9c2 71 for(Int_t ipta=0; ipta<7; ipta++){
72 fHistTrackEtaPhi[ipta]=0;
73 }
74
1101e468 75 for(Int_t icent = 0; icent<6; ++icent){
76 fHistJetPt[icent]=0;
81941f29 77 fHistJetPtBias[icent]=0;
41f87d59 78 fHistLeadJetPt[icent]=0;
79 fHistLeadJetPtBias[icent]=0;
81941f29 80 fHistJetPtTT[icent]=0;
55f64c2f 81 for(Int_t iptjet = 0; iptjet<5; ++iptjet){
1101e468 82 for(Int_t ieta = 0; ieta<3; ++ieta){
f569a5d2 83 fHistJetH[icent][iptjet][ieta]=0;
84 fHistJetHBias[icent][iptjet][ieta]=0;
85 fHistJetHTT[icent][iptjet][ieta]=0;
1101e468 86 }
87 }
88 }
89
90}
91//________________________________________________________________________
92AliAnalysisTaskEmcalJetHMEC::AliAnalysisTaskEmcalJetHMEC(const char *name) :
f569a5d2 93 AliAnalysisTaskEmcalJet(name,kTRUE),
94 fTracksName(""),
95 fJetsName(""),
1101e468 96 fPhimin(-10),
97 fPhimax(10),
98 fEtamin(-0.9),
99 fEtamax(0.9),
100 fAreacut(0.0),
965c985f 101 fTrkBias(5),
102 fClusBias(5),
22d3f9c2 103 fTrkEta(0.9),
55f64c2f 104 fDoEventMixing(0),
b1121dd6 105 fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
106 fTriggerEventType(AliVEvent::kEMCEJE), fMixingEventType(AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral),
107 fDoLessSparseAxes(0), fDoWiderTrackBin(0),
108 fCentBinSize(1),
f569a5d2 109 fESD(0),
110 fAOD(0),
55f64c2f 111 fPoolMgr(0x0),
1101e468 112 fHistTrackPt(0),
113 fHistCentrality(0),
114 fHistJetEtaPhi(0),
55f64c2f 115 fHistJetHEtaPhi(0),
22d3f9c2 116 fhnMixedEvents(0x0),
117 fhnJH(0x0)
1101e468 118{
119 // Constructor
22d3f9c2 120 for(Int_t ipta=0; ipta<7; ipta++){
121 fHistTrackEtaPhi[ipta]=0;
122 }
1101e468 123 for(Int_t icent = 0; icent<6; ++icent){
124 fHistJetPt[icent]=0;
125 fHistJetPtBias[icent]=0;
41f87d59 126 fHistLeadJetPt[icent]=0;
127 fHistLeadJetPtBias[icent]=0;
81941f29 128 fHistJetPtTT[icent]=0;
129 for(Int_t iptjet = 0; iptjet<5; ++iptjet){
1101e468 130 for(Int_t ieta = 0; ieta<3; ++ieta){
f569a5d2 131 fHistJetH[icent][iptjet][ieta]=0;
132 fHistJetHBias[icent][iptjet][ieta]=0;
133 fHistJetHTT[icent][iptjet][ieta]=0;
1101e468 134 }
135 }
136 }
137
1101e468 138}
139
140//________________________________________________________________________
f569a5d2 141void AliAnalysisTaskEmcalJetHMEC::UserCreateOutputObjects() {
142 // Called once
143 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
1101e468 144 OpenFile(1);
1101e468 145
55f64c2f 146 // Create histograms
1101e468 147 fHistTrackPt = new TH1F("fHistTrackPt", "P_{T} distribution", 1000, 0.0, 100.0);
1101e468 148 fHistCentrality = new TH1F("fHistCentrality","centrality",100,0,100);
9e4a0350 149 fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet eta-phi",900,-1.8,1.8,720,-3.2,3.2);
150 fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron deta-dphi",900,-1.8,1.8,720,-1.6,4.8);
1101e468 151
daf7dbf5 152 TString name;
1101e468 153
22d3f9c2 154 for(Int_t ipta=0; ipta<7; ++ipta){
daf7dbf5 155 name = Form("fHistTrackEtaPhi_%i", ipta);
9e4a0350 156 fHistTrackEtaPhi[ipta] = new TH2F(name,name,400,-1,1,720,0.0,2.0*TMath::Pi());
f569a5d2 157 fOutput->Add(fHistTrackEtaPhi[ipta]);
22d3f9c2 158 }
1101e468 159
160 for(Int_t icent = 0; icent<6; ++icent){
daf7dbf5 161 name = Form("fHistJetPt_%i",icent);
1101e468 162 fHistJetPt[icent] = new TH1F(name,name,200,0,200);
f569a5d2 163 fOutput->Add(fHistJetPt[icent]);
1101e468 164
daf7dbf5 165 name = Form("fHistJetPtBias_%i",icent);
1101e468 166 fHistJetPtBias[icent] = new TH1F(name,name,200,0,200);
f569a5d2 167 fOutput->Add(fHistJetPtBias[icent]);
1101e468 168
daf7dbf5 169 name = Form("fHistLeadJetPt_%i",icent);
41f87d59 170 fHistLeadJetPt[icent] = new TH1F(name,name,200,0,200);
f569a5d2 171 fOutput->Add(fHistLeadJetPt[icent]);
41f87d59 172
daf7dbf5 173 name = Form("fHistLeadJetPtBias_%i",icent);
41f87d59 174 fHistLeadJetPtBias[icent] = new TH1F(name,name,200,0,200);
f569a5d2 175 fOutput->Add(fHistLeadJetPtBias[icent]);
41f87d59 176
daf7dbf5 177 name = Form("fHistJetPtTT_%i",icent);
81941f29 178 fHistJetPtTT[icent] = new TH1F(name,name,200,0,200);
f569a5d2 179 fOutput->Add(fHistJetPtTT[icent]);
81941f29 180
181 for(Int_t iptjet = 0; iptjet<5; ++iptjet){
1101e468 182 for(Int_t ieta = 0; ieta<3; ++ieta){
f569a5d2 183 name = Form("fHistJetH_%i_%i_%i",icent,iptjet,ieta);
184 fHistJetH[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
185 fOutput->Add(fHistJetH[icent][iptjet][ieta]);
1101e468 186
f569a5d2 187 name = Form("fHistJetHBias_%i_%i_%i",icent,iptjet,ieta);
188 fHistJetHBias[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
189 fOutput->Add(fHistJetHBias[icent][iptjet][ieta]);
1101e468 190
f569a5d2 191 name = Form("fHistJetHTT_%i_%i_%i",icent,iptjet,ieta);
192 fHistJetHTT[icent][iptjet][ieta]=new TH2F(name,name,72,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
193 fOutput->Add(fHistJetHTT[icent][iptjet][ieta]);
81941f29 194
1101e468 195 }
196 }
197 }
198
22d3f9c2 199 UInt_t cifras = 0; // bit coded, see GetDimParams() below
b1121dd6 200 if(fDoLessSparseAxes) {
201 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
202 } else {
203 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
204 }
22d3f9c2 205 fhnJH = NewTHnSparseF("fhnJH", cifras);
22d3f9c2 206 fhnJH->Sumw2();
f569a5d2 207 fOutput->Add(fhnJH);
22d3f9c2 208
55f64c2f 209 if(fDoEventMixing){
b1121dd6 210 if(fDoLessSparseAxes) {
211 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5;
212 } else {
213 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
214 }
c1076e58 215 fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
c1076e58 216 fhnMixedEvents->Sumw2();
f569a5d2 217 fOutput->Add(fhnMixedEvents);
965c985f 218 }
219
f569a5d2 220 fOutput->Add(fHistTrackPt);
221 fOutput->Add(fHistCentrality);
222 fOutput->Add(fHistJetEtaPhi);
223 fOutput->Add(fHistJetHEtaPhi);
1101e468 224
f569a5d2 225 PostData(1, fOutput);
55f64c2f 226
227 //Event Mixing
228 Int_t trackDepth = fMixingTracks;
229 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
230
dc64e292 231 Int_t nZvtxBins = 5+1+5;
55f64c2f 232 // bins for second buffer are shifted by 100 cm
dc64e292 233 Double_t vertexBins[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, };
55f64c2f 234 Double_t* zvtxbin = vertexBins;
235
b1121dd6 236// Int_t nCentralityBins = 100;
237 Int_t nCentralityBins = 100;
238 Double_t mult = 1.0;
239 if(fCentBinSize==1) {
240 nCentralityBins = 100;
241 mult = 1.0;
242 } else if(fCentBinSize==2){
243 nCentralityBins = 50;
244 mult = 2.0;
245 } else if(fCentBinSize==5){
246 nCentralityBins = 20;
247 mult = 5.0;
248 } else if(fCentBinSize==10){
249 nCentralityBins = 10;
250 mult = 10.0;
251 }
55f64c2f 252 Double_t centralityBins[nCentralityBins];
253 for(Int_t ic=0; ic<nCentralityBins; ic++){
b1121dd6 254 centralityBins[ic]=mult*ic;
55f64c2f 255 }
55f64c2f 256
257 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
258
1101e468 259}
260
261//________________________________________________________________________
262
263Double_t AliAnalysisTaskEmcalJetHMEC:: RelativePhi(Double_t mphi,Double_t vphi) {
264
1101e468 265 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
266 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
267 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
268 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
269 Double_t dphi = vphi-mphi;
270 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
271 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
272
273 return dphi;//dphi in [-Pi, Pi]
274}
275
1101e468 276//________________________________________________________________________
f569a5d2 277Int_t AliAnalysisTaskEmcalJetHMEC::GetCentBin(Double_t cent) const {
1101e468 278 // Get centrality bin.
279
280 Int_t centbin = -1;
f569a5d2 281 if (cent>=0 && cent<10) centbin = 0;
282 else if (cent>=10 && cent<20) centbin = 1;
283 else if (cent>=20 && cent<30) centbin = 2;
284 else if (cent>=30 && cent<40) centbin = 3;
285 else if (cent>=40 && cent<50) centbin = 4;
286 else if (cent>=50 && cent<90) centbin = 5;
1101e468 287 return centbin;
288}
289
1101e468 290//________________________________________________________________________
f569a5d2 291Int_t AliAnalysisTaskEmcalJetHMEC::GetEtaBin(Double_t eta) const {
1101e468 292 // Get eta bin for histos.
293
294 Int_t etabin = -1;
f569a5d2 295 if (TMath::Abs(eta)<=0.4) etabin = 0;
296 else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8) etabin = 1;
297 else if (TMath::Abs(eta)>=0.8) etabin = 2;
1101e468 298 return etabin;
299}
300//________________________________________________________________________
81941f29 301Int_t AliAnalysisTaskEmcalJetHMEC::GetpTjetBin(Double_t pt) const
1101e468 302{
81941f29 303 // Get jet pt bin for histos.
304
305 Int_t ptbin = -1;
f569a5d2 306 if (pt>=15 && pt<20) ptbin = 0;
307 else if (pt>=20 && pt<25) ptbin = 1;
308 else if (pt>=25 && pt<30) ptbin = 2;
309 else if (pt>=30 && pt<60) ptbin = 3;
310 else if (pt>=60) ptbin = 4;
81941f29 311
312 return ptbin;
1101e468 313}
314
f569a5d2 315//________________________________________________________________________
316void AliAnalysisTaskEmcalJetHMEC::ExecOnce() {
317 AliAnalysisTaskEmcalJet::ExecOnce();
318
319}
1101e468 320
321//________________________________________________________________________
f569a5d2 322Bool_t AliAnalysisTaskEmcalJetHMEC::Run() {
323 // Main loop called for each event
324 if(!fTracks){
325 AliError(Form("No fTracks object!!\n"));
326 return kTRUE;
327 }
328 if(!fJets){
329 AliError(Form("No fJets object!!\n"));
330 return kTRUE;
331 }
1101e468 332
f569a5d2 333 // what kind of event do we have: AOD or ESD?
334 Bool_t esdMode = kTRUE;
335 if (dynamic_cast<AliAODEvent*>(InputEvent())) esdMode = kFALSE;
1101e468 336
f569a5d2 337 // if we have ESD event, set up ESD object
338 if(esdMode){
339 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
340 if (!fESD) {
341 AliError(Form("ERROR: fESD not available\n"));
342 return kTRUE;
343 }
344 }
1101e468 345
f569a5d2 346 // if we have AOD event, set up AOD object
347 if(!esdMode){
348 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
349 if(!fAOD) {
350 AliError(Form("ERROR: fAOD not available\n"));
351 return kTRUE;
352 }
353 }
1101e468 354
f569a5d2 355 TList *list = InputEvent()->GetList();
356 if(!list) {
357 AliError(Form("ERROR: list not attached\n"));
358 return kTRUE;
1101e468 359 }
1101e468 360
f569a5d2 361 // get centrality
362 if (fCent<0) {
363 AliError(Form("Centrality negative: %f", fCent));
364 return kTRUE;
1101e468 365 }
366
f569a5d2 367 Int_t centbin = GetCentBin(fCent);
368 if(centbin<0) return kTRUE;
369
9e4a0350 370 Double_t fvertex[3]={0,0,0};
371 InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
372 Double_t zVtx=fvertex[2];
373
f569a5d2 374 if(fabs(zVtx)>10.0) return kTRUE;
1101e468 375
f569a5d2 376 fHistCentrality->Fill(fCent);
1af24178 377
1101e468 378 TClonesArray *jets = 0;
379 TClonesArray *tracks = 0;
380
f569a5d2 381 tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
1101e468 382 if (!tracks) {
f569a5d2 383 AliError(Form("Pointer to tracks %s == 0", fTracks->GetName() ));
384 return kTRUE;
1101e468 385 }
386 const Int_t Ntracks=tracks->GetEntries();
387
f569a5d2 388 jets= dynamic_cast<TClonesArray*>(list->FindObject(fJets));
9e4a0350 389 if (!jets) {
f569a5d2 390 AliError(Form("Pointer to tracks %s == 0", fJets->GetName() ));
391 return kTRUE;
9e4a0350 392 }
1101e468 393 const Int_t Njets = jets->GetEntries();
394
81941f29 395 //Leticia's loop to find hardest track
81941f29 396 Int_t iTT=-1;
397 Double_t ptmax=-10;
398
f569a5d2 399 for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++) {
400 AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
401 if (!track) {
402 printf("ERROR: Could not receive track %d\n", iTracks);
403 continue;
404 }
81941f29 405
f569a5d2 406 if(TMath::Abs(track->Eta())>0.9) continue;
407 if(track->Pt()<0.15) continue;
408 //iCount++;
409 if(track->Pt()>ptmax){
410 ptmax=track->Pt();
411 iTT=iTracks;
81941f29 412 }
f569a5d2 413 }
1101e468 414
81941f29 415 Int_t ijethi=-1;
81941f29 416 Double_t highestjetpt=0.0;
81941f29 417 Int_t passedTTcut=0;
1101e468 418
f569a5d2 419 for (Int_t ijets = 0; ijets < Njets; ijets++){
420 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijets));
421 if (!jet) continue;
422 if(!AcceptthisJet(jet)) continue;
1101e468 423
1101e468 424 Double_t jetPt = jet->Pt();
425
59d5ae99 426 if(highestjetpt<jetPt){
f569a5d2 427 ijethi=ijets;
428 highestjetpt=jetPt;
59d5ae99 429 }
f569a5d2 430 }
59d5ae99 431
b1121dd6 432 // see if event is selected
433 UInt_t trig = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
434
f569a5d2 435 for (Int_t ijet = 0; ijet < Njets; ijet++){
c1076e58 436 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
f569a5d2 437 if (!jet) continue;
b1121dd6 438
439 // see if event is selected and our jet came from trigger event
440 if (!(trig & fTriggerEventType)) continue;
441 if (jet->Pt()<0.1) continue;
442
f569a5d2 443 if(!AcceptthisJet(jet)) continue;
c1076e58 444
55f64c2f 445 Double_t jetphi = jet->Phi();
59d5ae99 446 Double_t jetPt = jet->Pt();
55f64c2f 447 Double_t jeteta=jet->Eta();
1101e468 448
41f87d59 449 Double_t leadjet=0;
450 if (ijet==ijethi) leadjet=1;
451
965c985f 452 fHistJetPt[centbin]->Fill(jet->Pt());
41f87d59 453 fHistLeadJetPt[centbin]->Fill(jet->Pt());
965c985f 454
41f87d59 455 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
965c985f 456 fHistJetPtBias[centbin]->Fill(jet->Pt());
41f87d59 457 fHistLeadJetPtBias[centbin]->Fill(jet->Pt());
458 }
81941f29 459
f569a5d2 460 fHistJetEtaPhi->Fill(jet->Eta(),jetphi);
1101e468 461
f569a5d2 462 if(iTT>0){
463 AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
464 if(TMath::Abs(jetphi-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
465 else passedTTcut=0;
466 }
1101e468 467
f569a5d2 468 if(passedTTcut)
469 fHistJetPtTT[centbin]->Fill(jet->Pt());
1101e468 470
daf7dbf5 471 Int_t iptjet=-1;
472 iptjet=GetpTjetBin(jetPt);
473 if(iptjet<0) continue;
81941f29 474
b1121dd6 475 if (jetPt > 15) {
daf7dbf5 476
f569a5d2 477 for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++) {
478 AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
479 if (!track) {
480 printf("ERROR: Could not receive track %d\n", iTracks);
481 continue;
482 }
81941f29 483
f569a5d2 484 if(TMath::Abs(track->Eta())>fTrkEta) continue;
1101e468 485
f569a5d2 486 fHistTrackPt->Fill(track->Pt());
1101e468 487
f569a5d2 488 if (track->Pt()<0.15) continue;
1101e468 489
f569a5d2 490 Double_t trackphi = track->Phi();
491 if (trackphi > TMath::Pi())
492 trackphi = trackphi-2*TMath::Pi();
493
494 Double_t tracketa=track->Eta();
495 Double_t trackpt=track->Pt();
496 Double_t deta=tracketa-jeteta;
497 Int_t ieta=GetEtaBin(deta);
498 if (ieta<0) {
499 AliError(Form("Eta Bin negative: %f", deta));
500 continue;
501 }
daf7dbf5 502
f569a5d2 503 //Jet pt, track pt, dPhi,deta,fCent
504 Double_t dphijh = RelativePhi(jetphi,trackphi);
505 if (dphijh < -0.5*TMath::Pi())
506 dphijh+= 2*TMath::Pi();
507 if (dphijh > 1.5*TMath::Pi())
508 dphijh-=2.*TMath::Pi();
1101e468 509
f569a5d2 510 fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
511 fHistJetHEtaPhi->Fill(deta,dphijh);
22d3f9c2 512
f569a5d2 513 Double_t dR=sqrt(deta*deta+dphijh*dphijh);
1101e468 514
f569a5d2 515 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
516 fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,trackpt);
22d3f9c2 517
b1121dd6 518 if(fDoLessSparseAxes) { // check if we want all dimensions
519 Double_t triggerEntries[6] = {fCent,jetPt,trackpt,deta,dphijh,leadjet};
520 fhnJH->Fill(triggerEntries);
521 } else {
522 Double_t triggerEntries[8] = {fCent,jetPt,trackpt,deta,dphijh,leadjet,0.0,dR};
523 fhnJH->Fill(triggerEntries);
524 }
f569a5d2 525 }
22d3f9c2 526
f569a5d2 527 if(passedTTcut)
528 fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,trackpt);
55f64c2f 529
f569a5d2 530 } //track loop
531 }//jet pt cut
c1076e58 532 }//jet loop
55f64c2f 533
1af24178 534 //Prepare to do event mixing
535
55f64c2f 536 // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
b1121dd6 537 TObjArray* tracksClone = 0x0;
55f64c2f 538
55f64c2f 539 if(fDoEventMixing>0){
540
541 // event mixing
542
543 // 1. First get an event pool corresponding in mult (cent) and
544 // zvertex to the current event. Once initialized, the pool
545 // should contain nMix (reduced) events. This routine does not
546 // pre-scan the chain. The first several events of every chain
547 // will be skipped until the needed pools are filled to the
548 // specified depth. If the pool categories are not too rare, this
549 // should not be a problem. If they are rare, you could lose
550 // statistics.
551
552 // 2. Collect the whole pool's content of tracks into one TObjArray
553 // (bgTracks), which is effectively a single background super-event.
554
555 // 3. The reduced and bgTracks arrays must both be passed into
556 // FillCorrelations(). Also nMix should be passed in, so a weight
557 // of 1./nMix can be applied.
558
b1121dd6 559 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
560 // if event was not selected (triggered) for any reseason (should never happen) then return
561 if (trigger==0) return kTRUE;
562
f569a5d2 563 AliEventPool* pool = fPoolMgr->GetEventPool(fCent, zVtx);
55f64c2f 564
daf7dbf5 565 if (!pool){
f569a5d2 566 AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
567 return kTRUE;
daf7dbf5 568 }
55f64c2f 569
b1121dd6 570if(trigger & fTriggerEventType) {
1af24178 571 //check for a trigger jet
b1121dd6 572 if (pool->IsReady() || pool->NTracksInPool() > fNMIXtracks || pool->GetCurrentNEvents() >= fNMIXevents) {
1af24178 573
c1076e58 574 for (Int_t ijet = 0; ijet < Njets; ijet++){
c1076e58 575 Double_t leadjet=0;
576 if (ijet==ijethi) leadjet=1;
1af24178 577
c1076e58 578 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
b1121dd6 579 if (!jet) continue;
580
f569a5d2 581 if(!AcceptthisJet(jet)) continue;
c1076e58 582
583 Double_t jetPt = jet->Pt();
1af24178 584 Double_t jetphi = jet->Phi();
1af24178 585 Double_t jeteta=jet->Eta();
b1121dd6 586
587 // make sure event contains jet above our threshold (reduce stats of sparse)
588 if (jetPt < 15) continue;
589
1af24178 590 Int_t nMix = pool->GetCurrentNEvents();
591
965c985f 592 //Fill for biased jet triggers only
593 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
594
595 // Fill mixed-event histos here
f569a5d2 596 for (Int_t jMix=0; jMix<nMix; jMix++) {
597 TObjArray* bgTracks = pool->GetEvent(jMix);
598 const Int_t Nbgtrks = bgTracks->GetEntries();
965c985f 599
f569a5d2 600 for(Int_t ibg=0; ibg<Nbgtrks; ibg++){
601 AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
602 if(!part) continue;
b1121dd6 603 if(TMath::Abs(part->Eta())>0.9) continue;
604 if(part->Pt()<0.15) continue;
f569a5d2 605
606 Double_t DEta = part->Eta()-jeteta;
607 Double_t DPhi = RelativePhi(jetphi,part->Phi());
608
609 Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
610 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
611 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
b1121dd6 612 if(fDoLessSparseAxes) { // check if we want all the axis filled
613 Double_t triggerEntries[6] = {fCent,jetPt,part->Pt(),DEta,DPhi,leadjet};
614 fhnMixedEvents->Fill(triggerEntries,1./nMix);
615 } else {
616 Double_t triggerEntries[8] = {fCent,jetPt,part->Pt(),DEta,DPhi,leadjet,0.0,DR};
617 fhnMixedEvents->Fill(triggerEntries,1./nMix);
618 }
f569a5d2 619 }
1af24178 620 }
f569a5d2 621 }
965c985f 622 }
55f64c2f 623 }
b1121dd6 624}
625
626 if(trigger & fMixingEventType) {
627 tracksClone = CloneAndReduceTrackList(tracks);
628
629 //update pool if jet in event or not
630 pool->UpdatePool(tracksClone);
631 }
55f64c2f 632
b1121dd6 633 } // end of event mixing
59d5ae99 634
f569a5d2 635 return kTRUE;
1101e468 636}
637
638//________________________________________________________________________
639void AliAnalysisTaskEmcalJetHMEC::Terminate(Option_t *)
640{
641 //just terminate
1101e468 642}
643
c1076e58 644//________________________________________________________________________
f569a5d2 645Int_t AliAnalysisTaskEmcalJetHMEC::AcceptthisJet(AliEmcalJet *jet)
c1076e58 646{
647 //applies all jet cuts except pt
648 float jetphi = jet->Phi();
649 if (jetphi>TMath::Pi())
650 jetphi = jetphi-2*TMath::Pi();
f569a5d2 651
c1076e58 652 if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
653 return 0;
654 if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
655 return 0;
656 if (jet->Area()<fAreacut)
657 return 0;
658 //prevents 0 area jets from sneaking by when area cut == 0
659 if (jet->Area()==0)
f569a5d2 660 return 0;
c1076e58 661 //exclude jets with extremely high pt tracks which are likely misreconstructed
662 if(jet->MaxTrackPt()>100)
663 return 0;
664
665 //passed all above cuts
666 return 1;
c1076e58 667}
668
669//________________________________________________________________________
670
f569a5d2 671THnSparse* AliAnalysisTaskEmcalJetHMEC::NewTHnSparseF(const char* name, UInt_t entries){
55f64c2f 672 // generate new THnSparseF, axes are defined in GetDimParams()
673
674 Int_t count = 0;
675 UInt_t tmp = entries;
676 while(tmp!=0){
677 count++;
678 tmp = tmp &~ -tmp; // clear lowest bit
679 }
680
681 TString hnTitle(name);
682 const Int_t dim = count;
683 Int_t nbins[dim];
684 Double_t xmin[dim];
685 Double_t xmax[dim];
686
687 Int_t i=0;
688 Int_t c=0;
689 while(c<dim && i<32){
690 if(entries&(1<<i)){
691
692 TString label("");
693 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
694 hnTitle += Form(";%s",label.Data());
695 c++;
696 }
697
698 i++;
699 }
700 hnTitle += ";";
701
702 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
703}
704
705void AliAnalysisTaskEmcalJetHMEC::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
706{
707 // stores label and binning of axis for THnSparse
708
709 const Double_t pi = TMath::Pi();
710
711 switch(iEntry){
712
713 case 0:
714 label = "V0 centrality (%)";
55f64c2f 715 nbins = 10;
716 xmin = 0.;
717 xmax = 100.;
718 break;
719
55f64c2f 720 case 1:
721 label = "corrected jet pt";
f569a5d2 722 nbins = 20;
723 xmin = 0.;
724 xmax = 200.;
725 break;
55f64c2f 726
727 case 2:
b1121dd6 728 if(fDoWiderTrackBin) {
729 label = "track pT";
730 nbins = 40;
731 xmin = 0.;
732 xmax = 10.;
733 } else {
734 label = "track pT";
735 nbins = 100;
736 xmin = 0.;
737 xmax = 10;
738 }
55f64c2f 739 break;
b1121dd6 740
741 case 3:
55f64c2f 742 label = "deltaEta";
22d3f9c2 743 nbins = 24;
744 xmin = -1.2;
745 xmax = 1.2;
55f64c2f 746 break;
747
b1121dd6 748 case 4:
55f64c2f 749 label = "deltaPhi";
9e4a0350 750 nbins = 72;
55f64c2f 751 xmin = -0.5*pi;
752 xmax = 1.5*pi;
f569a5d2 753 break;
b1121dd6 754
755 case 5:
756 label = "leading jet";
757 nbins = 3;
758 xmin = -0.5;
759 xmax = 2.5;
55f64c2f 760 break;
761
b1121dd6 762 case 6:
55f64c2f 763 label = "trigger track";
764 nbins =10;
765 xmin = 0;
766 xmax = 50;
767 break;
768
b1121dd6 769 case 7:
770 label = "deltaR";
771 nbins = 10;
772 xmin = 0.;
773 xmax = 5.0;
774 break;
775
c1076e58 776 case 8:
b1121dd6 777 label = "leading track";
778 nbins = 13;
779 xmin = 0;
780 xmax = 50;
c1076e58 781 break;
55f64c2f 782 }
55f64c2f 783}
784
55f64c2f 785//_________________________________________________
786// From CF event mixing code PhiCorrelations
787TObjArray* AliAnalysisTaskEmcalJetHMEC::CloneAndReduceTrackList(TObjArray* tracks)
788{
22d3f9c2 789 // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
55f64c2f 790
791 TObjArray* tracksClone = new TObjArray;
792 tracksClone->SetOwner(kTRUE);
793
794 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
795 {
796 AliVParticle* particle = (AliVParticle*) tracks->At(i);
22d3f9c2 797 if(TMath::Abs(particle->Eta())>fTrkEta) continue;
f569a5d2 798 if(particle->Pt()<0.15) continue;
965c985f 799
22d3f9c2 800 Double_t trackpt=particle->Pt();
801
802 Int_t hadbin=-1;
803 if(trackpt<0.5) hadbin=0;
804 else if(trackpt<1) hadbin=1;
805 else if(trackpt<2) hadbin=2;
806 else if(trackpt<3) hadbin=3;
807 else if(trackpt<5) hadbin=4;
808 else if(trackpt<8) hadbin=5;
809 else if(trackpt<20) hadbin=6;
810
811 if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi());
812
5b38fc6d 813 tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
55f64c2f 814 }
815
816 return tracksClone;
817}