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