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