trigger and new tracking scenario
[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); \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 vertex contributors\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::kTPCSPDvtxUpdate || GetAnalysisMode() == AlidNdPtHelper::kMCRec) {\r
613     //multMBTracks = AlidNdPtHelper::GetSPDMBTrackMult(esdEvent,0.0);\r
614     multMBTracks = vtxESD->GetNContributors();\r
615   } \r
616   else {\r
617     AliDebug(AliLog::kError, Form("Found analysis type %s", GetAnalysisMode()));\r
618     return; \r
619   }\r
620 \r
621   //\r
622   // correct event and track histograms\r
623   //\r
624   TObjArray *allChargedTracks=0;\r
625   Int_t multRec=0, multRecTemp=0;\r
626   Int_t *labelsRec=0;\r
627 \r
628   if(isEventOK && isEventTriggered)\r
629   {\r
630     // get all charged tracks\r
631     //allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,vtxESD,GetAnalysisMode());\r
632     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());\r
633     if(!allChargedTracks) return;\r
634 \r
635     Int_t entries = allChargedTracks->GetEntries();\r
636     labelsRec = new Int_t[entries];\r
637 \r
638     // calculate mult of reconstructed tracks\r
639     for(Int_t i=0; i<entries;++i) \r
640     {\r
641       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);\r
642       if(!track) continue;\r
643       if(track->Charge()==0) continue;\r
644 \r
645       // only postive charged \r
646       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0) \r
647         continue;\r
648       \r
649       // only negative charged \r
650       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0) \r
651         continue;\r
652 \r
653       if(esdTrackCuts->AcceptTrack(track) && accCuts->AcceptTrack(track)) \r
654         multRecTemp++;\r
655     }  \r
656 \r
657     //\r
658     for(Int_t i=0; i<entries;++i) \r
659     {\r
660       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);\r
661       if(!track) continue;\r
662       if(track->Charge()==0) continue;\r
663 \r
664       // only postive charged \r
665       if(GetParticleMode() == AlidNdPtHelper::kPlus && track->Charge() < 0) \r
666       continue;\r
667       \r
668       // only negative charged \r
669       if(GetParticleMode() == AlidNdPtHelper::kMinus && track->Charge() > 0) \r
670       continue;\r
671         \r
672       // track-level corrections\r
673       if(esdTrackCuts->AcceptTrack(track)) \r
674       {\r
675          if (GetAnalysisMode()==AlidNdPtHelper::kTPCSPDvtxUpdate) \r
676          {\r
677            //\r
678            // update track parameters\r
679            //\r
680            AliExternalTrackParam cParam;\r
681            track->RelateToVertexTPC(esdEvent->GetPrimaryVertexSPD(),esdEvent->GetMagneticField(),kVeryBig,&cParam);\r
682            track->Set(cParam.GetX(),cParam.GetAlpha(),cParam.GetParameter(),cParam.GetCovariance());\r
683 \r
684            if(accCuts->AcceptTrack(track)) { \r
685              FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,vtxESD->GetZv(),multRecTemp); \r
686              labelsRec[multRec] = TMath::Abs(track->GetLabel());\r
687              multRec++;\r
688            }\r
689          }\r
690          else if (GetAnalysisMode()==AlidNdPtHelper::kMCRec && IsUseMCInfo()) \r
691          { \r
692            //\r
693            // Replace rec with MC\r
694            //\r
695            if(accCuts->AcceptTrack(track)) { \r
696              FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,vtxMC[2],multRecTemp); \r
697              labelsRec[multRec] = TMath::Abs(track->GetLabel());\r
698              multRec++;\r
699            }\r
700          } \r
701          else  {\r
702            //\r
703            // all the rest tracking scenarios\r
704            //\r
705            if(accCuts->AcceptTrack(track)) { \r
706              FillHistograms(track,stack,AlidNdPtHelper::kRecTracks,vtxESD->GetZv(),multRecTemp); \r
707              labelsRec[multRec] = TMath::Abs(track->GetLabel());\r
708              multRec++;\r
709            }\r
710          }\r
711        }\r
712     }\r
713     // event-level corrections\r
714     if(GetAnalysisMode()==AlidNdPtHelper::kMCRec && IsUseMCInfo()) { \r
715       FillHistograms(AlidNdPtHelper::kRecEvents,vtxMC[2],multMBTracks);\r
716     }\r
717     else {\r
718       FillHistograms(AlidNdPtHelper::kRecEvents,vtxESD->GetZv(),multMBTracks);\r
719     }\r
720 \r
721     // control event histograms\r
722     Double_t vRecEventHist1[3] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv()};\r
723     fRecEventHist1->Fill(vRecEventHist1);\r
724   } \r
725 \r
726   // empty events corrections\r
727   // no reconstructed zv\r
728   if(isEventTriggered && multMBTracks==0) {\r
729     if(GetAnalysisMode()==AlidNdPtHelper::kMCRec && IsUseMCInfo()) {\r
730       FillHistograms(AlidNdPtHelper::kTriggeredEvents,vtxMC[2],multMBTracks);\r
731     }\r
732     else {\r
733       Double_t zv = fZvNorm->GetRandom();\r
734       FillHistograms(AlidNdPtHelper::kTriggeredEvents,zv,multMBTracks);\r
735     }\r
736   }\r
737 \r
738   if(IsUseMCInfo())  \r
739   {\r
740     // select MC events \r
741     if(evtCuts->AcceptMCEvent(mcEvent))\r
742     {\r
743       //\r
744       // event histograms\r
745       //\r
746       Double_t vMCEventMatrix[2] = {vtxMC[2],multMBTracks};\r
747       fMCAllEventMultHist1->Fill(vMCEventMatrix);\r
748       if(evtType == AliPWG0Helper::kND) {\r
749         fMCAllNDEventMultHist1->Fill(vMCEventMatrix);\r
750       }\r
751       if(evtType != AliPWG0Helper::kSD) {\r
752         fMCAllNSDEventMultHist1->Fill(vMCEventMatrix);\r
753       }\r
754       if(isEventTriggered) fMCTriggerMultHist1->Fill(vMCEventMatrix);\r
755       if(isEventTriggered && isEventOK) fMCEventMultHist1->Fill(vMCEventMatrix);\r
756 \r
757       //\r
758       // MC histograms for efficiency studies\r
759       //\r
760       Int_t nPart  = stack->GetNtrack();\r
761       for (Int_t iMc = 0; iMc < nPart; ++iMc) \r
762       {\r
763         // print MC stack info\r
764         //AlidNdPtHelper::PrintMCInfo(stack,iMc);\r
765 \r
766         TParticle* particle = stack->Particle(iMc);\r
767         if (!particle)\r
768         continue;\r
769 \r
770         // only charged particles\r
771         Double_t charge = particle->GetPDG()->Charge()/3.;\r
772         if (charge == 0.0)\r
773           continue;\r
774 \r
775         // only postive charged \r
776         if(GetParticleMode() == AlidNdPtHelper::kPlus && charge  < 0.) \r
777         continue;\r
778       \r
779         // only negative charged \r
780         if(GetParticleMode() == AlidNdPtHelper::kMinus && charge > 0.) \r
781         continue;\r
782       \r
783         // physical primary\r
784         Bool_t prim = stack->IsPhysicalPrimary(iMc);\r
785 \r
786         // all primaries in acceptance\r
787         if(!accCuts->AcceptTrack(particle)) continue;\r
788         if(!prim) continue;\r
789 \r
790         Double_t gpt = particle->Pt();\r
791         Double_t geta = particle->Eta();\r
792 \r
793         Double_t valMCAllTrackMultHist1[3] = {gpt,geta,multRec};          \r
794         fMCAllPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
795         if(evtType == AliPWG0Helper::kND) {\r
796           fMCNDEventAllPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
797         }\r
798         if(evtType != AliPWG0Helper::kSD) {\r
799           fMCNSDEventAllPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
800         }\r
801         if(isEventTriggered) fMCTriggerPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
802         if(isEventTriggered && isEventOK) fMCEventPrimTrackMultHist1->Fill(valMCAllTrackMultHist1);\r
803       }\r
804     }\r
805   } // end bUseMC\r
806 \r
807   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;\r
808   if(labelsRec) delete [] labelsRec; labelsRec = 0;\r
809 \r
810   if(!evtCuts->IsRecVertexRequired() && vtxESD != NULL) delete vtxESD;\r
811 }\r
812 \r
813 //_____________________________________________________________________________\r
814 void AlidNdPtCorrection::FillHistograms(AlidNdPtHelper::EventObject eventObj, Double_t zv, Int_t multMBTracks)\r
815 {\r
816   //\r
817   // Fill corrected histograms\r
818   //\r
819 \r
820   Double_t vEventMatrix[2] = {zv,multMBTracks};\r
821   //\r
822   // Correct for efficiency \r
823   //\r
824   if(eventObj == AlidNdPtHelper::kRecEvents && multMBTracks>0)  \r
825   {\r
826     Double_t corrToMBF = GetCorrFactZvMult(fCorrEventMatrix,zv,multMBTracks);\r
827     Double_t corrToInelF = GetCorrFactZvMult(fCorrTriggerMBtoInelEventMatrix,zv,multMBTracks);\r
828     Double_t corrToNDF = GetCorrFactZvMult(fCorrTriggerMBtoNDEventMatrix,zv,multMBTracks);\r
829     Double_t corrToNSDF = GetCorrFactZvMult(fCorrTriggerMBtoNSDEventMatrix,zv,multMBTracks);\r
830     //printf("corrToMBF %f, corrToInelF %f, corrToNDF %f corrToNSDF %f \n",corrToMBF,corrToInelF,corrToNDF,corrToNSDF);\r
831 \r
832     fCorrRecEventHist1[0]->Fill(vEventMatrix);\r
833     fCorrRecEventHist1[1]->Fill(vEventMatrix,corrToMBF);\r
834     fCorrRecEventHist1[2]->Fill(vEventMatrix,corrToMBF*corrToInelF);\r
835     fCorrRecEventHist1[3]->Fill(vEventMatrix,corrToMBF*corrToNDF);\r
836     fCorrRecEventHist1[4]->Fill(vEventMatrix,corrToMBF*corrToNSDF);\r
837   }\r
838 \r
839   if(eventObj==AlidNdPtHelper::kTriggeredEvents && multMBTracks==0) // empty triggered events\r
840   {\r
841     Int_t bin = fZvEmptyEventsNorm->FindBin(zv);\r
842     Double_t Fz = fZvEmptyEventsNorm->GetBinContent(bin);\r
843     Double_t corrToInelF0 = GetCorrFactZvMult(fCorrTriggerMBtoInelEventMatrix,zv,multMBTracks);\r
844     Double_t corrToNDF0 = GetCorrFactZvMult(fCorrTriggerMBtoNDEventMatrix,zv,multMBTracks);\r
845     Double_t corrToNSDF0 = GetCorrFactZvMult(fCorrTriggerMBtoNSDEventMatrix,zv,multMBTracks);\r
846     //printf("Fz %f, corrToInelF0 %f, corrToNDF0 %f, corrToNSDF0 %f \n",Fz,corrToInelF0,corrToNDF0,corrToNSDF0);\r
847 \r
848     fCorrRecEventHist2[0]->Fill(vEventMatrix);\r
849     fCorrRecEventHist2[1]->Fill(vEventMatrix,Fz);\r
850     fCorrRecEventHist2[2]->Fill(vEventMatrix,Fz*corrToInelF0);\r
851     fCorrRecEventHist2[3]->Fill(vEventMatrix,Fz*corrToNDF0);\r
852     fCorrRecEventHist2[4]->Fill(vEventMatrix,Fz*corrToNSDF0);\r
853   }\r
854 }\r
855 \r
856 //_____________________________________________________________________________\r
857 void AlidNdPtCorrection::FillHistograms(AliESDtrack *esdTrack, AliStack *stack, AlidNdPtHelper::TrackObject trackObj, Double_t zv, Int_t mult)\r
858 {\r
859   //\r
860   // Fill ESD track and MC histograms \r
861   //\r
862   if(!esdTrack) return;\r
863 \r
864   //Float_t q = esdTrack->Charge();\r
865   Float_t pt = esdTrack->Pt();\r
866   Float_t eta = esdTrack->Eta();\r
867   Float_t phi = esdTrack->Phi();\r
868 \r
869   if(stack && GetAnalysisMode() == AlidNdPtHelper::kMCRec) \r
870   {\r
871     Int_t label = TMath::Abs(esdTrack->GetLabel());\r
872    \r
873     TParticle* particle = stack->Particle(label);\r
874     if(!particle) return;\r
875    \r
876     Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3\r
877     if(gq==0) return;\r
878     Float_t gpt = particle->Pt();\r
879     Float_t geta = particle->Eta();\r
880     Float_t gphi = particle->Phi();\r
881 \r
882     // replace reconstructed values with MC\r
883     pt = gpt;\r
884     eta = geta;\r
885     phi = gphi;\r
886   }\r
887 \r
888   //\r
889   // Fill histograms\r
890   //\r
891   Double_t values[3] = {pt,eta,phi};      \r
892   fRecTrackHist1[trackObj]->Fill(values);\r
893 \r
894   //\r
895   // Correct for contamination and efficiency \r
896   //\r
897   if(trackObj == AlidNdPtHelper::kRecTracks || GetAnalysisMode() == AlidNdPtHelper::kMCRec)  \r
898   {\r
899     // track level corrections\r
900     Double_t trackEffF = GetCorrFactZvPtEta(fCorrTrackMatrix,zv,pt,eta);\r
901     Double_t trackContF = GetContFactZvPtEta(fContTrackMatrix,zv,pt,eta);\r
902     Double_t multTrackContF = GetContFactZvPtEta(fContMultTrackMatrix,zv,pt,eta);\r
903     //printf("zv %f, pt %f, eta %f \n",zv,pt,eta);\r
904     //printf("trackEffF %f, trackContF %f, multTrackContF %f \n", trackEffF, trackContF, multTrackContF);\r
905    \r
906     // track-event level corrections\r
907     Double_t vertexEffF = GetCorrFactZvPtEta(fCorrTrackEventMatrix,zv,pt,eta);\r
908     Double_t trigMBToInel = GetCorrFactZvPtEta(fCorrTriggerMBtoInelTrackEventMatrix,zv,pt,eta);  \r
909     Double_t trigMBToND = GetCorrFactZvPtEta(fCorrTriggerMBtoNDTrackEventMatrix,zv,pt,eta);\r
910     Double_t trigMBToNSD = GetCorrFactZvPtEta(fCorrTriggerMBtoNSDTrackEventMatrix,zv,pt,eta);\r
911     //printf("vertexEffF %f, trigMBToInel %f, trigMBToNSD %f \n", vertexEffF, trigMBToInel, trigMBToNSD);\r
912     \r
913     Double_t corrF[8] = { 1.0, \r
914                           trackContF,\r
915                           trackContF*trackEffF,\r
916                           trackContF*trackEffF*multTrackContF,\r
917                           trackContF*trackEffF*multTrackContF*vertexEffF,\r
918                           trackContF*trackEffF*multTrackContF*vertexEffF*trigMBToInel,\r
919                           trackContF*trackEffF*multTrackContF*vertexEffF*trigMBToND,\r
920                           trackContF*trackEffF*multTrackContF*vertexEffF*trigMBToNSD\r
921                          }; \r
922  \r
923     // Fill histograms\r
924     Double_t valCorrRecTrackMultHist1[3] = {pt,eta,mult};         \r
925     Double_t valPtvsPt[2] = {pt,pt};      \r
926     for(Int_t i=0; i<8; i++) {\r
927       fCorrRecTrackMultHist1[i]->Fill(valCorrRecTrackMultHist1,corrF[i]);\r
928       fPtvsPt[i]->Fill(valPtvsPt,corrF[i]);\r
929     }\r
930   }\r
931 }\r
932 \r
933 //_____________________________________________________________________________\r
934 void AlidNdPtCorrection::FillHistograms(AliStack *stack, Int_t /*label*/, AlidNdPtHelper::TrackObject /*trackObj*/, Int_t /*mult*/)\r
935 {\r
936   // Fill MC histograms\r
937   if(!stack) return;\r
938 \r
939   /*\r
940   TParticle* particle = stack->Particle(label);\r
941   if(!particle) return;\r
942 \r
943   Int_t mother_pdg = -1;\r
944   TParticle* mother = 0;\r
945 \r
946   //TParticle* prim_mother = AlidNdPtHelper::FindPrimaryMother(stack,label);\r
947   Int_t motherLabel = particle->GetMother(0); \r
948   if(motherLabel>0) mother = stack->Particle(motherLabel);\r
949   if(mother) mother_pdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only\r
950   Int_t mech = particle->GetUniqueID(); // production mechanism\r
951 \r
952   Double_t gq = particle->GetPDG()->Charge()/3.0; // Charge units |e|/3 \r
953   Float_t gpt = particle->Pt();\r
954   Float_t qgpt = particle->Pt() * gq;\r
955   Float_t geta = particle->Eta();\r
956   Float_t gphi = particle->Phi();\r
957   Float_t gpz = particle->Pz();\r
958 \r
959   Bool_t prim = stack->IsPhysicalPrimary(label);\r
960   Float_t vx = particle->Vx(); Float_t vy = particle->Vy(); Float_t vz = particle->Vz();\r
961 \r
962   Int_t pid=-1;\r
963   if (TMath::Abs(particle->GetPdgCode()) == kElectron)         { pid = 0; }\r
964     else if (TMath::Abs(particle->GetPdgCode()) == kMuonMinus) { pid = 1; }\r
965     else if (TMath::Abs(particle->GetPdgCode()) == kPiPlus)    { pid = 2; }\r
966     else if (TMath::Abs(particle->GetPdgCode()) == kKPlus)     { pid = 3; }\r
967     else if (TMath::Abs(particle->GetPdgCode()) == kProton)    { pid = 4; }\r
968     else                                                       { pid = 5; }\r
969     */\r
970 \r
971   //if(!prim) printf("prim_mother %d, mother %d, particle %d, production mech %d\n",prim_mother->GetPdgCode(),mother->GetPdgCode(), particle->GetPdgCode(),mech);\r
972   \r
973 }\r
974 \r
975 //_____________________________________________________________________________\r
976 Double_t AlidNdPtCorrection::GetCorrFactZvPtEta(THnSparse *hist, Double_t zv, Double_t pt, Double_t eta) const {\r
977 // return correction factor F(zv,pt,eta)\r
978 \r
979  if(!hist) return 1.;\r
980 \r
981  //\r
982  TAxis *ax = hist->GetAxis(0);\r
983  TAxis *ay = hist->GetAxis(1);\r
984  TAxis *az = hist->GetAxis(2);\r
985 \r
986  Int_t binx = ax->FindBin(zv);\r
987  Int_t biny = ay->FindBin(pt);\r
988  Int_t binz = az->FindBin(eta);\r
989  Int_t dim[3] = {binx,biny,binz};\r
990 \r
991  Double_t fact  = hist->GetBinContent(dim);  \r
992 \r
993 return fact;\r
994 }\r
995 \r
996 //_____________________________________________________________________________\r
997 Double_t AlidNdPtCorrection::GetContFactZvPtEta(THnSparse *hist, Double_t zv, Double_t pt, Double_t eta) const {\r
998 // return contamination correction factor F(zv,pt,eta)\r
999 \r
1000  if(!hist) return 1.0;\r
1001 \r
1002  //\r
1003  TAxis *ax = hist->GetAxis(0);\r
1004  TAxis *ay = hist->GetAxis(1);\r
1005  TAxis *az = hist->GetAxis(2);\r
1006 \r
1007  Int_t binx = ax->FindBin(zv);\r
1008  Int_t biny = ay->FindBin(pt);\r
1009  Int_t binz = az->FindBin(eta);\r
1010  Int_t dim[3] = {binx,biny,binz};\r
1011 \r
1012  Double_t fact  = 1.0-hist->GetBinContent(dim);  \r
1013  //Double_t fact  = hist->GetBinContent(dim);  \r
1014 \r
1015 return fact;\r
1016 }\r
1017 \r
1018 //_____________________________________________________________________________\r
1019 Double_t AlidNdPtCorrection::GetCorrFactZvMult(THnSparse *hist, Double_t zv, Int_t mult) const {\r
1020 // return correction factor F(zv,mult)\r
1021 \r
1022  if(!hist) return 1.;\r
1023 \r
1024  TAxis *ax = hist->GetAxis(0);\r
1025  TAxis *ay = hist->GetAxis(1);\r
1026  Int_t binx = ax->FindBin(zv);\r
1027  Int_t biny = ay->FindBin(mult);\r
1028  Int_t dim[2] = {binx,biny};\r
1029 \r
1030  Double_t fact  = hist->GetBinContent(dim);  \r
1031 \r
1032 \r
1033 return fact;\r
1034 }\r
1035 \r
1036 //_____________________________________________________________________________\r
1037 Double_t AlidNdPtCorrection::GetContFactZvMult(THnSparse *hist, Double_t zv, Int_t mult) const {\r
1038 // return contamination correction factor F(zv,mult)\r
1039 \r
1040  if(!hist) return 1.;\r
1041 \r
1042  TAxis *ax = hist->GetAxis(0);\r
1043  TAxis *ay = hist->GetAxis(1);\r
1044  Int_t binx = ax->FindBin(zv);\r
1045  Int_t biny = ay->FindBin(mult);\r
1046  Int_t dim[2] = {binx,biny};\r
1047  Double_t fact  = 1.0-hist->GetBinContent(dim);  \r
1048 \r
1049 return fact;\r
1050 }\r
1051 \r
1052 //_____________________________________________________________________________\r
1053 Long64_t AlidNdPtCorrection::Merge(TCollection* list) \r
1054 {\r
1055   // Merge list of objects (needed by PROOF)\r
1056 \r
1057   if (!list)\r
1058   return 0;\r
1059 \r
1060   if (list->IsEmpty())\r
1061   return 1;\r
1062 \r
1063   TIterator* iter = list->MakeIterator();\r
1064   TObject* obj = 0;\r
1065 \r
1066   // collection of generated histograms\r
1067 \r
1068   Int_t count=0;\r
1069   while((obj = iter->Next()) != 0) {\r
1070     AlidNdPtCorrection* entry = dynamic_cast<AlidNdPtCorrection*>(obj);\r
1071     if (entry == 0) continue; \r
1072   \r
1073     fMCEventHist1->Add(entry->fMCEventHist1);\r
1074     fRecEventHist1->Add(entry->fRecEventHist1);\r
1075 \r
1076     fMCAllEventMultHist1->Add(entry->fMCAllEventMultHist1);\r
1077     fMCAllNDEventMultHist1->Add(entry->fMCAllNDEventMultHist1);\r
1078     fMCAllNSDEventMultHist1->Add(entry->fMCAllNSDEventMultHist1);\r
1079     fMCTriggerMultHist1->Add(entry->fMCTriggerMultHist1);\r
1080     fMCEventMultHist1->Add(entry->fMCEventMultHist1);\r
1081 \r
1082     fMCAllPrimTrackMultHist1->Add(entry->fMCAllPrimTrackMultHist1);\r
1083     fMCNDEventAllPrimTrackMultHist1->Add(entry->fMCNDEventAllPrimTrackMultHist1);\r
1084     fMCNSDEventAllPrimTrackMultHist1->Add(entry->fMCNSDEventAllPrimTrackMultHist1);\r
1085     fMCTriggerPrimTrackMultHist1->Add(entry->fMCTriggerPrimTrackMultHist1);\r
1086     fMCEventPrimTrackMultHist1->Add(entry->fMCEventPrimTrackMultHist1);\r
1087 \r
1088     for(Int_t i=0; i<AlidNdPtHelper::kCutSteps; i++) {\r
1089       fRecTrackHist1[i]->Add(entry->fRecTrackHist1[i]);\r
1090     }\r
1091 \r
1092     for(Int_t i=0; i<8; i++) {\r
1093       fCorrRecTrackMultHist1[i]->Add(entry->fCorrRecTrackMultHist1[i]);\r
1094       fPtvsPt[i]->Add(entry->fPtvsPt[i]);\r
1095     }\r
1096 \r
1097     for(Int_t i=0; i<5; i++) {\r
1098       fCorrRecEventHist1[i]->Add(entry->fCorrRecEventHist1[i]);\r
1099       fCorrRecEventHist2[i]->Add(entry->fCorrRecEventHist2[i]);\r
1100     }\r
1101 \r
1102   count++;\r
1103   }\r
1104 \r
1105 return count;\r
1106 }\r
1107  \r
1108 Int_t AlidNdPtCorrection::GetTrueMult(THnSparse *hist, Int_t mult)\r
1109 {\r
1110  if(!hist) return 0;\r
1111  Int_t true_mult = 0;\r
1112 \r
1113  // 0 bins exluded\r
1114  TAxis *ax = hist->GetAxis(0);\r
1115  TAxis *ay = hist->GetAxis(1);\r
1116  ax->SetRange(1,ax->GetNbins());\r
1117  ay->SetRange(1,ay->GetNbins());\r
1118 \r
1119  // measured mult\r
1120  ax->SetRangeUser((Float_t)mult,(Float_t)mult); \r
1121 \r
1122  // get true multiplicity\r
1123  TH1D *h1 = (TH1D *)hist->Projection(1);\r
1124  true_mult = (Int_t)h1->GetMean();\r
1125 \r
1126  return true_mult;\r
1127 }\r
1128 \r
1129 //_____________________________________________________________________________\r
1130 void AlidNdPtCorrection::Analyse() \r
1131 {\r
1132   // Analyse histograms\r
1133   //\r
1134   TH1::AddDirectory(kFALSE);\r
1135   TH1 *h = 0, *hs=0, *hsc=0; \r
1136   TH2 *h2D = 0; \r
1137 \r
1138   TObjArray *aFolderObj = new TObjArray;\r
1139 \r
1140   //\r
1141   // get cuts\r
1142   //\r
1143   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
1144   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
1145   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
1146 \r
1147   if(!evtCuts || !accCuts || !esdTrackCuts) {\r
1148     Error("AlidNdPtCutAnalysis::Analyse()", "cuts not available");\r
1149     return;\r
1150   }\r
1151 \r
1152   //\r
1153   // set min and max values\r
1154   //\r
1155   //Double_t minPt = accCuts->GetMinPt();\r
1156   //Double_t maxPt = accCuts->GetMaxPt();\r
1157   //Double_t minEta = accCuts->GetMinEta();\r
1158   //Double_t maxEta = accCuts->GetMaxEta()-0.00001;\r
1159 \r
1160   //\r
1161   // pt profile\r
1162   //\r
1163   char name[256];\r
1164   for(Int_t i=0; i<8; i++) {\r
1165     h2D = fPtvsPt[i]->Projection(1,0);\r
1166     sprintf(name,"PtvsMeanPt_%d",i);\r
1167     aFolderObj->Add(h2D);\r
1168   }\r
1169 \r
1170   //\r
1171   // event level \r
1172   //\r
1173   h = fCorrRecEventHist1[0]->Projection(1);\r
1174   h->SetName("mult_event_not_corrected");\r
1175   aFolderObj->Add(h);\r
1176 \r
1177   h = fCorrRecEventHist1[1]->Projection(1);\r
1178   h->SetName("mult_event_vertex_corrected");\r
1179   aFolderObj->Add(h);\r
1180 \r
1181   h = fCorrRecEventHist1[2]->Projection(1);\r
1182   h->SetName("mult_trigger_vertex_corrected");\r
1183   aFolderObj->Add(h);\r
1184 \r
1185   h = fCorrRecEventHist1[3]->Projection(1);\r
1186   h->SetName("mult_ND_trigger_vertex_corrected");\r
1187   aFolderObj->Add(h);\r
1188 \r
1189   h = fCorrRecEventHist1[4]->Projection(1);\r
1190   h->SetName("mult_NSD_trigger_vertex_corrected");\r
1191   aFolderObj->Add(h);\r
1192 \r
1193   // empty events\r
1194   h = fCorrRecEventHist2[0]->Projection(1);\r
1195   h->SetName("mult_empty_event_not_corrected");\r
1196   aFolderObj->Add(h);\r
1197 \r
1198   h = fCorrRecEventHist2[1]->Projection(1);\r
1199   h->SetName("mult_empty_event_vertex_corrected");\r
1200   aFolderObj->Add(h);\r
1201 \r
1202   h = fCorrRecEventHist2[2]->Projection(1);\r
1203   h->SetName("mult_empty_trigger_vertex_corrected");\r
1204   aFolderObj->Add(h);\r
1205 \r
1206   h = fCorrRecEventHist2[3]->Projection(1);\r
1207   h->SetName("mult_empty_ND_trigger_vertex_corrected");\r
1208   aFolderObj->Add(h);\r
1209 \r
1210   h = fCorrRecEventHist2[4]->Projection(1);\r
1211   h->SetName("mult_empty_NSD_trigger_vertex_corrected");\r
1212   aFolderObj->Add(h);\r
1213  \r
1214   //\r
1215   // MC available\r
1216   //\r
1217   if(IsUseMCInfo()) {\r
1218 \r
1219   // mc \r
1220   h = fMCAllEventMultHist1->Projection(1);\r
1221   h->SetName("mc_mult_event_acc_prim");\r
1222   aFolderObj->Add(h);\r
1223 \r
1224   h = fMCAllNDEventMultHist1->Projection(1);\r
1225   h->SetName("mc_mult_ND_event_acc_prim");\r
1226   aFolderObj->Add(h);\r
1227 \r
1228   h = fMCAllNSDEventMultHist1->Projection(1);\r
1229   h->SetName("mc_mult_NSD_event_acc_prim");\r
1230   aFolderObj->Add(h);\r
1231 \r
1232   h = fMCTriggerMultHist1->Projection(1);\r
1233   h->SetName("mc_mult_trigger_acc_prim");\r
1234   aFolderObj->Add(h);\r
1235 \r
1236   h = fMCEventMultHist1->Projection(1);\r
1237   h->SetName("mc_mult_trigger_event_acc_prim");\r
1238   aFolderObj->Add(h);\r
1239 \r
1240 \r
1241   //\r
1242   // track level\r
1243   //\r
1244   \r
1245   // limit eta range\r
1246   for(Int_t i=0;i<8;i++) { \r
1247       //fCorrRecTrackMultHist1[i]->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1248       //fCorrRecTrackMultHist1[i]->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1249   }\r
1250   //fMCAllPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1251   //fMCAllPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1252 \r
1253   //fMCNDEventAllPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1254   //fMCNDEventAllPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1255 \r
1256   //fMCNSDEventAllPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1257   //fMCNSDEventAllPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1258 \r
1259   //fMCTriggerPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1260   //fMCTriggerPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1261 \r
1262   //fMCEventPrimTrackMultHist1->GetAxis(0)->SetRangeUser(minPt,maxPt);\r
1263   //fMCEventPrimTrackMultHist1->GetAxis(1)->SetRangeUser(minEta,maxEta);\r
1264 \r
1265   } // end use MC info \r
1266   \r
1267   //\r
1268   h2D = fCorrRecTrackMultHist1[3]->Projection(1,0);\r
1269   h2D->SetName("pt_eta_rec_track_mult_eff_cont_corrected");\r
1270   aFolderObj->Add(h2D);\r
1271 \r
1272   h2D = fCorrRecTrackMultHist1[4]->Projection(1,0);\r
1273   h2D->SetName("pt_eta_rec_event_track_mult_eff_cont_corrected");\r
1274   aFolderObj->Add(h2D);\r
1275 \r
1276   h2D = fCorrRecTrackMultHist1[5]->Projection(1,0);\r
1277   h2D->SetName("pt_eta_rec_trig_event_track_mult_eff_cont_corrected");\r
1278   aFolderObj->Add(h2D);\r
1279 \r
1280   h2D = fCorrRecTrackMultHist1[6]->Projection(1,0);\r
1281   h2D->SetName("pt_eta_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1282   aFolderObj->Add(h2D);\r
1283 \r
1284   h2D = fCorrRecTrackMultHist1[7]->Projection(1,0);\r
1285   h2D->SetName("pt_eta_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1286   aFolderObj->Add(h2D);\r
1287 \r
1288 \r
1289   //\r
1290   h2D = fCorrRecTrackMultHist1[3]->Projection(2,0);\r
1291   h2D->SetName("pt_mult_rec_track_mult_eff_cont_corrected");\r
1292   aFolderObj->Add(h2D);\r
1293 \r
1294   h2D = fCorrRecTrackMultHist1[4]->Projection(2,0);\r
1295   h2D->SetName("pt_mult_rec_event_track_mult_eff_cont_corrected");\r
1296   aFolderObj->Add(h2D);\r
1297 \r
1298   h2D = fCorrRecTrackMultHist1[5]->Projection(2,0);\r
1299   h2D->SetName("pt_mult_rec_trig_event_track_mult_eff_cont_corrected");\r
1300   aFolderObj->Add(h2D);\r
1301 \r
1302   h2D = fCorrRecTrackMultHist1[6]->Projection(2,0);\r
1303   h2D->SetName("pt_mult_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1304   aFolderObj->Add(h2D);\r
1305 \r
1306   h2D = fCorrRecTrackMultHist1[7]->Projection(2,0);\r
1307   h2D->SetName("pt_mult_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1308   aFolderObj->Add(h2D);\r
1309 \r
1310   // pt axis\r
1311 \r
1312   h = fCorrRecTrackMultHist1[0]->Projection(0);\r
1313   h->SetName("pt_rec_track_not_corrected");\r
1314   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1315   hs->SetName("pt_rec_track_not_corrected_s");\r
1316   aFolderObj->Add(hs);\r
1317 \r
1318   //\r
1319   h = fCorrRecTrackMultHist1[1]->Projection(0);\r
1320   h->SetName("pt_rec_track_cont_corrected");\r
1321   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1322   hs->SetName("pt_rec_track_cont_corrected_s");\r
1323   aFolderObj->Add(hs);\r
1324 \r
1325   hsc = (TH1D*)hs->Clone("pt_rec_track_cont_corr_fact");\r
1326   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_not_corrected_s"));\r
1327   aFolderObj->Add(hsc);\r
1328 \r
1329   //\r
1330   h = fCorrRecTrackMultHist1[2]->Projection(0);\r
1331   h->SetName("pt_rec_track_eff_cont_corrected");\r
1332   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1333   hs->SetName("pt_rec_track_eff_cont_corrected_s");\r
1334   aFolderObj->Add(hs);\r
1335 \r
1336   hsc = (TH1D*)hs->Clone("pt_rec_track_eff_corr_fact");\r
1337   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_cont_corrected_s"));\r
1338   aFolderObj->Add(hsc);\r
1339 \r
1340   //\r
1341   h = fCorrRecTrackMultHist1[3]->Projection(0);\r
1342   h->SetName("pt_rec_track_mult_eff_cont_corrected");\r
1343   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1344   hs->SetName("pt_rec_track_mult_eff_cont_corrected_s");\r
1345   aFolderObj->Add(hs);\r
1346 \r
1347   hsc = (TH1D*)hs->Clone("pt_rec_track_mult_corr_fact");\r
1348   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_eff_cont_corrected_s"));\r
1349   aFolderObj->Add(hsc);\r
1350 \r
1351   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1352   hsc = (TH1D*)hs->Clone("pt_rec_track_mult_eff_cont_corr_fact");\r
1353   hsc->Divide((TH1D*)aFolderObj->FindObject("pt_rec_track_not_corrected_s"));\r
1354   aFolderObj->Add(hsc);\r
1355 \r
1356   hsc = (TH1D*)hs->Clone();\r
1357   hsc->SetName("pt_rec_track_mult_eff_cont_corrected_s_norm");\r
1358   hsc->Scale(1./(fCorrRecEventHist1[0]->Projection(1)->Integral()));\r
1359   aFolderObj->Add(hsc);\r
1360 \r
1361   //\r
1362   h = fCorrRecTrackMultHist1[4]->Projection(0);\r
1363   h->SetName("pt_rec_event_track_mult_eff_cont_corrected");\r
1364   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1365   hs->SetName("pt_rec_event_track_mult_eff_cont_corrected_s");\r
1366   aFolderObj->Add(hs);\r
1367 \r
1368   hsc = (TH1D*)hs->Clone();\r
1369   hsc->SetName("pt_rec_event_track_mult_eff_cont_corrected_s_norm");\r
1370   hsc->Scale(1./(fCorrRecEventHist1[1]->Projection(1)->Integral()+fCorrRecEventHist2[1]->Projection(1)->Integral()));\r
1371   aFolderObj->Add(hsc);\r
1372 \r
1373   //\r
1374   h = fCorrRecTrackMultHist1[5]->Projection(0);\r
1375   h->SetName("pt_rec_trig_event_track_mult_eff_cont_corrected");\r
1376   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1377   hs->SetName("pt_rec_trig_event_track_mult_eff_cont_corrected_s");\r
1378   aFolderObj->Add(hs);\r
1379 \r
1380   hsc = (TH1D*)hs->Clone();\r
1381   hsc->SetName("pt_rec_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1382   hsc->Scale(1./(fCorrRecEventHist1[2]->Projection(1)->Integral() + fCorrRecEventHist2[2]->Projection(1)->Integral()));\r
1383   aFolderObj->Add(hsc);\r
1384 \r
1385   //\r
1386   h = fCorrRecTrackMultHist1[6]->Projection(0);\r
1387   h->SetName("pt_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1388   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1389   hs->SetName("pt_rec_ND_trig_event_track_mult_eff_cont_corrected_s");\r
1390   aFolderObj->Add(hs);\r
1391 \r
1392   hsc = (TH1D*)hs->Clone();\r
1393   hsc->SetName("pt_rec_ND_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1394   hsc->Scale(1./(fCorrRecEventHist1[3]->Projection(1)->Integral()+fCorrRecEventHist2[3]->Projection(1)->Integral()));\r
1395   aFolderObj->Add(hsc);\r
1396 \r
1397   //\r
1398   h = fCorrRecTrackMultHist1[7]->Projection(0);\r
1399   h->SetName("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1400   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1401   hs->SetName("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected_s");\r
1402   aFolderObj->Add(hs);\r
1403 \r
1404   hsc = (TH1D*)hs->Clone();\r
1405   hsc->SetName("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1406   hsc->Scale(1./(fCorrRecEventHist1[4]->Projection(1)->Integral()+fCorrRecEventHist2[4]->Projection(1)->Integral()));\r
1407   aFolderObj->Add(hsc);\r
1408 \r
1409   // eta axis\r
1410   h = fCorrRecTrackMultHist1[0]->Projection(1);\r
1411   h->SetName("eta_rec_track_not_corrected");\r
1412   aFolderObj->Add(h);\r
1413   \r
1414   h = fCorrRecTrackMultHist1[1]->Projection(1);\r
1415   h->SetName("eta_rec_track_cont_corrected");\r
1416   aFolderObj->Add(h);\r
1417 \r
1418   h = fCorrRecTrackMultHist1[2]->Projection(1);\r
1419   h->SetName("eta_rec_track_eff_cont_corrected");\r
1420   aFolderObj->Add(h);\r
1421 \r
1422   h = fCorrRecTrackMultHist1[3]->Projection(1);\r
1423   h->SetName("eta_rec_track_mult_eff_cont_corrected");\r
1424   aFolderObj->Add(h);\r
1425 \r
1426   h = fCorrRecTrackMultHist1[4]->Projection(1);\r
1427   h->SetName("eta_rec_event_track_mult_eff_cont_corrected");\r
1428   aFolderObj->Add(h);\r
1429 \r
1430   h = fCorrRecTrackMultHist1[5]->Projection(1);\r
1431   h->SetName("eta_rec_trig_event_track_mult_eff_cont_corrected");\r
1432   aFolderObj->Add(h);\r
1433 \r
1434   h = fCorrRecTrackMultHist1[6]->Projection(1);\r
1435   h->SetName("eta_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1436   aFolderObj->Add(h);\r
1437 \r
1438   h = fCorrRecTrackMultHist1[7]->Projection(1);\r
1439   h->SetName("eta_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1440   aFolderObj->Add(h);\r
1441 \r
1442 \r
1443   //\r
1444   // MC available\r
1445   //\r
1446   if(IsUseMCInfo()) {\r
1447 \r
1448   //\r
1449   h2D = fMCAllPrimTrackMultHist1->Projection(2,0);\r
1450   h2D->SetName("mc_all_pt_mult_acc_prim");\r
1451   aFolderObj->Add(h2D);\r
1452 \r
1453   h2D = fMCAllPrimTrackMultHist1->Projection(1,0);\r
1454   h2D->SetName("mc_all_eta_pt_acc_prim");\r
1455   aFolderObj->Add(h2D);\r
1456 \r
1457   h2D = fMCNDEventAllPrimTrackMultHist1->Projection(2,0);\r
1458   h2D->SetName("mc_ND_all_pt_mult_acc_prim");\r
1459   aFolderObj->Add(h2D);\r
1460 \r
1461   h2D = fMCNDEventAllPrimTrackMultHist1->Projection(1,0);\r
1462   h2D->SetName("mc_ND_all_eta_pt_acc_prim");\r
1463   aFolderObj->Add(h2D);\r
1464 \r
1465   h2D = fMCNSDEventAllPrimTrackMultHist1->Projection(2,0);\r
1466   h2D->SetName("mc_NSD_all_pt_mult_acc_prim");\r
1467   aFolderObj->Add(h2D);\r
1468 \r
1469   h2D = fMCNSDEventAllPrimTrackMultHist1->Projection(1,0);\r
1470   h2D->SetName("mc_NSD_all_eta_pt_acc_prim");\r
1471   aFolderObj->Add(h2D);\r
1472 \r
1473   //\r
1474 \r
1475   h = fMCAllPrimTrackMultHist1->Projection(0);\r
1476   h->SetName("mc_all_pt_acc_prim");\r
1477   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1478   hs->SetName("mc_all_pt_acc_prim_s");\r
1479   aFolderObj->Add(hs);\r
1480 \r
1481   hsc = (TH1D*)hs->Clone();\r
1482   hsc->SetName("mc_all_pt_acc_prim_s_norm");\r
1483   hsc->Scale(1./fMCAllEventMultHist1->Projection(1)->Integral());\r
1484   aFolderObj->Add(hsc);\r
1485 \r
1486   h = fMCNDEventAllPrimTrackMultHist1->Projection(0);\r
1487   h->SetName("mc_ND_all_pt_acc_prim");\r
1488   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1489   hs->SetName("mc_ND_all_pt_acc_prim_s");\r
1490   aFolderObj->Add(hs);\r
1491 \r
1492   hsc = (TH1D*)hs->Clone();\r
1493   hsc->SetName("mc_ND_all_pt_acc_prim_s_norm");\r
1494   hsc->Scale(1./fMCAllNDEventMultHist1->Projection(1)->Integral());\r
1495   aFolderObj->Add(hsc);\r
1496 \r
1497   h = fMCNSDEventAllPrimTrackMultHist1->Projection(0);\r
1498   h->SetName("mc_NSD_all_pt_acc_prim");\r
1499   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1500   hs->SetName("mc_NSD_all_pt_acc_prim_s");\r
1501   aFolderObj->Add(hs);\r
1502 \r
1503   hsc = (TH1D*)hs->Clone();\r
1504   hsc->SetName("mc_NSD_all_pt_acc_prim_s_norm");\r
1505   hsc->Scale(1./fMCAllNSDEventMultHist1->Projection(1)->Integral());\r
1506   aFolderObj->Add(hsc);\r
1507 \r
1508   h = fMCTriggerPrimTrackMultHist1->Projection(0);\r
1509   h->SetName("mc_trigger_all_pt_acc_prim");\r
1510   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1511   hs->SetName("mc_trigger_all_pt_acc_prim_s");\r
1512   aFolderObj->Add(hs);\r
1513 \r
1514   hsc = (TH1D*)hs->Clone();\r
1515   hsc->SetName("mc_trigger_all_pt_acc_prim_s_norm");\r
1516   hsc->Scale(1./fMCTriggerMultHist1->Projection(1)->Integral());\r
1517   aFolderObj->Add(hsc);\r
1518 \r
1519   h = fMCEventPrimTrackMultHist1->Projection(0);\r
1520   h->SetName("mc_all_pt_acc_trig_event_prim");\r
1521   hs = AlidNdPtHelper::ScaleByBinWidth(h);\r
1522   hs->SetName("mc_all_pt_acc_trig_event_prim_s");\r
1523   aFolderObj->Add(hs);\r
1524 \r
1525   hsc = (TH1D*)hs->Clone();\r
1526   hsc->SetName("mc_all_pt_acc_trig_event_prim_s_norm");\r
1527   hsc->Scale(1./fMCEventMultHist1->Projection(1)->Integral());\r
1528   aFolderObj->Add(hsc);\r
1529 \r
1530   //\r
1531 \r
1532   h = fMCAllPrimTrackMultHist1->Projection(1);\r
1533   h->SetName("mc_all_eta_acc_prim");\r
1534   aFolderObj->Add(h);\r
1535 \r
1536   h = fMCNDEventAllPrimTrackMultHist1->Projection(1);\r
1537   h->SetName("mc_ND_all_eta_acc_prim");\r
1538   aFolderObj->Add(h);\r
1539 \r
1540   h = fMCNSDEventAllPrimTrackMultHist1->Projection(1);\r
1541   h->SetName("mc_NSD_all_eta_acc_prim");\r
1542   aFolderObj->Add(h);\r
1543 \r
1544   h = fMCTriggerPrimTrackMultHist1->Projection(1);\r
1545   h->SetName("mc_trigger_all_eta_acc_prim");\r
1546   aFolderObj->Add(h);\r
1547 \r
1548   h = fMCEventPrimTrackMultHist1->Projection(1);\r
1549   h->SetName("mc_all_eta_acc_trig_event_prim");\r
1550   aFolderObj->Add(h);\r
1551 \r
1552   //\r
1553   // calculate ratios (rec / mc)\r
1554   //\r
1555   \r
1556   hs = (TH1*)aFolderObj->FindObject("pt_rec_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1557   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_trig_event_track_mult_eff_cont_corrected"); \r
1558   hsc->Sumw2();\r
1559   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_pt_acc_prim_s_norm"));\r
1560   aFolderObj->Add(hsc);\r
1561 \r
1562   hs = (TH1*)aFolderObj->FindObject("eta_rec_trig_event_track_mult_eff_cont_corrected");\r
1563   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_trig_event_track_mult_eff_cont_corrected"); \r
1564   hsc->Sumw2();\r
1565   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_eta_acc_prim"));\r
1566   aFolderObj->Add(hsc);\r
1567 \r
1568   //\r
1569   hs = (TH1*)aFolderObj->FindObject("pt_rec_ND_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1570   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_ND_trig_event_track_mult_eff_cont_corrected"); \r
1571   hsc->Sumw2();\r
1572   hsc->Divide((TH1*)aFolderObj->FindObject("mc_ND_all_pt_acc_prim_s_norm"));\r
1573   aFolderObj->Add(hsc);\r
1574 \r
1575   hs = (TH1*)aFolderObj->FindObject("eta_rec_ND_trig_event_track_mult_eff_cont_corrected");\r
1576   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_ND_trig_event_track_mult_eff_cont_corrected"); \r
1577   hsc->Sumw2();\r
1578   hsc->Divide((TH1*)aFolderObj->FindObject("mc_ND_all_eta_acc_prim"));\r
1579   aFolderObj->Add(hsc);\r
1580 \r
1581   //\r
1582   hs = (TH1*)aFolderObj->FindObject("pt_rec_NSD_trig_event_track_mult_eff_cont_corrected_s_norm");\r
1583   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_NSD_trig_event_track_mult_eff_cont_corrected"); \r
1584   hsc->Sumw2();\r
1585   hsc->Divide((TH1*)aFolderObj->FindObject("mc_NSD_all_pt_acc_prim_s_norm"));\r
1586   aFolderObj->Add(hsc);\r
1587 \r
1588   hs = (TH1*)aFolderObj->FindObject("eta_rec_NSD_trig_event_track_mult_eff_cont_corrected");\r
1589   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_NSD_trig_event_track_mult_eff_cont_corrected"); \r
1590   hsc->Sumw2();\r
1591   hsc->Divide((TH1*)aFolderObj->FindObject("mc_NSD_all_eta_acc_prim"));\r
1592   aFolderObj->Add(hsc);\r
1593 \r
1594   //\r
1595   hs = (TH1*)aFolderObj->FindObject("pt_rec_event_track_mult_eff_cont_corrected_s_norm");\r
1596   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_event_track_mult_eff_cont_corrected"); \r
1597   hsc->Sumw2();\r
1598   hsc->Divide((TH1*)aFolderObj->FindObject("mc_trigger_all_pt_acc_prim_s_norm"));\r
1599   aFolderObj->Add(hsc);\r
1600 \r
1601   hs = (TH1*)aFolderObj->FindObject("eta_rec_event_track_mult_eff_cont_corrected");\r
1602   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_event_track_mult_eff_cont_corrected"); \r
1603   hsc->Sumw2();\r
1604   hsc->Divide((TH1*)aFolderObj->FindObject("mc_trigger_all_eta_acc_prim"));\r
1605   aFolderObj->Add(hsc);\r
1606 \r
1607   // track level\r
1608   hs = (TH1*)aFolderObj->FindObject("pt_rec_track_mult_eff_cont_corrected_s_norm");\r
1609   hsc = (TH1D*)hs->Clone("ratio_pt_rec_to_mc_track_mult_eff_cont_corrected"); \r
1610   hsc->Sumw2();\r
1611   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_pt_acc_trig_event_prim_s_norm"));\r
1612   aFolderObj->Add(hsc);\r
1613 \r
1614   hs = (TH1*)aFolderObj->FindObject("eta_rec_track_mult_eff_cont_corrected");\r
1615   hsc = (TH1D*)hs->Clone("ratio_eta_rec_to_mc_track_mult_eff_cont_corrected"); \r
1616   hsc->Sumw2();\r
1617   hsc->Divide((TH1*)aFolderObj->FindObject("mc_all_eta_acc_trig_event_prim"));\r
1618   aFolderObj->Add(hsc);\r
1619 \r
1620   } // end MC infor available\r
1621 \r
1622   // export objects to analysis folder\r
1623   fCorrectionFolder = ExportToFolder(aFolderObj);\r
1624 \r
1625   // delete only TObjArray\r
1626   if(aFolderObj) delete aFolderObj;\r
1627 }\r
1628 \r
1629 //_____________________________________________________________________________\r
1630 TFolder* AlidNdPtCorrection::ExportToFolder(TObjArray * array) \r
1631 {\r
1632   // recreate folder avery time and export objects to new one\r
1633   //\r
1634   AlidNdPtCorrection * comp=this;\r
1635   TFolder *folder = comp->GetCorrectionFolder();\r
1636 \r
1637   TString name, title;\r
1638   TFolder *newFolder = 0;\r
1639   Int_t i = 0;\r
1640   Int_t size = array->GetSize();\r
1641 \r
1642   if(folder) { \r
1643      // get name and title from old folder\r
1644      name = folder->GetName();  \r
1645      title = folder->GetTitle();  \r
1646 \r
1647          // delete old one\r
1648      delete folder;\r
1649 \r
1650          // create new one\r
1651      newFolder = CreateFolder(name.Data(),title.Data());\r
1652      newFolder->SetOwner();\r
1653 \r
1654          // add objects to folder\r
1655      while(i < size) {\r
1656            newFolder->Add(array->At(i));\r
1657            i++;\r
1658          }\r
1659   }\r
1660 \r
1661 return newFolder;\r
1662 }\r
1663 \r
1664 //_____________________________________________________________________________\r
1665 TFolder* AlidNdPtCorrection::CreateFolder(TString name,TString title) { \r
1666 // create folder for analysed histograms\r
1667 //\r
1668 TFolder *folder = 0;\r
1669   folder = new TFolder(name.Data(),title.Data());\r
1670 \r
1671   return folder;\r
1672 }\r