]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
cosmetics
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtAnalysisPbPbAOD.cxx
CommitLineData
d25bcbe6 1/**************************************************************************
2* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15//------------------------------------------------------------------------------
16// AlidNdPtAnalysisPbPbAOD class.
17//
18// Author: P. Luettig, 15.05.2013
19//------------------------------------------------------------------------------
20
21
22#include "AlidNdPtAnalysisPbPbAOD.h"
23
24
25using namespace std;
26
27ClassImp(AlidNdPtAnalysisPbPbAOD)
28
29// dummy constructor
30AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD() : AliAnalysisTaskSE(),
31fOutputList(0),
32// Histograms
33hPt(0),
34hMCPt(0),
35hnZvPtEtaCent(0),
36hnMCRecPrimZvPtEtaCent(0),
37hnMCGenZvPtEtaCent(0),
5747f2c6 38hnMCRecSecZvPtEtaCent(0),
d25bcbe6 39hEventStatistics(0),
40hEventStatisticsCentrality(0),
41hAllEventStatisticsCentrality(0),
42hnZvMultCent(0),
43hTriggerStatistics(0),
44hMCTrackPdgCode(0),
45hMCTrackStatusCode(0),
46hCharge(0),
47hMCCharge(0),
48hMCPdgPt(0),
49hMCHijingPrim(0),
50hAccNclsTPC(0),
51hAccCrossedRowsTPC(0),
5747f2c6 52hDCAPtAll(0),
53hDCAPtAccepted(0),
54hMCDCAPtSecondary(0),
55hMCDCAPtPrimary(0),
d25bcbe6 56//global
57bIsMonteCarlo(0),
58// event cut variables
59dCutMaxZVertex(10.),
60// track kinematic cut variables
61dCutPtMin(0.15),
62dCutPtMax(1000.),
63dCutEtaMin(-0.8),
64dCutEtaMax(0.8),
65// track quality cut variables
66bCutRequireTPCRefit(kTRUE),
67dCutMinNumberOfCrossedRows(120.),
68dCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
69dCutMaxChi2PerClusterTPC(4.),
70dCutMaxFractionSharedTPCClusters(0.4),
71dCutMaxDCAToVertexZ(3.0),
72dCutMaxDCAToVertexXY(3.0),
73bCutRequireITSRefit(kTRUE),
74dCutMaxChi2PerClusterITS(36.),
75dCutDCAToVertex2D(kFALSE),
76dCutRequireSigmaToVertex(kFALSE),
77dCutMaxDCAToVertexXYPtDepPar0(0.0182),
78dCutMaxDCAToVertexXYPtDepPar1(0.0350),
79dCutMaxDCAToVertexXYPtDepPar2(1.01),
80bCutAcceptKinkDaughters(kFALSE),
81dCutMaxChi2TPCConstrainedGlobal(36.),
82// binning for THnSparse
83fMultNbins(0),
84fPtNbins(0),
85fPtCorrNbins(0),
86fEtaNbins(0),
87fZvNbins(0),
88fCentralityNbins(0),
89fBinsMult(0),
90fBinsPt(0),
91fBinsPtCorr(0),
92fBinsEta(0),
93fBinsZv(0),
94fBinsCentrality(0)
95{
96
97 fMultNbins = 0;
98 fPtNbins = 0;
99 fPtCorrNbins = 0;
100 fEtaNbins = 0;
101 fZvNbins = 0;
102 fCentralityNbins = 0;
103 fBinsMult = 0;
104 fBinsPt = 0;
105 fBinsPtCorr = 0;
106 fBinsEta = 0;
107 fBinsEta = 0;
108 fBinsZv = 0;
109 fBinsCentrality = 0;
110
111}
112
113AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD(const char *name) : AliAnalysisTaskSE(name),
114fOutputList(0),
115// Histograms
116hPt(0),
117hMCPt(0),
118hnZvPtEtaCent(0),
119hnMCRecPrimZvPtEtaCent(0),
120hnMCGenZvPtEtaCent(0),
5747f2c6 121hnMCRecSecZvPtEtaCent(0),
d25bcbe6 122hEventStatistics(0),
123hEventStatisticsCentrality(0),
124hAllEventStatisticsCentrality(0),
125hnZvMultCent(0),
126hTriggerStatistics(0),
127hMCTrackPdgCode(0),
128hMCTrackStatusCode(0),
129hCharge(0),
130hMCCharge(0),
131hMCPdgPt(0),
132hMCHijingPrim(0),
133hAccNclsTPC(0),
134hAccCrossedRowsTPC(0),
5747f2c6 135hDCAPtAll(0),
136hDCAPtAccepted(0),
137hMCDCAPtSecondary(0),
138hMCDCAPtPrimary(0),
d25bcbe6 139//global
140bIsMonteCarlo(0),
141// event cut variables
142dCutMaxZVertex(10.),
143// track kinematic cut variables
144dCutPtMin(0.15),
145dCutPtMax(200.),
146dCutEtaMin(-0.8),
147dCutEtaMax(0.8),
148// track quality cut variables
149bCutRequireTPCRefit(kTRUE),
150dCutMinNumberOfCrossedRows(120.),
151dCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
152dCutMaxChi2PerClusterTPC(4.),
153dCutMaxFractionSharedTPCClusters(0.4),
154dCutMaxDCAToVertexZ(3.0),
155dCutMaxDCAToVertexXY(3.0),
156bCutRequireITSRefit(kTRUE),
157dCutMaxChi2PerClusterITS(36.),
158dCutDCAToVertex2D(kFALSE),
159dCutRequireSigmaToVertex(kFALSE),
160dCutMaxDCAToVertexXYPtDepPar0(0.0182),
161dCutMaxDCAToVertexXYPtDepPar1(0.0350),
162dCutMaxDCAToVertexXYPtDepPar2(1.01),
163bCutAcceptKinkDaughters(kFALSE),
164dCutMaxChi2TPCConstrainedGlobal(36.),
165// binning for THnSparse
166fMultNbins(0),
167fPtNbins(0),
168fPtCorrNbins(0),
169fEtaNbins(0),
170fZvNbins(0),
171fCentralityNbins(0),
172fBinsMult(0),
173fBinsPt(0),
174fBinsPtCorr(0),
175fBinsEta(0),
176fBinsZv(0),
177fBinsCentrality(0)
178{
179 fMultNbins = 0;
180 fPtNbins = 0;
181 fPtCorrNbins = 0;
182 fEtaNbins = 0;
183 fZvNbins = 0;
184 fCentralityNbins = 0;
185 fBinsMult = 0;
186 fBinsPt = 0;
187 fBinsPtCorr = 0;
188 fBinsEta = 0;
189 fBinsEta = 0;
190 fBinsZv = 0;
191 fBinsCentrality = 0;
192
193 DefineOutput(1, TList::Class());
194}
195
196// destructor
197AlidNdPtAnalysisPbPbAOD::~AlidNdPtAnalysisPbPbAOD()
198{
199 if(hnZvPtEtaCent) delete hnZvPtEtaCent; hnZvPtEtaCent = 0;
200 if(hPt) delete hPt; hPt = 0;
5747f2c6 201 if(hnMCRecPrimZvPtEtaCent) delete hnMCRecPrimZvPtEtaCent; hnMCRecPrimZvPtEtaCent = 0;
202 if(hnMCGenZvPtEtaCent) delete hnMCGenZvPtEtaCent; hnMCGenZvPtEtaCent = 0;
203 if(hnMCRecSecZvPtEtaCent) delete hnMCRecSecZvPtEtaCent; hnMCRecSecZvPtEtaCent = 0;
204 if(hMCPt) delete hMCPt; hMCPt = 0;
205 if(hEventStatistics) delete hEventStatistics; hEventStatistics = 0;
206 if(hEventStatisticsCentrality) delete hEventStatisticsCentrality; hEventStatisticsCentrality = 0;
207 if(hAllEventStatisticsCentrality) delete hAllEventStatisticsCentrality; hAllEventStatisticsCentrality = 0;
208 if(hnZvMultCent) delete hnZvMultCent; hnZvMultCent = 0;
209 if(hTriggerStatistics) delete hTriggerStatistics; hTriggerStatistics = 0;
210 if(hMCTrackPdgCode) delete hMCTrackPdgCode; hMCTrackPdgCode = 0;
211 if(hMCTrackStatusCode) delete hMCTrackStatusCode; hMCTrackStatusCode = 0;
212 if(hCharge) delete hCharge; hCharge = 0;
213 if(hMCCharge) delete hMCCharge; hMCCharge = 0;
214 if(hMCPdgPt) delete hMCPdgPt; hMCPdgPt = 0;
215 if(hMCHijingPrim) delete hMCHijingPrim; hMCHijingPrim = 0;
216 if(hAccNclsTPC) delete hAccNclsTPC; hAccNclsTPC = 0;
217 if(hAccCrossedRowsTPC) delete hAccCrossedRowsTPC; hAccCrossedRowsTPC = 0;
218 if(hDCAPtAll) delete hDCAPtAll; hDCAPtAll = 0;
219 if(hDCAPtAccepted) delete hDCAPtAccepted; hDCAPtAccepted = 0;
220 if(hMCDCAPtSecondary) delete hMCDCAPtSecondary; hMCDCAPtSecondary = 0;
221 if(hMCDCAPtPrimary) delete hMCDCAPtPrimary; hMCDCAPtPrimary = 0;
d25bcbe6 222}
223
224void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
225{
226 // create all output histograms here
227 OpenFile(1, "RECREATE");
228
229 fOutputList = new TList();
230 fOutputList->SetOwner();
231
232 //define default binning
233 Double_t binsMultDefault[48] = {-0.5, 0.5 , 1.5 , 2.5 , 3.5 , 4.5 , 5.5 , 6.5 , 7.5 , 8.5,9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5,19.5, 20.5, 30.5, 40.5 , 50.5 , 60.5 , 70.5 , 80.5 , 90.5 , 100.5,200.5, 300.5, 400.5, 500.5, 600.5, 700.5, 800.5, 900.5, 1000.5, 2000.5, 3000.5, 4000.5, 5000.5, 6000.5, 7000.5, 8000.5, 9000.5, 10000.5 };
234 Double_t binsPtDefault[82] = {0.0, 0.05, 0.1, 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.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0, 140.0, 150.0, 160.0, 180.0, 200.0};
235 Double_t binsPtCorrDefault[37] = {0.,0.05,0.1,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.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,3.0,4.0,200.0};
236 Double_t binsEtaDefault[31] = {-1.5,-1.4,-1.3,-1.2,-1.1,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5};
237 Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
238 Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
239
240 // if no binning is set, use the default
241 if (!fBinsMult) { SetBinsMult(48,binsMultDefault); }
242 if (!fBinsPt) { SetBinsPt(82,binsPtDefault); }
243 if (!fBinsPtCorr) { SetBinsPtCorr(37,binsPtCorrDefault); }
244 if (!fBinsEta) { SetBinsEta(31,binsEtaDefault); }
245 if (!fBinsZv) { SetBinsZv(13,binsZvDefault); }
246 if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
247
248 Int_t binsZvPtEtaCent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
249 Int_t binsZvMultCent[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
250
251 // define Histograms
252 hnZvPtEtaCent = new THnSparseF("hnZvPtEtaCent","Zv:Pt:Eta:Centrality",4,binsZvPtEtaCent);
253 hnZvPtEtaCent->SetBinEdges(0,fBinsZv);
254 hnZvPtEtaCent->SetBinEdges(1,fBinsPt);
255 hnZvPtEtaCent->SetBinEdges(2,fBinsEta);
256 hnZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
257 hnZvPtEtaCent->GetAxis(0)->SetTitle("Zv (cm)");
258 hnZvPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
259 hnZvPtEtaCent->GetAxis(2)->SetTitle("Eta");
260 hnZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
261 hnZvPtEtaCent->Sumw2();
262
263 hnMCRecPrimZvPtEtaCent = new THnSparseF("hnMCRecPrimZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
264 hnMCRecPrimZvPtEtaCent->SetBinEdges(0,fBinsZv);
265 hnMCRecPrimZvPtEtaCent->SetBinEdges(1,fBinsPt);
266 hnMCRecPrimZvPtEtaCent->SetBinEdges(2,fBinsEta);
267 hnMCRecPrimZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
268 hnMCRecPrimZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
269 hnMCRecPrimZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
270 hnMCRecPrimZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
271 hnMCRecPrimZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
272 hnMCRecPrimZvPtEtaCent->Sumw2();
273
274 hnMCGenZvPtEtaCent = new THnSparseF("hnMCGenZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
275 hnMCGenZvPtEtaCent->SetBinEdges(0,fBinsZv);
276 hnMCGenZvPtEtaCent->SetBinEdges(1,fBinsPt);
277 hnMCGenZvPtEtaCent->SetBinEdges(2,fBinsEta);
278 hnMCGenZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
279 hnMCGenZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
280 hnMCGenZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
281 hnMCGenZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
282 hnMCGenZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
283 hnMCGenZvPtEtaCent->Sumw2();
284
5747f2c6 285 hnMCRecSecZvPtEtaCent = new THnSparseF("hnMCRecSecZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
286 hnMCRecSecZvPtEtaCent->SetBinEdges(0,fBinsZv);
287 hnMCRecSecZvPtEtaCent->SetBinEdges(1,fBinsPt);
288 hnMCRecSecZvPtEtaCent->SetBinEdges(2,fBinsEta);
289 hnMCRecSecZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
290 hnMCRecSecZvPtEtaCent->GetAxis(0)->SetTitle("MC Sec Zv (cm)");
291 hnMCRecSecZvPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
292 hnMCRecSecZvPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
293 hnMCRecSecZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
294 hnMCRecSecZvPtEtaCent->Sumw2();
d25bcbe6 295
296 hPt = new TH1F("hPt","hPt",2000,0,200);
297 hPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
298 hPt->GetYaxis()->SetTitle("dN/dp_{T}");
299 hPt->Sumw2();
300
301 hMCPt = new TH1F("hMCPt","hMCPt",2000,0,200);
302 hMCPt->GetXaxis()->SetTitle("MC p_{T} (GeV/c)");
303 hMCPt->GetYaxis()->SetTitle("dN/dp_{T}");
304 hMCPt->Sumw2();
305
306 hEventStatistics = new TH1F("hEventStatistics","hEventStatistics",10,0,10);
307 hEventStatistics->GetYaxis()->SetTitle("number of events");
308 hEventStatistics->SetBit(TH1::kCanRebin);
309
310 hEventStatisticsCentrality = new TH1F("hEventStatisticsCentrality","hEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
311 hEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
312
313 hAllEventStatisticsCentrality = new TH1F("hAllEventStatisticsCentrality","hAllEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
314 hAllEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
315
316 hnZvMultCent = new THnSparseF("hnZvMultCent","Zv:mult:Centrality",3,binsZvMultCent);
317 hnZvMultCent->SetBinEdges(0,fBinsZv);
318 hnZvMultCent->SetBinEdges(1,fBinsMult);
319 hnZvMultCent->SetBinEdges(2,fBinsCentrality);
320 hnZvMultCent->GetAxis(0)->SetTitle("Zv (cm)");
321 hnZvMultCent->GetAxis(1)->SetTitle("N_{acc}");
322 hnZvMultCent->GetAxis(2)->SetTitle("Centrality");
323 hnZvMultCent->Sumw2();
324
325 hTriggerStatistics = new TH1F("hTriggerStatistics","hTriggerStatistics",10,0,10);
326 hTriggerStatistics->GetYaxis()->SetTitle("number of events");
327
328 hMCTrackPdgCode = new TH1F("hMCTrackPdgCode","hMCTrackPdgCode",100,0,10);
329 hMCTrackPdgCode->GetYaxis()->SetTitle("number of tracks");
330 hMCTrackPdgCode->SetBit(TH1::kCanRebin);
331
332 hMCTrackStatusCode = new TH1F("hMCTrackStatusCode","hMCTrackStatusCode",100,0,10);
333 hMCTrackStatusCode->GetYaxis()->SetTitle("number of tracks");
334 hMCTrackStatusCode->SetBit(TH1::kCanRebin);
335
336 hCharge = new TH1F("hCharge","hCharge",30, -5, 5);
337 hCharge->GetXaxis()->SetTitle("Charge");
338 hCharge->GetYaxis()->SetTitle("number of tracks");
339
340 hMCCharge = new TH1F("hMCCharge","hMCCharge",30, -5, 5);
341 hMCCharge->GetXaxis()->SetTitle("MC Charge");
342 hMCCharge->GetYaxis()->SetTitle("number of tracks");
343
344 hMCPdgPt = new TH2F("hMCPdgPt","hMCPdgPt",fPtNbins-1, fBinsPt, 100,0,100);
345 hMCPdgPt->GetYaxis()->SetTitle("particle");
346 hMCPdgPt->GetXaxis()->SetTitle("Pt (GeV/c)");
347
348 hMCHijingPrim = new TH1F("hMCHijingPrim","hMCHijingPrim",2,0,2);
349 hMCPdgPt->GetYaxis()->SetTitle("number of particles");
350
351 hAccNclsTPC = new TH1F("hAccNclsTPC","hAccNclsTPC",160,0,159);
352 hAccNclsTPC->GetXaxis()->SetTitle("number of clusters per track after cut");
353
354 hAccCrossedRowsTPC = new TH1F("hAccCrossedRowsTPC","hAccCrossedRowsTPC",160,0,159);
355 hAccCrossedRowsTPC->GetXaxis()->SetTitle("number of crossed rows per track after cut");
356
d0483ba3 357 Int_t binsDCAxyDCAzPt[3] = { 200,200, fPtNbins-1};
358 Double_t minDCAxyDCAzPt[3] = { -10, -10, 0};
359 Double_t maxDCAxyDCAzPt[3] = { 10., 10., 100};
360
361 hDCAPtAll = new THnSparseF("hDCAPtAll","hDCAPtAll;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
362 hDCAPtAccepted = new THnSparseF("hDCAPtAccepted","hDCAPtAccepted;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
363 hMCDCAPtSecondary = new THnSparseF("hMCDCAPtSecondary","hMCDCAPtSecondary;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
364 hMCDCAPtPrimary = new THnSparseF("hMCDCAPtPrimary","hMCDCAPtPrimary;p_{T} (GeV/c);DCAz",3, binsDCAxyDCAzPt, minDCAxyDCAzPt, maxDCAxyDCAzPt);
365
366 hDCAPtAll->SetBinEdges(2, fBinsPt);
367 hDCAPtAccepted->SetBinEdges(2, fBinsPt);
368 hMCDCAPtSecondary->SetBinEdges(2, fBinsPt);
369 hMCDCAPtPrimary->SetBinEdges(2, fBinsPt);
d25bcbe6 370
371 // Add Histos, Profiles etc to List
372 fOutputList->Add(hnZvPtEtaCent);
373 fOutputList->Add(hPt);
374 fOutputList->Add(hnMCRecPrimZvPtEtaCent);
375 fOutputList->Add(hnMCGenZvPtEtaCent);
5747f2c6 376 fOutputList->Add(hnMCRecSecZvPtEtaCent);
d25bcbe6 377 fOutputList->Add(hMCPt);
378 fOutputList->Add(hEventStatistics);
379 fOutputList->Add(hEventStatisticsCentrality);
380 fOutputList->Add(hAllEventStatisticsCentrality);
381 fOutputList->Add(hnZvMultCent);
382 fOutputList->Add(hTriggerStatistics);
383 fOutputList->Add(hMCTrackPdgCode);
384 fOutputList->Add(hMCTrackStatusCode);
385 fOutputList->Add(hCharge);
386 fOutputList->Add(hMCCharge);
387 fOutputList->Add(hMCPdgPt);
388 fOutputList->Add(hMCHijingPrim);
389 fOutputList->Add(hAccNclsTPC);
390 fOutputList->Add(hAccCrossedRowsTPC);
5747f2c6 391 fOutputList->Add(hDCAPtAll);
392 fOutputList->Add(hDCAPtAccepted);
393 fOutputList->Add(hMCDCAPtSecondary);
394 fOutputList->Add(hMCDCAPtPrimary);
395
d25bcbe6 396
397 PostData(1, fOutputList);
398}
399
400void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
401{
f856053f 402
d25bcbe6 403 // Main Loop
404 // called for each event
405 hEventStatistics->Fill("all events",1);
406
407 // set ZERO pointers:
408 AliInputEventHandler *inputHandler = NULL;
409 AliAODTrack *track = NULL;
410 AliAODMCParticle *mcPart = NULL;
411 AliAODMCHeader *mcHdr = NULL;
412 AliGenHijingEventHeader *genHijingHeader = NULL;
413 AliGenPythiaEventHeader *genPythiaHeader = NULL;
414
415 Bool_t bIsEventSelectedMB = kFALSE;
416 Bool_t bIsEventSelectedSemi = kFALSE;
417 Bool_t bIsEventSelectedCentral = kFALSE;
418 Bool_t bIsEventSelected = kFALSE;
419 Bool_t bIsPrimary = kFALSE;
420 Bool_t bIsHijingParticle = kFALSE;
421 Bool_t bIsPythiaParticle = kFALSE;
422 Bool_t bEventHasATrack = kFALSE;
423 Bool_t bEventHasATrackInRange = kFALSE;
424 Int_t nTriggerFired = 0;
425
426
427 Double_t dMCTrackZvPtEtaCent[4] = {0};
428 Double_t dTrackZvPtEtaCent[4] = {0};
429
430 Double_t dMCEventZv = -100;
431 Double_t dEventZv = -100;
432 Int_t iAcceptedMultiplicity = 0;
433
434 bIsMonteCarlo = kFALSE;
435
436 AliAODEvent *eventAOD = 0x0;
437 eventAOD = dynamic_cast<AliAODEvent*>( InputEvent() );
438 if (!eventAOD) {
439 AliWarning("ERROR: eventAOD not available \n");
440 return;
441 }
442
443 // check, which trigger has been fired
444 inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
445 bIsEventSelectedMB = ( inputHandler->IsEventSelected() & AliVEvent::kMB);
446 bIsEventSelectedSemi = ( inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
447 bIsEventSelectedCentral = ( inputHandler->IsEventSelected() & AliVEvent::kCentral);
448
449 if(bIsEventSelectedMB || bIsEventSelectedSemi || bIsEventSelectedCentral) hTriggerStatistics->Fill("all triggered events",1);
450 if(bIsEventSelectedMB) { hTriggerStatistics->Fill("MB trigger",1); nTriggerFired++; }
451 if(bIsEventSelectedSemi) { hTriggerStatistics->Fill("SemiCentral trigger",1); nTriggerFired++; }
452 if(bIsEventSelectedCentral) { hTriggerStatistics->Fill("Central trigger",1); nTriggerFired++; }
453 if(nTriggerFired == 0) { hTriggerStatistics->Fill("No trigger",1); }
454
455 bIsEventSelected = ( inputHandler->IsEventSelected() & GetCollisionCandidates() );
456
457 // only take tracks of events, which are triggered
458 if(nTriggerFired == 0) { return; }
459
460 // if( !bIsEventSelected || nTriggerFired>1 ) return;
461
462 // hEventStatistics->Fill("events with only coll. cand.", 1);
463
464
465
466 // check if there is a stack, if yes, then do MC loop
467 TList *list = eventAOD->GetList();
468 TClonesArray *stack = 0x0;
469 stack = (TClonesArray*)list->FindObject(AliAODMCParticle::StdBranchName());
470
471 if( stack )
472 {
473 bIsMonteCarlo = kTRUE;
474
475 mcHdr = (AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
476
477 genHijingHeader = GetHijingEventHeader(mcHdr);
478 // genPythiaHeader = GetPythiaEventHeader(mcHdr);
479
480 if(!genHijingHeader) { return; }
481
482 // if(!genPythiaHeader) { return; }
483
484 dMCEventZv = mcHdr->GetVtxZ();
485 dMCTrackZvPtEtaCent[0] = dMCEventZv;
486 hEventStatistics->Fill("MC all events",1);
487 }
488
489 AliCentrality* aCentrality = eventAOD->GetCentrality();
490 Double_t dCentrality = aCentrality->GetCentralityPercentile("V0M");
491
492 if( dCentrality < 0 ) return;
493 hEventStatistics->Fill("after centrality selection",1);
494
495
496
497 // start with MC truth analysis
498 if(bIsMonteCarlo)
499 {
500
501 if( dMCEventZv > dCutMaxZVertex ) { return; }
502
503 dMCTrackZvPtEtaCent[0] = dMCEventZv;
504
505 hEventStatistics->Fill("MC afterZv cut",1);
506
507 for(Int_t iMCtrack = 0; iMCtrack < stack->GetEntriesFast(); iMCtrack++)
508 {
509 mcPart =(AliAODMCParticle*)stack->At(iMCtrack);
510
511 // check for charge
512 if( !(IsMCTrackAccepted(mcPart)) ) continue;
513
514 if(!IsHijingParticle(mcPart, genHijingHeader)) { continue; }
515
516 if(mcPart->IsPhysicalPrimary() )
517 {
518 hMCHijingPrim->Fill("IsPhysicalPrimary",1);
519 }
520 else
521 {
522 hMCHijingPrim->Fill("NOT a primary",1);
523 continue;
524 }
525
d0483ba3 526
d25bcbe6 527 //
528 // ======================== fill histograms ========================
529 dMCTrackZvPtEtaCent[1] = mcPart->Pt();
530 dMCTrackZvPtEtaCent[2] = mcPart->Eta();
531 dMCTrackZvPtEtaCent[3] = dCentrality;
532
533 bEventHasATrack = kTRUE;
534
535 hnMCGenZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
536
537 if( (dMCTrackZvPtEtaCent[1] > dCutPtMin) &&
538 (dMCTrackZvPtEtaCent[1] < dCutPtMax) &&
539 (dMCTrackZvPtEtaCent[2] > dCutEtaMin) &&
540 (dMCTrackZvPtEtaCent[2] < dCutEtaMax) )
541 {
542 hMCPt->Fill(mcPart->Pt());
543 hMCCharge->Fill(mcPart->Charge()/3.);
544 bEventHasATrackInRange = kTRUE;
545 }
546
547 }
548 } // isMonteCarlo
549 if(bEventHasATrack) { hEventStatistics->Fill("MC events with tracks",1); }
550 if(bEventHasATrackInRange) { hEventStatistics->Fill("MC events with tracks in range",1); }
551 bEventHasATrack = kFALSE;
552 bEventHasATrackInRange = kFALSE;
553
554
555
556 // Loop over recontructed tracks
557
558 dEventZv = eventAOD->GetPrimaryVertex()->GetZ();
559 if( TMath::Abs(dEventZv) > dCutMaxZVertex ) return;
560
561 hAllEventStatisticsCentrality->Fill(dCentrality/*, nTriggerFired*/);
562
563 hEventStatistics->Fill("after Zv cut",1);
564
565 dTrackZvPtEtaCent[0] = dEventZv;
566
567 for(Int_t itrack = 0; itrack < eventAOD->GetNumberOfTracks(); itrack++)
568 {
569 track = eventAOD->GetTrack(itrack);
570 if(!track) continue;
571
572 mcPart = NULL;
573 dMCTrackZvPtEtaCent[1] = 0;
574 dMCTrackZvPtEtaCent[2] = 0;
575 dMCTrackZvPtEtaCent[3] = 0;
576
577 bIsPrimary = kFALSE;
578
d0483ba3 579 Double_t dDCAxyDCAzPt[3] = { GetDCAxy(track, eventAOD), GetDCAz(track, eventAOD), track->Pt() };
580
581 hDCAPtAll->Fill(dDCAxyDCAzPt);
582
d25bcbe6 583 if( !(IsTrackAccepted(track)) ) continue;
584
5747f2c6 585 dTrackZvPtEtaCent[1] = track->Pt();
586 dTrackZvPtEtaCent[2] = track->Eta();
587 dTrackZvPtEtaCent[3] = dCentrality;
588
d25bcbe6 589 if( bIsMonteCarlo )
590 {
591 mcPart = (AliAODMCParticle*)stack->At(TMath::Abs(track->GetLabel()));
592 if( !mcPart ) { continue; }
593
594 // check for charge
595 if( !(IsMCTrackAccepted(mcPart)) ) { continue; }
596
597 bIsHijingParticle = IsHijingParticle(mcPart, genHijingHeader);
598 // bIsPythiaParticle = IsPythiaParticle(mcPart, genPythiaHeader);
599
600 // if(!bIsHijingParticle) continue; // only take real tracks, not injected ones
601
602 bIsPrimary = mcPart->IsPhysicalPrimary();
603
604 dMCTrackZvPtEtaCent[1] = mcPart->Pt();
605 dMCTrackZvPtEtaCent[2] = mcPart->Eta();
606 dMCTrackZvPtEtaCent[3] = dCentrality;
607
608 if(bIsPrimary && bIsHijingParticle)
609 {
d0483ba3 610 hnMCRecPrimZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
611 hMCDCAPtPrimary->Fill(dDCAxyDCAzPt);
d25bcbe6 612 }
613
614 if(!bIsPrimary /*&& !bIsHijingParticle*/)
615 {
616 Int_t indexMoth = mcPart->GetMother();
617 if(indexMoth >= 0)
618 {
619 AliAODMCParticle* moth = (AliAODMCParticle*)stack->At(indexMoth);
620 Bool_t bMotherIsHijingParticle = IsHijingParticle(moth, genHijingHeader);
621
622 if(bMotherIsHijingParticle) // only store secondaries, which come from a not embedded signal!
623 {
624 hMCTrackStatusCode->Fill(Form("%d",mcPart->GetStatus()), 1);
625 if(TMath::Abs(mcPart->Eta()) < 0.8) { hMCPdgPt->Fill(mcPart->Pt(), Form("%s",GetParticleName(mcPart->GetPdgCode())), 1); }
626
d0483ba3 627 hnMCRecSecZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
628 hMCDCAPtSecondary->Fill(dDCAxyDCAzPt);
d25bcbe6 629 hMCTrackPdgCode->Fill(Form("%s_H%i_H%i",GetParticleName(moth->GetPdgCode()),bMotherIsHijingParticle, bIsHijingParticle), 1);
630 // delete moth;
631 }
632 }
633 }
634 } // end isMonteCarlo
635
636 // ======================== fill histograms ========================
637
d0483ba3 638 // if(bIsMonteCarlo && !bIsHijingParticle)
639 // {
640 // continue; //only store reco tracks, which do not come from embedded signal
641 // }
642
643
644
645 bEventHasATrack = kTRUE;
646
647 hnZvPtEtaCent->Fill(dTrackZvPtEtaCent);
648 hDCAPtAccepted->Fill(dDCAxyDCAzPt);
649
650 if( (dTrackZvPtEtaCent[1] > dCutPtMin) &&
651 (dTrackZvPtEtaCent[1] < dCutPtMax) &&
652 (dTrackZvPtEtaCent[2] > dCutEtaMin) &&
653 (dTrackZvPtEtaCent[2] < dCutEtaMax) )
654 {
655 iAcceptedMultiplicity++;
656 bEventHasATrackInRange = kTRUE;
657 hPt->Fill(track->Pt());
658 hCharge->Fill(track->Charge());
659 }
d25bcbe6 660 } // end track loop
661
662 if(bEventHasATrack) { hEventStatistics->Fill("events with tracks",1); bEventHasATrack = kFALSE; }
663
664 if(bEventHasATrackInRange)
665 {
666 hEventStatistics->Fill("events with tracks in range",1);
667 hEventStatisticsCentrality->Fill(dCentrality);
668 bEventHasATrackInRange = kFALSE;
669 }
670
671 Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
672 hnZvMultCent->Fill(dEventZvMultCent);
673
674
675
676 PostData(1, fOutputList);
677
678}
679
680Bool_t AlidNdPtAnalysisPbPbAOD::IsTrackAccepted(AliAODTrack *tr)
681{
682 if(!tr) return kFALSE;
683
684 if(tr->Charge()==0) { return kFALSE; }
685
686 if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { return kFALSE; }
687
d25bcbe6 688 Double_t dNClustersTPC = tr->GetTPCNcls();
689 Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
690
f856053f 691 if(dCrossedRowsTPC < GetCutMinNCrossedRowsTPC()) { return kFALSE; }
d0483ba3 692
d25bcbe6 693 hAccNclsTPC->Fill(dNClustersTPC);
694 hAccCrossedRowsTPC->Fill(dCrossedRowsTPC);
695 // Double_t dFindableClustersTPC = tr->GetTPCNclsF();
696 // Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
697 //
698 // Bool_t bIsFromKink = kFALSE;
699 // if(tr->GetProdVertex()->GetType() == AliAODVertex::kKink) bIsFromKink = kTRUE;
700 //
701 // // from AliAnalysisTaskPIDqa.cxx
702 // ULong_t uStatus = tr->GetStatus();
703 // Bool_t bHasRefitTPC = kFALSE;
704 // Bool_t bHasRefitITS = kFALSE;
705 //
706 // if ((uStatus & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) bHasRefitTPC = kTRUE;
707 // if ((uStatus & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) bHasRefitITS = kTRUE;
708 //
709 // // from AliDielectronVarManager.h
710 // Bool_t bHasHitInSPD = kFALSE;
711 // for (Int_t iC=0; iC<2; iC++)
712 // {
713 // if (((tr->GetITSClusterMap()) & (1<<(iC))) > 0) { bHasHitInSPD = kTRUE; }
714 // }
715 //
716 // Double_t dNClustersITS = tr->GetITSNcls();
717
718 // cuts to be done:
719 // TPC
720 // esdTrackCuts->SetMinNCrossedRowsTPC(70);
721 // esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
722 //
723 // esdTrackCuts->SetMaxChi2PerClusterTPC(4);
724 // esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
725 // esdTrackCuts->SetRequireTPCRefit(kTRUE);
726 // ITS
727 // esdTrackCuts->SetRequireITSRefit(kTRUE);
728 // esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
729 //
730 // esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
731 // esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
732 //
733 // esdTrackCuts->SetMaxDCAToVertexZ(2);
734 // esdTrackCuts->SetDCAToVertex2D(kFALSE);
735 // esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
736 //
737 // esdTrackCuts->SetMaxChi2PerClusterITS(36);
738
739
740 return kTRUE;
741}
742
d0483ba3 743Double_t AlidNdPtAnalysisPbPbAOD::GetDCAz(AliAODTrack *track, AliAODEvent *event)
744{
745 return GetDCA(track, event, kTRUE);
746}
747
748Double_t AlidNdPtAnalysisPbPbAOD::GetDCAxy(AliAODTrack *track, AliAODEvent *event)
749{
750 return GetDCA(track, event, kFALSE);
751}
752
753Double_t AlidNdPtAnalysisPbPbAOD::GetDCA(AliAODTrack *tr, AliAODEvent *evt, Bool_t bDCAz)
754{
755 if(!tr) return -999.;
756
757 if(tr->TestBit(AliAODTrack::kIsDCA))
758 {
759 if(bDCAz) return tr->ZAtDCA();
760 else return sqrt(tr->XAtDCA()*tr->XAtDCA() + tr->YAtDCA()*tr->YAtDCA());
761 }
762
763 Bool_t ok=kFALSE;
764 Double_t d0z0[2];
765 if(evt)
766 {
767 Double_t covd0z0[3];
768
769 AliExternalTrackParam etp;
770 etp.CopyFromVTrack(tr);
771
772 Float_t xstart = etp.GetX();
773 if(xstart>3.)
774 {
775 d0z0[0]=-999.;
776 d0z0[1]=-999.;
777 //printf("This method can be used only for propagation inside the beam pipe \n");
778 if(bDCAz) return d0z0[1];
779 else return d0z0[0];
780 }
781
782
783 AliAODVertex *vtx =(AliAODVertex*)(evt->GetPrimaryVertex());
784 Double_t fBzkG = evt->GetMagneticField(); // z componenent of field in kG
785 ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
786 }
787 if(!ok){
788 d0z0[0]=-999.;
789 d0z0[1]=-999.;
790 if(bDCAz) return d0z0[1];
791 else return d0z0[0];
792 }
793 if(bDCAz) return d0z0[1];
794 else return d0z0[0];
795}
796
d25bcbe6 797Bool_t AlidNdPtAnalysisPbPbAOD::IsMCTrackAccepted(AliAODMCParticle *part)
798{
799 if(!part) return kFALSE;
800
801 Double_t charge = part->Charge()/3.;
802 if (TMath::Abs(charge) < 0.001) return kFALSE;
803
804 return kTRUE;
805}
806
807const char * AlidNdPtAnalysisPbPbAOD::GetParticleName(Int_t pdg)
808{
809 TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
810 if(p1) return p1->GetName();
811 return Form("%d", pdg);
812}
813
814AliGenHijingEventHeader* AlidNdPtAnalysisPbPbAOD::GetHijingEventHeader(AliAODMCHeader *header)
815{
816 //
817 // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
818 //
819
820 if(!header) return 0x0;
821 AliGenHijingEventHeader* hijingGenHeader = NULL;
822
823 TList* headerList = header->GetCocktailHeaders();
824
825 for(Int_t i = 0; i < headerList->GetEntries(); i++)
826 {
827 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
828 if(hijingGenHeader) break;
829 }
830
831 if(!hijingGenHeader) return 0x0;
832
833 return hijingGenHeader;
834}
835
836AliGenPythiaEventHeader* AlidNdPtAnalysisPbPbAOD::GetPythiaEventHeader(AliAODMCHeader *header)
837{
838 //
839 // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
840 //
841
842 if(!header) return 0x0;
843 AliGenPythiaEventHeader* PythiaGenHeader = NULL;
844
845 TList* headerList = header->GetCocktailHeaders();
846
847 for(Int_t i = 0; i < headerList->GetEntries(); i++)
848 {
849 PythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
850 if(PythiaGenHeader) break;
851 }
852
853 if(!PythiaGenHeader) return 0x0;
854
855 return PythiaGenHeader;
856}
857
858//________________________________________________________________________
859Bool_t AlidNdPtAnalysisPbPbAOD::IsHijingParticle(const AliAODMCParticle *part, AliGenHijingEventHeader* hijingGenHeader){
860
861 // Check whether a particle is from Hijing or some injected
862
863 if(part->Label() > (hijingGenHeader->NProduced()-1)) return kFALSE;
864 return kTRUE;
865}
866
867//________________________________________________________________________
868Bool_t AlidNdPtAnalysisPbPbAOD::IsPythiaParticle(const AliAODMCParticle *part, AliGenPythiaEventHeader* pythiaGenHeader){
869
870 // Check whether a particle is from Pythia or some injected
871
872 if(part->Label() > (pythiaGenHeader->NProduced()-1)) return kFALSE;
873 return kTRUE;
874}
875
5747f2c6 876// Int_t AlidNdPtAnalysisPbPbAOD::IsMCSecondary(AliAODMCParticle *part, TClonesArray *arrayMC)
877// {
d0483ba3 878 // //
879 // // adapted from AliAnalysisTaskSpectraAOD.cxx
880 // //
881 // // returns
882 // // -1: no particle
883 // // 0: is primary
884 // // 1: is secondary from weak
885 // // 2: is secondary from material
886 //
887 // // usage for studies, currrently not implemented
888 //
889 // if(!part) return -1;
890 //
891 // if( part->IsPhysicalPrimary() ) return 0;
892 //
893 // Bool_t isSecondaryMaterial = kFALSE;
894 // Bool_t isSecondaryWeak = kFALSE;
895 // Int_t mfl = -999;
896 // Int_t codemoth = -999;
897 // Int_t indexMoth = part->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()
898 // if(indexMoth >= 0)
899 // {
900 // AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
901 // codemoth = TMath::Abs(moth->GetPdgCode());
902 // mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
903 // }
904 // // add if(partMC->GetStatus() & kPDecay)? FIXME
905 // if(mfl==3) isSecondaryWeak = kTRUE;
906 // else isSecondaryMaterial = kTRUE;
907 //
908 // if(isSecondaryWeak) return 1;
909 // if(isSecondaryMaterial) return 2;
910 //
911 // // if( isSecondaryMaterial || isSecondaryWeak ) return kTRUE;
912 //
913 // // return kFALSE; this line will not be reached, as either isSecondaryMaterial or isSecondaryWeak is true!
914 // // removed due to coverity
915 // }
916
917
918
919 void AlidNdPtAnalysisPbPbAOD::Terminate(Option_t *)
920 {
921
922 }
923
924 Double_t* AlidNdPtAnalysisPbPbAOD::GetArrayClone(Int_t n, Double_t* source)
925 {
926 if (!source || n==0) return 0;
927 Double_t* dest = new Double_t[n];
928 for (Int_t i=0; i<n ; i++) { dest[i] = source[i]; }
929 return dest;
930 }
931