]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadEPpid.cxx
Go from pointer to ifstream to ifstream.
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalJetHadEPpid.cxx
CommitLineData
fc392675 1// Jet-Hadron Correlations PID
2// Event plane dependence task.
3//
4// Author: J Mazer
5
6// task head include
7#include "AliAnalysisTaskEmcalJetHadEPpid.h"
8
9// general ROOT includes
10#include <TCanvas.h>
11#include <TChain.h>
12#include <TClonesArray.h>
13#include <TH1F.h>
14#include <TH2F.h>
15#include <TH3F.h>
16#include <THnSparse.h>
17#include <TList.h>
18#include <TLorentzVector.h>
19#include <TParameter.h>
20#include <TParticle.h>
21#include <TTree.h>
22#include <TVector3.h>
23#include <TObjArray.h>
24
25// AliROOT includes
26#include "AliAODEvent.h"
27#include "AliESDEvent.h"
28#include "AliAnalysisManager.h"
29#include "AliAnalysisTask.h"
30#include "AliCentrality.h"
31#include "AliEmcalJet.h"
32#include "AliAODJet.h"
33#include "AliVCluster.h"
34#include "AliVTrack.h"
35#include <AliVEvent.h>
36#include <AliVParticle.h>
37#include "AliRhoParameter.h"
38#include "AliEmcalParticle.h"
39
40// Localized Rho includes
41#include "AliLocalRhoParameter.h"
42#include "AliAnalysisTaskLocalRho.h"
43
44// event handler (and pico's) includes
45#include <AliInputEventHandler.h>
46#include <AliVEventHandler.h>
47#include "AliESDInputHandler.h"
48#include "AliPicoTrack.h"
49#include "AliEventPoolManager.h"
50
51// PID includes
52#include "AliPIDResponse.h"
53#include "AliTPCPIDResponse.h"
54#include "AliESDpid.h"
55
56// magnetic field includes
57#include "TGeoGlobalMagField.h"
58#include "AliMagF.h"
59
60ClassImp(AliAnalysisTaskEmcalJetHadEPpid)
61
62//________________________________________________________________________
63AliAnalysisTaskEmcalJetHadEPpid::AliAnalysisTaskEmcalJetHadEPpid() :
64 AliAnalysisTaskEmcalJet("correlations",kFALSE),
65 fPhimin(-10), fPhimax(10),
66 fEtamin(-0.9), fEtamax(0.9),
67 fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9), fJetPtcut(15.0),
68 fDoEventMixing(0), fMixingTracks(50000),
69 doPlotGlobalRho(kFALSE), doVariableBinning(kFALSE),
70 dovarbinTHnSparse(kFALSE), useAOD(kTRUE), doPID(kTRUE),
71 fLocalRhoVal(0),
72 fTracksName(""), fJetsName(""),
73 event(0),
74 isPItpc(0), isKtpc(0), isPtpc(0), nPIDtpc(0) , // pid TPC
75 isPIits(0), isKits(0), isPits(0), nPIDits(0) , // pid ITS
76 isPItof(0), isKtof(0), isPtof(0), nPIDtof(0) , // pid TOF
77 fPoolMgr(0x0),
78 fPIDResponse(0x0), fTPCResponse(),
79 fESD(0), fAOD(0),
80 fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
81 fHistRhovsCent(0), fHistNjetvsCent(0), fHistCentrality(0),
82 fHistZvtx(0), fHistMult(0),
83 fHistJetPhi(0), fHistTrackPhi(0),
84 fHistJetHaddPhiIN(0), fHistJetHaddPhiOUT(0), fHistJetHaddPhiMID(0),
85 fHistJetHaddPhiBias(0), fHistJetHaddPhiINBias(0), fHistJetHaddPhiOUTBias(0), fHistJetHaddPhiMIDBias(0),
86 fHistMEdPHI(0),
87 fHistTrackPtallcent(0),
88 fHistJetEtaPhi(0), fHistJetHEtaPhi(0),
89 fHistSEphieta(0), fHistMEphieta(0),
90 fHistJetHaddPHI(0),
91 fHistPID(0),
92 fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0)
93{
94 // Default Constructor
95 for(Int_t ilab=0; ilab<4; ilab++){
96 for(Int_t ipta=0; ipta<7; ipta++){
97 fHistTrackEtaPhi[ilab][ipta]=0;
98 } // end of pt-associated loop
99 } // end of lab loop
100
101 for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
102 fHistJetHadbindPhi[itrackpt]=0;
103 fHistJetHadbindPhiIN[itrackpt]=0;
104 fHistJetHadbindPhiMID[itrackpt]=0;
105 fHistJetHadbindPhiOUT[itrackpt]=0;
106 } // end of trackpt bin loop
107
108 for(Int_t icent = 0; icent<2; ++icent){
109 fHistJetPt[icent]=0;
110 fHistJetPtBias[icent]=0;
111 fHistJetPtTT[icent]=0;
112 fHistAreavsRawPt[icent]=0;
113
114 for(Int_t iptjet = 0; iptjet<5; ++iptjet){
115 for(Int_t ieta = 0; ieta<3; ++ieta){
116 fHistJetH[icent][iptjet][ieta]=0;
117 fHistJetHBias[icent][iptjet][ieta]=0;
118 fHistJetHTT[icent][iptjet][ieta]=0;
119 } // end of eta loop
120 } // end of pt-jet loop
121 } // end of centrality loop
122
123 // centrality dependent histo's
124 for (Int_t i = 0;i<6;++i){
125 fHistJetPtvsTrackPt[i] = 0;
126 fHistRawJetPtvsTrackPt[i] = 0;
127 fHistTrackPt[i] = 0;
128 fHistEP0[i] = 0;
129 fHistEP0A[i] = 0;
130 fHistEP0C[i] = 0;
131 fHistEPAvsC[i] = 0;
132 fHistJetPtvsdEP[i] = 0;
133 fHistJetPtvsdEPBias[i] = 0;
134 fHistRhovsdEP[i] = 0;
135 fHistJetEtaPhiPt[i] = 0;
136 fHistJetEtaPhiPtBias[i] = 0;
137 fHistJetPtArea[i] = 0;
138 fHistJetPtAreaBias[i] = 0;
139 fHistJetPtNcon[i] = 0;
140 fHistJetPtNconBias[i] = 0;
141 fHistJetPtNconCh[i] = 0;
142 fHistJetPtNconBiasCh[i] = 0;
143 fHistJetPtNconEm[i] = 0;
144 fHistJetPtNconBiasEm[i] = 0;
145 fHistJetHaddPhiINcent[i] = 0;
146 fHistJetHaddPhiOUTcent[i] = 0;
147 fHistJetHaddPhiMIDcent[i] = 0;
148 }
149
150 SetMakeGeneralHistograms(kTRUE);
151
152 // define input and output slots here
153 DefineInput(0, TChain::Class());
154 DefineOutput(1, TList::Class());
155}
156
157//________________________________________________________________________
158AliAnalysisTaskEmcalJetHadEPpid::AliAnalysisTaskEmcalJetHadEPpid(const char *name) :
159 AliAnalysisTaskEmcalJet(name,kTRUE),
160 fPhimin(-10), fPhimax(10),
161 fEtamin(-0.9), fEtamax(0.9),
162 fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9), fJetPtcut(15.0),
163 fDoEventMixing(0), fMixingTracks(50000),
164 doPlotGlobalRho(kFALSE), doVariableBinning(kFALSE),
165 dovarbinTHnSparse(kFALSE), useAOD(kTRUE), doPID(kTRUE),
166 fLocalRhoVal(0),
167 fTracksName(""), fJetsName(""),
168 event(0),
169 isPItpc(0), isKtpc(0), isPtpc(0), nPIDtpc(0), // pid TPC
170 isPIits(0), isKits(0), isPits(0), nPIDits(0), // pid ITS
171 isPItof(0), isKtof(0), isPtof(0), nPIDtof(0), // pid TOF
172 fPoolMgr(0x0),
173 fPIDResponse(0x0), fTPCResponse(),
174 fESD(0), fAOD(0),
175 fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
176 fHistRhovsCent(0), fHistNjetvsCent(0), fHistCentrality(0),
177 fHistZvtx(0), fHistMult(0),
178 fHistJetPhi(0), fHistTrackPhi(0),
179 fHistJetHaddPhiIN(0), fHistJetHaddPhiOUT(0), fHistJetHaddPhiMID(0),
180 fHistJetHaddPhiBias(0), fHistJetHaddPhiINBias(0), fHistJetHaddPhiOUTBias(0), fHistJetHaddPhiMIDBias(0),
181 fHistMEdPHI(0),
182 fHistTrackPtallcent(0),
183 fHistJetEtaPhi(0), fHistJetHEtaPhi(0),
184 fHistSEphieta(0), fHistMEphieta(0),
185 fHistJetHaddPHI(0),
186 fHistPID(0),
187 fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0)
188{
189 // Default Constructor
190 for(Int_t ilab=0; ilab<4; ilab++){
191 for(Int_t ipta=0; ipta<7; ipta++){
192 fHistTrackEtaPhi[ilab][ipta]=0;
193 } // end of pt-associated loop
194 } // end of lab loop
195
196 for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
197 fHistJetHadbindPhi[itrackpt]=0;
198 fHistJetHadbindPhiIN[itrackpt]=0;
199 fHistJetHadbindPhiMID[itrackpt]=0;
200 fHistJetHadbindPhiOUT[itrackpt]=0;
201 } // end of trackpt bin loop
202
203 for(Int_t icent = 0; icent<2; ++icent){
204 fHistJetPt[icent]=0;
205 fHistJetPtBias[icent]=0;
206 fHistJetPtTT[icent]=0;
207 fHistAreavsRawPt[icent]=0;
208
209 for(Int_t iptjet = 0; iptjet<5; ++iptjet){
210 for(Int_t ieta = 0; ieta<3; ++ieta){
211 fHistJetH[icent][iptjet][ieta]=0;
212 fHistJetHBias[icent][iptjet][ieta]=0;
213 fHistJetHTT[icent][iptjet][ieta]=0;
214 } // end of eta loop
215 } // end of pt-jet loop
216 } // end of centrality loop
217
218 // centrality dependent histo's
219 for (Int_t i = 0;i<6;++i){
220 fHistJetPtvsTrackPt[i] = 0;
221 fHistRawJetPtvsTrackPt[i] = 0;
222 fHistTrackPt[i] = 0;
223 fHistEP0[i] = 0;
224 fHistEP0A[i] = 0;
225 fHistEP0C[i] = 0;
226 fHistEPAvsC[i] = 0;
227 fHistJetPtvsdEP[i] = 0;
228 fHistJetPtvsdEPBias[i] = 0;
229 fHistRhovsdEP[i] = 0;
230 fHistJetEtaPhiPt[i] = 0;
231 fHistJetEtaPhiPtBias[i] = 0;
232 fHistJetPtArea[i] = 0;
233 fHistJetPtAreaBias[i] = 0;
234 fHistJetPtNcon[i] = 0;
235 fHistJetPtNconBias[i] = 0;
236 fHistJetPtNconCh[i] = 0;
237 fHistJetPtNconBiasCh[i] = 0;
238 fHistJetPtNconEm[i] = 0;
239 fHistJetPtNconBiasEm[i] = 0;
240 fHistJetHaddPhiINcent[i] = 0;
241 fHistJetHaddPhiOUTcent[i] = 0;
242 fHistJetHaddPhiMIDcent[i] = 0;
243 }
244
245 SetMakeGeneralHistograms(kTRUE);
246
247 // define input and output slots here
248 DefineInput(0, TChain::Class());
249 DefineOutput(1, TList::Class());
250}
251
252//_______________________________________________________________________
253AliAnalysisTaskEmcalJetHadEPpid::~AliAnalysisTaskEmcalJetHadEPpid()
254{
255 // destructor
256 if (fOutput) {
257 delete fOutput;
258 fOutput = 0;
259 }
260}
261
262//________________________________________________________________________
263void AliAnalysisTaskEmcalJetHadEPpid::UserCreateOutputObjects()
264{ // This is called ONCE!
265 if (!fCreateHisto) return;
266 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
267 OpenFile(1); // do I need the 1?
268
269 // create histograms that arn't array
270 fHistTPCdEdX = new TH2F("TPCdEdX", "TPCdEdX", 2000, 0.0, 100.0, 500, 0, 500);
271 fHistITSsignal = new TH2F("ITSsignal", "ITSsignal", 2000, 0.0, 100.0, 500, 0, 500);
272// fHistTOFsignal = new TH2F("TOFsignal", "TOFsignal", 2000, 0.0, 100.0, 500, 0, 500);
273 fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 500, 0, 500);
274 fHistNjetvsCent = new TH2F("NjetvsCent", "NjetvsCent", 100, 0.0, 100.0, 100, 0, 100);
275 fHistTrackPtallcent = new TH1F("fHistTrackPtallcent", "p_{T} distribution", 1000, 0.0, 100.0);
276 fHistCentrality = new TH1F("fHistCentrality","centrality",100,0,100);
277 fHistZvtx = new TH1F("fHistZvertex","z vertex",60,-30,30);
278 fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet #eta-#phi",900,-1.8,1.8,640,-3.2,3.2);
279 fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron #Delta#eta-#Delta#phi",900,-1.8,1.8,640,-1.6,4.8);
280 fHistSEphieta = new TH2F("fHistSEphieta", "Single Event #phi-#eta distribution", 64,-0.5*TMath::Pi(), 1.5*TMath::Pi(), 32,-1.8,1.8); // was 64 bins
281 fHistMEphieta = new TH2F("fHistMEphieta", "Mixed Event #phi-#eta distribution", 64, -0.5*TMath::Pi(), 1.5*TMath::Pi(), 32,-1.8,1.8); // was 64 bins
282 fHistJetHaddPHI = new TH1F("fHistJetHaddPHI", "Jet-Hadron #Delta#varphi", 128,-0.5*TMath::Pi(),1.5*TMath::Pi());
283 fHistJetPhi = new TH1F("fHistJetPhi", "Jet #phi Distribution", 128, -2.0*TMath::Pi(), 2.0*TMath::Pi());
284 fHistTrackPhi = new TH1F("fHistTrackPhi", "Track #phi Distribution", 128, -2.0*TMath::Pi(), 2.0*TMath::Pi());
285 fHistJetHaddPhiIN = new TH1F("fHistJetHaddPhiIN","Jet-Hadron #Delta#varphi IN PLANE", 128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
286 fHistJetHaddPhiOUT = new TH1F("fHistJetHaddPhiOUT","Jet-Hadron #Delta#varphi OUT PLANE",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
287 fHistJetHaddPhiMID = new TH1F("fHistJetHaddPhiMID","Jet-Hadron #Delta#varphi MIDDLE of PLANE",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
288 fHistJetHaddPhiBias = new TH1F("fHistJetHaddPhiBias","Jet-Hadron #Delta#varphi with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
289 fHistJetHaddPhiINBias = new TH1F("fHistJetHaddPhiINBias","Jet-Hadron #Delta#varphi IN PLANE with bias", 128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
290 fHistJetHaddPhiOUTBias = new TH1F("fHistJetHaddPhiOUTBias","Jet-Hadron #Delta#varphi OUT PLANE with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
291 fHistJetHaddPhiMIDBias = new TH1F("fHistJetHaddPhiMIDBias","Jet-Hadron #Delta#varphi MIDDLE of PLANE with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
292
293 fHistPID = new TH1F("fHistPID","PID Counter",12, -0.5, 11.5);
294 SetfHistPIDcounterLabels(fHistPID);
295 fOutput->Add(fHistPID);
296
297 char name[200];
298
299 for(Int_t ipta=0; ipta<7; ++ipta){
300 for(Int_t ilab=0; ilab<4; ++ilab){
81b592c4 301 snprintf(name,200, "fHistTrackEtaPhi_%i_%i", ilab,ipta);
fc392675 302 fHistTrackEtaPhi[ilab][ipta] = new TH2F(name,name,400,-1,1,640,0.0,2.*TMath::Pi());
303/// fOutput->Add(fHistTrackEtaPhi[ilab][ipta]);
304 } // end of lab loop
305 } // end of pt-associated loop
306
307
308 // jet hadron (binned) correlations in dPHI
309 for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
81b592c4 310 snprintf(name,200,"fHistJetHadbindPhi_%i",itrackpt);
fc392675 311 fHistJetHadbindPhi[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
312 fOutput->Add(fHistJetHadbindPhi[itrackpt]);
313
81b592c4 314 snprintf(name,200,"fHistJetHadbindPhiIN_%i",itrackpt);
fc392675 315 fHistJetHadbindPhiIN[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
316 fOutput->Add(fHistJetHadbindPhiIN[itrackpt]);
317
81b592c4 318 snprintf(name,200,"fHistJetHadbindPhiMID_%i",itrackpt);
fc392675 319 fHistJetHadbindPhiMID[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
320 fOutput->Add(fHistJetHadbindPhiMID[itrackpt]);
321
81b592c4 322 snprintf(name,200,"fHistJetHadbindPhiOUT_%i",itrackpt);
fc392675 323 fHistJetHadbindPhiOUT[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
324 fOutput->Add(fHistJetHadbindPhiOUT[itrackpt]);
325 } // end of trackpt bin loop
326
327 for(Int_t icent = 0; icent<2; ++icent){
81b592c4 328 snprintf(name,200,"fHistJetPt_%i",icent);
fc392675 329 fHistJetPt[icent] = new TH1F(name,name,200,0,200);
330 fOutput->Add(fHistJetPt[icent]);
331
81b592c4 332 snprintf(name,200,"fHistJetPtBias_%i",icent);
fc392675 333 fHistJetPtBias[icent] = new TH1F(name,name,200,0,200);
334 fOutput->Add(fHistJetPtBias[icent]);
335
81b592c4 336 snprintf(name,200,"fHistJetPtTT_%i",icent);
fc392675 337 fHistJetPtTT[icent] = new TH1F(name,name,200,0,200);
338 fOutput->Add(fHistJetPtTT[icent]);
339
81b592c4 340 snprintf(name,200,"fHistAreavsRawPt_%i",icent);
fc392675 341 fHistAreavsRawPt[icent] = new TH2F(name,name,100,0,1,200,0,200);
342 //fOutput->Add(fHistAreavsRawPt[icent]);
343
344 for(Int_t iptjet = 0; iptjet<5; ++iptjet){
345 for(Int_t ieta = 0; ieta<3; ++ieta){
81b592c4 346 snprintf(name,200,"fHistJetH_%i_%i_%i",icent,iptjet,ieta);
fc392675 347 fHistJetH[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
348 fOutput->Add(fHistJetH[icent][iptjet][ieta]);
349
81b592c4 350 snprintf(name,200,"fHistJetHBias_%i_%i_%i",icent,iptjet,ieta);
fc392675 351 fHistJetHBias[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
352 fOutput->Add(fHistJetHBias[icent][iptjet][ieta]);
353
81b592c4 354 snprintf(name,200,"fHistJetHTT_%i_%i_%i",icent,iptjet,ieta);
fc392675 355 fHistJetHTT[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
356 fOutput->Add(fHistJetHTT[icent][iptjet][ieta]);
357 } // end of eta loop
358 } // end of pt-jet loop
359 } // end of centrality loop
360
361 // strings for titles
362 TString name1;
363 TString title1;
364
365 for (Int_t i = 0;i<6;++i){
366 name1 = TString(Form("JetPtvsTrackPt_%i",i));
367 title1 = TString(Form("Jet p_{T} vs Leading Track p_{T} cent bin %i",i));
368 fHistJetPtvsTrackPt[i] = new TH2F(name1,title1,250,-50,200,250,0,50);
369 fOutput->Add(fHistJetPtvsTrackPt[i]);
370
371 name1 = TString(Form("RawJetPtvsTrackPt_%i",i));
372 title1 = TString(Form("Raw Jet p_{T} vs Leading Track p_{T} cent bin %i",i));
373 fHistRawJetPtvsTrackPt[i] = new TH2F(name1,title1,250,-50,200,250,0,50);
374 fOutput->Add(fHistRawJetPtvsTrackPt[i]);
375
376 name1 = TString(Form("TrackPt_%i",i));
377 title1 = TString(Form("Track p_{T} cent bin %i",i));
378 fHistTrackPt[i] = new TH1F(name1,title1,1000,0,100); // up to 200?
379 fOutput->Add(fHistTrackPt[i]);
380
381 name1 = TString(Form("EP0_%i",i));
382 title1 = TString(Form("EP VZero cent bin %i",i));
383 fHistEP0[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
384 fOutput->Add(fHistEP0[i]);
385
386/*
387 name1 = TString(Form("EP0A_%i",i));
388 title1 = TString(Form("EP VZero cent bin %i",i));
389 fHistEP0A[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
390 fOutput->Add(fHistEP0A[i]);
391
392 name1 = TString(Form("EP0C_%i",i));
393 title1 = TString(Form("EP VZero cent bin %i",i));
394 fHistEP0C[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
395 fOutput->Add(fHistEP0C[i]);
396
397 name1 = TString(Form("EPAvsC_%i",i));
398 title1 = TString(Form("EP VZero cent bin %i",i));
399 fHistEPAvsC[i] = new TH2F(name1,title1,144,-TMath::Pi(),TMath::Pi(),144,-TMath::Pi(),TMath::Pi());
400 fOutput->Add(fHistEPAvsC[i]);
401*/
402
403 name1 = TString(Form("JetPtvsdEP_%i",i));
404 title1 = TString(Form("Jet p_{T} vs #DeltaEP cent bin %i",i));
405 fHistJetPtvsdEP[i] = new TH2F(name1,title1,250,-50,200,288,-2*TMath::Pi(),2*TMath::Pi());
406 fOutput->Add(fHistJetPtvsdEP[i]);
407
408 name1 = TString(Form("JetPtvsdEPBias_%i",i));
409 title1 = TString(Form("Bias Jet p_{T} vs #DeltaEP cent bin %i",i));
410 fHistJetPtvsdEPBias[i] = new TH2F(name1,title1,250,-50,200,288,-2*TMath::Pi(),2*TMath::Pi());
411 fOutput->Add(fHistJetPtvsdEPBias[i]);
412
413 name1 = TString(Form("RhovsdEP_%i",i));
414 title1 = TString(Form("#rho vs #DeltaEP cent bin %i",i));
415 fHistRhovsdEP[i] = new TH2F(name1,title1,500,0,500,288,-2*TMath::Pi(),2*TMath::Pi());
416 fOutput->Add(fHistRhovsdEP[i]);
417
418 name1 = TString(Form("JetEtaPhiPt_%i",i));
419 title1 = TString(Form("Jet #eta-#phi p_{T} cent bin %i",i));
420 fHistJetEtaPhiPt[i] = new TH3F(name1,title1,250,-50,200,100,-1,1,64,-3.2,3.2);
421 fOutput->Add(fHistJetEtaPhiPt[i]);
422
423 name1 = TString(Form("JetEtaPhiPtBias_%i",i));
424 title1 = TString(Form("Jet #eta-#phi p_{T} Bias cent bin %i",i));
425 fHistJetEtaPhiPtBias[i] = new TH3F(name1,title1,250,-50,200,100,-1,1,64,-3.2,3.2);
426 fOutput->Add(fHistJetEtaPhiPtBias[i]);
427
428 name1 = TString(Form("JetPtArea_%i",i));
429 title1 = TString(Form("Jet p_{T} Area cent bin %i",i));
430 fHistJetPtArea[i] = new TH2F(name1,title1,250,-50,200,100,0,1);
431 fOutput->Add(fHistJetPtArea[i]);
432
433 name1 = TString(Form("JetPtAreaBias_%i",i));
434 title1 = TString(Form("Jet p_{T} Area Bias cent bin %i",i));
435 fHistJetPtAreaBias[i] = new TH2F(name1,title1,250,-50,200,100,0,1);
436 fOutput->Add(fHistJetPtAreaBias[i]);
437
438 name1 = TString(Form("JetPtNcon_%i",i));
439 title1 = TString(Form("Jet p_{T} Ncon cent bin %i",i));
440 fHistJetPtNcon[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
441 //fOutput->Add(fHistJetPtNcon[i]);
442
443 name1 = TString(Form("JetPtNconBias_%i",i));
444 title1 = TString(Form("Jet p_{T} NconBias cent bin %i",i));
445 fHistJetPtNconBias[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
446 fOutput->Add(fHistJetPtNconBias[i]);
447
448 name1 = TString(Form("JetPtNconCh_%i",i));
449 title1 = TString(Form("Jet p_{T} NconCh cent bin %i",i));
450 fHistJetPtNconCh[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
451 //fOutput->Add(fHistJetPtNconCh[i]);
452
453 name1 = TString(Form("JetPtNconBiasCh_%i",i));
454 title1 = TString(Form("Jet p_{T} NconBiasCh cent bin %i",i));
455 fHistJetPtNconBiasCh[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
456 //fOutput->Add(fHistJetPtNconBiasCh[i]);
457
458 name1 = TString(Form("JetPtNconEm_%i",i));
459 title1 = TString(Form("Jet p_{T} NconEm cent bin %i",i));
460 fHistJetPtNconEm[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
461 //fOutput->Add(fHistJetPtNconEm[i]);
462
463 name1 = TString(Form("JetPtNconBiasEm_%i",i));
464 title1 = TString(Form("Jet p_{T} NconBiasEm cent bin %i",i));
465 fHistJetPtNconBiasEm[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
466 //fOutput->Add(fHistJetPtNconBiasEm[i]);
467
468 name1 = TString(Form("JetHaddPhiINcent_%i",i));
469 title1 = TString(Form("Jet Hadron #Delta#varphi Distribution IN PLANE cent bin %i",i));
470 fHistJetHaddPhiINcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
471 fOutput->Add(fHistJetHaddPhiINcent[i]);
472
473 name1 = TString(Form("JetHaddPhiOUTcent_%i",i));
474 title1 = TString(Form("Jet Hadron #Delta#varphi Distribution OUT PLANE cent bin %i",i));
475 fHistJetHaddPhiOUTcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
476 fOutput->Add(fHistJetHaddPhiOUTcent[i]);
477
478 name1 = TString(Form("JetHaddPhiMIDcent_%i",i));
479 title1 = TString(Form("Jet Hadron #Delta#varphi Distribution MIDDLE of PLANE cent bin %i",i));
480 fHistJetHaddPhiMIDcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
481 fOutput->Add(fHistJetHaddPhiMIDcent[i]);
482 }
483
484 // add histo's to output list
485 fOutput->Add(fHistTPCdEdX);
486 fOutput->Add(fHistITSsignal);
487// fOutput->Add(fHistTOFsignal);
488 fOutput->Add(fHistRhovsCent);
489 fOutput->Add(fHistNjetvsCent);
490 fOutput->Add(fHistTrackPtallcent);
491 fOutput->Add(fHistJetEtaPhi);
492 //fOutput->Add(fHistTrackEtaPhi);
493 fOutput->Add(fHistJetHEtaPhi);
494 fOutput->Add(fHistCentrality);
495 fOutput->Add(fHistZvtx);
496 fOutput->Add(fHistJetPhi);
497 fOutput->Add(fHistTrackPhi);
498 fOutput->Add(fHistJetHaddPhiIN);
499 fOutput->Add(fHistJetHaddPhiOUT);
500 fOutput->Add(fHistJetHaddPhiMID);
501 fOutput->Add(fHistJetHaddPhiBias);
502 fOutput->Add(fHistJetHaddPhiINBias);
503 fOutput->Add(fHistJetHaddPhiOUTBias);
504 fOutput->Add(fHistJetHaddPhiMIDBias);
505 fOutput->Add(fHistSEphieta);
506 fOutput->Add(fHistMEphieta);
507 fOutput->Add(fHistJetHaddPHI);
508
509 // set up jet-hadron sparse
510 UInt_t bitcoded = 0; // bit coded, see GetDimParams() below
511 UInt_t cifras = 0;
512 UInt_t bitcode = 0; // bit coded, see GetDimParamsPID() below
513 bitcoded = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9;
514 fhnJH = NewTHnSparseF("fhnJH", bitcoded);
515 fhnJH->Sumw2();
516/// fOutput->Add(fhnJH);
517
518 // ****************************** PID *****************************************************
519 // set up PID handler
520 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
521 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
522 if(!inputHandler) {
81b592c4 523 AliError("Input handler needed");
524 return;
fc392675 525 }
526
527 // PID response object
528 fPIDResponse = inputHandler->GetPIDResponse();
529 if (!fPIDResponse) {
530 AliError("PIDResponse object was not created");
81b592c4 531 return;
fc392675 532 }
533
534 // ****************************************************************************************
535/*
536 // for pp we need mult bins for event mixing. Create binning here, to also make a histogram from it
537 Int_t nCentralityBins = 8;
538 Double_t centralityBins[9] = {0.0, 4., 9, 15, 25, 35, 55, 100.0,500.0};
539 Double_t centralityBins[nCentralityBins+1];
540 for(Int_t ic=0; ic<nCentralityBins+1; ic++){
541 if(ic==nCentralityBins) centralityBins[ic]=500;
542 else centralityBins[ic]=10.0*ic;
543 }
544*/
545 // setup for Pb-Pb collisions
546 Int_t nCentralityBins = 100;
547 Double_t centralityBins[nCentralityBins];
548 for(Int_t ic=0; ic<nCentralityBins; ic++){
549 centralityBins[ic]=1.0*ic;
550 }
551
552 fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBins,centralityBins);
553 fOutput->Add(fHistMult);
554
555 // Event Mixing
556 Int_t trackDepth = fMixingTracks;
557 Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
558 Int_t nZvtxBins = 5+1+5;
559 Double_t vertexBins[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10};
560 Double_t* zvtxbin = vertexBins;
561 fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
562
563 // set up event mixing sparse
564 if(fDoEventMixing){
565 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
566 fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
567 fhnMixedEvents->Sumw2();
568 fOutput->Add(fhnMixedEvents);
569 } // end of do-eventmixing
570
571 Double_t xlowjetPT[] = {-50,-45,-40,-35,-30,-25,-20,-18,-16,-14,-12,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,12,14,16,18,20,25,30,35,40,45,50,60,70,80,90,100,120,140,160,180,200,250,300,350,400};
572 Double_t xlowtrPT[] = {0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.25,2.50,2.75,3.0,3.25,3.5,3.75,4.0,4.25,4.50,4.75,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,22.0,24.0,26.0,28.0,30.0,35.0,40.0,45.0,50.0,60.0,70.0,80.0,90.0,100.0};
573
574 // number of bins you tell histogram should be (# in array - 1) because the last bin
575 // is the right-most edge of the histogram
576 // i.e. this is for PT and there are 57 numbers (bins) thus we have 56 bins in our histo
577 Int_t nbinsjetPT = sizeof(xlowjetPT)/sizeof(Double_t) - 1;
578 Int_t nbinstrPT = sizeof(xlowtrPT)/sizeof(Double_t) - 1;
579
580 // set up PID sparse
581 if(doPID){
582 bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
583 1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14 | 1<<15 | 1<<16 | 1<<17 | 1<<18 | 1<<19 |
584 1<<20 | 1<<21 | 1<<22;
585 fhnPID = NewTHnSparseFPID("fhnPID", bitcode);
586 fhnPID->Sumw2();
587 fOutput->Add(fhnPID);
588 } // end of do-PID
589
590 if(dovarbinTHnSparse){
591 fhnPID->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
592 fhnPID->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
593 }
594
595 // =========== Switch on Sumw2 for all histos ===========
596 for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
597 TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
598 if (h1){
599 h1->Sumw2();
600 continue;
601 }
602 TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
603 if (h2){
604 h2->Sumw2();
605 continue;
606 }
607 TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
608 if (h3){
609 h3->Sumw2();
610 continue;
611 }
612 THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
613 if(hn)hn->Sumw2();
614 }
615
616 PostData(1, fOutput);
617}
618
619//_________________________________________________________________________
620void AliAnalysisTaskEmcalJetHadEPpid::ExecOnce()
621{
622 AliAnalysisTaskEmcalJet::ExecOnce();
623}
624
625//_________________________________________________________________________
626Bool_t AliAnalysisTaskEmcalJetHadEPpid::Run()
627{ // Main loop called for each event
628 // TEST TEST TEST TEST for OBJECTS!
629 if(!fLocalRho){
630 AliError(Form("Couldn't get fLocalRho object, try to get it from Event based on name\n"));
631 fLocalRho = GetLocalRhoFromEvent(fLocalRhoName);
632 }
633 if(!fTracks){
634 AliError(Form("No fTracks object!!\n"));
635 return kTRUE;
636 }
637 if(!fJets){
638 AliError(Form("No fJets object!!\n"));
639 return kTRUE;
640 }
641
642 // what kind of event do we have: AOD or ESD?
643 if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
644 else useAOD = kFALSE;
645
646 // if we have ESD event, set up ESD object
647 if(!useAOD){
648 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
649 if (!fESD) {
650 AliError(Form("ERROR: fESD not available\n"));
651 return kTRUE;
652 }
653 }
654
655 // if we have AOD event, set up AOD object
656 if(useAOD){
657 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
658 if(!fAOD) {
659 AliError(Form("ERROR: fAOD not available\n"));
660 return kTRUE;
661 }
662 }
663
664 // get centrality
665 Int_t centbin = GetCentBin(fCent);
666 fHistCentrality->Fill(fCent); // won't be filled in pp collision (Keep this in mind!)
667
668 // for pp analyses we will just use the first centrality bin
669 if (centbin == -1) centbin = 0;
670
671 // get vertex information
672 Double_t fvertex[3]={0,0,0};
673 InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
674 Double_t zVtx=fvertex[2];
675 fHistZvtx->Fill(zVtx);
676
677 // get z-vertex bin
678 //Int_t zVbin = GetzVertexBin(zVtx);
679
680 // apply zVtx cut
681 if(fabs(zVtx)>10.0)
682 return kTRUE;
683
684 // apply cut to event on Centrality > 90%
685 if(fCent>90) return kTRUE;
686
687 // create pointer to list of input event
688 TList *list = InputEvent()->GetList();
689 if(!list) {
690 AliError(Form("ERROR: list not attached\n"));
81b592c4 691 return kFALSE;
fc392675 692 }
693
694 // initialize TClonesArray pointers to jets and tracks
695 TClonesArray *jets = 0;
696 TClonesArray *tracks = 0;
697
698 // get Tracks object
699 tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
700 if (!tracks) {
701 AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
702 return kTRUE;
703 } // verify existence of tracks
704
705 // get Jets object
706 jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
707 if(!jets){
708 AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
709 return kTRUE;
710 } // verify existence of jets
711
712 // get number of jets and tracks
713 const Int_t Njets = jets->GetEntries();
714 const Int_t Ntracks = tracks->GetEntries();
715 if(Ntracks<1) return kTRUE;
716 if(Njets<1) return kTRUE;
717
718 fHistMult->Fill(Ntracks); // fill multiplicity distribution
719
720 // initialize track parameters
721 Int_t iTT=-1;
722 Double_t ptmax=-10;
723
724 // loop over tracks - to get hardest track (highest pt)
725 for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++){
726 AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
727 if (!track) {
728 AliError(Form("Couldn't get VTrack track %d\n", iTracks));
729 continue;
730 } // verify existence of tracks
731
732 if(TMath::Abs(track->Eta())>0.9) continue;
733 if(track->Pt()<0.15) continue;
734 //iCount++;
735 if(track->Pt()>ptmax){
736 ptmax=track->Pt(); // max pT track
737 iTT=iTracks; // trigger tracks
738 } // check if Pt>maxpt
739
740 fHistTrackPhi->Fill(track->Phi());
741 fHistTrackPtallcent->Fill(track->Pt());
742 } // end of loop over tracks
743
744 // find tracks for event (alternative), and loop over
745 for (int i = 0;i<Ntracks;i++){
746 AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
747 if (!track) continue;
748
749 if(TMath::Abs(track->Eta())>0.9) continue;
750 if(track->Pt()<0.15) continue;
751
752 fHistTrackPt[centbin]->Fill(track->Pt());
753 }
754
755 // event plane histograms filled
756 fHistEP0[centbin]->Fill(fEPV0);
757 //fHistEP0A[centbin]->Fill(fEPV0A);
758 //fHistEP0C[centbin]->Fill(fEPV0C);
759 //fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
760
761 // get rho from event and fill relative histo's
762 fRho = GetRhoFromEvent(fRhoName);
763 fRhoVal = fRho->GetVal();
764 fHistRhovsCent->Fill(fCent,fRhoVal); // Global Rho vs Centrality
765 fHistRhovsdEP[centbin]->Fill(fRhoVal,fEPV0); // Global Rho vs delta Event Plane angle
766
767 // initialize jet parameters
768 Int_t ijethi=-1;
769 Double_t highestjetpt=0.0;
770 Int_t passedTTcut=0;
771 Int_t NjetAcc = 0;
772
773 // loop over jets in an event - to find highest jet pT and apply some cuts
774 for (Int_t ijet = 0; ijet < Njets; ijet++){
775 // get our jets
776 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
777 if (!jet) continue;
778
779 // cuts on jets
780 if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
781 continue;
782 if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
783 continue;
784 fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
785 if (jet->Area()<fAreacut)
786 continue;
787 //prevents 0 area jets from sneaking by when area cut == 0
788 if (jet->Area()==0)
789 continue;
790 //exclude jets with extremely high pt tracks which are likely misreconstructed
791 if(jet->MaxTrackPt()>100)
792 continue;
793 if(!AcceptJet(jet))
794 continue;
795
796 NjetAcc++; // # of accepted jets
797 fHistJetPhi->Fill(jet->Phi()); // Jet Phi histogram (filled)
798
799 // get highest Pt jet in event
800 if(highestjetpt<jet->Pt()){
801 ijethi=ijet;
802 highestjetpt=jet->Pt();
803 }
804 } // end of looping over jets
805
806 // accepted jets
807 fHistNjetvsCent->Fill(fCent,NjetAcc);
808 Int_t NJETAcc = 0;
809
810 // loop over jets in event and make appropriate cuts
811 for (Int_t ijet = 0; ijet < Njets; ++ijet) {
812 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ijet));
813 if (!jet) // see if we have a jet
814 continue;
815 if(!AcceptJet(jet)) // do we accept this jet?
816 continue;
817
818// Double_t leadjet=0;
819// if (ijet==ijethi) leadjet=1;
820
821 // cuts on jet
822 if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
823 continue;
824 if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
825 continue;
826 if (jet->Area()==0) // make sure jet has an area
827 continue;
828 if (jet->Pt()<0.1) // (should probably be higher..., but makes a cut on jet pT)
829 continue;
830 if (jet->MaxTrackPt()>100) // elminates fake tracks
831 continue;
832 if (!AcceptJet(jet)) // sees if jet is accepted
833 continue;
834 // jets.push_back(jet);
835
836 // check on lead jet
837 Double_t leadjet=0;
838 if (ijet==ijethi) leadjet=1;
839
840 // initialize and calculate various parameters: pt, eta, phi, rho, etc...
841 Double_t jetphi = jet->Phi(); // phi of jet
842 NJETAcc++; // # accepted jets
843 fLocalRhoVal = fLocalRho->GetLocalVal(jetphi, GetJetRadius(0)); //fmyJetRadius); // get local rho value
844 Double_t jeteta = jet->Eta(); // ETA of jet
845// Double_t jetptraw = jet->Pt(); // raw pT of jet
846 Double_t jetPt = -500; // initialize corr jet pt
847// jetPt = jet->Pt();
848// jetPt = jet->Pt()-jet->Area()*fRhoVal; // corrected pT of jet from rho value
849 jetPt = jet->Pt()-jet->Area()*fLocalRhoVal; // corrected pT of jet using Rho modulated for V2 and V3
850 Float_t dEP = -500; // initialize angle between jet and event plane
851 //Int_t ieta = -500; // initialize jet eta bin
852 //ieta = GetEtaBin(jeteta); // bin of eta
853 dEP = RelativeEPJET(jetphi,fEPV0); // angle betweeen jet and event plane
854
855 // make histo's
856 fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
857 fHistRawJetPtvsTrackPt[centbin]->Fill(jet->Pt(),jet->MaxTrackPt());
858// fHistJetPtvsdEP[centbin]->Fill(jetPt,RelativePhi((fEPV0+TMath::Pi()),jet->Phi()));
859 fHistJetPtvsdEP[centbin]->Fill(jetPt, dEP);
860 fHistJetEtaPhiPt[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
861 fHistJetPtArea[centbin]->Fill(jetPt,jet->Area());
862 fHistJetPtNcon[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
863 fHistJetPtNconCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
864 fHistJetPtNconEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
865
866 // make histo's with BIAS applied
867 if (jet->MaxTrackPt()>fTrkBias){
868 fHistJetPtvsdEPBias[centbin]->Fill(jetPt, dEP);
869 fHistJetEtaPhiPtBias[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
870 fHistJetPtAreaBias[centbin]->Fill(jetPt,jet->Area());
871 fHistJetPtNconBias[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
872 fHistJetPtNconBiasCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
873 fHistJetPtNconBiasEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
874 }
875
876 // fill some histo's
877 fHistJetPt[centbin]->Fill(jet->Pt()); // fill #jets vs pT histo
878 fHistJetEtaPhi->Fill(jet->Eta(),jet->Phi()); // fill jet eta-phi distribution histo
879 //fHistDeltaPtvsArea->Fill(jetPt,jet->Area());
880
881 if(leadjet && centbin==0) fHistJetPt[centbin+1]->Fill(jet->Pt());
882 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
883 fHistJetPtBias[centbin]->Fill(jet->Pt());
884 if(leadjet && centbin==0) fHistJetPtBias[centbin+1]->Fill(jet->Pt());
885 } // end of MaxTrackPt>ftrkBias or maxclusterPt>fclusBias
886
887 // do we have trigger tracks
888 if(iTT>0){
889 AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
890 if(TMath::Abs(jet->Phi()-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
891 else passedTTcut=0;
892 } // end of check on iTT > 0
893
894 if(passedTTcut) fHistJetPtTT[centbin]->Fill(jet->Pt());
895
896 // cut on HIGHEST jet pt in event (15 GeV default)
897 if (highestjetpt>fJetPtcut) {
898
899 // loop over all track for an event containing a jet with a pT>fJetPtCut (15)GeV
900 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
901 //AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
902 // check on even type to get object for PID info
903 //char eventtype[10];
904
905 AliVTrack *track = fAOD->GetTrack(iTracks);
906 if(!track) track = fESD->GetTrack(iTracks);
907// AliAODTrack* track = fAOD->GetTrack(iTracks);
908// //sprintf(eventtype,"fAOD");
909// if(!track) AliESDtrack* track = fESD->GetTrack(iTracks);
910// //sprintf(eventtype,"fESD");
911
912 if (!track) {
913 AliError(Form("Couldn't get AliESDtrack %d\n", iTracks));
914 continue;
915 }
916
917// AliESDtrack *track = GetAcceptTrack(esdtrack);
918// if(!track) continue;
919
920 // apply track cuts
921 if(TMath::Abs(track->Eta())>fTrkEta)
922 continue;
923 //fHistTrackPt->Fill(track->Pt()); // fill track pT histogram
924 if (track->Pt()<0.15) // make sure track pT > 0.15 GeV
925 continue;
926
927 // calculate and get some track parameters
928 track->GetTPCsignal(); // get the TPC signal of the track
929 Double_t tracketa=track->Eta(); // eta of track
930 Double_t deta=tracketa-jeteta; // dETA between track and jet
931 Int_t ieta=GetEtaBin(deta); // bin of eta
81b592c4 932 if(ieta<0) continue;
fc392675 933
934 // dPHI between jet and hadron
935 Double_t dphijh = RelativePhi(jet->Phi(), track->Phi()); // angle between jet and hadron
936
937 // jet pt bins
81b592c4 938 //Int_t iptjet=-1; // initialize jet pT bin
fc392675 939// iptjet=GetpTjetBin(jetPt); // bin of jet pT
81b592c4 940 Int_t iptjet=GetpTjetBin(jet->Pt());
941 if(iptjet<0) continue;
fc392675 942 Double_t dR=sqrt(deta*deta+dphijh*dphijh); // difference of R between jet and hadron track
943
944 // fill some jet-hadron histo's
945 fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt()); // fill jet-hadron dPHI--track pT distribution
946 fHistJetHEtaPhi->Fill(deta,dphijh); // fill jet-hadron eta--phi distribution
947 fHistJetHaddPHI->Fill(dphijh);
948
949 // does our max track or cluster pass the bias?
950 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
951 // set up and fill Jet-Hadron THnSparse
952 Double_t triggerEntries[9] = {fCent,jetPt,track->Pt(),dR,deta,dphijh,0,leadjet,zVtx};
953 fhnJH->Fill(triggerEntries); // fill Sparse Histo with trigger entries
954
955 // fill histo's
956 fHistSEphieta->Fill(dphijh, deta);
957 fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
958 fHistJetHaddPhiBias->Fill(dphijh);
959
960 // in plane and out of plane histo's
961 if( dEP>0 && dEP<=(TMath::Pi()/6) ){
962 // we are IN plane
963 fHistJetHaddPhiINBias->Fill(dphijh);
964 }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
965 // we are OUT of PLANE
966 fHistJetHaddPhiOUTBias->Fill(dphijh);
967 }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
968 // we are in middle of plane
969 fHistJetHaddPhiMIDBias->Fill(dphijh);
970 }
971 } // end of check maxtrackpt>ftrackbias or maxclusterpt>fclustbias
972
973 if(passedTTcut) fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
974
975// **************************************************************************************************************
976// *********************************** PID **********************************************************************
977// **************************************************************************************************************
978 if(doPID){
979 // data information for pid
980 Double_t eta, pt, EovP, fClsE = -1.;
981
982 EovP = 0;
983 fClsE = 0;
984
985 // nSigma of particles in TPC, TOF, and ITS
986 Double_t nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC = -1.;
987 Double_t nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF = -1.;
988 Double_t nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS = -1.;
989
990 // misc quantities
991 Double_t dEdx, ITSsig, TOFsig = -1.;
992 Double_t charge = -99;
993 //Double_t mom=track->GetTPCmomentum();
994
995 // get parameters of track
996 charge = track->Charge(); // charge of track
997 eta = track->Eta(); // ETA of track
998 pt = track->Pt(); // pT of track
999
1000 // extra attempt
1001 AliVEvent *vevent=InputEvent();
1002 if (!vevent||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
1003
1004 // get impact parameters and detector signals
1005 dEdx = track->GetTPCsignal();
1006 ITSsig = track->GetITSsignal();
1007 TOFsig = track->GetTOFsignal();
1008
1009 // fill detector signal histograms
1010 fHistTPCdEdX->Fill(pt, dEdx);
1011 fHistITSsignal->Fill(pt, ITSsig);
1012// fHistTOFsignal->Fill(pt, TOFsig);
1013
1014 // TPC nSigma's
1015 nSigmaPion_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
1016 nSigmaKaon_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kKaon);
1017 nSigmaProton_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
1018
1019 // TOF nSigma's
1020 nSigmaPion_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
1021 nSigmaKaon_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kKaon);
1022 nSigmaProton_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
1023
1024 // ITS nSigma's
1025 nSigmaPion_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kPion);
1026 nSigmaKaon_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kKaon);
1027 nSigmaProton_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kProton);
1028
1029 // Tests to PID pions, kaons, and protons, (default is undentified tracks)
1030 //Double_t nPIDtpc, nPIDits, nPIDtof = 0;
1031 Double_t nPID = -10;
1032
1033 // check track has pT < 0.900 GeV - use TPC pid
1034 if (pt<0.900 && dEdx>0) {
1035 nPIDtpc = 4;
1036
1037 // PION check - TPC
1038 if (TMath::Abs(nSigmaPion_TPC)<2 && TMath::Abs(nSigmaKaon_TPC)>2 && TMath::Abs(nSigmaProton_TPC)>2 ){
1039 isPItpc = kTRUE;
1040 nPIDtpc = 1;
1041 nPID=0.5;
1042 }else isPItpc = kFALSE;
1043
1044 // KAON check - TPC
1045 if (TMath::Abs(nSigmaKaon_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaProton_TPC)>2 ){
1046 isKtpc = kTRUE;
1047 nPIDtpc = 2;
1048 nPID=1.5;
1049 }else isKtpc = kFALSE;
1050
1051 // PROTON check - TPC
1052 if (TMath::Abs(nSigmaProton_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaKaon_TPC)>2 ){
1053 isPtpc = kTRUE;
1054 nPIDtpc = 3;
1055 nPID=2.5;
1056 }else isPtpc = kFALSE;
1057 } // cut on track pT for TPC
1058
1059 // check track has pT < 0.500 GeV - use ITS pid
1060 if (pt<0.500 && ITSsig>0) {
1061 nPIDits = 4;
1062
1063 // PION check - ITS
1064 if (TMath::Abs(nSigmaPion_ITS)<2 && TMath::Abs(nSigmaKaon_ITS)>2 && TMath::Abs(nSigmaProton_ITS)>2 ){
1065 isPIits = kTRUE;
1066 nPIDits = 1;
1067 nPID=4.5;
1068 }else isPIits = kFALSE;
1069
1070 // KAON check - ITS
1071 if (TMath::Abs(nSigmaKaon_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaProton_ITS)>2 ){
1072 isKits = kTRUE;
1073 nPIDits = 2;
1074 nPID=5.5;
1075 }else isKits = kFALSE;
1076
1077 // PROTON check - ITS
1078 if (TMath::Abs(nSigmaProton_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaKaon_ITS)>2 ){
1079 isPits = kTRUE;
1080 nPIDits = 3;
1081 nPID=6.5;
1082 }else isPits = kFALSE;
1083 } // cut on track pT for ITS
1084
1085 // check track has 0.900 GeV < pT < 2.500 GeV - use TOF pid
1086 if (pt>0.900 && pt<2.500 && TOFsig>0) {
1087 nPIDtof = 4;
1088
1089 // PION check - TOF
1090 if (TMath::Abs(nSigmaPion_TOF)<2 && TMath::Abs(nSigmaKaon_TOF)>2 && TMath::Abs(nSigmaProton_TOF)>2 ){
1091 isPItof = kTRUE;
1092 nPIDtof = 1;
1093 nPID=8.5;
1094 }else isPItof = kFALSE;
1095
1096 // KAON check - TOF
1097 if (TMath::Abs(nSigmaKaon_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaProton_TOF)>2 ){
1098 isKtof = kTRUE;
1099 nPIDtof = 2;
1100 nPID=9.5;
1101 }else isKtof = kFALSE;
1102
1103 // PROTON check - TOF
1104 if (TMath::Abs(nSigmaProton_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaKaon_TOF)>2 ){
1105 isPtof = kTRUE;
1106 nPIDtof = 3;
1107 nPID=10.5;
1108 }else isPtof = kFALSE;
1109 } // cut on track pT for TOF
1110
1111 if (nPIDtpc == 0) nPID = -0.5;
1112 if (nPIDits == 0) nPID = 3.5;
1113 if (nPIDtof == 0) nPID = 7.5;
1114 fHistPID->Fill(nPID);
1115
1116 // PID sparse getting filled
1117 Double_t nontrig_tracks_Entries[23] = {fCent,jetPt,pt,charge,eta,deta,dphijh,leadjet,zVtx,EovP,fClsE,
1118 nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC, //nSig in TPC
1119 nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS, //nSig in ITS
1120 nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF, //nSig in TOF
1121 nPIDtpc, nPIDits, nPIDtof //PID label for each detector
1122 }; //array for PID sparse
1123 fhnPID->Fill(nontrig_tracks_Entries); // fill Sparse histo of PID tracks
1124 } // end of doPID check
1125
1126 // get track pt bin
1127 Int_t itrackpt = -500; // initialize track pT bin
1128 itrackpt = GetpTtrackBin(track->Pt());
1129
1130 // all tracks: jet hadron correlations in hadron pt bins
1131 fHistJetHadbindPhi[itrackpt]->Fill(dphijh);
1132
1133 // in plane and out of plane jet-hadron histo's
1134 if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1135 // we are IN plane
1136 fHistJetHaddPhiINcent[centbin]->Fill(dphijh);
1137 fHistJetHaddPhiIN->Fill(dphijh);
1138 fHistJetHadbindPhiIN[itrackpt]->Fill(dphijh);
1139 //fHistJetHaddPhiPtcentbinIN[itrackpt][centbin]->Fill(dphijh);
1140 }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1141 // we are OUT of PLANE
1142 fHistJetHaddPhiOUTcent[centbin]->Fill(dphijh);
1143 fHistJetHaddPhiOUT->Fill(dphijh);
1144 fHistJetHadbindPhiOUT[itrackpt]->Fill(dphijh);
1145 //fHistJetHaddPhiPtcentbinOUT[itrackpt][centbin]->Fill(dphijh);
1146 }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1147 // we are in the middle of plane
1148 fHistJetHaddPhiMIDcent[centbin]->Fill(dphijh);
1149 fHistJetHaddPhiMID->Fill(dphijh);
1150 fHistJetHadbindPhiMID[itrackpt]->Fill(dphijh);
1151 }
1152 } // loop over tracks found in event with highest JET pT > 10.0 GeV (change)
1153 } // jet pt cut
1154 } // jet loop
1155
1156// ***************************************************************************************************************
1157// ******************************** Event MIXING *****************************************************************
1158// TObjArray* tracksClone = CloneAndReduceTrackList(tracks); // TEST
1159
1160 //Prepare to do event mixing
1161 if(fDoEventMixing>0){
1162 // event mixing
1163
1164 // 1. First get an event pool corresponding in mult (cent) and
1165 // zvertex to the current event. Once initialized, the pool
1166 // should contain nMix (reduced) events. This routine does not
1167 // pre-scan the chain. The first several events of every chain
1168 // will be skipped until the needed pools are filled to the
1169 // specified depth. If the pool categories are not too rare, this
1170 // should not be a problem. If they are rare, you could lose
1171 // statistics.
1172
1173 // 2. Collect the whole pool's content of tracks into one TObjArray
1174 // (bgTracks), which is effectively a single background super-event.
1175
1176 // 3. The reduced and bgTracks arrays must both be passed into
1177 // FillCorrelations(). Also nMix should be passed in, so a weight
1178 // of 1./nMix can be applied.
1179
1180 // mix jets from triggered events with tracks from MB events
1181 // get the trigger bit
1182 // need to change trigger bits between different runs
1183 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1184 if (trigger==0) return kTRUE; // return
1185
1186 //Double_t Ntrks=(Double_t)Ntracks*1.0;
1187 //cout<<"Test.. Ntrks: "<<fPoolMgr->GetEventPool(Ntrks);
1188
1189 AliEventPool* pool = fPoolMgr->GetEventPool(fCent, zVtx); // for PbPb? fcent
1190 //AliEventPool* pool = fPoolMgr->GetEventPool(Ntrks, zVtx); // for pp
1191
1192 if (!pool){
1193 AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
1194 //AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f", Ntrks, zVtx));
1195 return kTRUE;
1196 }
1197
1198 // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1199/// if (trigger & AliVEvent::kEMC1) {
1200 if (trigger & AliVEvent::kEMCEJE) { // TEST
1201 //check for a trigger jet
1202 // fmixingtrack/10 ??
1203 if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 /*100*/ || pool->GetCurrentNEvents() >= 5) {
1204 // loop over jets (passing cuts?)
1205 for (Int_t ijet = 0; ijet < Njets; ijet++) {
1206 //Double_t leadjet=0;
1207 //if (ijet==ijethi) leadjet=1;
1208
1209 // get jet object
1210 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1211 if(!AcceptJet(jet)) continue;
1212
1213 // get jet stats
1214 Double_t jetphi = jet->Phi(); // phi of jet
1215 Double_t jetPt = jet->Pt(); // pT of jet
1216
1217 Int_t nMix = pool->GetCurrentNEvents(); // how many particles in pool to mix
1218
1219 // Fill for biased jet triggers only
1220 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) {
1221 // Fill mixed-event histos here
1222 for (Int_t jMix=0; jMix<nMix; jMix++) {
1223 TObjArray* bgTracks = pool->GetEvent(jMix);
1224 const Int_t Nbgtrks = bgTracks->GetEntries();
1225 for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1226 AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1227 if(!part) continue;
1228
1229 //Double_t DPhi = jetphi - part->Phi();
1230 Double_t DEta = part->Eta()-jet->Eta(); // difference in eta
1231 Double_t DPhi = RelativePhi(jet->Phi(),part->Phi()); // difference in phi
1232 Float_t dEP = RelativeEPJET(jetphi,fEPV0);
1233 Double_t charge = part->Charge();
1234
1235 //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta); // difference in R
1236 if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
1237 if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
1238
1239 Double_t triggerEntries[8] = {fCent,jetPt,part->Pt(),DEta,DPhi,dEP,zVtx, charge}; //array for ME sparse
1240 fhnMixedEvents->Fill(triggerEntries,1./nMix); // fill Sparse histo of mixed events
1241
1242 // *********************
1243 //create phi distribution from mixed events
1244 //fHistMEdPHI->Fill(DPhi);
1245 fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1246
1247 } // end of background track loop
1248 } // end of filling mixed-event histo's
1249 } // end of check for biased jet triggers
1250 } // end of jet loop
1251 } // end of check for triggered jet
1252 } //end EMC triggered loop
1253
1254 // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1255/// if (trigger & AliVEvent::kMB) {
1256 if (trigger & AliVEvent::kAnyINT){ // test
1257 // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1258 TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
1259
1260 // update pool if jet in event or not
1261 pool->UpdatePool(tracksClone);
1262 } // check on track from MB events
1263 } // end of event mixing
1264
1265 // print some stats on the event
1266 event++;
1267
1268 //cout<<"Event #: "<<event<<endl;
1269 //cout<<"# of jets: "<<Njets<<" Highest jet pt: "<<highestjetpt<<endl;
1270 //cout<<"# tracks: "<<Ntracks<<" Highest track pt: "<<ptmax<<endl;
1271 //cout<<" =============================================== "<<endl;
1272
1273 return kTRUE; // used when the function is of type bool
1274} // end of RUN
1275
1276//________________________________________________________________________
1277Int_t AliAnalysisTaskEmcalJetHadEPpid::GetCentBin(Double_t cent) const
1278{ // Get centrality bin.
1279 Int_t centbin = -1;
1280 if (cent>=0 && cent<10) centbin = 0;
1281 else if (cent>=10 && cent<20) centbin = 1;
1282 else if (cent>=20 && cent<30) centbin = 2;
1283 else if (cent>=30 && cent<40) centbin = 3;
1284 else if (cent>=40 && cent<50) centbin = 4;
1285 else if (cent>=50 && cent<90) centbin = 5;
1286
1287 return centbin;
1288}
1289
1290//________________________________________________________________________
1291Float_t AliAnalysisTaskEmcalJetHadEPpid::RelativePhi(Double_t mphi,Double_t vphi) const
1292{ // function to calculate relative PHI
1293 double dphi = mphi-vphi;
1294
1295 // set dphi to operate on adjusted scale
1296 if(dphi<-0.5*TMath::Pi()) dphi+=2.*TMath::Pi();
1297 if(dphi>3./2.*TMath::Pi()) dphi-=2.*TMath::Pi();
1298
1299 // test
1300 if( dphi < -1.*TMath::Pi()/2 || dphi > 3.*TMath::Pi()/2 )
1301 AliWarning(Form("%s: dPHI not in range [-0.5*Pi, 1.5*Pi]!", GetName()));
1302
1303 return dphi; // dphi in [-0.5Pi, 1.5Pi]
1304}
1305
1306
1307
1308//_________________________________________________________________________
1309Float_t AliAnalysisTaskEmcalJetHadEPpid:: RelativeEPJET(Double_t jetAng, Double_t EPAng) const
1310{ // function to calculate angle between jet and EP in the 1st quadrant (0,Pi/2)
1311 Double_t dphi = (EPAng - jetAng);
1312
1313 // ran into trouble with a few dEP<-Pi so trying this...
1314 if( dphi<-1*TMath::Pi() ){
1315 dphi = dphi + 1*TMath::Pi();
1316 } // this assumes we are doing full jets currently
1317
1318 if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
1319 // Do nothing! we are in quadrant 1
1320 }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
1321 dphi = 1*TMath::Pi() - dphi;
1322 }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
1323 dphi = fabs(dphi);
1324 }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
1325 dphi = dphi + 1*TMath::Pi();
1326 }
1327
1328 // test
1329 if( dphi < 0 || dphi > TMath::Pi()/2 )
1330 AliWarning(Form("%s: dPHI not in range [0, 0.5*Pi]!", GetName()));
1331
1332 return dphi; // dphi in [0, Pi/2]
1333}
1334
1335//Int_t ieta=GetEtaBin(deta);
1336//________________________________________________________________________
1337Int_t AliAnalysisTaskEmcalJetHadEPpid::GetEtaBin(Double_t eta) const
1338{
1339 // Get eta bin for histos.
1340 Int_t etabin = -1;
1341 if (TMath::Abs(eta)<=0.4) etabin = 0;
1342 else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8) etabin = 1;
1343 else if (TMath::Abs(eta)>=0.8) etabin = 2;
1344
1345 return etabin;
1346} // end of get-eta-bin
1347
1348//________________________________________________________________________
1349Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTjetBin(Double_t pt) const
1350{
1351 // Get jet pt bin for histos.
1352 Int_t ptbin = -1;
1353 if (pt>=15 && pt<20) ptbin = 0;
1354 else if (pt>=20 && pt<25) ptbin = 1;
1355 else if (pt>=25 && pt<40) ptbin = 2;
1356 else if (pt>=40 && pt<60) ptbin = 3;
1357 else if (pt>=60) ptbin = 4;
1358
1359 return ptbin;
1360} // end of get-jet-pt-bin
1361
1362//________________________________________________________________________
1363Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTtrackBin(Double_t pt) const
1364{
1365 // May need to update bins for future runs... (testing locally)
1366
1367 // Get track pt bin for histos.
1368 Int_t ptbin = -1;
1369 if (pt < 0.5) ptbin = 0;
1370 else if (pt>=0.5 && pt<1.0) ptbin = 1;
1371 else if (pt>=1.0 && pt<1.5) ptbin = 2;
1372 else if (pt>=1.5 && pt<2.0) ptbin = 3;
1373 else if (pt>=2.0 && pt<2.5) ptbin = 4;
1374 else if (pt>=2.5 && pt<3.0) ptbin = 5;
1375 else if (pt>=3.0 && pt<4.0) ptbin = 6;
1376 else if (pt>=4.0 && pt<5.0) ptbin = 7;
1377 else if (pt>=5.0) ptbin = 8;
1378
1379 return ptbin;
1380} // end of get-jet-pt-bin
1381
1382
1383//___________________________________________________________________________
1384Int_t AliAnalysisTaskEmcalJetHadEPpid::GetzVertexBin(Double_t zVtx) const
1385{
1386 // get z-vertex bin for histo.
1387 int zVbin= -1;
1388 if (zVtx>=-10 && zVtx<-8) zVbin = 0;
1389 else if (zVtx>=-8 && zVtx<-6) zVbin = 1;
1390 else if (zVtx>=-6 && zVtx<-4) zVbin = 2;
1391 else if (zVtx>=-4 && zVtx<-2) zVbin = 3;
1392 else if (zVtx>=-2 && zVtx<0) zVbin = 4;
1393 else if (zVtx>=0 && zVtx<2) zVbin = 5;
1394 else if (zVtx>=2 && zVtx<4) zVbin = 6;
1395 else if (zVtx>=4 && zVtx<6) zVbin = 7;
1396 else if (zVtx>=6 && zVtx<8) zVbin = 8;
1397 else if (zVtx>=8 && zVtx<10) zVbin = 9;
1398 else zVbin = 10;
1399
1400 return zVbin;
1401} // end of get z-vertex bin
1402
1403//______________________________________________________________________
1404THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseF(const char* name, UInt_t entries)
1405{
1406 // generate new THnSparseF, axes are defined in GetDimParams()
1407 Int_t count = 0;
1408 UInt_t tmp = entries;
1409 while(tmp!=0){
1410 count++;
1411 tmp = tmp &~ -tmp; // clear lowest bit
1412 }
1413
1414 TString hnTitle(name);
1415 const Int_t dim = count;
1416 Int_t nbins[dim];
1417 Double_t xmin[dim];
1418 Double_t xmax[dim];
1419
1420 Int_t i=0;
1421 Int_t c=0;
1422 while(c<dim && i<32){
1423 if(entries&(1<<i)){
1424
1425 TString label("");
1426 GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1427 hnTitle += Form(";%s",label.Data());
1428 c++;
1429 }
1430
1431 i++;
1432 }
1433 hnTitle += ";";
1434
1435 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1436} // end of NewTHnSparseF
1437
1438void AliAnalysisTaskEmcalJetHadEPpid::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1439{
1440 // stores label and binning of axis for THnSparse
1441 const Double_t pi = TMath::Pi();
1442
1443 switch(iEntry){
1444
1445 case 0:
1446 label = "V0 centrality (%)";
1447 nbins = 10;
1448 xmin = 0.;
1449 xmax = 100.;
1450 break;
1451
1452 case 1:
1453 label = "Corrected Jet p_{T}";
1454 nbins = 250;
1455 xmin = -50.;
1456 xmax = 200.;
1457 break;
1458
1459 case 2:
1460 label = "Track p_{T}";
1461 nbins = 100;
1462 xmin = 0.;
1463 xmax = 100.;
1464 break;
1465
1466 case 3:
1467 label = "Relative Eta";
1468 nbins = 64;
1469 xmin = -1.6;
1470 xmax = 1.6;
1471 break;
1472
1473 case 4:
1474 label = "Relative Phi";
1475 nbins = 144;
1476 xmin = -0.5*pi;
1477 xmax = 1.5*pi;
1478 break;
1479
1480 case 5:
1481 label = "Relative angle of Jet and Reaction Plane";
1482 nbins = 144;
1483 xmin = -0.5*pi;
1484 xmax = 1.5*pi;
1485 break;
1486
1487 case 6:
1488 label = "z-vertex";
1489 nbins = 10;
1490 xmin = -10;
1491 xmax = 10;
1492 break;
1493
1494 case 7:
1495 label = "track charge";
1496 nbins =3;
1497 xmin = -1.5;
1498 xmax = 1.5;
1499 break;
1500
1501 case 8: // need to update
1502 label = "leading track";
1503 nbins = 10;
1504 xmin = 0;
1505 xmax = 50;
1506 break;
1507
1508 case 9: // need to update
1509 label = "leading jet";
1510 nbins = 3;
1511 xmin = -0.5;
1512 xmax = 2.5;
1513 break;
1514 } // end of switch
1515} // end of getting dim-params
1516
1517//_________________________________________________
1518// From CF event mixing code PhiCorrelations
1519TObjArray* AliAnalysisTaskEmcalJetHadEPpid::CloneAndReduceTrackList(TObjArray* tracks)
1520{
1521 // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
1522 TObjArray* tracksClone = new TObjArray;
1523 tracksClone->SetOwner(kTRUE);
1524
1525 /*
1526 Int_t nTrax = fESD->GetNumberOfTracks();
1527 //cout << "nTrax " << nTrax <<endl;
1528
1529 for (Int_t i = 0; i < nTrax; ++i)
1530 {
1531 AliESDtrack* esdtrack = fESD->GetTrack(i);
1532 if (!esdtrack)
1533 {
1534 AliError(Form("Couldn't get ESD track %d\n", i));
1535 continue;
1536 }
1537
1538 AliESDtrack *particle = GetAcceptTrack(esdtrack);
1539 if(!particle) continue;
1540
1541 cout << "RM Hybrid track : " << i << " " << particle->Pt() << endl;
1542 */
1543
1544 //cout << "nEntries " << tracks->GetEntriesFast() <<endl;
1545 for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1546 {
1547 AliVParticle* particle = (AliVParticle*) tracks->At(i);
1548 if(TMath::Abs(particle->Eta())>fTrkEta) continue;
1549 if(particle->Pt()<0.15)continue;
1550
1551 Double_t trackpt=particle->Pt(); // track pT
1552
1553 Int_t trklabel=-1;
1554 trklabel=particle->GetLabel();
1555 //cout << "TRACK_LABEL: " << particle->GetLabel()<<endl;
1556
1557 Int_t hadbin=-1;
1558 if(trackpt<0.5) hadbin=0;
1559 else if(trackpt<1) hadbin=1;
1560 else if(trackpt<2) hadbin=2;
1561 else if(trackpt<3) hadbin=3;
1562 else if(trackpt<5) hadbin=4;
1563 else if(trackpt<8) hadbin=5;
1564 else if(trackpt<20) hadbin=6;
1565
1566//feb10 comment out
1567 if(hadbin>-1 && trklabel>-1 && trklabel <3) fHistTrackEtaPhi[trklabel][hadbin]->Fill(particle->Eta(),particle->Phi());
1568 if(hadbin>-1) fHistTrackEtaPhi[3][hadbin]->Fill(particle->Eta(),particle->Phi());
1569
1570 //if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi()); // TEST
1571
1572 tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
1573 } // end of looping through tracks
1574
1575 return tracksClone;
1576}
1577
1578//____________________________________________________________________________________________
1579THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFPID(const char* name, UInt_t entries)
1580{
1581 // generate new THnSparseF PID, axes are defined in GetDimParams()
1582 Int_t count = 0;
1583 UInt_t tmp = entries;
1584 while(tmp!=0){
1585 count++;
1586 tmp = tmp &~ -tmp; // clear lowest bit
1587 }
1588
1589 TString hnTitle(name);
1590 const Int_t dim = count;
1591 Int_t nbins[dim];
1592 Double_t xmin[dim];
1593 Double_t xmax[dim];
1594
1595 Int_t i=0;
1596 Int_t c=0;
1597 while(c<dim && i<32){
1598 if(entries&(1<<i)){
1599
1600 TString label("");
1601 GetDimParamsPID(i, label, nbins[c], xmin[c], xmax[c]);
1602 hnTitle += Form(";%s",label.Data());
1603 c++;
1604 }
1605
1606 i++;
1607 }
1608 hnTitle += ";";
1609
1610 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1611} // end of NewTHnSparseF PID
1612
1613//________________________________________________________________________________
1614void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsPID(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1615{
1616 // stores label and binning of axis for THnSparse
1617 const Double_t pi = TMath::Pi();
1618
1619 switch(iEntry){
1620
1621 case 0:
1622 label = "V0 centrality (%)";
1623 nbins = 10;
1624 xmin = 0.;
1625 xmax = 100.;
1626 break;
1627
1628 case 1:
1629 label = "Corrected Jet p_{T}";
1630 nbins = 250;
1631 xmin = -50.;
1632 xmax = 200.;
1633 break;
1634
1635 case 2:
1636 label = "Track p_{T}";
1637 nbins = 500;
1638 xmin = 0.;
1639 xmax = 50.;
1640 break;
1641
1642 case 3:
1643 label = "Charge of Track";
1644 nbins = 3;
1645 xmin = -1.5;
1646 xmax = 1.5;
1647 break;
1648
1649 case 4:
1650 label = "Track Eta";
1651 nbins = 64;
1652 xmin = -1.2;
1653 xmax = 1.2;
1654 break;
1655
1656 case 5:
1657 label = "Relative Eta of Track and Jet";
1658 nbins = 64;
1659 xmin = -1.6;
1660 xmax = 1.6;
1661 break;
1662
1663 case 6:
1664 label = "Relative Phi of Track and Jet";
1665 nbins = 64;
1666 xmin = -0.5*pi;
1667 xmax = 1.5*pi;
1668 break;
1669
1670 case 7:
1671 label = "leading jet";
1672 nbins = 3;
1673 xmin = -.5;
1674 xmax = 2.5;
1675 break;
1676
1677 case 8:
1678 label = "Z-vertex";
1679 nbins = 10;
1680 xmin = -10.;
1681 xmax = 10.;
1682 break;
1683
1684 case 9: // don't need
1685 label = "E/p of track";
1686 nbins = 100;
1687 xmin = 0.;
1688 xmax = 1.;
1689 break;
1690
1691 case 10: // don't need
1692 label = "Cluster Energy";
1693 nbins = 50; //250;
1694 xmin = 0.;
1695 xmax = 50.; //250.;
1696 break;
1697
1698 case 11:
1699 label = "N-Sigma of pions in TPC";
1700 nbins = 100;
1701 xmin = -5.;
1702 xmax = 5.;
1703 break;
1704
1705 case 12:
1706 label = "N-Sigma of protons in TPC";
1707 nbins = 100;
1708 xmin = -5.;
1709 xmax = 5.;
1710 break;
1711
1712 case 13:
1713 label = "N-Sigma of kaons in TPC";
1714 nbins = 100;
1715 xmin = -5.;
1716 xmax = 5.;
1717 break;
1718
1719 case 14:
1720 label = "N-Sigma of pions in ITS";
1721 nbins = 100;
1722 xmin = -5.;
1723 xmax = 5.;
1724 break;
1725
1726 case 15:
1727 label = "N-Sigma of protons in ITS";
1728 nbins = 100;
1729 xmin = -5.;
1730 xmax = 5.;
1731 break;
1732
1733 case 16:
1734 label = "N-Sigma of kaons in ITS";
1735 nbins = 100;
1736 xmin = -5.;
1737 xmax = 5.;
1738 break;
1739
1740 case 17:
1741 label = "N-Sigma of pions in TOF";
1742 nbins = 100;
1743 xmin = -5.;
1744 xmax = 5.;
1745 break;
1746
1747 case 18:
1748 label = "N-Sigma of protons in TOF";
1749 nbins = 100;
1750 xmin = -5.;
1751 xmax = 5.;
1752 break;
1753
1754 case 19:
1755 label = "N-Sigma of kaons in TOF";
1756 nbins = 100;
1757 xmin = -5.;
1758 xmax = 5.;
1759 break;
1760
1761 case 20:
1762 label = "TPC PID determination";
1763 nbins = 4;
1764 xmin = 0.;
1765 xmax = 4.;
1766 break;
1767
1768 case 21:
1769 label = "ITS PID determination";
1770 nbins = 4;
1771 xmin = 0.;
1772 xmax = 4.;
1773 break;
1774
1775 case 22:
1776 label = "TOF PID determination";
1777 nbins = 4;
1778 xmin = 0.;
1779 xmax = 4.;
1780 break;
1781
1782 } // end of switch
1783} // end of get dimension parameters PID
1784
1785void AliAnalysisTaskEmcalJetHadEPpid::Terminate(Option_t *)
1786{
1787 cout<<"#########################"<<endl;
1788 cout<<"#### DONE RUNNING!!! ####"<<endl;
1789 cout<<"#########################"<<endl;
1790
1791} // end of terminate
1792
1793
1794//void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
1795void AliAnalysisTaskEmcalJetHadEPpid::SetfHistPIDcounterLabels(TH1* h) const
1796{
1797 // fill the analysis summary histrogram, saves all relevant analysis settigns
1798 h->GetXaxis()->SetBinLabel(1, "TPC: Unidentified"); // -0.5
1799 h->GetXaxis()->SetBinLabel(2, "TPC: Pion"); // 0.5
1800 h->GetXaxis()->SetBinLabel(3, "TPC: Kaon"); // 1.5
1801 h->GetXaxis()->SetBinLabel(4, "TPC: Proton"); // 2.5
1802 h->GetXaxis()->SetBinLabel(5, "ITS: Unidentified"); // 3.5
1803 h->GetXaxis()->SetBinLabel(6, "ITS: Pion"); // 4.5
1804 h->GetXaxis()->SetBinLabel(7, "ITS: Kaon"); // 5.5
1805 h->GetXaxis()->SetBinLabel(8, "ITS: Proton"); // 6.5
1806 h->GetXaxis()->SetBinLabel(9, "TOF: Unidentified"); // 7.5
1807 h->GetXaxis()->SetBinLabel(10, "TOF: Pion"); // 8.5
1808 h->GetXaxis()->SetBinLabel(11, "TOF: Kaon"); // 9.5
1809 h->GetXaxis()->SetBinLabel(12, "TOF: Proton"); // 10.5
1810
1811}