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