changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / AlidNdPtAnalysisPbPb2011.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 // AlidNdPtAnalysisPbPb2011 class. 
17 // 
18 // a. functionality:
19 // - fills analysis control histograms
20 // - fills generic correction matrices 
21 // - generates correction matrices 
22 //
23 // b. data members:
24 // - generic correction matrices
25 // - control histograms
26 //
27 // Author: J.Otwinowski 04/11/2008 
28 // last change: 2011-04-04 by M.Knichel
29 //------------------------------------------------------------------------------
30
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TCanvas.h"
34 #include "THnSparse.h"
35
36 #include "AliHeader.h"  
37 #include "AliGenEventHeader.h"  
38 #include "AliInputEventHandler.h"  
39 #include "AliAnalysisManager.h"  
40 #include "AliStack.h"  
41 #include "AliESDEvent.h"  
42 #include "AliMCEvent.h"  
43 #include "AliESDtrackCuts.h"  
44 #include "AliLog.h" 
45 #include "AliMultiplicity.h"
46 #include "AliTracker.h"
47
48 #include "AliCentrality.h"
49
50 #include "AlidNdPtEventCuts.h"
51 #include "AlidNdPtAcceptanceCuts.h"
52 #include "AliPhysicsSelection.h"
53 #include "AliTriggerAnalysis.h"
54
55 #include "AliPWG0Helper.h"
56 #include "AlidNdPtHelper.h"
57 #include "AlidNdPtAnalysisPbPb2011.h"
58
59 #include "TSystem.h"
60 #include "TROOT.h"
61
62 using namespace std;
63
64 ClassImp(AlidNdPtAnalysisPbPb2011)
65
66 //_____________________________________________________________________________
67 AlidNdPtAnalysisPbPb2011::AlidNdPtAnalysisPbPb2011(): AlidNdPt(),
68 fAnalysisFolder(0),
69 fHistogramsOn(kFALSE),
70
71 // rec. track pt vs true track pt correlation matrix 
72 fTrackPtCorrelationMatrix(0),
73
74 // event level correction
75 fGenEventMatrix(0),
76 fTriggerEventMatrix(0),
77 fRecEventMatrix(0),
78
79 //
80 // track-event level correction 
81 //
82 fGenTrackEventMatrix(0),
83
84 fTriggerTrackEventMatrix(0),
85
86 fRecTrackEventMatrix(0),
87
88 // track rec. efficiency correction (fRecPrimTrackMatrix / fGenPrimTrackMatrix)
89 fGenTrackMatrix(0),
90 fGenPrimTrackMatrix(0),
91 fRecPrimTrackMatrix(0),
92
93 // secondary track contamination correction (fRecSecTrackMatrix / fRecTrackMatrix)
94 fRecTrackMatrix(0),
95 fRecSecTrackMatrix(0),
96
97 // multiple rec. track contamination corrections (fRecMultTrackMatrix / fRecTrackMatrix)
98 fRecMultTrackMatrix(0),
99
100 // event control histograms
101 fMCEventHist1(0),
102 fRecEventHist1(0),
103 fRecEventHist2(0),
104 fRecMCEventHist1(0),
105 fRecMCEventHist2(0),
106
107 // rec. pt and eta resolution w.r.t MC
108 fRecMCTrackHist1(0),
109
110 //multple reconstructed tracks
111 fMCMultRecTrackHist1(0), 
112
113 // rec. track control histograms
114 fRecTrackHist3(0),
115
116 fTriggerAnalysis(0),
117
118 fCentralityEstimator(0),
119
120 fMultNbins(0),
121 fPtNbins(0),
122 fPtCorrNbins(0),
123 fEtaNbins(0),
124 fZvNbins(0),
125 fCentralityNbins(0),
126 fBinsMult(0),
127 fBinsPt(0),
128 fBinsPtCorr(0),
129 fBinsEta(0),
130 fBinsZv(0),
131 fBinsCentrality(0),
132
133 fIsInit(kFALSE)
134 {
135   // default constructor
136   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { 
137     fMCTrackHist1[i]=0;     
138     fMCPrimTrackHist1[i]=0;     
139     fMCPrimTrackHist2[i]=0;     
140     fMCSecTrackHist1[i]=0;     
141     fRecTrackHist1[i]=0;     
142     fRecTrackHist2[i]=0;     
143     fRecTrackMultHist1[i]=0;     
144   }
145   
146   fMultNbins = 0;
147   fPtNbins = 0;
148   fPtCorrNbins = 0;
149   fEtaNbins = 0;
150   fZvNbins = 0;
151   fCentralityNbins = 0;
152   fBinsMult = 0;
153   fBinsPt = 0;
154   fBinsPtCorr = 0;
155   fBinsEta = 0;
156   fBinsEta = 0;
157   fBinsZv = 0;
158   fBinsCentrality = 0;
159   
160   
161   //Init();
162   SetCentralityEstimator();
163 }
164
165 //_____________________________________________________________________________
166 AlidNdPtAnalysisPbPb2011::AlidNdPtAnalysisPbPb2011(Char_t* name, Char_t* title): AlidNdPt(name,title),
167 fAnalysisFolder(0),
168 fHistogramsOn(kFALSE),
169
170 // rec. track pt vs true track pt correlation matrix 
171 fTrackPtCorrelationMatrix(0),
172
173 // event level correction
174 fGenEventMatrix(0),
175 fTriggerEventMatrix(0),
176 fRecEventMatrix(0),
177
178 //
179 // track-event level correction 
180 //
181 fGenTrackEventMatrix(0),
182
183 fTriggerTrackEventMatrix(0),
184
185 fRecTrackEventMatrix(0),
186
187 // track rec. efficiency correction (fRecPrimTrackMatrix / fGenPrimTrackMatrix)
188 fGenTrackMatrix(0),
189 fGenPrimTrackMatrix(0),
190 fRecPrimTrackMatrix(0),
191
192 // secondary track contamination correction (fRecTrackMatrix - fRecSecTrackMatrix)
193 fRecTrackMatrix(0),
194 fRecSecTrackMatrix(0),
195
196 // multiple rec. track contamination corrections (fRecTrackMatrix - fRecMultTrackMatrix)
197 fRecMultTrackMatrix(0),
198
199
200
201
202 // event control histograms
203 fMCEventHist1(0),
204 fRecEventHist1(0),
205 fRecEventHist2(0),
206 fRecMCEventHist1(0),
207 fRecMCEventHist2(0),
208
209 // rec. pt and eta resolution w.r.t MC
210 fRecMCTrackHist1(0),
211
212 //multiple reconstructed tracks
213 fMCMultRecTrackHist1(0), 
214
215 // rec. track control histograms
216 fRecTrackHist3(0),
217
218 fTriggerAnalysis(0),
219
220 fCentralityEstimator(0),
221
222 fMultNbins(0),
223 fPtNbins(0),
224 fPtCorrNbins(0),
225 fEtaNbins(0),
226 fZvNbins(0),
227 fCentralityNbins(0),
228 fBinsMult(0),
229 fBinsPt(0),
230 fBinsPtCorr(0),
231 fBinsEta(0),
232 fBinsZv(0),
233 fBinsCentrality(0),
234
235 fIsInit(kFALSE)  
236
237 {
238   //
239   // constructor
240   //
241   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { 
242     fMCTrackHist1[i]=0;     
243     fMCPrimTrackHist1[i]=0;     
244     fMCPrimTrackHist2[i]=0;     
245     fMCSecTrackHist1[i]=0;     
246     fRecTrackHist1[i]=0;     
247     fRecTrackHist2[i]=0;     
248     fRecTrackMultHist1[i]=0; 
249   }
250   
251   fMultNbins = 0;
252   fPtNbins = 0;
253   fPtCorrNbins = 0;
254   fEtaNbins = 0;
255   fZvNbins = 0;
256   fCentralityNbins = 0;
257   fBinsMult = 0;
258   fBinsPt = 0;
259   fBinsPtCorr = 0;
260   fBinsEta = 0;
261   fBinsEta = 0;
262   fBinsZv = 0;
263   fBinsCentrality = 0;
264   
265   // Init();
266   SetCentralityEstimator();
267 }
268
269 //_____________________________________________________________________________
270 AlidNdPtAnalysisPbPb2011::~AlidNdPtAnalysisPbPb2011() {
271   //
272   // destructor
273   //
274   if(fTrackPtCorrelationMatrix) delete fTrackPtCorrelationMatrix; fTrackPtCorrelationMatrix=0;
275   //
276   if(fGenEventMatrix) delete fGenEventMatrix; fGenEventMatrix=0;
277   if(fTriggerEventMatrix) delete fTriggerEventMatrix; fTriggerEventMatrix=0;
278   if(fRecEventMatrix) delete fRecEventMatrix; fRecEventMatrix=0;
279   
280   //
281   if(fGenTrackEventMatrix) delete fGenTrackEventMatrix; fGenTrackEventMatrix=0;
282   if(fTriggerEventMatrix) delete fTriggerEventMatrix; fTriggerEventMatrix=0;
283   if(fRecTrackEventMatrix) delete fRecTrackEventMatrix; fRecTrackEventMatrix=0;
284   
285   //
286   if(fGenTrackMatrix) delete fGenTrackMatrix; fGenTrackMatrix=0;
287   if(fGenPrimTrackMatrix) delete fGenPrimTrackMatrix; fGenPrimTrackMatrix=0;
288   if(fRecPrimTrackMatrix) delete fRecPrimTrackMatrix; fRecPrimTrackMatrix=0;
289   //
290   if(fRecTrackMatrix) delete fRecTrackMatrix; fRecTrackMatrix=0;
291   if(fRecSecTrackMatrix) delete fRecSecTrackMatrix; fRecSecTrackMatrix=0;
292   // 
293   if(fRecMultTrackMatrix) delete fRecMultTrackMatrix; fRecMultTrackMatrix=0;
294   //
295   // Control histograms
296   //
297   if(fMCEventHist1) delete fMCEventHist1; fMCEventHist1=0;
298   if(fRecEventHist1) delete fRecEventHist1; fRecEventHist1=0;
299   if(fRecEventHist2) delete fRecEventHist2; fRecEventHist2=0;
300   if(fRecMCEventHist1) delete fRecMCEventHist1; fRecMCEventHist1=0;
301   if(fRecMCEventHist2) delete fRecMCEventHist2; fRecMCEventHist2=0;
302   
303   //
304   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { 
305     if(fMCTrackHist1[i]) delete fMCTrackHist1[i]; fMCTrackHist1[i]=0;
306     if(fMCPrimTrackHist1[i]) delete fMCPrimTrackHist1[i]; fMCPrimTrackHist1[i]=0;
307     if(fMCPrimTrackHist2[i]) delete fMCPrimTrackHist2[i]; fMCPrimTrackHist2[i]=0;
308     if(fMCSecTrackHist1[i]) delete fMCSecTrackHist1[i]; fMCSecTrackHist1[i]=0;
309     if(fRecTrackHist1[i]) delete fRecTrackHist1[i]; fRecTrackHist1[i]=0;
310     if(fRecTrackHist2[i]) delete fRecTrackHist2[i]; fRecTrackHist2[i]=0;
311     if(fRecTrackMultHist1[i]) delete fRecTrackMultHist1[i]; fRecTrackMultHist1[i]=0;
312   }
313   if(fRecMCTrackHist1) delete fRecMCTrackHist1; fRecMCTrackHist1=0;
314   if(fMCMultRecTrackHist1) delete fMCMultRecTrackHist1; fMCMultRecTrackHist1=0; 
315   if(fRecTrackHist3) delete fRecTrackHist3; fRecTrackHist3=0; 
316   //
317   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
318   
319   if (fTriggerAnalysis) delete fTriggerAnalysis;  fTriggerAnalysis = 0;
320   
321   if (fBinsMult) delete[] fBinsMult; fBinsMult=0;
322   if (fBinsPt) delete[] fBinsPt; fBinsPt=0;
323   if (fBinsPtCorr) delete[] fBinsPtCorr; fBinsPtCorr=0;
324   if (fBinsEta) delete[] fBinsEta; fBinsEta=0;
325   if (fBinsZv) delete[] fBinsZv; fBinsZv=0;
326   if (fBinsCentrality) delete[] fBinsCentrality; fBinsCentrality=0;  
327 }
328
329 //_____________________________________________________________________________
330 void AlidNdPtAnalysisPbPb2011::Init() {
331   
332   //define default binning
333   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 }; 
334   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};
335   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};    
336   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};
337   Double_t binsZvDefault[13] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};
338   Double_t binsCentralityDefault[12] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.};  
339   
340   // if no binning is set, use the default
341   if (!fBinsMult)                       { SetBinsMult(48,binsMultDefault); }
342   if (!fBinsPt)                         { SetBinsPt(82,binsPtDefault); }
343   if (!fBinsPtCorr)                     { SetBinsPtCorr(37,binsPtCorrDefault); }
344   if (!fBinsEta)                                { SetBinsEta(31,binsEtaDefault); }
345   if (!fBinsZv)                         { SetBinsZv(13,binsZvDefault); }  
346   if (!fBinsCentrality) { SetBinsCentrality(12,binsCentralityDefault); }
347   
348   for(Int_t i = 0; i < fPtNbins; i++)                           { AliDebug(AliLog::kDebug, Form("ptBin %i: %.2f", i, fBinsPt[i])); }
349   for(Int_t i = 0; i < fPtCorrNbins; i++)               { AliDebug(AliLog::kDebug, Form("ptCorrBin %i: %.2f", i, fBinsPtCorr[i])); }
350   for(Int_t i = 0; i < fEtaNbins; i++)                  { AliDebug(AliLog::kDebug, Form("etaBin %i: %.2f", i, fBinsEta[i])); }
351   for(Int_t i = 0; i < fZvNbins; i++)                           { AliDebug(AliLog::kDebug, Form("zvBin %i: %.2f", i, fBinsZv[i])); }
352   for(Int_t i = 0; i < fMultNbins; i++)                 { AliDebug(AliLog::kDebug, Form("multBin %i: %.2f", i, fBinsMult[i])); }
353   for(Int_t i = 0; i < fCentralityNbins; i++)   { AliDebug(AliLog::kDebug, Form("centBin %i: %.2f", i, fBinsCentrality[i])); }
354   
355   Int_t binsTrackEventCorrMatrix[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
356   Int_t binsTrackEvent[4]={fZvNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
357   Int_t binsTrackPtCorrelationMatrix[4]={fPtNbins-1,fPtNbins-1,fEtaNbins-1,fCentralityNbins-1};
358   
359   fTrackPtCorrelationMatrix = new THnSparseF("fTrackPtCorrelationMatrix","Pt:mcPt:mcEta:Centrality",4,binsTrackPtCorrelationMatrix);
360   fTrackPtCorrelationMatrix->SetBinEdges(0,fBinsPt);
361   fTrackPtCorrelationMatrix->SetBinEdges(1,fBinsPt);
362   fTrackPtCorrelationMatrix->SetBinEdges(2,fBinsEta);
363   fTrackPtCorrelationMatrix->SetBinEdges(3,fBinsCentrality);
364   fTrackPtCorrelationMatrix->GetAxis(0)->SetTitle("Pt (GeV/c)");
365   fTrackPtCorrelationMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
366   fTrackPtCorrelationMatrix->GetAxis(2)->SetTitle("mcEta");
367   fTrackPtCorrelationMatrix->GetAxis(3)->SetTitle("Centrality");
368   fTrackPtCorrelationMatrix->Sumw2();
369   
370   //
371   // Efficiency and contamination correction matrices
372   //
373   Int_t binsEventMatrix[3]={fZvNbins-1,fMultNbins-1,fCentralityNbins-1};
374   //   Double_t minEventMatrix[3]={-30.,-0.5,0.}; 
375   //   Double_t maxEventMatrix[3]={30.,10000.5,100.}; 
376   
377   fGenEventMatrix = new THnSparseF("fGenEventMatrix","mcZv:multMB:Centrality",3,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
378   fGenEventMatrix->SetBinEdges(0,fBinsZv);
379   fGenEventMatrix->SetBinEdges(1,fBinsMult);
380   fGenEventMatrix->SetBinEdges(2,fBinsCentrality);
381   fGenEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
382   fGenEventMatrix->GetAxis(1)->SetTitle("multiplicity MB");
383   fGenEventMatrix->GetAxis(2)->SetTitle("Centrality");
384   fGenEventMatrix->Sumw2();
385   //
386   fTriggerEventMatrix = new THnSparseF("fTriggerEventMatrix","mcZv:multMB:Centrality",3,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
387   fTriggerEventMatrix->SetBinEdges(0,fBinsZv);
388   fTriggerEventMatrix->SetBinEdges(1,fBinsMult);
389   fTriggerEventMatrix->SetBinEdges(2,fBinsCentrality);
390   fTriggerEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
391   fTriggerEventMatrix->GetAxis(1)->SetTitle("multiplicity MB");
392   fTriggerEventMatrix->GetAxis(2)->SetTitle("Centrality");
393   fTriggerEventMatrix->Sumw2();
394   //
395   fRecEventMatrix = new THnSparseF("fRecEventMatrix","mcZv:multMB:Centrality",3,binsEventMatrix); //,minEventMatrix,maxEventMatrix);
396   fRecEventMatrix->SetBinEdges(0,fBinsZv);
397   fRecEventMatrix->SetBinEdges(1,fBinsMult);
398   fRecEventMatrix->SetBinEdges(2,fBinsCentrality);    
399   fRecEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
400   fRecEventMatrix->GetAxis(1)->SetTitle("multiplicity MB");
401   fRecEventMatrix->GetAxis(2)->SetTitle("Centrality");
402   fRecEventMatrix->Sumw2();
403   
404   // 
405   // track to event corrections
406   //
407   fGenTrackEventMatrix = new THnSparseF("fGenTrackEventMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
408   fGenTrackEventMatrix->SetBinEdges(0,fBinsZv);
409   fGenTrackEventMatrix->SetBinEdges(1,fBinsPtCorr);
410   fGenTrackEventMatrix->SetBinEdges(2,fBinsEta);
411   fGenTrackEventMatrix->SetBinEdges(3,fBinsCentrality);
412   fGenTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
413   fGenTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
414   fGenTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
415   fGenTrackEventMatrix->GetAxis(3)->SetTitle("Centrality");
416   fGenTrackEventMatrix->Sumw2();
417   //
418   fTriggerTrackEventMatrix = new THnSparseF("fTriggerTrackEventMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
419   fTriggerTrackEventMatrix->SetBinEdges(0,fBinsZv);
420   fTriggerTrackEventMatrix->SetBinEdges(1,fBinsPtCorr);
421   fTriggerTrackEventMatrix->SetBinEdges(2,fBinsEta);
422   fTriggerTrackEventMatrix->SetBinEdges(3,fBinsCentrality);
423   fTriggerTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
424   fTriggerTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
425   fTriggerTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
426   fTriggerTrackEventMatrix->GetAxis(3)->SetTitle("Centrality");
427   fTriggerTrackEventMatrix->Sumw2();
428   //
429   fRecTrackEventMatrix = new THnSparseF("fRecTrackEventMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
430   fRecTrackEventMatrix->SetBinEdges(0,fBinsZv);
431   fRecTrackEventMatrix->SetBinEdges(1,fBinsPtCorr);
432   fRecTrackEventMatrix->SetBinEdges(2,fBinsEta);
433   fRecTrackEventMatrix->SetBinEdges(3,fBinsCentrality);
434   fRecTrackEventMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
435   fRecTrackEventMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
436   fRecTrackEventMatrix->GetAxis(2)->SetTitle("mcEta");
437   fRecTrackEventMatrix->GetAxis(3)->SetTitle("Centrality");
438   fRecTrackEventMatrix->Sumw2();
439   
440   //
441   // tracks correction matrices
442   //
443   fGenTrackMatrix = new THnSparseF("fGenTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
444   fGenTrackMatrix->SetBinEdges(0,fBinsZv);
445   fGenTrackMatrix->SetBinEdges(1,fBinsPtCorr);
446   fGenTrackMatrix->SetBinEdges(2,fBinsEta);
447   fGenTrackMatrix->SetBinEdges(3,fBinsCentrality);
448   fGenTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
449   fGenTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
450   fGenTrackMatrix->GetAxis(2)->SetTitle("mcEta");
451   fGenTrackMatrix->GetAxis(3)->SetTitle("Centrality");
452   fGenTrackMatrix->Sumw2();
453   
454   fGenPrimTrackMatrix = new THnSparseF("fGenPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
455   fGenPrimTrackMatrix->SetBinEdges(0,fBinsZv);
456   fGenPrimTrackMatrix->SetBinEdges(1,fBinsPtCorr);
457   fGenPrimTrackMatrix->SetBinEdges(2,fBinsEta);
458   fGenPrimTrackMatrix->SetBinEdges(3,fBinsCentrality);
459   fGenPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
460   fGenPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
461   fGenPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
462   fGenPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
463   fGenPrimTrackMatrix->Sumw2();
464   
465   
466   fRecPrimTrackMatrix = new THnSparseF("fRecPrimTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
467   fRecPrimTrackMatrix->SetBinEdges(0,fBinsZv);
468   fRecPrimTrackMatrix->SetBinEdges(1,fBinsPtCorr);
469   fRecPrimTrackMatrix->SetBinEdges(2,fBinsEta);
470   fRecPrimTrackMatrix->SetBinEdges(3,fBinsCentrality);
471   fRecPrimTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
472   fRecPrimTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
473   fRecPrimTrackMatrix->GetAxis(2)->SetTitle("mcEta");
474   fRecPrimTrackMatrix->GetAxis(3)->SetTitle("Centrality");
475   fRecPrimTrackMatrix->Sumw2();
476   
477   //
478   fRecTrackMatrix = new THnSparseF("fRecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
479   fRecTrackMatrix->SetBinEdges(0,fBinsZv);
480   fRecTrackMatrix->SetBinEdges(1,fBinsPtCorr);
481   fRecTrackMatrix->SetBinEdges(2,fBinsEta);
482   fRecTrackMatrix->SetBinEdges(3,fBinsCentrality);
483   fRecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
484   fRecTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
485   fRecTrackMatrix->GetAxis(2)->SetTitle("mcEta");
486   fRecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
487   fRecTrackMatrix->Sumw2();
488   
489   fRecSecTrackMatrix = new THnSparseF("fRecSecTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
490   fRecSecTrackMatrix->SetBinEdges(0,fBinsZv);
491   fRecSecTrackMatrix->SetBinEdges(1,fBinsPtCorr);
492   fRecSecTrackMatrix->SetBinEdges(2,fBinsEta);
493   fRecSecTrackMatrix->SetBinEdges(3,fBinsCentrality);
494   fRecSecTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
495   fRecSecTrackMatrix->GetAxis(1)->SetTitle("Pt (GeV/c)");
496   fRecSecTrackMatrix->GetAxis(2)->SetTitle("Eta");
497   fRecSecTrackMatrix->GetAxis(3)->SetTitle("Centrality");
498   fRecSecTrackMatrix->Sumw2();
499   
500   //
501   fRecMultTrackMatrix = new THnSparseF("fRecMultTrackMatrix","mcZv:mcPt:mcEta:Centrality",4,binsTrackEventCorrMatrix);
502   fRecMultTrackMatrix->SetBinEdges(0,fBinsZv);
503   fRecMultTrackMatrix->SetBinEdges(1,fBinsPtCorr);
504   fRecMultTrackMatrix->SetBinEdges(2,fBinsEta);
505   fRecMultTrackMatrix->SetBinEdges(3,fBinsCentrality);
506   fRecMultTrackMatrix->GetAxis(0)->SetTitle("mcZv (cm)");
507   fRecMultTrackMatrix->GetAxis(1)->SetTitle("mcPt (GeV/c)");
508   fRecMultTrackMatrix->GetAxis(2)->SetTitle("mcEta");
509   fRecMultTrackMatrix->GetAxis(3)->SetTitle("Centrality");
510   fRecMultTrackMatrix->Sumw2();
511   
512   //
513   // Control analysis histograms
514   //
515   Int_t binsMCEventHist1[4]={100,100,fZvNbins-1,fCentralityNbins-1};
516   Double_t minMCEventHist1[4]={-0.1,-0.1,-30.,0.}; 
517   Double_t maxMCEventHist1[4]={0.1,0.1,30.,100.}; 
518   fMCEventHist1 = new THnSparseF("fMCEventHist1","mcXv:mcYv:mcZv:Centrality",4,binsMCEventHist1,minMCEventHist1,maxMCEventHist1);
519   fMCEventHist1->SetBinEdges(2,fBinsZv);
520   fMCEventHist1->SetBinEdges(3,fBinsCentrality);
521   fMCEventHist1->GetAxis(0)->SetTitle("mcXv (cm)");
522   fMCEventHist1->GetAxis(1)->SetTitle("mcYv (cm)");
523   fMCEventHist1->GetAxis(2)->SetTitle("mcZv (cm)");
524   fMCEventHist1->GetAxis(3)->SetTitle("Centrality");  
525   fMCEventHist1->Sumw2();
526   
527   //
528   Int_t binsRecEventHist1[4]={100,100,fZvNbins-1,fCentralityNbins-1};
529   Double_t minRecEventHist1[4]={-3.,-3.,-30.,0.}; 
530   Double_t maxRecEventHist1[4]={3.,3.,30.,100.}; 
531   
532   fRecEventHist1 = new THnSparseF("fRecEventHist1","Xv:Yv:Zv:Centrality",4,binsRecEventHist1,minRecEventHist1,maxRecEventHist1);
533   fRecEventHist1->SetBinEdges(2,fBinsZv);
534   fRecEventHist1->SetBinEdges(3,fBinsCentrality);
535   fRecEventHist1->GetAxis(0)->SetTitle("Xv (cm)");
536   fRecEventHist1->GetAxis(1)->SetTitle("Yv (cm)");
537   fRecEventHist1->GetAxis(2)->SetTitle("Zv (cm)");
538   fRecEventHist1->GetAxis(3)->SetTitle("Centrality");  
539   fRecEventHist1->Sumw2();
540   
541   //
542   Int_t binsRecEventHist2[3]={fZvNbins-1, fMultNbins-1, fCentralityNbins-1};
543   Double_t minRecEventHist2[3]={-30., -0.5, 0.}; 
544   Double_t maxRecEventHist2[3]={30., 10000.5, 100.}; 
545   
546   fRecEventHist2 = new THnSparseF("fRecEventHist2","Zv:multMB:Centrality",3,binsRecEventHist2,minRecEventHist2,maxRecEventHist2);
547   fRecEventHist2->SetBinEdges(0,fBinsZv);
548   fRecEventHist2->SetBinEdges(1,fBinsMult);
549   fRecEventHist2->SetBinEdges(2,fBinsCentrality);
550   fRecEventHist2->GetAxis(0)->SetTitle("Zv (cm)");
551   fRecEventHist2->GetAxis(1)->SetTitle("multiplicity MB");
552   fRecEventHist2->GetAxis(2)->SetTitle("Centrality"); 
553   fRecEventHist2->Sumw2();
554   
555   //
556   Double_t kFact = 0.1;
557   Int_t binsRecMCEventHist1[4]={100,100,100, fCentralityNbins-1};
558   Double_t minRecMCEventHist1[4]={-10.0*kFact,-10.0*kFact,-10.0*kFact, 0.}; 
559   Double_t maxRecMCEventHist1[4]={10.0*kFact,10.0*kFact,10.0*kFact, 100.}; 
560   
561   fRecMCEventHist1 = new THnSparseF("fRecMCEventHist1","Xv-mcXv:Yv-mcYv:Zv-mcZv:Centrality",4,binsRecMCEventHist1,minRecMCEventHist1,maxRecMCEventHist1);
562   fRecMCEventHist1->SetBinEdges(3,fBinsCentrality);
563   fRecMCEventHist1->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
564   fRecMCEventHist1->GetAxis(1)->SetTitle("Yv-mcYv (cm)");
565   fRecMCEventHist1->GetAxis(2)->SetTitle("Zv-mcZv (cm)");
566   fRecMCEventHist1->GetAxis(3)->SetTitle("Centrality"); 
567   
568   fRecMCEventHist1->Sumw2();
569   
570   //
571   Int_t binsRecMCEventHist2[4]={100,100,fMultNbins-1, fCentralityNbins-1};
572   Double_t minRecMCEventHist2[4]={-10.0*kFact,-10.0*kFact,-0.5, 0.}; 
573   Double_t maxRecMCEventHist2[4]={10.0*kFact,10.0*kFact,10000.5, 100.}; 
574   
575   fRecMCEventHist2 = new THnSparseF("fRecMCEventHist2","Xv-mcXv:Zv-mcZv:mult:Centrality",4,binsRecMCEventHist2,minRecMCEventHist2,maxRecMCEventHist2);
576   fRecMCEventHist2->SetBinEdges(2,fBinsMult);
577   fRecMCEventHist2->SetBinEdges(3,fBinsCentrality);
578   fRecMCEventHist2->GetAxis(0)->SetTitle("Xv-mcXv (cm)");
579   fRecMCEventHist2->GetAxis(1)->SetTitle("Zv-mcZv (cm)");
580   fRecMCEventHist2->GetAxis(2)->SetTitle("multiplicity");
581   fRecMCEventHist2->GetAxis(3)->SetTitle("Centrality"); 
582   fRecMCEventHist2->Sumw2();
583   //
584   char name[256];
585   char title[256];
586   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) 
587   {
588     // THnSparse track histograms
589     
590     Int_t binsMCTrackHist1[4]=  {fPtNbins-1, fEtaNbins-1, 90, fCentralityNbins-1};
591     Double_t minMCTrackHist1[4]={0.,-1.5,0., 0.}; 
592     Double_t maxMCTrackHist1[4]={50,1.5,2.*TMath::Pi(), 100.}; 
593     snprintf(name,256,"fMCTrackHist1_%d",i);
594     snprintf(title,256,"mcPt:mcEta:mcPhi:Centrality");
595     
596     fMCTrackHist1[i] = new THnSparseF(name,title,4,binsMCTrackHist1,minMCTrackHist1,maxMCTrackHist1);
597     fMCTrackHist1[i]->SetBinEdges(0,fBinsPt);
598     fMCTrackHist1[i]->SetBinEdges(1,fBinsEta);
599     fMCTrackHist1[i]->SetBinEdges(3,fBinsCentrality);
600     fMCTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
601     fMCTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
602     fMCTrackHist1[i]->GetAxis(2)->SetTitle("mcPhi (rad)");
603     fMCTrackHist1[i]->GetAxis(3)->SetTitle("Centrality"); 
604     fMCTrackHist1[i]->Sumw2();
605     
606     Int_t binsMCPrimTrackHist1[6]=  {fPtNbins-1,fEtaNbins-1,6,20,4000, fCentralityNbins-1};
607     Double_t minMCPrimTrackHist1[6]={0.,-1.5,0.,0.,0., 0.}; 
608     Double_t maxMCPrimTrackHist1[6]={50.,1.5,6.,20.,4000., 100.}; 
609     snprintf(name,256,"fMCPrimTrackHist1_%d",i);
610     snprintf(title,256,"mcPt:mcEta:pid:mech:mother:Centrality");
611     
612     fMCPrimTrackHist1[i] = new THnSparseF(name,title,6,binsMCPrimTrackHist1,minMCPrimTrackHist1,maxMCPrimTrackHist1);
613     fMCPrimTrackHist1[i]->SetBinEdges(0,fBinsPt);
614     fMCPrimTrackHist1[i]->SetBinEdges(1,fBinsEta);
615     fMCPrimTrackHist1[i]->SetBinEdges(5,fBinsCentrality);
616     fMCPrimTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
617     fMCPrimTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
618     fMCPrimTrackHist1[i]->GetAxis(2)->SetTitle("pid");
619     fMCPrimTrackHist1[i]->GetAxis(3)->SetTitle("mech");
620     fMCPrimTrackHist1[i]->GetAxis(4)->SetTitle("mother");
621     fMCPrimTrackHist1[i]->GetAxis(5)->SetTitle("Centrality");
622     fMCPrimTrackHist1[i]->Sumw2();
623     
624     Int_t binsMCPrimTrackHist2[4]=  {4000,20,4000,fCentralityNbins-1};
625     Double_t minMCPrimTrackHist2[4]={0.,0.,0., 0.}; 
626     Double_t maxMCPrimTrackHist2[4]={4000.,20.,4000., 100.}; 
627     snprintf(name,256,"fMCPrimTrackHist2_%d",i);
628     snprintf(title,256,"pdg:mech:mother:Centrality");
629     
630     fMCPrimTrackHist2[i] = new THnSparseF(name,title,4,binsMCPrimTrackHist2,minMCPrimTrackHist2,maxMCPrimTrackHist2);
631     fMCPrimTrackHist2[i]->SetBinEdges(3,fBinsCentrality);
632     fMCPrimTrackHist2[i]->GetAxis(0)->SetTitle("pdg");
633     fMCPrimTrackHist2[i]->GetAxis(1)->SetTitle("mech");
634     fMCPrimTrackHist2[i]->GetAxis(2)->SetTitle("mother");
635     fMCPrimTrackHist2[i]->GetAxis(3)->SetTitle("Centrality");
636     fMCPrimTrackHist2[i]->Sumw2();
637     
638     Int_t binsMCSecTrackHist1[6]=  {fPtNbins-1,fEtaNbins-1,6,20,4000, fCentralityNbins-1};
639     Double_t minMCSecTrackHist1[6]={0.,-1.5,0.,0.,0., 0.}; 
640     Double_t maxMCSecTrackHist1[6]={50.,1.5,6.,20.,4000., 100.}; 
641     snprintf(name,256,"fMCSecTrackHist1_%d",i);
642     snprintf(title,256,"mcPt:mcEta:pid:mech:mother:Centrality");
643     
644     fMCSecTrackHist1[i] = new THnSparseF(name,title,6,binsMCSecTrackHist1,minMCSecTrackHist1,maxMCSecTrackHist1);
645     fMCSecTrackHist1[i]->SetBinEdges(0,fBinsPt);
646     fMCSecTrackHist1[i]->SetBinEdges(1,fBinsEta);
647     fMCSecTrackHist1[i]->SetBinEdges(5,fBinsCentrality);
648     fMCSecTrackHist1[i]->GetAxis(0)->SetTitle("mcPt (GeV/c)");
649     fMCSecTrackHist1[i]->GetAxis(1)->SetTitle("mcEta");
650     fMCSecTrackHist1[i]->GetAxis(2)->SetTitle("pid");
651     fMCSecTrackHist1[i]->GetAxis(3)->SetTitle("mech");
652     fMCSecTrackHist1[i]->GetAxis(4)->SetTitle("mother");
653     fMCSecTrackHist1[i]->GetAxis(5)->SetTitle("Centrality");
654     fMCSecTrackHist1[i]->Sumw2();
655     
656     Int_t binsRecTrackHist1[4]={fPtNbins-1,fEtaNbins-1,90, fCentralityNbins-1};
657     Double_t minRecTrackHist1[4]={0.,-1.5,0., 0.}; 
658     Double_t maxRecTrackHist1[4]={50.,1.5,2.*TMath::Pi(), 100.};
659     snprintf(name,256,"fRecTrackHist1_%d",i);
660     snprintf(title,256,"Pt:Eta:Phi:Centrality");
661     fRecTrackHist1[i] = new THnSparseF(name,title,4,binsRecTrackHist1,minRecTrackHist1,maxRecTrackHist1);
662     fRecTrackHist1[i]->SetBinEdges(0,fBinsPt);
663     fRecTrackHist1[i]->SetBinEdges(1,fBinsEta);
664     fRecTrackHist1[i]->SetBinEdges(3,fBinsCentrality);
665     fRecTrackHist1[i]->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
666     fRecTrackHist1[i]->GetAxis(1)->SetTitle("#eta");
667     fRecTrackHist1[i]->GetAxis(2)->SetTitle("#phi (rad)");
668     fRecTrackHist1[i]->GetAxis(3)->SetTitle("Centrality");
669     fRecTrackHist1[i]->Sumw2();
670     
671     snprintf(name,256,"fRecTrackHist2_%d",i);
672     snprintf(title,256,"Zv:Pt:Eta:Centrality");
673     fRecTrackHist2[i] = new THnSparseF(name,title,4,binsTrackEvent);
674     fRecTrackHist2[i]->SetBinEdges(0,fBinsZv);
675     fRecTrackHist2[i]->SetBinEdges(1,fBinsPt);
676     fRecTrackHist2[i]->SetBinEdges(2,fBinsEta);
677     fRecTrackHist2[i]->SetBinEdges(3,fBinsCentrality); 
678     fRecTrackHist2[i]->GetAxis(0)->SetTitle("Zv (cm)");
679     fRecTrackHist2[i]->GetAxis(1)->SetTitle("p_{T} (GeV/c)");
680     fRecTrackHist2[i]->GetAxis(2)->SetTitle("#eta");
681     fRecTrackHist2[i]->GetAxis(3)->SetTitle("Centrality");
682     fRecTrackHist2[i]->Sumw2();
683     
684     // 
685     Int_t binsRecTrackMultHist1[3]={fPtNbins-1,fMultNbins-1, fCentralityNbins-1};
686     Double_t minRecTrackMultHist1[3]={0.,-0.5, -0.}; 
687     Double_t maxRecTrackMultHist1[3]={50.,10000.5, 100.};
688     snprintf(name,256,"fRecTrackMultHist_%d",i);
689     snprintf(title,256,"Pt:Mult:Centrality");
690     fRecTrackMultHist1[i] = new THnSparseF(name,title,3,binsRecTrackMultHist1,minRecTrackMultHist1,maxRecTrackMultHist1);
691     fRecTrackMultHist1[i]->SetBinEdges(0,fBinsPt);
692     fRecTrackMultHist1[i]->SetBinEdges(1,fBinsMult);
693     fRecTrackMultHist1[i]->SetBinEdges(2,fBinsCentrality);   
694     fRecTrackMultHist1[i]->GetAxis(0)->SetTitle("Pt (GeV/c)");
695     fRecTrackMultHist1[i]->GetAxis(1)->SetTitle("multiplicity");
696     fRecTrackMultHist1[i]->GetAxis(2)->SetTitle("Centrality");  
697     fRecTrackMultHist1[i]->Sumw2();
698   }
699   
700   Int_t binsRecMCTrackHist1[5] = {fPtNbins-1,fEtaNbins-1,100,100, fCentralityNbins-1};
701   Double_t minRecMCTrackHist1[5]={0.,-1.5,-0.5,-0.5, 0.}; 
702   Double_t maxRecMCTrackHist1[5]={50.,1.5,0.5,0.5, 100.}; 
703   
704   snprintf(name,256,"fRecMCTrackHist1");
705   snprintf(title,256,"mcPt:mcEta:(Pt-mcPt)/mcPt:(Eta-mcEta):Centrality");
706   fRecMCTrackHist1 = new THnSparseF(name,title,5,binsRecMCTrackHist1,minRecMCTrackHist1,maxRecMCTrackHist1);
707   fRecMCTrackHist1->SetBinEdges(0,fBinsPt);
708   fRecMCTrackHist1->SetBinEdges(1,fBinsEta);
709   fRecMCTrackHist1->SetBinEdges(4,fBinsCentrality);  
710   fRecMCTrackHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");
711   fRecMCTrackHist1->GetAxis(1)->SetTitle("mcEta");
712   fRecMCTrackHist1->GetAxis(2)->SetTitle("(Pt-mcPt)/mcPt");
713   fRecMCTrackHist1->GetAxis(3)->SetTitle("Eta-mcEta");
714   fRecMCTrackHist1->GetAxis(4)->SetTitle("Centrality"); 
715   
716   Int_t binsMCMultRecTrackHist1[4] = {fPtNbins-1,fEtaNbins-1,6, fCentralityNbins-1};
717   Double_t minMCMultRecTrackHist1[4]={0.,-1.5,0., 0.}; 
718   Double_t maxMCMultRecTrackHist1[4]={50.,1.5,6., 100.}; 
719   snprintf(name,256,"fMCMultRecTrackHist1");
720   snprintf(title,256,"mcPt:mcEta:pid:Centrality");
721   fMCMultRecTrackHist1 = new THnSparseF(name,title,4,binsMCMultRecTrackHist1,minMCMultRecTrackHist1,maxMCMultRecTrackHist1);
722   fMCMultRecTrackHist1->SetBinEdges(0,fBinsPt);
723   fMCMultRecTrackHist1->SetBinEdges(1,fBinsEta);
724   fMCMultRecTrackHist1->SetBinEdges(3,fBinsCentrality);  
725   fMCMultRecTrackHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");
726   fMCMultRecTrackHist1->GetAxis(1)->SetTitle("mcEta");
727   fMCMultRecTrackHist1->GetAxis(2)->SetTitle("pid");
728   fMCMultRecTrackHist1->GetAxis(3)->SetTitle("Centrality"); 
729   
730   //nClust:chi2PerClust:pt:eta:phi:Centrality
731   Int_t binsRecTrackHist3[6]={160,100,fPtNbins-1,fEtaNbins-1,90, fCentralityNbins-1};
732   Double_t minRecTrackHist3[6]={0., 0., 0., -1.5, 0., 0.};
733   Double_t maxRecRecTrackHist3[6]={160.,10., 50., 1.5, 2.*TMath::Pi(), 100.};
734   
735   fRecTrackHist3 = new THnSparseF("fRecTrackHist3","nClust:chi2PerClust:pt:eta:phi:Centrality",6,binsRecTrackHist3,minRecTrackHist3,maxRecRecTrackHist3);
736   fRecTrackHist3->SetBinEdges(2,fBinsPt);
737   fRecTrackHist3->SetBinEdges(3,fBinsEta);
738   fRecTrackHist3->SetBinEdges(5,fBinsCentrality);  
739   fRecTrackHist3->GetAxis(0)->SetTitle("nClust");
740   fRecTrackHist3->GetAxis(1)->SetTitle("chi2PerClust");
741   fRecTrackHist3->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
742   fRecTrackHist3->GetAxis(3)->SetTitle("#eta");
743   fRecTrackHist3->GetAxis(4)->SetTitle("#phi (rad)");
744   fRecTrackHist3->GetAxis(5)->SetTitle("Centrality");
745   fRecTrackHist3->Sumw2();
746   
747   
748   // init folder
749   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");
750   
751   // init trigger analysis (for zdc cut)
752   fTriggerAnalysis = new AliTriggerAnalysis;
753   
754   // set init flag
755   fIsInit = kTRUE;
756 }
757
758 //_____________________________________________________________________________
759 void AlidNdPtAnalysisPbPb2011::Process(AliESDEvent *const esdEvent, AliMCEvent *const mcEvent)
760 {
761   //  init if not done already
762   if (!fIsInit) { Init(); }
763   //
764   // Process real and/or simulated events
765   //
766   if(!esdEvent) {
767     AliDebug(AliLog::kError, "esdEvent not available");
768     return;
769   }
770   
771   // get selection cuts
772   AlidNdPtEventCuts *evtCuts = GetEventCuts(); 
773   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); 
774   AlidNdPtAcceptanceCuts *recCuts = GetRecAcceptanceCuts();   
775   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); 
776   
777   if(!evtCuts || !accCuts  || !esdTrackCuts) {
778     AliDebug(AliLog::kError, "cuts not available");
779     return;
780   }
781   if (0 == recCuts) { recCuts = accCuts;}
782   
783   // trigger selection
784   Bool_t isEventTriggered = kTRUE;
785   AliPhysicsSelection *physicsSelection = NULL;
786   AliTriggerAnalysis* triggerAnalysis = NULL;
787   
788   // 
789   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
790   if (!inputHandler)
791   {
792     Printf("ERROR: Could not receive input handler");
793     return;
794   }
795   
796   if(evtCuts->IsTriggerRequired())  
797   {
798     // always MB
799     //isEventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
800     isEventTriggered = inputHandler->IsEventSelected() & GetTriggerMask();
801     
802     physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
803     if(!physicsSelection) return;
804     //SetPhysicsTriggerSelection(physicsSelection);
805     
806     if (isEventTriggered && (GetTrigger() == AliTriggerAnalysis::kV0AND)) {
807       // set trigger (V0AND)
808       triggerAnalysis = physicsSelection->GetTriggerAnalysis();
809       if(!triggerAnalysis) return;
810       isEventTriggered = triggerAnalysis->IsOfflineTriggerFired(esdEvent, GetTrigger());
811     }
812   }
813   
814   
815   // centrality determination
816   Float_t centralityF = -2.;
817   AliCentrality *esdCentrality = esdEvent->GetCentrality();
818   centralityF = esdCentrality->GetCentralityPercentile((const char*)GetCentralityEstimator());
819   
820   // if centrality percentile is larger then the highest bin, return
821   if(centralityF > fBinsCentrality[fCentralityNbins-1]) 
822   {
823     AliDebug(AliLog::kError, Form("Event out of centrality binning (event: %.2f, maxallowed: %.2f)",centralityF, fBinsCentrality[fCentralityNbins-1]));
824     return;
825   }
826   
827   // use MC information
828   AliHeader* header = 0;
829   AliGenEventHeader* genHeader = 0;
830   AliStack* stack = 0;
831   TArrayF vtxMC(3);
832   
833   Int_t multMCTrueTracks = 0;
834   if(IsUseMCInfo())
835   {
836     //
837     if(!mcEvent) {
838       AliDebug(AliLog::kError, "mcEvent not available");
839       return;
840     }
841     // get MC event header
842     header = mcEvent->Header();
843     if (!header) {
844       AliDebug(AliLog::kError, "Header not available");
845       return;
846     }
847     // MC particle stack
848     stack = mcEvent->Stack();
849     if (!stack) {
850       AliDebug(AliLog::kError, "Stack not available");
851       return;
852     }
853     
854     // get MC vertex
855     genHeader = header->GenEventHeader();
856     if (!genHeader) {
857       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
858       return;
859     }
860     genHeader->PrimaryVertex(vtxMC);
861     
862     Double_t vMCEventHist1[4]={vtxMC[0],vtxMC[1],vtxMC[2],centralityF};
863     fMCEventHist1->Fill(vMCEventHist1);
864     
865     // multipliticy of all MC primary tracks
866     // in Zv, pt and eta ranges)
867     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);
868     
869   } // end bUseMC
870   
871   // get reconstructed vertex  
872   const AliESDVertex* vtxESD = 0; 
873   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {
874     vtxESD = esdEvent->GetPrimaryVertexTPC();
875   }
876   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {
877     vtxESD = esdEvent->GetPrimaryVertexTracks();
878   }
879   else {
880     return;
881   }
882   
883   if(!vtxESD) return;
884   
885   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); 
886   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);
887   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());
888   
889   // vertex contributors
890   Int_t multMBTracks = 0; 
891   if(GetAnalysisMode() == AlidNdPtHelper::kTPC || GetAnalysisMode() == AlidNdPtHelper::kTPCITS) 
892   {
893     if(vtxESD->GetStatus()) {
894       multMBTracks = vtxESD->GetNContributors();
895     }
896   } 
897   else {
898     AliDebug(AliLog::kError, Form("Found analysis type %d", GetAnalysisMode()));
899     return; 
900   }
901   
902   TObjArray *allChargedTracks=0;
903   //Int_t multAll=0, multAcc=0, multRec=0;
904   Int_t multAll=0, multRec=0;
905   Int_t *labelsAll=0, *labelsAcc=0, *labelsRec=0;
906   
907   // check event cuts
908   if(isEventOK && isEventTriggered)
909   {
910     
911     // get all charged tracks
912     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());
913     if(!allChargedTracks) return;
914     
915     Int_t entries = allChargedTracks->GetEntries();
916     
917     labelsAll = new Int_t[entries];
918     labelsAcc = new Int_t[entries];
919     labelsRec = new Int_t[entries];
920     for(Int_t i=0; i<entries;++i) 
921     {
922       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
923       
924       if(!track) continue;
925       if(track->Charge()==0) continue;
926       
927       if(IsUseMCInfo() && ( mcEvent->IsFromBGEvent(track->GetLabel()) ) && (track->GetLabel() >= 0) ) 
928       {
929         // Printf("i = %i, label = %i", i, track->GetLabel());
930         continue; 
931       }
932       
933       // only postive charged 
934       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0) 
935         continue;
936       
937       // only negative charged 
938       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0) 
939         continue;
940       
941       //
942       Double_t values[4] = {vtxESD->GetZ(),track->Pt(),track->Eta(), centralityF};        
943       
944       fRecTrackHist2[AlidNdPtHelper::kAllTracks]->Fill(values);
945       FillHistograms(track,stack,AlidNdPtHelper::kAllTracks,centralityF); 
946       labelsAll[multAll] = TMath::Abs(track->GetLabel());
947       
948       multAll++;      
949       if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track) && recCuts->AcceptTrackLocalTPC(track)) {
950         
951         fRecTrackHist2[AlidNdPtHelper::kRecTracks]->Fill(values);
952         FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,centralityF); 
953         labelsRec[multRec] = TMath::Abs(track->GetLabel());
954         
955         multRec++;
956         
957       }
958     }//loop over entries
959     
960     
961     // fill track multiplicity histograms
962     //FillHistograms(allChargedTracks,labelsAll,multAll,labelsAcc,multAcc,labelsRec,multRec);
963     
964     Double_t vRecEventHist1[4] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ(),centralityF};
965     fRecEventHist1->Fill(vRecEventHist1);
966     
967     Double_t vRecEventHist2[3] = {vtxESD->GetZ(),Double_t(multMBTracks),centralityF};
968     fRecEventHist2->Fill(vRecEventHist2);
969     
970   } // triggered and event vertex
971   
972   if(IsUseMCInfo())  
973   {
974     // 
975     // event level corrections (zv,N_MB)
976     //
977     // all inelastic
978     Double_t vEventMatrix[3] = {vtxMC[2],Double_t(multMBTracks),centralityF};
979     fGenEventMatrix->Fill(vEventMatrix); 
980     if(isEventTriggered) fTriggerEventMatrix->Fill(vEventMatrix);
981     if(isEventOK && isEventTriggered) fRecEventMatrix->Fill(vEventMatrix);
982     
983     //
984     // track-event level corrections (zv,pt,eta)
985     //
986     for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) 
987     {
988       TParticle* particle = stack->Particle(iMc);
989       if (!particle)
990         continue;
991       
992       // only charged particles
993       if(!particle->GetPDG()) continue;
994       Double_t charge = particle->GetPDG()->Charge()/3.;
995       if ( TMath::Abs(charge) < 0.001 )
996         continue;
997       
998       // only postive charged 
999       if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.) 
1000         continue;
1001       
1002       // only negative charged 
1003       if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.) 
1004         continue;
1005       
1006       // physical primary
1007       Bool_t prim = stack->IsPhysicalPrimary(iMc);
1008       if(!prim) continue;
1009       
1010       // injected signal
1011       Bool_t isInjected = mcEvent->IsFromBGEvent(iMc);
1012       if(isInjected) continue;
1013       
1014       // checked accepted
1015       if(accCuts->AcceptTrack(particle)) 
1016       {
1017         Double_t vTrackEventMatrix[4] = {vtxMC[2], particle->Pt(), particle->Eta(), centralityF}; 
1018         fGenTrackEventMatrix->Fill(vTrackEventMatrix);
1019         
1020         if(!isEventTriggered) continue;  
1021         fTriggerTrackEventMatrix->Fill(vTrackEventMatrix);
1022         
1023         if(!isEventOK) continue;
1024         fRecTrackEventMatrix->Fill(vTrackEventMatrix);
1025         
1026       }// if(accCuts->AcceptTrack(particle))
1027     }// for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc)
1028     
1029     // 
1030     // track-level corrections (zv,pt,eta)
1031     //
1032     if(isEventOK && isEventTriggered)
1033     {
1034       
1035       // fill MC and rec event control histograms
1036       if(fHistogramsOn) {
1037         Double_t vRecMCEventHist1[4] = {vtxESD->GetX()-vtxMC[0],vtxESD->GetY()-vtxMC[1],vtxESD->GetZ()-vtxMC[2], centralityF};
1038         fRecMCEventHist1->Fill(vRecMCEventHist1);//
1039         
1040         Double_t vRecMCEventHist2[4] = {vtxESD->GetX()-vtxMC[0],vtxESD->GetZ()-vtxMC[2],Double_t(multMBTracks), centralityF};
1041         fRecMCEventHist2->Fill(vRecMCEventHist2);
1042         
1043       }//
1044       
1045       //
1046       // MC histograms for track efficiency studies
1047       //
1048       for (Int_t iMc = 0; iMc < stack->GetNtrack(); ++iMc) 
1049       {
1050         TParticle* particle = stack->Particle(iMc);
1051         if (!particle)
1052           continue;
1053         
1054         Double_t vTrackMatrix[4] = {vtxMC[2],particle->Pt(),particle->Eta(),centralityF}; 
1055         
1056         // only charged particles
1057         if(!particle->GetPDG()) continue;
1058         Double_t charge = particle->GetPDG()->Charge()/3.;
1059         if (TMath::Abs(charge) < 0.001)
1060           continue;
1061         
1062         // only postive charged 
1063         if(GetParticleMode() == AlidNdPtHelper::kPlus && charge < 0.) 
1064           continue;
1065         
1066         // only negative charged 
1067         if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.) 
1068           continue;
1069         
1070         // physical primary
1071         Bool_t prim = stack->IsPhysicalPrimary(iMc);
1072         
1073         // injected signal
1074         Bool_t isInjected = mcEvent->IsFromBGEvent(iMc);
1075         if(isInjected) continue;
1076         
1077         // check accepted
1078         if(accCuts->AcceptTrack(particle)) 
1079         {
1080           
1081           if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) 
1082             fGenPrimTrackMatrix->Fill(vTrackMatrix);
1083           
1084           // fill control histograms
1085           if(fHistogramsOn) 
1086             FillHistograms(stack,iMc,AlidNdPtHelper::kAccTracks, centralityF); 
1087           
1088           // check multiple found tracks
1089           Int_t multCount = 0;
1090           for(Int_t iRec=0; iRec<multRec; ++iRec)
1091           {
1092             if(iMc == labelsRec[iRec]) 
1093             {
1094               multCount++;
1095               if(multCount>1)
1096               {  
1097                 fRecMultTrackMatrix->Fill(vTrackMatrix);
1098                 
1099                 // fill control histogram
1100                 if(fHistogramsOn) {
1101                   Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1102                   Double_t vMCMultRecTrackHist1[4] = {particle->Pt(), particle->Eta(), Double_t(pid), centralityF};
1103                   fMCMultRecTrackHist1->Fill(vMCMultRecTrackHist1);
1104                 }
1105               }
1106             }
1107           }
1108           
1109           // check reconstructed
1110           for(Int_t iRec=0; iRec<multRec; ++iRec)
1111           {
1112             if(iMc == labelsRec[iRec]) 
1113             {
1114               fRecTrackMatrix->Fill(vTrackMatrix);
1115               
1116               if( AlidNdPtHelper::IsPrimaryParticle(stack, iMc, GetParticleMode()) ) {
1117                 fRecPrimTrackMatrix->Fill(vTrackMatrix);
1118               }
1119               if(!prim) fRecSecTrackMatrix->Fill(vTrackMatrix);
1120               
1121               // fill control histograms
1122               if(fHistogramsOn) 
1123                 FillHistograms(stack,iMc,AlidNdPtHelper::kRecTracks, centralityF); 
1124               
1125               break;
1126             }//if(iMc == labelsRec[iRec])
1127           }//reco tracks
1128         }//accepted tracks
1129       }//stack loop
1130     }//is triggered
1131   } // end bUseMC
1132   
1133   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;
1134   if(labelsAll) delete [] labelsAll; labelsAll = 0;
1135   if(labelsAcc) delete [] labelsAcc; labelsAcc = 0;
1136   if(labelsRec) delete [] labelsRec; labelsRec = 0;
1137   
1138 }
1139
1140 //_____________________________________________________________________________
1141 void AlidNdPtAnalysisPbPb2011::FillHistograms(TObjArray *const allChargedTracks,Int_t *const labelsAll,Int_t multAll,Int_t *const labelsAcc,Int_t multAcc,Int_t *const labelsRec,Int_t multRec, Float_t centralityF) {
1142   // multiplicity  histograms
1143   
1144   
1145   if(!allChargedTracks) return;
1146   if(!labelsAll) return;
1147   if(!labelsAcc) return;
1148   if(!labelsRec) return;
1149   
1150   Int_t entries = allChargedTracks->GetEntries();
1151   for(Int_t i=0; i<entries; ++i)
1152   {
1153     AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);
1154     if(!track) continue;
1155     if(track->Charge() == 0) continue;
1156     
1157     Int_t label = TMath::Abs(track->GetLabel());
1158     for(Int_t iAll=0; iAll<multAll; ++iAll) {
1159       if(label == labelsAll[iAll]) {
1160         Double_t v1[3] = {track->Pt(), Double_t(multAll), centralityF}; 
1161         fRecTrackMultHist1[AlidNdPtHelper::kAllTracks]->Fill(v1);
1162       }
1163     }
1164     for(Int_t iAcc=0; iAcc<multAcc; ++iAcc) {
1165       if(label == labelsAcc[iAcc]) {
1166         Double_t v2[3] = {track->Pt(), Double_t(multAcc), centralityF}; 
1167         fRecTrackMultHist1[AlidNdPtHelper::kAccTracks]->Fill(v2);
1168       }
1169     }
1170     for(Int_t iRec=0; iRec<multRec; ++iRec) {
1171       if(label == labelsRec[iRec]) {
1172         Double_t v3[3] = {track->Pt(), Double_t(multRec), centralityF}; 
1173         fRecTrackMultHist1[AlidNdPtHelper::kRecTracks]->Fill(v3);
1174       }//out
1175     }
1176   }
1177 }
1178
1179 //_____________________________________________________________________________
1180 void AlidNdPtAnalysisPbPb2011::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1181 {
1182   
1183   //
1184   // Fill ESD track and MC histograms 
1185   //
1186   if(!esdTrack) return;
1187   
1188   Float_t q = esdTrack->Charge();
1189   if(TMath::Abs(q) < 0.001) return;
1190   
1191   Float_t pt = esdTrack->Pt();
1192   Float_t eta = esdTrack->Eta();
1193   Float_t phi = esdTrack->Phi();
1194   
1195   Float_t dca[2], bCov[3];
1196   esdTrack->GetImpactParameters(dca,bCov);
1197   
1198   Int_t nClust = esdTrack->GetTPCclusters(0);
1199   Float_t chi2PerCluster = 0.;
1200   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
1201   
1202   
1203   // fill histograms
1204   Double_t values[4] = {pt,eta,phi,centralityF};          
1205   fRecTrackHist1[trackObj]->Fill(values);
1206   
1207   //
1208   // Fill rec vs MC information
1209   //
1210   if(!stack) return;
1211   
1212   Int_t label = TMath::Abs(esdTrack->GetLabel()); 
1213   //if(label == 0) return;
1214   
1215   if(label > stack->GetNtrack()) return;
1216   TParticle* particle = stack->Particle(label);
1217   if(!particle) return;
1218   
1219   Int_t motherPdg = -1;
1220   TParticle* mother = 0;
1221   
1222   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1223   Int_t motherLabel = particle->GetMother(0); 
1224   if(motherLabel>0) mother = stack->Particle(motherLabel);
1225   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1226     //Int_t mech = particle->GetUniqueID(); // production mechanism
1227     
1228     if(!particle->GetPDG()) return;
1229     Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3 
1230     if(TMath::Abs(gq)<0.001) return;
1231     Float_t gpt = particle->Pt();
1232   Float_t geta = particle->Eta();
1233   
1234   Double_t dpt=0;
1235   //printf("pt %f, gpt %f \n",pt,gpt);
1236   if(gpt) dpt = (pt-gpt)/gpt;
1237   Double_t deta = (eta-geta);
1238   
1239   // fill histograms
1240   if(trackObj == AlidNdPtHelper::kRecTracks)  //RecTracks???
1241   {
1242     Double_t vTrackPtCorrelationMatrix[4]={pt,gpt,geta,centralityF};
1243     fTrackPtCorrelationMatrix->Fill(vTrackPtCorrelationMatrix);
1244     
1245     Double_t vRecMCTrackHist1[5]={gpt,geta,dpt,deta,centralityF};
1246     fRecMCTrackHist1->Fill(vRecMCTrackHist1);
1247   }
1248 }
1249
1250 //_____________________________________________________________________________
1251 void AlidNdPtAnalysisPbPb2011::FillHistograms(AliStack *const stack, Int_t label, AlidNdPtHelper::TrackObject trackObj, Float_t centralityF)
1252 {
1253   
1254   // Fill MC histograms
1255   if(!stack) return;
1256   
1257   if(label > stack->GetNtrack()) return;
1258   TParticle* particle = stack->Particle(label);
1259   if(!particle) return;
1260   
1261   Int_t motherPdg = -1;
1262   TParticle* mother = 0;
1263   
1264   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);
1265   Int_t motherLabel = particle->GetMother(0); 
1266   if(motherLabel>0) mother = stack->Particle(motherLabel);
1267   if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only
1268     Int_t mech = particle->GetUniqueID(); // production mechanism
1269     
1270     if(!particle->GetPDG()) return;
1271     Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3 
1272     if(TMath::Abs(gq) < 0.001) return;
1273     
1274     Float_t gpt = particle->Pt();
1275   //Float_t qgpt = particle->Pt() * gq;
1276   Float_t geta = particle->Eta();
1277   Float_t gphi = particle->Phi();
1278   //Float_t gpz = particle->Pz();
1279   
1280   Bool_t prim = stack->IsPhysicalPrimary(label);
1281   //Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();
1282   
1283   Int_t pid = AlidNdPtHelper::ConvertPdgToPid(particle);
1284   
1285   //if(prim&&pid==5) printf("pdgcode %d, production mech %d \n",particle->GetPdgCode(),mech);
1286   //if(!prim) printf("motherPdg %d, particle %d, production mech %d\n",motherPdg, particle->GetPdgCode(),mech);
1287   
1288   //
1289   // fill histogram
1290   //
1291   Double_t vMCTrackHist1[4] = {gpt,geta,gphi,centralityF};
1292   fMCTrackHist1[trackObj]->Fill(vMCTrackHist1);
1293   
1294   Double_t vMCPrimTrackHist1[6] = {gpt,geta,Double_t(pid),Double_t(mech),Double_t(motherPdg),centralityF};
1295   Double_t vMCPrimTrackHist2[4] = {Double_t(TMath::Abs(particle->GetPdgCode())),Double_t(mech),Double_t(motherPdg),centralityF};
1296   
1297   //if(prim && AliPWG0Helper::IsPrimaryCharged(particle, label)) fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1298   
1299   if(prim) { 
1300     fMCPrimTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1301     if(pid == 5) fMCPrimTrackHist2[trackObj]->Fill(vMCPrimTrackHist2);
1302   }
1303   else { 
1304     fMCSecTrackHist1[trackObj]->Fill(vMCPrimTrackHist1);
1305   }
1306   
1307 }
1308
1309 //_____________________________________________________________________________
1310 Long64_t AlidNdPtAnalysisPbPb2011::Merge(TCollection* const list) 
1311 {
1312   // Merge list of objects (needed by PROOF)
1313   
1314   //  init if not done already
1315   if (!fIsInit) { Init(); }
1316   
1317   if (!list)
1318     return 0;
1319   
1320   if (list->IsEmpty())
1321     return 1;
1322   
1323   TIterator* iter = list->MakeIterator();
1324   TObject* obj = 0;
1325   
1326   //
1327   //TList *collPhysSelection = new TList;
1328   
1329   // collection of generated histograms
1330   
1331   Int_t count=0;
1332   while((obj = iter->Next()) != 0) {
1333     AlidNdPtAnalysisPbPb2011* entry = dynamic_cast<AlidNdPtAnalysisPbPb2011*>(obj);
1334     if (entry == 0) continue; 
1335     
1336     // physics selection
1337     //printf("entry->GetPhysicsTriggerSelection() %p \n", entry->GetPhysicsTriggerSelection());
1338     //AliPhysicsSelection *physSel = entry->GetPhysicsTriggerSelection();
1339     //if( physSel ) collPhysSelection->Add(physSel); 
1340     
1341     //
1342     fTrackPtCorrelationMatrix->Add(entry->fTrackPtCorrelationMatrix);
1343     
1344     //
1345     fGenEventMatrix->Add(entry->fGenEventMatrix);
1346     
1347     fTriggerEventMatrix->Add(entry->fTriggerEventMatrix);
1348     
1349     fRecEventMatrix->Add(entry->fRecEventMatrix);
1350     //
1351     fGenTrackEventMatrix->Add(entry->fGenTrackEventMatrix);
1352     
1353     fTriggerTrackEventMatrix->Add(entry->fTriggerTrackEventMatrix);
1354     
1355     fRecTrackEventMatrix->Add(entry->fRecTrackEventMatrix);
1356     
1357     //
1358     fGenTrackMatrix->Add(entry->fGenTrackMatrix);
1359     fGenPrimTrackMatrix->Add(entry->fGenPrimTrackMatrix);
1360     fRecPrimTrackMatrix->Add(entry->fRecPrimTrackMatrix);
1361     //
1362     fRecTrackMatrix->Add(entry->fRecTrackMatrix);
1363     fRecSecTrackMatrix->Add(entry->fRecSecTrackMatrix);
1364     //
1365     fRecMultTrackMatrix->Add(entry->fRecMultTrackMatrix);
1366     
1367     //
1368     // control analysis histograms
1369     //
1370     fMCEventHist1->Add(entry->fMCEventHist1);
1371     fRecEventHist1->Add(entry->fRecEventHist1);
1372     fRecEventHist2->Add(entry->fRecEventHist2);
1373     fRecMCEventHist1->Add(entry->fRecMCEventHist1);
1374     fRecMCEventHist2->Add(entry->fRecMCEventHist2);
1375     
1376     
1377     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {
1378       fMCTrackHist1[i]->Add(entry->fMCTrackHist1[i]);
1379       
1380       fMCPrimTrackHist1[i]->Add(entry->fMCPrimTrackHist1[i]);
1381       fMCPrimTrackHist2[i]->Add(entry->fMCPrimTrackHist2[i]);
1382       fMCSecTrackHist1[i]->Add(entry->fMCSecTrackHist1[i]);
1383       
1384       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);
1385       fRecTrackHist2[i]->Add(entry->fRecTrackHist2[i]);
1386       fRecTrackMultHist1[i]->Add(entry->fRecTrackMultHist1[i]);
1387     }
1388     fRecMCTrackHist1->Add(entry->fRecMCTrackHist1);
1389     fMCMultRecTrackHist1->Add(entry->fMCMultRecTrackHist1);
1390     fRecTrackHist3->Add(entry->fRecTrackHist3);
1391     
1392     count++;
1393   }
1394   
1395   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();
1396   //if( trigSelection ) trigSelection->Merge(collPhysSelection);
1397   //if(collPhysSelection) delete collPhysSelection;
1398   
1399   return count;
1400 }
1401
1402
1403
1404 //_____________________________________________________________________________
1405 void AlidNdPtAnalysisPbPb2011::Analyse() 
1406 {
1407   
1408   //  init if not done already
1409   if (!fIsInit) { Init(); }
1410   
1411   // Analyse histograms
1412   //
1413   TH1::AddDirectory(kFALSE);
1414   TH1 *h=0, *h1=0, *h2=0, *h2c = 0; 
1415   THnSparse *hs=0; 
1416   TH2 *h2D=0; 
1417   
1418   char name[256];
1419   TObjArray *aFolderObj = new TObjArray;
1420   if(!aFolderObj) return;
1421   
1422   //
1423   // LHC backgraund in all and 0-bins
1424   //
1425   //AliPhysicsSelection *trigSel = GetPhysicsTriggerSelection();
1426   //trigSel->SaveHistograms("physics_selection");
1427   
1428   //
1429   // Reconstructed event vertex
1430   //
1431   h = fRecEventHist1->Projection(0);
1432   h->SetName("Xv");
1433   aFolderObj->Add(h);
1434   
1435   h = fRecEventHist1->Projection(1);
1436   h->SetName("Yv");
1437   aFolderObj->Add(h);
1438   
1439   h = fRecEventHist1->Projection(2);
1440   h->SetName("Zv");
1441   aFolderObj->Add(h);
1442   
1443   //
1444   // multiplicity
1445   //
1446   h = fRecEventHist2->Projection(1);
1447   h->SetName("multMB");
1448   aFolderObj->Add(h);
1449   
1450   h2D = fRecEventHist2->Projection(0,1); 
1451   h2D->SetName("Zv_vs_multiplicity_MB");
1452   aFolderObj->Add(h2D);
1453   
1454   //
1455   // reconstructed pt histograms
1456   //
1457   h = fRecTrackHist1[0]->Projection(0);
1458   h->Scale(1.,"width");
1459   h->SetName("pt_all_ch");
1460   aFolderObj->Add(h);
1461   
1462   h = fRecTrackHist1[1]->Projection(0);
1463   h->Scale(1.,"width");
1464   h->SetName("pt_acc");
1465   aFolderObj->Add(h);
1466   
1467   h = fRecTrackHist1[2]->Projection(0);
1468   h->Scale(1.,"width");
1469   h->SetName("pt_rec");
1470   aFolderObj->Add(h);
1471   
1472   //
1473   // reconstructed eta histograms
1474   //
1475   h = fRecTrackHist1[0]->Projection(1);
1476   h->SetName("eta_all_ch");
1477   aFolderObj->Add(h);
1478   
1479   h = fRecTrackHist1[1]->Projection(1);
1480   h->SetName("eta_acc");
1481   aFolderObj->Add(h);
1482   
1483   h = fRecTrackHist1[2]->Projection(1);
1484   h->SetName("eta_rec");
1485   aFolderObj->Add(h);
1486   
1487   //
1488   // reconstructed phi histograms
1489   //
1490   h = fRecTrackHist1[0]->Projection(2);
1491   h->SetName("phi_all_ch");
1492   aFolderObj->Add(h);
1493   
1494   h = fRecTrackHist1[1]->Projection(2);
1495   h->SetName("phi_acc");
1496   aFolderObj->Add(h);
1497   
1498   h = fRecTrackHist1[2]->Projection(2);
1499   h->SetName("phi_rec");
1500   aFolderObj->Add(h);
1501   
1502   //
1503   // reconstructed eta:pt histograms
1504   //
1505   h2D = fRecTrackHist1[0]->Projection(1,0);
1506   h2D->SetName("pt_eta_all_ch");
1507   aFolderObj->Add(h2D);
1508   
1509   h2D = fRecTrackHist1[1]->Projection(1,0);
1510   h2D->SetName("pt_eta_acc");
1511   aFolderObj->Add(h2D);
1512   
1513   h2D = fRecTrackHist1[2]->Projection(1,0);
1514   h2D->SetName("pt_eta_rec");
1515   aFolderObj->Add(h2D);
1516   
1517   //
1518   // reconstructed phi:pt histograms
1519   //
1520   h2D = fRecTrackHist1[0]->Projection(2,0);
1521   h2D->SetName("pt_phi_all_ch");
1522   aFolderObj->Add(h2D);
1523   
1524   h2D = fRecTrackHist1[1]->Projection(2,0);
1525   h2D->SetName("pt_phi_acc");
1526   aFolderObj->Add(h2D);
1527   
1528   h2D = fRecTrackHist1[2]->Projection(2,0);
1529   h2D->SetName("pt_phi_rec");
1530   aFolderObj->Add(h2D);
1531   
1532   //
1533   // reconstructed phi:eta histograms
1534   //
1535   h2D = fRecTrackHist1[0]->Projection(2,1);
1536   h2D->SetName("eta_phi_all_ch");
1537   aFolderObj->Add(h2D);
1538   
1539   h2D = fRecTrackHist1[1]->Projection(2,1);
1540   h2D->SetName("eta_phi_acc");
1541   aFolderObj->Add(h2D);
1542   
1543   h2D = fRecTrackHist1[2]->Projection(2,1);
1544   h2D->SetName("eta_phi_rec");
1545   aFolderObj->Add(h2D);
1546   
1547   //
1548   // reconstructed nClust, chi2 vs pt, eta, phi
1549   //
1550   if(fHistogramsOn) {
1551     
1552     h2D = fRecTrackHist3->Projection(0,1);
1553     h2D->SetName("nClust_chi2_rec");
1554     aFolderObj->Add(h2D);
1555     
1556     h2D = fRecTrackHist3->Projection(0,2);
1557     h2D->SetName("nClust_pt_rec");
1558     aFolderObj->Add(h2D);
1559     
1560     h2D = fRecTrackHist3->Projection(0,3);
1561     h2D->SetName("nClust_eta_rec");
1562     aFolderObj->Add(h2D);
1563     
1564     h2D = fRecTrackHist3->Projection(0,4);
1565     h2D->SetName("nClust_phi_rec");
1566     aFolderObj->Add(h2D);
1567     
1568     h2D = fRecTrackHist3->Projection(1,2);
1569     h2D->SetName("chi2_pt_rec");
1570     aFolderObj->Add(h2D);
1571     
1572     h2D = fRecTrackHist3->Projection(1,3);
1573     h2D->SetName("chi2_eta_rec");
1574     aFolderObj->Add(h2D);
1575     
1576     h2D = fRecTrackHist3->Projection(1,4);
1577     h2D->SetName("chi2_phi_rec");
1578     aFolderObj->Add(h2D);
1579     
1580   }
1581   
1582   //
1583   // calculate corrections for empty events
1584   // with multMB==0 
1585   //
1586   
1587   //
1588   // normalised zv to generate zv for triggered events
1589   //
1590   h = fRecEventHist2->Projection(0);
1591   if( h->Integral() ) h->Scale(1./h->Integral());
1592   h->SetName("zv_distribution_norm");
1593   aFolderObj->Add(h);
1594   
1595   //
1596   // MC available
1597   //
1598   if(IsUseMCInfo()) {
1599     
1600     //
1601     // Event vertex resolution
1602     //
1603     h2D = fRecMCEventHist2->Projection(0,2);
1604     h2D->SetName("DeltaXv_vs_mult");
1605     aFolderObj->Add(h2D);
1606     
1607     h2D = fRecMCEventHist2->Projection(1,2);
1608     h2D->SetName("DeltaZv_vs_mult");
1609     aFolderObj->Add(h2D);
1610     
1611     //
1612     // normalised zv to get trigger/trigger+vertex event differences
1613     // F(zv) = E_trig(zv,0)/Int(E_trig(zv,0) / Sum(E_trigvtx(zv,n))/Sum(Int(E_trigvtx(zv,n))dzv)
1614     //
1615     fTriggerEventMatrix->GetAxis(1)->SetRangeUser(0.,0.);
1616     h = fTriggerEventMatrix->Projection(0);
1617     h2D = fTriggerEventMatrix->Projection(0,1);
1618     if(h2D->Integral()) h->Scale(1./h2D->Integral());
1619     
1620     h1 = fRecEventMatrix->Projection(0);
1621     h2D = fRecEventMatrix->Projection(0,1);
1622     if(h2D->Integral()) h1->Scale(1./h2D->Integral());
1623     
1624     h->Divide(h1);
1625     h->SetName("zv_empty_events_norm");
1626     aFolderObj->Add(h);
1627     
1628     fTriggerEventMatrix->GetAxis(1)->SetRange(1,fTriggerEventMatrix->GetAxis(1)->GetNbins());
1629     
1630     //
1631     // rec. vs true track pt correlation matrix
1632     //
1633     hs = (THnSparse*)fTrackPtCorrelationMatrix->Clone("track_pt_correlation_matrix");
1634     aFolderObj->Add(hs);
1635     
1636     //
1637     // trigger efficiency for INEL
1638     //
1639     h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fGenEventMatrix->Projection(0),"zv_trig_INEL_eff_matrix");
1640     aFolderObj->Add(h);
1641     
1642     
1643     //
1644     // trigger bias correction (MB to INEL)
1645     //
1646     hs = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix,fTriggerEventMatrix,"zv_mult_trig_MBtoInel_corr_matrix");
1647     aFolderObj->Add(hs);
1648     
1649     h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0),fTriggerEventMatrix->Projection(0),"zv_trig_MBtoInel_corr_matrix");
1650     aFolderObj->Add(h);
1651     
1652     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(0,1),fTriggerEventMatrix->Projection(0,1),"zv_mult_trig_MBtoInel_corr_matrix_2D");
1653     aFolderObj->Add(h2D);
1654     
1655     
1656     h = AlidNdPtHelper::GenerateCorrMatrix(fGenEventMatrix->Projection(1),fTriggerEventMatrix->Projection(1),"mult_trig_MBtoInel_corr_matrix");
1657     aFolderObj->Add(h);
1658     
1659     
1660     //
1661     // event vertex reconstruction correction (MB)
1662     //
1663     hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix,fRecEventMatrix,"zv_mult_event_corr_matrix");
1664     aFolderObj->Add(hs);
1665     
1666     h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0,1),fRecEventMatrix->Projection(0,1),"zv_mult_event_corr_matrix_2D");
1667     aFolderObj->Add(h2D);
1668     
1669     
1670     h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(1),fRecEventMatrix->Projection(1),"mult_event_corr_matrix");
1671     aFolderObj->Add(h);
1672     
1673     
1674     h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerEventMatrix->Projection(0),fRecEventMatrix->Projection(0),"zv_event_corr_matrix");
1675     aFolderObj->Add(h);
1676     
1677     
1678     
1679     //
1680     // track-event trigger bias correction (MB to INEL)
1681     //
1682     hs = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix,fTriggerTrackEventMatrix,"zv_pt_eta_track_trig_MBtoInel_corr_matrix");
1683     aFolderObj->Add(hs);
1684     
1685     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,2),fTriggerTrackEventMatrix->Projection(1,2),"eta_pt_track_trig_MBtoInel_corr_matrix");
1686     aFolderObj->Add(h2D);
1687     
1688     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(1,0),fTriggerTrackEventMatrix->Projection(1,0),"pt_zv_track_trig_MBtoInel_corr_matrix");
1689     aFolderObj->Add(h2D);
1690     
1691     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenTrackEventMatrix->Projection(2,0),fTriggerTrackEventMatrix->Projection(2,0),"zv_eta_track_trig_MBtoInel_corr_matrix");
1692     aFolderObj->Add(h2D);
1693     
1694     // efficiency
1695     
1696     h = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1),fGenTrackEventMatrix->Projection(1),"pt_track_trig_MBtoInel_eff_matrix");
1697     aFolderObj->Add(h);
1698     
1699     
1700     //
1701     // track-event vertex reconstruction correction (MB)
1702     //
1703     hs = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix,fRecTrackEventMatrix,"zv_pt_eta_track_event_corr_matrix");
1704     aFolderObj->Add(hs);
1705     
1706     h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,2),fRecTrackEventMatrix->Projection(1,2),"eta_pt_track_event_corr_matrix");
1707     aFolderObj->Add(h2D);
1708     
1709     h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(1,0),fRecTrackEventMatrix->Projection(1,0),"pt_zv_track_event_corr_matrix");
1710     aFolderObj->Add(h2D);
1711     
1712     h2D = AlidNdPtHelper::GenerateCorrMatrix(fTriggerTrackEventMatrix->Projection(2,0),fRecTrackEventMatrix->Projection(2,0),"zv_eta_track_event_corr_matrix");
1713     aFolderObj->Add(h2D);
1714     
1715     // efficiency
1716     
1717     h = AlidNdPtHelper::GenerateCorrMatrix(fRecTrackEventMatrix->Projection(1),fTriggerTrackEventMatrix->Projection(1),"pt_track_event_eff_matrix");
1718     aFolderObj->Add(h);
1719     
1720     
1721     //
1722     // track rec. efficiency correction
1723     //
1724     hs = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix,fRecPrimTrackMatrix,"zv_pt_eta_track_corr_matrix");
1725     aFolderObj->Add(hs);
1726     
1727     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,2),fRecPrimTrackMatrix->Projection(1,2),"eta_pt_track_corr_matrix");
1728     aFolderObj->Add(h2D);
1729     
1730     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1,0),fRecPrimTrackMatrix->Projection(1,0),"pt_zv_track_corr_matrix");
1731     aFolderObj->Add(h2D);
1732     
1733     h2D = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2,0),fRecPrimTrackMatrix->Projection(2,0),"zv_eta_track_corr_matrix");
1734     aFolderObj->Add(h2D);
1735     
1736     
1737     h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(0),fRecPrimTrackMatrix->Projection(0),"zv_track_corr_matrix");
1738     aFolderObj->Add(h);
1739     
1740     h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(1),fRecPrimTrackMatrix->Projection(1),"pt_track_corr_matrix");
1741     aFolderObj->Add(h);
1742     
1743     // efficiency
1744     
1745     h = AlidNdPtHelper::GenerateCorrMatrix(fRecPrimTrackMatrix->Projection(1), fGenPrimTrackMatrix->Projection(1),"pt_track_eff_matrix");
1746     aFolderObj->Add(h);
1747     
1748     h = AlidNdPtHelper::GenerateCorrMatrix(fGenPrimTrackMatrix->Projection(2),fRecPrimTrackMatrix->Projection(2),"eta_track_corr_matrix");
1749     aFolderObj->Add(h);
1750     
1751     //
1752     // secondary track contamination correction
1753     //
1754     //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1755     hs = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix,fRecTrackMatrix,"zv_pt_eta_track_cont_matrix");
1756     aFolderObj->Add(hs);
1757     
1758     h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_track_cont_matrix");
1759     aFolderObj->Add(h2D);
1760     
1761     h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_track_cont_matrix");
1762     aFolderObj->Add(h2D);
1763     
1764     h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_track_cont_matrix");
1765     aFolderObj->Add(h2D);
1766     
1767     h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_track_cont_matrix");
1768     aFolderObj->Add(h);
1769     
1770     
1771     h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_track_cont_matrix");
1772     aFolderObj->Add(h);
1773     
1774     h = AlidNdPtHelper::GenerateCorrMatrix(fRecSecTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_track_cont_matrix");
1775     aFolderObj->Add(h);
1776     
1777     //
1778     // multiple track reconstruction correction
1779     //
1780     //hs = AlidNdPtHelper::GenerateContCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1781     hs = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix,fRecTrackMatrix,"zv_pt_eta_mult_track_cont_matrix");
1782     aFolderObj->Add(hs);
1783     
1784     h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,2),fRecTrackMatrix->Projection(1,2),"eta_pt_mult_track_cont_matrix");
1785     aFolderObj->Add(h2D);
1786     
1787     h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1,0),fRecTrackMatrix->Projection(1,0),"pt_zv_mult_track_cont_matrix");
1788     aFolderObj->Add(h2D);
1789     
1790     h2D = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2,0),fRecTrackMatrix->Projection(2,0),"zv_eta_mult_track_cont_matrix");
1791     aFolderObj->Add(h2D);
1792     
1793     h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(0),fRecTrackMatrix->Projection(0),"zv_mult_track_cont_matrix");
1794     aFolderObj->Add(h);
1795     
1796     h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(1),fRecTrackMatrix->Projection(1),"pt_mult_track_cont_matrix");
1797     aFolderObj->Add(h);
1798     
1799     h = AlidNdPtHelper::GenerateCorrMatrix(fRecMultTrackMatrix->Projection(2),fRecTrackMatrix->Projection(2),"eta_mult_track_cont_matrix");
1800     aFolderObj->Add(h);
1801     
1802     //
1803     // Control histograms
1804     //
1805     
1806     if(fHistogramsOn) {
1807       
1808       // Efficiency electrons, muons, pions, kaons, protons, all
1809       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,1); 
1810       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,1); 
1811       h1 = fMCPrimTrackHist1[1]->Projection(0);
1812       h2 = fMCPrimTrackHist1[2]->Projection(0);
1813       h2c = (TH1D *)h2->Clone();
1814       h2c->Divide(h1);
1815       h2c->SetName("eff_pt_electrons");
1816       aFolderObj->Add(h2c);
1817       
1818       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(2,2); 
1819       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(2,2); 
1820       h1 = fMCPrimTrackHist1[1]->Projection(0);
1821       h2 = fMCPrimTrackHist1[2]->Projection(0);
1822       h2c = (TH1D *)h2->Clone();
1823       h2c->Divide(h1);
1824       h2c->SetName("eff_pt_muons");
1825       aFolderObj->Add(h2c);
1826       
1827       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(3,3); 
1828       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(3,3); 
1829       h1 = fMCPrimTrackHist1[1]->Projection(0);
1830       h2 = fMCPrimTrackHist1[2]->Projection(0);
1831       h2c = (TH1D *)h2->Clone();
1832       h2c->Divide(h1);
1833       h2c->SetName("eff_pt_pions");
1834       aFolderObj->Add(h2c);
1835       
1836       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(4,4); 
1837       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(4,4); 
1838       h1 = fMCPrimTrackHist1[1]->Projection(0);
1839       h2 = fMCPrimTrackHist1[2]->Projection(0);
1840       h2c = (TH1D *)h2->Clone();
1841       h2c->Divide(h1);
1842       h2c->SetName("eff_pt_kaons");
1843       aFolderObj->Add(h2c);
1844       
1845       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(5,5); 
1846       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(5,5); 
1847       h1 = fMCPrimTrackHist1[1]->Projection(0);
1848       h2 = fMCPrimTrackHist1[2]->Projection(0);
1849       h2c = (TH1D *)h2->Clone();
1850       h2c->Divide(h1);
1851       h2c->SetName("eff_pt_protons");
1852       aFolderObj->Add(h2c);
1853       
1854       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,5); 
1855       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,5); 
1856       h1 = fMCPrimTrackHist1[1]->Projection(0);
1857       h2 = fMCPrimTrackHist1[2]->Projection(0);
1858       h2c = (TH1D *)h2->Clone();
1859       h2c->Divide(h1);
1860       h2c->SetName("eff_pt_selected");
1861       aFolderObj->Add(h2c);
1862       
1863       fMCPrimTrackHist1[1]->GetAxis(2)->SetRange(1,6); 
1864       fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(1,6); 
1865       h1 = fMCPrimTrackHist1[1]->Projection(0);
1866       h2 = fMCPrimTrackHist1[2]->Projection(0);
1867       h2c = (TH1D *)h2->Clone();
1868       h2c->Divide(h1);
1869       h2c->SetName("eff_pt_all");
1870       aFolderObj->Add(h2c);
1871       
1872       fMCPrimTrackHist1[1]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[1]->GetAxis(1)->GetNbins()); 
1873       fMCPrimTrackHist1[2]->GetAxis(1)->SetRange(1,fMCPrimTrackHist1[2]->GetAxis(1)->GetNbins());
1874       
1875       //  pt spetra
1876       // - rec, primaries, secondaries
1877       // - primaries (pid) 
1878       // - secondaries (pid)
1879       // - secondaries (mech)
1880       // - secondaries (mother)
1881       //
1882       
1883       TH1D *mcPtAccall = fMCTrackHist1[1]->Projection(0);
1884       mcPtAccall->SetName("mc_pt_acc_all");
1885       aFolderObj->Add(mcPtAccall);
1886       
1887       TH1D *mcPtAccprim = fMCPrimTrackHist1[1]->Projection(0);
1888       mcPtAccprim->SetName("mc_pt_acc_prim");
1889       aFolderObj->Add(mcPtAccprim);
1890       
1891       TH1D *mcPtRecall = fMCTrackHist1[2]->Projection(0);
1892       mcPtRecall->SetName("mc_pt_rec_all");
1893       aFolderObj->Add(mcPtRecall);
1894       
1895       TH1D *mcPtRecprim = fMCPrimTrackHist1[2]->Projection(0);
1896       mcPtRecprim->SetName("mc_pt_rec_prim");
1897       aFolderObj->Add(mcPtRecprim);
1898       
1899       TH1D *mcPtRecsec = fMCSecTrackHist1[2]->Projection(0);
1900       mcPtRecsec->SetName("mc_pt_rec_sec");
1901       aFolderObj->Add(mcPtRecsec);
1902       
1903       for(Int_t i = 0; i<6; i++) 
1904       { 
1905         snprintf(name,256,"mc_pt_rec_prim_pid_%d",i); 
1906         fMCPrimTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1907         h = fMCPrimTrackHist1[2]->Projection(0);
1908         h->SetName(name);
1909         aFolderObj->Add(h);
1910         
1911         snprintf(name,256,"mc_pt_rec_sec_pid_%d",i); 
1912         fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1913         h = fMCSecTrackHist1[2]->Projection(0);
1914         h->SetName(name);
1915         aFolderObj->Add(h);
1916         
1917         // production mechanisms for given pid
1918         fMCSecTrackHist1[2]->GetAxis(2)->SetRange(i+1,i+1);
1919         
1920         for(Int_t j=0; j<20; j++) {
1921           if(j == 4) {
1922             // decay
1923             
1924             snprintf(name,256,"mc_pt_rec_sec_pid_%d_decay",i); 
1925             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1926             h = fMCSecTrackHist1[2]->Projection(0);
1927             h->SetName(name);
1928             aFolderObj->Add(h);
1929             
1930             snprintf(name,256,"mc_eta_rec_sec_pid_%d_decay",i); 
1931             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1932             h = fMCSecTrackHist1[2]->Projection(1);
1933             h->SetName(name);
1934             aFolderObj->Add(h);
1935             
1936             snprintf(name,256,"mc_mother_rec_sec_pid_%d_decay",i); 
1937             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1938             h = fMCSecTrackHist1[2]->Projection(4);
1939             h->SetName(name);
1940             aFolderObj->Add(h);
1941             
1942           } else if (j == 5) {
1943             // conversion
1944             
1945             snprintf(name,256,"mc_pt_rec_sec_pid_%d_conv",i); 
1946             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1947             h = fMCSecTrackHist1[2]->Projection(0);
1948             h->SetName(name);
1949             aFolderObj->Add(h);
1950             
1951             snprintf(name,256,"mc_eta_rec_sec_pid_%d_conv",i); 
1952             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1953             h = fMCSecTrackHist1[2]->Projection(1);
1954             h->SetName(name);
1955             aFolderObj->Add(h);
1956             
1957             snprintf(name,256,"mc_mother_rec_sec_pid_%d_conv",i); 
1958             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1959             h = fMCSecTrackHist1[2]->Projection(4);
1960             h->SetName(name);
1961             aFolderObj->Add(h);
1962             
1963           } else if (j == 13) {
1964             // mat
1965             
1966             snprintf(name,256,"mc_pt_rec_sec_pid_%d_mat",i); 
1967             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1968             h = fMCSecTrackHist1[2]->Projection(0);
1969             h->SetName(name);
1970             aFolderObj->Add(h);
1971             
1972             snprintf(name,256,"mc_eta_rec_sec_pid_%d_mat",i); 
1973             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1974             h = fMCSecTrackHist1[2]->Projection(1);
1975             h->SetName(name);
1976             aFolderObj->Add(h);
1977             
1978             snprintf(name,256,"mc_eta_mother_rec_sec_pid_%d_mat",i); 
1979             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1980             h = fMCSecTrackHist1[2]->Projection(4,1);
1981             h->SetName(name);
1982             aFolderObj->Add(h);
1983             
1984             snprintf(name,256,"mc_mother_rec_sec_pid_%d_mat",i); 
1985             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1986             h = fMCSecTrackHist1[2]->Projection(4);
1987             h->SetName(name);
1988             aFolderObj->Add(h);
1989             
1990             snprintf(name,256,"mc_pt_mother_rec_sec_pid_%d_mat",i); 
1991             fMCSecTrackHist1[2]->GetAxis(3)->SetRange(j+1,j+1);
1992             h = fMCSecTrackHist1[2]->Projection(4,0);
1993             h->SetName(name);
1994             aFolderObj->Add(h);
1995             
1996           } else {
1997             continue;
1998           }
1999         }
2000         
2001       }
2002     } // end fHistogramOn
2003     
2004     //
2005     //  resolution histograms
2006     //  only for reconstructed tracks
2007     //
2008     
2009     TH2F *h2F=0;
2010     TCanvas * c = new TCanvas("resol","resol");
2011     c->cd();
2012     
2013     //
2014     fRecMCTrackHist1->GetAxis(1)->SetRangeUser(-0.8,0.79); 
2015     
2016     h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
2017     h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2018     h->SetXTitle("p_{tmc} (GeV/c)");
2019     h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
2020     h->Draw();
2021     h->SetName("pt_resolution_vs_mcpt");
2022     aFolderObj->Add(h);
2023     
2024     h2F = (TH2F*)fRecMCTrackHist1->Projection(2,0);
2025     h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2026     h->SetXTitle("p_{tmc} (GeV/c)");
2027     h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2028     h->Draw();
2029     h->SetName("dpt_mean_vs_mcpt");
2030     aFolderObj->Add(h);
2031     
2032     //
2033     h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
2034     h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2035     h->SetXTitle("p_{tmc} (GeV/c)");
2036     h->SetYTitle("(#eta-#eta_{mc}) resolution");
2037     h->Draw();
2038     h->SetName("eta_resolution_vs_mcpt");
2039     aFolderObj->Add(h);
2040     
2041     h2F = (TH2F*)fRecMCTrackHist1->Projection(3,0);
2042     h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2043     h->SetXTitle("p_{tmc} (GeV/c)");
2044     h->SetYTitle("(#eta-mc#eta) mean");
2045     h->Draw();
2046     h->SetName("deta_mean_vs_mcpt");
2047     aFolderObj->Add(h);
2048     
2049     // 
2050     fRecMCTrackHist1->GetAxis(1)->SetRange(1,fRecMCTrackHist1->GetAxis(1)->GetNbins()); 
2051     
2052     h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2053     h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2054     h->SetXTitle("#eta_{mc}");
2055     h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} resolution");
2056     h->Draw();
2057     h->SetName("pt_resolution_vs_mceta");
2058     aFolderObj->Add(h);
2059     
2060     h2F = (TH2F*)fRecMCTrackHist1->Projection(2,1);
2061     h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2062     h->SetXTitle("#eta_{mc}");
2063     h->SetYTitle("(p_{t}-p_{tmc})/p_{tmc} mean");
2064     h->Draw();
2065     h->SetName("dpt_mean_vs_mceta");
2066     aFolderObj->Add(h);
2067     
2068     //
2069     h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2070     h = AlidNdPtHelper::MakeResol(h2F,1,0,kTRUE,10);
2071     h->SetXTitle("#eta_{mc}");
2072     h->SetYTitle("(#eta-#eta_{mc}) resolution");
2073     h->Draw();
2074     h->SetName("eta_resolution_vs_mceta");
2075     aFolderObj->Add(h);
2076     
2077     h2F = (TH2F*)fRecMCTrackHist1->Projection(3,1);
2078     h = AlidNdPtHelper::MakeResol(h2F,1,1,kTRUE,10);
2079     h->SetXTitle("#eta_{mc}");
2080     h->SetYTitle("(#eta-mc#eta) mean");
2081     h->Draw();
2082     h->SetName("deta_mean_vs_mceta");
2083     aFolderObj->Add(h);
2084     
2085     fRecMCTrackHist1->GetAxis(0)->SetRange(1,fRecMCTrackHist1->GetAxis(0)->GetNbins()); 
2086     
2087   } // end use MC info
2088   
2089   // export objects to analysis folder
2090   fAnalysisFolder = ExportToFolder(aFolderObj);
2091   
2092   // delete only TObjArray
2093   if(aFolderObj) delete aFolderObj;
2094 }
2095
2096 //_____________________________________________________________________________
2097 TFolder* AlidNdPtAnalysisPbPb2011::ExportToFolder(TObjArray * const array) 
2098 {
2099   // recreate folder avery time and export objects to new one
2100   //
2101   AlidNdPtAnalysisPbPb2011 * comp=this;
2102   TFolder *folder = comp->GetAnalysisFolder();
2103   
2104   TString name, title;
2105   TFolder *newFolder = 0;
2106   Int_t i = 0;
2107   Int_t size = array->GetSize();
2108   
2109   if(folder) { 
2110     // get name and title from old folder
2111     name = folder->GetName();  
2112     title = folder->GetTitle();  
2113     
2114     // delete old one
2115     delete folder;
2116     
2117     // create new one
2118     newFolder = CreateFolder(name.Data(),title.Data());
2119     newFolder->SetOwner();
2120     
2121     // add objects to folder
2122     while(i < size) {
2123       newFolder->Add(array->At(i));
2124       i++;
2125     }
2126   }
2127   
2128   return newFolder;
2129 }
2130
2131 //_____________________________________________________________________________
2132 TFolder* AlidNdPtAnalysisPbPb2011::CreateFolder(TString name,TString title) { 
2133   // create folder for analysed histograms
2134   //
2135   TFolder *folder = 0;
2136   folder = new TFolder(name.Data(),title.Data());
2137   
2138   return folder;
2139 }