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