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