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