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