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