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