743dad890bca38a9557a38c769ce96e7a00c7029
[u/mrichter/AliRoot.git] / PWG0 / dNdPt / AlidNdPtCorrection.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 #include <iostream>\r
17 \r
18 #include "TFile.h"\r
19 #include "TCint.h"\r
20 #include "TH1.h"\r
21 #include "TH2.h"\r
22 \r
23 #include "AliHeader.h"  \r
24 #include "AliGenEventHeader.h"  \r
25 #include "AliStack.h"  \r
26 #include "AliESDEvent.h"  \r
27 #include "AliMCEvent.h"  \r
28 #include "AliESDtrackCuts.h"  \r
29 #include "AliLog.h" \r
30 \r
31 #include "AlidNdPtEventCuts.h"\r
32 #include "AlidNdPtAcceptanceCuts.h"\r
33 \r
34 #include "AliPWG0Helper.h"\r
35 #include "AlidNdPtHelper.h"\r
36 #include "AlidNdPtCorrection.h"\r
37 \r
38 using namespace std;\r
39 \r
40 ClassImp(AlidNdPtCorrection)\r
41 \r
42 //_____________________________________________________________________________\r
43 //AlidNdPtCorrection::AlidNdPtCorrection(): TNamed(),\r
44   AlidNdPtCorrection::AlidNdPtCorrection(): AlidNdPt(),\r
45   fCorrectionFolder(0),\r
46   fMCEventHist1(0),\r
47   fRecEventHist1(0),\r
48   fMCAllEventMultHist1(0),\r
49   fMCAllNDEventMultHist1(0),\r
50   fMCAllNSDEventMultHist1(0),\r
51   fMCTriggerMultHist1(0),\r
52   fMCEventMultHist1(0),\r
53   fMCAllPrimTrackMultHist1(0),\r
54   fMCNDEventAllPrimTrackMultHist1(0),\r
55   fMCNSDEventAllPrimTrackMultHist1(0),\r
56   fMCTriggerPrimTrackMultHist1(0),\r
57   fMCEventPrimTrackMultHist1(0),\r
58   fEventMultCorrelationMatrix(0),\r
59   fZvNorm(0),\r
60   fZvEmptyEventsNorm(0),\r
61   fCorrTriggerMBtoInelEventMatrix(0),\r
62   fCorrTriggerMBtoNDEventMatrix(0),\r
63   fCorrTriggerMBtoNSDEventMatrix(0),\r
64   fCorrEventMatrix(0),\r
65   fCorrTriggerMBtoInelTrackEventMatrix(0),\r
66   fCorrTriggerMBtoNDTrackEventMatrix(0),\r
67   fCorrTriggerMBtoNSDTrackEventMatrix(0),\r
68   fCorrTrackEventMatrix(0),\r
69   fCorrTrackMatrix(0),\r
70   fContTrackMatrix(0),\r
71   fContMultTrackMatrix(0),\r
72   fCorrMatrixFileName("")\r
73 {\r
74   // default constructor\r
75   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { \r
76     fRecTrackHist1[i]=0;     \r
77   }\r
78 \r
79   for(Int_t i=0; i<8; i++) { \r
80     fCorrRecTrackMultHist1[i] = 0;\r
81   }\r
82 \r
83   for(Int_t i=0; i<5; i++) { \r
84     fCorrRecEventHist1[i] = 0;\r
85     fCorrRecEventHist2[i] = 0;\r
86   }\r
87 \r
88   Init();\r
89 }\r
90 \r
91 //_____________________________________________________________________________\r
92 AlidNdPtCorrection::AlidNdPtCorrection(Char_t* name, Char_t* title, TString corrMatrixFileName): AlidNdPt(name,title),\r
93   fCorrectionFolder(0),\r
94   fMCEventHist1(0),\r
95   fRecEventHist1(0),\r
96   fMCAllEventMultHist1(0),\r
97   fMCAllNDEventMultHist1(0),\r
98   fMCAllNSDEventMultHist1(0),\r
99   fMCTriggerMultHist1(0),\r
100   fMCEventMultHist1(0),\r
101   fMCAllPrimTrackMultHist1(0),\r
102   fMCNDEventAllPrimTrackMultHist1(0),\r
103   fMCNSDEventAllPrimTrackMultHist1(0),\r
104   fMCTriggerPrimTrackMultHist1(0),\r
105   fMCEventPrimTrackMultHist1(0),\r
106   fEventMultCorrelationMatrix(0),\r
107   fZvNorm(0),\r
108   fZvEmptyEventsNorm(0),\r
109   fCorrTriggerMBtoInelEventMatrix(0),\r
110   fCorrTriggerMBtoNDEventMatrix(0),\r
111   fCorrTriggerMBtoNSDEventMatrix(0),\r
112   fCorrEventMatrix(0),\r
113   fCorrTriggerMBtoInelTrackEventMatrix(0),\r
114   fCorrTriggerMBtoNDTrackEventMatrix(0),\r
115   fCorrTriggerMBtoNSDTrackEventMatrix(0),\r
116   fCorrTrackEventMatrix(0),\r
117   fCorrTrackMatrix(0),\r
118   fContTrackMatrix(0),\r
119   fContMultTrackMatrix(0),\r
120   fCorrMatrixFileName(corrMatrixFileName)\r
121 {\r
122   // constructor\r
123   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { \r
124     fRecTrackHist1[i]=0;     \r
125   }\r
126 \r
127   for(Int_t i=0; i<8; i++) { \r
128     fCorrRecTrackMultHist1[i] = 0;\r
129     fPtvsPt[i] = 0;\r
130   }\r
131 \r
132   for(Int_t i=0; i<5; i++) { \r
133     fCorrRecEventHist1[i] = 0;\r
134     fCorrRecEventHist2[i] = 0;\r
135   }\r
136 \r
137   Init();\r
138 }\r
139 \r
140 //_____________________________________________________________________________\r
141 AlidNdPtCorrection::~AlidNdPtCorrection() {\r
142   // \r
143   // destructor\r
144   //\r
145   if(fMCEventHist1) delete fMCEventHist1; fMCEventHist1=0;\r
146   if(fRecEventHist1) delete fRecEventHist1; fRecEventHist1=0;\r
147 \r
148   if(fMCAllEventMultHist1) delete fMCAllEventMultHist1; fMCAllEventMultHist1=0;\r
149   if(fMCAllNDEventMultHist1) delete fMCAllNDEventMultHist1; fMCAllNDEventMultHist1=0;\r
150   if(fMCAllNSDEventMultHist1) delete fMCAllNSDEventMultHist1; fMCAllNSDEventMultHist1=0;\r
151   if(fMCTriggerMultHist1) delete fMCTriggerMultHist1; fMCTriggerMultHist1=0;\r
152   if(fMCEventMultHist1) delete fMCEventMultHist1; fMCEventMultHist1=0;\r
153 \r
154   if(fMCAllPrimTrackMultHist1) delete fMCAllPrimTrackMultHist1; fMCAllPrimTrackMultHist1=0;\r
155   if(fMCNDEventAllPrimTrackMultHist1) delete fMCNDEventAllPrimTrackMultHist1; fMCNDEventAllPrimTrackMultHist1=0;\r
156   if(fMCNSDEventAllPrimTrackMultHist1) delete fMCNSDEventAllPrimTrackMultHist1; fMCNSDEventAllPrimTrackMultHist1=0;\r
157   if(fMCTriggerPrimTrackMultHist1) delete fMCTriggerPrimTrackMultHist1; fMCTriggerPrimTrackMultHist1=0;\r
158   if(fMCEventPrimTrackMultHist1) delete fMCEventPrimTrackMultHist1; fMCEventPrimTrackMultHist1=0;\r
159 \r
160   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) { \r
161     if(fRecTrackHist1[i]) delete fRecTrackHist1[i]; fRecTrackHist1[i]=0;\r
162   }\r
163 \r
164   for(Int_t i=0; i<8; i++) { \r
165     if(fCorrRecTrackMultHist1[i]) delete fCorrRecTrackMultHist1[i]; fCorrRecTrackMultHist1[i]=0;\r
166     if(fPtvsPt[i]) delete fPtvsPt[i]; fPtvsPt[i]=0;\r
167   }\r
168 \r
169   for(Int_t i=0; i<5; i++) { \r
170     if(fCorrRecEventHist1[i]) delete fCorrRecEventHist1[i]; fCorrRecEventHist1[i]=0;\r
171     if(fCorrRecEventHist2[i]) delete fCorrRecEventHist2[i]; fCorrRecEventHist2[i]=0;\r
172   }\r
173 \r
174   if(fCorrectionFolder) delete fCorrectionFolder; fCorrectionFolder=0;\r
175 }\r
176 \r
177 //_____________________________________________________________________________\r
178 void AlidNdPtCorrection::Init(){\r
179   //\r
180   // Init histograms\r
181   //\r
182   const Int_t etaNbins = 30; \r
183   const Int_t zvNbins = 12;\r
184 \r
185   // UA1 bining\r
186   //const Int_t ptNbins = 52; \r
187   //Double_t binsPt[ptNbins+1] = { 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.10, 3.20, 3.30, 3.40, 3.50, 3.60, 3.70, 3.80, 3.90, 4.00, 4.20, 4.40, 4.60, 4.80, 5.00, 5.20, 5.40, 5.60, 5.80, 6.00, 7.00, 7.60, 8.80, 9.60 }; \r
188 \r
189   const Int_t ptNbins = 56; \r
190   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,7.5,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0};\r
191   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};\r
192   Double_t binsZv[zvNbins+1] = {-30.,-25.,-20.,-15.,-10.,-5.,0.,5.,10.,15.,20.,25.,30.};\r
193 \r
194 \r
195   //\r
196   Int_t binsMCEventHist1[3]={100,100,140};\r
197   Double_t minMCEventHist1[3]={-0.1,-0.1,-35.}; \r
198   Double_t maxMCEventHist1[3]={0.1,0.1,35.}; \r
199   fMCEventHist1 = new THnSparseF("fMCEventHist1","mcXv:mcYv:mcZv",3,binsMCEventHist1,minMCEventHist1,maxMCEventHist1);\r
200   fMCEventHist1->GetAxis(0)->SetTitle("mcXv (cm)");\r
201   fMCEventHist1->GetAxis(1)->SetTitle("mcYv (cm)");\r
202   fMCEventHist1->GetAxis(2)->SetTitle("mcZv (cm)");\r
203   fMCEventHist1->Sumw2();\r
204 \r
205   //\r
206   Int_t binsRecEventHist1[3]={100,100,140};\r
207   Double_t minRecEventHist1[3]={-3.,-3.,-35.}; \r
208   Double_t maxRecEventHist1[3]={3.,3.,35.}; \r
209   fRecEventHist1 = new THnSparseF("fRecEventHist1","Xv:Yv:Zv",3,binsRecEventHist1,minRecEventHist1,maxRecEventHist1);\r
210   fRecEventHist1->GetAxis(0)->SetTitle("Xv (cm)");\r
211   fRecEventHist1->GetAxis(1)->SetTitle("Yv (cm)");\r
212   fRecEventHist1->GetAxis(2)->SetTitle("Zv (cm)");\r
213   fRecEventHist1->Sumw2();\r
214 \r
215   //\r
216   char name[256];\r
217   char title[256];\r
218 \r
219   Int_t binsMCAllPrimTrackMultHist1[3]={ptNbins,etaNbins,150};\r
220   Double_t minMCAllPrimTrackMultHist1[3]={0.,-1.,-0.5}; \r
221   Double_t maxMCAllPrimTrackMultHist1[3]={20.,1.,149.5}; \r
222   sprintf(name,"fMCAllPrimTrackMultHist1");\r
223   sprintf(title,"mcPt:mcEta:multiplicity");\r
224   \r
225   fMCAllPrimTrackMultHist1 = new THnSparseF(name,title,3,binsMCAllPrimTrackMultHist1,minMCAllPrimTrackMultHist1,maxMCAllPrimTrackMultHist1);\r
226   fMCAllPrimTrackMultHist1->SetBinEdges(0,binsPt);\r
227   fMCAllPrimTrackMultHist1->SetBinEdges(1,binsEta);\r
228   fMCAllPrimTrackMultHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");\r
229   fMCAllPrimTrackMultHist1->GetAxis(1)->SetTitle("mcEta");\r
230   fMCAllPrimTrackMultHist1->GetAxis(2)->SetTitle("multiplicity");\r
231   fMCAllPrimTrackMultHist1->Sumw2();\r
232 \r
233   Int_t binsMCNDEventAllPrimTrackMultHist1[3]={ptNbins,etaNbins,150};\r
234   Double_t minMCNDEventAllPrimTrackMultHist1[3]={0.,-1.,-0.5}; \r
235   Double_t maxMCNDEventAllPrimTrackMultHist1[3]={20.,1.,149.5}; \r
236   sprintf(name,"fMCNDEventAllPrimTrackMultHist1");\r
237   sprintf(title,"mcPt:mcEta:multiplicity");\r
238   \r
239   fMCNDEventAllPrimTrackMultHist1 = new THnSparseF(name,title,3,binsMCNDEventAllPrimTrackMultHist1,minMCNDEventAllPrimTrackMultHist1,maxMCNDEventAllPrimTrackMultHist1);\r
240   fMCNDEventAllPrimTrackMultHist1->SetBinEdges(0,binsPt);\r
241   fMCNDEventAllPrimTrackMultHist1->SetBinEdges(1,binsEta);\r
242   fMCNDEventAllPrimTrackMultHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");\r
243   fMCNDEventAllPrimTrackMultHist1->GetAxis(1)->SetTitle("mcEta");\r
244   fMCNDEventAllPrimTrackMultHist1->GetAxis(2)->SetTitle("multiplicity");\r
245   fMCNDEventAllPrimTrackMultHist1->Sumw2();\r
246 \r
247   Int_t binsMCNSDEventAllPrimTrackMultHist1[3]={ptNbins,etaNbins,150};\r
248   Double_t minMCNSDEventAllPrimTrackMultHist1[3]={0.,-1.,-0.5}; \r
249   Double_t maxMCNSDEventAllPrimTrackMultHist1[3]={20.,1.,140.9}; \r
250   sprintf(name,"fMCNSDEventAllPrimTrackMultHist1");\r
251   sprintf(title,"mcPt:mcEta:multiplicity");\r
252   \r
253   fMCNSDEventAllPrimTrackMultHist1 = new THnSparseF(name,title,3,binsMCNSDEventAllPrimTrackMultHist1,minMCNSDEventAllPrimTrackMultHist1,maxMCNSDEventAllPrimTrackMultHist1);\r
254   fMCNSDEventAllPrimTrackMultHist1->SetBinEdges(0,binsPt);\r
255   fMCNSDEventAllPrimTrackMultHist1->SetBinEdges(1,binsEta);\r
256   fMCNSDEventAllPrimTrackMultHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");\r
257   fMCNSDEventAllPrimTrackMultHist1->GetAxis(1)->SetTitle("mcEta");\r
258   fMCNSDEventAllPrimTrackMultHist1->GetAxis(2)->SetTitle("multiplicity");\r
259   fMCNSDEventAllPrimTrackMultHist1->Sumw2();\r
260 \r
261   Int_t binsMCEventTriggerPrimTrackMultHist1[3]={ptNbins,etaNbins,150};\r
262   Double_t minMCEventTriggerPrimTrackMultHist1[3]={0.,-1.,-0.5}; \r
263   Double_t maxMCEventTriggerPrimTrackMultHist1[3]={20.,1.,149.5}; \r
264   sprintf(name,"fMCTriggerPrimTrackMultHist1");\r
265   sprintf(title,"mcPt:mcEta:multiplicity");\r
266   \r
267   fMCTriggerPrimTrackMultHist1 = new THnSparseF(name,title,3,binsMCEventTriggerPrimTrackMultHist1,minMCEventTriggerPrimTrackMultHist1,maxMCEventTriggerPrimTrackMultHist1);\r
268   fMCTriggerPrimTrackMultHist1->SetBinEdges(0,binsPt);\r
269   fMCTriggerPrimTrackMultHist1->SetBinEdges(1,binsEta);\r
270   fMCTriggerPrimTrackMultHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");\r
271   fMCTriggerPrimTrackMultHist1->GetAxis(1)->SetTitle("mcEta");\r
272   fMCTriggerPrimTrackMultHist1->GetAxis(2)->SetTitle("multiplicity");\r
273   fMCTriggerPrimTrackMultHist1->Sumw2();\r
274 \r
275   Int_t binsMCEventPrimTrackMultHist1[3]={ptNbins,etaNbins,150};\r
276   Double_t minMCEventPrimTrackMultHist1[3]={0.,-1.,-0.5}; \r
277   Double_t maxMCEventPrimTrackMultHist1[3]={20.,1.,149.5}; \r
278   sprintf(name,"fMCEventPrimTrackMultHist1");\r
279   sprintf(title,"mcPt:mcEta:multiplicity");\r
280   \r
281   fMCEventPrimTrackMultHist1 = new THnSparseF(name,title,3,binsMCEventPrimTrackMultHist1,minMCEventPrimTrackMultHist1,maxMCEventPrimTrackMultHist1);\r
282   fMCEventPrimTrackMultHist1->SetBinEdges(0,binsPt);\r
283   fMCEventPrimTrackMultHist1->SetBinEdges(1,binsEta);\r
284   fMCEventPrimTrackMultHist1->GetAxis(0)->SetTitle("mcPt (GeV/c)");\r
285   fMCEventPrimTrackMultHist1->GetAxis(1)->SetTitle("mcEta");\r
286   fMCEventPrimTrackMultHist1->GetAxis(2)->SetTitle("multiplicity");\r
287   fMCEventPrimTrackMultHist1->Sumw2();\r
288 \r
289   for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) \r
290   {\r
291     // THnSparse track histograms\r
292     //\r
293     Int_t binsRecTrackHist1[3]={ptNbins,etaNbins,90};\r
294     Double_t minRecTrackHist1[3]={0.,-1.,0.}; \r
295     Double_t maxRecTrackHist1[3]={10.,1.,2.*TMath::Pi()};\r
296     sprintf(name,"fRecTrackHist1_%d",i);\r
297     sprintf(title,"Pt:Eta:Phi");\r
298   \r
299     fRecTrackHist1[i] = new THnSparseF(name,title,3,binsRecTrackHist1,minRecTrackHist1,maxRecTrackHist1);\r
300     fRecTrackHist1[i]->SetBinEdges(0,binsPt);\r
301     fRecTrackHist1[i]->SetBinEdges(1,binsEta);\r
302     fRecTrackHist1[i]->GetAxis(0)->SetTitle("Pt (GeV/c)");\r
303     fRecTrackHist1[i]->GetAxis(1)->SetTitle("Eta");\r
304     fRecTrackHist1[i]->GetAxis(2)->SetTitle("Phi (rad)");\r
305     fRecTrackHist1[i]->Sumw2();\r
306   }\r
307 \r
308   Int_t binsCorrRecTrackMultHist1[3]={ptNbins,etaNbins,150};\r
309   Double_t minCorrRecTrackMultHist1[3]={0.,-1.,-0.5}; \r
310   Double_t maxCorrRecTrackMultHist1[3]={20.,1.,149.5};\r
311 \r
312   Int_t binsPtvsPt[3]={ptNbins,320};\r
313   Double_t minPtvsPt[3]={0.,0.}; \r
314   Double_t maxPtvsPt[3]={20.,16.};\r
315 \r
316   for(Int_t i=0; i<8; i++) \r
317   {\r
318     // THnSparse track histograms\r
319     sprintf(name,"fCorrRecTrackMultHist1_%d",i);\r
320     sprintf(title,"Pt:Eta:mult");\r
321     fCorrRecTrackMultHist1[i] = new THnSparseF(name,title,3,binsCorrRecTrackMultHist1,minCorrRecTrackMultHist1,maxCorrRecTrackMultHist1);\r
322     fCorrRecTrackMultHist1[i]->SetBinEdges(0,binsPt);\r
323     fCorrRecTrackMultHist1[i]->SetBinEdges(1,binsEta);\r
324     fCorrRecTrackMultHist1[i]->GetAxis(0)->SetTitle("Pt (GeV/c)");\r
325     fCorrRecTrackMultHist1[i]->GetAxis(1)->SetTitle("Eta");\r
326     fCorrRecTrackMultHist1[i]->GetAxis(2)->SetTitle("multiplicity");\r
327     fCorrRecTrackMultHist1[i]->Sumw2();\r
328 \r
329     sprintf(name,"fPtvsPt_%d",i);\r
330     sprintf(title,"Pt:Pt");\r
331     fPtvsPt[i] = new THnSparseF(name,title,2,binsPtvsPt,minPtvsPt,maxPtvsPt);\r
332     fPtvsPt[i]->SetBinEdges(0,binsPt);\r
333     fPtvsPt[i]->GetAxis(0)->SetTitle("Pt (GeV/c)"); \r
334     fPtvsPt[i]->GetAxis(1)->SetTitle("Pt (GeV/c)"); \r
335     fPtvsPt[i]->Sumw2();\r
336   }\r
337 \r
338   Int_t binsEventMatrix[2]={zvNbins,150};\r
339   Double_t minEventMatrix[2]={-25.,-0.5};\r
340   Double_t maxEventMatrix[2]={25.,149.5};\r
341 \r
342   fMCAllEventMultHist1 = new THnSparseF("fMCAllEventMultHist1","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
343   fMCAllEventMultHist1->SetBinEdges(0,binsZv);\r
344   fMCAllEventMultHist1->GetAxis(0)->SetTitle("mcZv (cm)");\r
345   fMCAllEventMultHist1->GetAxis(1)->SetTitle("multiplicity");\r
346   fMCAllEventMultHist1->Sumw2();\r
347 \r
348   fMCAllNDEventMultHist1 = new THnSparseF("fMCAllNDEventMultHist1","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
349   fMCAllNDEventMultHist1->SetBinEdges(0,binsZv);\r
350   fMCAllNDEventMultHist1->GetAxis(0)->SetTitle("mcZv (cm)");\r
351   fMCAllNDEventMultHist1->GetAxis(1)->SetTitle("multiplicity");\r
352   fMCAllNDEventMultHist1->Sumw2();\r
353 \r
354   fMCAllNSDEventMultHist1 = new THnSparseF("fMCAllNSDEventMultHist1","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
355   fMCAllNSDEventMultHist1->SetBinEdges(0,binsZv);\r
356   fMCAllNSDEventMultHist1->GetAxis(0)->SetTitle("mcZv (cm)");\r
357   fMCAllNSDEventMultHist1->GetAxis(1)->SetTitle("multiplicity");\r
358   fMCAllNSDEventMultHist1->Sumw2();\r
359 \r
360   fMCTriggerMultHist1 = new THnSparseF("fMCTriggerMultHist1","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
361   fMCTriggerMultHist1->SetBinEdges(0,binsZv);\r
362   fMCTriggerMultHist1->GetAxis(0)->SetTitle("mcZv (cm)");\r
363   fMCTriggerMultHist1->GetAxis(1)->SetTitle("multiplicity");\r
364   fMCTriggerMultHist1->Sumw2();\r
365 \r
366   fMCEventMultHist1 = new THnSparseF("fMCEventMultHist1","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
367   fMCEventMultHist1->SetBinEdges(0,binsZv);\r
368   fMCEventMultHist1->GetAxis(0)->SetTitle("mcZv (cm)");\r
369   fMCEventMultHist1->GetAxis(1)->SetTitle("multiplicity");\r
370   fMCEventMultHist1->Sumw2();\r
371 \r
372   for(Int_t i=0; i<5; i++) \r
373   {\r
374     // event corrected histograms\r
375     sprintf(name,"fCorrRecEventHist1_%d",i);\r
376     sprintf(title,"mcZv:mult");\r
377     fCorrRecEventHist1[i] = new THnSparseF("fCorrRecEventHist1","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
378     fCorrRecEventHist1[i]->SetBinEdges(0,binsZv);\r
379     fCorrRecEventHist1[i]->GetAxis(0)->SetTitle("Zv (cm)");\r
380     fCorrRecEventHist1[i]->GetAxis(1)->SetTitle("multiplicity");\r
381     fCorrRecEventHist1[i]->Sumw2();\r
382 \r
383     // empty event corrected histograms\r
384     sprintf(name,"fCorrRecEventHist2_%d",i);\r
385     sprintf(title,"mcZv:mult");\r
386     fCorrRecEventHist2[i] = new THnSparseF("fCorrRecEventHist2","mcZv:mult",2,binsEventMatrix,minEventMatrix,maxEventMatrix);\r
387     fCorrRecEventHist2[i]->SetBinEdges(0,binsZv);\r
388     fCorrRecEventHist2[i]->GetAxis(0)->SetTitle("Zv (cm)");\r
389     fCorrRecEventHist2[i]->GetAxis(1)->SetTitle("multiplicity");\r
390     fCorrRecEventHist2[i]->Sumw2();\r
391   }\r
392 \r
393   // init output folder\r
394   fCorrectionFolder = CreateFolder("folderdNdPt","Correction dNdPt Folder");\r
395 \r
396   // init correction matrices\r
397   TFile *file = (TFile *)TFile::Open(fCorrMatrixFileName.Data()); \r
398   if(!file) { \r
399     AliDebug(AliLog::kError, "file with efficiency matrices not available");\r
400     printf("file with efficiency matrices not available \n");\r
401   } else {\r
402     TFolder *folder = (TFolder *)file->FindObjectAny("folderdNdPt");\r
403     if(!folder) { \r
404       AliDebug(AliLog::kError, "file without folderdNdPt");\r
405       printf("file without folderdNdPt \n");\r
406     } else {\r
407       // rec. event mult vs true multiplicity \r
408       fEventMultCorrelationMatrix = (THnSparseF*)folder->FindObject("event_mult_correlation_matrix");\r
409       if(!fEventMultCorrelationMatrix) {\r
410          Printf("No %s matrix \n", "event_mult_correlation_matrix");\r
411       }\r
412 \r
413       //\r
414       // event level corrections (zv,mult_MB)\r
415       //\r
416  \r
417       // trigger bias correction (MBtoND) \r
418       fCorrTriggerMBtoNDEventMatrix = (THnSparseF*)folder->FindObject("zv_mult_trig_MBtoND_corr_matrix");\r
419       if(!fCorrTriggerMBtoNDEventMatrix) {\r
420          Printf("No %s matrix \n", "zv_mult_trig_MBtoND_corr_matrix");\r
421       }\r
422 \r
423       // trigger bias correction (MBtoNSD)\r
424       fCorrTriggerMBtoNSDEventMatrix = (THnSparseF*)folder->FindObject("zv_mult_trig_MBtoNSD_corr_matrix");\r
425       if(!fCorrTriggerMBtoNSDEventMatrix) {\r
426          Printf("No %s matrix \n", "zv_mult_trig_MBtoNSD_corr_matrix");\r
427       }\r
428 \r
429       // trigger bias correction (MBtoInel)\r
430       fCorrTriggerMBtoInelEventMatrix = (THnSparseF*)folder->FindObject("zv_mult_trig_MBtoInel_corr_matrix");\r
431       if(!fCorrTriggerMBtoInelEventMatrix) {\r
432          Printf("No %s matrix \n", "zv_mult_trig_MBtoInel_corr_matrix"); \r
433       }\r
434      \r
435       // vertex reconstruction efficiency correction\r
436       fCorrEventMatrix = (THnSparseF*)folder->FindObject("zv_mult_event_corr_matrix");\r
437       if(!fCorrEventMatrix) {\r
438          Printf("No %s matrix \n", "zv_mult_event_corr_matrix");\r
439       }\r
440 \r
441       //\r
442       // histogram needed for empty events corrections\r
443       //\r
444       fZvNorm = (TH1D*)folder->FindObject("zv_distribution_norm");\r
445       if(!fZvNorm) {\r
446          Printf("No %s matrix \n", "fZvNorm");\r
447       }\r
448 \r
449       fZvEmptyEventsNorm = (TH1D*)folder->FindObject("zv_empty_events_norm");\r
450       if(!fZvEmptyEventsNorm) {\r
451          Printf("No %s matrix \n", "fZvEmptyEventsNorm");\r
452       }\r
453 \r
454       //\r
455       // track-event level corrections (zv,pt,eta)\r
456       //\r
457 \r
458       // trigger bias correction (MBtoND) \r
459       fCorrTriggerMBtoNDTrackEventMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_track_trig_MBtoND_corr_matrix");\r
460       if(!fCorrTriggerMBtoNDTrackEventMatrix) {\r
461          Printf("No %s matrix \n", "zv_pt_eta_track_trig_MBtoND_corr_matrix");\r
462       }\r
463 \r
464       // trigger bias correction (MBtoNSD)\r
465       fCorrTriggerMBtoNSDTrackEventMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_track_trig_MBtoNSD_corr_matrix");\r
466       if(!fCorrTriggerMBtoNSDTrackEventMatrix) {\r
467          Printf("No %s matrix \n", "zv_pt_eta_track_trig_MBtoNSD_corr_matrix");\r
468       }\r
469 \r
470       // trigger bias correction (MBtoInel) \r
471       fCorrTriggerMBtoInelTrackEventMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_track_trig_MBtoInel_corr_matrix");\r
472       if(!fCorrTriggerMBtoInelTrackEventMatrix) {\r
473          Printf("No %s matrix \n", "zv_pt_eta_track_trig_MBtoInel_corr_matrix");\r
474       }\r
475     \r
476       // vertex reconstruction efficiency correction (zv,pt,eta)\r
477       fCorrTrackEventMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_track_event_corr_matrix");\r
478       if(!fCorrTrackEventMatrix) {\r
479          Printf("No %s matrix \n", "zv_pt_eta_track_event_corr_matrix");\r
480       }\r
481 \r
482       // track reconstruction efficiency correction (zv,pt,eta)\r
483       fCorrTrackMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_track_corr_matrix");\r
484       if(!fCorrTrackMatrix) {\r
485          Printf("No %s matrix \n", "zv_pt_eta_track_corr_matrix");\r
486       }\r
487 \r
488       // secondary tracks contamination correction (zv,pt,eta)\r
489       fContTrackMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_track_cont_matrix");\r
490       if(!fContTrackMatrix) {\r
491          Printf("No %s matrix \n", "zv_pt_eta_track_cont_matrix");\r
492       }\r
493 \r
494       // multiply reconstructed tracks correction\r
495       fContMultTrackMatrix = (THnSparseF*)folder->FindObject("zv_pt_eta_mult_track_cont_matrix");\r
496       if(!fContMultTrackMatrix) {\r
497          Printf("No %s matrix \n", "zv_pt_eta_mult_track_cont_matrix");\r
498       }\r
499     }\r
500   }\r
501 \r
502 }\r
503 \r
504 //_____________________________________________________________________________\r
505 void AlidNdPtCorrection::Process(AliESDEvent *esdEvent, AliMCEvent *mcEvent)\r
506 {\r
507   //\r
508   // Process real and/or simulated events\r
509   //\r
510   if(!esdEvent) {\r
511     AliDebug(AliLog::kError, "esdEvent not available");\r
512     return;\r
513   }\r
514 \r
515   // get selection cuts\r
516   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
517   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
518   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
519 \r
520   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
521     AliDebug(AliLog::kError, "cuts not available");\r
522     return;\r
523   }\r
524 \r
525   // trigger selection\r
526   Bool_t isEventTriggered = kTRUE;\r
527   if(evtCuts->IsTriggerRequired())  {\r
528     if(IsUseMCInfo()) { \r
529       static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;\r
530       isEventTriggered = triggerAnalysis->IsTriggerFired(esdEvent, GetTrigger());\r
531     }\r
532     else {\r
533       isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());\r
534     }\r
535   }\r
536 \r
537   // use MC information\r
538   AliHeader* header = 0;\r
539   AliGenEventHeader* genHeader = 0;\r
540   AliStack* stack = 0;\r
541   TArrayF vtxMC(3);\r
542   AliPWG0Helper::MCProcessType evtType = AliPWG0Helper::kInvalidProcess;\r
543   Int_t multMCTrueTracks = 0;\r
544 \r
545   if(IsUseMCInfo())\r
546   {\r
547     if(!mcEvent) {\r
548       AliDebug(AliLog::kError, "mcEvent not available");\r
549       return;\r
550     }\r
551 \r
552     // get MC event header\r
553     header = mcEvent->Header();\r
554     if (!header) {\r
555       AliDebug(AliLog::kError, "Header not available");\r
556       return;\r
557     }\r
558     // MC particle stack\r
559     stack = mcEvent->Stack();\r
560     if (!stack) {\r
561       AliDebug(AliLog::kError, "Stack not available");\r
562       return;\r
563     }\r
564 \r
565     // get event type (ND=0x1, DD=0x2, SD=0x4)\r
566     evtType = AliPWG0Helper::GetEventProcessType(header);\r
567     //Printf("evtType %d \n", evtType);\r
568     AliDebug(AliLog::kDebug+1, Form("Found process type %d", evtType));\r
569 \r
570     // get MC vertex\r
571     genHeader = header->GenEventHeader();\r
572     if (!genHeader) {\r
573       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
574       return;\r
575     }\r
576     genHeader->PrimaryVertex(vtxMC);\r
577 \r
578     // Fill MC event histogram\r
579     Double_t vMCEventHist1[3]={vtxMC[0],vtxMC[1],vtxMC[2]};\r
580     fMCEventHist1->Fill(vMCEventHist1);\r
581 \r
582     // multipliticy of all MC primary tracks\r
583     // in Zvtx, pt and eta ranges\r
584     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
585 \r
586   } // end bUseMC\r
587 \r
588   // get reconstructed vertex  \r
589   const AliESDVertex* vtxESD = 0; \r
590   Bool_t isRecVertex = kFALSE;\r
591   if(evtCuts->IsRecVertexRequired()) \r
592   {\r
593     Bool_t bRedoTPCVertex = evtCuts->IsRedoTPCVertex();\r
594     Bool_t bUseConstraints = evtCuts->IsUseBeamSpotConstraint();\r
595     vtxESD = AlidNdPtHelper::GetVertex(esdEvent,evtCuts,accCuts,esdTrackCuts,GetAnalysisMode(),kFALSE,bRedoTPCVertex,bUseConstraints); \r
596     isRecVertex = AlidNdPtHelper::TestRecVertex(vtxESD, GetAnalysisMode(), kFALSE); // should be moved to AcceptEvent\r
597   }\r
598   if( IsUseMCInfo() &&  !evtCuts->IsRecVertexRequired() ) {\r
599     vtxESD = new AliESDVertex(vtxMC[2],10.,genHeader->NProduced());\r
600     isRecVertex = kTRUE;\r
601   }\r
602   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD) && isRecVertex; \r
603   //printf("isEventOK %d \n",isEventOK);\r
604 \r
605   //\r
606   // get multiplicity of min. bias tracks\r
607   //\r
608   Int_t multMBTracks = 0; \r
609   if(GetAnalysisMode() == AlidNdPtHelper::kTPC || GetAnalysisMode() == AlidNdPtHelper::kMCRec) {  \r
610     multMBTracks = AlidNdPtHelper::GetTPCMBTrackMult(esdEvent,evtCuts,accCuts,esdTrackCuts);\r
611   } \r
612   else if(GetAnalysisMode() == AlidNdPtHelper::kTPCSPDvtx || GetAnalysisMode() == AlidNdPtHelper::kMCRec) {\r
613     multMBTracks = AlidNdPtHelper::GetSPDMBTrackMult(esdEvent,0.0);\r
614   } \r
615   else {\r
616     AliDebug(AliLog::kError, Form("Found analysis type %s", GetAnalysisMode()));\r
617     return; \r
618   }\r
619 \r
620   //\r
621   // correct event and track histograms\r
622   //\r
623   TObjArray *allChargedTracks=0;\r
624   Int_t multRec=0, multRecTemp=0;\r
625   Int_t *labelsRec=0;\r
626 \r
627   if(isEventOK && isEventTriggered)\r
628   {\r
629     // get all charged tracks\r
630     //allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,vtxESD,GetAnalysisMode());\r
631     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());\r
632     if(!allChargedTracks) return;\r
633 \r
634     Int_t entries = allChargedTracks->GetEntries();\r
635     labelsRec = new Int_t[entries];\r
636 \r
637     // calculate mult of reconstructed tracks\r
638     for(Int_t i=0; i<entries;++i) \r
639     {\r
640       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);\r
641       if(!track) continue;\r
642       if(accCuts->AcceptTrack(track)) \r
643       {\r
644         if(esdTrackCuts->AcceptTrack(track)) \r
645         {\r
646           multRecTemp++;\r
647         }\r
648       }\r
649     }  \r
650 \r
651     for(Int_t i=0; i<entries;++i) \r
652     {\r
653       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);\r
654       if(!track) continue;\r
655         \r
656       if(accCuts->AcceptTrack(track)) \r
657       {\r
658         if(esdTrackCuts->AcceptTrack(track)) \r
659         {\r
660           // track-level corrections\r
661           if(GetAnalysisMode()==AlidNdPtHelper::kMCRec && IsUseMCInfo()) { \r
662             FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,vtxMC[2],multRecTemp); \r
663           } else {\r
664             FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,vtxESD->GetZv(),multRecTemp); \r
665           }\r
666           labelsRec[multRec] = TMath::Abs(track->GetLabel());\r
667           multRec++;\r
668         }\r
669       }\r
670     }\r
671     // event-level corrections\r
672     if(GetAnalysisMode()==AlidNdPtHelper::kMCRec && IsUseMCInfo()) { \r
673       FillHistograms(AlidNdPtHelper::kRecEvents,vtxMC[2],multMBTracks);\r
674     }\r
675     else {\r
676       FillHistograms(AlidNdPtHelper::kRecEvents,vtxESD->GetZv(),multMBTracks);\r
677     }\r
678 \r
679     // control event histograms\r
680     Double_t vRecEventHist1[3] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv()};\r
681     fRecEventHist1->Fill(vRecEventHist1);\r
682   } \r
683 \r
684   // empty events corrections\r
685   // no reconstructed zv\r
686   if(isEventTriggered && multMBTracks==0) {\r
687     if(GetAnalysisMode()==AlidNdPtHelper::kMCRec && IsUseMCInfo()) {\r
688       FillHistograms(AlidNdPtHelper::kTriggeredEvents,vtxMC[2],multMBTracks);\r
689     }\r
690     else {\r
691       Double_t zv = fZvNorm->GetRandom();\r
692       FillHistograms(AlidNdPtHelper::kTriggeredEvents,zv,multMBTracks);\r
693     }\r
694   }\r
695 \r
696   if(IsUseMCInfo())  \r
697   {\r
698     // select MC events \r
699     if(evtCuts->AcceptMCEvent(mcEvent))\r
700     {\r
701       //\r
702       // event histograms\r
703       //\r
704       Double_t vMCEventMatrix[2] = {vtxMC[2],multMBTracks};\r
705       fMCAllEventMultHist1->Fill(vMCEventMatrix);\r
706       if(evtType == AliPWG0Helper::kND) {\r
707         fMCAllNDEventMultHist1->Fill(vMCEventMatrix);\r
708       }\r
709       if(evtType != AliPWG0Helper::kSD) {\r
710         fMCAllNSDEventMultHist1->Fill(vMCEventMatrix);\r
711       }\r
712       if(isEventTriggered) fMCTriggerMultHist1->Fill(vMCEventMatrix);\r
713       if(isEventTriggered && isEventOK) fMCEventMultHist1->Fill(vMCEventMatrix);\r
714 \r
715       //\r
716       // MC histograms for efficiency studies\r
717       //\r
718       Int_t nPart  = stack->GetNtrack();\r
719       for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
720       {\r
721         // print MC stack info\r
722         //AlidNdPtHelper::PrintMCInfo(stack,iMc);\r
723 \r
724         TParticle* particle = stack->Particle(iMc);\r
725         if (!particle)\r
726         continue;\r
727 \r
728         // only charged particles\r
729         Double_t charge = particle->GetPDG()->Charge()/3.;\r
730         if (charge == 0.0)\r
731           continue;\r
732       \r
733         // physical primary\r
734         Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
735 \r
736         // all primaries in acceptance\r
737         if(!accCuts->AcceptTrack(particle)) continue;\r
738         if(!prim) continue;\r
739 \r
740         Double_t gpt = particle->Pt();\r
741         Double_t geta = particle->Eta();\r
742 \r
743         Double_t valMCAllTrackMultHist1[3] = {gpt,geta,multRec};          \r
744         fMCAllPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
745         if(evtType == AliPWG0Helper::kND) {\r
746           fMCNDEventAllPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
747         }\r
748         if(evtType != AliPWG0Helper::kSD) {\r
749           fMCNSDEventAllPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
750         }\r
751         if(isEventTriggered) fMCTriggerPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
752         if(isEventTriggered && isEventOK) fMCEventPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
753       }\r
754     }\r
755   } // end bUseMC\r
756 \r
757   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;\r
758   if(labelsRec) delete [] labelsRec; labelsRec = 0;\r
759 \r
760   if(!evtCuts->IsRecVertexRequired() && vtxESD != NULL) delete vtxESD;\r
761 }\r
762 \r
763 //_____________________________________________________________________________\r
764 void AlidNdPtCorrection::FillHistograms(AlidNdPtHelper::EventObject eventObj, Double_t zv, Int_t multMBTracks)\r
765 {\r
766   //\r
767   // Fill corrected histograms\r
768   //\r
769 \r
770   Double_t vEventMatrix[2] = {zv,multMBTracks};\r
771   //\r
772   // Correct for efficiency \r
773   //\r
774   if(eventObj == AlidNdPtHelper::kRecEvents && multMBTracks>0)  \r
775   {\r
776     Double_t corrToMBF = GetCorrFactZvMult(fCorrEventMatrix,zv,multMBTracks);\r
777     Double_t corrToInelF = GetCorrFactZvMult(fCorrTriggerMBtoInelEventMatrix,zv,multMBTracks);\r
778     Double_t corrToNDF = GetCorrFactZvMult(fCorrTriggerMBtoNDEventMatrix,zv,multMBTracks);\r
779     Double_t corrToNSDF = GetCorrFactZvMult(fCorrTriggerMBtoNSDEventMatrix,zv,multMBTracks);\r
780     //printf("corrToMBF %f, corrToInelF %f, corrToNDF %f corrToNSDF %f \n",corrToMBF,corrToInelF,corrToNDF,corrToNSDF);\r
781 \r
782     fCorrRecEventHist1[0]->Fill(vEventMatrix);\r
783     fCorrRecEventHist1[1]->Fill(vEventMatrix,corrToMBF);\r
784     fCorrRecEventHist1[2]->Fill(vEventMatrix,corrToMBF*corrToInelF);\r
785     fCorrRecEventHist1[3]->Fill(vEventMatrix,corrToMBF*corrToNDF);\r
786     fCorrRecEventHist1[4]->Fill(vEventMatrix,corrToMBF*corrToNSDF);\r
787   }\r
788 \r
789   if(eventObj==AlidNdPtHelper::kTriggeredEvents && multMBTracks==0) // empty triggered events\r
790   {\r
791     Int_t bin = fZvEmptyEventsNorm->FindBin(zv);\r
792     Double_t Fz = fZvEmptyEventsNorm->GetBinContent(bin);\r
793     Double_t corrToInelF0 = GetCorrFactZvMult(fCorrTriggerMBtoInelEventMatrix,zv,multMBTracks);\r
794     Double_t corrToNDF0 = GetCorrFactZvMult(fCorrTriggerMBtoNDEventMatrix,zv,multMBTracks);\r
795     Double_t corrToNSDF0 = GetCorrFactZvMult(fCorrTriggerMBtoNSDEventMatrix,zv,multMBTracks);\r
796     //printf("Fz %f, corrToInelF0 %f, corrToNDF0 %f, corrToNSDF0 %f \n",Fz,corrToInelF0,corrToNDF0,corrToNSDF0);\r
797 \r
798     fCorrRecEventHist2[0]->Fill(vEventMatrix);\r
799     fCorrRecEventHist2[1]->Fill(vEventMatrix,Fz);\r
800     fCorrRecEventHist2[2]->Fill(vEventMatrix,Fz*corrToInelF0);\r
801     fCorrRecEventHist2[3]->Fill(vEventMatrix,Fz*corrToNDF0);\r
802     fCorrRecEventHist2[4]->Fill(vEventMatrix,Fz*corrToNSDF0);\r
803   }\r
804 }\r
805 \r
806 //_____________________________________________________________________________\r
807 void AlidNdPtCorrection::FillHistograms(AliESDtrack *esdTrack, AliStack *stack, AlidNdPtHelper::TrackObject trackObj, Double_t zv, Int_t mult)\r
808 {\r
809   //\r
810   // Fill ESD track and MC histograms \r
811   //\r
812   if(!esdTrack) return;\r
813 \r
814   //Float_t q = esdTrack->Charge();\r
815   Float_t pt = esdTrack->Pt();\r
816   Float_t eta = esdTrack->Eta();\r
817   Float_t phi = esdTrack->Phi();\r
818 \r
819   if(stack && GetAnalysisMode() == AlidNdPtHelper::kMCRec) \r
820   {\r
821     Int_t label = TMath::Abs(esdTrack->GetLabel());\r
822    \r
823     TParticle* particle = stack->Particle(label);\r
824     if(!particle) return;\r
825    \r
826     Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3\r
827     if(gq==0) return;\r
828     Float_t gpt = particle->Pt();\r
829     Float_t geta = particle->Eta();\r
830     Float_t gphi = particle->Phi();\r
831 \r
832     // replace reconstructed values with MC\r
833     pt = gpt;\r
834     eta = geta;\r
835     phi = gphi;\r
836   }\r
837 \r
838   //\r
839   // Fill histograms\r
840   //\r
841   Double_t values[3] = {pt,eta,phi};      \r
842   fRecTrackHist1[trackObj]->Fill(values);\r
843 \r
844   //\r
845   // Correct for contamination and efficiency \r
846   //\r
847   if(trackObj == AlidNdPtHelper::kRecTracks || GetAnalysisMode() == AlidNdPtHelper::kMCRec)  \r
848   {\r
849     // track level corrections\r
850     Double_t trackEffF = GetCorrFactZvPtEta(fCorrTrackMatrix,zv,pt,eta);\r
851     Double_t trackContF = GetContFactZvPtEta(fContTrackMatrix,zv,pt,eta);\r
852     Double_t multTrackContF = GetContFactZvPtEta(fContMultTrackMatrix,zv,pt,eta);\r
853     //printf("zv %f, pt %f, eta %f \n",zv,pt,eta);\r
854     //printf("trackEffF %f, trackContF %f, multTrackContF %f \n", trackEffF, trackContF, multTrackContF);\r
855    \r
856     // track-event level corrections\r
857     Double_t vertexEffF = GetCorrFactZvPtEta(fCorrTrackEventMatrix,zv,pt,eta);\r
858     Double_t trigMBToInel = GetCorrFactZvPtEta(fCorrTriggerMBtoInelTrackEventMatrix,zv,pt,eta);  \r
859     Double_t trigMBToND = GetCorrFactZvPtEta(fCorrTriggerMBtoNDTrackEventMatrix,zv,pt,eta);\r
860     Double_t trigMBToNSD = GetCorrFactZvPtEta(fCorrTriggerMBtoNSDTrackEventMatrix,zv,pt,eta);\r
861     //printf("vertexEffF %f, trigMBToInel %f, trigMBToNSD %f \n", vertexEffF, trigMBToInel, trigMBToNSD);\r
862     \r
863     Double_t corrF[8] = { 1.0, \r
864                           trackContF,\r
865                           trackContF*trackEffF,\r
866                           trackContF*trackEffF*multTrackContF,\r
867                           trackContF*trackEffF*multTrackContF*vertexEffF,\r
868                           trackContF*trackEffF*multTrackContF*vertexEffF*trigMBToInel,\r
869                           trackContF*trackEffF*multTrackContF*vertexEffF*trigMBToND,\r
870                           trackContF*trackEffF*multTrackContF*vertexEffF*trigMBToNSD\r
871                          }; \r
872  \r
873     // Fill histograms\r
874     Double_t valCorrRecTrackMultHist1[3] = {pt,eta,mult};         \r
875     Double_t valPtvsPt[2] = {pt,pt};      \r
876     for(Int_t i=0; i<8; i++) {\r
877       fCorrRecTrackMultHist1[i]->Fill(valCorrRecTrackMultHist1,corrF[i]);\r
878       fPtvsPt[i]->Fill(valPtvsPt,corrF[i]);\r
879     }\r
880   }\r
881 }\r
882 \r
883 //_____________________________________________________________________________\r
884 void AlidNdPtCorrection::FillHistograms(AliStack *stack, Int_t /*label*/, AlidNdPtHelper::TrackObject /*trackObj*/, Int_t /*mult*/)\r
885 {\r
886   // Fill MC histograms\r
887   if(!stack) return;\r
888 \r
889   /*\r
890   TParticle* particle = stack->Particle(label);\r
891   if(!particle) return;\r
892 \r
893   Int_t mother_pdg = -1;\r
894   TParticle* mother = 0;\r
895 \r
896   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);\r
897   Int_t motherLabel = particle->GetMother(0); \r
898   if(motherLabel>0) mother = stack->Particle(motherLabel);\r
899   if(mother) mother_pdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only\r
900   Int_t mech = particle->GetUniqueID(); // production mechanism\r
901 \r
902   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3 \r
903   Float_t gpt = particle->Pt();\r
904   Float_t qgpt = particle->Pt() * gq;\r
905   Float_t geta = particle->Eta();\r
906   Float_t gphi = particle->Phi();\r
907   Float_t gpz = particle->Pz();\r
908 \r
909   Bool_t prim = stack->IsPhysicalPrimary(label);\r
910   Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();\r
911 \r
912   Int_t pid=-1;\r
913   if (TMath::Abs(particle->GetPdgCode()) == kElectron)         { pid = 0; }\r
914     else if (TMath::Abs(particle->GetPdgCode()) == kMuonMinus) { pid = 1; }\r
915     else if (TMath::Abs(particle->GetPdgCode()) == kPiPlus)    { pid = 2; }\r
916     else if (TMath::Abs(particle->GetPdgCode()) == kKPlus)     { pid = 3; }\r
917     else if (TMath::Abs(particle->GetPdgCode()) == kProton)    { pid = 4; }\r
918     else                                                       { pid = 5; }\r
919     */\r
920 \r
921   //if(!prim) printf("prim_mother %d, mother %d, particle %d, production mech %d\n",prim_mother->GetPdgCode(),mother->GetPdgCode(), particle->GetPdgCode(),mech);\r
922   \r
923 }\r
924 \r
925 //_____________________________________________________________________________\r
926 Double_t AlidNdPtCorrection::GetCorrFactZvPtEta(THnSparse *hist, Double_t zv, Double_t pt, Double_t eta) const {\r
927 // return correction factor F(zv,pt,eta)\r
928 \r
929  if(!hist) return 1.;\r
930 \r
931  //\r
932  TAxis *ax = hist->GetAxis(0);\r
933  TAxis *ay = hist->GetAxis(1);\r
934  TAxis *az = hist->GetAxis(2);\r
935 \r
936  Int_t binx = ax->FindBin(zv);\r
937  Int_t biny = ay->FindBin(pt);\r
938  Int_t binz = az->FindBin(eta);\r
939  Int_t dim[3] = {binx,biny,binz};\r
940 \r
941  Double_t fact  = hist->GetBinContent(dim);  \r
942 \r
943 return fact;\r
944 }\r
945 \r
946 //_____________________________________________________________________________\r
947 Double_t AlidNdPtCorrection::GetContFactZvPtEta(THnSparse *hist, Double_t zv, Double_t pt, Double_t eta) const {\r
948 // return contamination correction factor F(zv,pt,eta)\r
949 \r
950  if(!hist) return 1.0;\r
951 \r
952  //\r
953  TAxis *ax = hist->GetAxis(0);\r
954  TAxis *ay = hist->GetAxis(1);\r
955  TAxis *az = hist->GetAxis(2);\r
956 \r
957  Int_t binx = ax->FindBin(zv);\r
958  Int_t biny = ay->FindBin(pt);\r
959  Int_t binz = az->FindBin(eta);\r
960  Int_t dim[3] = {binx,biny,binz};\r
961 \r
962  Double_t fact  = 1.0-hist->GetBinContent(dim);  \r
963  //Double_t fact  = hist->GetBinContent(dim);  \r
964 \r
965 return fact;\r
966 }\r
967 \r
968 //_____________________________________________________________________________\r
969 Double_t AlidNdPtCorrection::GetCorrFactZvMult(THnSparse *hist, Double_t zv, Int_t mult) const {\r
970 // return correction factor F(zv,mult)\r
971 \r
972  if(!hist) return 1.;\r
973 \r
974  TAxis *ax = hist->GetAxis(0);\r
975  TAxis *ay = hist->GetAxis(1);\r
976  Int_t binx = ax->FindBin(zv);\r
977  Int_t biny = ay->FindBin(mult);\r
978  Int_t dim[2] = {binx,biny};\r
979 \r
980  Double_t fact  = hist->GetBinContent(dim);  \r
981 \r
982 \r
983 return fact;\r
984 }\r
985 \r
986 //_____________________________________________________________________________\r
987 Double_t AlidNdPtCorrection::GetContFactZvMult(THnSparse *hist, Double_t zv, Int_t mult) const {\r
988 // return contamination correction factor F(zv,mult)\r
989 \r
990  if(!hist) return 1.;\r
991 \r
992  TAxis *ax = hist->GetAxis(0);\r
993  TAxis *ay = hist->GetAxis(1);\r
994  Int_t binx = ax->FindBin(zv);\r
995  Int_t biny = ay->FindBin(mult);\r
996  Int_t dim[2] = {binx,biny};\r
997  Double_t fact  = 1.0-hist->GetBinContent(dim);  \r
998 \r
999 return fact;\r
1000 }\r
1001 \r
1002 //_____________________________________________________________________________\r
1003 Long64_t AlidNdPtCorrection::Merge(TCollection* list) \r
1004 {\r
1005   // Merge list of objects (needed by PROOF)\r
1006 \r
1007   if (!list)\r
1008   return 0;\r
1009 \r
1010   if (list->IsEmpty())\r
1011   return 1;\r
1012 \r
1013   TIterator* iter = list->MakeIterator();\r
1014   TObject* obj = 0;\r
1015 \r
1016   // collection of generated histograms\r
1017 \r
1018   Int_t count=0;\r
1019   while((obj = iter->Next()) != 0) {\r
1020     AlidNdPtCorrection* entry = dynamic_cast<AlidNdPtCorrection*>(obj);\r
1021     if (entry == 0) continue; \r
1022   \r
1023     fMCEventHist1->Add(entry->fMCEventHist1);\r
1024     fRecEventHist1->Add(entry->fRecEventHist1);\r
1025 \r
1026     fMCAllEventMultHist1->Add(entry->fMCAllEventMultHist1);\r
1027     fMCAllNDEventMultHist1->Add(entry->fMCAllNDEventMultHist1);\r
1028     fMCAllNSDEventMultHist1->Add(entry->fMCAllNSDEventMultHist1);\r
1029     fMCTriggerMultHist1->Add(entry->fMCTriggerMultHist1);\r
1030     fMCEventMultHist1->Add(entry->fMCEventMultHist1);\r
1031 \r
1032     fMCAllPrimTrackMultHist1->Add(entry->fMCAllPrimTrackMultHist1);\r
1033     fMCNDEventAllPrimTrackMultHist1->Add(entry->fMCNDEventAllPrimTrackMultHist1);\r
1034     fMCNSDEventAllPrimTrackMultHist1->Add(entry->fMCNSDEventAllPrimTrackMultHist1);\r
1035     fMCTriggerPrimTrackMultHist1->Add(entry->fMCTriggerPrimTrackMultHist1);\r
1036     fMCEventPrimTrackMultHist1->Add(entry->fMCEventPrimTrackMultHist1);\r
1037 \r
1038     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {\r
1039       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);\r
1040     }\r
1041 \r
1042     for(Int_t i=0; i<8; i++) {\r
1043       fCorrRecTrackMultHist1[i]->Add(entry->fCorrRecTrackMultHist1[i]);\r
1044       fPtvsPt[i]->Add(entry->fPtvsPt[i]);\r
1045     }\r
1046 \r
1047     for(Int_t i=0; i<5; i++) {\r
1048       fCorrRecEventHist1[i]->Add(entry->fCorrRecEventHist1[i]);\r
1049       fCorrRecEventHist2[i]->Add(entry->fCorrRecEventHist2[i]);\r
1050     }\r
1051 \r
1052   count++;\r
1053   }\r
1054 \r
1055 return count;\r
1056 }\r
1057  \r
1058 Int_t AlidNdPtCorrection::GetTrueMult(THnSparse *hist, Int_t mult)\r
1059 {\r
1060  if(!hist) return 0;\r
1061  Int_t true_mult = 0;\r
1062 \r
1063  // 0 bins exluded\r
1064  TAxis *ax = hist->GetAxis(0);\r
1065  TAxis *ay = hist->GetAxis(1);\r
1066  ax->SetRange(1,ax->GetNbins());\r
1067  ay->SetRange(1,ay->GetNbins());\r
1068 \r
1069  // measured mult\r
1070  ax->SetRangeUser((Float_t)mult,(Float_t)mult); \r
1071 \r
1072  // get true multiplicity\r
1073  TH1D *h1 = (TH1D *)hist->Projection(1);\r
1074  true_mult = (Int_t)h1->GetMean();\r
1075 \r
1076  return true_mult;\r
1077 }\r
1078 \r
1079 //_____________________________________________________________________________\r
1080 void AlidNdPtCorrection::Analyse() \r
1081 {\r
1082   // Analyse histograms\r
1083   //\r
1084   TH1::AddDirectory(kFALSE);\r
1085   TH1 *h = 0, *hs=0, *hsc=0; \r
1086   TH2 *h2D = 0; \r
1087 \r
1088   TObjArray *aFolderObj = new TObjArray;\r
1089 \r
1090   //\r
1091   // get cuts\r
1092   //\r
1093   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
1094   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1095   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1096 \r
1097   if(!evtCuts || !accCuts || !esdTrackCuts) {\r
1098     Error("AlidNdPtCutAnalysis::Analyse()", "cuts not available");\r
1099     return;\r
1100   }\r
1101 \r
1102   //\r
1103   // set min and max values\r
1104   //\r
1105   Double_t minPt = accCuts->GetMinPt();\r
1106   Double_t maxPt = accCuts->GetMaxPt();\r
1107   Double_t minEta = accCuts->GetMinEta();\r
1108   Double_t maxEta = accCuts->GetMaxEta()-0.00001;\r
1109 \r
1110   //\r
1111   // pt profile\r
1112   //\r
1113   char name[256];\r
1114   for(Int_t i=0; i<8; i++) {\r
1115     h2D = fPtvsPt[i]->Projection(1,0);\r
1116     sprintf(name,"PtvsMeanPt_%d",i);\r
1117     aFolderObj->Add(h2D);\r
1118   }\r
1119 \r
1120   //\r
1121   // event level \r
1122   //\r
1123   h = fCorrRecEventHist1[0]->Projection(1);\r
1124   h->SetName("mult_event_not_corrected");\r
1125   aFolderObj->Add(h);\r
1126 \r
1127   h = fCorrRecEventHist1[1]->Projection(1);\r
1128   h->SetName("mult_event_vertex_corrected");\r
1129   aFolderObj->Add(h);\r
1130 \r
1131   h = fCorrRecEventHist1[2]->Projection(1);\r
1132   h->SetName("mult_trigger_vertex_corrected");\r
1133   aFolderObj->Add(h);\r
1134 \r
1135   h = fCorrRecEventHist1[3]->Projection(1);\r
1136   h->SetName("mult_ND_trigger_vertex_corrected");\r
1137   aFolderObj->Add(h);\r
1138 \r
1139   h = fCorrRecEventHist1[4]->Projection(1);\r
1140   h->SetName("mult_NSD_trigger_vertex_corrected");\r
1141   aFolderObj->Add(h);\r
1142 \r
1143   // empty events\r
1144   h = fCorrRecEventHist2[0]->Projection(1);\r
1145   h->SetName("mult_empty_event_not_corrected");\r
1146   aFolderObj->Add(h);\r
1147 \r
1148   h = fCorrRecEventHist2[1]->Projection(1);\r
1149   h->SetName("mult_empty_event_vertex_corrected");\r
1150   aFolderObj->Add(h);\r
1151 \r
1152   h = fCorrRecEventHist2[2]->Projection(1);\r
1153   h->SetName("mult_empty_trigger_vertex_corrected");\r
1154   aFolderObj->Add(h);\r
1155 \r
1156   h = fCorrRecEventHist2[3]->Projection(1);\r
1157   h->SetName("mult_empty_ND_trigger_vertex_corrected");\r
1158   aFolderObj->Add(h);\r
1159 \r
1160   h = fCorrRecEventHist2[4]->Projection(1);\r
1161   h->SetName("mult_empty_NSD_trigger_vertex_corrected");\r
1162   aFolderObj->Add(h);\r
1163  \r
1164   //\r
1165   // MC available\r
1166   //\r
1167   if(IsUseMCInfo()) {\r
1168 \r
1169   // mc \r
1170   h = fMCAllEventMultHist1->Projection(1);\r
1171   h->SetName("mc_mult_event_acc_prim");\r
1172   aFolderObj->Add(h);\r
1173 \r
1174   h = fMCAllNDEventMultHist1->Projection(1);\r
1175   h->SetName("mc_mult_ND_event_acc_prim");\r
1176   aFolderObj->Add(h);\r
1177 \r
1178   h = fMCAllNSDEventMultHist1->Projection(1);\r
1179   h->SetName("mc_mult_NSD_event_acc_prim");\r
1180   aFolderObj->Add(h);\r
1181 \r
1182   h = fMCTriggerMultHist1->Projection(1);\r
1183   h->SetName("mc_mult_trigger_acc_prim");\r
1184   aFolderObj->Add(h);\r
1185 \r
1186   h = fMCEventMultHist1->Projection(1);\r
1187   h->SetName("mc_mult_trigger_event_acc_prim");\r
1188   aFolderObj->Add(h);\r
1189 \r
1190 \r
1191   //\r
1192   // track level\r
1193   //\r
1194   \r
1195   // limit eta range\r
1196   for(Int_t i=0;i<8;i++) { \r
1197       //fCorrRecTrackMultHist1[i]->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1198       //fCorrRecTrackMultHist1[i]->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1199   }\r
1200   //fMCAllPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1201   //fMCAllPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1202 \r
1203   //fMCNDEventAllPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1204   //fMCNDEventAllPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1205 \r
1206   //fMCNSDEventAllPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1207   //fMCNSDEventAllPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1208 \r
1209   //fMCTriggerPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1210   //fMCTriggerPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1211 \r
1212   //fMCEventPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1213   //fMCEventPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1214 \r
1215   } // end use MC info \r
1216   \r
1217   //\r
1218   h2D = fCorrRecTrackMultHist1[3]->Projection(1,0);\r
1219   h2D->SetName("pt_eta_rec_track_mult_eff_cont_corrected");\r
1220   aFolderObj->Add(h2D);\r
1221 \r
1222   h2D = fCorrRecTrackMultHist1[4]->Projection(1,0);\r
1223   h2D->SetName("pt_eta_rec_event_track_mult_eff_cont_corrected");\r
1224   aFolderObj->Add(h2D);\r
1225 \r
1226   h2D = fCorrRecTrackMultHist1[5]->Projection(1,0);\r
1227   h2D->SetName("pt_eta_rec_trig_event_track_mult_eff_cont_corrected");\r
1228   aFolderObj->Add(h2D);\r
1229 \r
1230   h2D = fCorrRecTrackMultHist1[6]->Projection(1,0);\r
1231   h2D->SetName("pt_eta_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1232   aFolderObj->Add(h2D);\r
1233 \r
1234   h2D = fCorrRecTrackMultHist1[7]->Projection(1,0);\r
1235   h2D->SetName("pt_eta_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1236   aFolderObj->Add(h2D);\r
1237 \r
1238 \r
1239   //\r
1240   h2D = fCorrRecTrackMultHist1[3]->Projection(2,0);\r
1241   h2D->SetName("pt_mult_rec_track_mult_eff_cont_corrected");\r
1242   aFolderObj->Add(h2D);\r
1243 \r
1244   h2D = fCorrRecTrackMultHist1[4]->Projection(2,0);\r
1245   h2D->SetName("pt_mult_rec_event_track_mult_eff_cont_corrected");\r
1246   aFolderObj->Add(h2D);\r
1247 \r
1248   h2D = fCorrRecTrackMultHist1[5]->Projection(2,0);\r
1249   h2D->SetName("pt_mult_rec_trig_event_track_mult_eff_cont_corrected");\r
1250   aFolderObj->Add(h2D);\r
1251 \r
1252   h2D = fCorrRecTrackMultHist1[6]->Projection(2,0);\r
1253   h2D->SetName("pt_mult_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1254   aFolderObj->Add(h2D);\r
1255 \r
1256   h2D = fCorrRecTrackMultHist1[7]->Projection(2,0);\r
1257   h2D->SetName("pt_mult_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1258   aFolderObj->Add(h2D);\r
1259 \r
1260   // pt axis\r
1261 \r
1262   h = fCorrRecTrackMultHist1[0]->Projection(0);\r
1263   h->SetName("pt_rec_track_not_corrected");\r
1264   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1265   hs->SetName("pt_rec_track_not_corrected_s");\r
1266   aFolderObj->Add(hs);\r
1267 \r
1268   //\r
1269   h = fCorrRecTrackMultHist1[1]->Projection(0);\r
1270   h->SetName("pt_rec_track_cont_corrected");\r
1271   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1272   hs->SetName("pt_rec_track_cont_corrected_s");\r
1273   aFolderObj->Add(hs);\r
1274 \r
1275   hsc = (TH1D*)hs->Clone("pt_rec_track_cont_corr_fact");\r
1276   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_not_corrected_s"));\r
1277   aFolderObj->Add(hsc);\r
1278 \r
1279   //\r
1280   h = fCorrRecTrackMultHist1[2]->Projection(0);\r
1281   h->SetName("pt_rec_track_eff_cont_corrected");\r
1282   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1283   hs->SetName("pt_rec_track_eff_cont_corrected_s");\r
1284   aFolderObj->Add(hs);\r
1285 \r
1286   hsc = (TH1D*)hs->Clone("pt_rec_track_eff_corr_fact");\r
1287   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_cont_corrected_s"));\r
1288   aFolderObj->Add(hsc);\r
1289 \r
1290   //\r
1291   h = fCorrRecTrackMultHist1[3]->Projection(0);\r
1292   h->SetName("pt_rec_track_mult_eff_cont_corrected");\r
1293   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1294   hs->SetName("pt_rec_track_mult_eff_cont_corrected_s");\r
1295   aFolderObj->Add(hs);\r
1296 \r
1297   hsc = (TH1D*)hs->Clone("pt_rec_track_mult_corr_fact");\r
1298   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_eff_cont_corrected_s"));\r
1299   aFolderObj->Add(hsc);\r
1300 \r
1301   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1302   hsc = (TH1D*)hs->Clone("pt_rec_track_mult_eff_cont_corr_fact");\r
1303   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_not_corrected_s"));\r
1304   aFolderObj->Add(hsc);\r
1305 \r
1306   hsc = (TH1D*)hs->Clone();\r
1307   hsc->SetName("pt_rec_track_mult_eff_cont_corrected_s_norm");\r
1308   hsc->Scale(1./(fCorrRecEventHist1[0]->Projection(1)->Integral()));\r
1309   aFolderObj->Add(hsc);\r
1310 \r
1311   //\r
1312   h = fCorrRecTrackMultHist1[4]->Projection(0);\r
1313   h->SetName("pt_rec_event_track_mult_eff_cont_corrected");\r
1314   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1315   hs->SetName("pt_rec_event_track_mult_eff_cont_corrected_s");\r
1316   aFolderObj->Add(hs);\r
1317 \r
1318   hsc = (TH1D*)hs->Clone();\r
1319   hsc->SetName("pt_rec_event_track_mult_eff_cont_corrected_s_norm");\r
1320   hsc->Scale(1./(fCorrRecEventHist1[1]->Projection(1)->Integral()+fCorrRecEventHist2[1]->Projection(1)->Integral()));\r
1321   aFolderObj->Add(hsc);\r
1322 \r
1323   //\r
1324   h = fCorrRecTrackMultHist1[5]->Projection(0);\r
1325   h->SetName("pt_rec_trig_event_track_mult_eff_cont_corrected");\r
1326   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1327   hs->SetName("pt_rec_trig_event_track_mult_eff_cont_corrected_s");\r
1328   aFolderObj->Add(hs);\r
1329 \r
1330   hsc = (TH1D*)hs->Clone();\r
1331   hsc->SetName("pt_rec_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1332   hsc->Scale(1./(fCorrRecEventHist1[2]->Projection(1)->Integral() + fCorrRecEventHist2[2]->Projection(1)->Integral()));\r
1333   aFolderObj->Add(hsc);\r
1334 \r
1335   //\r
1336   h = fCorrRecTrackMultHist1[6]->Projection(0);\r
1337   h->SetName("pt_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1338   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1339   hs->SetName("pt_rec_ND_trig_event_track_mult_eff_cont_corrected_s");\r
1340   aFolderObj->Add(hs);\r
1341 \r
1342   hsc = (TH1D*)hs->Clone();\r
1343   hsc->SetName("pt_rec_ND_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1344   hsc->Scale(1./(fCorrRecEventHist1[3]->Projection(1)->Integral()+fCorrRecEventHist2[3]->Projection(1)->Integral()));\r
1345   aFolderObj->Add(hsc);\r
1346 \r
1347   //\r
1348   h = fCorrRecTrackMultHist1[7]->Projection(0);\r
1349   h->SetName("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1350   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1351   hs->SetName("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected_s");\r
1352   aFolderObj->Add(hs);\r
1353 \r
1354   hsc = (TH1D*)hs->Clone();\r
1355   hsc->SetName("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1356   hsc->Scale(1./(fCorrRecEventHist1[4]->Projection(1)->Integral()+fCorrRecEventHist2[4]->Projection(1)->Integral()));\r
1357   aFolderObj->Add(hsc);\r
1358 \r
1359   // eta axis\r
1360   h = fCorrRecTrackMultHist1[0]->Projection(1);\r
1361   h->SetName("eta_rec_track_not_corrected");\r
1362   aFolderObj->Add(h);\r
1363   \r
1364   h = fCorrRecTrackMultHist1[1]->Projection(1);\r
1365   h->SetName("eta_rec_track_cont_corrected");\r
1366   aFolderObj->Add(h);\r
1367 \r
1368   h = fCorrRecTrackMultHist1[2]->Projection(1);\r
1369   h->SetName("eta_rec_track_eff_cont_corrected");\r
1370   aFolderObj->Add(h);\r
1371 \r
1372   h = fCorrRecTrackMultHist1[3]->Projection(1);\r
1373   h->SetName("eta_rec_track_mult_eff_cont_corrected");\r
1374   aFolderObj->Add(h);\r
1375 \r
1376   h = fCorrRecTrackMultHist1[4]->Projection(1);\r
1377   h->SetName("eta_rec_event_track_mult_eff_cont_corrected");\r
1378   aFolderObj->Add(h);\r
1379 \r
1380   h = fCorrRecTrackMultHist1[5]->Projection(1);\r
1381   h->SetName("eta_rec_trig_event_track_mult_eff_cont_corrected");\r
1382   aFolderObj->Add(h);\r
1383 \r
1384   h = fCorrRecTrackMultHist1[6]->Projection(1);\r
1385   h->SetName("eta_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1386   aFolderObj->Add(h);\r
1387 \r
1388   h = fCorrRecTrackMultHist1[7]->Projection(1);\r
1389   h->SetName("eta_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1390   aFolderObj->Add(h);\r
1391 \r
1392 \r
1393   //\r
1394   // MC available\r
1395   //\r
1396   if(IsUseMCInfo()) {\r
1397 \r
1398   //\r
1399   h2D = fMCAllPrimTrackMultHist1->Projection(2,0);\r
1400   h2D->SetName("mc_all_pt_mult_acc_prim");\r
1401   aFolderObj->Add(h2D);\r
1402 \r
1403   h2D = fMCAllPrimTrackMultHist1->Projection(1,0);\r
1404   h2D->SetName("mc_all_eta_pt_acc_prim");\r
1405   aFolderObj->Add(h2D);\r
1406 \r
1407   h2D = fMCNDEventAllPrimTrackMultHist1->Projection(2,0);\r
1408   h2D->SetName("mc_ND_all_pt_mult_acc_prim");\r
1409   aFolderObj->Add(h2D);\r
1410 \r
1411   h2D = fMCNDEventAllPrimTrackMultHist1->Projection(1,0);\r
1412   h2D->SetName("mc_ND_all_eta_pt_acc_prim");\r
1413   aFolderObj->Add(h2D);\r
1414 \r
1415   h2D = fMCNSDEventAllPrimTrackMultHist1->Projection(2,0);\r
1416   h2D->SetName("mc_NSD_all_pt_mult_acc_prim");\r
1417   aFolderObj->Add(h2D);\r
1418 \r
1419   h2D = fMCNSDEventAllPrimTrackMultHist1->Projection(1,0);\r
1420   h2D->SetName("mc_NSD_all_eta_pt_acc_prim");\r
1421   aFolderObj->Add(h2D);\r
1422 \r
1423   //\r
1424 \r
1425   h = fMCAllPrimTrackMultHist1->Projection(0);\r
1426   h->SetName("mc_all_pt_acc_prim");\r
1427   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1428   hs->SetName("mc_all_pt_acc_prim_s");\r
1429   aFolderObj->Add(hs);\r
1430 \r
1431   hsc = (TH1D*)hs->Clone();\r
1432   hsc->SetName("mc_all_pt_acc_prim_s_norm");\r
1433   hsc->Scale(1./fMCAllEventMultHist1->Projection(1)->Integral());\r
1434   aFolderObj->Add(hsc);\r
1435 \r
1436   h = fMCNDEventAllPrimTrackMultHist1->Projection(0);\r
1437   h->SetName("mc_ND_all_pt_acc_prim");\r
1438   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1439   hs->SetName("mc_ND_all_pt_acc_prim_s");\r
1440   aFolderObj->Add(hs);\r
1441 \r
1442   hsc = (TH1D*)hs->Clone();\r
1443   hsc->SetName("mc_ND_all_pt_acc_prim_s_norm");\r
1444   hsc->Scale(1./fMCAllNDEventMultHist1->Projection(1)->Integral());\r
1445   aFolderObj->Add(hsc);\r
1446 \r
1447   h = fMCNSDEventAllPrimTrackMultHist1->Projection(0);\r
1448   h->SetName("mc_NSD_all_pt_acc_prim");\r
1449   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1450   hs->SetName("mc_NSD_all_pt_acc_prim_s");\r
1451   aFolderObj->Add(hs);\r
1452 \r
1453   hsc = (TH1D*)hs->Clone();\r
1454   hsc->SetName("mc_NSD_all_pt_acc_prim_s_norm");\r
1455   hsc->Scale(1./fMCAllNSDEventMultHist1->Projection(1)->Integral());\r
1456   aFolderObj->Add(hsc);\r
1457 \r
1458   h = fMCTriggerPrimTrackMultHist1->Projection(0);\r
1459   h->SetName("mc_trigger_all_pt_acc_prim");\r
1460   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1461   hs->SetName("mc_trigger_all_pt_acc_prim_s");\r
1462   aFolderObj->Add(hs);\r
1463 \r
1464   hsc = (TH1D*)hs->Clone();\r
1465   hsc->SetName("mc_trigger_all_pt_acc_prim_s_norm");\r
1466   hsc->Scale(1./fMCTriggerMultHist1->Projection(1)->Integral());\r
1467   aFolderObj->Add(hsc);\r
1468 \r
1469   h = fMCEventPrimTrackMultHist1->Projection(0);\r
1470   h->SetName("mc_all_pt_acc_trig_event_prim");\r
1471   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1472   hs->SetName("mc_all_pt_acc_trig_event_prim_s");\r
1473   aFolderObj->Add(hs);\r
1474 \r
1475   hsc = (TH1D*)hs->Clone();\r
1476   hsc->SetName("mc_all_pt_acc_trig_event_prim_s_norm");\r
1477   hsc->Scale(1./fMCEventMultHist1->Projection(1)->Integral());\r
1478   aFolderObj->Add(hsc);\r
1479 \r
1480   //\r
1481 \r
1482   h = fMCAllPrimTrackMultHist1->Projection(1);\r
1483   h->SetName("mc_all_eta_acc_prim");\r
1484   aFolderObj->Add(h);\r
1485 \r
1486   h = fMCNDEventAllPrimTrackMultHist1->Projection(1);\r
1487   h->SetName("mc_ND_all_eta_acc_prim");\r
1488   aFolderObj->Add(h);\r
1489 \r
1490   h = fMCNSDEventAllPrimTrackMultHist1->Projection(1);\r
1491   h->SetName("mc_NSD_all_eta_acc_prim");\r
1492   aFolderObj->Add(h);\r
1493 \r
1494   h = fMCTriggerPrimTrackMultHist1->Projection(1);\r
1495   h->SetName("mc_trigger_all_eta_acc_prim");\r
1496   aFolderObj->Add(h);\r
1497 \r
1498   h = fMCEventPrimTrackMultHist1->Projection(1);\r
1499   h->SetName("mc_all_eta_acc_trig_event_prim");\r
1500   aFolderObj->Add(h);\r
1501 \r
1502   //\r
1503   // calculate ratios (rec / mc)\r
1504   //\r
1505   \r
1506   hs = (TH1*)aFolderObj->FindObject("pt_rec_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1507   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_trig_event_track_mult_eff_cont_corrected"); \r
1508   hsc->Sumw2();\r
1509   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_pt_acc_prim_s_norm"));\r
1510   aFolderObj->Add(hsc);\r
1511 \r
1512   hs = (TH1*)aFolderObj->FindObject("eta_rec_trig_event_track_mult_eff_cont_corrected");\r
1513   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_trig_event_track_mult_eff_cont_corrected"); \r
1514   hsc->Sumw2();\r
1515   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_eta_acc_prim"));\r
1516   aFolderObj->Add(hsc);\r
1517 \r
1518   //\r
1519   hs = (TH1*)aFolderObj->FindObject("pt_rec_ND_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1520   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_ND_trig_event_track_mult_eff_cont_corrected"); \r
1521   hsc->Sumw2();\r
1522   hsc->Divide((TH1*)aFolderObj->FindObject("mc_ND_all_pt_acc_prim_s_norm"));\r
1523   aFolderObj->Add(hsc);\r
1524 \r
1525   hs = (TH1*)aFolderObj->FindObject("eta_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1526   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_ND_trig_event_track_mult_eff_cont_corrected"); \r
1527   hsc->Sumw2();\r
1528   hsc->Divide((TH1*)aFolderObj->FindObject("mc_ND_all_eta_acc_prim"));\r
1529   aFolderObj->Add(hsc);\r
1530 \r
1531   //\r
1532   hs = (TH1*)aFolderObj->FindObject("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1533   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_NSD_trig_event_track_mult_eff_cont_corrected"); \r
1534   hsc->Sumw2();\r
1535   hsc->Divide((TH1*)aFolderObj->FindObject("mc_NSD_all_pt_acc_prim_s_norm"));\r
1536   aFolderObj->Add(hsc);\r
1537 \r
1538   hs = (TH1*)aFolderObj->FindObject("eta_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1539   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_NSD_trig_event_track_mult_eff_cont_corrected"); \r
1540   hsc->Sumw2();\r
1541   hsc->Divide((TH1*)aFolderObj->FindObject("mc_NSD_all_eta_acc_prim"));\r
1542   aFolderObj->Add(hsc);\r
1543 \r
1544   //\r
1545   hs = (TH1*)aFolderObj->FindObject("pt_rec_event_track_mult_eff_cont_corrected_s_norm");\r
1546   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_event_track_mult_eff_cont_corrected"); \r
1547   hsc->Sumw2();\r
1548   hsc->Divide((TH1*)aFolderObj->FindObject("mc_trigger_all_pt_acc_prim_s_norm"));\r
1549   aFolderObj->Add(hsc);\r
1550 \r
1551   hs = (TH1*)aFolderObj->FindObject("eta_rec_event_track_mult_eff_cont_corrected");\r
1552   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_event_track_mult_eff_cont_corrected"); \r
1553   hsc->Sumw2();\r
1554   hsc->Divide((TH1*)aFolderObj->FindObject("mc_trigger_all_eta_acc_prim"));\r
1555   aFolderObj->Add(hsc);\r
1556 \r
1557   // track level\r
1558   hs = (TH1*)aFolderObj->FindObject("pt_rec_track_mult_eff_cont_corrected_s_norm");\r
1559   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_track_mult_eff_cont_corrected"); \r
1560   hsc->Sumw2();\r
1561   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_pt_acc_trig_event_prim_s_norm"));\r
1562   aFolderObj->Add(hsc);\r
1563 \r
1564   hs = (TH1*)aFolderObj->FindObject("eta_rec_track_mult_eff_cont_corrected");\r
1565   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_track_mult_eff_cont_corrected"); \r
1566   hsc->Sumw2();\r
1567   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_eta_acc_trig_event_prim"));\r
1568   aFolderObj->Add(hsc);\r
1569 \r
1570   } // end MC infor available\r
1571 \r
1572   // export objects to analysis folder\r
1573   fCorrectionFolder = ExportToFolder(aFolderObj);\r
1574 \r
1575   // delete only TObjArray\r
1576   if(aFolderObj) delete aFolderObj;\r
1577 }\r
1578 \r
1579 //_____________________________________________________________________________\r
1580 TFolder* AlidNdPtCorrection::ExportToFolder(TObjArray * array) \r
1581 {\r
1582   // recreate folder avery time and export objects to new one\r
1583   //\r
1584   AlidNdPtCorrection * comp=this;\r
1585   TFolder *folder = comp->GetCorrectionFolder();\r
1586 \r
1587   TString name, title;\r
1588   TFolder *newFolder = 0;\r
1589   Int_t i = 0;\r
1590   Int_t size = array->GetSize();\r
1591 \r
1592   if(folder) { \r
1593      // get name and title from old folder\r
1594      name = folder->GetName();  \r
1595      title = folder->GetTitle();  \r
1596 \r
1597          // delete old one\r
1598      delete folder;\r
1599 \r
1600          // create new one\r
1601      newFolder = CreateFolder(name.Data(),title.Data());\r
1602      newFolder->SetOwner();\r
1603 \r
1604          // add objects to folder\r
1605      while(i < size) {\r
1606            newFolder->Add(array->At(i));\r
1607            i++;\r
1608          }\r
1609   }\r
1610 \r
1611 return newFolder;\r
1612 }\r
1613 \r
1614 //_____________________________________________________________________________\r
1615 TFolder* AlidNdPtCorrection::CreateFolder(TString name,TString title) { \r
1616 // create folder for analysed histograms\r
1617 //\r
1618 TFolder *folder = 0;\r
1619   folder = new TFolder(name.Data(),title.Data());\r
1620 \r
1621   return folder;\r
1622 }\r