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