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