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