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