a5c5f5e0957ccc6c6119262206ef4eaa62d7cd90
[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     hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1499     
1500     if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1501     if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1502     if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1503
1504     // backfolding: apply response matrix to unfolded spectrum
1505     TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); 
1506     if(hist) 
1507       hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1508
1509     if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1510     if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
1511     if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
1512     
1513     // ratio unfolded to original histo 
1514     TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1515     hRatioUnfolded->Reset();
1516     if (hist) 
1517       hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1518
1519     if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1520     if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
1521     if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;
1522
1523
1524     // ratio backfolded to original histo
1525     TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1526     hRatioBackFolded->Reset();
1527     hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1528
1529     if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1530     if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
1531     if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
1532     
1533     delete hnHist;
1534     delete hnFlatEfficiency;
1535
1536  }
1537 }
1538
1539 //_____________________________________________________________________________________________________
1540 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1541 {
1542   
1543   Int_t type = kFlagPt;
1544   UnfoldHistos(nIter, useCorrelatedErrors, type);
1545 }
1546
1547 //_____________________________________________________________________________________________________
1548 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1549 {
1550   
1551   Int_t type = kFlagZ;
1552   UnfoldHistos(nIter, useCorrelatedErrors, type);
1553 }
1554
1555 //_____________________________________________________________________________________________________
1556 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1557 {
1558   
1559   Int_t type = kFlagXi;
1560   UnfoldHistos(nIter, useCorrelatedErrors, type);
1561 }
1562
1563 //______________________________________________________________________________________________
1564 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1565 {
1566   // apply (multiply) response matrix to hist 
1567
1568   TH1F* hOut = new TH1F(*hist);
1569   hOut->Reset();
1570
1571   const Int_t axisM = 0; 
1572   const Int_t axisT = 1;
1573  
1574   Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1575   Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1576
1577   // response matrix normalization
1578   // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1579   // so for unfolding proper we leave it to CORRFW ...
1580
1581   Double_t normFacResponse[nBinsT];
1582
1583   for(Int_t iT=1; iT<=nBinsT; iT++){
1584
1585     Double_t sumResp = 0;
1586     
1587     for(Int_t iM=1; iM<=nBinsM; iM++){
1588       
1589       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1590       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1591       
1592       Double_t binCoord[] = {coordM,coordT};
1593       
1594       Long64_t binIndex = hnResponse->GetBin(binCoord);
1595       
1596       Double_t resp = hnResponse->GetBinContent(binIndex); 
1597       
1598       sumResp += resp;
1599     }
1600     
1601     normFacResponse[iT] = 0;
1602     if(sumResp) normFacResponse[iT] = 1/sumResp;
1603   }
1604   
1605   
1606   
1607   for(Int_t iM=1; iM<=nBinsM; iM++){
1608     
1609     Double_t contM   = 0;
1610
1611     for(Int_t iT=1; iT<=nBinsT; iT++){
1612
1613       Double_t contT = hist->GetBinContent(iT);
1614       
1615       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1616       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1617
1618       Double_t binCoord[] = {coordM,coordT};
1619       
1620       Long64_t binIndex = hnResponse->GetBin(binCoord);
1621       
1622       Double_t resp = hnResponse->GetBinContent(binIndex); 
1623       
1624       contM   += resp*normFacResponse[iT]*contT; 
1625     }
1626
1627     hOut->SetBinContent(iM,contM);
1628   }
1629
1630   return hOut;
1631 }
1632
1633 //_______________________________________________________________________________________________________
1634 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1635 {
1636   // read reconstruction efficiency from file
1637   // argument strlist optional - read from directory strdir if not specified
1638
1639   // temporary histos to hold histos from file
1640   TH1F* hEffPt[fNJetPtSlices]; 
1641   TH1F* hEffZ[fNJetPtSlices];
1642   TH1F* hEffXi[fNJetPtSlices];
1643   
1644   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1645   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
1646   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1647
1648   TFile f(strfile,"READ");
1649
1650   if(!f.IsOpen()){
1651     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1652     return;
1653   }
1654
1655   if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1656  
1657   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1658
1659   TList* list = 0;
1660
1661   if(strlist && strlist.Length()){
1662    
1663     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1664       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1665       return;
1666     }
1667   }  
1668
1669   for(Int_t i=0; i<fNJetPtSlices; i++){
1670     
1671     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1672     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1673     
1674     TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1675     TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1676     TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1677      
1678    
1679     if(list){
1680       hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt); 
1681       hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ); 
1682       hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi); 
1683     }
1684     else{
1685       hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt); 
1686       hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ); 
1687       hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi); 
1688     }
1689     
1690     if(!hEffPt[i]){
1691       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1692     }
1693   
1694     if(!hEffZ[i]){
1695       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1696     }    
1697
1698     if(!hEffXi[i]){
1699       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1700     }
1701
1702
1703     if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1704     if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1705     if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1706
1707     if(hEffPt[i]) hEffPt[i]->SetDirectory(0); 
1708     if(hEffZ[i])  hEffZ[i]->SetDirectory(0); 
1709     if(hEffXi[i]) hEffXi[i]->SetDirectory(0); 
1710
1711   } // jet slices loop
1712
1713   f.Close();
1714
1715   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1716     if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1717     if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
1718     if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1719   }
1720 }
1721
1722 //___________________________________________________________________________________________________________
1723 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1724 {
1725   // read bgr eff from file
1726   // argument strlist optional - read from directory strdir if not specified
1727  
1728   TH1F* hEffPtBgr[fNJetPtSlices]; 
1729   TH1F* hEffZBgr [fNJetPtSlices];
1730   TH1F* hEffXiBgr[fNJetPtSlices];
1731
1732   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1733   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
1734   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1735
1736
1737   TFile f(strfile,"READ");
1738
1739   if(!f.IsOpen()){
1740     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1741     return;
1742   }
1743
1744   if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1745  
1746   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1747
1748   TList* list = 0;
1749
1750   if(strlist && strlist.Length()){
1751        
1752     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1753       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1754       return;
1755     }
1756   }  
1757
1758   for(Int_t i=0; i<fNJetPtSlices; i++){
1759     
1760     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1761     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1762     
1763     TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1764     TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1765     TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1766   
1767    
1768     if(list){
1769       hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr); 
1770       hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr); 
1771       hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr); 
1772     }
1773     else{
1774       hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr); 
1775       hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr); 
1776       hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr); 
1777     }
1778     
1779     if(!hEffPtBgr[i]){
1780       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1781     }
1782   
1783     if(!hEffZBgr[i]){
1784       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1785     }    
1786
1787     if(!hEffXiBgr[i]){
1788       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1789     }
1790
1791
1792     if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1793     if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1794     if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1795
1796     if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0); 
1797     if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0); 
1798     if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0); 
1799
1800   } // jet slices loop
1801
1802   f.Close();
1803
1804   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1805     if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1806     if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
1807     if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1808   }
1809 }
1810
1811 // ________________________________________________
1812 void AliFragmentationFunctionCorrections::EffCorr()
1813 {
1814   // apply efficiency correction
1815
1816   AddCorrectionLevel("eff");
1817
1818   for(Int_t i=0; i<fNJetPtSlices; i++){
1819
1820     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1821     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1822     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1823
1824     TString histNamePt = histPt->GetName();
1825     TString histNameZ  = histZ->GetName();
1826     TString histNameXi = histXi->GetName();
1827
1828     
1829     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1830     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1831     
1832     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1833     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1834     
1835     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1836     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1837     
1838     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1839   }
1840 }
1841
1842 //___________________________________________________
1843 void AliFragmentationFunctionCorrections::EffCorrBgr()
1844 {
1845   // apply efficiency correction to bgr distributions
1846
1847   AddCorrectionLevelBgr("eff");
1848
1849   Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1850
1851   for(Int_t i=0; i<fNJetPtSlices; i++){
1852
1853     TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1854     TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1855     TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1856     
1857     TString histNamePt = histPt->GetName();
1858     TString histNameZ  = histZ->GetName();
1859     TString histNameXi = histXi->GetName();
1860
1861     
1862     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1863     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1864     
1865     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1866     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1867     
1868     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1869     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1870     
1871     fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1872   }
1873 }
1874
1875 //______________________________________________________________________
1876 void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1877 {
1878   // re-evaluate jet energy after FF corrections from dN/dpt distribution
1879   // 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'
1880   // argument corrlevel: which level of correction to be corrected/shifted to 
1881
1882   if(corrLevel>=fNCorrectionLevels){ 
1883     Printf(" calc xi shift: corrLevel exceeded - do nothing");
1884     return;
1885   }
1886
1887   Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1888   Double_t* jetPtCorr   = new Double_t[fNJetPtSlices];
1889   Double_t* deltaXi     = new Double_t[fNJetPtSlices];
1890
1891   for(Int_t i=0; i<fNJetPtSlices; i++){
1892     
1893     TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1894     TH1F* histPt    = fCorrFF[corrLevel]->GetTrackPt(i);
1895
1896     Double_t ptUncorr = 0;
1897     Double_t ptCorr   = 0;
1898
1899     for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1900
1901       Double_t cont   = histPtRaw->GetBinContent(bin);
1902       Double_t width  = histPtRaw->GetBinWidth(bin);
1903       Double_t meanPt = histPtRaw->GetBinCenter(bin);
1904
1905       ptUncorr += meanPt*cont*width;
1906     }
1907     
1908     for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1909       
1910       Double_t cont   = histPt->GetBinContent(bin);
1911       Double_t width  = histPt->GetBinWidth(bin);
1912       Double_t meanPt = histPt->GetBinCenter(bin);
1913       
1914       ptCorr += meanPt*cont*width;
1915     }
1916
1917     jetPtUncorr[i] = ptUncorr; 
1918     jetPtCorr[i]   = ptCorr;   
1919   }
1920
1921   // calc dXi from dN/dpt distribution : 
1922   // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt 
1923
1924   for(Int_t i=0; i<fNJetPtSlices; i++){
1925
1926     Float_t jetPtLoLim = fJetPtSlices->At(i);
1927     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1928
1929     Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1930     
1931     Double_t ptUncorr = jetPtUncorr[i]; 
1932     Double_t ptCorr   = jetPtCorr[i]; 
1933
1934     Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1935     
1936     Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1937            ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1938     
1939     deltaXi[i] = dXi;
1940   }
1941   
1942   // book & fill new dN/dxi histos
1943
1944   for(Int_t i=0; i<fNJetPtSlices; i++){
1945
1946     TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1947             
1948     Double_t dXi = deltaXi[i];
1949
1950     Int_t nBins  = histXi->GetNbinsX();
1951     const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1952     Float_t binsVecNew[nBins+1];
1953     
1954     TString strName = histXi->GetName(); 
1955     strName.Append("_shift");
1956     TString strTit  = histXi->GetTitle();  
1957
1958     TH1F* hXiCorr; 
1959
1960     // create shifted histo ...
1961
1962     if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1963
1964       for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;    
1965       hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1966     }
1967     else{ // uniform bin size
1968       
1969       Double_t xMin  = histXi->GetXaxis()->GetXmin();
1970       Double_t xMax  = histXi->GetXaxis()->GetXmax();
1971
1972       xMin += dXi;
1973       xMax += dXi;
1974       
1975       hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1976     }
1977
1978     // ... and fill
1979
1980     for(Int_t bin=1; bin<nBins+1; bin++){
1981       Double_t cont = histXi->GetBinContent(bin);
1982       Double_t err  = histXi->GetBinError(bin);
1983       
1984       hXiCorr->SetBinContent(bin,cont);
1985       hXiCorr->SetBinError(bin,err);
1986     }
1987     
1988     new(fh1FFXiShift[i]) TH1F(*hXiCorr);
1989     delete hXiCorr;
1990   }
1991
1992   delete[] jetPtUncorr;
1993   delete[] jetPtCorr;
1994   delete[] deltaXi;
1995 }
1996
1997 //_____________________________________________________
1998 void AliFragmentationFunctionCorrections::SubtractBgr()
1999 {
2000   // subtract bgr distribution from FF
2001   // the current corr level is used for both 
2002   
2003   AddCorrectionLevel("bgrSub");
2004
2005   for(Int_t i=0; i<fNJetPtSlices; i++){
2006
2007     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2008     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2009     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2010     
2011     TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2012     TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2013     TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2014
2015     TString histNamePt = histPt->GetName();
2016     TString histNameZ  = histZ->GetName();
2017     TString histNameXi = histXi->GetName();
2018     
2019
2020     TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2021     hFFTrackPtBgrSub->Add(histPtBgr,-1);
2022     
2023     TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
2024     hFFZBgrSub->Add(histZBgr,-1);
2025     
2026     TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2027     hFFXiBgrSub->Add(histXiBgr,-1);
2028
2029     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2030   }
2031 }
2032
2033 //________________________________________________________________________________________________________________
2034 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile, 
2035                                                               Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2036
2037   // read task ouput from MC and write single track eff - standard dir/list 
2038      
2039   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2040   TString strlist = "fracfunc_" + strID;
2041     
2042   WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2043 }
2044
2045 //___________________________________________________________________________________________________________________________________
2046 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, 
2047                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2048 {
2049   // read task output from MC and write single track eff as function of pt, eta, phi
2050   // argument strlist optional - read from directory strdir if not specified
2051   // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2052
2053
2054   TH1D* hdNdptTracksMCPrimGen;
2055   TH2D* hdNdetadphiTracksMCPrimGen;
2056   
2057   TH1D* hdNdptTracksMCPrimRec;
2058   TH2D* hdNdetadphiTracksMCPrimRec;
2059     
2060
2061   TFile f(strInfile,"READ");
2062
2063   if(!f.IsOpen()){
2064     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2065     return;
2066   }
2067   
2068   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2069  
2070   if(strdir && strdir.Length()){
2071     if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2072     gDirectory->cd(strdir);
2073   }
2074
2075   TList* list = 0;
2076
2077   if(strlist && strlist.Length()){
2078
2079     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2080       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2081       return;
2082     }
2083   }
2084
2085
2086   TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2087   if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data()); 
2088
2089   TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2090   if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2091
2092   TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2093   if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2094   
2095   TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2096   if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2097   
2098
2099   if(list){
2100     hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
2101     hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2102     
2103     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
2104     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2105   }
2106   else{
2107     hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2108     hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2109
2110     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2111     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2112   }
2113
2114   hdNdptTracksMCPrimGen->SetDirectory(0);
2115   hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2116   hdNdptTracksMCPrimRec->SetDirectory(0);
2117   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2118   
2119   f.Close();
2120
2121   // projections: dN/deta, dN/dphi 
2122
2123   TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2124   TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2125  
2126   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2127   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2128
2129   // rebin
2130
2131   TString strNamePtGen = "hTrackPtGenPrim";
2132   TString strNamePtRec = "hTrackPtRecPrim";
2133
2134   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2135   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2136  
2137     // efficiency: divide 
2138
2139   TString hNameTrackEffPt = "hSingleTrackEffPt";
2140   if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2141                                                
2142   TString hNameTrackEffEta = "hSingleTrackEffEta";
2143   if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2144
2145   TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2146   if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2147
2148
2149   TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2150   hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2151
2152   TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2153   hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2154   
2155   TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2156   hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2157   
2158   
2159   TString outfileOption = "RECREATE";
2160   if(updateOutfile)  outfileOption = "UPDATE";
2161
2162   TFile out(strOutfile,outfileOption);
2163
2164   if(!out.IsOpen()){
2165     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2166     return;
2167   }
2168
2169   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2170
2171   if(strOutDir && strOutDir.Length()){
2172     
2173     TDirectory* dir;
2174     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2175     else{
2176       dir = out.mkdir(strOutDir);
2177       dir->cd(); 
2178     } 
2179   }
2180
2181   hSingleTrackEffPt->Write();
2182   hSingleTrackEffEta->Write();
2183   hSingleTrackEffPhi->Write();
2184   
2185   out.Close();
2186 }
2187
2188 //________________________________________________________________________________________________________________
2189 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile, 
2190                                                                   Bool_t updateOutfile, TString strOutDir)
2191
2192   // read task ouput from MC and write single track eff - standard dir/list 
2193      
2194   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2195   TString strlist = "fracfunc_" + strID;
2196     
2197   WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2198 }
2199
2200 //___________________________________________________________________________________________________________________________________
2201 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, 
2202                                                                   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2203 {
2204   // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2205   // argument strlist optional - read from directory strdir if not specified
2206   // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2207
2208   TH1D* hdNdptTracksMCPrimRec;
2209   TH2D* hdNdetadphiTracksMCPrimRec;
2210   
2211   TH1D* hdNdptTracksMCSecRec;
2212   TH2D* hdNdetadphiTracksMCSecRec;
2213
2214   TFile f(strInfile,"READ");
2215
2216   if(!f.IsOpen()){
2217     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2218     return;
2219   }
2220   
2221   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2222  
2223   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2224
2225   TList* list = 0;
2226
2227   if(strlist && strlist.Length()){
2228
2229     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2230       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2231       return;
2232     }
2233   }
2234
2235
2236   if(list){
2237     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2238     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2239     
2240     hdNdptTracksMCSecRec       = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2241     hdNdetadphiTracksMCSecRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2242   }
2243   else{
2244     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2245     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2246
2247     hdNdptTracksMCSecRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2248     hdNdetadphiTracksMCSecRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2249   }
2250   
2251   hdNdptTracksMCPrimRec->SetDirectory(0);
2252   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2253   hdNdptTracksMCSecRec->SetDirectory(0);
2254   hdNdetadphiTracksMCSecRec->SetDirectory(0);
2255   
2256   f.Close();
2257
2258   // projections: dN/deta, dN/dphi 
2259
2260   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2261   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2262  
2263   TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2264   TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2265
2266  
2267   // rebin
2268
2269   TString strNamePtPrim = "hTrackPtPrim";
2270   TString strNamePtSec  = "hTrackPtSec";
2271
2272   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2273   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec  = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2274
2275
2276   // secondary correction factor: divide prim/(prim+sec)
2277
2278   TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2279   TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2280   hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2281   hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2282
2283   TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2284   TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2285   hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2286   hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2287
2288   TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2289   TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2290   hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2291   hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2292   
2293   
2294   TString outfileOption = "RECREATE";
2295   if(updateOutfile)  outfileOption = "UPDATE";
2296
2297   TFile out(strOutfile,outfileOption);
2298
2299   if(!out.IsOpen()){
2300     Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2301     return;
2302   }
2303
2304   if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2305
2306   if(strOutDir && strOutDir.Length()){  
2307
2308     TDirectory* dir;
2309     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2310     else{
2311       dir = out.mkdir(strOutDir);
2312       dir->cd(); 
2313     } 
2314   } 
2315
2316   hdNdptTracksMCSecRec->Write();
2317   hdNdetaTracksMCSecRec->Write();
2318   hdNdphiTracksMCSecRec->Write();
2319
2320   hSingleTrackSecCorrPt->Write();
2321   hSingleTrackSecCorrEta->Write();
2322   hSingleTrackSecCorrPhi->Write();
2323   
2324   out.Close();
2325 }
2326
2327 //________________________________________________________________________________________________________________
2328 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile, 
2329                                                               Bool_t updateOutfile, TString strOutDir)
2330
2331   // read task ouput from MC and write single track eff - standard dir/list 
2332      
2333   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2334   TString strlist = "fracfunc_" + strID;
2335     
2336   WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2337 }
2338
2339
2340 //_____________________________________________________________________________________________________________________________________
2341 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2342                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2343 {
2344   // read 2d THnSparse response matrices in pt from file
2345   // project TH2 
2346   // write to strOutfile 
2347
2348   THnSparse* hnResponseSinglePt; 
2349   TH2F*      h2ResponseSinglePt; 
2350  
2351   TFile f(strInfile,"READ");
2352
2353   if(!f.IsOpen()){
2354     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2355     return;
2356   }
2357
2358   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2359  
2360   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2361
2362   TList* list = 0;
2363
2364   if(strlist && strlist.Length()){
2365     
2366     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2367       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2368       return;
2369     }
2370   }  
2371   
2372   if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2373   else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2374   
2375
2376   if(!hnResponseSinglePt){
2377     Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2378     return;
2379   }
2380
2381   f.Close();
2382
2383
2384   // project 2d histo 
2385   h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2386   h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt",""); 
2387     
2388   
2389   // write 
2390
2391   TString outfileOption = "RECREATE";
2392   if(updateOutfile)  outfileOption = "UPDATE";
2393   
2394   TFile out(strOutfile,outfileOption);
2395   
2396   if(!out.IsOpen()){
2397     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2398     return;
2399   }
2400
2401   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2402   
2403   if(strOutDir && strOutDir.Length()){  
2404
2405     TDirectory* dir;
2406     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
2407     else{
2408       dir = out.mkdir(strOutDir);
2409       dir->cd(); 
2410     } 
2411   }
2412   
2413   hnResponseSinglePt->Write();
2414   h2ResponseSinglePt->Write();
2415   
2416   out.Close();  
2417 }
2418
2419 //________________________________________________________________________________________________________________
2420 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile, 
2421                                                            Bool_t updateOutfile)
2422
2423   // read task ouput from MC and write single track eff - standard dir/list 
2424      
2425   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2426   TString strlist = "fracfunc_" + strID;
2427     
2428   WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2429 }
2430
2431 //___________________________________________________________________________________________________________________________________
2432 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, 
2433                                                            TString strOutfile, Bool_t updateOutfile)
2434 {
2435   // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2436   // argument strlist optional - read from directory strdir if not specified
2437   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2438
2439   TH1F* hEffPt[fNJetPtSlices];
2440   TH1F* hEffXi[fNJetPtSlices];
2441   TH1F* hEffZ[fNJetPtSlices];
2442
2443   TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2444   TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2445   TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2446     
2447   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2448   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2449   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2450
2451
2452   TH1F* fh1FFJetPtRecEffGen;
2453
2454   TH2F* fh2FFTrackPtRecEffGen;
2455   TH2F* fh2FFZRecEffGen;
2456   TH2F* fh2FFXiRecEffGen;
2457   
2458   TH2F* fh2FFTrackPtRecEffRec;
2459   TH2F* fh2FFZRecEffRec;
2460   TH2F* fh2FFXiRecEffRec;
2461  
2462
2463   TFile f(strInfile,"READ");
2464
2465   if(!f.IsOpen()){
2466     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2467     return;
2468   }
2469   
2470   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2471  
2472   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2473
2474   TList* list = 0;
2475
2476   if(strlist && strlist.Length()){
2477
2478     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2479       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2480       return;
2481     }
2482   }
2483
2484   if(list){
2485     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2486
2487     fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2488     fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2489     fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2490     
2491     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2492     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2493     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2494   }
2495   else{
2496     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2497
2498     fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2499     fh2FFZRecEffGen        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2500     fh2FFXiRecEffGen       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2501     
2502     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2503     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2504     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2505   }
2506   
2507   fh1FFJetPtRecEffGen->SetDirectory(0); 
2508
2509   fh2FFTrackPtRecEffGen->SetDirectory(0);
2510   fh2FFZRecEffGen->SetDirectory(0);
2511   fh2FFXiRecEffGen->SetDirectory(0);
2512   
2513   fh2FFTrackPtRecEffRec->SetDirectory(0);
2514   fh2FFZRecEffRec->SetDirectory(0);
2515   fh2FFXiRecEffRec->SetDirectory(0);
2516   
2517   f.Close();
2518
2519
2520   // projections: FF for generated and reconstructed primaries 
2521   
2522   for(Int_t i=0; i<fNJetPtSlices; i++){
2523     
2524     Float_t jetPtLoLim = fJetPtSlices->At(i);
2525     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2526
2527     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2528     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2529
2530     if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2531       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2532       return; 
2533     }
2534     
2535     TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2536     TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2537     TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2538     
2539     TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2540     TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2541     TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2542     
2543     // project 
2544     // appendix 'unbinned' to avoid histos with same name after rebinning
2545
2546     hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2547     hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2548     hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2549     
2550     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2551     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2552     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2553     
2554     // rebin
2555
2556     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2557     if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2558     if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2559
2560     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2561     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2562     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2563
2564     hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2565     hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2566     hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2567     
2568     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2569     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2570     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2571  
2572     // normalize
2573     
2574     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2575
2576     NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin); 
2577     NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin); 
2578     NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin); 
2579
2580     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2581     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2582     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2583     
2584     // divide rec/gen : efficiency
2585
2586     TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2587     TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2588     TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2589  
2590     hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2591     hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2592     
2593     hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2594     hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2595     
2596     hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2597     hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2598   } 
2599   
2600   // write 
2601
2602   TString outfileOption = "RECREATE";
2603   if(updateOutfile)  outfileOption = "UPDATE";
2604
2605   TFile out(strOutfile,outfileOption);
2606   
2607   if(!out.IsOpen()){
2608     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2609     return;
2610   }
2611
2612   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2613
2614 //   if(strdir && strdir.Length()){ 
2615 //     TDirectory* dir = out.mkdir(strdir);
2616 //     dir->cd(); 
2617 //   }
2618
2619   for(Int_t i=0; i<fNJetPtSlices; i++){
2620
2621     hdNdptTracksMCPrimGen[i]->Write();
2622     hdNdxiMCPrimGen[i]->Write();
2623     hdNdzMCPrimGen[i]->Write();
2624     
2625     hdNdptTracksMCPrimRec[i]->Write();
2626     hdNdxiMCPrimRec[i]->Write();
2627     hdNdzMCPrimRec[i]->Write();
2628   
2629     hEffPt[i]->Write();
2630     hEffXi[i]->Write();
2631     hEffZ[i]->Write();
2632   }
2633
2634   out.Close();
2635
2636 }
2637
2638 //________________________________________________________________________________________________________________
2639 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile, 
2640                                                            Bool_t updateOutfile)
2641
2642   // read task ouput from MC and write secondary correction - standard dir/list 
2643      
2644   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2645   TString strlist = "fracfunc_" + strID;
2646     
2647   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2648 }
2649
2650 //___________________________________________________________________________________________________________________________________
2651 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, 
2652                                                            TString strOutfile, Bool_t updateOutfile)
2653 {
2654   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2655   // argument strlist optional - read from directory strdir if not specified
2656   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2657
2658   TH1F* hSecCorrPt[fNJetPtSlices];
2659   TH1F* hSecCorrXi[fNJetPtSlices];
2660   TH1F* hSecCorrZ[fNJetPtSlices];
2661
2662   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2663   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2664   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2665     
2666   TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2667   TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2668   TH1F* hdNdzMCSecRec[fNJetPtSlices];
2669
2670   TH1F* fh1FFJetPtRecEffGen;
2671
2672   TH2F* fh2FFTrackPtRecEffRec;
2673   TH2F* fh2FFZRecEffRec;
2674   TH2F* fh2FFXiRecEffRec;
2675   
2676   TH2F* fh2FFTrackPtSecRec;
2677   TH2F* fh2FFZSecRec;
2678   TH2F* fh2FFXiSecRec;
2679
2680   TFile f(strInfile,"READ");
2681
2682   if(!f.IsOpen()){
2683     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2684     return;
2685   }
2686   
2687   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2688  
2689   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2690
2691   TList* list = 0;
2692
2693   if(strlist && strlist.Length()){
2694
2695     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2696       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2697       return;
2698     }
2699   }
2700
2701   if(list){
2702     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2703
2704     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2705     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2706     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2707     
2708     fh2FFTrackPtSecRec  = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2709     fh2FFZSecRec        = (TH2F*) list->FindObject("fh2FFZSecRec");
2710     fh2FFXiSecRec       = (TH2F*) list->FindObject("fh2FFXiSecRec");
2711   }
2712   else{
2713     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2714
2715     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2716     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2717     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2718     
2719     fh2FFTrackPtSecRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2720     fh2FFZSecRec        = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2721     fh2FFXiSecRec       = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2722   }
2723   
2724   fh1FFJetPtRecEffGen->SetDirectory(0); 
2725
2726   fh2FFTrackPtRecEffRec->SetDirectory(0);
2727   fh2FFZRecEffRec->SetDirectory(0);
2728   fh2FFXiRecEffRec->SetDirectory(0);
2729   
2730   fh2FFTrackPtSecRec->SetDirectory(0);
2731   fh2FFZSecRec->SetDirectory(0);
2732   fh2FFXiSecRec->SetDirectory(0);
2733   
2734   f.Close();
2735
2736
2737   // projections: FF for generated and reconstructed primaries 
2738   
2739   for(Int_t i=0; i<fNJetPtSlices; i++){
2740     
2741     Float_t jetPtLoLim = fJetPtSlices->At(i);
2742     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2743
2744     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2745     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2746
2747     if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2748       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2749       return; 
2750     }
2751     
2752     TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2753     TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2754     TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2755     
2756     TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2757     TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2758     TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2759     
2760     // project 
2761     // appendix 'unbinned' to avoid histos with same name after rebinning
2762
2763     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2764     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2765     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2766     
2767     hdNdptTracksMCSecRec[i]  = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2768     hdNdzMCSecRec[i]         = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2769     hdNdxiMCSecRec[i]        = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2770     
2771     // rebin
2772
2773     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2774     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2775     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2776
2777     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2778     if(fNHistoBinsZ[i])  hdNdzMCSecRec[i]        = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2779     if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i]       = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2780
2781     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2782     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2783     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2784     
2785     hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2786     hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2787     hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2788     
2789     // normalize
2790     
2791     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2792
2793     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2794     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2795     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2796
2797     NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin); 
2798     NormalizeTH1(hdNdzMCSecRec[i],nJetsBin); 
2799     NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin); 
2800     
2801     // divide rec/gen : efficiency
2802
2803     TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2804     TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2805     TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2806  
2807     hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2808     TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2809     hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2810     hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2811
2812     hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2813     TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2814     hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2815     hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2816
2817     hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2818     TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2819     hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2820     hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2821   } 
2822   
2823   // write 
2824
2825   TString outfileOption = "RECREATE";
2826   if(updateOutfile)  outfileOption = "UPDATE";
2827
2828   TFile out(strOutfile,outfileOption);
2829   
2830   if(!out.IsOpen()){
2831     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2832     return;
2833   }
2834   
2835   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2836
2837   for(Int_t i=0; i<fNJetPtSlices; i++){
2838
2839     //  hdNdptTracksMCSecRec[i]->Write();
2840     //  hdNdxiMCSecRec[i]->Write();
2841     //  hdNdzMCSecRec[i]->Write();
2842   
2843     hSecCorrPt[i]->Write();
2844     hSecCorrXi[i]->Write();
2845     hSecCorrZ[i]->Write();
2846   }
2847
2848   out.Close();
2849 }
2850
2851 //________________________________________________________________________________________________________________
2852 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile, 
2853                                                            Bool_t updateOutfile)
2854
2855   // read task ouput from MC and write single track eff - standard dir/list 
2856      
2857   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2858   TString strlist = "fracfunc_" + strID;
2859     
2860   WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2861 }
2862
2863 //_____________________________________________________________________________________________________________________________________
2864 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2865                                                            TString strOutfile, Bool_t updateOutfile)
2866 {
2867   // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2868   // project THnSparse + TH2 in jet pt slices 
2869   // write to strOutfile 
2870
2871   THnSparse* hn3ResponseJetPt;
2872   THnSparse* hn3ResponseJetZ;
2873   THnSparse* hn3ResponseJetXi;
2874
2875   // 2D response matrices 
2876
2877   THnSparse* hnResponsePt[fNJetPtSlices];
2878   THnSparse* hnResponseZ[fNJetPtSlices];
2879   THnSparse* hnResponseXi[fNJetPtSlices];
2880
2881   TH2F* h2ResponsePt[fNJetPtSlices];
2882   TH2F* h2ResponseZ[fNJetPtSlices];
2883   TH2F* h2ResponseXi[fNJetPtSlices];
2884
2885   // 1D projections on gen pt / rec pt axes
2886
2887   TH1F* h1FFPtRec[fNJetPtSlices]; 
2888   TH1F* h1FFZRec[fNJetPtSlices];
2889   TH1F* h1FFXiRec[fNJetPtSlices];
2890
2891   TH1F* h1FFPtGen[fNJetPtSlices]; 
2892   TH1F* h1FFZGen[fNJetPtSlices];
2893   TH1F* h1FFXiGen[fNJetPtSlices];
2894
2895   TH1F* h1RatioPt[fNJetPtSlices]; 
2896   TH1F* h1RatioZ[fNJetPtSlices];
2897   TH1F* h1RatioXi[fNJetPtSlices];
2898
2899
2900
2901   TFile f(strInfile,"READ");
2902
2903   if(!f.IsOpen()){
2904     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2905     return;
2906   }
2907
2908   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2909  
2910   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2911
2912   TList* list = 0;
2913
2914   if(strlist && strlist.Length()){
2915     
2916     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2917       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2918       return;
2919     }
2920   }  
2921   
2922   if(list){
2923     hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2924     hn3ResponseJetZ  = (THnSparse*) list->FindObject("fhnResponseJetZ");
2925     hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2926   }
2927   else{
2928     hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2929     hn3ResponseJetZ  = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2930     hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2931   }
2932
2933   
2934   if(!hn3ResponseJetPt){
2935     Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2936     return;
2937   }
2938
2939   if(!hn3ResponseJetZ){
2940     Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2941     return;
2942   }
2943
2944   if(!hn3ResponseJetXi){
2945     Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2946     return;
2947   }
2948
2949   f.Close();  
2950
2951   // axes 
2952
2953   Int_t axisJetPtTHn3 = -1;
2954   Int_t axisGenPtTHn3 = -1;
2955   Int_t axisRecPtTHn3 = -1;
2956
2957   for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2958     
2959     TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle(); 
2960
2961     if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i; 
2962     if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i; 
2963     if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i; 
2964   }
2965
2966   if(axisJetPtTHn3 == -1){
2967     Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2968     return;
2969   }
2970
2971   if(axisGenPtTHn3 == -1){
2972     Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2973     return;
2974   }
2975
2976   if(axisRecPtTHn3 == -1){
2977     Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2978     return;
2979   }
2980
2981   for(Int_t i=0; i<fNJetPtSlices; i++){
2982     
2983     Float_t jetPtLoLim = fJetPtSlices->At(i);
2984     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2985
2986     Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2987     Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2988
2989     if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
2990       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2991       return; 
2992     }
2993     
2994     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2995     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2996     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2997
2998     TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2999     TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3000     TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3001          
3002     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3003     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3004     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3005  
3006     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3007     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3008     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3009  
3010     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3011     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3012     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3013  
3014     
3015     // 2D projections in jet pt range
3016   
3017     hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
3018     hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
3019     hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
3020
3021     Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};  
3022
3023     hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
3024     hnResponseZ[i]  = hn3ResponseJetZ->Projection(2,axesProj);
3025     hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
3026   
3027     hnResponsePt[i]->SetNameTitle(strNameRespPt,""); 
3028     hnResponseZ[i]->SetNameTitle(strNameRespZ,""); 
3029     hnResponseXi[i]->SetNameTitle(strNameRespXi,""); 
3030
3031     h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3032     h2ResponseZ[i]  = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3033     h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3034  
3035     h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,""); 
3036     h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,""); 
3037     h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,""); 
3038
3039
3040     // 1D projections
3041
3042     Int_t axisGenPtTHn2 = -1;
3043     Int_t axisRecPtTHn2 = -1;
3044
3045     for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3046     
3047       TString title = hnResponsePt[i]->GetAxis(d)->GetTitle(); 
3048
3049       if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d; 
3050       if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d; 
3051     }
3052
3053     
3054     if(axisGenPtTHn2 == -1){
3055       Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3056       return;
3057     }
3058     
3059     if(axisRecPtTHn2 == -1){
3060       Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3061       return;
3062     }
3063     
3064
3065     h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3066     h1FFZRec[i]  = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3067     h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3068     
3069     h1FFPtRec[i]->SetNameTitle(strNameRecPt,""); 
3070     h1FFZRec[i]->SetNameTitle(strNameRecZ,""); 
3071     h1FFXiRec[i]->SetNameTitle(strNameRecXi,""); 
3072
3073     h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3074     h1FFZGen[i]  = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3075     h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3076     
3077     h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3078     h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3079     h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3080
3081     // normalize 1D projections
3082
3083     if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
3084     if(fNHistoBinsZ[i])  h1FFZRec[i]  = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3085     if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3086     
3087     if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3088     if(fNHistoBinsZ[i])  h1FFZGen[i]  = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3089     if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3090        
3091     NormalizeTH1(h1FFPtRec[i],fNJets->At(i)); 
3092     NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3093     NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3094
3095     NormalizeTH1(h1FFPtGen[i],fNJets->At(i)); 
3096     NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3097     NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3098
3099     // ratios 1D projections
3100
3101     h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3102     h1RatioPt[i]->Reset();
3103     h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3104
3105     h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3106     h1RatioZ[i]->Reset();
3107     h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3108     
3109     h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3110     h1RatioXi[i]->Reset();
3111     h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3112   }
3113   
3114   
3115   // write 
3116
3117   TString outfileOption = "RECREATE";
3118   if(updateOutfile)  outfileOption = "UPDATE";
3119   
3120   TFile out(strOutfile,outfileOption);
3121   
3122   if(!out.IsOpen()){
3123     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3124     return;
3125   }
3126
3127   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3128
3129   //if(strdir && strdir.Length()){ 
3130   //  TDirectory* dir = out.mkdir(strdir);
3131   //  dir->cd(); 
3132   //}
3133
3134   for(Int_t i=0; i<fNJetPtSlices; i++){
3135
3136     hnResponsePt[i]->Write();
3137     hnResponseXi[i]->Write();
3138     hnResponseZ[i]->Write();
3139
3140     h2ResponsePt[i]->Write();
3141     h2ResponseXi[i]->Write();
3142     h2ResponseZ[i]->Write();
3143
3144     h1FFPtRec[i]->Write(); 
3145     h1FFZRec[i]->Write();
3146     h1FFXiRec[i]->Write();
3147     
3148     h1FFPtGen[i]->Write(); 
3149     h1FFZGen[i]->Write();
3150     h1FFXiGen[i]->Write();
3151     
3152     h1RatioPt[i]->Write(); 
3153     h1RatioZ[i]->Write();
3154     h1RatioXi[i]->Write();
3155
3156   }
3157
3158   out.Close();  
3159 }
3160
3161 //______________________________________________________________________________________________________
3162 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3163 {
3164   // read response matrices from file
3165   // argument strlist optional - read from directory strdir if not specified
3166   // note: THnSparse are not rebinned 
3167
3168   THnSparse* hRespPt[fNJetPtSlices]; 
3169   THnSparse* hRespZ[fNJetPtSlices];
3170   THnSparse* hRespXi[fNJetPtSlices];
3171   
3172   for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3173   for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i]  = 0;
3174   for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3175
3176   TFile f(strfile,"READ");
3177
3178   if(!f.IsOpen()){
3179     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3180     return;
3181   }
3182
3183   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3184  
3185   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3186
3187   TList* list = 0;
3188
3189   if(strlist && strlist.Length()){
3190    
3191     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3192       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3193       return;
3194     }
3195   }  
3196
3197   for(Int_t i=0; i<fNJetPtSlices; i++){
3198     
3199     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3200     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3201     
3202     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3203     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3204     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3205         
3206     if(list){
3207       hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt); 
3208       hRespZ[i]  = (THnSparse*) list->FindObject(strNameRespZ); 
3209       hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi); 
3210     }
3211     else{
3212       hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt); 
3213       hRespZ[i]  = (THnSparse*) gDirectory->Get(strNameRespZ); 
3214       hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi); 
3215     }
3216     
3217     if(!hRespPt[i]){
3218       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3219     }
3220   
3221     if(!hRespZ[i]){
3222       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3223     }    
3224
3225     if(!hRespXi[i]){
3226       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3227     }
3228    
3229     //    if(0){ // can't rebin THnSparse ...
3230     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3231     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3232     
3233     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3234     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3235     
3236     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3237     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3238     //     }
3239    
3240  
3241   } // jet slices loop
3242
3243   f.Close();  // THnSparse pointers still valid even if file closed
3244
3245 //   for(Int_t i=0; i<fNJetPtSlices; i++){  // no copy c'tor ...
3246 //     if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3247 //     if(hRespZ[i])  new(fhnResponseZ[i])  THnSparseF(*hRespZ[i]);
3248 //     if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3249 //   }
3250
3251   for(Int_t i=0; i<fNJetPtSlices; i++){ 
3252     fhnResponsePt[i] = hRespPt[i];
3253     fhnResponseZ[i]  = hRespZ[i];
3254     fhnResponseXi[i] = hRespXi[i];
3255   }
3256 }
3257
3258 //______________________________________________________________________________________________________________________
3259 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3260 {
3261   // read priors from file: rec primaries, gen pt dist
3262
3263   if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3264   
3265   // temporary histos to store pointers from file
3266   TH1F* hist[fNJetPtSlices]; 
3267
3268   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3269
3270   TFile f(strfile,"READ");
3271
3272   if(!f.IsOpen()){
3273     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3274     return;
3275   }
3276
3277   for(Int_t i=0; i<fNJetPtSlices; i++){
3278     
3279     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3280     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3281     
3282     TString strName;
3283     
3284     if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); 
3285     if(type == kFlagZ)  strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);       
3286     if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);      
3287     
3288     hist[i] = (TH1F*) gDirectory->Get(strName); 
3289     
3290     if(!hist[i]){
3291       Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3292     }
3293     
3294
3295     //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3296
3297     if(hist[i]) hist[i]->SetDirectory(0); 
3298     
3299   } // jet slices loop
3300   
3301   f.Close();
3302
3303   
3304   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3305     if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]); 
3306     if(hist[i] && type == kFlagZ)  new(fh1FFZPrior[i])       TH1F(*hist[i]); 
3307     if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i])      TH1F(*hist[i]); 
3308   }
3309 }
3310
3311 //_____________________________________________________
3312 // void AliFragmentationFunctionCorrections::RatioRecGen()
3313 // {
3314 //   // create ratio reconstructed over generated FF
3315 //   // use current highest corrLevel 
3316   
3317 //   Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3318
3319 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3320  
3321 //     TH1F* histPtRec =  fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3322 //     TH1F* histZRec  =  fCorrFF[fNCorrectionLevels-1]->GetZ(i);       // levels -1: latest corr level
3323 //     TH1F* histXiRec =  fCorrFF[fNCorrectionLevels-1]->GetXi(i);      // levels -1: latest corr level
3324     
3325 //     TH1F* histPtGen = fh1FFTrackPtGenPrim[i]; 
3326 //     TH1F* histZGen  = fh1FFZGenPrim[i];
3327 //     TH1F* histXiGen = fh1FFXiGenPrim[i];
3328
3329 //     TString histNamePt = histPtRec->GetName();
3330 //     TString histNameZ  = histZRec->GetName();
3331 //     TString histNameXi = histXiRec->GetName();
3332
3333 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3334 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3335 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3336
3337 //     // ratio 
3338 //     TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3339 //     hRatioRecGenPt->Reset();
3340 //     hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3341
3342 //     TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3343 //     hRatioRecGenZ->Reset();
3344 //     hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3345     
3346 //     TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3347 //     hRatioRecGenXi->Reset();
3348 //     hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3349
3350 //     new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3351 //     new(fh1FFRatioRecGenZ[i])  TH1F(*hRatioRecGenZ);
3352 //     new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3353 //   }
3354 // }
3355
3356 // //___________________________________________________________
3357 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3358 // {
3359 //   // create ratio reconstructed tracks over reconstructed primaries 
3360 //   // use raw FF (corrLevel 0)
3361   
3362 //   Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3363
3364 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3365  
3366 //     const Int_t corrLevel = 0; 
3367
3368 //     TH1F* histPtRec =  fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3369 //     TH1F* histZRec  =  fCorrFF[corrLevel]->GetZ(i);       // levels -1: latest corr level
3370 //     TH1F* histXiRec =  fCorrFF[corrLevel]->GetXi(i);      // levels -1: latest corr level
3371     
3372 //     TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i]; 
3373 //     TH1F* histZRecPrim  = fh1FFZRecPrim[i];
3374 //     TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3375
3376 //     TString histNamePt = histPtRec->GetName();
3377 //     TString histNameZ  = histZRec->GetName();
3378 //     TString histNameXi = histXiRec->GetName();
3379
3380 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3381 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3382 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3383
3384 //     // ratio 
3385 //     TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3386 //     hRatioRecPrimPt->Reset();
3387 //     hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3388
3389 //     TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3390 //     hRatioRecPrimZ->Reset();
3391 //     hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3392     
3393 //     TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3394 //     hRatioRecPrimXi->Reset();
3395 //     hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3396
3397
3398 //     new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3399 //     new(fh1FFRatioRecPrimZ[i])  TH1F(*hRatioRecPrimZ);
3400 //     new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3401 //   }
3402 // }
3403
3404 // __________________________________________________________________________________
3405 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3406 {
3407
3408   // project response matrices on both axes: 
3409   // FF for rec primaries, in terms of generated and reconstructed momentum
3410   // write FF and ratios to outFile
3411
3412   Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3413
3414   TH1F* hFFPtRec[fNJetPtSlices]; 
3415   TH1F* hFFZRec[fNJetPtSlices];
3416   TH1F* hFFXiRec[fNJetPtSlices];
3417
3418   TH1F* hFFPtGen[fNJetPtSlices]; 
3419   TH1F* hFFZGen[fNJetPtSlices];
3420   TH1F* hFFXiGen[fNJetPtSlices];
3421
3422   TH1F* hRatioPt[fNJetPtSlices]; 
3423   TH1F* hRatioZ[fNJetPtSlices];
3424   TH1F* hRatioXi[fNJetPtSlices];
3425
3426
3427   Int_t axisGenPt = 1;
3428   Int_t axisRecPt = 0;
3429
3430   for(Int_t i=0; i<fNJetPtSlices; i++){ 
3431
3432     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3433     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3434  
3435     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3436     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3437     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3438  
3439     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3440     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3441     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3442  
3443     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3444     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3445     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3446  
3447
3448     hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3449     hFFZRec[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3450     hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3451     
3452     hFFPtRec[i]->SetNameTitle(strNameRecPt,""); 
3453     hFFZRec[i]->SetNameTitle(strNameRecZ,""); 
3454     hFFXiRec[i]->SetNameTitle(strNameRecXi,""); 
3455
3456
3457     hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3458     hFFZGen[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3459     hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3460     
3461     hFFPtGen[i]->SetNameTitle(strNameGenPt,""); 
3462     hFFZGen[i]->SetNameTitle(strNameGenZ,""); 
3463     hFFXiGen[i]->SetNameTitle(strNameGenXi,""); 
3464    
3465
3466     if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
3467     if(fNHistoBinsZ[i])  hFFZRec[i]   = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3468     if(fNHistoBinsXi[i]) hFFXiRec[i]  = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3469     
3470     if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3471     if(fNHistoBinsZ[i])  hFFZGen[i]  = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3472     if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3473     
3474     NormalizeTH1(hFFPtGen[i],fNJets->At(i)); 
3475     NormalizeTH1(hFFZGen[i],fNJets->At(i));
3476     NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3477    
3478     NormalizeTH1(hFFPtRec[i],fNJets->At(i)); 
3479     NormalizeTH1(hFFZRec[i],fNJets->At(i));
3480     NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3481
3482
3483     hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3484     hRatioPt[i]->Reset();
3485     hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3486
3487     hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3488     hRatioZ[i]->Reset();
3489     hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3490     
3491     hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3492     hRatioXi[i]->Reset();
3493     hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3494   }
3495     
3496
3497
3498   // write 
3499   
3500   TFile out(strOutfile,"RECREATE");
3501   
3502   if(!out.IsOpen()){
3503     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3504     return;
3505   }
3506   
3507   for(Int_t i=0; i<fNJetPtSlices; i++){
3508     
3509     hFFPtRec[i]->Write(); 
3510     hFFZRec[i]->Write();
3511     hFFXiRec[i]->Write();
3512     
3513     hFFPtGen[i]->Write(); 
3514     hFFZGen[i]->Write();
3515     hFFXiGen[i]->Write();
3516     
3517     hRatioPt[i]->Write(); 
3518     hRatioZ[i]->Write();
3519     hRatioXi[i]->Write();
3520   }
3521
3522   out.Close();  
3523 }
3524
3525 // ____________________________________________________________________________________________________________________________
3526 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir )
3527
3528   // project response matrix on both axes: 
3529   // pt spec for rec primaries, in terms of generated and reconstructed momentum
3530   // write spec and ratios to outFile
3531   
3532   Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3533   
3534   TH1F* hSpecPtRec;  
3535   TH1F* hSpecPtGen; 
3536   TH1F* hRatioPt; 
3537   
3538   Int_t axisGenPt = 1;
3539   Int_t axisRecPt = 0;
3540   
3541   TString strNameRecPt   = "h1SpecTrackPtRecPrim_recPt";
3542   TString strNameGenPt   = "h1SpecTrackPtRecPrim_genPt";
3543   TString strNameRatioPt = "h1RatioTrackPtRecPrim";
3544     
3545   hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
3546   hSpecPtRec->SetNameTitle(strNameRecPt,""); 
3547   
3548   hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
3549   hSpecPtGen->SetNameTitle(strNameGenPt,""); 
3550   
3551   hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
3552   hRatioPt->Reset();
3553   hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
3554   
3555   TString outfileOption = "RECREATE";
3556   if(updateOutfile) outfileOption = "UPDATE";
3557
3558   TFile out(strOutfile,outfileOption);
3559   
3560   if(!out.IsOpen()){
3561     Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3562     return;
3563   }
3564       
3565
3566   if(strOutDir && strOutDir.Length()){  
3567     
3568     TDirectory* dir;
3569     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
3570     else{
3571       dir = out.mkdir(strOutDir);
3572       dir->cd(); 
3573     } 
3574   }
3575   
3576   hSpecPtRec->Write(); 
3577   hSpecPtGen->Write(); 
3578   hRatioPt->Write(); 
3579   
3580   out.Close();  
3581 }
3582
3583
3584 //__________________________________________________________________________________________________________________________________________________________________
3585 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
3586
3587   // rebin histo, rescale bins according to new width
3588   // only correct for input histos with equal bin size
3589
3590   // args: jetPtSlice, type, use current corr level
3591  
3592   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
3593   // array size of binsLimits: nBinsLimits 
3594   // array size of binsWidth: nBinsLimits-1 
3595   // binsLimits have to be in increasing order
3596   // if binning undefined for any slice, original binning will be kept
3597
3598   if(!fNJetPtSlices){
3599     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
3600     return;
3601   }
3602
3603   if(jetPtSlice>=fNJetPtSlices){
3604     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
3605     return;
3606   }
3607
3608
3609   Double_t binLimitMin = binsLimits[0];
3610   Double_t binLimitMax = binsLimits[nBinsLimits-1];
3611
3612   Double_t binLimit = binLimitMin; // start value 
3613   
3614   Int_t sizeUpperLim = 1000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1; -  only if first bin has smallest width, but not case for dN/dxi ...
3615   TArrayD binsArray(sizeUpperLim);
3616   Int_t nBins = 0; 
3617   binsArray.SetAt(binLimitMin,nBins++);
3618
3619   while(binLimit<binLimitMax && nBins<sizeUpperLim){
3620
3621     Int_t currentSlice = -1;
3622     for(Int_t i=0; i<nBinsLimits; i++){
3623       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
3624     }
3625     
3626     Double_t currentBinWidth = binsWidth[currentSlice];
3627     binLimit += currentBinWidth;
3628
3629     binsArray.SetAt(binLimit,nBins++);
3630   }
3631   
3632
3633   TH1F* hist = 0;
3634   if(type == kFlagPt)            hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice); 
3635   else if(type == kFlagZ)        hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);       
3636   else if(type == kFlagXi)       hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);      
3637   else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0); 
3638   else{ 
3639     Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
3640     return;
3641   }
3642   
3643   Double_t binWidthNoRebin = hist->GetBinWidth(1);
3644
3645   Double_t* bins = binsArray.GetArray();
3646
3647   hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
3648  
3649   for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
3650
3651     Double_t binWidthRebin = hist->GetBinWidth(bin);
3652     Double_t scaleF = binWidthNoRebin / binWidthRebin;
3653
3654     Double_t binCont  = hist->GetBinContent(bin);
3655     Double_t binErr   = hist->GetBinError(bin);
3656     
3657     binCont *= scaleF;
3658     binErr  *= scaleF;
3659
3660     hist->SetBinContent(bin,binCont);
3661     hist->SetBinError(bin,binErr);
3662   }
3663
3664
3665
3666   TH1F* temp = new TH1F(*hist);
3667
3668   if(type == kFlagPt)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
3669   if(type == kFlagZ)        fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
3670   if(type == kFlagXi)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
3671   if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
3672
3673
3674   delete temp;
3675 }
3676 //__________________________________________________________________________________________________________________________________________________________________
3677 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
3678
3679  
3680  if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3681  
3682   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3683
3684   TList* list = 0;
3685
3686   if(strlist && strlist.Length()){ 
3687     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3688       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3689       return;
3690     }
3691   }  
3692   if(list == 0)return; // catch strlist.Lenght() == 0;
3693    
3694   THnSparse* hn6ResponseJetPt  = (THnSparse*) list->FindObject("fhnCorrelation");
3695  
3696   Int_t axis6RecJetPt = 0;
3697   Int_t axis6GenJetPt = 3;
3698
3699   hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)"); 
3700   hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)"); 
3701   
3702   Int_t nBinsRecPt    = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins(); 
3703   Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
3704   Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
3705
3706   Int_t nBinsGenPt    = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins(); 
3707   Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
3708   Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
3709
3710   Int_t nBinsTrackPt = 200;
3711   Int_t loLimTrackPt = 0;
3712   Int_t upLimTrackPt = 200;
3713   
3714
3715   Int_t    nBinsResponse[4]  = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
3716   Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
3717   Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
3718
3719   const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"};
3720   
3721   THnSparseD* hn4ResponseTrackPtJetPt =  new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
3722   
3723   for(Int_t i=0; i<4; i++){
3724     hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]); 
3725   }
3726   
3727
3728   // fill 
3729   
3730
3731 }
3732
3733 //_____________________________________________________________________________________________________________________________________
3734 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
3735 {
3736
3737   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
3738
3739 }
3740
3741 //_____________________________________________________________________________________________________________________________________
3742 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
3743 {
3744
3745   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
3746
3747 }
3748
3749 //_____________________________________________________________________________________________________________________________________
3750 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
3751 {
3752
3753   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
3754
3755 }
3756
3757 //______________________________________________________________________________________________________________________________________________________
3758 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
3759 {
3760   // read single track correction (pt) from file
3761   // type: efficiency / response / secondaries correction 
3762
3763   if(!((type == kFlagEfficiency)  || (type == kFlagResponse) || (type == kFlagSecondaries))){
3764     Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
3765     return;
3766   }
3767
3768   TFile f(strfile,"READ");
3769
3770   if(!f.IsOpen()){
3771     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
3772     return;
3773   }
3774
3775   if(fDebug>0 && type==kFlagEfficiency)  Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3776   if(fDebug>0 && type==kFlagResponse)    Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3777   if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3778  
3779   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3780
3781   TList* list = 0;
3782
3783   if(strlist && strlist.Length()){
3784     
3785     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3786       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3787       return;
3788     }
3789   }  
3790   
3791   TH1F*      h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
3792   THnSparse* hnCorrHist = 0;
3793   
3794   if(type == kFlagEfficiency || type == kFlagSecondaries){
3795
3796     if(list) h1CorrHist = (TH1F*) list->FindObject(strname); 
3797     else     h1CorrHist = (TH1F*) gDirectory->Get(strname);
3798
3799     if(!h1CorrHist){
3800       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3801       return;
3802     }
3803
3804   }
3805   else if(type == kFlagResponse){
3806     
3807     if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
3808     else     hnCorrHist = (THnSparse*) gDirectory->Get(strname);
3809
3810     if(!hnCorrHist){
3811       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3812       return;
3813     }
3814  
3815   }
3816
3817   if(h1CorrHist) h1CorrHist->SetDirectory(0);
3818   //if(hnCorrHist) hnCorrHist->SetDirectory(0);
3819   
3820   f.Close();  
3821
3822   if(type == kFlagEfficiency)       fh1EffSinglePt      = h1CorrHist; 
3823   else if(type == kFlagResponse)    fhnResponseSinglePt = hnCorrHist;
3824   else if(type == kFlagSecondaries) fh1SecCorrSinglePt  = h1CorrHist; 
3825  
3826 }
3827
3828 //________________________________________________________________________________________________________________
3829 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
3830
3831   // read track pt spec from task ouput - standard dir/list 
3832      
3833   TString strdir  = "PWG4_FragmentationFunction_" + strID;
3834   TString strlist = "fracfunc_" + strID;
3835     
3836   ReadRawPtSpec(strInfile,strdir,strlist);
3837 }
3838
3839 //_______________________________________________________________________________________________________
3840 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
3841 {
3842   // get raw pt spectra from file 
3843
3844
3845   // book histos
3846   fNCorrectionLevelsSinglePt = 0; 
3847   fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3848   AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3849
3850   // get raw pt spec from input file, normalize
3851
3852   TFile f(strfile,"READ");
3853
3854   if(!f.IsOpen()){
3855     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3856     return;
3857   }
3858
3859   if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3860
3861   gDirectory->cd(strdir);
3862