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