]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.cxx
Added distributions in phi, updated AddTask accordingly
[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
ea3cfeda 19// last modified: 17.10.2013
d25bcbe6 20//------------------------------------------------------------------------------
21
22
23#include "AlidNdPtAnalysisPbPbAOD.h"
24
25
26using namespace std;
27
28ClassImp(AlidNdPtAnalysisPbPbAOD)
29
ea3cfeda 30
d25bcbe6 31
32AlidNdPtAnalysisPbPbAOD::AlidNdPtAnalysisPbPbAOD(const char *name) : AliAnalysisTaskSE(name),
33fOutputList(0),
34// Histograms
ea3cfeda
PL
35fPt(0),
36fMCPt(0),
37fZvPtEtaCent(0),
bc2a9da9 38fPhiPtEtaCent(0),
ea3cfeda
PL
39fMCRecPrimZvPtEtaCent(0),
40fMCGenZvPtEtaCent(0),
41fMCRecSecZvPtEtaCent(0),
bc2a9da9
PL
42fMCRecPrimPhiPtEtaCent(0),
43fMCGenPhiPtEtaCent(0),
44fMCRecSecPhiPtEtaCent(0),
ea3cfeda
PL
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),
d25bcbe6 68//global
ea3cfeda 69fIsMonteCarlo(0),
d25bcbe6 70// event cut variables
ea3cfeda 71fCutMaxZVertex(10.),
d25bcbe6 72// track kinematic cut variables
ea3cfeda
PL
73fCutPtMin(0.15),
74fCutPtMax(200.),
75fCutEtaMin(-0.8),
76fCutEtaMax(0.8),
d25bcbe6 77// track quality cut variables
ea3cfeda
PL
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),
d25bcbe6 100// binning for THnSparse
101fMultNbins(0),
102fPtNbins(0),
103fPtCorrNbins(0),
72bb4ceb 104fPtCheckNbins(0),
d25bcbe6 105fEtaNbins(0),
72bb4ceb 106fEtaCheckNbins(0),
d25bcbe6 107fZvNbins(0),
108fCentralityNbins(0),
72bb4ceb 109fPhiNbins(0),
d25bcbe6 110fBinsMult(0),
111fBinsPt(0),
112fBinsPtCorr(0),
72bb4ceb 113fBinsPtCheck(0),
d25bcbe6 114fBinsEta(0),
72bb4ceb 115fBinsEtaCheck(0),
d25bcbe6 116fBinsZv(0),
72bb4ceb 117fBinsCentrality(0),
118fBinsPhi(0)
d25bcbe6 119{
72bb4ceb 120
121 for(Int_t i = 0; i < cqMax; i++)
122 {
ea3cfeda
PL
123 fCrossCheckAll[i] = 0;
124 fCrossCheckAcc[i] = 0;
72bb4ceb 125 }
126
d25bcbe6 127 fMultNbins = 0;
128 fPtNbins = 0;
129 fPtCorrNbins = 0;
72bb4ceb 130 fPtCheckNbins = 0;
d25bcbe6 131 fEtaNbins = 0;
72bb4ceb 132 fEtaCheckNbins = 0;
d25bcbe6 133 fZvNbins = 0;
134 fCentralityNbins = 0;
135 fBinsMult = 0;
136 fBinsPt = 0;
137 fBinsPtCorr = 0;
72bb4ceb 138 fBinsPtCheck = 0;
d25bcbe6 139 fBinsEta = 0;
72bb4ceb 140 fBinsEtaCheck = 0;
d25bcbe6 141 fBinsZv = 0;
142 fBinsCentrality = 0;
72bb4ceb 143 fBinsPhi = 0;
d25bcbe6 144
145 DefineOutput(1, TList::Class());
146}
147
148// destructor
149AlidNdPtAnalysisPbPbAOD::~AlidNdPtAnalysisPbPbAOD()
150{
9db7eb94 151
ea3cfeda
PL
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;
72bb4ceb 177
178 for(Int_t i = 0; i < cqMax; i++)
179 {
ea3cfeda
PL
180 if(fCrossCheckAll[i]) delete fCrossCheckAll[i]; fCrossCheckAll[i] = 0;
181 if(fCrossCheckAcc[i]) delete fCrossCheckAcc[i]; fCrossCheckAcc[i] = 0;
72bb4ceb 182 }
183
d25bcbe6 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};
72bb4ceb 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};
d25bcbe6 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
72bb4ceb 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
d25bcbe6 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); }
72bb4ceb 211 if (!fBinsPtCheck) { SetBinsPtCheck(20,binsPtCheckDefault); }
d25bcbe6 212 if (!fBinsEta) { SetBinsEta(31,binsEtaDefault); }
72bb4ceb 213 if (!fBinsEtaCheck) { SetBinsEtaCheck(7,binsEtaCheckDefault); }
d25bcbe6 214 if (!fBinsZv) { SetBinsZv(13,binsZvDefault); }
215 if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
72bb4ceb 216 if (!fBinsPhi) { SetBinsPhi(37,binsPhiDefault); }
d25bcbe6 217
218 Int_t binsZvPtEtaCent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
bc2a9da9 219 Int_t binsPhiPtEtaCent[4]={fPhiNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
d25bcbe6 220 Int_t binsZvMultCent[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
221
222 // define Histograms
ea3cfeda
PL
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
bc2a9da9
PL
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
ea3cfeda
PL
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
bc2a9da9
PL
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
ea3cfeda
PL
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");
d25bcbe6 373
95f26ffa 374
95f26ffa 375
72bb4ceb 376 Int_t binsDCAxyDCAzPtEtaPhi[6] = { 200,200, fPtCheckNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
9db7eb94 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};
d25bcbe6 379
ea3cfeda
PL
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");
9db7eb94 432
72bb4ceb 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;
ea3cfeda 445 const Int_t iNbinLength = 165;
72bb4ceb 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 {
bc2a9da9
PL
458 snprintf(cTempTitleAxis0All,255, "NcrossedRows before Cut");
459 snprintf(cTempTitleAxis0Acc,255, "NcrossedRows after Cut");
460 snprintf(cTempNameAxis0,255, "CrossedRows");
72bb4ceb 461 iNbin = iNbinRowsClusters;
462 dBinMin = 0;
463 dBinMax = 159.;
464 }
465 else if(iCheckQuant == cqNcluster)
466 {
bc2a9da9
PL
467 snprintf(cTempTitleAxis0All,255, "Nclusters before Cut");
468 snprintf(cTempTitleAxis0Acc,255, "Nclusters after Cut");
469 snprintf(cTempNameAxis0,255, "Clusters");
72bb4ceb 470 iNbin = iNbinRowsClusters;
471 dBinMin = 0;
472 dBinMax = 159.;
473 }
474 else if(iCheckQuant == cqChi)
475 {
bc2a9da9
PL
476 snprintf(cTempTitleAxis0All,255, "#Chi^{2}/cluster before Cut");
477 snprintf(cTempTitleAxis0Acc,255, "#Chi^{2}/cluster after Cut");
478 snprintf(cTempNameAxis0,255, "Chi");
72bb4ceb 479 iNbin = iNbinChi;
480 dBinMin = 0;
481 dBinMax = 10.;
482 }
ea3cfeda
PL
483 else if(iCheckQuant == cqLength)
484 {
bc2a9da9
PL
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");
ea3cfeda
PL
488 iNbin = iNbinLength;
489 dBinMin = 0;
490 dBinMax = 165.;
491 }
72bb4ceb 492
493 Int_t binsCheckPtEtaPhi[5] = { iNbin, fPtCheckNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
ea3cfeda 494// Int_t binsCheckPtEtaPhi[5] = { iNbin, fPtNbins-1, fEtaCheckNbins-1, 18, fCentralityNbins-1};
72bb4ceb 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
bc2a9da9
PL
498 snprintf(cFullTempName, 255, "f%sPtEtaPhiAll",cTempNameAxis0);
499 snprintf(cFullTempTitle, 255,"%s;%s;p_{T} (GeV/c);#eta;#phi;Centrality", cFullTempName, cTempTitleAxis0All);
ea3cfeda
PL
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();
72bb4ceb 504
bc2a9da9
PL
505 snprintf(cFullTempName, 255, "f%sPtEtaPhiAcc",cTempNameAxis0);
506 snprintf(cFullTempTitle, 255,"%s;%s;p_{T} (GeV/c);#eta;#phi;Centrality", cFullTempName, cTempTitleAxis0Acc);
ea3cfeda
PL
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();
72bb4ceb 511 } // end iCheckQuant
fad9b70b 512
ea3cfeda
PL
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
d25bcbe6 531 // Add Histos, Profiles etc to List
ea3cfeda 532 fOutputList->Add(fZvPtEtaCent);
bc2a9da9 533 fOutputList->Add(fPhiPtEtaCent);
ea3cfeda
PL
534 fOutputList->Add(fPt);
535 fOutputList->Add(fMCRecPrimZvPtEtaCent);
536 fOutputList->Add(fMCGenZvPtEtaCent);
537 fOutputList->Add(fMCRecSecZvPtEtaCent);
bc2a9da9
PL
538 fOutputList->Add(fMCRecPrimPhiPtEtaCent);
539 fOutputList->Add(fMCGenPhiPtEtaCent);
540 fOutputList->Add(fMCRecSecPhiPtEtaCent);
ea3cfeda
PL
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);
72bb4ceb 559 for(Int_t i = 0; i < cqMax; i++)
560 {
ea3cfeda
PL
561 fOutputList->Add(fCrossCheckAll[i]);
562 fOutputList->Add(fCrossCheckAcc[i]);
72bb4ceb 563 }
ea3cfeda
PL
564 fOutputList->Add(fCutPercClusters);
565 fOutputList->Add(fCutPercCrossed);
566 fOutputList->Add(fCrossCheckRowsLength);
567 fOutputList->Add(fCrossCheckClusterLength);
568 fOutputList->Add(fCrossCheckRowsLengthAcc);
569 fOutputList->Add(fCrossCheckClusterLengthAcc);
d25bcbe6 570
571 PostData(1, fOutputList);
572}
573
574void AlidNdPtAnalysisPbPbAOD::UserExec(Option_t *option)
575{
f856053f 576
d25bcbe6 577 // Main Loop
578 // called for each event
ea3cfeda 579 fEventStatistics->Fill("all events",1);
d25bcbe6 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;
b3341d37 587 //AliGenPythiaEventHeader *genPythiaHeader = NULL;
d25bcbe6 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;
9db7eb94 595 Bool_t bMotherIsHijingParticle = kFALSE;
b3341d37 596 //Bool_t bIsPythiaParticle = kFALSE;
d25bcbe6 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
bc2a9da9
PL
605 Double_t dMCTrackPhiPtEtaCent[4] = {0};
606 Double_t dTrackPhiPtEtaCent[4] = {0};
607
b3341d37 608 Double_t dDCA[2] = {0};
609
d25bcbe6 610 Double_t dMCEventZv = -100;
611 Double_t dEventZv = -100;
612 Int_t iAcceptedMultiplicity = 0;
613
ea3cfeda 614 fIsMonteCarlo = kFALSE;
d25bcbe6 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
ea3cfeda
PL
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); }
d25bcbe6 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
ea3cfeda 642 // fEventStatistics->Fill("events with only coll. cand.", 1);
d25bcbe6 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 {
ea3cfeda 653 fIsMonteCarlo = kTRUE;
d25bcbe6 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;
ea3cfeda 666 fEventStatistics->Fill("MC all events",1);
d25bcbe6 667 }
668
669 AliCentrality* aCentrality = eventAOD->GetCentrality();
670 Double_t dCentrality = aCentrality->GetCentralityPercentile("V0M");
671
672 if( dCentrality < 0 ) return;
ea3cfeda 673 fEventStatistics->Fill("after centrality selection",1);
d25bcbe6 674
675
676
677 // start with MC truth analysis
ea3cfeda 678 if(fIsMonteCarlo)
d25bcbe6 679 {
680
ea3cfeda 681 if( dMCEventZv > GetCutMaxZVertex() ) { return; }
d25bcbe6 682
683 dMCTrackZvPtEtaCent[0] = dMCEventZv;
684
ea3cfeda 685 fEventStatistics->Fill("MC afterZv cut",1);
d25bcbe6 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 {
ea3cfeda 698 fMCHijingPrim->Fill("IsPhysicalPrimary",1);
d25bcbe6 699 }
700 else
701 {
ea3cfeda 702 fMCHijingPrim->Fill("NOT a primary",1);
d25bcbe6 703 continue;
704 }
705
d0483ba3 706
d25bcbe6 707 //
708 // ======================== fill histograms ========================
709 dMCTrackZvPtEtaCent[1] = mcPart->Pt();
710 dMCTrackZvPtEtaCent[2] = mcPart->Eta();
711 dMCTrackZvPtEtaCent[3] = dCentrality;
bc2a9da9
PL
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);
d25bcbe6 719
720 bEventHasATrack = kTRUE;
721
d25bcbe6 722
ea3cfeda
PL
723 if( (dMCTrackZvPtEtaCent[1] > GetCutPtMin() ) &&
724 (dMCTrackZvPtEtaCent[1] < GetCutPtMax() ) &&
725 (dMCTrackZvPtEtaCent[2] > GetCutEtaMin() ) &&
726 (dMCTrackZvPtEtaCent[2] < GetCutEtaMax() ) )
d25bcbe6 727 {
ea3cfeda
PL
728 fMCPt->Fill(mcPart->Pt());
729 fMCCharge->Fill(mcPart->Charge()/3.);
d25bcbe6 730 bEventHasATrackInRange = kTRUE;
731 }
732
733 }
734 } // isMonteCarlo
ea3cfeda
PL
735
736 if(bEventHasATrack) { fEventStatistics->Fill("MC events with tracks",1); }
9db7eb94 737 if(bEventHasATrackInRange)
738 {
ea3cfeda
PL
739 fEventStatistics->Fill("MC events with tracks in range",1);
740 fMCEventStatisticsCentrality->Fill(dCentrality);
9db7eb94 741 }
d25bcbe6 742 bEventHasATrack = kFALSE;
743 bEventHasATrackInRange = kFALSE;
744
745
746
747 // Loop over recontructed tracks
748
749 dEventZv = eventAOD->GetPrimaryVertex()->GetZ();
ea3cfeda 750 if( TMath::Abs(dEventZv) > GetCutMaxZVertex() ) return;
d25bcbe6 751
ea3cfeda 752 fAllEventStatisticsCentrality->Fill(dCentrality/*, nTriggerFired*/);
d25bcbe6 753
ea3cfeda 754 fEventStatistics->Fill("after Zv cut",1);
d25bcbe6 755
756 dTrackZvPtEtaCent[0] = dEventZv;
757
ea3cfeda
PL
758 if(AreRelativeCutsEnabled())
759 {
760 if(!SetRelativeCuts(eventAOD)) return;
761 }
762
d25bcbe6 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
bc2a9da9
PL
773 dMCTrackPhiPtEtaCent[0] = 0;
774 dMCTrackPhiPtEtaCent[1] = 0;
775 dMCTrackPhiPtEtaCent[2] = 0;
776 dMCTrackPhiPtEtaCent[3] = 0;
777
d25bcbe6 778 bIsPrimary = kFALSE;
779
b3341d37 780 GetDCA(track, eventAOD, dDCA);
781
9db7eb94 782 Double_t dDCAxyDCAzPt[5] = { dDCA[0], dDCA[1], track->Pt(), track->Eta(), track->Phi() };
d0483ba3 783
ea3cfeda 784 fDCAPtAll->Fill(dDCAxyDCAzPt);
d0483ba3 785
ea3cfeda 786 if( !(IsTrackAccepted(track, dCentrality, eventAOD->GetMagneticField())) ) continue;
d25bcbe6 787
5747f2c6 788 dTrackZvPtEtaCent[1] = track->Pt();
789 dTrackZvPtEtaCent[2] = track->Eta();
790 dTrackZvPtEtaCent[3] = dCentrality;
791
bc2a9da9
PL
792 dTrackPhiPtEtaCent[0] = track->Phi();
793 dTrackPhiPtEtaCent[1] = track->Pt();
794 dTrackPhiPtEtaCent[2] = track->Eta();
795 dTrackPhiPtEtaCent[3] = dCentrality;
796
ea3cfeda 797 if( fIsMonteCarlo )
d25bcbe6 798 {
799 mcPart = (AliAODMCParticle*)stack->At(TMath::Abs(track->GetLabel()));
800 if( !mcPart ) { continue; }
801
802 // check for charge
ea3cfeda 803 // if( !(IsMCTrackAccepted(mcPart)) ) { continue; }
d25bcbe6 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
bc2a9da9
PL
816 dMCTrackPhiPtEtaCent[0] = mcPart->Phi();
817 dMCTrackPhiPtEtaCent[1] = mcPart->Pt();
818 dMCTrackPhiPtEtaCent[2] = mcPart->Eta();
819 dMCTrackPhiPtEtaCent[3] = dCentrality;
820
d25bcbe6 821 if(bIsPrimary && bIsHijingParticle)
822 {
ea3cfeda 823 fMCRecPrimZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
bc2a9da9 824 fMCRecPrimPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
ea3cfeda
PL
825 if( (TMath::Abs(mcPart->Eta()) < 0.8) && (dCentrality<5.) ) { fMCPdgPt->Fill(mcPart->Pt(), Form("%s",GetParticleName(mcPart->GetPdgCode())), 1); }
826 fMCDCAPtPrimary->Fill(dDCAxyDCAzPt);
d25bcbe6 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);
9db7eb94 835 bMotherIsHijingParticle = IsHijingParticle(moth, genHijingHeader);
d25bcbe6 836
837 if(bMotherIsHijingParticle) // only store secondaries, which come from a not embedded signal!
838 {
ea3cfeda
PL
839 fMCTrackStatusCode->Fill(Form("%d",mcPart->GetStatus()), 1);
840
d25bcbe6 841
ea3cfeda 842 fMCRecSecZvPtEtaCent->Fill(dMCTrackZvPtEtaCent);
bc2a9da9 843 fMCRecSecPhiPtEtaCent->Fill(dMCTrackPhiPtEtaCent);
ea3cfeda
PL
844 fMCDCAPtSecondary->Fill(dDCAxyDCAzPt);
845 fMCTrackPdgCode->Fill(Form("%s_H%i_H%i",GetParticleName(moth->GetPdgCode()),bMotherIsHijingParticle, bIsHijingParticle), 1);
d25bcbe6 846 // delete moth;
847 }
848 }
849 }
850 } // end isMonteCarlo
851
852 // ======================== fill histograms ========================
853
9db7eb94 854 // only keep prim and sec from not embedded signal
855 Bool_t bKeepMCTrack = kFALSE;
ea3cfeda 856 if(fIsMonteCarlo)
9db7eb94 857 {
858 if( (bIsHijingParticle && bIsPrimary) ^ (bMotherIsHijingParticle && !bIsPrimary) )
d0483ba3 859 {
9db7eb94 860 bKeepMCTrack = kTRUE;
861 }
862 else
863 {
864 continue;
d0483ba3 865 }
9db7eb94 866 }
867
868 bEventHasATrack = kTRUE;
869
ea3cfeda 870 fZvPtEtaCent->Fill(dTrackZvPtEtaCent);
bc2a9da9
PL
871 fPhiPtEtaCent->Fill(dTrackPhiPtEtaCent);
872
ea3cfeda 873 fDCAPtAccepted->Fill(dDCAxyDCAzPt);
9db7eb94 874
ea3cfeda
PL
875 if( (dTrackZvPtEtaCent[1] > GetCutPtMin()) &&
876 (dTrackZvPtEtaCent[1] < GetCutPtMax()) &&
877 (dTrackZvPtEtaCent[2] > GetCutEtaMin()) &&
878 (dTrackZvPtEtaCent[2] < GetCutEtaMax()) )
9db7eb94 879 {
880 iAcceptedMultiplicity++;
881 bEventHasATrackInRange = kTRUE;
ea3cfeda
PL
882 fPt->Fill(track->Pt());
883 fCharge->Fill(track->Charge());
9db7eb94 884 }
d25bcbe6 885 } // end track loop
886
ea3cfeda 887 if(bEventHasATrack) { fEventStatistics->Fill("events with tracks",1); bEventHasATrack = kFALSE; }
d25bcbe6 888
889 if(bEventHasATrackInRange)
890 {
ea3cfeda
PL
891 fEventStatistics->Fill("events with tracks in range",1);
892 fEventStatisticsCentrality->Fill(dCentrality);
d25bcbe6 893 bEventHasATrackInRange = kFALSE;
fad9b70b 894
ea3cfeda
PL
895 if(bIsEventSelectedMB) fEventStatisticsCentralityTrigger->Fill(dCentrality, 0);
896 if(bIsEventSelectedSemi) fEventStatisticsCentralityTrigger->Fill(dCentrality, 1);
897 if(bIsEventSelectedCentral) fEventStatisticsCentralityTrigger->Fill(dCentrality, 2);
d25bcbe6 898 }
899
900 Double_t dEventZvMultCent[3] = {dEventZv, iAcceptedMultiplicity, dCentrality};
ea3cfeda 901 fZvMultCent->Fill(dEventZvMultCent);
d25bcbe6 902
d25bcbe6 903 PostData(1, fOutputList);
904
905}
906
ea3cfeda
PL
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)
d25bcbe6 973{
974 if(!tr) return kFALSE;
975
976 if(tr->Charge()==0) { return kFALSE; }
977
ea3cfeda
PL
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);
d25bcbe6 1010 Double_t dNClustersTPC = tr->GetTPCNcls();
1011 Double_t dCrossedRowsTPC = tr->GetTPCClusterInfo(2,1);
ea3cfeda 1012 Double_t dFindableClustersTPC = tr->GetTPCNclsF();
9db7eb94 1013 Double_t dChi2PerClusterTPC = (dNClustersTPC>0)?tr->Chi2perNDF()*(dNClustersTPC-5)/dNClustersTPC:-1.; // see AliDielectronVarManager.h
1014
1015 // hAllCrossedRowsTPC->Fill(dCrossedRowsTPC);
d25bcbe6 1016
ea3cfeda
PL
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;
72bb4ceb 1026
1027
1028 FillDebugHisto(dCheck, dKine, dCentrality, kFALSE);
1029
ea3cfeda
PL
1030 // first cut on length
1031
1032 if( DoCutLengthInTPCPtDependent() && ( dLengthInTPC < GetPrefactorLengthInTPCPtDependent()*(130-5*TMath::Abs(1./tr->Pt())) ) ) { return kFALSE; }
95f26ffa 1033
9db7eb94 1034 // filter bit 5
1035 if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobal)) ) { return kFALSE; }
95f26ffa 1036
9db7eb94 1037 // filter bit 4
1038 // if(!(tr->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) ) { return kFALSE; }
1039
1040 // hFilterCrossedRowsTPC->Fill(dCrossedRowsTPC);
95f26ffa 1041
d0483ba3 1042
ea3cfeda
PL
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; }
72bb4ceb 1047
ea3cfeda 1048 if (!(tr->GetStatus() & AliVTrack::kITSrefit)) { return kFALSE; } // no ITS refit
72bb4ceb 1049
ea3cfeda
PL
1050 // do a relativ cut in Nclusters, both time at 80% of mean
1051 // if(fIsMonteCarlo)
1052 // {
1053 // if(dNClustersTPC < 88) { return kFALSE; }
d25bcbe6 1054 // }
ea3cfeda
PL
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;
d25bcbe6 1113}
1114
72bb4ceb 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};
ea3cfeda 1122 fCrossCheckAcc[iCrossCheck]->Fill(dFillIt);
72bb4ceb 1123 }
ea3cfeda
PL
1124
1125 fCrossCheckRowsLengthAcc->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqCrossedRows]);
1126 fCrossCheckClusterLengthAcc->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqNcluster]);
72bb4ceb 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};
ea3cfeda 1133 fCrossCheckAll[iCrossCheck]->Fill(dFillIt);
72bb4ceb 1134 }
ea3cfeda
PL
1135
1136 fCrossCheckRowsLength->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqCrossedRows]);
1137 fCrossCheckClusterLength->Fill(dCrossCheckVar[cqLength], dCrossCheckVar[cqNcluster]);
72bb4ceb 1138 }
1139
1140 return kTRUE;
1141
1142}
1143
b3341d37 1144Bool_t AlidNdPtAnalysisPbPbAOD::GetDCA(const AliAODTrack *track, AliAODEvent *evt, Double_t d0z0[2])
d0483ba3 1145{
b3341d37 1146 // function adapted from AliDielectronVarManager.h
d0483ba3 1147
b3341d37 1148 if(track->TestBit(AliAODTrack::kIsDCA)){
1149 d0z0[0]=track->DCA();
1150 d0z0[1]=track->ZAtDCA();
1151 return kTRUE;
d0483ba3 1152 }
1153
1154 Bool_t ok=kFALSE;
b3341d37 1155 if(evt) {
d0483ba3 1156 Double_t covd0z0[3];
b3341d37 1157 //AliAODTrack copy(*track);
1158 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
d0483ba3 1159
1160 Float_t xstart = etp.GetX();
b3341d37 1161 if(xstart>3.) {
d0483ba3 1162 d0z0[0]=-999.;
1163 d0z0[1]=-999.;
1164 //printf("This method can be used only for propagation inside the beam pipe \n");
b3341d37 1165 return kFALSE;
d0483ba3 1166 }
b3341d37 1167
1168
d0483ba3 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);
b3341d37 1172 //ok = copy.PropagateToDCA(vtx,fBzkG,kVeryBig,d0z0,covd0z0);
d0483ba3 1173 }
1174 if(!ok){
1175 d0z0[0]=-999.;
1176 d0z0[1]=-999.;
d0483ba3 1177 }
b3341d37 1178 return ok;
d0483ba3 1179}
1180
b3341d37 1181
d25bcbe6 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
ea3cfeda 1247 // returns kFALSE if particle is injected
d25bcbe6 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;
b3341d37 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;
d25bcbe6 1268}
1269
b3341d37 1270void AlidNdPtAnalysisPbPbAOD::Terminate(Option_t *)
1271{
1272
1273}
1274
1275