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