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