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