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