]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Added distributions in phi, updated AddTask accordingly
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtAnalysisPbPbAOD.cxx
... / ...
CommitLineData
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// last modified: 17.10.2013
20//------------------------------------------------------------------------------
21
22
23#include "AlidNdPtAnalysisPbPbAOD.h"
24
25
26using namespace std;
27
28ClassImp(AlidNdPtAnalysisPbPbAOD)
29
30
31
32AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD(const char *name) : AliAnalysisTaskSE(name),
33fOutputList(0),
34// Histograms
35fPt(0),
36fMCPt(0),
37fZvPtEtaCent(0),
38fPhiPtEtaCent(0),
39fMCRecPrimZvPtEtaCent(0),
40fMCGenZvPtEtaCent(0),
41fMCRecSecZvPtEtaCent(0),
42fMCRecPrimPhiPtEtaCent(0),
43fMCGenPhiPtEtaCent(0),
44fMCRecSecPhiPtEtaCent(0),
45fEventStatistics(0),
46fEventStatisticsCentrality(0),
47fMCEventStatisticsCentrality(0),
48fAllEventStatisticsCentrality(0),
49fEventStatisticsCentralityTrigger(0),
50fZvMultCent(0),
51fTriggerStatistics(0),
52fMCTrackPdgCode(0),
53fMCTrackStatusCode(0),
54fCharge(0),
55fMCCharge(0),
56fMCPdgPt(0),
57fMCHijingPrim(0),
58fDCAPtAll(0),
59fDCAPtAccepted(0),
60fMCDCAPtSecondary(0),
61fMCDCAPtPrimary(0),
62fCutPercClusters(0),
63fCutPercCrossed(0),
64fCrossCheckRowsLength(0),
65fCrossCheckClusterLength(0),
66fCrossCheckRowsLengthAcc(0),
67fCrossCheckClusterLengthAcc(0),
68//global
69fIsMonteCarlo(0),
70// event cut variables
71fCutMaxZVertex(10.),
72// track kinematic cut variables
73fCutPtMin(0.15),
74fCutPtMax(200.),
75fCutEtaMin(-0.8),
76fCutEtaMax(0.8),
77// track quality cut variables
78fUseRelativeCuts(kFALSE),
79fCutRequireTPCRefit(kTRUE),
80fCutMinNumberOfClusters(60),
81fCutPercMinNumberOfClusters(0.2),
82fCutMinNumberOfCrossedRows(120.),
83fCutPercMinNumberOfCrossedRows(0.2),
84fCutMinRatioCrossedRowsOverFindableClustersTPC(0.8),
85fCutMaxChi2PerClusterTPC(4.),
86fCutMaxFractionSharedTPCClusters(0.4),
87fCutMaxDCAToVertexZ(3.0),
88fCutMaxDCAToVertexXY(3.0),
89fCutRequireITSRefit(kTRUE),
90fCutMaxChi2PerClusterITS(36.),
91fCutDCAToVertex2D(kFALSE),
92fCutRequireSigmaToVertex(kFALSE),
93fCutMaxDCAToVertexXYPtDepPar0(0.0182),
94fCutMaxDCAToVertexXYPtDepPar1(0.0350),
95fCutMaxDCAToVertexXYPtDepPar2(1.01),
96fCutAcceptKinkDaughters(kFALSE),
97fCutMaxChi2TPCConstrainedGlobal(36.),
98fCutLengthInTPCPtDependent(kFALSE),
99fPrefactorLengthInTPCPtDependent(1),
100// binning for THnSparse
101fMultNbins(0),
102fPtNbins(0),
103fPtCorrNbins(0),
104fPtCheckNbins(0),
105fEtaNbins(0),
106fEtaCheckNbins(0),
107fZvNbins(0),
108fCentralityNbins(0),
109fPhiNbins(0),
110fBinsMult(0),
111fBinsPt(0),
112fBinsPtCorr(0),
113fBinsPtCheck(0),
114fBinsEta(0),
115fBinsEtaCheck(0),
116fBinsZv(0),
117fBinsCentrality(0),
118fBinsPhi(0)
119{
120
121 for(Int_t i = 0; i < cqMax; i++)
122 {
123 fCrossCheckAll[i] = 0;
124 fCrossCheckAcc[i] = 0;
125 }
126
127 fMultNbins = 0;
128 fPtNbins = 0;
129 fPtCorrNbins = 0;
130 fPtCheckNbins = 0;
131 fEtaNbins = 0;
132 fEtaCheckNbins = 0;
133 fZvNbins = 0;
134 fCentralityNbins = 0;
135 fBinsMult = 0;
136 fBinsPt = 0;
137 fBinsPtCorr = 0;
138 fBinsPtCheck = 0;
139 fBinsEta = 0;
140 fBinsEtaCheck = 0;
141 fBinsZv = 0;
142 fBinsCentrality = 0;
143 fBinsPhi = 0;
144
145 DefineOutput(1, TList::Class());
146}
147
148// destructor
149AlidNdPtAnalysisPbPbAOD::~AlidNdPtAnalysisPbPbAOD()
150{
151
152 if(fZvPtEtaCent) delete fZvPtEtaCent; fZvPtEtaCent = 0;
153 if(fPt) delete fPt; fPt = 0;
154 if(fMCRecPrimZvPtEtaCent) delete fMCRecPrimZvPtEtaCent; fMCRecPrimZvPtEtaCent = 0;
155 if(fMCGenZvPtEtaCent) delete fMCGenZvPtEtaCent; fMCGenZvPtEtaCent = 0;
156 if(fMCRecSecZvPtEtaCent) delete fMCRecSecZvPtEtaCent; fMCRecSecZvPtEtaCent = 0;
157 if(fMCPt) delete fMCPt; fMCPt = 0;
158 if(fEventStatistics) delete fEventStatistics; fEventStatistics = 0;
159 if(fEventStatisticsCentrality) delete fEventStatisticsCentrality; fEventStatisticsCentrality = 0;
160 if(fMCEventStatisticsCentrality) delete fMCEventStatisticsCentrality; fMCEventStatisticsCentrality = 0;
161 if(fAllEventStatisticsCentrality) delete fAllEventStatisticsCentrality; fAllEventStatisticsCentrality = 0;
162 if(fEventStatisticsCentralityTrigger) delete fEventStatisticsCentralityTrigger; fEventStatisticsCentralityTrigger = 0;
163 if(fZvMultCent) delete fZvMultCent; fZvMultCent = 0;
164 if(fTriggerStatistics) delete fTriggerStatistics; fTriggerStatistics = 0;
165 if(fMCTrackPdgCode) delete fMCTrackPdgCode; fMCTrackPdgCode = 0;
166 if(fMCTrackStatusCode) delete fMCTrackStatusCode; fMCTrackStatusCode = 0;
167 if(fCharge) delete fCharge; fCharge = 0;
168 if(fMCCharge) delete fMCCharge; fMCCharge = 0;
169 if(fMCPdgPt) delete fMCPdgPt; fMCPdgPt = 0;
170 if(fMCHijingPrim) delete fMCHijingPrim; fMCHijingPrim = 0;
171 if(fDCAPtAll) delete fDCAPtAll; fDCAPtAll = 0;
172 if(fDCAPtAccepted) delete fDCAPtAccepted; fDCAPtAccepted = 0;
173 if(fMCDCAPtSecondary) delete fMCDCAPtSecondary; fMCDCAPtSecondary = 0;
174 if(fMCDCAPtPrimary) delete fMCDCAPtPrimary; fMCDCAPtPrimary = 0;
175 if(fMCDCAPtSecondary) delete fMCDCAPtSecondary; fMCDCAPtSecondary = 0;
176 if(fMCDCAPtPrimary) delete fMCDCAPtPrimary; fMCDCAPtPrimary = 0;
177
178 for(Int_t i = 0; i < cqMax; i++)
179 {
180 if(fCrossCheckAll[i]) delete fCrossCheckAll[i]; fCrossCheckAll[i] = 0;
181 if(fCrossCheckAcc[i]) delete fCrossCheckAcc[i]; fCrossCheckAcc[i] = 0;
182 }
183
184}
185
186void AlidNdPtAnalysisPbPbAOD::UserCreateOutputObjects()
187{
188 // create all output histograms here
189 OpenFile(1, "RECREATE");
190
191 fOutputList = new TList();
192 fOutputList->SetOwner();
193
194 //define default binning
195 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 };
196 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};
197 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};
198 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};
199 Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
200 Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};
201
202 Double_t binsPhiDefault[37] = { 0, 0.174533, 0.349066, 0.523599, 0.698132, 0.872665, 1.0472, 1.22173, 1.39626, 1.5708, 1.74533, 1.91986, 2.0944, 2.26893, 2.44346, 2.61799, 2.79253, 2.96706, 3.14159, 3.31613, 3.49066, 3.66519, 3.83972, 4.01426, 4.18879, 4.36332, 4.53786, 4.71239, 4.88692, 5.06145, 5.23599, 5.41052, 5.58505, 5.75959, 5.93412, 6.10865, 2.*TMath::Pi()};
203
204 Double_t binsPtCheckDefault[20] = {0.,0.15,0.5,1.0,2.0,3.0,4.0, 5.0, 10.0, 13.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 70.0, 100.0, 150.0, 200.0};
205 Double_t binsEtaCheckDefault[7] = {-1.0,-0.8,-0.4,0.,0.4,0.8,1.0};
206
207 // if no binning is set, use the default
208 if (!fBinsMult) { SetBinsMult(48,binsMultDefault); }
209 if (!fBinsPt) { SetBinsPt(82,binsPtDefault); }
210 if (!fBinsPtCorr) { SetBinsPtCorr(37,binsPtCorrDefault); }
211 if (!fBinsPtCheck) { SetBinsPtCheck(20,binsPtCheckDefault); }
212 if (!fBinsEta) { SetBinsEta(31,binsEtaDefault); }
213 if (!fBinsEtaCheck) { SetBinsEtaCheck(7,binsEtaCheckDefault); }
214 if (!fBinsZv) { SetBinsZv(13,binsZvDefault); }
215 if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
216 if (!fBinsPhi) { SetBinsPhi(37,binsPhiDefault); }
217
218 Int_t binsZvPtEtaCent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
219 Int_t binsPhiPtEtaCent[4]={fPhiNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
220 Int_t binsZvMultCent[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
221
222 // define Histograms
223 fZvPtEtaCent = new THnSparseF("fZvPtEtaCent","Zv:Pt:Eta:Centrality",4,binsZvPtEtaCent);
224 fZvPtEtaCent->SetBinEdges(0,fBinsZv);
225 fZvPtEtaCent->SetBinEdges(1,fBinsPt);
226 fZvPtEtaCent->SetBinEdges(2,fBinsEta);
227 fZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
228 fZvPtEtaCent->GetAxis(0)->SetTitle("Zv (cm)");
229 fZvPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
230 fZvPtEtaCent->GetAxis(2)->SetTitle("Eta");
231 fZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
232 fZvPtEtaCent->Sumw2();
233
234 fPhiPtEtaCent = new THnSparseF("fPhiPtEtaCent","Phi:Pt:Eta:Centrality",4,binsPhiPtEtaCent);
235 fPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
236 fPhiPtEtaCent->SetBinEdges(1,fBinsPt);
237 fPhiPtEtaCent->SetBinEdges(2,fBinsEta);
238 fPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
239 fPhiPtEtaCent->GetAxis(0)->SetTitle("Phi (cm)");
240 fPhiPtEtaCent->GetAxis(1)->SetTitle("Pt (GeV/c)");
241 fPhiPtEtaCent->GetAxis(2)->SetTitle("Eta");
242 fPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
243 fPhiPtEtaCent->Sumw2();
244
245
246
247 fMCRecPrimZvPtEtaCent = new THnSparseF("fMCRecPrimZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
248 fMCRecPrimZvPtEtaCent->SetBinEdges(0,fBinsZv);
249 fMCRecPrimZvPtEtaCent->SetBinEdges(1,fBinsPt);
250 fMCRecPrimZvPtEtaCent->SetBinEdges(2,fBinsEta);
251 fMCRecPrimZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
252 fMCRecPrimZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
253 fMCRecPrimZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
254 fMCRecPrimZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
255 fMCRecPrimZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
256 fMCRecPrimZvPtEtaCent->Sumw2();
257
258 fMCGenZvPtEtaCent = new THnSparseF("fMCGenZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
259 fMCGenZvPtEtaCent->SetBinEdges(0,fBinsZv);
260 fMCGenZvPtEtaCent->SetBinEdges(1,fBinsPt);
261 fMCGenZvPtEtaCent->SetBinEdges(2,fBinsEta);
262 fMCGenZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
263 fMCGenZvPtEtaCent->GetAxis(0)->SetTitle("MC Zv (cm)");
264 fMCGenZvPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
265 fMCGenZvPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
266 fMCGenZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
267 fMCGenZvPtEtaCent->Sumw2();
268
269 fMCRecSecZvPtEtaCent = new THnSparseF("fMCRecSecZvPtEtaCent","mcZv:mcPt:mcEta:Centrality",4,binsZvPtEtaCent);
270 fMCRecSecZvPtEtaCent->SetBinEdges(0,fBinsZv);
271 fMCRecSecZvPtEtaCent->SetBinEdges(1,fBinsPt);
272 fMCRecSecZvPtEtaCent->SetBinEdges(2,fBinsEta);
273 fMCRecSecZvPtEtaCent->SetBinEdges(3,fBinsCentrality);
274 fMCRecSecZvPtEtaCent->GetAxis(0)->SetTitle("MC Sec Zv (cm)");
275 fMCRecSecZvPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
276 fMCRecSecZvPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
277 fMCRecSecZvPtEtaCent->GetAxis(3)->SetTitle("Centrality");
278 fMCRecSecZvPtEtaCent->Sumw2();
279
280 fMCRecPrimPhiPtEtaCent = new THnSparseF("fMCRecPrimPhiPtEtaCent","mcPhi:mcPt:mcEta:Centrality",4,binsPhiPtEtaCent);
281 fMCRecPrimPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
282 fMCRecPrimPhiPtEtaCent->SetBinEdges(1,fBinsPt);
283 fMCRecPrimPhiPtEtaCent->SetBinEdges(2,fBinsEta);
284 fMCRecPrimPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
285 fMCRecPrimPhiPtEtaCent->GetAxis(0)->SetTitle("MC Phi (cm)");
286 fMCRecPrimPhiPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
287 fMCRecPrimPhiPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
288 fMCRecPrimPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
289 fMCRecPrimPhiPtEtaCent->Sumw2();
290
291 fMCGenPhiPtEtaCent = new THnSparseF("fMCGenPhiPtEtaCent","mcPhi:mcPt:mcEta:Centrality",4,binsPhiPtEtaCent);
292 fMCGenPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
293 fMCGenPhiPtEtaCent->SetBinEdges(1,fBinsPt);
294 fMCGenPhiPtEtaCent->SetBinEdges(2,fBinsEta);
295 fMCGenPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
296 fMCGenPhiPtEtaCent->GetAxis(0)->SetTitle("MC Phi (cm)");
297 fMCGenPhiPtEtaCent->GetAxis(1)->SetTitle("MC Pt (GeV/c)");
298 fMCGenPhiPtEtaCent->GetAxis(2)->SetTitle("MC Eta");
299 fMCGenPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
300 fMCGenPhiPtEtaCent->Sumw2();
301
302 fMCRecSecPhiPtEtaCent = new THnSparseF("fMCRecSecPhiPtEtaCent","mcPhi:mcPt:mcEta:Centrality",4,binsPhiPtEtaCent);
303 fMCRecSecPhiPtEtaCent->SetBinEdges(0,fBinsPhi);
304 fMCRecSecPhiPtEtaCent->SetBinEdges(1,fBinsPt);
305 fMCRecSecPhiPtEtaCent->SetBinEdges(2,fBinsEta);
306 fMCRecSecPhiPtEtaCent->SetBinEdges(3,fBinsCentrality);
307 fMCRecSecPhiPtEtaCent->GetAxis(0)->SetTitle("MC Sec Phi (cm)");
308 fMCRecSecPhiPtEtaCent->GetAxis(1)->SetTitle("MC Sec Pt (GeV/c)");
309 fMCRecSecPhiPtEtaCent->GetAxis(2)->SetTitle("MC Sec Eta");
310 fMCRecSecPhiPtEtaCent->GetAxis(3)->SetTitle("Centrality");
311 fMCRecSecPhiPtEtaCent->Sumw2();
312
313 fPt = new TH1F("fPt","fPt",2000,0,200);
314 fPt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
315 fPt->GetYaxis()->SetTitle("dN/dp_{T}");
316 fPt->Sumw2();
317
318 fMCPt = new TH1F("fMCPt","fMCPt",2000,0,200);
319 fMCPt->GetXaxis()->SetTitle("MC p_{T} (GeV/c)");
320 fMCPt->GetYaxis()->SetTitle("dN/dp_{T}");
321 fMCPt->Sumw2();
322
323 fEventStatistics = new TH1F("fEventStatistics","fEventStatistics",10,0,10);
324 fEventStatistics->GetYaxis()->SetTitle("number of events");
325 fEventStatistics->SetBit(TH1::kCanRebin);
326
327 fEventStatisticsCentrality = new TH1F("fEventStatisticsCentrality","fEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
328 fEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
329
330 fMCEventStatisticsCentrality = new TH1F("fMCEventStatisticsCentrality","fMCEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
331 fMCEventStatisticsCentrality->GetYaxis()->SetTitle("number of MC events");
332
333 fAllEventStatisticsCentrality = new TH1F("fAllEventStatisticsCentrality","fAllEventStatisticsCentrality",fCentralityNbins-1, fBinsCentrality);
334 fAllEventStatisticsCentrality->GetYaxis()->SetTitle("number of events");
335
336 fEventStatisticsCentralityTrigger = new TH2F("fEventStatisticsCentralityTrigger","fEventStatisticsCentralityTrigger;centrality;trigger",100,0,100,3,0,3);
337 fEventStatisticsCentralityTrigger->Sumw2();
338
339 fZvMultCent = new THnSparseF("fZvMultCent","Zv:mult:Centrality",3,binsZvMultCent);
340 fZvMultCent->SetBinEdges(0,fBinsZv);
341 fZvMultCent->SetBinEdges(1,fBinsMult);
342 fZvMultCent->SetBinEdges(2,fBinsCentrality);
343 fZvMultCent->GetAxis(0)->SetTitle("Zv (cm)");
344 fZvMultCent->GetAxis(1)->SetTitle("N_{acc}");
345 fZvMultCent->GetAxis(2)->SetTitle("Centrality");
346 fZvMultCent->Sumw2();
347
348 fTriggerStatistics = new TH1F("fTriggerStatistics","fTriggerStatistics",10,0,10);
349 fTriggerStatistics->GetYaxis()->SetTitle("number of events");
350
351 fMCTrackPdgCode = new TH1F("fMCTrackPdgCode","fMCTrackPdgCode",100,0,10);
352 fMCTrackPdgCode->GetYaxis()->SetTitle("number of tracks");
353 fMCTrackPdgCode->SetBit(TH1::kCanRebin);
354
355 fMCTrackStatusCode = new TH1F("fMCTrackStatusCode","fMCTrackStatusCode",100,0,10);
356 fMCTrackStatusCode->GetYaxis()->SetTitle("number of tracks");
357 fMCTrackStatusCode->SetBit(TH1::kCanRebin);
358
359 fCharge = new TH1F("fCharge","fCharge",30, -5, 5);
360 fCharge->GetXaxis()->SetTitle("Charge");
361 fCharge->GetYaxis()->SetTitle("number of tracks");
362
363 fMCCharge = new TH1F("fMCCharge","fMCCharge",30, -5, 5);
364 fMCCharge->GetXaxis()->SetTitle("MC Charge");
365 fMCCharge->GetYaxis()->SetTitle("number of tracks");
366
367 fMCPdgPt = new TH2F("fMCPdgPt","fMCPdgPt",fPtNbins-1, fBinsPt, 100,0,100);
368 fMCPdgPt->GetYaxis()->SetTitle("particle");
369 fMCPdgPt->GetXaxis()->SetTitle("Pt (GeV/c)");
370
371 fMCHijingPrim = new TH1F("fMCHijingPrim","fMCHijingPrim",2,0,2);
372 fMCHijingPrim->GetYaxis()->SetTitle("number of particles");
373
374
375
376 Int_t binsDCAxyDCAzPtEtaPhi[6] = { 200,200, fPtCheckNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
377 Double_t minDCAxyDCAzPtEtaPhi[6] = { -5, -5, 0, -1.5, 0., 0, };
378 Double_t maxDCAxyDCAzPtEtaPhi[6] = { 5., 5., 100, 1.5, 2.*TMath::Pi(), 100};
379
380 fDCAPtAll = new THnSparseF("fDCAPtAll","fDCAPtAll",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
381 fDCAPtAccepted = new THnSparseF("fDCAPtAccepted","fDCAPtAccepted",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
382 fMCDCAPtSecondary = new THnSparseF("fMCDCAPtSecondary","fMCDCAPtSecondary",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
383 fMCDCAPtPrimary = new THnSparseF("fMCDCAPtPrimary","fMCDCAPtPrimary",6, binsDCAxyDCAzPtEtaPhi, minDCAxyDCAzPtEtaPhi, maxDCAxyDCAzPtEtaPhi);
384
385 fDCAPtAll->SetBinEdges(2, fBinsPtCheck);
386 fDCAPtAccepted->SetBinEdges(2, fBinsPtCheck);
387 fMCDCAPtSecondary->SetBinEdges(2, fBinsPtCheck);
388 fMCDCAPtPrimary->SetBinEdges(2, fBinsPtCheck);
389
390 fDCAPtAll->SetBinEdges(3, fBinsEtaCheck);
391 fDCAPtAccepted->SetBinEdges(3, fBinsEtaCheck);
392 fMCDCAPtSecondary->SetBinEdges(3, fBinsEtaCheck);
393 fMCDCAPtPrimary->SetBinEdges(3, fBinsEtaCheck);
394
395 fDCAPtAll->SetBinEdges(5, fBinsCentrality);
396 fDCAPtAccepted->SetBinEdges(5, fBinsCentrality);
397 fMCDCAPtSecondary->SetBinEdges(5, fBinsCentrality);
398 fMCDCAPtPrimary->SetBinEdges(5, fBinsCentrality);
399
400 fDCAPtAll->Sumw2();
401 fDCAPtAccepted->Sumw2();
402 fMCDCAPtSecondary->Sumw2();
403 fMCDCAPtPrimary->Sumw2();
404
405 fDCAPtAll->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
406 fDCAPtAll->GetAxis(1)->SetTitle("DCA_{z} (cm)");
407 fDCAPtAll->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
408 fDCAPtAll->GetAxis(3)->SetTitle("#eta");
409 fDCAPtAll->GetAxis(4)->SetTitle("#phi");
410 fDCAPtAll->GetAxis(5)->SetTitle("Centrality");
411
412 fDCAPtAccepted->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
413 fDCAPtAccepted->GetAxis(1)->SetTitle("DCA_{z} (cm)");
414 fDCAPtAccepted->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
415 fDCAPtAccepted->GetAxis(3)->SetTitle("#eta");
416 fDCAPtAccepted->GetAxis(4)->SetTitle("#phi");
417 fDCAPtAccepted->GetAxis(5)->SetTitle("Centrality");
418
419 fMCDCAPtSecondary->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
420 fMCDCAPtSecondary->GetAxis(1)->SetTitle("DCA_{z} (cm)");
421 fMCDCAPtSecondary->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
422 fMCDCAPtSecondary->GetAxis(3)->SetTitle("#eta");
423 fMCDCAPtSecondary->GetAxis(4)->SetTitle("#phi");
424 fMCDCAPtSecondary->GetAxis(5)->SetTitle("Centrality");
425
426 fMCDCAPtPrimary->GetAxis(0)->SetTitle("DCA_{xy} (cm)");
427 fMCDCAPtPrimary->GetAxis(1)->SetTitle("DCA_{z} (cm)");
428 fMCDCAPtPrimary->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
429 fMCDCAPtPrimary->GetAxis(3)->SetTitle("#eta");
430 fMCDCAPtPrimary->GetAxis(4)->SetTitle("#phi");
431 fMCDCAPtPrimary->GetAxis(5)->SetTitle("Centrality");
432
433
434 char cFullTempTitle[255];
435 char cTempTitleAxis0All[255];
436 char cTempTitleAxis0Acc[255];
437 // char cTempTitleAxis1[255];
438 char cFullTempName[255];
439 char cTempNameAxis0[255];
440 // char cTempNameAxis1[255];
441 const Int_t iNbinRowsClusters = 21;
442 // Double_t dBinsRowsClusters[iNbinRowsClusters] = {0, 7.95, 15.9, 23.85, 31.8, 39.75, 47.7, 55.65, 63.6, 71.55, 79.5, 87.45, 95.4, 103.35, 111.3, 119.25, 127.2, 135.15, 143.1, 151.05, 159.};
443
444 const Int_t iNbinChi = 51;
445 const Int_t iNbinLength = 165;
446 // Double_t dBinsChi[iNbinChi] = {0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.2, 3.4, 3.6, 3.8, 4, 4.2, 4.4, 4.6, 4.8, 5, 5.2, 5.4, 5.6, 5.8, 6, 6.2, 6.4, 6.6, 6.8, 7, 7.2, 7.4, 7.6, 7.8, 8, 8.2, 8.4, 8.6, 8.8, 9, 9.2, 9.4, 9.6, 9.8,10.};
447
448 Int_t iNbin = 0;
449 // Double_t *dBins = 0x0;
450 Double_t dBinMin = 0;
451 Double_t dBinMax = 0;
452
453 for(Int_t iCheckQuant = 0; iCheckQuant < cqMax; iCheckQuant++)
454 {
455 // iCheckQuant: 0 = CrossedRows, 1 = Nclusters, 2 = Chi^2/clusterTPC
456 if(iCheckQuant == cqCrossedRows)
457 {
458 snprintf(cTempTitleAxis0All,255, "NcrossedRows before Cut");
459 snprintf(cTempTitleAxis0Acc,255, "NcrossedRows after Cut");
460 snprintf(cTempNameAxis0,255, "CrossedRows");
461 iNbin = iNbinRowsClusters;
462 dBinMin = 0;
463 dBinMax = 159.;
464 }
465 else if(iCheckQuant == cqNcluster)
466 {
467 snprintf(cTempTitleAxis0All,255, "Nclusters before Cut");
468 snprintf(cTempTitleAxis0Acc,255, "Nclusters after Cut");
469 snprintf(cTempNameAxis0,255, "Clusters");
470 iNbin = iNbinRowsClusters;
471 dBinMin = 0;
472 dBinMax = 159.;
473 }
474 else if(iCheckQuant == cqChi)
475 {
476 snprintf(cTempTitleAxis0All,255, "#Chi^{2}/cluster before Cut");
477 snprintf(cTempTitleAxis0Acc,255, "#Chi^{2}/cluster after Cut");
478 snprintf(cTempNameAxis0,255, "Chi");
479 iNbin = iNbinChi;
480 dBinMin = 0;
481 dBinMax = 10.;
482 }
483 else if(iCheckQuant == cqLength)
484 {
485 snprintf(cTempTitleAxis0All,255, "Length in TPC before Cut (cm)");
486 snprintf(cTempTitleAxis0Acc,255, "Length in TPC after Cut (cm)");
487 snprintf(cTempNameAxis0,255, "Length");
488 iNbin = iNbinLength;
489 dBinMin = 0;
490 dBinMax = 165.;
491 }
492
493 Int_t binsCheckPtEtaPhi[5] = { iNbin, fPtCheckNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
494// Int_t binsCheckPtEtaPhi[5] = { iNbin, fPtNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
495 Double_t minCheckPtEtaPhi[5] = { dBinMin, 0, -1.5, 0., 0, };
496 Double_t maxCheckPtEtaPhi[5] = { dBinMax, 100, 1.5, 2.*TMath::Pi(), 100};
497
498 snprintf(cFullTempName, 255, "f%sPtEtaPhiAll",cTempNameAxis0);
499 snprintf(cFullTempTitle, 255,"%s;%s;p_{T} (GeV/c);#eta;#phi;Centrality", cFullTempName, cTempTitleAxis0All);
500 fCrossCheckAll[iCheckQuant] = new THnF(cFullTempName, cFullTempTitle, 5, binsCheckPtEtaPhi, minCheckPtEtaPhi, maxCheckPtEtaPhi);
501 fCrossCheckAll[iCheckQuant]->SetBinEdges(1, fBinsPtCheck);
502 fCrossCheckAll[iCheckQuant]->SetBinEdges(2, fBinsEtaCheck);
503 fCrossCheckAll[iCheckQuant]->Sumw2();
504
505 snprintf(cFullTempName, 255, "f%sPtEtaPhiAcc",cTempNameAxis0);
506 snprintf(cFullTempTitle, 255,"%s;%s;p_{T} (GeV/c);#eta;#phi;Centrality", cFullTempName, cTempTitleAxis0Acc);
507 fCrossCheckAcc[iCheckQuant] = new THnF(cFullTempName, cFullTempTitle, 5, binsCheckPtEtaPhi, minCheckPtEtaPhi, maxCheckPtEtaPhi);
508 fCrossCheckAcc[iCheckQuant]->SetBinEdges(1, fBinsPtCheck);
509 fCrossCheckAcc[iCheckQuant]->SetBinEdges(2, fBinsEtaCheck);
510 fCrossCheckAcc[iCheckQuant]->Sumw2();
511 } // end iCheckQuant
512
513 fCutPercClusters = new TH1F("fCutPercClusters","fCutPercClusters;NclustersTPC;counts",160,0,160);
514 fCutPercClusters->Sumw2();
515 fCutPercCrossed = new TH1F("fCutPercCrossed","fCutPercCrossed;NcrossedRowsTPC;counts",160,0,160);
516 fCutPercCrossed->Sumw2();
517
518 fCrossCheckRowsLength = new TH2F("fCrossCheckRowsLength","fCrossCheckRowsLength;Length in TPC;NcrossedRows",170,0,170,170,0,170);
519 fCrossCheckRowsLength->Sumw2();
520
521 fCrossCheckClusterLength = new TH2F("fCrossCheckClusterLength","fCrossCheckClusterLength;Length in TPC;NClusters",170,0,170,170,0,170);
522 fCrossCheckClusterLength->Sumw2();
523
524 fCrossCheckRowsLengthAcc = new TH2F("fCrossCheckRowsLengthAcc","fCrossCheckRowsLengthAcc;Length in TPC;NcrossedRows",170,0,170,170,0,170);
525 fCrossCheckRowsLengthAcc->Sumw2();
526
527 fCrossCheckClusterLengthAcc = new TH2F("fCrossCheckClusterLengthAcc","fCrossCheckClusterLengthAcc;Length in TPC;NClusters",170,0,170,170,0,170);
528 fCrossCheckClusterLengthAcc->Sumw2();
529
530
531 // Add Histos, Profiles etc to List
532 fOutputList->Add(fZvPtEtaCent);
533 fOutputList->Add(fPhiPtEtaCent);
534 fOutputList->Add(fPt);
535 fOutputList->Add(fMCRecPrimZvPtEtaCent);
536 fOutputList->Add(fMCGenZvPtEtaCent);
537 fOutputList->Add(fMCRecSecZvPtEtaCent);
538 fOutputList->Add(fMCRecPrimPhiPtEtaCent);
539 fOutputList->Add(fMCGenPhiPtEtaCent);
540 fOutputList->Add(fMCRecSecPhiPtEtaCent);
541 fOutputList->Add(fMCPt);
542 fOutputList->Add(fEventStatistics);
543 fOutputList->Add(fEventStatisticsCentrality);
544 fOutputList->Add(fMCEventStatisticsCentrality);
545 fOutputList->Add(fAllEventStatisticsCentrality);
546 fOutputList->Add(fEventStatisticsCentralityTrigger);
547 fOutputList->Add(fZvMultCent);
548 fOutputList->Add(fTriggerStatistics);
549 fOutputList->Add(fMCTrackPdgCode);
550 fOutputList->Add(fMCTrackStatusCode);
551 fOutputList->Add(fCharge);
552 fOutputList->Add(fMCCharge);
553 fOutputList->Add(fMCPdgPt);
554 fOutputList->Add(fMCHijingPrim);
555 fOutputList->Add(fDCAPtAll);
556 fOutputList->Add(fDCAPtAccepted);
557 fOutputList->Add(fMCDCAPtSecondary);
558 fOutputList->Add(fMCDCAPtPrimary);
559 for(Int_t i = 0; i < cqMax; i++)
560 {
561 fOutputList->Add(fCrossCheckAll[i]);
562 fOutputList->Add(fCrossCheckAcc[i]);
563 }
564 fOutputList->Add(fCutPercClusters);
565 fOutputList->Add(fCutPercCrossed);
566 fOutputList->Add(fCrossCheckRowsLength);
567 fOutputList->Add(fCrossCheckClusterLength);
568 fOutputList->Add(fCrossCheckRowsLengthAcc);
569 fOutputList->Add(fCrossCheckClusterLengthAcc);
570
571 PostData(1, fOutputList);
572}
573
574void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
575{
576
577 // Main Loop
578 // called for each event
579 fEventStatistics->Fill("all events",1);
580
581 // set ZERO pointers:
582 AliInputEventHandler *inputHandler = NULL;
583 AliAODTrack *track = NULL;
584 AliAODMCParticle *mcPart = NULL;
585 AliAODMCHeader *mcHdr = NULL;
586 AliGenHijingEventHeader *genHijingHeader = NULL;
587 //AliGenPythiaEventHeader *genPythiaHeader = NULL;
588
589 Bool_t bIsEventSelectedMB = kFALSE;
590 Bool_t bIsEventSelectedSemi = kFALSE;
591 Bool_t bIsEventSelectedCentral = kFALSE;
592 Bool_t bIsEventSelected = kFALSE;
593 Bool_t bIsPrimary = kFALSE;
594 Bool_t bIsHijingParticle = kFALSE;
595 Bool_t bMotherIsHijingParticle = kFALSE;
596 //Bool_t bIsPythiaParticle = kFALSE;
597 Bool_t bEventHasATrack = kFALSE;
598 Bool_t bEventHasATrackInRange = kFALSE;
599 Int_t nTriggerFired = 0;
600
601
602 Double_t dMCTrackZvPtEtaCent[4] = {0};
603 Double_t dTrackZvPtEtaCent[4] = {0};
604
605 Double_t dMCTrackPhiPtEtaCent[4] = {0};
606 Double_t dTrackPhiPtEtaCent[4] = {0};
607
608 Double_t dDCA[2] = {0};
609
610 Double_t dMCEventZv = -100;
611 Double_t dEventZv = -100;
612 Int_t iAcceptedMultiplicity = 0;
613
614 fIsMonteCarlo = kFALSE;
615
616 AliAODEvent *eventAOD = 0x0;
617 eventAOD = dynamic_cast<AliAODEvent*>( InputEvent() );
618 if (!eventAOD) {
619 AliWarning("ERROR: eventAOD not available \n");
620 return;
621 }
622
623 // check, which trigger has been fired
624 inputHandler = (AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
625 bIsEventSelectedMB = ( inputHandler->IsEventSelected() & AliVEvent::kMB);
626 bIsEventSelectedSemi = ( inputHandler->IsEventSelected() & AliVEvent::kSemiCentral);
627 bIsEventSelectedCentral = ( inputHandler->IsEventSelected() & AliVEvent::kCentral);
628
629 if(bIsEventSelectedMB || bIsEventSelectedSemi || bIsEventSelectedCentral) fTriggerStatistics->Fill("all triggered events",1);
630 if(bIsEventSelectedMB) { fTriggerStatistics->Fill("MB trigger",1); nTriggerFired++; }
631 if(bIsEventSelectedSemi) { fTriggerStatistics->Fill("SemiCentral trigger",1); nTriggerFired++; }
632 if(bIsEventSelectedCentral) { fTriggerStatistics->Fill("Central trigger",1); nTriggerFired++; }
633 if(nTriggerFired == 0) { fTriggerStatistics->Fill("No trigger",1); }
634
635 bIsEventSelected = ( inputHandler->IsEventSelected() & GetCollisionCandidates() );
636
637 // only take tracks of events, which are triggered
638 if(nTriggerFired == 0) { return; }
639
640 // if( !bIsEventSelected || nTriggerFired>1 ) return;
641
642 // fEventStatistics->Fill("events with only coll. cand.", 1);
643
644
645
646 // check if there is a stack, if yes, then do MC loop
647 TList *list = eventAOD->GetList();
648 TClonesArray *stack = 0x0;
649 stack = (TClonesArray*)list->FindObject(AliAODMCParticle::StdBranchName());
650
651 if( stack )
652 {
653 fIsMonteCarlo = kTRUE;
654
655 mcHdr = (AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
656
657 genHijingHeader = GetHijingEventHeader(mcHdr);
658 // genPythiaHeader = GetPythiaEventHeader(mcHdr);
659
660 if(!genHijingHeader) { return; }
661
662 // if(!genPythiaHeader) { return; }
663
664 dMCEventZv = mcHdr->GetVtxZ();
665 dMCTrackZvPtEtaCent[0] = dMCEventZv;
666 fEventStatistics->Fill("MC all events",1);
667 }
668
669 AliCentrality* aCentrality = eventAOD->GetCentrality();
670 Double_t dCentrality = aCentrality->GetCentralityPercentile("V0M");
671
672 if( dCentrality < 0 ) return;
673 fEventStatistics->Fill("after centrality selection",1);
674
675
676
677 // start with MC truth analysis
678 if(fIsMonteCarlo)
679 {
680
681 if( dMCEventZv > GetCutMaxZVertex() ) { return; }
682
683 dMCTrackZvPtEtaCent[0] = dMCEventZv;
684
685 fEventStatistics->Fill("MC afterZv cut",1);
686
687 for(Int_t iMCtrack = 0; iMCtrack < stack->GetEntriesFast(); iMCtrack++)
688 {
689 mcPart =(AliAODMCParticle*)stack->At(iMCtrack);
690
691 // check for charge
692 if( !(IsMCTrackAccepted(mcPart)) ) continue;
693
694 if(!IsHijingParticle(mcPart, genHijingHeader)) { continue; }
695
696 if(mcPart->IsPhysicalPrimary() )
697 {
698 fMCHijingPrim->Fill("IsPhysicalPrimary",1);
699 }
700 else
701 {
702 fMCHijingPrim->Fill("NOT a primary",1);
703 continue;
704 }
705
706
707 //
708 // ======================== fill histograms ========================
709 dMCTrackZvPtEtaCent[1] = mcPart->Pt();
710 dMCTrackZvPtEtaCent[2] = mcPart->Eta();
711 dMCTrackZvPtEtaCent[3] = dCentrality;
712 fMCGenZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
713
714 dMCTrackPhiPtEtaCent[0] = mcPart->Phi();
715 dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
716 dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
717 dMCTrackPhiPtEtaCent[3] = dCentrality;
718 fMCGenPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
719
720 bEventHasATrack = kTRUE;
721
722
723 if( (dMCTrackZvPtEtaCent[1] > GetCutPtMin() ) &&
724 (dMCTrackZvPtEtaCent[1] < GetCutPtMax() ) &&
725 (dMCTrackZvPtEtaCent[2] > GetCutEtaMin() ) &&
726 (dMCTrackZvPtEtaCent[2] < GetCutEtaMax() ) )
727 {
728 fMCPt->Fill(mcPart->Pt());
729 fMCCharge->Fill(mcPart->Charge()/3.);
730 bEventHasATrackInRange = kTRUE;
731 }
732
733 }
734 } // isMonteCarlo
735
736 if(bEventHasATrack) { fEventStatistics->Fill("MC events with tracks",1); }
737 if(bEventHasATrackInRange)
738 {
739 fEventStatistics->Fill("MC events with tracks in range",1);
740 fMCEventStatisticsCentrality->Fill(dCentrality);
741 }
742 bEventHasATrack = kFALSE;
743 bEventHasATrackInRange = kFALSE;
744
745
746
747 // Loop over recontructed tracks
748
749 dEventZv = eventAOD->GetPrimaryVertex()->GetZ();
750 if( TMath::Abs(dEventZv) > GetCutMaxZVertex() ) return;
751
752 fAllEventStatisticsCentrality->Fill(dCentrality/*, nTriggerFired*/);
753
754 fEventStatistics->Fill("after Zv cut",1);
755
756 dTrackZvPtEtaCent[0] = dEventZv;
757
758 if(AreRelativeCutsEnabled())
759 {
760 if(!SetRelativeCuts(eventAOD)) return;
761 }
762
763 for(Int_t itrack = 0; itrack < eventAOD->GetNumberOfTracks(); itrack++)
764 {
765 track = eventAOD->GetTrack(itrack);
766 if(!track) continue;
767
768 mcPart = NULL;
769 dMCTrackZvPtEtaCent[1] = 0;
770 dMCTrackZvPtEtaCent[2] = 0;
771 dMCTrackZvPtEtaCent[3] = 0;
772
773 dMCTrackPhiPtEtaCent[0] = 0;
774 dMCTrackPhiPtEtaCent[1] = 0;
775 dMCTrackPhiPtEtaCent[2] = 0;
776 dMCTrackPhiPtEtaCent[3] = 0;
777
778 bIsPrimary = kFALSE;
779
780 GetDCA(track, eventAOD, dDCA);
781
782 Double_t dDCAxyDCAzPt[5] = { dDCA[0], dDCA[1], track->Pt(), track->Eta(), track->Phi() };
783
784 fDCAPtAll->Fill(dDCAxyDCAzPt);
785
786 if( !(IsTrackAccepted(track, dCentrality, eventAOD->GetMagneticField())) ) continue;
787
788 dTrackZvPtEtaCent[1] = track->Pt();
789 dTrackZvPtEtaCent[2] = track->Eta();
790 dTrackZvPtEtaCent[3] = dCentrality;
791
792 dTrackPhiPtEtaCent[0] = track->Phi();
793 dTrackPhiPtEtaCent[1] = track->Pt();
794 dTrackPhiPtEtaCent[2] = track->Eta();
795 dTrackPhiPtEtaCent[3] = dCentrality;
796
797 if( fIsMonteCarlo )
798 {
799 mcPart = (AliAODMCParticle*)stack->At(TMath::Abs(track->GetLabel()));
800 if( !mcPart ) { continue; }
801
802 // check for charge
803 // if( !(IsMCTrackAccepted(mcPart)) ) { continue; }
804
805 bIsHijingParticle = IsHijingParticle(mcPart, genHijingHeader);
806 // bIsPythiaParticle = IsPythiaParticle(mcPart, genPythiaHeader);
807
808 // if(!bIsHijingParticle) continue; // only take real tracks, not injected ones
809
810 bIsPrimary = mcPart->IsPhysicalPrimary();
811
812 dMCTrackZvPtEtaCent[1] = mcPart->Pt();
813 dMCTrackZvPtEtaCent[2] = mcPart->Eta();
814 dMCTrackZvPtEtaCent[3] = dCentrality;
815
816 dMCTrackPhiPtEtaCent[0] = mcPart->Phi();
817 dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
818 dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
819 dMCTrackPhiPtEtaCent[3] = dCentrality;
820
821 if(bIsPrimary && bIsHijingParticle)
822 {
823 fMCRecPrimZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
824 fMCRecPrimPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
825 if( (TMath::Abs(mcPart->Eta()) < 0.8) && (dCentrality<5.) ) { fMCPdgPt->Fill(mcPart->Pt(), Form("%s",GetParticleName(mcPart->GetPdgCode())), 1); }
826 fMCDCAPtPrimary->Fill(dDCAxyDCAzPt);
827 }
828
829 if(!bIsPrimary /*&& !bIsHijingParticle*/)
830 {
831 Int_t indexMoth = mcPart->GetMother();
832 if(indexMoth >= 0)
833 {
834 AliAODMCParticle* moth = (AliAODMCParticle*)stack->At(indexMoth);
835 bMotherIsHijingParticle = IsHijingParticle(moth, genHijingHeader);
836
837 if(bMotherIsHijingParticle) // only store secondaries, which come from a not embedded signal!
838 {
839 fMCTrackStatusCode->Fill(Form("%d",mcPart->GetStatus()), 1);
840
841
842 fMCRecSecZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
843 fMCRecSecPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
844 fMCDCAPtSecondary->Fill(dDCAxyDCAzPt);
845 fMCTrackPdgCode->Fill(Form("%s_H%i_H%i",GetParticleName(moth->GetPdgCode()),bMotherIsHijingParticle, bIsHijingParticle), 1);
846 // delete moth;
847 }
848 }
849 }
850 } // end isMonteCarlo
851
852 // ======================== fill histograms ========================
853
854 // only keep prim and sec from not embedded signal
855 Bool_t bKeepMCTrack = kFALSE;
856 if(fIsMonteCarlo)
857 {
858 if( (bIsHijingParticle && bIsPrimary) ^ (bMotherIsHijingParticle && !bIsPrimary) )
859 {
860 bKeepMCTrack = kTRUE;
861 }
862 else
863 {
864 continue;
865 }
866 }
867
868 bEventHasATrack = kTRUE;
869
870 fZvPtEtaCent->Fill(dTrackZvPtEtaCent);
871 fPhiPtEtaCent->Fill(dTrackPhiPtEtaCent);
872
873 fDCAPtAccepted->Fill(dDCAxyDCAzPt);
874
875 if( (dTrackZvPtEtaCent[1] > GetCutPtMin()) &&
876 (dTrackZvPtEtaCent[1] < GetCutPtMax()) &&
877 (dTrackZvPtEtaCent[2] > GetCutEtaMin()) &&
878 (dTrackZvPtEtaCent[2] < GetCutEtaMax()) )
879 {
880 iAcceptedMultiplicity++;
881 bEventHasATrackInRange = kTRUE;
882 fPt->Fill(track->Pt());
883 fCharge->Fill(track->Charge());
884 }
885 } // end track loop
886
887 if(bEventHasATrack) { fEventStatistics->Fill("events with tracks",1); bEventHasATrack = kFALSE; }
888
889 if(bEventHasATrackInRange)
890 {
891 fEventStatistics->Fill("events with tracks in range",1);
892 fEventStatisticsCentrality->Fill(dCentrality);
893 bEventHasATrackInRange = kFALSE;
894
895 if(bIsEventSelectedMB) fEventStatisticsCentralityTrigger->Fill(dCentrality, 0);
896 if(bIsEventSelectedSemi) fEventStatisticsCentralityTrigger->Fill(dCentrality, 1);
897 if(bIsEventSelectedCentral) fEventStatisticsCentralityTrigger->Fill(dCentrality, 2);
898 }
899
900 Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
901 fZvMultCent->Fill(dEventZvMultCent);
902
903 PostData(1, fOutputList);
904
905}
906
907Bool_t AlidNdPtAnalysisPbPbAOD::SetRelativeCuts(AliAODEvent *event)
908{
909 if(!event) return kFALSE;
910
911 AliAODTrack *tr = 0x0;
912 TH1F *hCluster = new TH1F("hCluster","hCluster",160,0,160);
913 TH1F *hCrossed = new TH1F("hCrossed","hCrossed",160,0,160);
914
915 for(Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++)
916 {
917 tr = event->GetTrack(itrack);
918 if(!tr) continue;
919
920 // do some selection already
921 //if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { continue; }
922
923 Double_t dNClustersTPC = tr->GetTPCNcls();
924 Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
925
926 hCluster->Fill(dNClustersTPC);
927 hCrossed->Fill(dCrossedRowsTPC);
928 }
929
930 // loop trough histogram to check, where percentage is reach
931 Double_t dTotIntCluster = hCluster->Integral();
932 Double_t dTotIntCrossed = hCrossed->Integral();
933 Float_t dIntCluster = 0;
934 Float_t dIntCrossed = 0;
935
936 if(dTotIntCluster)
937 {
938 for(Int_t i = 0; i < hCluster->GetNbinsX(); i++)
939 {
940 if(hCluster->GetBinCenter(i) < 0) continue;
941 dIntCluster += hCluster->GetBinContent(i);
942 if(dIntCluster/dTotIntCluster > (1-GetCutPercMinNClustersTPC()))
943 {
944 SetCutMinNClustersTPC(hCluster->GetBinCenter(i));
945 fCutPercClusters->Fill(hCluster->GetBinCenter(i));
946 break;
947 }
948 }
949 }
950
951 if(dTotIntCrossed)
952 {
953 for(Int_t i = 0; i < hCrossed->GetNbinsX(); i++)
954 {
955 if(hCrossed->GetBinCenter(i) < 0) continue;
956 dIntCrossed += hCrossed->GetBinContent(i);
957 if(dIntCrossed/dTotIntCrossed > (1-GetCutPercMinNCrossedRowsTPC()))
958 {
959 SetCutMinNClustersTPC(hCrossed->GetBinCenter(i));
960 fCutPercCrossed->Fill(hCrossed->GetBinCenter(i));
961 break;
962 }
963 }
964 }
965
966 delete hCrossed;
967 delete hCluster;
968 return kTRUE;
969
970}
971
972Bool_t AlidNdPtAnalysisPbPbAOD::IsTrackAccepted(AliAODTrack *tr, Double_t dCentrality, Double_t bMagZ)
973{
974 if(!tr) return kFALSE;
975
976 if(tr->Charge()==0) { return kFALSE; }
977
978 //
979 // as done in AliAnalysisTaskFragmentationFunction
980 //
981
982 Short_t sign = tr->Charge();
983 Double_t xyz[50];
984 Double_t pxpypz[50];
985 Double_t cv[100];
986
987 for(Int_t i = 0; i < 100; i++) cv[i] = 0;
988 for(Int_t i = 0; i < 50; i++) xyz[i] = 0;
989 for(Int_t i = 0; i < 50; i++) pxpypz[i] = 0;
990
991 tr->GetXYZ(xyz);
992 tr->GetPxPyPz(pxpypz);
993
994 // similar error occured as this one:
995 // See https://savannah.cern.ch/bugs/?102721
996 // which is one of the two 11h re-filtering follow-ups:
997 // Andrea Dainese now first does the beam pipe
998 // check and then copies from the vtrack (was the other
999 // way around) to avoid the crash in the etp::Set()
1000
1001// if(xyz[0]*xyz[0]+xyz[1]*xyz[1] > 3.*3.) { return kFALSE; }
1002
1003 AliExternalTrackParam * par = new AliExternalTrackParam(xyz, pxpypz, cv, sign);
1004 if(!par) { return kFALSE; }
1005 AliESDtrack dummy;
1006// Double_t dLength = dummy.GetLengthInActiveZone(par,3,236, -5 ,0,0);
1007// Double_t dLengthInTPC = GetLengthInTPC(tr, 1.8, 220, bMagZ);
1008
1009 Double_t dLengthInTPC = dummy.GetLengthInActiveZone(par,3,236, bMagZ ,0,0);
1010 Double_t dNClustersTPC = tr->GetTPCNcls();
1011 Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
1012 Double_t dFindableClustersTPC = tr->GetTPCNclsF();
1013 Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
1014
1015 // hAllCrossedRowsTPC->Fill(dCrossedRowsTPC);
1016
1017 Double_t dCheck[cqMax] = {dCrossedRowsTPC, dNClustersTPC, dChi2PerClusterTPC, dLengthInTPC};// = new Double_t[cqMax];
1018 Double_t dKine[kqMax] = {tr->Pt(), tr->Eta(), tr->Phi()};// = new Double_t[kqMax];
1019 // dKine[0] = tr->Pt();
1020 // dKine[1] = tr->Eta();
1021 // dKine[2] = tr->Phi();
1022 //
1023 // dCheck[0] = dCrossedRowsTPC;
1024 // dCheck[1] = dNClustersTPC;
1025 // dCheck[2] = dChi2PerClusterTPC;
1026
1027
1028 FillDebugHisto(dCheck, dKine, dCentrality, kFALSE);
1029
1030 // first cut on length
1031
1032 if( DoCutLengthInTPCPtDependent() && ( dLengthInTPC < GetPrefactorLengthInTPCPtDependent()*(130-5*TMath::Abs(1./tr->Pt())) ) ) { return kFALSE; }
1033
1034 // filter bit 5
1035 if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { return kFALSE; }
1036
1037 // filter bit 4
1038 // if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) ) { return kFALSE; }
1039
1040 // hFilterCrossedRowsTPC->Fill(dCrossedRowsTPC);
1041
1042
1043 if(dFindableClustersTPC == 0) {return kFALSE; }
1044 if(dCrossedRowsTPC < GetCutMinNCrossedRowsTPC()) { return kFALSE; }
1045 if( (dCrossedRowsTPC/dFindableClustersTPC) < GetCutMinRatioCrossedRowsOverFindableClustersTPC() ) { return kFALSE; }
1046 if(dNClustersTPC < GetCutMinNClustersTPC()) { return kFALSE; }
1047
1048 if (!(tr->GetStatus() & AliVTrack::kITSrefit)) { return kFALSE; } // no ITS refit
1049
1050 // do a relativ cut in Nclusters, both time at 80% of mean
1051 // if(fIsMonteCarlo)
1052 // {
1053 // if(dNClustersTPC < 88) { return kFALSE; }
1054 // }
1055 // else
1056 // {
1057 // if(dNClustersTPC < 76) { return kFALSE; }
1058 // }
1059
1060
1061 FillDebugHisto(dCheck, dKine, dCentrality, kTRUE);
1062
1063
1064 // hAccNclsTPC->Fill(dNClustersTPC);
1065 // hAccCrossedRowsTPC->Fill(dCrossedRowsTPC);
1066 // Double_t dFindableClustersTPC = tr->GetTPCNclsF();
1067 // Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
1068 //
1069 // Bool_t bIsFromKink = kFALSE;
1070 // if(tr->GetProdVertex()->GetType() == AliAODVertex::kKink) bIsFromKink = kTRUE;
1071 //
1072 // // from AliAnalysisTaskPIDqa.cxx
1073 // ULong_t uStatus = tr->GetStatus();
1074 // Bool_t bHasRefitTPC = kFALSE;
1075 // Bool_t bHasRefitITS = kFALSE;
1076 //
1077 // if ((uStatus & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) bHasRefitTPC = kTRUE;
1078 // if ((uStatus & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) bHasRefitITS = kTRUE;
1079 //
1080 // // from AliDielectronVarManager.h
1081 // Bool_t bHasHitInSPD = kFALSE;
1082 // for (Int_t iC=0; iC<2; iC++)
1083 // {
1084 // if (((tr->GetITSClusterMap()) & (1<<(iC))) > 0) { bHasHitInSPD = kTRUE; }
1085 // }
1086 //
1087 // Double_t dNClustersITS = tr->GetITSNcls();
1088
1089 // cuts to be done:
1090 // TPC
1091 // esdTrackCuts->SetMinNCrossedRowsTPC(70);
1092 // esdTrackCuts->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);
1093 //
1094 // esdTrackCuts->SetMaxChi2PerClusterTPC(4);
1095 // esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
1096 // esdTrackCuts->SetRequireTPCRefit(kTRUE);
1097 // ITS
1098 // esdTrackCuts->SetRequireITSRefit(kTRUE);
1099 // esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1100 //
1101 // esdTrackCuts->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.1");
1102 // esdTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);
1103 //
1104 // esdTrackCuts->SetMaxDCAToVertexZ(2);
1105 // esdTrackCuts->SetDCAToVertex2D(kFALSE);
1106 // esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1107 //
1108 // esdTrackCuts->SetMaxChi2PerClusterITS(36);
1109
1110 // delete [] dKine;
1111 // delete [] dCheck;
1112 return kTRUE;
1113}
1114
1115Bool_t AlidNdPtAnalysisPbPbAOD::FillDebugHisto(Double_t *dCrossCheckVar, Double_t *dKineVar, Double_t dCentrality, Bool_t bIsAccepted)
1116{
1117 if(bIsAccepted)
1118 {
1119 for(Int_t iCrossCheck = 0; iCrossCheck < cqMax; iCrossCheck++)
1120 {
1121 Double_t dFillIt[5] = {dCrossCheckVar[iCrossCheck], dKineVar[0], dKineVar[1], dKineVar[2], dCentrality};
1122 fCrossCheckAcc[iCrossCheck]->Fill(dFillIt);
1123 }
1124
1125 fCrossCheckRowsLengthAcc->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqCrossedRows]);
1126 fCrossCheckClusterLengthAcc->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqNcluster]);
1127 }
1128 else
1129 {
1130 for(Int_t iCrossCheck = 0; iCrossCheck < cqMax; iCrossCheck++)
1131 {
1132 Double_t dFillIt[5] = {dCrossCheckVar[iCrossCheck], dKineVar[0], dKineVar[1], dKineVar[2], dCentrality};
1133 fCrossCheckAll[iCrossCheck]->Fill(dFillIt);
1134 }
1135
1136 fCrossCheckRowsLength->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqCrossedRows]);
1137 fCrossCheckClusterLength->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqNcluster]);
1138 }
1139
1140 return kTRUE;
1141
1142}
1143
1144Bool_t AlidNdPtAnalysisPbPbAOD::GetDCA(const AliAODTrack *track, AliAODEvent *evt, Double_t d0z0[2])
1145{
1146 // function adapted from AliDielectronVarManager.h
1147
1148 if(track->TestBit(AliAODTrack::kIsDCA)){
1149 d0z0[0]=track->DCA();
1150 d0z0[1]=track->ZAtDCA();
1151 return kTRUE;
1152 }
1153
1154 Bool_t ok=kFALSE;
1155 if(evt) {
1156 Double_t covd0z0[3];
1157 //AliAODTrack copy(*track);
1158 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1159
1160 Float_t xstart = etp.GetX();
1161 if(xstart>3.) {
1162 d0z0[0]=-999.;
1163 d0z0[1]=-999.;
1164 //printf("This method can be used only for propagation inside the beam pipe \n");
1165 return kFALSE;
1166 }
1167
1168
1169 AliAODVertex *vtx =(AliAODVertex*)(evt->GetPrimaryVertex());
1170 Double_t fBzkG = evt->GetMagneticField(); // z componenent of field in kG
1171 ok = etp.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
1172 //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
1173 }
1174 if(!ok){
1175 d0z0[0]=-999.;
1176 d0z0[1]=-999.;
1177 }
1178 return ok;
1179}
1180
1181
1182Bool_t AlidNdPtAnalysisPbPbAOD::IsMCTrackAccepted(AliAODMCParticle *part)
1183{
1184 if(!part) return kFALSE;
1185
1186 Double_t charge = part->Charge()/3.;
1187 if (TMath::Abs(charge) < 0.001) return kFALSE;
1188
1189 return kTRUE;
1190}
1191
1192const char * AlidNdPtAnalysisPbPbAOD::GetParticleName(Int_t pdg)
1193{
1194 TParticlePDG * p1 = TDatabasePDG::Instance()->GetParticle(pdg);
1195 if(p1) return p1->GetName();
1196 return Form("%d", pdg);
1197}
1198
1199AliGenHijingEventHeader* AlidNdPtAnalysisPbPbAOD::GetHijingEventHeader(AliAODMCHeader *header)
1200{
1201 //
1202 // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
1203 //
1204
1205 if(!header) return 0x0;
1206 AliGenHijingEventHeader* hijingGenHeader = NULL;
1207
1208 TList* headerList = header->GetCocktailHeaders();
1209
1210 for(Int_t i = 0; i < headerList->GetEntries(); i++)
1211 {
1212 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(headerList->At(i));
1213 if(hijingGenHeader) break;
1214 }
1215
1216 if(!hijingGenHeader) return 0x0;
1217
1218 return hijingGenHeader;
1219}
1220
1221AliGenPythiaEventHeader* AlidNdPtAnalysisPbPbAOD::GetPythiaEventHeader(AliAODMCHeader *header)
1222{
1223 //
1224 // inspired by PWGJE/AliPWG4HighPtSpectra.cxx
1225 //
1226
1227 if(!header) return 0x0;
1228 AliGenPythiaEventHeader* PythiaGenHeader = NULL;
1229
1230 TList* headerList = header->GetCocktailHeaders();
1231
1232 for(Int_t i = 0; i < headerList->GetEntries(); i++)
1233 {
1234 PythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1235 if(PythiaGenHeader) break;
1236 }
1237
1238 if(!PythiaGenHeader) return 0x0;
1239
1240 return PythiaGenHeader;
1241}
1242
1243//________________________________________________________________________
1244Bool_t AlidNdPtAnalysisPbPbAOD::IsHijingParticle(const AliAODMCParticle *part, AliGenHijingEventHeader* hijingGenHeader){
1245
1246 // Check whether a particle is from Hijing or some injected
1247 // returns kFALSE if particle is injected
1248
1249 if(part->Label() > (hijingGenHeader->NProduced()-1)) return kFALSE;
1250 return kTRUE;
1251}
1252
1253//________________________________________________________________________
1254Bool_t AlidNdPtAnalysisPbPbAOD::IsPythiaParticle(const AliAODMCParticle *part, AliGenPythiaEventHeader* pythiaGenHeader){
1255
1256 // Check whether a particle is from Pythia or some injected
1257
1258 if(part->Label() > (pythiaGenHeader->NProduced()-1)) return kFALSE;
1259 return kTRUE;
1260}
1261
1262Double_t* AlidNdPtAnalysisPbPbAOD::GetArrayClone(Int_t n, Double_t* source)
1263{
1264 if (!source || n==0) return 0;
1265 Double_t* dest = new Double_t[n];
1266 for (Int_t i=0; i<n ; i++) { dest[i] = source[i]; }
1267 return dest;
1268}
1269
1270void AlidNdPtAnalysisPbPbAOD::Terminate(Option_t *)
1271{
1272
1273}
1274
1275