change method name
[u/mrichter/AliRoot.git] / PWGJE / AliFragmentationFunctionCorrections.cxx
1 // *************************************************************************
2 // *                                                                       *
3 // *               corrections to Fragmentation Functions                  *
4 // *                                                                       *
5 // *************************************************************************
6
7
8 /**************************************************************************
9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10  *                                                                        *
11  * Author: The ALICE Off-line Project.                                    *
12  * Contributors are mentioned in the code where appropriate.              *
13  *                                                                        *
14  * Permission to use, copy, modify and distribute this software and its   *
15  * documentation strictly for non-commercial purposes is hereby granted   *
16  * without fee, provided that the above copyright notice appears in all   *
17  * copies and that both the copyright notice and this permission notice   *
18  * appear in the supporting documentation. The authors make no claims     *
19  * about the suitability of this software for any purpose. It is          *
20  * provided "as is" without express or implied warranty.                  *
21  **************************************************************************/
22
23 /* $Id: */
24
25 #include "TMath.h"
26 #include "TH2F.h"
27 #include "THnSparse.h"
28 #include "TFile.h"
29 #include "TDirectory.h"
30 #include "AliCFUnfolding.h"
31 #include "AliFragmentationFunctionCorrections.h"
32 #include <iostream> // OB TEST!!!
33 #include <fstream>
34 using std::cout;
35 using std::endl;
36 using std::cerr;
37 using std::fstream;
38
39 ///////////////////////////////////////////////////////////////////////////////
40 //                                                                           //
41 //  correction class for ouput of AliAnalysisTaskFragmentationFunction       //     
42 //                                                                           //
43 //  corrections for: reconstruction efficiency, momentum resolution,         //
44 //                   secondaries, UE / HI background, fluctuations           // 
45 //                   back-correction on jet energy on dN/dxi                 //
46 //                                                                           //
47 //  read MC ouput and write out efficiency histos, response matrices etc.    //
48 //  read measured distributions and apply efficiency, response matrices etc. //
49 //                                                                           //
50 //  contact: o.busch@gsi.de                                                  //
51 //                                                                           //
52 ///////////////////////////////////////////////////////////////////////////////
53
54
55 ClassImp(AliFragmentationFunctionCorrections)
56
57 //________________________________________________________________________
58 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
59    : TObject()
60    ,fDebug(0)
61    ,fNJetPtSlices(0)
62    ,fJetPtSlices(0)   
63    ,fNJets(0)
64    ,fNJetsBgr(0)
65    ,fNHistoBinsSinglePt(0)
66    ,fHistoBinsSinglePt(0)
67    ,fNHistoBinsPt(0)
68    ,fNHistoBinsZ(0)
69    ,fNHistoBinsXi(0)
70    ,fHistoBinsPt(0)
71    ,fHistoBinsZ(0)
72    ,fHistoBinsXi(0)
73    ,fNCorrectionLevels(0)
74    ,fCorrFF(0)
75    ,fNCorrectionLevelsBgr(0)
76    ,fCorrBgr(0)
77    ,fNCorrectionLevelsSinglePt(0)
78    ,fCorrSinglePt(0)
79    ,fh1EffSinglePt(0)
80    ,fh1EffPt(0) 
81    ,fh1EffZ(0) 
82    ,fh1EffXi(0) 
83    ,fh1EffBgrPt(0)
84    ,fh1EffBgrZ(0)
85    ,fh1EffBgrXi(0)
86    ,fh1BbBCorrSinglePt(0) 
87    ,fh1BbBPt(0)
88    ,fh1BbBZ(0)
89    ,fh1BbBXi(0)
90    ,fh1BbBBgrPt(0)
91    ,fh1BbBBgrZ(0)
92    ,fh1BbBBgrXi(0)
93    ,fh1FoldingCorrPt(0)
94    ,fh1FoldingCorrZ(0)
95    ,fh1FoldingCorrXi(0)
96    ,fh1SecCorrPt(0)
97    ,fh1SecCorrZ(0)
98    ,fh1SecCorrXi(0)
99    ,fh1SecCorrBgrPt(0)
100    ,fh1SecCorrBgrZ(0)
101    ,fh1SecCorrBgrXi(0)
102    ,fh1FFTrackPtBackFolded(0)   
103    ,fh1FFZBackFolded(0)
104    ,fh1FFXiBackFolded(0)        
105    ,fh1FFRatioTrackPtFolded(0)
106    ,fh1FFRatioZFolded(0)
107    ,fh1FFRatioXiFolded(0)
108    ,fh1FFRatioTrackPtBackFolded(0)
109    ,fh1FFRatioZBackFolded(0)
110    ,fh1FFRatioXiBackFolded(0)
111    ,fh1SingleTrackPtBackFolded(0)   
112    ,fh1RatioSingleTrackPtFolded(0) 
113    ,fh1RatioSingleTrackPtBackFolded(0)
114    ,fhnResponseSinglePt(0)
115    ,fhnResponsePt(0)
116    ,fhnResponseZ(0) 
117    ,fhnResponseXi(0)       
118    ,fh1FFTrackPtPrior(0)
119    ,fh1FFZPrior(0)
120    ,fh1FFXiPrior(0)
121    ,fh1SecCorrSinglePt(0)
122 {
123    // default constructor
124 }
125
126 //________________________________________________________________________________________________________________________
127 AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const  AliFragmentationFunctionCorrections &copy)
128   : TObject()
129   ,fDebug(copy.fDebug)                        
130   ,fNJetPtSlices(copy.fNJetPtSlices)                 
131   ,fJetPtSlices(copy.fJetPtSlices)                  
132   ,fNJets(copy.fNJets)                        
133   ,fNJetsBgr(copy.fNJetsBgr)
134   ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)  
135   ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt) 
136   ,fNHistoBinsPt(copy.fNHistoBinsPt)                 
137   ,fNHistoBinsZ(copy.fNHistoBinsZ)                  
138   ,fNHistoBinsXi(copy.fNHistoBinsXi)                 
139   ,fHistoBinsPt(copy.fHistoBinsPt)                  
140   ,fHistoBinsZ(copy.fHistoBinsZ)                   
141   ,fHistoBinsXi(copy.fHistoBinsXi)                  
142   ,fNCorrectionLevels(copy.fNCorrectionLevels)            
143   ,fCorrFF(copy.fCorrFF)                       
144   ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)         
145   ,fCorrBgr(copy.fCorrBgr)    
146   ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
147   ,fCorrSinglePt(copy.fCorrSinglePt)   
148   ,fh1EffSinglePt(copy.fh1EffSinglePt)
149   ,fh1EffPt(copy.fh1EffPt)                      
150   ,fh1EffZ(copy.fh1EffZ)                       
151   ,fh1EffXi(copy.fh1EffXi)                      
152   ,fh1EffBgrPt(copy.fh1EffBgrPt)                   
153   ,fh1EffBgrZ(copy.fh1EffBgrZ)                    
154   ,fh1EffBgrXi(copy.fh1EffBgrXi)
155   ,fh1BbBCorrSinglePt(copy.fh1BbBCorrSinglePt) 
156   ,fh1BbBPt(copy.fh1BbBPt)
157   ,fh1BbBZ(copy.fh1BbBZ)
158   ,fh1BbBXi(copy.fh1BbBXi)
159   ,fh1BbBBgrPt(copy.fh1BbBBgrPt)
160   ,fh1BbBBgrZ(copy.fh1BbBBgrZ)
161   ,fh1BbBBgrXi(copy.fh1BbBBgrXi)
162   ,fh1FoldingCorrPt(copy.fh1FoldingCorrPt)
163   ,fh1FoldingCorrZ(copy.fh1FoldingCorrZ)
164   ,fh1FoldingCorrXi(copy.fh1FoldingCorrXi)
165   ,fh1SecCorrPt(copy.fh1SecCorrPt)
166   ,fh1SecCorrZ(copy.fh1SecCorrZ)
167   ,fh1SecCorrXi(copy.fh1SecCorrXi)
168   ,fh1SecCorrBgrPt(copy.fh1SecCorrBgrPt)
169   ,fh1SecCorrBgrZ(copy.fh1SecCorrBgrZ) 
170   ,fh1SecCorrBgrXi(copy.fh1SecCorrBgrXi)
171   ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)        
172   ,fh1FFZBackFolded(copy.fh1FFZBackFolded)              
173   ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)             
174   ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)       
175   ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)             
176   ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)            
177   ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)   
178   ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)         
179   ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)        
180   ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)     
181   ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)    
182   ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
183   ,fhnResponseSinglePt(copy.fhnResponseSinglePt)                 
184   ,fhnResponsePt(copy.fhnResponsePt)                 
185   ,fhnResponseZ(copy.fhnResponseZ)                  
186   ,fhnResponseXi(copy.fhnResponseXi)                 
187   ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)   
188   ,fh1FFZPrior(copy.fh1FFZPrior)
189   ,fh1FFXiPrior(copy.fh1FFXiPrior)
190   ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
191 {
192   // copy constructor
193   
194 }
195
196 // ______________________________________________________________________________________________________________________________
197 AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
198 {
199   // assignment
200   
201   if(this!=&o){
202     TObject::operator=(o);
203     fDebug                          = o.fDebug;                        
204     fNJetPtSlices                   = o.fNJetPtSlices;                 
205     fJetPtSlices                    = o.fJetPtSlices;                  
206     fNJets                          = o.fNJets;                        
207     fNJetsBgr                       = o.fNJetsBgr;                     
208     fNHistoBinsSinglePt             = o.fNHistoBinsSinglePt;  
209     fHistoBinsSinglePt              = o.fHistoBinsSinglePt; 
210     fNHistoBinsPt                   = o.fNHistoBinsPt;                 
211     fNHistoBinsZ                    = o.fNHistoBinsZ;                  
212     fNHistoBinsXi                   = o.fNHistoBinsXi;                 
213     fHistoBinsPt                    = o.fHistoBinsPt;                  
214     fHistoBinsZ                     = o.fHistoBinsZ;                   
215     fHistoBinsXi                    = o.fHistoBinsXi;                  
216     fNCorrectionLevels              = o.fNCorrectionLevels;            
217     fCorrFF                         = o.fCorrFF;                       
218     fNCorrectionLevelsBgr           = o.fNCorrectionLevelsBgr;         
219     fCorrBgr                        = o.fCorrBgr;                      
220     fNCorrectionLevelsSinglePt      = o.fNCorrectionLevelsSinglePt;
221     fCorrSinglePt                   = o.fCorrSinglePt;
222     fh1EffSinglePt                  = o.fh1EffSinglePt;
223     fh1EffPt                        = o.fh1EffPt;                      
224     fh1EffZ                         = o.fh1EffZ;                       
225     fh1EffXi                        = o.fh1EffXi;                      
226     fh1EffBgrPt                     = o.fh1EffBgrPt;                   
227     fh1EffBgrZ                      = o.fh1EffBgrZ;                    
228     fh1EffBgrXi                     = o.fh1EffBgrXi;                   
229     fh1BbBCorrSinglePt              = o.fh1BbBCorrSinglePt; 
230     fh1BbBPt                        = o.fh1BbBPt;
231     fh1BbBZ                         = o.fh1BbBZ;
232     fh1BbBXi                        = o.fh1BbBXi;
233     fh1BbBBgrPt                     = o.fh1BbBBgrPt;
234     fh1BbBBgrZ                      = o.fh1BbBBgrZ;
235     fh1BbBBgrXi                     = o.fh1BbBBgrXi;
236     fh1FoldingCorrPt                = o.fh1FoldingCorrPt;
237     fh1FoldingCorrZ                 = o.fh1FoldingCorrZ;
238     fh1FoldingCorrXi                = o.fh1FoldingCorrXi;
239     fh1SecCorrPt                    = o.fh1SecCorrPt;
240     fh1SecCorrZ                     = o.fh1SecCorrZ;
241     fh1SecCorrXi                    = o.fh1SecCorrXi;
242     fh1SecCorrBgrPt                 = o.fh1SecCorrBgrPt;
243     fh1SecCorrBgrZ                  = o.fh1SecCorrBgrZ;
244     fh1SecCorrBgrXi                 = o.fh1SecCorrBgrXi;
245     fh1FFTrackPtBackFolded          = o.fh1FFTrackPtBackFolded;        
246     fh1FFZBackFolded                = o.fh1FFZBackFolded;              
247     fh1FFXiBackFolded               = o.fh1FFXiBackFolded;             
248     fh1FFRatioTrackPtFolded         = o.fh1FFRatioTrackPtFolded;       
249     fh1FFRatioZFolded               = o.fh1FFRatioZFolded;             
250     fh1FFRatioXiFolded              = o.fh1FFRatioXiFolded;            
251     fh1FFRatioTrackPtBackFolded     = o.fh1FFRatioTrackPtBackFolded;   
252     fh1FFRatioZBackFolded           = o.fh1FFRatioZBackFolded;         
253     fh1FFRatioXiBackFolded          = o.fh1FFRatioXiBackFolded;        
254     fh1SingleTrackPtBackFolded      = o.fh1SingleTrackPtBackFolded;     
255     fh1RatioSingleTrackPtFolded     = o.fh1RatioSingleTrackPtFolded;    
256     fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded;
257     fhnResponseSinglePt             = o.fhnResponseSinglePt;                 
258     fhnResponsePt                   = o.fhnResponsePt;                 
259     fhnResponseZ                    = o.fhnResponseZ;                  
260     fhnResponseXi                   = o.fhnResponseXi;                 
261     fh1FFTrackPtPrior               = o.fh1FFTrackPtPrior;
262     fh1FFZPrior                     = o.fh1FFZPrior;
263     fh1FFXiPrior                    = o.fh1FFXiPrior;
264     fh1SecCorrSinglePt              = o.fh1SecCorrSinglePt;
265   }
266   
267   return *this;
268 }
269
270 //_________________________________________________________________________
271 AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
272 {
273   // destructor  
274
275   if(fJetPtSlices) delete fJetPtSlices;
276   if(fNJets)       delete fNJets;
277   if(fNJetsBgr)    delete fNJetsBgr;
278
279   DeleteHistoArray(fh1EffPt);
280   DeleteHistoArray(fh1EffXi);
281   DeleteHistoArray(fh1EffZ );
282
283   DeleteHistoArray(fh1EffBgrPt);
284   DeleteHistoArray(fh1EffBgrXi);
285   DeleteHistoArray(fh1EffBgrZ);
286
287   // BbB
288   DeleteHistoArray(fh1BbBPt);
289   DeleteHistoArray(fh1BbBXi);
290   DeleteHistoArray(fh1BbBZ);
291
292   DeleteHistoArray(fh1BbBBgrPt);
293   DeleteHistoArray(fh1BbBBgrXi);
294   DeleteHistoArray(fh1BbBBgrZ);
295
296   DeleteHistoArray(fh1FoldingCorrPt);
297   DeleteHistoArray(fh1FoldingCorrXi);
298   DeleteHistoArray(fh1FoldingCorrZ);
299
300   // sec corr
301   DeleteHistoArray(fh1SecCorrPt);
302   DeleteHistoArray(fh1SecCorrXi);
303   DeleteHistoArray(fh1SecCorrZ);
304
305   DeleteHistoArray(fh1SecCorrBgrPt);
306   DeleteHistoArray(fh1SecCorrBgrXi);
307   DeleteHistoArray(fh1SecCorrBgrZ);
308
309   // unfolding
310
311   DeleteHistoArray(fh1FFTrackPtBackFolded);
312   DeleteHistoArray(fh1FFZBackFolded);
313   DeleteHistoArray(fh1FFXiBackFolded);
314
315   DeleteHistoArray(fh1FFRatioTrackPtFolded);
316   DeleteHistoArray(fh1FFRatioZFolded);
317   DeleteHistoArray(fh1FFRatioXiFolded);
318
319   DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
320   DeleteHistoArray(fh1FFRatioZBackFolded);
321   DeleteHistoArray(fh1FFRatioXiBackFolded);
322
323   DeleteTHnSparseArray(fhnResponsePt);
324   DeleteTHnSparseArray(fhnResponseZ);
325   DeleteTHnSparseArray(fhnResponseXi);
326
327   DeleteHistoArray(fh1FFTrackPtPrior);
328   DeleteHistoArray(fh1FFZPrior);
329   DeleteHistoArray(fh1FFXiPrior);
330
331   // clean up corrected FF 
332   
333   for(Int_t c=0; c<fNCorrectionLevels; c++) delete  fCorrFF[c];
334   delete[] fCorrFF;
335
336   // clean up bgr
337
338   for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
339   delete[] fCorrBgr;
340
341   // clean up inclusive pt 
342   for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
343   delete[] fCorrSinglePt;
344
345   delete[] fNHistoBinsPt;
346   delete[] fNHistoBinsZ;
347   delete[] fNHistoBinsXi;
348
349   // clean up histo bins
350
351   if(fHistoBinsSinglePt) delete fHistoBinsSinglePt; 
352
353   for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsPt[i];
354   for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsZ[i];
355   for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsXi[i];
356
357   delete[] fHistoBinsPt;
358   delete[] fHistoBinsZ;
359   delete[] fHistoBinsXi;
360 }
361
362 //_________________________________________________________________________________
363 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
364   : TObject()
365   ,fArraySize(0)
366   ,fh1CorrFFTrackPt(0)
367   ,fh1CorrFFZ(0)
368   ,fh1CorrFFXi(0)
369   ,fCorrLabel(0)
370 {
371   // default constructor
372   
373 }
374
375 //__________________________________________________________________________________________________________________
376 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& copy)
377   : TObject()
378   ,fArraySize(copy.fArraySize)
379   ,fh1CorrFFTrackPt(0)
380   ,fh1CorrFFZ(0)
381   ,fh1CorrFFXi(0)
382   ,fCorrLabel(copy.fCorrLabel)
383 {
384   // copy constructor
385
386   fh1CorrFFTrackPt = new TH1F*[copy.fArraySize];
387   fh1CorrFFZ       = new TH1F*[copy.fArraySize];
388   fh1CorrFFXi      = new TH1F*[copy.fArraySize];
389
390   for(Int_t i=0; i<copy.fArraySize; i++){
391     fh1CorrFFTrackPt[i] = new TH1F(*(copy.fh1CorrFFTrackPt[i]));
392     fh1CorrFFZ[i]       = new TH1F(*(copy.fh1CorrFFZ[i]));
393     fh1CorrFFXi[i]      = new TH1F(*(copy.fh1CorrFFXi[i]));
394   }
395 }
396
397 //_______________________________________________________________________________________________________________________________________________________________
398 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
399 {
400   // assignment
401   
402   if(this!=&o){
403     TObject::operator=(o);
404     Int_t fArraySize_new = o.fArraySize;
405     fCorrLabel       = o.fCorrLabel;
406  
407     TH1F** fh1CorrFFTrackPt_new = new TH1F*[fArraySize_new];
408     TH1F** fh1CorrFFZ_new       = new TH1F*[fArraySize_new];
409     TH1F** fh1CorrFFXi_new      = new TH1F*[fArraySize_new];
410
411     for(Int_t i=0; i<fArraySize_new; i++){
412       fh1CorrFFTrackPt_new[i] = new TH1F(*(o.fh1CorrFFTrackPt[i]));
413       fh1CorrFFZ_new[i]       = new TH1F(*(o.fh1CorrFFZ[i]));
414       fh1CorrFFXi_new[i]      = new TH1F(*(o.fh1CorrFFXi[i]));
415     }
416     
417     // ---
418
419     if(fArraySize){
420       for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
421       for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
422       for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];          
423     }
424
425     if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
426     if(fh1CorrFFZ)       delete[] fh1CorrFFZ;
427     if(fh1CorrFFXi)      delete[] fh1CorrFFXi;
428         
429     // ---
430
431     fArraySize = fArraySize_new;
432
433     fh1CorrFFTrackPt = fh1CorrFFTrackPt_new;
434     fh1CorrFFZ       = fh1CorrFFZ_new;
435     fh1CorrFFXi      = fh1CorrFFXi_new;
436     
437     for(Int_t i=0; i<fArraySize; i++){
438       fh1CorrFFTrackPt[i] = fh1CorrFFTrackPt_new[i];
439       fh1CorrFFZ[i]       = fh1CorrFFZ_new[i];
440       fh1CorrFFXi[i]      = fh1CorrFFXi_new[i];
441     }
442   }
443   
444   return *this;
445 }
446
447 //__________________________________________________________________________________
448 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
449 {
450   // destructor 
451
452   if(fArraySize){
453     for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
454     for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
455     for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];          
456   }
457
458   if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
459   if(fh1CorrFFZ)       delete[] fh1CorrFFZ;
460   if(fh1CorrFFXi)      delete[] fh1CorrFFXi;
461
462 }
463
464 //___________________________________________________________________________________________________________________
465 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
466   : TObject()
467   ,fArraySize(0)
468   ,fh1CorrFFTrackPt(0)
469   ,fh1CorrFFZ(0)
470   ,fh1CorrFFXi(0)
471   ,fCorrLabel(label)
472 {
473   // constructor
474
475   fArraySize = arraySize;
476   fh1CorrFFTrackPt = new TH1F*[arraySize];
477   fh1CorrFFZ       = new TH1F*[arraySize];
478   fh1CorrFFXi      = new TH1F*[arraySize];
479
480   for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
481   for(int i=0; i<arraySize; i++) fh1CorrFFZ[i]       = new TH1F;
482   for(int i=0; i<arraySize; i++) fh1CorrFFXi[i]      = new TH1F;
483 }
484
485 //_______________________________________________________________________________________________________________________________
486 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
487 {
488   // place histo in array, add corrLabel to histo name
489
490   if(slice>=fArraySize){
491     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
492     return;
493   }
494
495   if(histPt){
496     TString name = histPt->GetName();
497     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
498     histPt->SetName(name);
499     
500     if(!(histPt->GetSumw2())) histPt->Sumw2(); 
501
502     new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
503   }
504   
505   if(histZ){
506     TString name = histZ->GetName();
507     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
508     histZ->SetName(name);
509
510     if(!(histZ->GetSumw2())) histZ->Sumw2(); 
511
512     new(fh1CorrFFZ[slice]) TH1F(*histZ);
513   }
514
515   if(histXi){
516     TString name = histXi->GetName();
517     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
518     histXi->SetName(name);
519
520     if(!(histXi->GetSumw2())) histXi->Sumw2(); 
521
522     new(fh1CorrFFXi[slice]) TH1F(*histXi);
523   }
524 }
525
526 //___________________________________________________________________________________________________________________________________
527 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
528 {
529   // replace histo in array 
530
531   if(slice>=fArraySize){
532     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
533     return;
534   }
535
536   if(histPt){
537     if(!(histPt->GetSumw2())) histPt->Sumw2(); 
538
539     TString name = histPt->GetName();
540     histPt->SetName(name);
541
542     delete fh1CorrFFTrackPt[slice];
543     fh1CorrFFTrackPt[slice] =  new TH1F;
544     new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
545   }
546   
547   if(histZ){
548     if(!(histZ->GetSumw2())) histZ->Sumw2(); 
549
550     TString name = histZ->GetName();
551     histZ->SetName(name);
552
553     delete fh1CorrFFZ[slice];
554     fh1CorrFFZ[slice] = new TH1F;
555     new(fh1CorrFFZ[slice]) TH1F(*histZ);
556   }
557
558   if(histXi){
559     if(!(histXi->GetSumw2())) histXi->Sumw2(); 
560     
561     TString name = histXi->GetName();
562     histXi->SetName(name);
563
564     delete fh1CorrFFXi[slice];
565     fh1CorrFFXi[slice] = new TH1F;
566     new(fh1CorrFFXi[slice]) TH1F(*histXi);
567   }
568 }
569
570 // ___________________________________________________________________________________________
571 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
572
573   // return pt histo 
574   
575   if(slice>=fArraySize){
576     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
577     return 0;
578   }
579
580   return fh1CorrFFTrackPt[slice]; 
581
582 }
583
584 // ______________________________________________________________________________________
585 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
586
587   // return z histo 
588   
589   if(slice>=fArraySize){
590     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
591     return 0;
592   }
593
594   return fh1CorrFFZ[slice]; 
595 }
596
597 // ________________________________________________________________________________________
598 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
599
600   // return xi histo
601
602   if(slice>=fArraySize){
603     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
604     return 0;
605   }
606
607   return fh1CorrFFXi[slice]; 
608 }
609
610 // __________________________________________________________________________
611 void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
612 {
613   // delete array of TH1 
614  
615   if(hist){
616     for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
617     delete[] hist;
618   }
619 }
620
621 // ____________________________________________________________________________________
622 void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
623 {
624   // delete array of THnSparse 
625  
626   if(hist){
627     for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
628     delete[] hist;
629   }
630 }
631
632 // _________________________________________________________
633 TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
634 {
635   // book array of TH1
636
637   if(!fNJetPtSlices){
638     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
639     return 0;
640   }
641
642   TH1F** hist = new TH1F*[fNJetPtSlices];
643   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
644   
645   return hist;
646 }
647
648 //__________________________________________________________________
649 THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
650 {
651   // book array of THnSparse
652
653   if(!fNJetPtSlices){
654     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
655     return 0;
656   }
657
658   THnSparse** hist = new THnSparse*[fNJetPtSlices];
659   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF(); 
660   
661   return hist;
662 }
663
664 //_____________________________________________________________________________
665 void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
666 {
667   // increase corr level 
668
669   if(!fNJetPtSlices){
670     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
671     return;
672   }
673
674   if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
675     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
676     return;
677   }
678
679   fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
680   fNCorrectionLevels++;
681 }
682
683
684 //________________________________________________________________________________
685 void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
686 {
687   // increase corr level for bgr FF
688
689   if(!fNJetPtSlices){
690     if(fDebug>0)  Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
691     return;
692   }
693   
694   if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
695     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
696     return;
697   }
698
699   fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
700   fNCorrectionLevelsBgr++;
701 }
702
703 //_____________________________________________________________________________________
704 void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
705 {
706   // increase corr level single pt spec
707
708   Int_t nJetPtSlicesSingle = 1; // pro forma 
709
710   if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
711     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
712     return;
713   }
714
715   fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
716   fNCorrectionLevelsSinglePt++;
717 }
718
719 //_____________________________________________________________________________________________
720 void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
721
722   // set jet pt slices
723   // once slices are known, can book all other histos 
724
725   fNJetPtSlices = nJetPtSlices;
726
727   const Int_t size = nJetPtSlices+1;
728   fJetPtSlices  = new TArrayF(size,bins);
729  
730   // nJets array
731
732   fNJets = new TArrayF(size);
733   fNJets->Reset(0);
734
735   fNJetsBgr = new TArrayF(size);
736   fNJetsBgr->Reset(0);
737
738   // histos 
739
740   fNCorrectionLevels = 0; 
741   fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
742   AddCorrectionLevel(); // first 'correction' level = raw FF
743
744   fNCorrectionLevelsBgr = 0; 
745   fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
746   AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
747
748   // eff histos
749
750   fh1EffPt = BookHistoArray();
751   fh1EffXi = BookHistoArray();
752   fh1EffZ  = BookHistoArray();
753
754   fh1EffBgrPt = BookHistoArray();
755   fh1EffBgrXi = BookHistoArray();
756   fh1EffBgrZ  = BookHistoArray();
757
758
759   // BbB
760
761   fh1BbBPt = BookHistoArray();
762   fh1BbBXi = BookHistoArray();
763   fh1BbBZ  = BookHistoArray();
764
765   fh1BbBBgrPt = BookHistoArray();
766   fh1BbBBgrXi = BookHistoArray();
767   fh1BbBBgrZ  = BookHistoArray();
768
769   fh1FoldingCorrPt = BookHistoArray();
770   fh1FoldingCorrXi = BookHistoArray();
771   fh1FoldingCorrZ  = BookHistoArray();
772
773   // SecCorr
774
775   fh1SecCorrPt = BookHistoArray();
776   fh1SecCorrXi = BookHistoArray();
777   fh1SecCorrZ  = BookHistoArray();
778
779   fh1SecCorrBgrPt = BookHistoArray();
780   fh1SecCorrBgrXi = BookHistoArray();
781   fh1SecCorrBgrZ  = BookHistoArray();
782
783   // unfolding
784
785   fh1FFTrackPtBackFolded = BookHistoArray();
786   fh1FFXiBackFolded      = BookHistoArray();
787   fh1FFZBackFolded       = BookHistoArray();
788
789   fh1FFRatioTrackPtFolded = BookHistoArray();
790   fh1FFRatioXiFolded      = BookHistoArray();
791   fh1FFRatioZFolded       = BookHistoArray();
792   
793   fh1FFRatioTrackPtBackFolded = BookHistoArray();
794   fh1FFRatioXiBackFolded      = BookHistoArray();
795   fh1FFRatioZBackFolded       = BookHistoArray();
796
797   //
798
799   fhnResponsePt      = BookTHnSparseArray(); 
800   fhnResponseZ       = BookTHnSparseArray();
801   fhnResponseXi      = BookTHnSparseArray();
802
803   fh1FFTrackPtPrior =  BookHistoArray();
804   fh1FFZPrior       =  BookHistoArray();
805   fh1FFXiPrior      =  BookHistoArray();
806
807   // histos bins 
808
809   fNHistoBinsPt = new Int_t[fNJetPtSlices];
810   fNHistoBinsXi = new Int_t[fNJetPtSlices];
811   fNHistoBinsZ  = new Int_t[fNJetPtSlices];
812
813   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
814   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
815   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i]  = 0;
816   
817   fHistoBinsPt = new TArrayD*[fNJetPtSlices];
818   fHistoBinsXi = new TArrayD*[fNJetPtSlices];
819   fHistoBinsZ  = new TArrayD*[fNJetPtSlices];
820
821   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
822   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
823   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i]  = new TArrayD(0);
824 }
825
826 //_____________________________________________________________________________________________________________________________________
827 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
828
829   // set histo bins for jet pt slice
830   // if binning undefined for any slice, original binning will be used
831
832   if(!fNJetPtSlices){
833     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
834     return;
835   }
836   
837   if(jetPtSlice>=fNJetPtSlices){
838     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
839     return;
840   }
841   
842   if(type == kFlagPt){
843     fNHistoBinsPt[jetPtSlice] = sizeBins-1;
844     fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins); 
845   }
846   else if(type==kFlagZ){
847     fNHistoBinsZ[jetPtSlice] = sizeBins-1;
848     fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins); 
849   }
850
851   else if(type==kFlagXi){
852     fNHistoBinsXi[jetPtSlice] = sizeBins-1;
853     fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins); 
854   }
855 }
856
857 //__________________________________________________________________________________________________________________________________________________________________
858 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
859
860   // set histo bins for jet pt slice
861   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
862   // array size of binsLimits: nBinsLimits 
863   // array size of binsWidth: nBinsLimits-1 
864   // binsLimits have to be in increasing order
865   // if binning undefined for any slice, original binning will be used
866
867   if(!fNJetPtSlices){
868     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
869     return;
870   }
871
872   if(jetPtSlice>=fNJetPtSlices){
873     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
874     return;
875   }
876
877
878   Double_t binLimitMin = binsLimits[0];
879   Double_t binLimitMax = binsLimits[nBinsLimits-1];
880
881   Double_t binLimit = binLimitMin; // start value 
882   
883   Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
884   TArrayD binsArray(sizeUpperLim);
885   Int_t nBins = 0; 
886   binsArray.SetAt(binLimitMin,nBins++);
887
888   while(binLimit<0.999999*binLimitMax && nBins<sizeUpperLim){ // 0.999 numerical safety factor
889
890     Int_t currentSlice = -1;
891     for(Int_t i=0; i<nBinsLimits; i++){
892       if(binLimit >= 0.999999*binsLimits[i]) currentSlice = i;  
893     }
894     
895     Double_t currentBinWidth = binsWidth[currentSlice];
896     binLimit += currentBinWidth;
897
898     // if(type == kFlagZ) std::cout<<" SetHistoBins: nBins "<<nBins<<" binLimit "<<binLimit<<" binLimitMax "<<binLimitMax
899     //                          <<" currentSlice "<<currentSlice<<std::endl;
900
901     binsArray.SetAt(binLimit,nBins++);
902   }
903   
904   Double_t* bins = binsArray.GetArray();
905
906   SetHistoBins(jetPtSlice,nBins,bins,type); 
907 }
908
909 //_____________________________________________________________________________________________________________________________________
910 TArrayD* AliFragmentationFunctionCorrections::GetHistoBins(const Int_t& jetPtSlice, const Int_t& type)
911
912   // set histo bins for jet pt slice
913   // if binning undefined for any slice, original binning will be used
914
915   if(!fNJetPtSlices){
916     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
917     return 0;
918   }
919   
920   if(jetPtSlice>=fNJetPtSlices){
921     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
922     return 0;
923   }
924   
925   if(type == kFlagPt){
926     return fHistoBinsPt[jetPtSlice]; 
927   }
928   else if(type==kFlagZ){
929     return fHistoBinsZ[jetPtSlice]; 
930   }
931   
932   else if(type==kFlagXi){
933     return fHistoBinsXi[jetPtSlice]; 
934   }
935   else{ 
936     Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
937     return 0;
938   }
939 }
940
941 //__________________________________________________________________________________________________
942 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
943
944   // set histo bins for inclusive pt spectra
945   // if binning undefined, original binning will be used
946
947   fNHistoBinsSinglePt = sizeBins-1;
948   
949   fHistoBinsSinglePt = new TArrayD(sizeBins);
950   fHistoBinsSinglePt->Set(sizeBins,bins);  
951 }
952
953 //__________________________________________________________________________________________________________________________________________________________________
954 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
955
956   // set histo bins for inclusive pt spectra 
957   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
958   // array size of binsLimits: nBinsLimits 
959   // array size of binsWidth: nBinsLimits-1 
960   // binsLimits have to be in increasing order
961   // if binning undefined for any slice, original binning will be used
962
963
964   Double_t binLimitMin = binsLimits[0];
965   Double_t binLimitMax = binsLimits[nBinsLimits-1];
966
967   Double_t binLimit = binLimitMin; // start value 
968   
969   Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
970   TArrayD binsArray(sizeUpperLim);
971   Int_t nBins = 0; 
972   binsArray.SetAt(binLimitMin,nBins++);
973
974   while(binLimit<binLimitMax && nBins<sizeUpperLim){
975
976     Int_t currentSlice = -1;
977     for(Int_t i=0; i<nBinsLimits; i++){
978       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
979     }
980     
981     Double_t currentBinWidth = binsWidth[currentSlice];
982     binLimit += currentBinWidth;
983
984     binsArray.SetAt(binLimit,nBins++);
985   }
986   
987   Double_t* bins = binsArray.GetArray();
988   
989   SetHistoBinsSinglePt(nBins,bins); 
990 }
991
992 //____________________________________________________________________________________
993 void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
994 {
995   // FF normalization: divide by bin width and normalize to entries/jets
996   // should also work for TH2, but to be tested !!!
997
998   if(nJets == 0){ // nothing to do
999     if(fDebug>0)  Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
1000     return; 
1001   }
1002   
1003   Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
1004
1005   for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
1006
1007     Double_t binWidth = hist->GetBinWidth(bin);
1008     Double_t binCont  = hist->GetBinContent(bin);
1009     Double_t binErr   = hist->GetBinError(bin);
1010     
1011     binCont /= binWidth;
1012     binErr  /= binWidth;
1013
1014     hist->SetBinContent(bin,binCont);
1015     hist->SetBinError(bin,binErr);
1016   }
1017
1018   hist->Scale(1/nJets);
1019 }
1020
1021 //_____________________________________________________
1022 void AliFragmentationFunctionCorrections::NormalizeFF()
1023 {
1024   // normalize FF
1025
1026   if(fNCorrectionLevels>1){
1027     Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
1028   }
1029
1030   for(Int_t i=0; i<fNJetPtSlices; i++){
1031     
1032     if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
1033
1034     NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
1035     NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
1036     NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
1037   } 
1038 }
1039
1040 //______________________________________________________
1041 void AliFragmentationFunctionCorrections::NormalizeBgr()
1042 {
1043   // normalize bgr/UE FF
1044
1045   if(fNCorrectionLevelsBgr>1){
1046     Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
1047   }
1048
1049   for(Int_t i=0; i<fNJetPtSlices; i++){
1050     NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
1051     NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
1052     NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
1053   } 
1054
1055 }
1056
1057 //__________________________________________________________________________________________________
1058 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
1059
1060   // read raw FF - standard dir/list name
1061      
1062   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
1063   TString strlist = "fracfunc_" + strID;
1064     
1065   ReadRawFF(strfile,strdir,strlist,strFFID);
1066 }
1067
1068 //____________________________________________________________________________________________________________________
1069 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
1070 {
1071   // get raw FF from input file, project in jet pt slice
1072   // normalization done separately 
1073
1074   TFile f(strfile,"READ");
1075
1076   if(!f.IsOpen()){
1077     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1078     return;
1079   }
1080
1081   if(fDebug>0) Printf("%s:%d -- read FF from file %s, dir %s ",(char*)__FILE__,__LINE__,strfile.Data(), strdir.Data());
1082
1083   gDirectory->cd(strdir);
1084
1085   TList* list = 0;
1086   
1087   if(strlist && strlist.Length()){
1088     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1089       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1090       return;
1091     }
1092   }
1093
1094   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
1095   TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1096   TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
1097   TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
1098
1099   TH1F* fh1FFJetPt   = 0;
1100   TH2F* fh2FFTrackPt = 0;
1101   TH2F* fh2FFZ       = 0;
1102   TH2F* fh2FFXi      = 0;
1103   
1104   if(list){
1105     fh1FFJetPt   = (TH1F*) list->FindObject(hnameJetPt);
1106     fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
1107     fh2FFZ       = (TH2F*) list->FindObject(hnameZ);  
1108     fh2FFXi      = (TH2F*) list->FindObject(hnameXi); 
1109   }
1110   else{
1111     fh1FFJetPt   = (TH1F*) gDirectory->Get(hnameJetPt);
1112     fh2FFTrackPt = (TH2F*) gDirectory->Get(hnameTrackPt);
1113     fh2FFZ       = (TH2F*) gDirectory->Get(hnameZ);  
1114     fh2FFXi      = (TH2F*) gDirectory->Get(hnameXi);
1115   }
1116
1117
1118   if(!fh1FFJetPt)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());   return; }
1119   if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
1120   if(!fh2FFZ)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data());       return; }
1121   if(!fh2FFXi)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data());      return; }
1122
1123   fh1FFJetPt->SetDirectory(0);
1124   fh2FFTrackPt->SetDirectory(0);
1125   fh2FFZ->SetDirectory(0);  
1126   fh2FFXi->SetDirectory(0); 
1127
1128   f.Close();  
1129
1130
1131   // nJets per bin
1132
1133   for(Int_t i=0; i<fNJetPtSlices; i++){
1134     
1135     Float_t jetPtLoLim = fJetPtSlices->At(i);
1136     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1137     
1138     Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
1139     Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
1140                                      
1141     Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
1142     
1143     fNJets->SetAt(nJetsBin,i); 
1144
1145     if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
1146   }
1147   
1148   // projections: FF 
1149   
1150   for(Int_t i=0; i<fNJetPtSlices; i++){
1151     
1152     Float_t jetPtLoLim = fJetPtSlices->At(i);
1153     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1154
1155     Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
1156     Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
1157
1158     if(binUp > fh2FFTrackPt->GetNbinsX()){
1159       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
1160       return; 
1161     }
1162     
1163     TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1164     TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1165     TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1166     
1167     // appendix 'unbinned' to avoid histos with same name after rebinning
1168     TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
1169     TH1F* projZ  = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
1170     TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
1171     
1172     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
1173     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
1174     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
1175
1176     projPt->SetNameTitle(strNameFFPt,"");
1177     projZ->SetNameTitle(strNameFFZ,"");
1178     projXi->SetNameTitle(strNameFFXi,"");
1179
1180     // raw FF = corr level 0
1181     fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
1182   }
1183
1184
1185   delete fh1FFJetPt;
1186   delete fh2FFTrackPt;
1187   delete fh2FFZ;  
1188   delete fh2FFXi; 
1189 }
1190
1191 //_____________________________________________________________________________________________________________________
1192 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1193
1194   // read raw FF - standard dir/list name
1195      
1196   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
1197   TString strlist = "fracfunc_" + strID;
1198     
1199   ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1200 }
1201
1202 //_______________________________________________________________________________________________________________________________________
1203 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1204 {
1205   // get raw FF from input file, project in jet pt slice
1206   // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1207   // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1208   // normalization done separately 
1209
1210   TString strID = strBgrID + strFFID;
1211   
1212   TFile f(strfile,"READ");
1213
1214   if(!f.IsOpen()){
1215     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1216     return;
1217   }
1218
1219   if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s, dir %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data(),strdir.Data());
1220
1221   gDirectory->cd(strdir);
1222
1223   TList* list = 0;
1224   
1225   if(strlist && strlist.Length()){
1226     if(!(list = (TList*) gDirectory->Get(strlist))){
1227       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1228       return;
1229     }
1230   }
1231
1232   TString hnameNJets = "fh1nRecJetsCuts"; 
1233   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1234   TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1235   TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1236   TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1237
1238   TH1F* fh1NJets;
1239   TH1F* fh1FFJetPtBgr;  
1240   TH2F* fh2FFTrackPtBgr;
1241   TH2F* fh2FFZBgr;
1242   TH2F* fh2FFXiBgr;     
1243   
1244   if(list){
1245     fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1246     fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
1247     fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1248     fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);  
1249     fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi); 
1250   }
1251   else{
1252     fh1NJets        = (TH1F*) gDirectory->Get(hnameNJets); // needed for normalization of bgr out of 2 jets
1253     fh1FFJetPtBgr   = (TH1F*) gDirectory->Get(hnameJetPt);
1254     fh2FFTrackPtBgr = (TH2F*) gDirectory->Get(hnameBgrTrackPt);
1255     fh2FFZBgr       = (TH2F*) gDirectory->Get(hnameBgrZ);  
1256     fh2FFXiBgr      = (TH2F*) gDirectory->Get(hnameBgrXi);  
1257   }
1258
1259
1260   if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
1261   if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());              }
1262   if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1263   if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
1264   if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }
1265
1266   fh1FFJetPtBgr->SetDirectory(0);
1267   fh2FFTrackPtBgr->SetDirectory(0);
1268   fh2FFZBgr->SetDirectory(0);  
1269   fh2FFXiBgr->SetDirectory(0); 
1270   if(fh1NJets) fh1NJets->SetDirectory(0);
1271
1272   f.Close();  
1273
1274   // nJets per bin
1275
1276   for(Int_t i=0; i<fNJetPtSlices; i++){
1277     
1278     Float_t jetPtLoLim = fJetPtSlices->At(i);
1279     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1280     
1281     Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1282     Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1283                                      
1284     Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1285     Double_t scaleF = 1;
1286
1287     //if(strBgrID.Contains("Out2Jets")){  // scale by ratio 2 jets events / all events
1288     //  scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX()) 
1289     //  / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1290
1291
1292     if(strBgrID.Contains("OutAllJets") ){  // scale by ratio >3 jets events / all events
1293
1294       if(fh1NJets){
1295         scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX()) 
1296           / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1297       }
1298       else{
1299         Printf("%s:%d -- use bgr OutAllJets, but histo fhn1NJets not found",(char*)__FILE__,__LINE__);    
1300       }
1301     }
1302     
1303  
1304     fNJetsBgr->SetAt(nJetsBin*scaleF,i); 
1305
1306     if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1307                         static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1308
1309   }
1310   
1311   // projections: FF 
1312   
1313   for(Int_t i=0; i<fNJetPtSlices; i++){
1314     
1315     Float_t jetPtLoLim = fJetPtSlices->At(i);
1316     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1317
1318     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1319     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1320
1321     if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1322       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
1323       return; 
1324     }
1325     
1326     TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1327     TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1328     TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1329     
1330     // appendix 'unbinned' to avoid histos with same name after rebinning
1331     TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
1332     TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1333     TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1334     
1335     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1336     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1337     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1338
1339     projPt->SetNameTitle(strNameBgrPt,"");
1340     projZ->SetNameTitle(strNameBgrZ,"");
1341     projXi->SetNameTitle(strNameBgrXi,"");
1342     
1343     // raw bgr = corr level 0
1344     fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1345   }
1346
1347   delete fh1FFJetPtBgr;
1348   if(fh1NJets) delete fh1NJets;
1349   delete fh2FFTrackPtBgr;
1350   delete fh2FFZBgr;  
1351   delete fh2FFXiBgr; 
1352 }
1353
1354 //_____________________________________________________________________________________________________________________
1355 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1356
1357   // read raw FF - standard dir/list name
1358      
1359   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
1360   TString strlist = "fracfunc_" + strID;
1361     
1362   ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1363 }
1364
1365 //_______________________________________________________________________________________________________________________________________
1366 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1367 {
1368   // get raw FF from input file, project in jet pt slice
1369   // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1370   // normalization done separately 
1371
1372   TString strBgrID = "BckgEmbed";
1373   TString strID = strBgrID + strFFID;
1374
1375   TFile f(strfile,"READ");
1376
1377   if(!f.IsOpen()){
1378     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1379     return;
1380   }
1381
1382   if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1383
1384   gDirectory->cd(strdir);
1385
1386   TList* list = 0;
1387   
1388   if(!(list = (TList*) gDirectory->Get(strlist))){
1389     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1390     return;
1391   }
1392
1393   TString hnameNJets = "fh1nRecJetsCuts"; 
1394   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); 
1395   TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1396   TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1397   TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1398
1399   TH1F* fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1400   TH1F* fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
1401   TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1402   TH2F* fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);  
1403   TH2F* fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi); 
1404
1405   if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
1406   if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());      return; }
1407   if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1408   if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
1409   if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }
1410
1411   fh1FFJetPtBgr->SetDirectory(0);
1412   fh1NJets->SetDirectory(0);
1413   fh2FFTrackPtBgr->SetDirectory(0);
1414   fh2FFZBgr->SetDirectory(0);  
1415   fh2FFXiBgr->SetDirectory(0); 
1416
1417   f.Close();  
1418
1419   // nJets per bin
1420
1421   for(Int_t i=0; i<fNJetPtSlices; i++){
1422     
1423     Float_t jetPtLoLim = fJetPtSlices->At(i);
1424     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1425     
1426     Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1427     Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1428                                      
1429     Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1430     Double_t scaleF = 1;
1431
1432     fNJetsBgr->SetAt(nJetsBin*scaleF,i); 
1433
1434     if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1435                         static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1436
1437   }
1438   
1439   // projections: FF 
1440   
1441   for(Int_t i=0; i<fNJetPtSlices; i++){
1442     
1443     Float_t jetPtLoLim = fJetPtSlices->At(i);
1444     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1445
1446     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1447     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1448
1449     if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1450       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
1451       return; 
1452     }
1453     
1454     TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1455     TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1456     TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1457     
1458     // appendix 'unbinned' to avoid histos with same name after rebinning
1459     TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
1460     TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1461     TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1462     
1463     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1464     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1465     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1466
1467     projPt->SetNameTitle(strNameBgrPt,"");
1468     projZ->SetNameTitle(strNameBgrZ,"");
1469     projXi->SetNameTitle(strNameBgrXi,"");
1470     
1471     // raw bgr = corr level 0
1472     fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1473   }
1474
1475   delete fh1FFJetPtBgr;
1476   delete fh1NJets;
1477   delete fh2FFTrackPtBgr;
1478   delete fh2FFZBgr;  
1479   delete fh2FFXiBgr; 
1480 }
1481
1482
1483 //__________________________________________________________________________________________________________
1484 void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1485 {
1486   // write histos to file
1487   // skip histos with 0 entries
1488
1489   TString outfileOption = "RECREATE";
1490   if(updateOutfile) outfileOption = "UPDATE";
1491
1492   TFile f(strfile,outfileOption); 
1493
1494   if(!f.IsOpen()){
1495     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1496     return;
1497   }
1498
1499   if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1500
1501   if(strdir && strdir.Length()){
1502     TDirectory* dir = f.mkdir(strdir);
1503     dir->cd(); 
1504   } 
1505
1506   for(Int_t i=0; i<fNJetPtSlices; i++){
1507     
1508     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1509     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries())       fCorrFF[c]->GetZ(i)->Write();
1510     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries())      fCorrFF[c]->GetXi(i)->Write();
1511
1512     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1513     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries())       fCorrBgr[c]->GetZ(i)->Write();
1514     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries())      fCorrBgr[c]->GetXi(i)->Write();
1515
1516     
1517     if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1518     if(fh1FFZBackFolded[i]       && fh1FFZBackFolded[i]->GetEntries())       fh1FFZBackFolded[i]->Write();
1519     if(fh1FFXiBackFolded[i]      && fh1FFXiBackFolded[i]->GetEntries())      fh1FFXiBackFolded[i]->Write();
1520
1521
1522     if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1523     if(fh1FFRatioZFolded[i]       && fh1FFRatioZFolded[i]->GetEntries())       fh1FFRatioZFolded[i]->Write();
1524     if(fh1FFRatioXiFolded[i]      && fh1FFRatioXiFolded[i]->GetEntries())      fh1FFRatioXiFolded[i]->Write();
1525
1526     if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1527     if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries())             fh1FFRatioZBackFolded[i]->Write();
1528     if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries())            fh1FFRatioXiBackFolded[i]->Write();
1529
1530     if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetEntries()) fh1FoldingCorrPt[i]->Write(); // TEST!!!
1531     if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetEntries())     fh1FoldingCorrZ[i]->Write();  // TEST!!!
1532     if(fh1FoldingCorrXi[i] && fh1FoldingCorrXi[i]->GetEntries()) fh1FoldingCorrXi[i]->Write(); // TEST!!!
1533   }
1534   
1535   // inclusive track pt
1536
1537   for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1538   if(fh1SingleTrackPtBackFolded)      fh1SingleTrackPtBackFolded->Write();  
1539   if(fh1RatioSingleTrackPtFolded)     fh1RatioSingleTrackPtFolded->Write();  
1540   if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();  
1541
1542   
1543
1544
1545   f.Close();  
1546 }
1547
1548 //____________________________________________________________________________________________________________________________________
1549 THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1550 {
1551   // copy 1-dimensional histo to THnSparse 
1552   // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1553   // histos with variable bin size are supported
1554
1555   // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1556
1557   THnSparseF* fhnHist;
1558
1559   Int_t nBins       = hist->GetXaxis()->GetNbins();
1560   Int_t nBinsVec[1] = { nBins };
1561
1562   const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1563   
1564   if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1565
1566     fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1567     fhnHist->SetBinEdges(0,binsVec);
1568   }
1569   else{ // uniform bin size
1570     
1571     Double_t xMin  = hist->GetXaxis()->GetXmin();
1572     Double_t xMax  = hist->GetXaxis()->GetXmax();
1573     
1574     Double_t binMinVec[1] = { xMin };
1575     Double_t binMaxVec[1] = { xMax };
1576     
1577     fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1578   }
1579  
1580
1581   for(Int_t bin=1; bin<=nBins; bin++){ 
1582
1583     Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1584
1585     Double_t binCoord[] = {binCenter}; 
1586     fhnHist->Fill(binCoord,1); // initially need to create the bin 
1587
1588     Long64_t binIndex = fhnHist->GetBin(binCoord);
1589
1590     Double_t cont =  hist->GetBinContent(bin);
1591     Double_t err  =  hist->GetBinError(bin);
1592     
1593     if(fillConst){
1594       cont = 1;
1595       err  = 0;
1596     }
1597
1598     fhnHist->SetBinContent(binIndex,cont);
1599     fhnHist->SetBinError(binIndex,err);
1600  }
1601
1602  return fhnHist;
1603 }
1604
1605 //______________________________________________________________________________________________________________________________________________
1606 TH1F* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1607                                                   const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
1608 {
1609   // unfold input histo 
1610
1611   AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior);  // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1612   unfolding.SetMaxNumberOfIterations(nIter);
1613   // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1614   // if(useSmoothing) unfolding.UseSmoothing(); 
1615
1616   if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1617   unfolding.Unfold();
1618
1619   THnSparse* unfolded = unfolding.GetUnfolded();
1620
1621   TString hnameUnf = hnHist->GetName();
1622   hnameUnf.Append("_unf");
1623   unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
1624   
1625   TH1F* h1Unfolded = (TH1F*) unfolded->Projection(0); 
1626   h1Unfolded->SetNameTitle(hnHist->GetName(),hnHist->GetTitle());
1627   
1628   return h1Unfolded;
1629 }
1630
1631 //___________________________________________________________________________________________________________________________
1632 void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1633 {
1634   // loop over jet pt slices and unfold dN/dpt spectra
1635   
1636   TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel(); 
1637   if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1638   
1639   for(Int_t i=0; i<fNJetPtSlices; i++){
1640  
1641     TH1F* hist = 0;
1642     if(type == kFlagPt)      hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1643     else if(type == kFlagZ)  hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i);       // level -2: before unfolding, level -1: unfolded
1644     else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i);      // level -2: before unfolding, level -1: unfolded
1645     else{ 
1646       Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
1647       return;
1648     }
1649     
1650     if(!hist){
1651       Printf("%s%d no histo found ",(char*)__FILE__,__LINE__);
1652       return;
1653     }
1654
1655
1656     THnSparse* hnResponse = 0;
1657     if(type == kFlagPt)      hnResponse = fhnResponsePt[i];
1658     else if(type == kFlagZ)  hnResponse = fhnResponseZ[i];
1659     else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1660
1661     TH1F* hPrior = 0;
1662     if(type == kFlagPt && fh1FFTrackPtPrior[i]  && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1663     else if(type == kFlagZ  && fh1FFZPrior[i]   && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0)       ) hPrior = fh1FFZPrior[i];
1664     else if(type == kFlagXi && fh1FFXiPrior[i]  && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0)      ) hPrior = fh1FFXiPrior[i];
1665
1666
1667     TString histNameTHn;
1668     histNameTHn = hist->GetName();
1669     histNameTHn.ReplaceAll("TH1","THn");
1670     
1671
1672     TString priorNameTHn; 
1673     if(hPrior){
1674       priorNameTHn = hPrior->GetName();
1675       priorNameTHn.ReplaceAll("TH1","THn");
1676     }
1677
1678     TString histNameBackFolded;
1679     histNameBackFolded = hist->GetName();
1680     histNameBackFolded.Append("_backfold");
1681     
1682     TString histNameRatioFolded;
1683     histNameRatioFolded = hist->GetName();
1684     histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1685     
1686     histNameRatioFolded.Append("_unfold");
1687
1688     TString histNameRatioBackFolded;
1689     histNameRatioBackFolded = hist->GetName();
1690     histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1691     histNameRatioBackFolded.Append("_backfold");
1692  
1693     THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1694     THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff 
1695     THnSparse* hnPrior          = 0;
1696     if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1697     
1698     //THnSparse* hnUnfolded 
1699     //  = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);  
1700      
1701     //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
1702     //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1703     
1704     TH1F* hUnfolded 
1705       = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);  
1706      
1707     //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
1708     //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1709
1710     // errors
1711     for(Int_t bin=1; bin<=hist->GetNbinsX(); bin++){
1712       
1713       Double_t contOrg = hist->GetBinContent(bin);
1714       Double_t errOrg  = hist->GetBinError(bin);
1715       
1716       Double_t contUnf = hUnfolded->GetBinContent(bin);
1717       //Double_t errUnf  = hUnfolded->GetBinError(bin);
1718
1719       Double_t relErrOrg = 0;
1720       if(contOrg) relErrOrg = errOrg/contOrg;
1721       
1722       //       std::cout<<" cont original "<<contOrg<<" err original "<<errOrg<<std::endl;
1723       //       std::cout<<" cont unf "<<contUnf<<" errUnf"<<errUnf<<std::endl;
1724       //       std::cout<<" err/cont original "<<relErrOrg<<std::endl;
1725       //       if(contUnf)  std::cout<<" err/conf unf "<<errUnf/contUnf<<std::endl;
1726       
1727       hUnfolded->SetBinError(bin,relErrOrg*contUnf);
1728       
1729       //       if(hUnfolded->GetBinContent(bin)){
1730       //        std::cout<<" modified err unfolded "<<hUnfolded->GetBinError(bin)/hUnfolded->GetBinContent(bin)<<std::endl; 
1731       //       }
1732     }
1733     
1734
1735     if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1736     if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1737     if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1738
1739     // backfolding: apply response matrix to unfolded spectrum
1740     TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); 
1741     hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1742
1743     if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1744     if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
1745     if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
1746     
1747     // ratio unfolded to original histo 
1748     TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1749     hRatioUnfolded->Reset();
1750     hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1751
1752     if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1753     if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
1754     if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;
1755
1756
1757     // ratio backfolded to original histo
1758     TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1759     hRatioBackFolded->Reset();
1760     hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1761
1762     if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1763     if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
1764     if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
1765     
1766     delete hnHist;
1767     delete hnFlatEfficiency;
1768  }
1769 }
1770
1771 //_____________________________________________________________________________________________________
1772 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1773 {
1774   
1775   Int_t type = kFlagPt;
1776   UnfoldHistos(nIter, useCorrelatedErrors, type);
1777 }
1778
1779 //_____________________________________________________________________________________________________
1780 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1781 {
1782   
1783   Int_t type = kFlagZ;
1784   UnfoldHistos(nIter, useCorrelatedErrors, type);
1785 }
1786
1787 //_____________________________________________________________________________________________________
1788 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1789 {
1790   
1791   Int_t type = kFlagXi;
1792   UnfoldHistos(nIter, useCorrelatedErrors, type);
1793 }
1794
1795
1796 //______________________________________________________________________________________________
1797 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1798 {
1799   // apply (multiply) response matrix to hist 
1800
1801   TH1F* hOut = new TH1F(*hist);
1802   hOut->Reset();
1803
1804   const Int_t axisM = 0; 
1805   const Int_t axisT = 1;
1806  
1807   Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1808   Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1809
1810   // response matrix normalization
1811   // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1812   // so for unfolding proper we leave it to CORRFW ...
1813
1814   Double_t normFacResponse[nBinsT+1];
1815
1816   for(Int_t iT=1; iT<=nBinsT; iT++){
1817
1818     Double_t sumResp = 0;
1819     
1820     for(Int_t iM=1; iM<=nBinsM; iM++){
1821       
1822       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1823       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1824       
1825       Double_t binCoord[] = {coordM,coordT};
1826       
1827       Long64_t binIndex = hnResponse->GetBin(binCoord);
1828       
1829       Double_t resp = hnResponse->GetBinContent(binIndex); 
1830       
1831       sumResp += resp;
1832     }
1833     
1834     normFacResponse[iT] = 0;
1835     if(sumResp) normFacResponse[iT] = 1/sumResp;
1836   }
1837   
1838   
1839   
1840   for(Int_t iM=1; iM<=nBinsM; iM++){
1841     
1842     Double_t contM   = 0;
1843
1844     for(Int_t iT=1; iT<=nBinsT; iT++){
1845
1846       Double_t contT = hist->GetBinContent(iT);
1847       
1848       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1849       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1850
1851       Double_t binCoord[] = {coordM,coordT};
1852       
1853       Long64_t binIndex = hnResponse->GetBin(binCoord);
1854       
1855       Double_t resp = hnResponse->GetBinContent(binIndex); 
1856       
1857       contM   += resp*normFacResponse[iT]*contT; 
1858     }
1859
1860     hOut->SetBinContent(iM,contM);
1861   }
1862
1863   return hOut;
1864 }
1865
1866 //_______________________________________________________________________________________________________
1867 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1868 {
1869   // read reconstruction efficiency from file
1870   // argument strlist optional - read from directory strdir if not specified
1871
1872   // temporary histos to hold histos from file
1873   TH1F* hEffPt[fNJetPtSlices]; 
1874   TH1F* hEffZ[fNJetPtSlices];
1875   TH1F* hEffXi[fNJetPtSlices];
1876   
1877   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1878   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
1879   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1880
1881   TFile f(strfile,"READ");
1882
1883   if(!f.IsOpen()){
1884     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1885     return;
1886   }
1887
1888   if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1889  
1890   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1891
1892   TList* list = 0;
1893
1894   if(strlist && strlist.Length()){
1895    
1896     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1897       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1898       return;
1899     }
1900   }  
1901
1902   for(Int_t i=0; i<fNJetPtSlices; i++){
1903     
1904     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1905     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1906     
1907     TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1908     TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1909     TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1910      
1911    
1912     if(list){
1913       hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt); 
1914       hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ); 
1915       hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi); 
1916     }
1917     else{
1918       hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt); 
1919       hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ); 
1920       hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi); 
1921     }
1922     
1923     if(!hEffPt[i]){
1924       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1925     }
1926   
1927     if(!hEffZ[i]){
1928       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1929     }    
1930
1931     if(!hEffXi[i]){
1932       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1933     }
1934
1935
1936     if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1937     if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1938     if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1939
1940     if(hEffPt[i]) hEffPt[i]->SetDirectory(0); 
1941     if(hEffZ[i])  hEffZ[i]->SetDirectory(0); 
1942     if(hEffXi[i]) hEffXi[i]->SetDirectory(0); 
1943
1944   } // jet slices loop
1945
1946   f.Close();
1947
1948   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1949     if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1950     if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
1951     if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1952   }
1953 }
1954
1955 //___________________________________________________________________________________________________________
1956 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1957 {
1958   // read bgr eff from file
1959   // argument strlist optional - read from directory strdir if not specified
1960  
1961   TH1F* hEffPtBgr[fNJetPtSlices]; 
1962   TH1F* hEffZBgr [fNJetPtSlices];
1963   TH1F* hEffXiBgr[fNJetPtSlices];
1964
1965   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1966   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
1967   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1968
1969
1970   TFile f(strfile,"READ");
1971
1972   if(!f.IsOpen()){
1973     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1974     return;
1975   }
1976
1977   if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1978  
1979   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1980
1981   TList* list = 0;
1982
1983   if(strlist && strlist.Length()){
1984        
1985     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1986       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1987       return;
1988     }
1989   }  
1990
1991   for(Int_t i=0; i<fNJetPtSlices; i++){
1992     
1993     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1994     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1995     
1996     TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1997     TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1998     TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1999   
2000    
2001     if(list){
2002       hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr); 
2003       hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr); 
2004       hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr); 
2005     }
2006     else{
2007       hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr); 
2008       hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr); 
2009       hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr); 
2010     }
2011     
2012     if(!hEffPtBgr[i]){
2013       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
2014     }
2015   
2016     if(!hEffZBgr[i]){
2017       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
2018     }    
2019
2020     if(!hEffXiBgr[i]){
2021       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
2022     }
2023
2024
2025     if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
2026     if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
2027     if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
2028
2029     if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0); 
2030     if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0); 
2031     if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0); 
2032
2033   } // jet slices loop
2034
2035   f.Close();
2036
2037   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
2038     if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
2039     if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
2040     if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
2041   }
2042 }
2043
2044 // ________________________________________________
2045 void AliFragmentationFunctionCorrections::EffCorr()
2046 {
2047   // apply efficiency correction
2048
2049   AddCorrectionLevel("eff");
2050
2051   for(Int_t i=0; i<fNJetPtSlices; i++){
2052
2053     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2054     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2055     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2056
2057     TString histNamePt = histPt->GetName();
2058     TString histNameZ  = histZ->GetName();
2059     TString histNameXi = histXi->GetName();
2060
2061     
2062     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
2063     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
2064     
2065     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
2066     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
2067     
2068     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
2069     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
2070     
2071     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
2072   }
2073 }
2074
2075 //___________________________________________________
2076 void AliFragmentationFunctionCorrections::EffCorrBgr()
2077 {
2078   // apply efficiency correction to bgr distributions
2079
2080   AddCorrectionLevelBgr("eff");
2081
2082   Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
2083
2084   for(Int_t i=0; i<fNJetPtSlices; i++){
2085
2086     TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
2087     TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
2088     TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
2089     
2090     TString histNamePt = histPt->GetName();
2091     TString histNameZ  = histZ->GetName();
2092     TString histNameXi = histXi->GetName();
2093
2094     
2095     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
2096     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
2097     
2098     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
2099     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
2100     
2101     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
2102     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
2103     
2104     fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
2105   }
2106 }
2107
2108 //_____________________________________________________
2109 void AliFragmentationFunctionCorrections::SubtractBgr(Double_t sysErr)
2110 {
2111   // subtract bgr distribution from FF
2112   // the current corr level is used for both 
2113   
2114   AddCorrectionLevel("bgrSub");
2115
2116   fstream ascii_out_pt;
2117   fstream ascii_out_z;
2118   fstream ascii_out_xi;
2119
2120   if(sysErr) ascii_out_pt.open("sysErrUE_pt.txt",std::ios::out);
2121   if(sysErr) ascii_out_z.open("sysErrUE_z.txt",std::ios::out);
2122   if(sysErr) ascii_out_xi.open("sysErrUE_xi.txt",std::ios::out);
2123
2124   for(Int_t i=0; i<fNJetPtSlices; i++){ // jet slices
2125
2126     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2127     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2128     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2129     
2130     TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2131     TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2132     TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2133
2134     TString histNamePt = histPt->GetName();
2135     TString histNameZ  = histZ->GetName();
2136     TString histNameXi = histXi->GetName();
2137     
2138
2139     TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2140     hFFTrackPtBgrSub->Add(histPtBgr,-1);
2141     
2142     TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
2143     hFFZBgrSub->Add(histZBgr,-1);
2144     
2145     TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2146     hFFXiBgrSub->Add(histXiBgr,-1);
2147
2148     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2149
2150     if(sysErr){
2151
2152       for(Int_t bin=1; bin<=histPt->GetNbinsX(); bin++){
2153
2154         //Double_t sigPlBgr = histPt->GetBinContent(bin);
2155         Double_t sig      = hFFTrackPtBgrSub->GetBinContent(bin);
2156         Double_t bgr      = histPtBgr->GetBinContent(bin);
2157         
2158         Double_t relErr = 0; 
2159         if(sig>0)  relErr = sysErr*bgr/sig;
2160
2161 //      std::cout<<" pt bin "<<bin<<" mean "<<histPt->GetBinCenter(bin)
2162 //               <<" sigPlBgr "<<sigPlBgr<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2163         
2164         ascii_out_pt<<i<<" "<<bin<<" "<<relErr<<std::endl;
2165       }
2166
2167
2168       for(Int_t bin=1; bin<=histZ->GetNbinsX(); bin++){
2169         
2170         //Double_t sigPlBgr = histZ->GetBinContent(bin);
2171         Double_t sig      = hFFZBgrSub->GetBinContent(bin);
2172         Double_t bgr      = histZBgr->GetBinContent(bin);
2173         
2174         Double_t relErr = 0; 
2175         if(sig>0)  relErr = sysErr*bgr/sig;
2176         
2177         std::cout<<" z bin "<<bin<<" mean "<<histZ->GetBinCenter(bin)
2178                  <<" sigPlBgr "<<histZ->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2179         
2180         ascii_out_z<<i<<" "<<bin<<" "<<relErr<<std::endl;
2181       }
2182
2183
2184       for(Int_t bin=1; bin<=histXi->GetNbinsX(); bin++){
2185         
2186         //Double_t sigPlBgr = histXi->GetBinContent(bin);
2187         Double_t sig      = hFFXiBgrSub->GetBinContent(bin);
2188         Double_t bgr      = histXiBgr->GetBinContent(bin);
2189         
2190         Double_t relErr = 0; 
2191         if(sig>0)  relErr = sysErr*bgr/sig;
2192         
2193         std::cout<<" xi bin "<<bin<<" mean "<<histXi->GetBinCenter(bin)
2194                  <<" sigPlBgr "<<histXi->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
2195         
2196         ascii_out_xi<<i<<" "<<bin<<" "<<relErr<<std::endl;
2197       }
2198     }
2199   }
2200
2201   if(sysErr) ascii_out_pt.close(); 
2202   if(sysErr) ascii_out_xi.close(); 
2203
2204
2205 }
2206
2207 //________________________________________________________________________________________________________________
2208 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile, 
2209                                                               Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2210
2211   // read task ouput from MC and write single track eff - standard dir/list 
2212      
2213   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2214   TString strlist = "fracfunc_" + strID;
2215     
2216   WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2217 }
2218
2219 //___________________________________________________________________________________________________________________________________
2220 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, 
2221                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2222 {
2223   // read task output from MC and write single track eff as function of pt, eta, phi
2224   // argument strlist optional - read from directory strdir if not specified
2225   // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2226
2227
2228   TH1D* hdNdptTracksMCPrimGen;
2229   TH2D* hdNdetadphiTracksMCPrimGen;
2230   
2231   TH1D* hdNdptTracksMCPrimRec;
2232   TH2D* hdNdetadphiTracksMCPrimRec;
2233     
2234
2235   TFile f(strInfile,"READ");
2236
2237   if(!f.IsOpen()){
2238     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2239     return;
2240   }
2241   
2242   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2243  
2244   if(strdir && strdir.Length()){
2245     if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2246     gDirectory->cd(strdir);
2247   }
2248
2249   TList* list = 0;
2250
2251   if(strlist && strlist.Length()){
2252
2253     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2254       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2255       return;
2256     }
2257   }
2258
2259
2260   TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2261   if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data()); 
2262
2263   TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2264   if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2265
2266   TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2267   if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2268   
2269   TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2270   if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2271   
2272
2273   if(list){
2274     hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
2275     hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2276     
2277     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
2278     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2279   }
2280   else{
2281     hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2282     hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2283
2284     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2285     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2286   }
2287
2288   hdNdptTracksMCPrimGen->SetDirectory(0);
2289   hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2290   hdNdptTracksMCPrimRec->SetDirectory(0);
2291   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2292   
2293   f.Close();
2294
2295   // projections: dN/deta, dN/dphi 
2296
2297   TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2298   TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2299  
2300   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2301   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2302
2303   // rebin
2304
2305   TString strNamePtGen = "hTrackPtGenPrim";
2306   TString strNamePtRec = "hTrackPtRecPrim";
2307
2308   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2309   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2310  
2311     // efficiency: divide 
2312
2313   TString hNameTrackEffPt = "hSingleTrackEffPt";
2314   if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2315                                                
2316   TString hNameTrackEffEta = "hSingleTrackEffEta";
2317   if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2318
2319   TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2320   if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2321
2322
2323   TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2324   hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2325
2326   TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2327   hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2328   
2329   TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2330   hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2331   
2332   
2333   TString outfileOption = "RECREATE";
2334   if(updateOutfile)  outfileOption = "UPDATE";
2335
2336   TFile out(strOutfile,outfileOption);
2337
2338   if(!out.IsOpen()){
2339     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2340     return;
2341   }
2342
2343   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2344
2345   if(strOutDir && strOutDir.Length()){
2346     
2347     TDirectory* dir;
2348     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2349     else{
2350       dir = out.mkdir(strOutDir);
2351       dir->cd(); 
2352     } 
2353   }
2354
2355   hSingleTrackEffPt->Write();
2356   hSingleTrackEffEta->Write();
2357   hSingleTrackEffPhi->Write();
2358   
2359   out.Close();
2360
2361   delete hdNdptTracksMCPrimGen;
2362   delete hdNdetadphiTracksMCPrimGen;
2363   delete hdNdptTracksMCPrimRec;
2364   delete hdNdetadphiTracksMCPrimRec;
2365 }
2366
2367 //________________________________________________________________________________________________________________
2368 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile, 
2369                                                                   Bool_t updateOutfile, TString strOutDir)
2370
2371   // read task ouput from MC and write single track eff - standard dir/list 
2372      
2373   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2374   TString strlist = "fracfunc_" + strID;
2375     
2376   WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2377 }
2378
2379 //___________________________________________________________________________________________________________________________________
2380 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, 
2381                                                                   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2382 {
2383   // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2384   // argument strlist optional - read from directory strdir if not specified
2385   // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2386
2387   TH1D* hdNdptTracksMCPrimRec;
2388   TH2D* hdNdetadphiTracksMCPrimRec;
2389   
2390   TH1D* hdNdptTracksMCSecRecNS;
2391   TH2D* hdNdetadphiTracksMCSecRecNS;
2392
2393   TH1D* hdNdptTracksMCSecRecSsc;
2394   TH2D* hdNdetadphiTracksMCSecRecSsc;
2395
2396   TFile f(strInfile,"READ");
2397
2398   if(!f.IsOpen()){
2399     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2400     return;
2401   }
2402   
2403   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2404  
2405   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2406
2407   TList* list = 0;
2408
2409   if(strlist && strlist.Length()){
2410
2411     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2412       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2413       return;
2414     }
2415   }
2416
2417
2418   if(list){
2419     hdNdptTracksMCPrimRec        = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2420     hdNdetadphiTracksMCPrimRec   = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2421     
2422     hdNdptTracksMCSecRecNS       = (TH1D*) list->FindObject("fh1TrackQAPtSecRecNS");
2423     hdNdetadphiTracksMCSecRecNS  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecNS");
2424
2425     hdNdptTracksMCSecRecSsc      = (TH1D*) list->FindObject("fh1TrackQAPtSecRecSsc");
2426     hdNdetadphiTracksMCSecRecSsc = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecSsc");
2427
2428   }
2429   else{
2430     hdNdptTracksMCPrimRec        = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2431     hdNdetadphiTracksMCPrimRec   = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2432
2433     hdNdptTracksMCSecRecNS       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecNS");
2434     hdNdetadphiTracksMCSecRecNS  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecNS");
2435
2436     hdNdptTracksMCSecRecSsc      = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecSsc");
2437     hdNdetadphiTracksMCSecRecSsc = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecSsc");
2438   }
2439   
2440   hdNdptTracksMCPrimRec->SetDirectory(0);
2441   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2442
2443   hdNdptTracksMCSecRecNS->SetDirectory(0);
2444   hdNdetadphiTracksMCSecRecNS->SetDirectory(0);
2445
2446   hdNdptTracksMCSecRecSsc->SetDirectory(0);
2447   hdNdetadphiTracksMCSecRecSsc->SetDirectory(0);
2448   
2449   f.Close();
2450
2451   // projections: dN/deta, dN/dphi 
2452
2453   TH1D* hdNdetaTracksMCPrimRec   = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2454   TH1D* hdNdphiTracksMCPrimRec   = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2455  
2456   TH1D* hdNdetaTracksMCSecRecNS  = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionX("hdNdetaTracksMcSecRecNS");
2457   TH1D* hdNdphiTracksMCSecRecNS  = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionY("hdNdphiTracksMcSecRecNS");
2458
2459   TH1D* hdNdetaTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionX("hdNdetaTracksMcSecRecSsc");
2460   TH1D* hdNdphiTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionY("hdNdphiTracksMcSecRecSsc");
2461  
2462   // rebin
2463
2464   TString strNamePtPrim = "hTrackPtPrim";
2465   TString strNamePtSec  = "hTrackPtSec";
2466
2467   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec   = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2468   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecNS  = (TH1D*) hdNdptTracksMCSecRecNS->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2469   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecSsc = (TH1D*) hdNdptTracksMCSecRecSsc->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2470
2471
2472   // secondary correction factor: divide prim/(prim+sec)
2473
2474   TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSingleTrackSecCorrPt");
2475   TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSumPrimSecPt");
2476   hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS);
2477   hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc);
2478   hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2479
2480
2481   TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSingleTrackSecCorrEta");
2482   TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSumPrimSecEta");
2483   hSumPrimSecEta->Add(hdNdetaTracksMCSecRecNS);
2484   hSumPrimSecEta->Add(hdNdetaTracksMCSecRecSsc);
2485   hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2486
2487   TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSingleTrackSecCorrPhi");
2488   TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSumPrimSecPhi");
2489   hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecNS);
2490   hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecSsc);
2491   hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2492
2493   // 
2494   
2495   TString outfileOption = "RECREATE";
2496   if(updateOutfile)  outfileOption = "UPDATE";
2497
2498   TFile out(strOutfile,outfileOption);
2499
2500   if(!out.IsOpen()){
2501     Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2502     return;
2503   }
2504
2505   if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2506
2507   if(strOutDir && strOutDir.Length()){  
2508
2509     TDirectory* dir;
2510     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2511     else{
2512       dir = out.mkdir(strOutDir);
2513       dir->cd(); 
2514     } 
2515   } 
2516
2517   hdNdptTracksMCSecRecNS->Write();
2518   hdNdetaTracksMCSecRecNS->Write();
2519   hdNdphiTracksMCSecRecNS->Write();
2520
2521   hdNdptTracksMCSecRecSsc->Write();
2522   hdNdetaTracksMCSecRecSsc->Write();
2523   hdNdphiTracksMCSecRecSsc->Write();
2524
2525   hSingleTrackSecCorrPt->Write();
2526   hSingleTrackSecCorrEta->Write();
2527   hSingleTrackSecCorrPhi->Write();
2528   
2529   out.Close();
2530
2531   delete hdNdptTracksMCPrimRec;       
2532   delete hdNdetadphiTracksMCPrimRec;  
2533   delete hdNdptTracksMCSecRecNS;
2534   delete hdNdetadphiTracksMCSecRecNS;
2535   delete hdNdptTracksMCSecRecSsc;
2536   delete hdNdetadphiTracksMCSecRecSsc;
2537 }
2538
2539 //________________________________________________________________________________________________________________
2540 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile, 
2541                                                               Bool_t updateOutfile, TString strOutDir)
2542
2543   // read task ouput from MC and write single track eff - standard dir/list 
2544      
2545   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2546   TString strlist = "fracfunc_" + strID;
2547     
2548   WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2549 }
2550
2551
2552 //_____________________________________________________________________________________________________________________________________
2553 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2554                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2555 {
2556   // read 2d THnSparse response matrices in pt from file
2557   // project TH2 
2558   // write to strOutfile 
2559
2560   THnSparse* hnResponseSinglePt; 
2561   TH2F*      h2ResponseSinglePt; 
2562  
2563   TFile f(strInfile,"READ");
2564
2565   if(!f.IsOpen()){
2566     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2567     return;
2568   }
2569
2570   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2571  
2572   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2573
2574   TList* list = 0;
2575
2576   if(strlist && strlist.Length()){
2577     
2578     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2579       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2580       return;
2581     }
2582   }  
2583   
2584   if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2585   else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2586   
2587
2588   if(!hnResponseSinglePt){
2589     Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2590     return;
2591   }
2592
2593   f.Close();
2594
2595
2596   // project 2d histo 
2597   h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2598   h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt",""); 
2599     
2600   
2601   // write 
2602
2603   TString outfileOption = "RECREATE";
2604   if(updateOutfile)  outfileOption = "UPDATE";
2605   
2606   TFile out(strOutfile,outfileOption);
2607   
2608   if(!out.IsOpen()){
2609     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2610     return;
2611   }
2612
2613   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2614   
2615   if(strOutDir && strOutDir.Length()){  
2616
2617     TDirectory* dir;
2618     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
2619     else{
2620       dir = out.mkdir(strOutDir);
2621       dir->cd(); 
2622     } 
2623   }
2624   
2625   hnResponseSinglePt->Write();
2626   h2ResponseSinglePt->Write();
2627   
2628   out.Close();  
2629 }
2630
2631 //________________________________________________________________________________________________________________
2632 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile, 
2633                                                            Bool_t updateOutfile)
2634
2635   // read task ouput from MC and write single track eff - standard dir/list 
2636      
2637   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2638   TString strlist = "fracfunc_" + strID;
2639     
2640   WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2641 }
2642
2643 //___________________________________________________________________________________________________________________________________
2644 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, 
2645                                                            TString strOutfile, Bool_t updateOutfile)
2646 {
2647   // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2648   // argument strlist optional - read from directory strdir if not specified
2649   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2650
2651   TH1F* hEffPt[fNJetPtSlices];
2652   TH1F* hEffXi[fNJetPtSlices];
2653   TH1F* hEffZ[fNJetPtSlices];
2654
2655   TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2656   TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2657   TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2658     
2659   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2660   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2661   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2662
2663
2664   TH1F* fh1FFJetPtRecEffGen;
2665
2666   TH2F* fh2FFTrackPtRecEffGen;
2667   TH2F* fh2FFZRecEffGen;
2668   TH2F* fh2FFXiRecEffGen;
2669   
2670   TH2F* fh2FFTrackPtRecEffRec;
2671   TH2F* fh2FFZRecEffRec;
2672   TH2F* fh2FFXiRecEffRec;
2673  
2674
2675   TFile f(strInfile,"READ");
2676
2677   if(!f.IsOpen()){
2678     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2679     return;
2680   }
2681   
2682   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2683  
2684   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2685
2686   TList* list = 0;
2687
2688   if(strlist && strlist.Length()){
2689
2690     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2691       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2692       return;
2693     }
2694   }
2695
2696   if(list){
2697     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2698
2699     fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2700     fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2701     fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2702     
2703     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2704     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2705     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2706   }
2707   else{
2708     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffGen");
2709
2710     fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffGen");
2711     fh2FFZRecEffGen        = (TH2F*) gDirectory->Get("fh2FFZRecEffGen");
2712     fh2FFXiRecEffGen       = (TH2F*) gDirectory->Get("fh2FFXiRecEffGen");
2713     
2714     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
2715     fh2FFZRecEffRec        = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
2716     fh2FFXiRecEffRec       = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
2717   }
2718   
2719   fh1FFJetPtRecEffGen->SetDirectory(0); 
2720
2721   fh2FFTrackPtRecEffGen->SetDirectory(0);
2722   fh2FFZRecEffGen->SetDirectory(0);
2723   fh2FFXiRecEffGen->SetDirectory(0);
2724   
2725   fh2FFTrackPtRecEffRec->SetDirectory(0);
2726   fh2FFZRecEffRec->SetDirectory(0);
2727   fh2FFXiRecEffRec->SetDirectory(0);
2728   
2729   f.Close();
2730
2731
2732   // projections: FF for generated and reconstructed primaries 
2733   
2734   for(Int_t i=0; i<fNJetPtSlices; i++){
2735     
2736     Float_t jetPtLoLim = fJetPtSlices->At(i);
2737     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2738
2739     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2740     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2741
2742     if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2743       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2744       return; 
2745     }
2746     
2747     TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2748     TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2749     TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2750     
2751     TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2752     TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2753     TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2754     
2755     // project 
2756     // appendix 'unbinned' to avoid histos with same name after rebinning
2757
2758     hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2759     hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2760     hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2761     
2762     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2763     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2764     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2765     
2766     // rebin
2767
2768     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2769     if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2770     if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2771
2772     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2773     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2774     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2775
2776     hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2777     hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2778     hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2779     
2780     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2781     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2782     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2783  
2784     // normalize
2785     
2786     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2787
2788     NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin); 
2789     NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin); 
2790     NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin); 
2791
2792     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2793     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2794     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2795     
2796     // divide rec/gen : efficiency
2797
2798     TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2799     TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2800     TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2801  
2802     hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2803     hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2804     
2805     hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2806     hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2807     
2808     hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2809     hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2810   } 
2811   
2812   // write 
2813
2814   TString outfileOption = "RECREATE";
2815   if(updateOutfile)  outfileOption = "UPDATE";
2816
2817   TFile out(strOutfile,outfileOption);
2818   
2819   if(!out.IsOpen()){
2820     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2821     return;
2822   }
2823
2824   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2825
2826 //   if(strdir && strdir.Length()){ 
2827 //     TDirectory* dir = out.mkdir(strdir);
2828 //     dir->cd(); 
2829 //   }
2830
2831   for(Int_t i=0; i<fNJetPtSlices; i++){
2832
2833     hdNdptTracksMCPrimGen[i]->Write();
2834     hdNdxiMCPrimGen[i]->Write();
2835     hdNdzMCPrimGen[i]->Write();
2836     
2837     hdNdptTracksMCPrimRec[i]->Write();
2838     hdNdxiMCPrimRec[i]->Write();
2839     hdNdzMCPrimRec[i]->Write();
2840   
2841     hEffPt[i]->Write();
2842     hEffXi[i]->Write();
2843     hEffZ[i]->Write();
2844   }
2845
2846   out.Close();
2847
2848  
2849   delete fh1FFJetPtRecEffGen; 
2850
2851   delete fh2FFTrackPtRecEffGen;
2852   delete fh2FFZRecEffGen;
2853   delete fh2FFXiRecEffGen;
2854   
2855   delete fh2FFTrackPtRecEffRec;
2856   delete fh2FFZRecEffRec;
2857   delete fh2FFXiRecEffRec;
2858 }
2859
2860 //________________________________________________________________________________________________________________
2861 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile, 
2862                                                           Bool_t updateOutfile, TString strOutDir)
2863
2864   // read task ouput from MC and write secondary correction - standard dir/list 
2865      
2866   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2867   TString strlist = "fracfunc_" + strID;
2868     
2869   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kFALSE,"");
2870 }
2871
2872 //________________________________________________________________________________________________________________
2873 void AliFragmentationFunctionCorrections::WriteBgrJetSecCorr(TString strInfile, TString strBgrID, TString strID, TString strOutfile, 
2874                                                              Bool_t updateOutfile, TString strOutDir,Double_t scaleFacBgrRec)
2875
2876   // read task ouput from MC and write secondary correction - standard dir/list 
2877      
2878   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2879   TString strlist = "fracfunc_" + strID;
2880     
2881   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kTRUE,strBgrID,scaleFacBgrRec);
2882 }
2883
2884
2885 //___________________________________________________________________________________________________________________________________
2886 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, 
2887                                                           TString strOutfile, Bool_t updateOutfile, TString strOutDir,
2888                                                           Bool_t writeBgr, TString strBgrID,Double_t scaleFacBgrRec)
2889 {
2890   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2891   // argument strlist optional - read from directory strdir if not specified
2892   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2893
2894   TH1F* hSecCorrPt[fNJetPtSlices];
2895   TH1F* hSecCorrXi[fNJetPtSlices];
2896   TH1F* hSecCorrZ[fNJetPtSlices];
2897
2898   // corr factors using naive Pythia strangeness - for sys uncertainties
2899   TH1F* hSecCorrPt_nonSc[fNJetPtSlices];
2900   TH1F* hSecCorrXi_nonSc[fNJetPtSlices];
2901   TH1F* hSecCorrZ_nonSc[fNJetPtSlices];
2902
2903   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2904   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2905   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2906     
2907   TH1F* hdNdptTracksMCSecRecNS[fNJetPtSlices];
2908   TH1F* hdNdxiMCSecRecNS[fNJetPtSlices];
2909   TH1F* hdNdzMCSecRecNS[fNJetPtSlices];
2910
2911   TH1F* hdNdptTracksMCSecRecS[fNJetPtSlices];
2912   TH1F* hdNdxiMCSecRecS[fNJetPtSlices];
2913   TH1F* hdNdzMCSecRecS[fNJetPtSlices];
2914
2915   TH1F* hdNdptTracksMCSecRecSsc[fNJetPtSlices];
2916   TH1F* hdNdxiMCSecRecSsc[fNJetPtSlices];
2917   TH1F* hdNdzMCSecRecSsc[fNJetPtSlices];
2918   
2919   // ---
2920
2921   TH1F* fh1FFJetPtRecEffRec;
2922
2923   TH2F* fh2FFTrackPtRecEffRec;
2924   TH2F* fh2FFZRecEffRec;
2925   TH2F* fh2FFXiRecEffRec;
2926   
2927   TH2F* fh2FFTrackPtSecRecNS;
2928   TH2F* fh2FFZSecRecNS;
2929   TH2F* fh2FFXiSecRecNS;
2930
2931   TH2F* fh2FFTrackPtSecRecS;
2932   TH2F* fh2FFZSecRecS;
2933   TH2F* fh2FFXiSecRecS;
2934
2935   TH2F* fh2FFTrackPtSecRecSsc;
2936   TH2F* fh2FFZSecRecSsc;
2937   TH2F* fh2FFXiSecRecSsc;
2938   
2939   // ---
2940
2941   TFile f(strInfile,"READ");
2942
2943   if(!f.IsOpen()){
2944     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2945     return;
2946   }
2947   
2948   if(fDebug>0) Printf("%s:%d -- writeJetTrackSecCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2949  
2950   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2951
2952   TList* list = 0;
2953
2954   if(strlist && strlist.Length()){
2955
2956     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2957       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2958       return;
2959     }
2960   }
2961
2962   if(list){
2963     fh1FFJetPtRecEffRec    = (TH1F*) list->FindObject("fh1FFJetPtRecEffRec");
2964
2965     if(writeBgr){ 
2966       fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
2967       fh2FFZRecEffRec        = (TH2F*) list->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
2968       fh2FFXiRecEffRec       = (TH2F*) list->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data())); 
2969
2970       fh2FFTrackPtSecRecNS   = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
2971       fh2FFZSecRecNS         = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
2972       fh2FFXiSecRecNS        = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));
2973
2974       fh2FFTrackPtSecRecS    = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
2975       fh2FFZSecRecS          = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
2976       fh2FFXiSecRecS         = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
2977       
2978       fh2FFTrackPtSecRecSsc  = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
2979       fh2FFZSecRecSsc        = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
2980       fh2FFXiSecRecSsc       = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
2981     }
2982     else{
2983       fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2984       fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2985       fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2986       
2987       fh2FFTrackPtSecRecNS   = (TH2F*) list->FindObject("fh2FFTrackPtSecRecNS");
2988       fh2FFZSecRecNS         = (TH2F*) list->FindObject("fh2FFZSecRecNS");
2989       fh2FFXiSecRecNS        = (TH2F*) list->FindObject("fh2FFXiSecRecNS");
2990       
2991       fh2FFTrackPtSecRecS    = (TH2F*) list->FindObject("fh2FFTrackPtSecRecS");
2992       fh2FFZSecRecS          = (TH2F*) list->FindObject("fh2FFZSecRecS");
2993       fh2FFXiSecRecS         = (TH2F*) list->FindObject("fh2FFXiSecRecS");
2994       
2995       fh2FFTrackPtSecRecSsc  = (TH2F*) list->FindObject("fh2FFTrackPtSecRecSsc");
2996       fh2FFZSecRecSsc        = (TH2F*) list->FindObject("fh2FFZSecRecSsc");
2997       fh2FFXiSecRecSsc       = (TH2F*) list->FindObject("fh2FFXiSecRecSsc");
2998     }
2999   }
3000   else{
3001     fh1FFJetPtRecEffRec    = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");
3002
3003     if(writeBgr){
3004       fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
3005       fh2FFZRecEffRec        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
3006       fh2FFXiRecEffRec       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data())); 
3007
3008       fh2FFTrackPtSecRecNS   = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
3009       fh2FFZSecRecNS         = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
3010       fh2FFXiSecRecNS        = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));
3011
3012       fh2FFTrackPtSecRecS    = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
3013       fh2FFZSecRecS          = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
3014       fh2FFXiSecRecS         = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
3015       
3016       fh2FFTrackPtSecRecSsc  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
3017       fh2FFZSecRecSsc        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
3018       fh2FFXiSecRecSsc       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
3019     }
3020     else{
3021       fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
3022       fh2FFZRecEffRec        = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
3023       fh2FFXiRecEffRec       = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
3024     
3025       fh2FFTrackPtSecRecNS   = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecNS");
3026       fh2FFZSecRecNS         = (TH2F*) gDirectory->Get("fh2FFZSecRecNS");
3027       fh2FFXiSecRecNS        = (TH2F*) gDirectory->Get("fh2FFXiSecRecNS");
3028       
3029       fh2FFTrackPtSecRecS    = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecS");
3030       fh2FFZSecRecS          = (TH2F*) gDirectory->Get("fh2FFZSecRecS");
3031       fh2FFXiSecRecS         = (TH2F*) gDirectory->Get("fh2FFXiSecRecS");
3032       
3033       fh2FFTrackPtSecRecSsc  = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecSsc");
3034       fh2FFZSecRecSsc        = (TH2F*) gDirectory->Get("fh2FFZSecRecSsc");
3035       fh2FFXiSecRecSsc       = (TH2F*) gDirectory->Get("fh2FFXiSecRecSsc");
3036     }
3037   }
3038   
3039   fh1FFJetPtRecEffRec->SetDirectory(0); 
3040    
3041   fh2FFTrackPtRecEffRec->SetDirectory(0);
3042   fh2FFZRecEffRec->SetDirectory(0);
3043   fh2FFXiRecEffRec->SetDirectory(0);
3044   
3045   fh2FFTrackPtSecRecNS->SetDirectory(0);
3046   fh2FFZSecRecNS->SetDirectory(0);
3047   fh2FFXiSecRecNS->SetDirectory(0);
3048   
3049   fh2FFTrackPtSecRecS->SetDirectory(0);
3050   fh2FFZSecRecS->SetDirectory(0);
3051   fh2FFXiSecRecS->SetDirectory(0);
3052   
3053   fh2FFTrackPtSecRecSsc->SetDirectory(0);
3054   fh2FFZSecRecSsc->SetDirectory(0);
3055   fh2FFXiSecRecSsc->SetDirectory(0);
3056
3057   f.Close();
3058
3059
3060   // projections: FF for generated and reconstructed primaries 
3061   
3062   for(Int_t i=0; i<fNJetPtSlices; i++){
3063     
3064     Float_t jetPtLoLim = fJetPtSlices->At(i);
3065     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
3066
3067     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
3068     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
3069
3070     if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
3071       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
3072       return; 
3073     }
3074     
3075     TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3076     TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3077     TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3078     
3079     TString strNameFFPtSecRecNS(Form("fh1FFTrackPtRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3080     TString strNameFFZSecRecNS(Form("fh1FFZRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3081     TString strNameFFXiSecRecNS(Form("fh1FFXiRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3082      
3083     TString strNameFFPtSecRecS(Form("fh1FFTrackPtRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3084     TString strNameFFZSecRecS(Form("fh1FFZRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3085     TString strNameFFXiSecRecS(Form("fh1FFXiRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3086      
3087     TString strNameFFPtSecRecSsc(Form("fh1FFTrackPtRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3088     TString strNameFFZSecRecSsc(Form("fh1FFZRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3089     TString strNameFFXiSecRecSsc(Form("fh1FFXiRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3090
3091     if(writeBgr){
3092       strNameFFPtPrimRec.Form("fh1BgrTrackPt%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3093       strNameFFZPrimRec.Form("fh1BgrZ%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3094       strNameFFXiPrimRec.Form("fh1BgrXi%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3095       
3096       strNameFFPtSecRecNS.Form("fh1BgrTrackPt%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3097       strNameFFZSecRecNS.Form("fh1BgrZ%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3098       strNameFFXiSecRecNS.Form("fh1BgrXi%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3099       
3100       strNameFFPtSecRecS.Form("fh1BgrTrackPt%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3101       strNameFFZSecRecS.Form("fh1BgrZ%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3102       strNameFFXiSecRecS.Form("fh1BgrXi%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3103       
3104       strNameFFPtSecRecSsc.Form("fh1BgrTrackPt%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3105       strNameFFZSecRecSsc.Form("fh1BgrZ%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3106       strNameFFXiSecRecSsc.Form("fh1BgrXi%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
3107     }
3108     
3109     // project 
3110     // appendix 'unbinned' to avoid histos with same name after rebinning
3111
3112     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
3113     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
3114     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
3115     
3116     hdNdptTracksMCSecRecNS[i]  = (TH1F*) fh2FFTrackPtSecRecNS->ProjectionY(strNameFFPtSecRecNS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
3117     hdNdzMCSecRecNS[i]         = (TH1F*) fh2FFZSecRecNS->ProjectionY(strNameFFZSecRecNS+"_unBinned",binLo,binUp,"o");
3118     hdNdxiMCSecRecNS[i]        = (TH1F*) fh2FFXiSecRecNS->ProjectionY(strNameFFXiSecRecNS+"_unBinned",binLo,binUp,"o");
3119     
3120     hdNdptTracksMCSecRecS[i]  = (TH1F*) fh2FFTrackPtSecRecS->ProjectionY(strNameFFPtSecRecS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
3121     hdNdzMCSecRecS[i]         = (TH1F*) fh2FFZSecRecS->ProjectionY(strNameFFZSecRecS+"_unBinned",binLo,binUp,"o");
3122     hdNdxiMCSecRecS[i]        = (TH1F*) fh2FFXiSecRecS->ProjectionY(strNameFFXiSecRecS+"_unBinned",binLo,binUp,"o");
3123
3124     hdNdptTracksMCSecRecSsc[i]  = (TH1F*) fh2FFTrackPtSecRecSsc->ProjectionY(strNameFFPtSecRecSsc+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
3125     hdNdzMCSecRecSsc[i]         = (TH1F*) fh2FFZSecRecSsc->ProjectionY(strNameFFZSecRecSsc+"_unBinned",binLo,binUp,"o");
3126     hdNdxiMCSecRecSsc[i]        = (TH1F*) fh2FFXiSecRecSsc->ProjectionY(strNameFFXiSecRecSsc+"_unBinned",binLo,binUp,"o");
3127
3128
3129     // rebin
3130     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
3131     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
3132     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
3133
3134     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecNS[i] = (TH1F*) hdNdptTracksMCSecRecNS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecNS,fHistoBinsPt[i]->GetArray());
3135     if(fNHistoBinsZ[i])  hdNdzMCSecRecNS[i]        = (TH1F*) hdNdzMCSecRecNS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecNS,fHistoBinsZ[i]->GetArray());
3136     if(fNHistoBinsXi[i]) hdNdxiMCSecRecNS[i]       = (TH1F*) hdNdxiMCSecRecNS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecNS,fHistoBinsXi[i]->GetArray());
3137
3138     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecS[i] = (TH1F*) hdNdptTracksMCSecRecS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecS,fHistoBinsPt[i]->GetArray());
3139     if(fNHistoBinsZ[i])  hdNdzMCSecRecS[i]        = (TH1F*) hdNdzMCSecRecS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecS,fHistoBinsZ[i]->GetArray());
3140     if(fNHistoBinsXi[i]) hdNdxiMCSecRecS[i]       = (TH1F*) hdNdxiMCSecRecS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecS,fHistoBinsXi[i]->GetArray());
3141
3142     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecSsc[i] = (TH1F*) hdNdptTracksMCSecRecSsc[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecSsc,fHistoBinsPt[i]->GetArray());
3143     if(fNHistoBinsZ[i])  hdNdzMCSecRecSsc[i]        = (TH1F*) hdNdzMCSecRecSsc[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecSsc,fHistoBinsZ[i]->GetArray());
3144     if(fNHistoBinsXi[i]) hdNdxiMCSecRecSsc[i]       = (TH1F*) hdNdxiMCSecRecSsc[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecSsc,fHistoBinsXi[i]->GetArray());
3145
3146
3147     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
3148     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
3149     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
3150     
3151     hdNdptTracksMCSecRecNS[i]->SetNameTitle(strNameFFPtSecRecNS,"");
3152     hdNdzMCSecRecNS[i]->SetNameTitle(strNameFFZSecRecNS,"");
3153     hdNdxiMCSecRecNS[i]->SetNameTitle(strNameFFXiSecRecNS,"");
3154     
3155     hdNdptTracksMCSecRecS[i]->SetNameTitle(strNameFFPtSecRecS,"");
3156     hdNdzMCSecRecS[i]->SetNameTitle(strNameFFZSecRecS,"");
3157     hdNdxiMCSecRecS[i]->SetNameTitle(strNameFFXiSecRecS,"");
3158
3159     hdNdptTracksMCSecRecSsc[i]->SetNameTitle(strNameFFPtSecRecSsc,"");
3160     hdNdzMCSecRecSsc[i]->SetNameTitle(strNameFFZSecRecSsc,"");
3161     hdNdxiMCSecRecSsc[i]->SetNameTitle(strNameFFXiSecRecSsc,"");
3162     
3163     // normalize 
3164     Double_t nJetsBin = fh1FFJetPtRecEffRec->Integral(binLo,binUp);
3165
3166     // scale fac for perp2 bgr
3167     if(writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBin /= scaleFacBgrRec;
3168
3169     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
3170     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
3171     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
3172
3173     NormalizeTH1(hdNdptTracksMCSecRecNS[i],nJetsBin); 
3174     NormalizeTH1(hdNdzMCSecRecNS[i],nJetsBin); 
3175     NormalizeTH1(hdNdxiMCSecRecNS[i],nJetsBin); 
3176
3177     NormalizeTH1(hdNdptTracksMCSecRecS[i],nJetsBin); 
3178     NormalizeTH1(hdNdzMCSecRecS[i],nJetsBin); 
3179     NormalizeTH1(hdNdxiMCSecRecS[i],nJetsBin); 
3180
3181     NormalizeTH1(hdNdptTracksMCSecRecSsc[i],nJetsBin); 
3182     NormalizeTH1(hdNdzMCSecRecSsc[i],nJetsBin); 
3183     NormalizeTH1(hdNdxiMCSecRecSsc[i],nJetsBin); 
3184
3185     
3186     // divide prim / (prim+sec) : corr factor
3187     TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3188     TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3189     TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3190     
3191     if(writeBgr){
3192       strNameSecCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3193       strNameSecCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3194       strNameSecCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3195     }
3196      
3197     hSecCorrPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt);
3198     TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt");
3199     hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS[i]);
3200     hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc[i]);
3201     hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
3202
3203     hSecCorrXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi);
3204     TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi");
3205     hSumPrimSecXi->Add(hdNdxiMCSecRecNS[i]);
3206     hSumPrimSecXi->Add(hdNdxiMCSecRecSsc[i]);
3207     hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
3208
3209     hSecCorrZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ);
3210     TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ");
3211     hSumPrimSecZ->Add(hdNdzMCSecRecNS[i]);
3212     hSumPrimSecZ->Add(hdNdzMCSecRecSsc[i]);
3213     hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
3214
3215     // the same using unscaled strangeness
3216     TString strNameSecCorrPt_nonSc(Form("hSecCorrPt_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3217     TString strNameSecCorrZ_nonSc(Form("hSecCorrZ_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3218     TString strNameSecCorrXi_nonSc(Form("hSecCorrXi_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3219
3220     if(writeBgr){
3221       strNameSecCorrPt_nonSc.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3222       strNameSecCorrZ_nonSc.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3223       strNameSecCorrXi_nonSc.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
3224     }
3225
3226     hSecCorrPt_nonSc[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt_nonSc);
3227     TH1F* hSumPrimSecPt_nonSc = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt_nonSc");
3228     hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecNS[i]);
3229     hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecS[i]); // non-scaled secondaries from strangeness
3230     hSecCorrPt_nonSc[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt_nonSc,1,1,"B"); // binominal errors
3231     
3232     hSecCorrZ_nonSc[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ_nonSc);
3233     TH1F* hSumPrimSecZ_nonSc = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ_nonSc");
3234     hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecNS[i]);
3235     hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecS[i]); // non-scaled secondaries from strangeness
3236     hSecCorrZ_nonSc[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ_nonSc,1,1,"B"); // binominal errors
3237     
3238     hSecCorrXi_nonSc[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi_nonSc);
3239     TH1F* hSumPrimSecXi_nonSc = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi_nonSc");
3240     hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecNS[i]);
3241     hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecS[i]); // non-scaled secondaries from strangeness
3242     hSecCorrXi_nonSc[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi_nonSc,1,1,"B"); // binominal errors
3243   } 
3244   
3245   // write 
3246
3247   TString outfileOption = "RECREATE";
3248   if(updateOutfile)  outfileOption = "UPDATE";
3249
3250   TFile out(strOutfile,outfileOption);
3251   
3252   if(!out.IsOpen()){
3253     Printf("%s:%d -- error opening sec corr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3254     return;
3255   }
3256   
3257   if(fDebug>0) Printf("%s:%d -- write jet sec corr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3258
3259   if(strOutDir && strOutDir.Length()){  
3260
3261     TDirectory* dir;
3262     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
3263     else{
3264       dir = out.mkdir(strOutDir);
3265       dir->cd(); 
3266     } 
3267   }
3268   
3269
3270   for(Int_t i=0; i<fNJetPtSlices; i++){
3271   
3272     hSecCorrPt[i]->Write();
3273     hSecCorrXi[i]->Write();
3274     hSecCorrZ[i]->Write();
3275   
3276     hSecCorrPt_nonSc[i]->Write();
3277     hSecCorrXi_nonSc[i]->Write();