Adding ability to use multiple test functions
[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  = "PWGJE_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  = "PWGJE_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  = "PWGJE_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     if(!hist){
1459       Printf("%s%d no histo found ",(char*)__FILE__,__LINE__);
1460       return;
1461     }
1462
1463
1464     THnSparse* hnResponse = 0;
1465     if(type == kFlagPt)      hnResponse = fhnResponsePt[i];
1466     else if(type == kFlagZ)  hnResponse = fhnResponseZ[i];
1467     else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1468
1469     TH1F* hPrior = 0;
1470     if(type == kFlagPt && fh1FFTrackPtPrior[i]  && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1471     else if(type == kFlagZ  && fh1FFZPrior[i]   && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0)       ) hPrior = fh1FFZPrior[i];
1472     else if(type == kFlagXi && fh1FFXiPrior[i]  && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0)      ) hPrior = fh1FFXiPrior[i];
1473
1474
1475     TString histNameTHn;
1476     histNameTHn = hist->GetName();
1477     histNameTHn.ReplaceAll("TH1","THn");
1478     
1479
1480     TString priorNameTHn; 
1481     if(hPrior){
1482       priorNameTHn = hPrior->GetName();
1483       priorNameTHn.ReplaceAll("TH1","THn");
1484     }
1485
1486     TString histNameBackFolded;
1487     histNameBackFolded = hist->GetName();
1488     histNameBackFolded.Append("_backfold");
1489     
1490     TString histNameRatioFolded;
1491     histNameRatioFolded = hist->GetName();
1492     histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1493     
1494     histNameRatioFolded.Append("_unfold");
1495
1496     TString histNameRatioBackFolded;
1497     histNameRatioBackFolded = hist->GetName();
1498     histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1499     histNameRatioBackFolded.Append("_backfold");
1500  
1501     THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1502     THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff 
1503     THnSparse* hnPrior          = 0;
1504     if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1505     
1506     THnSparse* hnUnfolded 
1507       = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);  
1508      
1509     TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
1510     hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1511     
1512     if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1513     if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1514     if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1515
1516     // backfolding: apply response matrix to unfolded spectrum
1517     TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); 
1518     hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1519
1520     if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1521     if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
1522     if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
1523     
1524     // ratio unfolded to original histo 
1525     TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1526     hRatioUnfolded->Reset();
1527     hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1528
1529     if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1530     if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
1531     if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;
1532
1533
1534     // ratio backfolded to original histo
1535     TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1536     hRatioBackFolded->Reset();
1537     hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1538
1539     if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1540     if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
1541     if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
1542     
1543     delete hnHist;
1544     delete hnFlatEfficiency;
1545
1546  }
1547 }
1548
1549 //_____________________________________________________________________________________________________
1550 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1551 {
1552   
1553   Int_t type = kFlagPt;
1554   UnfoldHistos(nIter, useCorrelatedErrors, type);
1555 }
1556
1557 //_____________________________________________________________________________________________________
1558 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1559 {
1560   
1561   Int_t type = kFlagZ;
1562   UnfoldHistos(nIter, useCorrelatedErrors, type);
1563 }
1564
1565 //_____________________________________________________________________________________________________
1566 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1567 {
1568   
1569   Int_t type = kFlagXi;
1570   UnfoldHistos(nIter, useCorrelatedErrors, type);
1571 }
1572
1573 //______________________________________________________________________________________________
1574 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1575 {
1576   // apply (multiply) response matrix to hist 
1577
1578   TH1F* hOut = new TH1F(*hist);
1579   hOut->Reset();
1580
1581   const Int_t axisM = 0; 
1582   const Int_t axisT = 1;
1583  
1584   Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1585   Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1586
1587   // response matrix normalization
1588   // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1589   // so for unfolding proper we leave it to CORRFW ...
1590
1591   Double_t normFacResponse[nBinsT];
1592
1593   for(Int_t iT=1; iT<=nBinsT; iT++){
1594
1595     Double_t sumResp = 0;
1596     
1597     for(Int_t iM=1; iM<=nBinsM; iM++){
1598       
1599       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1600       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1601       
1602       Double_t binCoord[] = {coordM,coordT};
1603       
1604       Long64_t binIndex = hnResponse->GetBin(binCoord);
1605       
1606       Double_t resp = hnResponse->GetBinContent(binIndex); 
1607       
1608       sumResp += resp;
1609     }
1610     
1611     normFacResponse[iT] = 0;
1612     if(sumResp) normFacResponse[iT] = 1/sumResp;
1613   }
1614   
1615   
1616   
1617   for(Int_t iM=1; iM<=nBinsM; iM++){
1618     
1619     Double_t contM   = 0;
1620
1621     for(Int_t iT=1; iT<=nBinsT; iT++){
1622
1623       Double_t contT = hist->GetBinContent(iT);
1624       
1625       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1626       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1627
1628       Double_t binCoord[] = {coordM,coordT};
1629       
1630       Long64_t binIndex = hnResponse->GetBin(binCoord);
1631       
1632       Double_t resp = hnResponse->GetBinContent(binIndex); 
1633       
1634       contM   += resp*normFacResponse[iT]*contT; 
1635     }
1636
1637     hOut->SetBinContent(iM,contM);
1638   }
1639
1640   return hOut;
1641 }
1642
1643 //_______________________________________________________________________________________________________
1644 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1645 {
1646   // read reconstruction efficiency from file
1647   // argument strlist optional - read from directory strdir if not specified
1648
1649   // temporary histos to hold histos from file
1650   TH1F* hEffPt[fNJetPtSlices]; 
1651   TH1F* hEffZ[fNJetPtSlices];
1652   TH1F* hEffXi[fNJetPtSlices];
1653   
1654   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1655   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
1656   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1657
1658   TFile f(strfile,"READ");
1659
1660   if(!f.IsOpen()){
1661     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1662     return;
1663   }
1664
1665   if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1666  
1667   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1668
1669   TList* list = 0;
1670
1671   if(strlist && strlist.Length()){
1672    
1673     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1674       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1675       return;
1676     }
1677   }  
1678
1679   for(Int_t i=0; i<fNJetPtSlices; i++){
1680     
1681     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1682     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1683     
1684     TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1685     TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1686     TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1687      
1688    
1689     if(list){
1690       hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt); 
1691       hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ); 
1692       hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi); 
1693     }
1694     else{
1695       hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt); 
1696       hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ); 
1697       hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi); 
1698     }
1699     
1700     if(!hEffPt[i]){
1701       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1702     }
1703   
1704     if(!hEffZ[i]){
1705       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1706     }    
1707
1708     if(!hEffXi[i]){
1709       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1710     }
1711
1712
1713     if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1714     if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1715     if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1716
1717     if(hEffPt[i]) hEffPt[i]->SetDirectory(0); 
1718     if(hEffZ[i])  hEffZ[i]->SetDirectory(0); 
1719     if(hEffXi[i]) hEffXi[i]->SetDirectory(0); 
1720
1721   } // jet slices loop
1722
1723   f.Close();
1724
1725   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1726     if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1727     if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
1728     if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1729   }
1730 }
1731
1732 //___________________________________________________________________________________________________________
1733 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1734 {
1735   // read bgr eff from file
1736   // argument strlist optional - read from directory strdir if not specified
1737  
1738   TH1F* hEffPtBgr[fNJetPtSlices]; 
1739   TH1F* hEffZBgr [fNJetPtSlices];
1740   TH1F* hEffXiBgr[fNJetPtSlices];
1741
1742   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1743   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
1744   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1745
1746
1747   TFile f(strfile,"READ");
1748
1749   if(!f.IsOpen()){
1750     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1751     return;
1752   }
1753
1754   if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1755  
1756   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1757
1758   TList* list = 0;
1759
1760   if(strlist && strlist.Length()){
1761        
1762     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1763       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1764       return;
1765     }
1766   }  
1767
1768   for(Int_t i=0; i<fNJetPtSlices; i++){
1769     
1770     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1771     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1772     
1773     TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1774     TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1775     TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1776   
1777    
1778     if(list){
1779       hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr); 
1780       hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr); 
1781       hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr); 
1782     }
1783     else{
1784       hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr); 
1785       hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr); 
1786       hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr); 
1787     }
1788     
1789     if(!hEffPtBgr[i]){
1790       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1791     }
1792   
1793     if(!hEffZBgr[i]){
1794       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1795     }    
1796
1797     if(!hEffXiBgr[i]){
1798       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1799     }
1800
1801
1802     if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1803     if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1804     if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1805
1806     if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0); 
1807     if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0); 
1808     if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0); 
1809
1810   } // jet slices loop
1811
1812   f.Close();
1813
1814   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1815     if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1816     if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
1817     if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1818   }
1819 }
1820
1821 // ________________________________________________
1822 void AliFragmentationFunctionCorrections::EffCorr()
1823 {
1824   // apply efficiency correction
1825
1826   AddCorrectionLevel("eff");
1827
1828   for(Int_t i=0; i<fNJetPtSlices; i++){
1829
1830     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1831     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1832     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1833
1834     TString histNamePt = histPt->GetName();
1835     TString histNameZ  = histZ->GetName();
1836     TString histNameXi = histXi->GetName();
1837
1838     
1839     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1840     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1841     
1842     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1843     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1844     
1845     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1846     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1847     
1848     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1849   }
1850 }
1851
1852 //___________________________________________________
1853 void AliFragmentationFunctionCorrections::EffCorrBgr()
1854 {
1855   // apply efficiency correction to bgr distributions
1856
1857   AddCorrectionLevelBgr("eff");
1858
1859   Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1860
1861   for(Int_t i=0; i<fNJetPtSlices; i++){
1862
1863     TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1864     TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1865     TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1866     
1867     TString histNamePt = histPt->GetName();
1868     TString histNameZ  = histZ->GetName();
1869     TString histNameXi = histXi->GetName();
1870
1871     
1872     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1873     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1874     
1875     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1876     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1877     
1878     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1879     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1880     
1881     fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1882   }
1883 }
1884
1885 //______________________________________________________________________
1886 void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1887 {
1888   // re-evaluate jet energy after FF corrections from dN/dpt distribution
1889   // 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'
1890   // argument corrlevel: which level of correction to be corrected/shifted to 
1891
1892   if(corrLevel>=fNCorrectionLevels){ 
1893     Printf(" calc xi shift: corrLevel exceeded - do nothing");
1894     return;
1895   }
1896
1897   Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1898   Double_t* jetPtCorr   = new Double_t[fNJetPtSlices];
1899   Double_t* deltaXi     = new Double_t[fNJetPtSlices];
1900
1901   for(Int_t i=0; i<fNJetPtSlices; i++){
1902     
1903     TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1904     TH1F* histPt    = fCorrFF[corrLevel]->GetTrackPt(i);
1905
1906     Double_t ptUncorr = 0;
1907     Double_t ptCorr   = 0;
1908
1909     for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1910
1911       Double_t cont   = histPtRaw->GetBinContent(bin);
1912       Double_t width  = histPtRaw->GetBinWidth(bin);
1913       Double_t meanPt = histPtRaw->GetBinCenter(bin);
1914
1915       ptUncorr += meanPt*cont*width;
1916     }
1917     
1918     for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1919       
1920       Double_t cont   = histPt->GetBinContent(bin);
1921       Double_t width  = histPt->GetBinWidth(bin);
1922       Double_t meanPt = histPt->GetBinCenter(bin);
1923       
1924       ptCorr += meanPt*cont*width;
1925     }
1926
1927     jetPtUncorr[i] = ptUncorr; 
1928     jetPtCorr[i]   = ptCorr;   
1929   }
1930
1931   // calc dXi from dN/dpt distribution : 
1932   // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt 
1933
1934   for(Int_t i=0; i<fNJetPtSlices; i++){
1935
1936     Float_t jetPtLoLim = fJetPtSlices->At(i);
1937     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1938
1939     Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1940     
1941     Double_t ptUncorr = jetPtUncorr[i]; 
1942     Double_t ptCorr   = jetPtCorr[i]; 
1943
1944     Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1945     
1946     Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1947            ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1948     
1949     deltaXi[i] = dXi;
1950   }
1951   
1952   // book & fill new dN/dxi histos
1953
1954   for(Int_t i=0; i<fNJetPtSlices; i++){
1955
1956     TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1957             
1958     Double_t dXi = deltaXi[i];
1959
1960     Int_t nBins  = histXi->GetNbinsX();
1961     const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1962     Float_t binsVecNew[nBins+1];
1963     
1964     TString strName = histXi->GetName(); 
1965     strName.Append("_shift");
1966     TString strTit  = histXi->GetTitle();  
1967
1968     TH1F* hXiCorr; 
1969
1970     // create shifted histo ...
1971
1972     if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1973
1974       for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;    
1975       hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1976     }
1977     else{ // uniform bin size
1978       
1979       Double_t xMin  = histXi->GetXaxis()->GetXmin();
1980       Double_t xMax  = histXi->GetXaxis()->GetXmax();
1981
1982       xMin += dXi;
1983       xMax += dXi;
1984       
1985       hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1986     }
1987
1988     // ... and fill
1989
1990     for(Int_t bin=1; bin<nBins+1; bin++){
1991       Double_t cont = histXi->GetBinContent(bin);
1992       Double_t err  = histXi->GetBinError(bin);
1993       
1994       hXiCorr->SetBinContent(bin,cont);
1995       hXiCorr->SetBinError(bin,err);
1996     }
1997     
1998     new(fh1FFXiShift[i]) TH1F(*hXiCorr);
1999     delete hXiCorr;
2000   }
2001
2002   delete[] jetPtUncorr;
2003   delete[] jetPtCorr;
2004   delete[] deltaXi;
2005 }
2006
2007 //_____________________________________________________
2008 void AliFragmentationFunctionCorrections::SubtractBgr()
2009 {
2010   // subtract bgr distribution from FF
2011   // the current corr level is used for both 
2012   
2013   AddCorrectionLevel("bgrSub");
2014
2015   for(Int_t i=0; i<fNJetPtSlices; i++){
2016
2017     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
2018     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
2019     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
2020     
2021     TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
2022     TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
2023     TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
2024
2025     TString histNamePt = histPt->GetName();
2026     TString histNameZ  = histZ->GetName();
2027     TString histNameXi = histXi->GetName();
2028     
2029
2030     TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
2031     hFFTrackPtBgrSub->Add(histPtBgr,-1);
2032     
2033     TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
2034     hFFZBgrSub->Add(histZBgr,-1);
2035     
2036     TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
2037     hFFXiBgrSub->Add(histXiBgr,-1);
2038
2039     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
2040   }
2041 }
2042
2043 //________________________________________________________________________________________________________________
2044 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile, 
2045                                                               Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2046
2047   // read task ouput from MC and write single track eff - standard dir/list 
2048      
2049   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2050   TString strlist = "fracfunc_" + strID;
2051     
2052   WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2053 }
2054
2055 //___________________________________________________________________________________________________________________________________
2056 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, 
2057                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2058 {
2059   // read task output from MC and write single track eff as function of pt, eta, phi
2060   // argument strlist optional - read from directory strdir if not specified
2061   // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2062
2063
2064   TH1D* hdNdptTracksMCPrimGen;
2065   TH2D* hdNdetadphiTracksMCPrimGen;
2066   
2067   TH1D* hdNdptTracksMCPrimRec;
2068   TH2D* hdNdetadphiTracksMCPrimRec;
2069     
2070
2071   TFile f(strInfile,"READ");
2072
2073   if(!f.IsOpen()){
2074     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2075     return;
2076   }
2077   
2078   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2079  
2080   if(strdir && strdir.Length()){
2081     if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2082     gDirectory->cd(strdir);
2083   }
2084
2085   TList* list = 0;
2086
2087   if(strlist && strlist.Length()){
2088
2089     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2090       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2091       return;
2092     }
2093   }
2094
2095
2096   TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2097   if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data()); 
2098
2099   TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2100   if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2101
2102   TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2103   if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2104   
2105   TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2106   if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2107   
2108
2109   if(list){
2110     hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
2111     hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2112     
2113     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
2114     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2115   }
2116   else{
2117     hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2118     hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2119
2120     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2121     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2122   }
2123
2124   hdNdptTracksMCPrimGen->SetDirectory(0);
2125   hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2126   hdNdptTracksMCPrimRec->SetDirectory(0);
2127   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2128   
2129   f.Close();
2130
2131   // projections: dN/deta, dN/dphi 
2132
2133   TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2134   TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2135  
2136   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2137   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2138
2139   // rebin
2140
2141   TString strNamePtGen = "hTrackPtGenPrim";
2142   TString strNamePtRec = "hTrackPtRecPrim";
2143
2144   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2145   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2146  
2147     // efficiency: divide 
2148
2149   TString hNameTrackEffPt = "hSingleTrackEffPt";
2150   if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2151                                                
2152   TString hNameTrackEffEta = "hSingleTrackEffEta";
2153   if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2154
2155   TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2156   if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2157
2158
2159   TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2160   hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2161
2162   TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2163   hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2164   
2165   TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2166   hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2167   
2168   
2169   TString outfileOption = "RECREATE";
2170   if(updateOutfile)  outfileOption = "UPDATE";
2171
2172   TFile out(strOutfile,outfileOption);
2173
2174   if(!out.IsOpen()){
2175     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2176     return;
2177   }
2178
2179   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2180
2181   if(strOutDir && strOutDir.Length()){
2182     
2183     TDirectory* dir;
2184     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2185     else{
2186       dir = out.mkdir(strOutDir);
2187       dir->cd(); 
2188     } 
2189   }
2190
2191   hSingleTrackEffPt->Write();
2192   hSingleTrackEffEta->Write();
2193   hSingleTrackEffPhi->Write();
2194   
2195   out.Close();
2196 }
2197
2198 //________________________________________________________________________________________________________________
2199 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile, 
2200                                                                   Bool_t updateOutfile, TString strOutDir)
2201
2202   // read task ouput from MC and write single track eff - standard dir/list 
2203      
2204   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2205   TString strlist = "fracfunc_" + strID;
2206     
2207   WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2208 }
2209
2210 //___________________________________________________________________________________________________________________________________
2211 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, 
2212                                                                   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2213 {
2214   // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2215   // argument strlist optional - read from directory strdir if not specified
2216   // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2217
2218   TH1D* hdNdptTracksMCPrimRec;
2219   TH2D* hdNdetadphiTracksMCPrimRec;
2220   
2221   TH1D* hdNdptTracksMCSecRec;
2222   TH2D* hdNdetadphiTracksMCSecRec;
2223
2224   TFile f(strInfile,"READ");
2225
2226   if(!f.IsOpen()){
2227     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2228     return;
2229   }
2230   
2231   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2232  
2233   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2234
2235   TList* list = 0;
2236
2237   if(strlist && strlist.Length()){
2238
2239     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2240       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2241       return;
2242     }
2243   }
2244
2245
2246   if(list){
2247     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2248     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2249     
2250     hdNdptTracksMCSecRec       = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2251     hdNdetadphiTracksMCSecRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2252   }
2253   else{
2254     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2255     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2256
2257     hdNdptTracksMCSecRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2258     hdNdetadphiTracksMCSecRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2259   }
2260   
2261   hdNdptTracksMCPrimRec->SetDirectory(0);
2262   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2263   hdNdptTracksMCSecRec->SetDirectory(0);
2264   hdNdetadphiTracksMCSecRec->SetDirectory(0);
2265   
2266   f.Close();
2267
2268   // projections: dN/deta, dN/dphi 
2269
2270   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2271   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2272  
2273   TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2274   TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2275
2276  
2277   // rebin
2278
2279   TString strNamePtPrim = "hTrackPtPrim";
2280   TString strNamePtSec  = "hTrackPtSec";
2281
2282   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2283   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec  = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2284
2285
2286   // secondary correction factor: divide prim/(prim+sec)
2287
2288   TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2289   TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2290   hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2291   hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2292
2293   TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2294   TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2295   hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2296   hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2297
2298   TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2299   TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2300   hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2301   hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2302   
2303   
2304   TString outfileOption = "RECREATE";
2305   if(updateOutfile)  outfileOption = "UPDATE";
2306
2307   TFile out(strOutfile,outfileOption);
2308
2309   if(!out.IsOpen()){
2310     Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2311     return;
2312   }
2313
2314   if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2315
2316   if(strOutDir && strOutDir.Length()){  
2317
2318     TDirectory* dir;
2319     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2320     else{
2321       dir = out.mkdir(strOutDir);
2322       dir->cd(); 
2323     } 
2324   } 
2325
2326   hdNdptTracksMCSecRec->Write();
2327   hdNdetaTracksMCSecRec->Write();
2328   hdNdphiTracksMCSecRec->Write();
2329
2330   hSingleTrackSecCorrPt->Write();
2331   hSingleTrackSecCorrEta->Write();
2332   hSingleTrackSecCorrPhi->Write();
2333   
2334   out.Close();
2335 }
2336
2337 //________________________________________________________________________________________________________________
2338 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile, 
2339                                                               Bool_t updateOutfile, TString strOutDir)
2340
2341   // read task ouput from MC and write single track eff - standard dir/list 
2342      
2343   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2344   TString strlist = "fracfunc_" + strID;
2345     
2346   WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2347 }
2348
2349
2350 //_____________________________________________________________________________________________________________________________________
2351 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2352                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2353 {
2354   // read 2d THnSparse response matrices in pt from file
2355   // project TH2 
2356   // write to strOutfile 
2357
2358   THnSparse* hnResponseSinglePt; 
2359   TH2F*      h2ResponseSinglePt; 
2360  
2361   TFile f(strInfile,"READ");
2362
2363   if(!f.IsOpen()){
2364     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2365     return;
2366   }
2367
2368   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2369  
2370   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2371
2372   TList* list = 0;
2373
2374   if(strlist && strlist.Length()){
2375     
2376     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2377       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2378       return;
2379     }
2380   }  
2381   
2382   if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2383   else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2384   
2385
2386   if(!hnResponseSinglePt){
2387     Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2388     return;
2389   }
2390
2391   f.Close();
2392
2393
2394   // project 2d histo 
2395   h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2396   h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt",""); 
2397     
2398   
2399   // write 
2400
2401   TString outfileOption = "RECREATE";
2402   if(updateOutfile)  outfileOption = "UPDATE";
2403   
2404   TFile out(strOutfile,outfileOption);
2405   
2406   if(!out.IsOpen()){
2407     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2408     return;
2409   }
2410
2411   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2412   
2413   if(strOutDir && strOutDir.Length()){  
2414
2415     TDirectory* dir;
2416     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
2417     else{
2418       dir = out.mkdir(strOutDir);
2419       dir->cd(); 
2420     } 
2421   }
2422   
2423   hnResponseSinglePt->Write();
2424   h2ResponseSinglePt->Write();
2425   
2426   out.Close();  
2427 }
2428
2429 //________________________________________________________________________________________________________________
2430 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile, 
2431                                                            Bool_t updateOutfile)
2432
2433   // read task ouput from MC and write single track eff - standard dir/list 
2434      
2435   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2436   TString strlist = "fracfunc_" + strID;
2437     
2438   WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2439 }
2440
2441 //___________________________________________________________________________________________________________________________________
2442 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, 
2443                                                            TString strOutfile, Bool_t updateOutfile)
2444 {
2445   // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2446   // argument strlist optional - read from directory strdir if not specified
2447   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2448
2449   TH1F* hEffPt[fNJetPtSlices];
2450   TH1F* hEffXi[fNJetPtSlices];
2451   TH1F* hEffZ[fNJetPtSlices];
2452
2453   TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2454   TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2455   TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2456     
2457   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2458   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2459   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2460
2461
2462   TH1F* fh1FFJetPtRecEffGen;
2463
2464   TH2F* fh2FFTrackPtRecEffGen;
2465   TH2F* fh2FFZRecEffGen;
2466   TH2F* fh2FFXiRecEffGen;
2467   
2468   TH2F* fh2FFTrackPtRecEffRec;
2469   TH2F* fh2FFZRecEffRec;
2470   TH2F* fh2FFXiRecEffRec;
2471  
2472
2473   TFile f(strInfile,"READ");
2474
2475   if(!f.IsOpen()){
2476     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2477     return;
2478   }
2479   
2480   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2481  
2482   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2483
2484   TList* list = 0;
2485
2486   if(strlist && strlist.Length()){
2487
2488     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2489       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2490       return;
2491     }
2492   }
2493
2494   if(list){
2495     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2496
2497     fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2498     fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2499     fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2500     
2501     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2502     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2503     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2504   }
2505   else{
2506     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2507
2508     fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2509     fh2FFZRecEffGen        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2510     fh2FFXiRecEffGen       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2511     
2512     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2513     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2514     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2515   }
2516   
2517   fh1FFJetPtRecEffGen->SetDirectory(0); 
2518
2519   fh2FFTrackPtRecEffGen->SetDirectory(0);
2520   fh2FFZRecEffGen->SetDirectory(0);
2521   fh2FFXiRecEffGen->SetDirectory(0);
2522   
2523   fh2FFTrackPtRecEffRec->SetDirectory(0);
2524   fh2FFZRecEffRec->SetDirectory(0);
2525   fh2FFXiRecEffRec->SetDirectory(0);
2526   
2527   f.Close();
2528
2529
2530   // projections: FF for generated and reconstructed primaries 
2531   
2532   for(Int_t i=0; i<fNJetPtSlices; i++){
2533     
2534     Float_t jetPtLoLim = fJetPtSlices->At(i);
2535     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2536
2537     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2538     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2539
2540     if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2541       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2542       return; 
2543     }
2544     
2545     TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2546     TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2547     TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2548     
2549     TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2550     TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2551     TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2552     
2553     // project 
2554     // appendix 'unbinned' to avoid histos with same name after rebinning
2555
2556     hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2557     hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2558     hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2559     
2560     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2561     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2562     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2563     
2564     // rebin
2565
2566     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2567     if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2568     if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2569
2570     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2571     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2572     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2573
2574     hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2575     hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2576     hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2577     
2578     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2579     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2580     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2581  
2582     // normalize
2583     
2584     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2585
2586     NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin); 
2587     NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin); 
2588     NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin); 
2589
2590     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2591     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2592     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2593     
2594     // divide rec/gen : efficiency
2595
2596     TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2597     TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2598     TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2599  
2600     hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2601     hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2602     
2603     hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2604     hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2605     
2606     hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2607     hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2608   } 
2609   
2610   // write 
2611
2612   TString outfileOption = "RECREATE";
2613   if(updateOutfile)  outfileOption = "UPDATE";
2614
2615   TFile out(strOutfile,outfileOption);
2616   
2617   if(!out.IsOpen()){
2618     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2619     return;
2620   }
2621
2622   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2623
2624 //   if(strdir && strdir.Length()){ 
2625 //     TDirectory* dir = out.mkdir(strdir);
2626 //     dir->cd(); 
2627 //   }
2628
2629   for(Int_t i=0; i<fNJetPtSlices; i++){
2630
2631     hdNdptTracksMCPrimGen[i]->Write();
2632     hdNdxiMCPrimGen[i]->Write();
2633     hdNdzMCPrimGen[i]->Write();
2634     
2635     hdNdptTracksMCPrimRec[i]->Write();
2636     hdNdxiMCPrimRec[i]->Write();
2637     hdNdzMCPrimRec[i]->Write();
2638   
2639     hEffPt[i]->Write();
2640     hEffXi[i]->Write();
2641     hEffZ[i]->Write();
2642   }
2643
2644   out.Close();
2645
2646 }
2647
2648 //________________________________________________________________________________________________________________
2649 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile, 
2650                                                            Bool_t updateOutfile)
2651
2652   // read task ouput from MC and write secondary correction - standard dir/list 
2653      
2654   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2655   TString strlist = "fracfunc_" + strID;
2656     
2657   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2658 }
2659
2660 //___________________________________________________________________________________________________________________________________
2661 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, 
2662                                                            TString strOutfile, Bool_t updateOutfile)
2663 {
2664   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2665   // argument strlist optional - read from directory strdir if not specified
2666   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2667
2668   TH1F* hSecCorrPt[fNJetPtSlices];
2669   TH1F* hSecCorrXi[fNJetPtSlices];
2670   TH1F* hSecCorrZ[fNJetPtSlices];
2671
2672   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2673   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2674   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2675     
2676   TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2677   TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2678   TH1F* hdNdzMCSecRec[fNJetPtSlices];
2679
2680   TH1F* fh1FFJetPtRecEffGen;
2681
2682   TH2F* fh2FFTrackPtRecEffRec;
2683   TH2F* fh2FFZRecEffRec;
2684   TH2F* fh2FFXiRecEffRec;
2685   
2686   TH2F* fh2FFTrackPtSecRec;
2687   TH2F* fh2FFZSecRec;
2688   TH2F* fh2FFXiSecRec;
2689
2690   TFile f(strInfile,"READ");
2691
2692   if(!f.IsOpen()){
2693     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2694     return;
2695   }
2696   
2697   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2698  
2699   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2700
2701   TList* list = 0;
2702
2703   if(strlist && strlist.Length()){
2704
2705     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2706       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2707       return;
2708     }
2709   }
2710
2711   if(list){
2712     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2713
2714     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2715     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2716     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2717     
2718     fh2FFTrackPtSecRec  = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2719     fh2FFZSecRec        = (TH2F*) list->FindObject("fh2FFZSecRec");
2720     fh2FFXiSecRec       = (TH2F*) list->FindObject("fh2FFXiSecRec");
2721   }
2722   else{
2723     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2724
2725     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2726     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2727     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2728     
2729     fh2FFTrackPtSecRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2730     fh2FFZSecRec        = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2731     fh2FFXiSecRec       = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2732   }
2733   
2734   fh1FFJetPtRecEffGen->SetDirectory(0); 
2735
2736   fh2FFTrackPtRecEffRec->SetDirectory(0);
2737   fh2FFZRecEffRec->SetDirectory(0);
2738   fh2FFXiRecEffRec->SetDirectory(0);
2739   
2740   fh2FFTrackPtSecRec->SetDirectory(0);
2741   fh2FFZSecRec->SetDirectory(0);
2742   fh2FFXiSecRec->SetDirectory(0);
2743   
2744   f.Close();
2745
2746
2747   // projections: FF for generated and reconstructed primaries 
2748   
2749   for(Int_t i=0; i<fNJetPtSlices; i++){
2750     
2751     Float_t jetPtLoLim = fJetPtSlices->At(i);
2752     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2753
2754     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2755     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2756
2757     if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2758       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2759       return; 
2760     }
2761     
2762     TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2763     TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2764     TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2765     
2766     TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2767     TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2768     TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2769     
2770     // project 
2771     // appendix 'unbinned' to avoid histos with same name after rebinning
2772
2773     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2774     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2775     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2776     
2777     hdNdptTracksMCSecRec[i]  = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2778     hdNdzMCSecRec[i]         = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2779     hdNdxiMCSecRec[i]        = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2780     
2781     // rebin
2782
2783     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2784     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2785     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2786
2787     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2788     if(fNHistoBinsZ[i])  hdNdzMCSecRec[i]        = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2789     if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i]       = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2790
2791     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2792     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2793     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2794     
2795     hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2796     hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2797     hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2798     
2799     // normalize
2800     
2801     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2802
2803     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2804     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2805     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2806
2807     NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin); 
2808     NormalizeTH1(hdNdzMCSecRec[i],nJetsBin); 
2809     NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin); 
2810     
2811     // divide rec/gen : efficiency
2812
2813     TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2814     TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2815     TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2816  
2817     hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2818     TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2819     hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2820     hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2821
2822     hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2823     TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2824     hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2825     hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2826
2827     hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2828     TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2829     hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2830     hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2831   } 
2832   
2833   // write 
2834
2835   TString outfileOption = "RECREATE";
2836   if(updateOutfile)  outfileOption = "UPDATE";
2837
2838   TFile out(strOutfile,outfileOption);
2839   
2840   if(!out.IsOpen()){
2841     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2842     return;
2843   }
2844   
2845   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2846
2847   for(Int_t i=0; i<fNJetPtSlices; i++){
2848
2849     //  hdNdptTracksMCSecRec[i]->Write();
2850     //  hdNdxiMCSecRec[i]->Write();
2851     //  hdNdzMCSecRec[i]->Write();
2852   
2853     hSecCorrPt[i]->Write();
2854     hSecCorrXi[i]->Write();
2855     hSecCorrZ[i]->Write();
2856   }
2857
2858   out.Close();
2859 }
2860
2861 //________________________________________________________________________________________________________________
2862 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile, 
2863                                                            Bool_t updateOutfile)
2864
2865   // read task ouput from MC and write single track eff - standard dir/list 
2866      
2867   TString strdir  = "PWGJE_FragmentationFunction_" + strID;
2868   TString strlist = "fracfunc_" + strID;
2869     
2870   WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2871 }
2872
2873 //_____________________________________________________________________________________________________________________________________
2874 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2875                                                            TString strOutfile, Bool_t updateOutfile)
2876 {
2877   // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2878   // project THnSparse + TH2 in jet pt slices 
2879   // write to strOutfile 
2880
2881   THnSparse* hn3ResponseJetPt;
2882   THnSparse* hn3ResponseJetZ;
2883   THnSparse* hn3ResponseJetXi;
2884
2885   // 2D response matrices 
2886
2887   THnSparse* hnResponsePt[fNJetPtSlices];
2888   THnSparse* hnResponseZ[fNJetPtSlices];
2889   THnSparse* hnResponseXi[fNJetPtSlices];
2890
2891   TH2F* h2ResponsePt[fNJetPtSlices];
2892   TH2F* h2ResponseZ[fNJetPtSlices];
2893   TH2F* h2ResponseXi[fNJetPtSlices];
2894
2895   // 1D projections on gen pt / rec pt axes
2896
2897   TH1F* h1FFPtRec[fNJetPtSlices]; 
2898   TH1F* h1FFZRec[fNJetPtSlices];
2899   TH1F* h1FFXiRec[fNJetPtSlices];
2900
2901   TH1F* h1FFPtGen[fNJetPtSlices]; 
2902   TH1F* h1FFZGen[fNJetPtSlices];
2903   TH1F* h1FFXiGen[fNJetPtSlices];
2904
2905   TH1F* h1RatioPt[fNJetPtSlices]; 
2906   TH1F* h1RatioZ[fNJetPtSlices];
2907   TH1F* h1RatioXi[fNJetPtSlices];
2908
2909
2910
2911   TFile f(strInfile,"READ");
2912
2913   if(!f.IsOpen()){
2914     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2915     return;
2916   }
2917
2918   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2919  
2920   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2921
2922   TList* list = 0;
2923
2924   if(strlist && strlist.Length()){
2925     
2926     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2927       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2928       return;
2929     }
2930   }  
2931   
2932   if(list){
2933     hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2934     hn3ResponseJetZ  = (THnSparse*) list->FindObject("fhnResponseJetZ");
2935     hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2936   }
2937   else{
2938     hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2939     hn3ResponseJetZ  = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2940     hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2941   }
2942
2943   
2944   if(!hn3ResponseJetPt){
2945     Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2946     return;
2947   }
2948
2949   if(!hn3ResponseJetZ){
2950     Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2951     return;
2952   }
2953
2954   if(!hn3ResponseJetXi){
2955     Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2956     return;
2957   }
2958
2959   f.Close();  
2960
2961   // axes 
2962
2963   Int_t axisJetPtTHn3 = -1;
2964   Int_t axisGenPtTHn3 = -1;
2965   Int_t axisRecPtTHn3 = -1;
2966
2967   for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2968     
2969     TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle(); 
2970
2971     if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i; 
2972     if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i; 
2973     if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i; 
2974   }
2975
2976   if(axisJetPtTHn3 == -1){
2977     Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2978     return;
2979   }
2980
2981   if(axisGenPtTHn3 == -1){
2982     Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2983     return;
2984   }
2985
2986   if(axisRecPtTHn3 == -1){
2987     Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2988     return;
2989   }
2990
2991   for(Int_t i=0; i<fNJetPtSlices; i++){
2992     
2993     Float_t jetPtLoLim = fJetPtSlices->At(i);
2994     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2995
2996     Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2997     Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2998
2999     if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
3000       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
3001       return; 
3002     }
3003     
3004     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3005     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3006     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3007
3008     TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3009     TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3010     TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
3011          
3012     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3013     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3014     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3015  
3016     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3017     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3018     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3019  
3020     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3021     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3022     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
3023  
3024     
3025     // 2D projections in jet pt range
3026   
3027     hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
3028     hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
3029     hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
3030
3031     Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};  
3032
3033     hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
3034     hnResponseZ[i]  = hn3ResponseJetZ->Projection(2,axesProj);
3035     hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
3036   
3037     hnResponsePt[i]->SetNameTitle(strNameRespPt,""); 
3038     hnResponseZ[i]->SetNameTitle(strNameRespZ,""); 
3039     hnResponseXi[i]->SetNameTitle(strNameRespXi,""); 
3040
3041     h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
3042     h2ResponseZ[i]  = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
3043     h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
3044  
3045     h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,""); 
3046     h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,""); 
3047     h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,""); 
3048
3049
3050     // 1D projections
3051
3052     Int_t axisGenPtTHn2 = -1;
3053     Int_t axisRecPtTHn2 = -1;
3054
3055     for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3056     
3057       TString title = hnResponsePt[i]->GetAxis(d)->GetTitle(); 
3058
3059       if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d; 
3060       if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d; 
3061     }
3062
3063     
3064     if(axisGenPtTHn2 == -1){
3065       Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3066       return;
3067     }
3068     
3069     if(axisRecPtTHn2 == -1){
3070       Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3071       return;
3072     }
3073     
3074
3075     h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3076     h1FFZRec[i]  = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3077     h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3078     
3079     h1FFPtRec[i]->SetNameTitle(strNameRecPt,""); 
3080     h1FFZRec[i]->SetNameTitle(strNameRecZ,""); 
3081     h1FFXiRec[i]->SetNameTitle(strNameRecXi,""); 
3082
3083     h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3084     h1FFZGen[i]  = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3085     h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3086     
3087     h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3088     h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3089     h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3090
3091     // normalize 1D projections
3092
3093     if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
3094     if(fNHistoBinsZ[i])  h1FFZRec[i]  = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3095     if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3096     
3097     if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3098     if(fNHistoBinsZ[i])  h1FFZGen[i]  = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3099     if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3100        
3101     NormalizeTH1(h1FFPtRec[i],fNJets->At(i)); 
3102     NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3103     NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3104
3105     NormalizeTH1(h1FFPtGen[i],fNJets->At(i)); 
3106     NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3107     NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3108
3109     // ratios 1D projections
3110
3111     h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3112     h1RatioPt[i]->Reset();
3113     h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3114
3115     h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3116     h1RatioZ[i]->Reset();
3117     h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3118     
3119     h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3120     h1RatioXi[i]->Reset();
3121     h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3122   }
3123   
3124   
3125   // write 
3126
3127   TString outfileOption = "RECREATE";
3128   if(updateOutfile)  outfileOption = "UPDATE";
3129   
3130   TFile out(strOutfile,outfileOption);
3131   
3132   if(!out.IsOpen()){
3133     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3134     return;
3135   }
3136
3137   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3138
3139   //if(strdir && strdir.Length()){ 
3140   //  TDirectory* dir = out.mkdir(strdir);
3141   //  dir->cd(); 
3142   //}
3143
3144   for(Int_t i=0; i<fNJetPtSlices; i++){
3145
3146     hnResponsePt[i]->Write();
3147     hnResponseXi[i]->Write();
3148     hnResponseZ[i]->Write();
3149
3150     h2ResponsePt[i]->Write();
3151     h2ResponseXi[i]->Write();
3152     h2ResponseZ[i]->Write();
3153
3154     h1FFPtRec[i]->Write(); 
3155     h1FFZRec[i]->Write();
3156     h1FFXiRec[i]->Write();
3157     
3158     h1FFPtGen[i]->Write(); 
3159     h1FFZGen[i]->Write();
3160     h1FFXiGen[i]->Write();
3161     
3162     h1RatioPt[i]->Write(); 
3163     h1RatioZ[i]->Write();
3164     h1RatioXi[i]->Write();
3165
3166   }
3167
3168   out.Close();  
3169 }
3170
3171 //______________________________________________________________________________________________________
3172 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3173 {
3174   // read response matrices from file
3175   // argument strlist optional - read from directory strdir if not specified
3176   // note: THnSparse are not rebinned 
3177
3178   THnSparse* hRespPt[fNJetPtSlices]; 
3179   THnSparse* hRespZ[fNJetPtSlices];
3180   THnSparse* hRespXi[fNJetPtSlices];
3181   
3182   for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3183   for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i]  = 0;
3184   for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3185
3186   TFile f(strfile,"READ");
3187
3188   if(!f.IsOpen()){
3189     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3190     return;
3191   }
3192
3193   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3194  
3195   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3196
3197   TList* list = 0;
3198
3199   if(strlist && strlist.Length()){
3200    
3201     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3202       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3203       return;
3204     }
3205   }  
3206
3207   for(Int_t i=0; i<fNJetPtSlices; i++){
3208     
3209     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3210     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3211     
3212     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3213     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3214     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3215         
3216     if(list){
3217       hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt); 
3218       hRespZ[i]  = (THnSparse*) list->FindObject(strNameRespZ); 
3219       hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi); 
3220     }
3221     else{
3222       hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt); 
3223       hRespZ[i]  = (THnSparse*) gDirectory->Get(strNameRespZ); 
3224       hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi); 
3225     }
3226     
3227     if(!hRespPt[i]){
3228       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3229     }
3230   
3231     if(!hRespZ[i]){
3232       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3233     }    
3234
3235     if(!hRespXi[i]){
3236       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3237     }
3238    
3239     //    if(0){ // can't rebin THnSparse ...
3240     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3241     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3242     
3243     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3244     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3245     
3246     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3247     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3248     //     }
3249    
3250  
3251   } // jet slices loop
3252
3253   f.Close();  // THnSparse pointers still valid even if file closed
3254
3255 //   for(Int_t i=0; i<fNJetPtSlices; i++){  // no copy c'tor ...
3256 //     if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3257 //     if(hRespZ[i])  new(fhnResponseZ[i])  THnSparseF(*hRespZ[i]);
3258 //     if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3259 //   }
3260
3261   for(Int_t i=0; i<fNJetPtSlices; i++){ 
3262     fhnResponsePt[i] = hRespPt[i];
3263     fhnResponseZ[i]  = hRespZ[i];
3264     fhnResponseXi[i] = hRespXi[i];
3265   }
3266 }
3267
3268 //______________________________________________________________________________________________________________________
3269 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3270 {
3271   // read priors from file: rec primaries, gen pt dist
3272
3273   if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3274   
3275   // temporary histos to store pointers from file
3276   TH1F* hist[fNJetPtSlices]; 
3277
3278   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3279
3280   TFile f(strfile,"READ");
3281
3282   if(!f.IsOpen()){
3283     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3284     return;
3285   }
3286
3287   for(Int_t i=0; i<fNJetPtSlices; i++){
3288     
3289     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3290     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3291     
3292     TString strName;
3293     
3294     if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); 
3295     if(type == kFlagZ)  strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);       
3296     if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);      
3297     
3298     hist[i] = (TH1F*) gDirectory->Get(strName); 
3299     
3300     if(!hist[i]){
3301       Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3302     }
3303     
3304
3305     //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3306
3307     if(hist[i]) hist[i]->SetDirectory(0); 
3308     
3309   } // jet slices loop
3310