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