]> 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
changed binning for event pool, added more trigger QA histos, added more mixed event...
[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),
a409c5ba 79 fDoEventMixing(0), fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
d60e724e 80 fTriggerEventType(AliVEvent::kAny), fMixingEventType(AliVEvent::kAny),
8ee2d316 81 doPlotGlobalRho(0), doVariableBinning(0), dovarbinTHnSparse(0),
daa728b0 82 makeQAhistos(0), makeBIAShistos(0), makeextraCORRhistos(0), makeoldJEThadhistos(0),
139e13f4 83 allpidAXIS(0), fcutType("EMCAL"), doPID(0), doPIDtrackBIAS(0),
56523939 84 doComments(0),
85 doFlavourJetAnalysis(0), fJetFlavTag(-99),
86 fBeam(kNA),
fc392675 87 fLocalRhoVal(0),
56523939 88 fTracksName(""), fTracksNameME(""), fJetsName(""),
fc392675 89 event(0),
aa5f3717 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
fc392675 93 fPoolMgr(0x0),
94 fPIDResponse(0x0), fTPCResponse(),
d60e724e 95 fESD(0), fAOD(0), fVevent(0),
96 fHistEventQA(0), fHistEventSelectionQA(0),
a409c5ba 97 fHistCentZvertGA(0), fHistCentZvertJE(0), fHistCentZvertMB(0), fHistCentZvertAny(0),
fc392675 98 fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
99 fHistRhovsCent(0), fHistNjetvsCent(0), fHistCentrality(0),
100 fHistZvtx(0), fHistMult(0),
101 fHistJetPhi(0), fHistTrackPhi(0),
d60e724e 102 fHistLocalRhoJetpt(0),
fc392675 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),
daa728b0 111 fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0), fhnCorr(0x0),
8ee2d316 112 fJetsCont(0), fTracksCont(0), fCaloClustersCont(0),
113 fContainerAllJets(0), fContainerPIDJets(1)
fc392675 114{
115 // Default Constructor
116 for(Int_t ilab=0; ilab<4; ilab++){
117 for(Int_t ipta=0; ipta<7; ipta++){
72a4d320 118 //fHistTrackEtaPhi[ilab][ipta]=0; // keep out for now
fc392675 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
daa728b0 129 for(Int_t icent = 0; icent<6; ++icent){
fc392675 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){
daa728b0 141 fHistJetPt[i] = 0;
142 fHistJetPtBias[i] = 0;
143 fHistJetPtTT[i] = 0;
144 fHistAreavsRawPt[i] = 0;
fc392675 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;
8ee2d316 152 fHistJetPtcorrGlRho[i] = 0;
fc392675 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
fc392675 173}
174
175//________________________________________________________________________
176AliAnalysisTaskEmcalJetHadEPpid::AliAnalysisTaskEmcalJetHadEPpid(const char *name) :
177 AliAnalysisTaskEmcalJet(name,kTRUE),
178 fPhimin(-10), fPhimax(10),
179 fEtamin(-0.9), fEtamax(0.9),
f1ce4cc3 180 fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9),
181 fJetPtcut(15.0), fJetRad(0.4), fConstituentCut(0.15),
3263f852 182 fesdTrackCuts(0),
a409c5ba 183 fDoEventMixing(0), fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
d60e724e 184 fTriggerEventType(AliVEvent::kAny), fMixingEventType(AliVEvent::kAny),
8ee2d316 185 doPlotGlobalRho(0), doVariableBinning(0), dovarbinTHnSparse(0),
daa728b0 186 makeQAhistos(0), makeBIAShistos(0), makeextraCORRhistos(0), makeoldJEThadhistos(0),
139e13f4 187 allpidAXIS(0), fcutType("EMCAL"), doPID(0), doPIDtrackBIAS(0),
56523939 188 doComments(0),
189 doFlavourJetAnalysis(0), fJetFlavTag(-99),
190 fBeam(kNA),
fc392675 191 fLocalRhoVal(0),
56523939 192 fTracksName(""), fTracksNameME(""), fJetsName(""),
fc392675 193 event(0),
aa5f3717 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
fc392675 197 fPoolMgr(0x0),
198 fPIDResponse(0x0), fTPCResponse(),
d60e724e 199 fESD(0), fAOD(0), fVevent(0),
200 fHistEventQA(0), fHistEventSelectionQA(0),
a409c5ba 201 fHistCentZvertGA(0), fHistCentZvertJE(0), fHistCentZvertMB(0), fHistCentZvertAny(0),
fc392675 202 fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
203 fHistRhovsCent(0), fHistNjetvsCent(0), fHistCentrality(0),
204 fHistZvtx(0), fHistMult(0),
205 fHistJetPhi(0), fHistTrackPhi(0),
d60e724e 206 fHistLocalRhoJetpt(0),
fc392675 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),
daa728b0 215 fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0), fhnCorr(0x0),
8ee2d316 216 fJetsCont(0), fTracksCont(0), fCaloClustersCont(0),
217 fContainerAllJets(0), fContainerPIDJets(1)
fc392675 218{
219 // Default Constructor
220 for(Int_t ilab=0; ilab<4; ilab++){
221 for(Int_t ipta=0; ipta<7; ipta++){
72a4d320 222 //fHistTrackEtaPhi[ilab][ipta]=0; //keep out for now
fc392675 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
fc392675 232
daa728b0 233 for(Int_t icent = 0; icent<6; ++icent){
fc392675 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){
daa728b0 245 fHistJetPt[i] = 0;
246 fHistJetPtBias[i] = 0;
247 fHistJetPtTT[i] = 0;
248 fHistAreavsRawPt[i] = 0;
fc392675 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;
8ee2d316 256 fHistJetPtcorrGlRho[i] = 0;
fc392675 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
fc392675 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
8ee2d316 296 // char array for naming histograms
297 int nchar = 200;
298 char name[nchar];
299
daa728b0 300 // strings for titles
301 TString name1;
302 TString title1;
303
fc392675 304 // create histograms that arn't array
daa728b0 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());
d60e724e 310 fHistLocalRhoJetpt = new TH1F("fHistLocalRhoJetpt","Local Rho corrected Jet p_{T}", 50, -50, 200);
311
a409c5ba 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
d60e724e 322 // Event QA histo
323 fHistEventQA = new TH1F("fHistEventQA", "Event Counter at checkpoints in code", 20, 0.5, 20.5);
f9b4eee4 324 SetfHistQAcounterLabels(fHistEventQA);
325 fOutput->Add(fHistEventQA);
d60e724e 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
daa728b0 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);
d60e724e 338 fOutput->Add(fHistLocalRhoJetpt);
daa728b0 339
340 // create histo's used for general QA
8ee2d316 341 if (makeQAhistos) {
a409c5ba 342 fHistTPCdEdX = new TH2F("TPCdEdX", "TPCdEdX", 2000, 0.0, 100.0, 500, 0, 500);
8ee2d316 343 fHistITSsignal = new TH2F("ITSsignal", "ITSsignal", 2000, 0.0, 100.0, 500, 0, 500);
daa728b0 344 // fHistTOFsignal = new TH2F("TOFsignal", "TOFsignal", 2000, 0.0, 100.0, 500, 0, 500);
8ee2d316 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);
daa728b0 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
8ee2d316 355
356 // add to output list
d60e724e 357 fOutput->Add(fHistTPCdEdX);
8ee2d316 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);
daa728b0 365 fOutput->Add(fHistTrackPtallcent);
366 fOutput->Add(fHistJetEtaPhi);
367 fOutput->Add(fHistJetHEtaPhi);
368 fOutput->Add(fHistSEphieta);
369 fOutput->Add(fHistMEphieta);
8ee2d316 370
72a4d320 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
8ee2d316 378
daa728b0 379 for (Int_t i = 0;i<6;++i){
8ee2d316 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]);
daa728b0 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
fc392675 444
daa728b0 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]);
fc392675 464
daa728b0 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]);
fc392675 469
8ee2d316 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]);
daa728b0 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]);
8ee2d316 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 }
3263f852 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};
e56c3f6e 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};
3263f852 581
e56c3f6e 582 // tracks: 51, jets: 26
3263f852 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;
8ee2d316 588
fc392675 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
daa728b0 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;
3badfe33 595 //if(fDoEventMixing) {
e56c3f6e 596 fhnJH = NewTHnSparseF("fhnJH", bitcoded);
3263f852 597
598 if(dovarbinTHnSparse){
599 fhnJH->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
600 fhnJH->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
601 }
602
603 fOutput->Add(fhnJH);
3badfe33 604 //}
fc392675 605
f9b4eee4 606 bitcodeCorr = 1<<0 | 1<<1 | 1<<2 | 1<<3; // | 1<<4 | 1<<5;
e56c3f6e 607 fhnCorr = NewTHnSparseFCorr("fhnCorr", bitcodeCorr);
3263f852 608 if(dovarbinTHnSparse) fhnCorr->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
3badfe33 609 fOutput->Add(fhnCorr);
f1ce4cc3 610
3badfe33 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
56523939 621 Int_t nCentralityBinspp = 8;
3badfe33 622 //Double_t centralityBinspp[nCentralityBinspp+1];
56523939 623 Double_t centralityBinspp[9] = {0.0, 4., 9, 15, 25, 35, 55, 100.0, 500.0};
3badfe33 624
625 // Setup for Pb-Pb collisions
a409c5ba 626 Int_t nCentralityBinsPbPb = 10; //100;
3badfe33 627 Double_t centralityBinsPbPb[nCentralityBinsPbPb+1];
628 for(Int_t ic=0; ic<nCentralityBinsPbPb; ic++){
a409c5ba 629 centralityBinsPbPb[ic]=10.0*ic; //1.0*ic;
fc392675 630 }
631
56523939 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);
8ee2d316 636// fOutput->Add(fHistMult);
fc392675 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;
56523939 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);
3badfe33 646// if(GetBeamType() == 0) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinspp, centralityBinspp, nZvtxBins, zvtxbin);
56523939 647// if(GetBeamType() == 1) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinsPbPb, centralityBinsPbPb, nZvtxBins, zvtxbin);
fc392675 648
649 // set up event mixing sparse
650 if(fDoEventMixing){
daa728b0 651 cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8; // | 1<<9;
e56c3f6e 652 fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
fc392675 653
3263f852 654 if(dovarbinTHnSparse){
655 fhnMixedEvents->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
656 fhnMixedEvents->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
657 }
fc392675 658
3263f852 659 fOutput->Add(fhnMixedEvents);
660 } // end of do-eventmixing
fc392675 661
662 // set up PID sparse
663 if(doPID){
8ee2d316 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");
f1ce4cc3 670 return;
8ee2d316 671 }
672
673 // PID response object
674 fPIDResponse = inputHandler->GetPIDResponse();
675 if (!fPIDResponse) {
676 AliError("PIDResponse object was not created");
f1ce4cc3 677 return;
8ee2d316 678 }
679 // *****************************************************************************************
680
681 // PID counter
d60e724e 682 fHistPID = new TH1F("fHistPID","PID Counter", 15, 0.5, 15.5);
8ee2d316 683 SetfHistPIDcounterLabels(fHistPID);
684 fOutput->Add(fHistPID);
685
139e13f4 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 |
daa728b0 688 1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14 | 1<<15 | 1<<16 | 1<<17 | 1<<18 | 1<<19 |
139e13f4 689 1<<20;
e56c3f6e 690 fhnPID = NewTHnSparseFPID("fhnPID", bitcode);
139e13f4 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;
e56c3f6e 694 fhnPID = NewTHnSparseFPID("fhnPID", bitcode);
139e13f4 695 }
8ee2d316 696
697 if(dovarbinTHnSparse){
139e13f4 698 fhnPID->GetAxis(1)->Set(nbinstrPT, xlowtrPT);
3263f852 699 fhnPID->GetAxis(8)->Set(nbinsjetPT, xlowjetPT);
8ee2d316 700 }
701
fc392675 702 fOutput->Add(fhnPID);
703 } // end of do-PID
704
fc392675 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();
8ee2d316 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;
fc392675 737}
738
739//_________________________________________________________________________
740Bool_t AliAnalysisTaskEmcalJetHadEPpid::Run()
741{ // Main loop called for each event
742 // TEST TEST TEST TEST for OBJECTS!
3badfe33 743
f9b4eee4 744 fHistEventQA->Fill(1); // All Events that get entered
745
d60e724e 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
fc392675 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
f9b4eee4 786 fHistEventQA->Fill(2); // events after object check
787
fc392675 788 // what kind of event do we have: AOD or ESD?
139e13f4 789 Bool_t useAOD;
790 if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
791 else useAOD = kFALSE;
792
fc392675 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
f9b4eee4 811 fHistEventQA->Fill(3); // events after Aod/esd check
812
fc392675 813 // get centrality
814 Int_t centbin = GetCentBin(fCent);
8ee2d316 815 if (makeQAhistos) fHistCentrality->Fill(fCent); // won't be filled in pp collision (Keep this in mind!)
fc392675 816
3badfe33 817 // BEAM TYPE enumerator: kNA = -1, kpp = 0, kAA = 1, kpA = 2
fc392675 818 // for pp analyses we will just use the first centrality bin
3badfe33 819 if(GetBeamType() == 0) if (centbin == -1) centbin = 0;
820 if(GetBeamType() == 1) if (centbin == -1) return kTRUE;
fc392675 821
3badfe33 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 }
8ee2d316 827
f9b4eee4 828 fHistEventQA->Fill(4); // events after centrality check
829
fc392675 830 // get vertex information
831 Double_t fvertex[3]={0,0,0};
832 InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
833 Double_t zVtx=fvertex[2];
8ee2d316 834 if (makeQAhistos) fHistZvtx->Fill(zVtx);
fc392675 835
836 // get z-vertex bin
837 //Int_t zVbin = GetzVertexBin(zVtx);
838
839 // apply zVtx cut
8ee2d316 840 if(fabs(zVtx)>10.0) return kTRUE;
fc392675 841
f9b4eee4 842 fHistEventQA->Fill(5); // events after zvertex check
843
fc392675 844 // create pointer to list of input event
845 TList *list = InputEvent()->GetList();
846 if(!list) {
847 AliError(Form("ERROR: list not attached\n"));
8ee2d316 848 return kTRUE;
fc392675 849 }
850
f9b4eee4 851 fHistEventQA->Fill(6); // events after list check
852
fc392675 853 // initialize TClonesArray pointers to jets and tracks
854 TClonesArray *jets = 0;
855 TClonesArray *tracks = 0;
3badfe33 856 TClonesArray *tracksME = 0;
fc392675 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
3badfe33 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
fc392675 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
f9b4eee4 879 fHistEventQA->Fill(7); // events after track/jet pointer check
880
fc392675 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
f9b4eee4 887 fHistEventQA->Fill(8); // events after #track and jets < 1 check
888
8ee2d316 889 if (makeQAhistos) fHistMult->Fill(Ntracks); // fill multiplicity distribution
fc392675 890
891 // initialize track parameters
892 Int_t iTT=-1;
893 Double_t ptmax=-10;
a409c5ba 894 Int_t NtrackAcc = 0;
fc392675 895
3badfe33 896 fVevent = dynamic_cast<AliVEvent*>(InputEvent());
897 if (!fVevent) {
898 printf("ERROR: fVEvent not available\n");
899 return kTRUE;
900 }
901
a409c5ba 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);
3badfe33 907
fc392675 908 // loop over tracks - to get hardest track (highest pt)
909 for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++){
3badfe33 910 AliVParticle* Vtrack = static_cast<AliVParticle*>(tracks->At(iTracks));
911 if (!Vtrack) {
912 printf("ERROR: Could not receive track %d\n", iTracks);
e56c3f6e 913 continue;
3263f852 914 }
3badfe33 915
916 AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
917 if(!track) continue;
3263f852 918
daa728b0 919 // track cuts
fc392675 920 if(TMath::Abs(track->Eta())>0.9) continue;
921 if(track->Pt()<0.15) continue;
922 //iCount++;
56523939 923 NtrackAcc++;
924
fc392675 925 if(track->Pt()>ptmax){
926 ptmax=track->Pt(); // max pT track
927 iTT=iTracks; // trigger tracks
928 } // check if Pt>maxpt
929
8ee2d316 930 if (makeQAhistos) fHistTrackPhi->Fill(track->Phi());
931 if (makeQAhistos) fHistTrackPt[centbin]->Fill(track->Pt());
daa728b0 932 if (makeQAhistos) fHistTrackPtallcent->Fill(track->Pt());
fc392675 933 } // end of loop over tracks
934
fc392675 935 // get rho from event and fill relative histo's
936 fRho = GetRhoFromEvent(fRhoName);
937 fRhoVal = fRho->GetVal();
fc392675 938
8ee2d316 939 if (makeQAhistos) {
daa728b0 940 fHistRhovsdEP[centbin]->Fill(fRhoVal,fEPV0); // Global Rho vs delta Event Plane angle
8ee2d316 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
fc392675 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;
daa728b0 953 Double_t leadhadronPT = 0;
fc392675 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
f1ce4cc3 961 // apply jet cuts
8ee2d316 962 if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) continue;
963 if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) continue;
daa728b0 964 if (makeQAhistos) fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
8ee2d316 965 if(!AcceptMyJet(jet)) continue;
fc392675 966
967 NjetAcc++; // # of accepted jets
56523939 968
969 // if FlavourJetAnalysis, get desired FlavTag and check against Jet
970 if(doFlavourJetAnalysis) { if(!AcceptFlavourJet(jet, fJetFlavTag)) continue;}
971
8ee2d316 972 // use this to get total # of jets passing cuts in events!!!!!!!!
973 if (makeQAhistos) fHistJetPhi->Fill(jet->Phi()); // Jet Phi histogram (filled)
fc392675 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;
f9b4eee4 985 fHistEventQA->Fill(9); // events after track/jet loop to get highest pt
986
fc392675 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));
8ee2d316 990 if (!jet) continue;
991
d60e724e 992 if (!(trig & fTriggerEventType)) continue;
993
8ee2d316 994 // (should probably be higher..., but makes a cut on jet pT)
995 if (jet->Pt()<0.1) continue;
f1ce4cc3 996 // do we accept jet? apply jet cuts
8ee2d316 997 if (!AcceptMyJet(jet)) continue;
fc392675 998
56523939 999 // if FlavourJetAnalysis, get desired FlavTag and check against Jet
1000 if(doFlavourJetAnalysis) { if(!AcceptFlavourJet(jet, fJetFlavTag)) continue;}
1001
f9b4eee4 1002 fHistEventQA->Fill(10); // accepted jets
1003
fc392675 1004 // check on lead jet
1005 Double_t leadjet=0;
1006 if (ijet==ijethi) leadjet=1;
1007
daa728b0 1008 // check on leading hadron pt
139e13f4 1009 if (ijet==ijethi) leadhadronPT = GetLeadingHadronPt(jet);
daa728b0 1010
fc392675 1011 // initialize and calculate various parameters: pt, eta, phi, rho, etc...
1012 Double_t jetphi = jet->Phi(); // phi of jet
1013 NJETAcc++; // # accepted jets
fc392675 1014 Double_t jeteta = jet->Eta(); // ETA of jet
aa5f3717 1015 Double_t jetPt = -500;
1016 Double_t jetPtGlobal = -500;
d60e724e 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 }
8ee2d316 1022 jetPt = jet->Pt();
1023 jetPtGlobal = jet->Pt()-jet->Area()*fRhoVal; // corrected pT of jet from rho value
a8150e70 1024 Double_t dEP = -500; // initialize angle between jet and event plane
fc392675 1025 dEP = RelativeEPJET(jetphi,fEPV0); // angle betweeen jet and event plane
1026
1027 // make histo's
8ee2d316 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);
daa728b0 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
8ee2d316 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());
daa728b0 1038 if (makeoldJEThadhistos) fHistJetPt[centbin]->Fill(jet->Pt()); // fill #jets vs pT histo
8ee2d316 1039 //fHistDeltaPtvsArea->Fill(jetPt,jet->Area());
fc392675 1040
1041 // make histo's with BIAS applied
1042 if (jet->MaxTrackPt()>fTrkBias){
8ee2d316 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());
fc392675 1049 }
1050
8ee2d316 1051 //if(leadjet && centbin==0){
1052 // if(makeextraCORRhistos) fHistJetPt[centbin+1]->Fill(jet->Pt());
1053 //}
fc392675 1054 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
daa728b0 1055 if (makeoldJEThadhistos){
8ee2d316 1056 fHistJetPtBias[centbin]->Fill(jet->Pt());
1057 //if(leadjet && centbin==0) fHistJetPtBias[centbin+1]->Fill(jet->Pt());
1058 }
fc392675 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
8ee2d316 1068 if(passedTTcut) {
daa728b0 1069 if (makeoldJEThadhistos) fHistJetPtTT[centbin]->Fill(jet->Pt());
8ee2d316 1070 }
fc392675 1071
1072 // cut on HIGHEST jet pt in event (15 GeV default)
139e13f4 1073 //if (highestjetpt>fJetPtcut) {
1074 if (jet->Pt() > fJetPtcut) {
f9b4eee4 1075 fHistEventQA->Fill(11); // jets meeting pt threshold
72a4d320 1076
daa728b0 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
f9b4eee4 1080 Double_t CorrEntries[4] = {fCent, jet->Pt(), dEP, zVtx};
daa728b0 1081 fhnCorr->Fill(CorrEntries); // fill Sparse Histo with Correction entries
d60e724e 1082
1083 if(GetBeamType() == 1) fHistLocalRhoJetpt->Fill(jetPtLocal);
daa728b0 1084 }
1085
fc392675 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) {
3badfe33 1088 AliVParticle* Vtrackass = static_cast<AliVParticle*>(tracks->At(iTracks));
1089 if (!Vtrackass) {
1090 printf("ERROR: Could not receive track %d\n", iTracks);
1091 continue;
3263f852 1092 }
1093
3badfe33 1094 AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrackass);
fc392675 1095 if (!track) {
8ee2d316 1096 AliError(Form("Couldn't get AliVtrack %d\n", iTracks));
fc392675 1097 continue;
1098 }
1099
8ee2d316 1100 // apply track cuts
daa728b0 1101 if(TMath::Abs(track->Eta())>fTrkEta) continue;
1102 if (track->Pt()<0.15) continue;
fc392675 1103
f9b4eee4 1104 fHistEventQA->Fill(12); // accepted tracks in events from trigger jets
1105
fc392675 1106 // calculate and get some track parameters
a8150e70 1107 Double_t trCharge = -99;
1108 trCharge = track->Charge();
fc392675 1109 Double_t tracketa=track->Eta(); // eta of track
1110 Double_t deta=tracketa-jeteta; // dETA between track and jet
daa728b0 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
f1ce4cc3 1117 if(ieta<0) continue; // double check we don't have a negative array index
daa728b0 1118 iptjet=GetpTjetBin(jet->Pt()); // bin of jet pt
1119 if(iptjet<0) continue; // double check we don't have a negative array index
f1ce4cc3 1120 }
fc392675 1121
1122 // dPHI between jet and hadron
1123 Double_t dphijh = RelativePhi(jet->Phi(), track->Phi()); // angle between jet and hadron
1124
fc392675 1125 // fill some jet-hadron histo's
daa728b0 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
8ee2d316 1128 fHistJetHaddPHI->Fill(dphijh);
1129 if(passedTTcut){
daa728b0 1130 if (makeoldJEThadhistos) fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
8ee2d316 1131 }
fc392675 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
daa728b0 1136 Double_t triggerEntries[9] = {fCent, jet->Pt(), track->Pt(), deta, dphijh, dEP, zVtx, trCharge, leadjet};
72a4d320 1137 fhnJH->Fill(triggerEntries); // fill Sparse Histo with trigger entries
d60e724e 1138
fc392675 1139 // fill histo's
daa728b0 1140 if(makeQAhistos) fHistSEphieta->Fill(dphijh, deta); // single event distribution
1141 if (makeoldJEThadhistos) fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
8ee2d316 1142
1143 if (makeBIAShistos) {
1144 fHistJetHaddPhiBias->Fill(dphijh);
fc392675 1145
8ee2d316 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
fc392675 1158 } // end of check maxtrackpt>ftrackbias or maxclusterpt>fclustbias
1159
daa728b0 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 }
8ee2d316 1167
daa728b0 1168 // some variables for PID
daa728b0 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
daa728b0 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
f9b4eee4 1189 fHistEventQA->Fill(13); // check for AliVEvent and fPIDresponse objects
1190
e56c3f6e 1191///////////////////////////////////////
daa728b0 1192
3badfe33 1193 // get detector signals
1194 dEdx = track->GetTPCsignal();
1195 ITSsig = track->GetITSsignal();
1196 TOFsig = track->GetTOFsignal();
daa728b0 1197
3badfe33 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);
e56c3f6e 1212
daa728b0 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;
d60e724e 1227 nPID = 1;
daa728b0 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;
d60e724e 1233 nPID=2;
daa728b0 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;
d60e724e 1240 nPID=3;
daa728b0 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;
d60e724e 1247 nPID=4;
daa728b0 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;
d60e724e 1254 nPID = 5;
daa728b0 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;
d60e724e 1260 nPID=6;
daa728b0 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;
d60e724e 1267 nPID=7;
daa728b0 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;
d60e724e 1274 nPID=8;
daa728b0 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;
d60e724e 1281 nPID = 9;
daa728b0 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;
d60e724e 1287 nPID=10;
daa728b0 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;
d60e724e 1294 nPID=11;
daa728b0 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;
d60e724e 1301 nPID=12;
daa728b0 1302 }else isPtof = kFALSE;
1303 } // cut on track pT for TOF
1304
d60e724e 1305 if (nPID == -99) nPID = 14;
daa728b0 1306 fHistPID->Fill(nPID);
fc392675 1307
1308 // PID sparse getting filled
139e13f4 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
fc392675 1323 }; //array for PID sparse
139e13f4 1324 fhnPID->Fill(pid_Entries); // fill Sparse histo of PID tracks
1325 } // minimal pid sparse filling
1326
daa728b0 1327 } // end of doPID check
fc392675 1328
8ee2d316 1329 // get track pt bin
fc392675 1330 Int_t itrackpt = -500; // initialize track pT bin
1331 itrackpt = GetpTtrackBin(track->Pt());
1332
8ee2d316 1333 // all tracks: jet hadron correlations in hadron pt bins
daa728b0 1334 if(makeextraCORRhistos) fHistJetHadbindPhi[itrackpt]->Fill(dphijh);
fc392675 1335
1336 // in plane and out of plane jet-hadron histo's
1337 if( dEP>0 && dEP<=(TMath::Pi()/6) ){
daa728b0 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);
fc392675 1343 }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
daa728b0 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);
fc392675 1349 }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
daa728b0 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);
fc392675 1354 }
1355 } // loop over tracks found in event with highest JET pT > 10.0 GeV (change)
1356 } // jet pt cut
1357 } // jet loop
1358
f9b4eee4 1359 fHistEventQA->Fill(14); // events right before event mixing
1360
fc392675 1361// ***************************************************************************************************************
1362// ******************************** Event MIXING *****************************************************************
3badfe33 1363// ***************************************************************************************************************
1364
1365 // initialize object array of cloned picotracks
1366 TObjArray* tracksClone = 0x0;
1367
1368 // PbPb collisions - create cloned picotracks
d60e724e 1369 //if(GetBeamType() == 1) tracksClone = CloneAndReduceTrackList(tracks); // TEST
fc392675 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
3badfe33 1392 // get the trigger bit, need to change trigger bits between different runs
1393 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
d60e724e 1394 // if event was not selected (triggered) for any reseason (should never happen) then return
1395 if (trigger==0) return kTRUE;
fc392675 1396
3badfe33 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);
3badfe33 1406 poolpp = fPoolMgr->GetEventPool(Ntrks, zVtx); // for pp
1407 }
fc392675 1408
3badfe33 1409 // PbPb collisions - get event pool
1410 if(GetBeamType() == 1) pool = fPoolMgr->GetEventPool(fCent, zVtx); // for PbPb? fcent
fc392675 1411
d60e724e 1412 // if we don't have a pool, return
3badfe33 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));
fc392675 1416 return kTRUE;
1417 }
1418
f9b4eee4 1419 fHistEventQA->Fill(15); // mixed events cases that have pool
1420
3badfe33 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 ??
a409c5ba 1428 if(GetBeamType() == 1) if(trigger & fTriggerEventType) { //kEMCEJE)) {
1429 if (pool->IsReady() || pool->NTracksInPool() > fNMIXtracks || pool->GetCurrentNEvents() >= fNMIXevents) {
3badfe33 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
d60e724e 1475 fhnMixedEvents->Fill(triggerEntries,1./nMix); // fill Sparse histo of mixed events
3badfe33 1476
d60e724e 1477 fHistEventQA->Fill(18); // event mixing - nbgtracks
3badfe33 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
d60e724e 1483 } // end of check for pool being ready
1484 } // end EMC triggered loop
3badfe33 1485
1486//=============================================================================================================
1487
fc392675 1488 // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1489/// if (trigger & AliVEvent::kEMC1) {
3badfe33 1490 // pp collisions
d60e724e 1491 if(GetBeamType() == 0) if(trigger & fTriggerEventType) { //kEMC1)) {
a409c5ba 1492 if (poolpp->IsReady() || poolpp->NTracksInPool() > fNMIXtracks || poolpp->GetCurrentNEvents() >= fNMIXevents) {
3badfe33 1493
fc392675 1494 // loop over jets (passing cuts?)
1495 for (Int_t ijet = 0; ijet < Njets; ijet++) {
8ee2d316 1496 Double_t leadjet=0;
1497 if (ijet==ijethi) leadjet=1;
fc392675 1498
1499 // get jet object
1500 AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
8ee2d316 1501 if (!jet) continue;
fc392675 1502
8ee2d316 1503 // (should probably be higher..., but makes a cut on jet pT)
1504 if (jet->Pt()<0.1) continue;
1505 if (!AcceptMyJet(jet)) continue;
fc392675 1506
f9b4eee4 1507 fHistEventQA->Fill(16); // event mixing jets
1508
3badfe33 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;
fc392675 1511
3badfe33 1512 // get number of current events in pool
1513 Int_t nMix = poolpp->GetCurrentNEvents(); // how many particles in pool to mix
e56c3f6e 1514
fc392675 1515 // Fill for biased jet triggers only
e56c3f6e 1516 if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) { // && jet->Pt() > fJetPtcut) {
fc392675 1517 // Fill mixed-event histos here
1518 for (Int_t jMix=0; jMix<nMix; jMix++) {
3badfe33 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
d60e724e 1535
3badfe33 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
d60e724e 1539
1540 fHistEventQA->Fill(18); // event mixing - nbgtracks
3badfe33 1541 if(makeextraCORRhistos) fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1542 } // end of background track loop
1543 } // end of filling mixed-event histo's
fc392675 1544 } // end of check for biased jet triggers
1545 } // end of jet loop
d60e724e 1546 } // end of check for pool being ready
3badfe33 1547 } //end EMC triggered loop
fc392675 1548
3badfe33 1549 // pp collisions
fc392675 1550 // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1551/// if (trigger & AliVEvent::kMB) {
8ee2d316 1552//T if (trigger & AliVEvent::kAnyINT){ // test
3badfe33 1553 if(GetBeamType() == 0) {
1554
d60e724e 1555 // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1556 if(trigger & fMixingEventType) { //kMB) {
3badfe33 1557
1558 // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
d60e724e 1559 tracksClone = CloneAndReduceTrackList(tracks);
1560
3badfe33 1561 // update pool if jet in event or not
1562 poolpp->UpdatePool(tracksClone);
1563
1564 } // check on track from MB events
1565 }
fc392675 1566
3badfe33 1567 // PbPb collisions
1568 if(GetBeamType() == 1) {
d60e724e 1569
1570 // use only tracks from MB events
1571 if(trigger & fMixingEventType) { //kMB) {
3badfe33 1572
d60e724e 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
fc392675 1581 } // end of event mixing
1582
1583 // print some stats on the event
1584 event++;
f9b4eee4 1585 fHistEventQA->Fill(19); // events making it to end
1586
f1ce4cc3 1587 if (doComments) {
1588 cout<<"Event #: "<<event<<" Jet Radius: "<<fJetRad<<" Constituent Pt Cut: "<<fConstituentCut<<endl;
56523939 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;
f1ce4cc3 1591 cout<<" =============================================== "<<endl;
1592 }
fc392675 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//________________________________________________________________________
a8150e70 1612Double_t AliAnalysisTaskEmcalJetHadEPpid::RelativePhi(Double_t mphi,Double_t vphi) const
fc392675 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
fc392675 1627//_________________________________________________________________________
a8150e70 1628Double_t AliAnalysisTaskEmcalJetHadEPpid:: RelativeEPJET(Double_t jetAng, Double_t EPAng) const
fc392675 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;
f1ce4cc3 1660 if (TMath::Abs(eta)<=0.4) etabin = 0;
fc392675 1661 else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8) etabin = 1;
f1ce4cc3 1662 else if (TMath::Abs(eta)>=0.8) etabin = 2;
fc392675 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;
f1ce4cc3 1676 else if (pt>=60) ptbin = 4;
fc392675 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
fc392675 1701//___________________________________________________________________________
1702Int_t AliAnalysisTaskEmcalJetHadEPpid::GetzVertexBin(Double_t zVtx) const
1703{
1704 // get z-vertex bin for histo.
1705 int zVbin= -1;
f1ce4cc3 1706 if (zVtx>=-10 && zVtx<-8) zVbin = 0;
fc392675 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//______________________________________________________________________
e56c3f6e 1722THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseF(const char* name, UInt_t entries)
fc392675 1723{
e56c3f6e 1724 // generate new THnSparseF, axes are defined in GetDimParams()
fc392675 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
e56c3f6e 1753 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1754} // end of NewTHnSparseF
fc392675 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:
a8150e70 1771 label = "Jet p_{T}";
e56c3f6e 1772 nbins = 50; // 216
a8150e70 1773 xmin = 0.;
e56c3f6e 1774 xmax = 250.; // 216
fc392675 1775 break;
1776
1777 case 2:
1778 label = "Track p_{T}";
e56c3f6e 1779 nbins = 80; //300; // 750 pid
fc392675 1780 xmin = 0.;
e56c3f6e 1781 xmax = 20.; //75.;
fc392675 1782 break;
1783
1784 case 3:
1785 label = "Relative Eta";
daa728b0 1786 nbins = 48;
fc392675 1787 xmin = -1.6;
1788 xmax = 1.6;
1789 break;
1790
1791 case 4:
1792 label = "Relative Phi";
daa728b0 1793 nbins = 72;
fc392675 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";
e56c3f6e 1800 nbins = 12; // 72
1801 xmin = 0;
1802 xmax = 0.5*pi;
fc392675 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";
a8150e70 1814 nbins = 3;
fc392675 1815 xmin = -1.5;
1816 xmax = 1.5;
1817 break;
1818
daa728b0 1819 case 8:
fc392675 1820 label = "leading jet";
1821 nbins = 3;
1822 xmin = -0.5;
1823 xmax = 2.5;
1824 break;
a8150e70 1825
1826 case 9: // need to update
1827 label = "leading track";
1828 nbins = 10;
1829 xmin = 0;
1830 xmax = 50;
1831 break;
1832
fc392675 1833 } // end of switch
1834} // end of getting dim-params
1835
1836//_________________________________________________
1837// From CF event mixing code PhiCorrelations
3badfe33 1838TObjArray* AliAnalysisTaskEmcalJetHadEPpid::CloneAndReduceTrackList(TObjArray* tracksME)
fc392675 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
3263f852 1844// ===============================
1845
1846// cout << "RM Hybrid track : " << i << " " << particle->Pt() << endl;
fc392675 1847
1848 //cout << "nEntries " << tracks->GetEntriesFast() <<endl;
3badfe33 1849 for (Int_t i=0; i<tracksME->GetEntriesFast(); i++) { // AOD/general case
1850 AliVParticle* particle = (AliVParticle*) tracksME->At(i); // AOD/general case
fc392675 1851 if(TMath::Abs(particle->Eta())>fTrkEta) continue;
1852 if(particle->Pt()<0.15)continue;
1853
8ee2d316 1854/*
1855// DON'T USE
fc392675 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;
8ee2d316 1870// end of DON'T USE
fc392675 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
8ee2d316 1876 if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi()); // TEST
1877*/
fc392675 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//____________________________________________________________________________________________
e56c3f6e 1886THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFPID(const char* name, UInt_t entries)
fc392675 1887{
e56c3f6e 1888 // generate new THnSparseF PID, axes are defined in GetDimParams()
fc392675 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
e56c3f6e 1917 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1918} // end of NewTHnSparseF PID
fc392675 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:
fc392675 1936 label = "Track p_{T}";
e56c3f6e 1937 nbins = 80; //300; // 750
fc392675 1938 xmin = 0.;
e56c3f6e 1939 xmax = 20.; //75.;
fc392675 1940 break;
1941
139e13f4 1942 case 2:
fc392675 1943 label = "Charge of Track";
1944 nbins = 3;
1945 xmin = -1.5;
1946 xmax = 1.5;
1947 break;
1948
139e13f4 1949 case 3:
fc392675 1950 label = "Relative Eta of Track and Jet";
72a4d320 1951 nbins = 48;
fc392675 1952 xmin = -1.6;
1953 xmax = 1.6;
1954 break;
1955
139e13f4 1956 case 4:
fc392675 1957 label = "Relative Phi of Track and Jet";
daa728b0 1958 nbins = 72;
fc392675 1959 xmin = -0.5*pi;
1960 xmax = 1.5*pi;
1961 break;
1962
139e13f4 1963 case 5:
fc392675 1964 label = "leading jet";
1965 nbins = 3;
1966 xmin = -.5;
1967 xmax = 2.5;
1968 break;
1969
139e13f4 1970 case 6:
fc392675 1971 label = "Z-vertex";
1972 nbins = 10;
1973 xmin = -10.;
1974 xmax = 10.;
1975 break;
1976
139e13f4 1977 case 7:
8ee2d316 1978 label = "Relative angle: Jet and Reaction Plane";
e56c3f6e 1979 nbins = 12; // 48
fc392675 1980 xmin = 0.;
8ee2d316 1981 xmax = 0.5*pi;
fc392675 1982 break;
1983
139e13f4 1984 case 8:
1985 label = "Jet p_{T}";
e56c3f6e 1986 nbins = 50; // 216
fc392675 1987 xmin = 0.;
e56c3f6e 1988 xmax = 250.; // 216
fc392675 1989 break;
1990
139e13f4 1991 case 9:
fc392675 1992 label = "N-Sigma of pions in TPC";
8ee2d316 1993 nbins = 200;
daa728b0 1994 xmin = -10.0;
1995 xmax = 10.0;
fc392675 1996 break;
1997
139e13f4 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
fc392675 2012 case 12:
139e13f4 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:
fc392675 2027 label = "N-Sigma of protons in TPC";
daa728b0 2028 nbins = 200;
2029 xmin = -10.;
2030 xmax = 10.;
fc392675 2031 break;
2032
139e13f4 2033 case 15:
fc392675 2034 label = "N-Sigma of kaons in TPC";
daa728b0 2035 nbins = 200;
2036 xmin = -10.;
2037 xmax = 10.;
fc392675 2038 break;
2039
139e13f4 2040 case 16:
fc392675 2041 label = "N-Sigma of pions in ITS";
8ee2d316 2042 nbins = 200;
daa728b0 2043 xmin = -10.0;
2044 xmax = 10.0;
fc392675 2045 break;
2046
139e13f4 2047 case 17:
fc392675 2048 label = "N-Sigma of protons in ITS";
daa728b0 2049 nbins = 200;
2050 xmin = -10.;
2051 xmax = 10.;
fc392675 2052 break;
2053
139e13f4 2054 case 18:
fc392675 2055 label = "N-Sigma of kaons in ITS";
daa728b0 2056 nbins = 200;
2057 xmin = -10.;
2058 xmax = 10.;
fc392675 2059 break;
2060
139e13f4 2061 case 19:
fc392675 2062 label = "N-Sigma of protons in TOF";
daa728b0 2063 nbins = 200;
2064 xmin = -10.;
2065 xmax = 10.;
fc392675 2066 break;
2067
139e13f4 2068 case 20:
fc392675 2069 label = "N-Sigma of kaons in TOF";
daa728b0 2070 nbins = 200;
2071 xmin = -10.;
2072 xmax = 10.;
fc392675 2073 break;
2074
fc392675 2075 } // end of switch
2076} // end of get dimension parameters PID
2077
8ee2d316 2078void AliAnalysisTaskEmcalJetHadEPpid::Terminate(Option_t *) {
fc392675 2079 cout<<"#########################"<<endl;
2080 cout<<"#### DONE RUNNING!!! ####"<<endl;
2081 cout<<"#########################"<<endl;
fc392675 2082} // end of terminate
2083
8ee2d316 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;
8ee2d316 2094
2095 //passed all above cuts
2096 return 1;
2097}
fc392675 2098
2099//void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
2100void AliAnalysisTaskEmcalJetHadEPpid::SetfHistPIDcounterLabels(TH1* h) const
2101{
2102 // fill the analysis summary histrogram, saves all relevant analysis settigns
d60e724e 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");
fc392675 2119}
daa728b0 2120
f9b4eee4 2121//void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
2122void AliAnalysisTaskEmcalJetHadEPpid::SetfHistQAcounterLabels(TH1* h) const
2123{
d60e724e 2124 // label bins of the analysis event summary
f9b4eee4 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");
e56c3f6e 2132 h->GetXaxis()->SetBinLabel(8, "tracks & jets < than 1 check");
f9b4eee4 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");
d60e724e 2136 h->GetXaxis()->SetBinLabel(12, "accepted tracks in events w/ trigger jet");
2137 h->GetXaxis()->SetBinLabel(13, "after AliVEvent & fPIDResponse");
f9b4eee4 2138 h->GetXaxis()->SetBinLabel(14, "events before event mixing");
d60e724e 2139 h->GetXaxis()->SetBinLabel(15, "mixed events w/ pool");
f9b4eee4 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");
d60e724e 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");
f9b4eee4 2176}
2177
daa728b0 2178//______________________________________________________________________
e56c3f6e 2179THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFCorr(const char* name, UInt_t entries) {
daa728b0 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
e56c3f6e 2209 return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
2210} // end of NewTHnSparseF
daa728b0 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}";
e56c3f6e 2229 nbins = 50; // 216
daa728b0 2230 xmin = 0.;
e56c3f6e 2231 xmax = 250.; // 216
daa728b0 2232 break;
2233
2234 case 2:
2235 label = "Relative angle: Jet and Reaction Plane";
e56c3f6e 2236 nbins = 12; // 48
daa728b0 2237 xmin = 0.;
2238 xmax = 0.5*pi;
2239 break;
2240
f9b4eee4 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";
e56c3f6e 2250 nbins = 50; // 250
f9b4eee4 2251 xmin = -50.;
2252 xmax = 200.;
2253 break;
2254
2255 case 5:
2256 label = "Jet p_{T} corrected with Global Rho";
e56c3f6e 2257 nbins = 50; // 250
f9b4eee4 2258 xmin = -50.;
2259 xmax = 200.;
2260 break;
2261
daa728b0 2262 }// end of switch
2263} // end of Correction (ME) sparse
2264
56523939 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
fdf5fcfe 2275 Int_t flavNUM = -99;//, flavNUM2 = -99, flavNUM3 = -99; FIXME commented out to avoid compiler warning
56523939 2276 flavNUM = NUM;
2277 //flavNUM2 = NUM2;
2278 //flavNUM3 = NUM3;
daa728b0 2279
56523939 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}