coverity and coding conventions (O. Busch)
[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 AliFragFuncCorrHistos& copy)
307   : TObject()
308   ,fArraySize(copy.fArraySize)
309   ,fh1CorrFFTrackPt(copy.fh1CorrFFTrackPt)
310   ,fh1CorrFFZ(copy.fh1CorrFFZ)
311   ,fh1CorrFFXi(copy.fh1CorrFFXi)
312   ,fCorrLabel(copy.fCorrLabel)
313 {
314   // copy constructor
315   for(Int_t i=0; i<copy.fArraySize; i++){
316     fh1CorrFFTrackPt[i] = copy.fh1CorrFFTrackPt[i];
317     fh1CorrFFZ[i]       = copy.fh1CorrFFZ[i];
318     fh1CorrFFXi[i]      = copy.fh1CorrFFXi[i];
319   }
320 }
321
322 //_______________________________________________________________________________________________________________________________________________________________
323 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
324 {
325   // assignment
326   
327   if(this!=&o){
328     TObject::operator=(o);
329     fArraySize       = o.fArraySize;
330     fCorrLabel       = o.fCorrLabel;
331  
332     for(Int_t i=0; i<o.fArraySize; i++){
333       fh1CorrFFTrackPt[i] = o.fh1CorrFFTrackPt[i];
334       fh1CorrFFZ[i]       = o.fh1CorrFFZ[i];
335       fh1CorrFFXi[i]      = o.fh1CorrFFXi[i];
336     }
337   }
338     
339   return *this;
340 }
341
342 //__________________________________________________________________________________
343 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
344 {
345   // destructor 
346
347   if(fArraySize){
348     for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
349     for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
350     for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];          
351   }
352
353   if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
354   if(fh1CorrFFZ)       delete[] fh1CorrFFZ;
355   if(fh1CorrFFXi)      delete[] fh1CorrFFXi;
356
357 }
358
359 //___________________________________________________________________________________________________________________
360 AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
361   : TObject()
362   ,fArraySize(0)
363   ,fh1CorrFFTrackPt(0)
364   ,fh1CorrFFZ(0)
365   ,fh1CorrFFXi(0)
366   ,fCorrLabel(label)
367 {
368   // constructor
369
370   fArraySize = arraySize;
371   fh1CorrFFTrackPt = new TH1F*[arraySize];
372   fh1CorrFFZ       = new TH1F*[arraySize];
373   fh1CorrFFXi      = new TH1F*[arraySize];
374
375   for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
376   for(int i=0; i<arraySize; i++) fh1CorrFFZ[i]       = new TH1F;
377   for(int i=0; i<arraySize; i++) fh1CorrFFXi[i]      = new TH1F;
378 }
379
380 //_______________________________________________________________________________________________________________________________
381 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
382 {
383   // place histo in array, add corrLabel to histo name
384
385   if(slice>=fArraySize){
386     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
387     return;
388   }
389
390   if(histPt){
391     TString name = histPt->GetName();
392     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
393     histPt->SetName(name);
394     
395     if(!(histPt->GetSumw2())) histPt->Sumw2(); 
396
397     new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
398   }
399   
400   if(histZ){
401     TString name = histZ->GetName();
402     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
403     histZ->SetName(name);
404
405     if(!(histZ->GetSumw2())) histZ->Sumw2(); 
406
407     new(fh1CorrFFZ[slice]) TH1F(*histZ);
408   }
409
410   if(histXi){
411     TString name = histXi->GetName();
412     if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
413     histXi->SetName(name);
414
415     if(!(histXi->GetSumw2())) histXi->Sumw2(); 
416
417     new(fh1CorrFFXi[slice]) TH1F(*histXi);
418   }
419 }
420
421 //___________________________________________________________________________________________________________________________________
422 void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
423 {
424   // replace histo in array 
425
426   if(slice>=fArraySize){
427     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
428     return;
429   }
430
431   if(histPt){
432     if(!(histPt->GetSumw2())) histPt->Sumw2(); 
433
434     TString name = histPt->GetName();
435     histPt->SetName(name);
436
437     delete fh1CorrFFTrackPt[slice];
438     fh1CorrFFTrackPt[slice] =  new TH1F;
439     new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
440   }
441   
442   if(histZ){
443     if(!(histZ->GetSumw2())) histZ->Sumw2(); 
444
445     TString name = histZ->GetName();
446     histZ->SetName(name);
447
448     delete fh1CorrFFZ[slice];
449     fh1CorrFFZ[slice] = new TH1F;
450     new(fh1CorrFFZ[slice]) TH1F(*histZ);
451   }
452
453   if(histXi){
454     if(!(histXi->GetSumw2())) histXi->Sumw2(); 
455     
456     TString name = histXi->GetName();
457     histXi->SetName(name);
458
459     delete fh1CorrFFXi[slice];
460     fh1CorrFFXi[slice] = new TH1F;
461     new(fh1CorrFFXi[slice]) TH1F(*histXi);
462   }
463 }
464
465 // ___________________________________________________________________________________________
466 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
467
468   // return pt histo 
469   
470   if(slice>=fArraySize){
471     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
472     return 0;
473   }
474
475   return fh1CorrFFTrackPt[slice]; 
476
477 }
478
479 // ______________________________________________________________________________________
480 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
481
482   // return z histo 
483   
484   if(slice>=fArraySize){
485     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
486     return 0;
487   }
488
489   return fh1CorrFFZ[slice]; 
490 }
491
492 // ________________________________________________________________________________________
493 TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
494
495   // return xi histo
496
497   if(slice>=fArraySize){
498     Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
499     return 0;
500   }
501
502   return fh1CorrFFXi[slice]; 
503 }
504
505 // __________________________________________________________________________
506 void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
507 {
508   // delete array of TH1 
509  
510   if(hist){
511     for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
512     delete[] hist;
513   }
514 }
515
516 // ____________________________________________________________________________________
517 void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
518 {
519   // delete array of THnSparse 
520  
521   if(hist){
522     for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
523     delete[] hist;
524   }
525 }
526
527 // _________________________________________________________
528 TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
529 {
530   // book array of TH1
531
532   if(!fNJetPtSlices){
533     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
534     return 0;
535   }
536
537   TH1F** hist = new TH1F*[fNJetPtSlices];
538   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
539   
540   return hist;
541 }
542
543 //__________________________________________________________________
544 THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
545 {
546   // book array of THnSparse
547
548   if(!fNJetPtSlices){
549     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
550     return 0;
551   }
552
553   THnSparse** hist = new THnSparse*[fNJetPtSlices];
554   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF(); 
555   
556   return hist;
557 }
558
559 //_____________________________________________________________________________
560 void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
561 {
562   // increase corr level 
563
564   if(!fNJetPtSlices){
565     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
566     return;
567   }
568
569   if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
570     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
571     return;
572   }
573
574   fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
575   fNCorrectionLevels++;
576 }
577
578
579 //________________________________________________________________________________
580 void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
581 {
582   // increase corr level for bgr FF
583
584   if(!fNJetPtSlices){
585     if(fDebug>0)  Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
586     return;
587   }
588   
589   if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
590     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
591     return;
592   }
593
594   fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
595   fNCorrectionLevelsBgr++;
596 }
597
598 //_____________________________________________________________________________________
599 void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
600 {
601   // increase corr level single pt spec
602
603   Int_t nJetPtSlicesSingle = 1; // pro forma 
604
605   if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
606     Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
607     return;
608   }
609
610   fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
611   fNCorrectionLevelsSinglePt++;
612 }
613
614 //_____________________________________________________________________________________________
615 void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
616
617   // set jet pt slices
618   // once slices are known, can book all other histos 
619
620   fNJetPtSlices = nJetPtSlices;
621
622   const Int_t size = nJetPtSlices+1;
623   fJetPtSlices  = new TArrayF(size,bins);
624  
625   // nJets array
626
627   fNJets = new TArrayF(size);
628   fNJets->Reset(0);
629
630   fNJetsBgr = new TArrayF(size);
631   fNJetsBgr->Reset(0);
632
633   // histos 
634
635   fNCorrectionLevels = 0; 
636   fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
637   AddCorrectionLevel(); // first 'correction' level = raw FF
638
639   fNCorrectionLevelsBgr = 0; 
640   fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
641   AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist
642
643   fh1FFXiShift = BookHistoArray();
644
645   // eff histos
646
647   fh1EffPt = BookHistoArray();
648   fh1EffXi = BookHistoArray();
649   fh1EffZ  = BookHistoArray();
650
651   fh1EffBgrPt = BookHistoArray();
652   fh1EffBgrXi = BookHistoArray();
653   fh1EffBgrZ  = BookHistoArray();
654
655
656   // unfolding
657
658   fh1FFTrackPtBackFolded = BookHistoArray();
659   fh1FFXiBackFolded      = BookHistoArray();
660   fh1FFZBackFolded       = BookHistoArray();
661
662   fh1FFRatioTrackPtFolded = BookHistoArray();
663   fh1FFRatioXiFolded      = BookHistoArray();
664   fh1FFRatioZFolded       = BookHistoArray();
665   
666   fh1FFRatioTrackPtBackFolded = BookHistoArray();
667   fh1FFRatioXiBackFolded      = BookHistoArray();
668   fh1FFRatioZBackFolded       = BookHistoArray();
669
670   //
671
672   fhnResponsePt      = BookTHnSparseArray(); 
673   fhnResponseZ       = BookTHnSparseArray();
674   fhnResponseXi      = BookTHnSparseArray();
675
676   fh1FFTrackPtPrior =  BookHistoArray();
677   fh1FFZPrior       =  BookHistoArray();
678   fh1FFXiPrior      =  BookHistoArray();
679
680
681   // histos bins 
682
683   fNHistoBinsPt = new Int_t[fNJetPtSlices];
684   fNHistoBinsXi = new Int_t[fNJetPtSlices];
685   fNHistoBinsZ  = new Int_t[fNJetPtSlices];
686
687   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
688   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
689   for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i]  = 0;
690   
691   fHistoBinsPt = new TArrayD*[fNJetPtSlices];
692   fHistoBinsXi = new TArrayD*[fNJetPtSlices];
693   fHistoBinsZ  = new TArrayD*[fNJetPtSlices];
694
695   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
696   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
697   for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i]  = new TArrayD(0);
698 }
699
700 //_____________________________________________________________________________________________________________________________________
701 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
702
703   // set histo bins for jet pt slice
704   // if binning undefined for any slice, original binning will be used
705
706   if(!fNJetPtSlices){
707     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
708     return;
709   }
710   
711   if(jetPtSlice>=fNJetPtSlices){
712     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
713     return;
714   }
715   
716   if(type == kFlagPt){
717     fNHistoBinsPt[jetPtSlice] = sizeBins-1;
718     fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins); 
719   }
720   else if(type==kFlagZ){
721     fNHistoBinsZ[jetPtSlice] = sizeBins-1;
722     fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins); 
723   }
724
725   else if(type==kFlagXi){
726     fNHistoBinsXi[jetPtSlice] = sizeBins-1;
727     fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins); 
728   }
729 }
730
731 //__________________________________________________________________________________________________________________________________________________________________
732 void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
733
734   // set histo bins for jet pt slice
735   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
736   // array size of binsLimits: nBinsLimits 
737   // array size of binsWidth: nBinsLimits-1 
738   // binsLimits have to be in increasing order
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
752   Double_t binLimitMin = binsLimits[0];
753   Double_t binLimitMax = binsLimits[nBinsLimits-1];
754
755   Double_t binLimit = binLimitMin; // start value 
756   
757   Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
758   TArrayD binsArray(sizeUpperLim);
759   Int_t nBins = 0; 
760   binsArray.SetAt(binLimitMin,nBins++);
761
762   while(binLimit<binLimitMax && nBins<sizeUpperLim){
763
764     Int_t currentSlice = -1;
765     for(Int_t i=0; i<nBinsLimits; i++){
766       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
767     }
768     
769     Double_t currentBinWidth = binsWidth[currentSlice];
770     binLimit += currentBinWidth;
771
772     binsArray.SetAt(binLimit,nBins++);
773   }
774   
775   Double_t* bins = binsArray.GetArray();
776
777   SetHistoBins(jetPtSlice,nBins,bins,type); 
778 }
779
780 //__________________________________________________________________________________________________
781 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
782
783   // set histo bins for inclusive pt spectra
784   // if binning undefined, original binning will be used
785
786   fNHistoBinsSinglePt = sizeBins-1;
787   
788   fHistoBinsSinglePt = new TArrayD(sizeBins);
789   fHistoBinsSinglePt->Set(sizeBins,bins);  
790 }
791
792 //__________________________________________________________________________________________________________________________________________________________________
793 void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
794
795   // set histo bins for inclusive pt spectra 
796   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
797   // array size of binsLimits: nBinsLimits 
798   // array size of binsWidth: nBinsLimits-1 
799   // binsLimits have to be in increasing order
800   // if binning undefined for any slice, original binning will be used
801
802
803   Double_t binLimitMin = binsLimits[0];
804   Double_t binLimitMax = binsLimits[nBinsLimits-1];
805
806   Double_t binLimit = binLimitMin; // start value 
807   
808   Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
809   TArrayD binsArray(sizeUpperLim);
810   Int_t nBins = 0; 
811   binsArray.SetAt(binLimitMin,nBins++);
812
813   while(binLimit<binLimitMax && nBins<sizeUpperLim){
814
815     Int_t currentSlice = -1;
816     for(Int_t i=0; i<nBinsLimits; i++){
817       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
818     }
819     
820     Double_t currentBinWidth = binsWidth[currentSlice];
821     binLimit += currentBinWidth;
822
823     binsArray.SetAt(binLimit,nBins++);
824   }
825   
826   Double_t* bins = binsArray.GetArray();
827   
828   SetHistoBinsSinglePt(nBins,bins); 
829 }
830
831 //____________________________________________________________________________________
832 void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
833 {
834   // FF normalization: divide by bin width and normalize to entries/jets
835   // should also work for TH2, but to be tested !!!
836
837   if(nJets == 0){ // nothing to do
838     if(fDebug>0)  Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
839     return; 
840   }
841   
842   Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();
843
844   for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)
845
846     Double_t binWidth = hist->GetBinWidth(bin);
847     Double_t binCont  = hist->GetBinContent(bin);
848     Double_t binErr   = hist->GetBinError(bin);
849     
850     binCont /= binWidth;
851     binErr  /= binWidth;
852
853     hist->SetBinContent(bin,binCont);
854     hist->SetBinError(bin,binErr);
855   }
856
857   hist->Scale(1/nJets);
858 }
859
860 //_____________________________________________________
861 void AliFragmentationFunctionCorrections::NormalizeFF()
862 {
863   // normalize FF
864
865   if(fNCorrectionLevels>1){
866     Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
867   }
868
869   for(Int_t i=0; i<fNJetPtSlices; i++){
870     
871     if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));
872
873     NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
874     NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
875     NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
876   } 
877 }
878
879 //______________________________________________________
880 void AliFragmentationFunctionCorrections::NormalizeBgr()
881 {
882   // normalize bgr/UE FF
883
884   if(fNCorrectionLevelsBgr>1){
885     Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
886   }
887
888   for(Int_t i=0; i<fNJetPtSlices; i++){
889     NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
890     NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
891     NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
892   } 
893
894 }
895
896 //__________________________________________________________________________________________________
897 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
898
899   // read raw FF - standard dir/list name
900      
901   TString strdir  = "PWG4_FragmentationFunction_" + strID;
902   TString strlist = "fracfunc_" + strID;
903     
904   ReadRawFF(strfile,strdir,strlist,strFFID);
905 }
906
907 //____________________________________________________________________________________________________________________
908 void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
909 {
910   // get raw FF from input file, project in jet pt slice
911   // normalization done separately 
912
913   TFile f(strfile,"READ");
914
915   if(!f.IsOpen()){
916     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
917     return;
918   }
919
920   if(fDebug>0) Printf("%s:%d -- read FF from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
921
922   gDirectory->cd(strdir);
923
924   TList* list = 0;
925   
926   if(!(list = (TList*) gDirectory->Get(strlist))){ 
927     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
928     return;
929   }
930
931   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
932   TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
933   TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
934   TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));
935
936   TH1F* fh1FFJetPt   = (TH1F*) list->FindObject(hnameJetPt);
937   TH2F* fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
938   TH2F* fh2FFZ       = (TH2F*) list->FindObject(hnameZ);  
939   TH2F* fh2FFXi      = (TH2F*) list->FindObject(hnameXi); 
940
941   if(!fh1FFJetPt)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());   return; }
942   if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
943   if(!fh2FFZ)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data());       return; }
944   if(!fh2FFXi)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data());      return; }
945
946   fh1FFJetPt->SetDirectory(0);
947   fh2FFTrackPt->SetDirectory(0);
948   fh2FFZ->SetDirectory(0);  
949   fh2FFXi->SetDirectory(0); 
950
951   f.Close();  
952
953
954   // nJets per bin
955
956   for(Int_t i=0; i<fNJetPtSlices; i++){
957     
958     Float_t jetPtLoLim = fJetPtSlices->At(i);
959     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
960     
961     Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
962     Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
963                                      
964     Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
965     
966     fNJets->SetAt(nJetsBin,i); 
967
968     if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
969   }
970   
971   // projections: FF 
972   
973   for(Int_t i=0; i<fNJetPtSlices; i++){
974     
975     Float_t jetPtLoLim = fJetPtSlices->At(i);
976     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
977
978     Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
979     Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;
980
981     if(binUp > fh2FFTrackPt->GetNbinsX()){
982       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
983       return; 
984     }
985     
986     TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
987     TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
988     TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
989     
990     // appendix 'unbinned' to avoid histos with same name after rebinning
991     TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
992     TH1F* projZ  = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
993     TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
994     
995     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
996     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
997     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());
998
999     projPt->SetNameTitle(strNameFFPt,"");
1000     projZ->SetNameTitle(strNameFFZ,"");
1001     projXi->SetNameTitle(strNameFFXi,"");
1002
1003     // raw FF = corr level 0
1004     fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
1005   }  
1006 }
1007
1008 //_____________________________________________________________________________________________________________________
1009 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
1010
1011   // read raw FF - standard dir/list name
1012      
1013   TString strdir  = "PWG4_FragmentationFunction_" + strID;
1014   TString strlist = "fracfunc_" + strID;
1015     
1016   ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
1017 }
1018
1019 //_______________________________________________________________________________________________________________________________________
1020 void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
1021 {
1022   // get raw FF from input file, project in jet pt slice
1023   // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
1024   // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
1025   // normalization done separately 
1026
1027   TString strID = strBgrID + strFFID;
1028   
1029   TFile f(strfile,"READ");
1030
1031   if(!f.IsOpen()){
1032     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1033     return;
1034   }
1035
1036   if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data());
1037
1038   gDirectory->cd(strdir);
1039
1040   TList* list = 0;
1041   
1042   if(!(list = (TList*) gDirectory->Get(strlist))){
1043     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1044     return;
1045   }
1046
1047   TString hnameNJets = "fh1nRecJetsCuts"; 
1048   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
1049   TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
1050   TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
1051   TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));
1052
1053   TH1F* fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1054   TH1F* fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
1055   TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1056   TH2F* fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);  
1057   TH2F* fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi); 
1058
1059   if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
1060   if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());      return; }
1061   if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1062   if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
1063   if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }
1064
1065   fh1FFJetPtBgr->SetDirectory(0);
1066   fh1NJets->SetDirectory(0);
1067   fh2FFTrackPtBgr->SetDirectory(0);
1068   fh2FFZBgr->SetDirectory(0);  
1069   fh2FFXiBgr->SetDirectory(0); 
1070
1071   f.Close();  
1072
1073   // nJets per bin
1074
1075   for(Int_t i=0; i<fNJetPtSlices; i++){
1076     
1077     Float_t jetPtLoLim = fJetPtSlices->At(i);
1078     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1079     
1080     Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1081     Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1082                                      
1083     Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1084     Double_t scaleF = 1;
1085
1086     //if(strBgrID.Contains("Out2Jets")){  // scale by ratio 2 jets events / all events
1087     //  scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX()) 
1088     //  / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}
1089
1090
1091     if(strBgrID.Contains("OutAllJets")){  // scale by ratio >3 jets events / all events
1092       scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX()) 
1093         / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
1094     }
1095   
1096     fNJetsBgr->SetAt(nJetsBin*scaleF,i); 
1097
1098     if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1099                         static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1100
1101   }
1102   
1103   // projections: FF 
1104   
1105   for(Int_t i=0; i<fNJetPtSlices; i++){
1106     
1107     Float_t jetPtLoLim = fJetPtSlices->At(i);
1108     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1109
1110     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1111     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1112
1113     if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1114       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
1115       return; 
1116     }
1117     
1118     TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1119     TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1120     TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1121     
1122     // appendix 'unbinned' to avoid histos with same name after rebinning
1123     TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
1124     TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1125     TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1126     
1127     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1128     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1129     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1130
1131     projPt->SetNameTitle(strNameBgrPt,"");
1132     projZ->SetNameTitle(strNameBgrZ,"");
1133     projXi->SetNameTitle(strNameBgrXi,"");
1134     
1135     // raw bgr = corr level 0
1136     fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1137   }  
1138 }
1139
1140 //_____________________________________________________________________________________________________________________
1141 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
1142
1143   // read raw FF - standard dir/list name
1144      
1145   TString strdir  = "PWG4_FragmentationFunction_" + strID;
1146   TString strlist = "fracfunc_" + strID;
1147     
1148   ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
1149 }
1150
1151 //_______________________________________________________________________________________________________________________________________
1152 void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
1153 {
1154   // get raw FF from input file, project in jet pt slice
1155   // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
1156   // normalization done separately 
1157
1158   TString strBgrID = "BckgEmbed";
1159   TString strID = strBgrID + strFFID;
1160
1161   TFile f(strfile,"READ");
1162
1163   if(!f.IsOpen()){
1164     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1165     return;
1166   }
1167
1168   if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());
1169
1170   gDirectory->cd(strdir);
1171
1172   TList* list = 0;
1173   
1174   if(!(list = (TList*) gDirectory->Get(strlist))){
1175     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1176     return;
1177   }
1178
1179   TString hnameNJets = "fh1nRecJetsCuts"; 
1180   TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); 
1181   TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
1182   TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
1183   TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));
1184
1185   TH1F* fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
1186   TH1F* fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
1187   TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
1188   TH2F* fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);  
1189   TH2F* fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi); 
1190
1191   if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
1192   if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());      return; }
1193   if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
1194   if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
1195   if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }
1196
1197   fh1FFJetPtBgr->SetDirectory(0);
1198   fh1NJets->SetDirectory(0);
1199   fh2FFTrackPtBgr->SetDirectory(0);
1200   fh2FFZBgr->SetDirectory(0);  
1201   fh2FFXiBgr->SetDirectory(0); 
1202
1203   f.Close();  
1204
1205   // nJets per bin
1206
1207   for(Int_t i=0; i<fNJetPtSlices; i++){
1208     
1209     Float_t jetPtLoLim = fJetPtSlices->At(i);
1210     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1211     
1212     Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
1213     Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
1214                                      
1215     Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
1216     Double_t scaleF = 1;
1217
1218     fNJetsBgr->SetAt(nJetsBin*scaleF,i); 
1219
1220     if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
1221                         static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);
1222
1223   }
1224   
1225   // projections: FF 
1226   
1227   for(Int_t i=0; i<fNJetPtSlices; i++){
1228     
1229     Float_t jetPtLoLim = fJetPtSlices->At(i);
1230     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1231
1232     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
1233     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;
1234
1235     if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
1236       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
1237       return; 
1238     }
1239     
1240     TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1241     TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1242     TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
1243     
1244     // appendix 'unbinned' to avoid histos with same name after rebinning
1245     TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
1246     TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
1247     TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
1248     
1249     if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
1250     if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
1251     if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());
1252
1253     projPt->SetNameTitle(strNameBgrPt,"");
1254     projZ->SetNameTitle(strNameBgrZ,"");
1255     projXi->SetNameTitle(strNameBgrXi,"");
1256     
1257     // raw bgr = corr level 0
1258     fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
1259   }  
1260 }
1261
1262
1263 //__________________________________________________________________________________________________________
1264 void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
1265 {
1266   // write histos to file
1267   // skip histos with 0 entries
1268
1269   TString outfileOption = "RECREATE";
1270   if(updateOutfile) outfileOption = "UPDATE";
1271
1272   TFile f(strfile,outfileOption); 
1273
1274   if(!f.IsOpen()){
1275     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1276     return;
1277   }
1278
1279   if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1280
1281   if(strdir && strdir.Length()){
1282     TDirectory* dir = f.mkdir(strdir);
1283     dir->cd(); 
1284   } 
1285
1286   for(Int_t i=0; i<fNJetPtSlices; i++){
1287     
1288     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
1289     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries())       fCorrFF[c]->GetZ(i)->Write();
1290     for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries())      fCorrFF[c]->GetXi(i)->Write();
1291
1292     if(fh1FFXiShift[i]->GetEntries()) fh1FFXiShift[i]->Write();
1293
1294     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
1295     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries())       fCorrBgr[c]->GetZ(i)->Write();
1296     for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries())      fCorrBgr[c]->GetXi(i)->Write();
1297
1298     
1299     if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
1300     if(fh1FFZBackFolded[i]       && fh1FFZBackFolded[i]->GetEntries())       fh1FFZBackFolded[i]->Write();
1301     if(fh1FFXiBackFolded[i]      && fh1FFXiBackFolded[i]->GetEntries())      fh1FFXiBackFolded[i]->Write();
1302
1303
1304     if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
1305     if(fh1FFRatioZFolded[i]       && fh1FFRatioZFolded[i]->GetEntries())       fh1FFRatioZFolded[i]->Write();
1306     if(fh1FFRatioXiFolded[i]      && fh1FFRatioXiFolded[i]->GetEntries())      fh1FFRatioXiFolded[i]->Write();
1307
1308     if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
1309     if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries())             fh1FFRatioZBackFolded[i]->Write();
1310     if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries())            fh1FFRatioXiBackFolded[i]->Write();
1311
1312   }
1313   
1314   // inclusive track pt
1315
1316   for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
1317   if(fh1SingleTrackPtBackFolded)      fh1SingleTrackPtBackFolded->Write();  
1318   if(fh1RatioSingleTrackPtFolded)     fh1RatioSingleTrackPtFolded->Write();  
1319   if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();  
1320
1321   f.Close();  
1322 }
1323
1324 //____________________________________________________________________________________________________________________________________
1325 THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
1326 {
1327   // copy 1-dimensional histo to THnSparse 
1328   // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
1329   // histos with variable bin size are supported
1330
1331   // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user
1332
1333   THnSparseF* fhnHist;
1334
1335   Int_t nBins       = hist->GetXaxis()->GetNbins();
1336   Int_t nBinsVec[1] = { nBins };
1337
1338   const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
1339   
1340   if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1341
1342     fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
1343     fhnHist->SetBinEdges(0,binsVec);
1344   }
1345   else{ // uniform bin size
1346     
1347     Double_t xMin  = hist->GetXaxis()->GetXmin();
1348     Double_t xMax  = hist->GetXaxis()->GetXmax();
1349     
1350     Double_t binMinVec[1] = { xMin };
1351     Double_t binMaxVec[1] = { xMax };
1352     
1353     fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
1354   }
1355  
1356
1357   for(Int_t bin=1; bin<=nBins; bin++){ 
1358
1359     Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);
1360
1361     Double_t binCoord[] = {binCenter}; 
1362     fhnHist->Fill(binCoord,1); // initially need to create the bin 
1363
1364     Long64_t binIndex = fhnHist->GetBin(binCoord);
1365
1366     Double_t cont =  hist->GetBinContent(bin);
1367     Double_t err  =  hist->GetBinError(bin);
1368     
1369     if(fillConst){
1370       cont = 1;
1371       err  = 0;
1372     }
1373
1374     fhnHist->SetBinContent(binIndex,cont);
1375     fhnHist->SetBinError(binIndex,err);
1376  }
1377
1378  return fhnHist;
1379 }
1380
1381 //______________________________________________________________________________________________________________________________________________
1382 THnSparse* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
1383                                                        const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
1384 {
1385   // unfold input histo 
1386
1387   AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior);  // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
1388   unfolding.SetMaxNumberOfIterations(nIter);
1389   // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
1390   // if(useSmoothing) unfolding.UseSmoothing(); 
1391
1392   if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
1393   unfolding.Unfold();
1394
1395   THnSparse* unfolded = unfolding.GetUnfolded();
1396
1397   TString hnameUnf = hnHist->GetName();
1398   hnameUnf.Append("_unf");
1399   unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
1400
1401   return unfolded;
1402 }
1403
1404 //___________________________________________________________________________________________________________________________
1405 void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
1406 {
1407   // loop over jet pt slices and unfold dN/dpt spectra
1408   
1409   TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel(); 
1410   if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
1411   
1412   for(Int_t i=0; i<fNJetPtSlices; i++){
1413  
1414     TH1F* hist = 0;
1415     if(type == kFlagPt)      hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
1416     else if(type == kFlagZ)  hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i);       // level -2: before unfolding, level -1: unfolded
1417     else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i);      // level -2: before unfolding, level -1: unfolded
1418
1419     THnSparse* hnResponse = 0;
1420     if(type == kFlagPt) hnResponse = fhnResponsePt[i];
1421     else if(type == kFlagZ)  hnResponse = fhnResponseZ[i];
1422     else if(type == kFlagXi) hnResponse = fhnResponseXi[i];
1423
1424
1425     TH1F* hPrior = 0;
1426     if(type == kFlagPt && fh1FFTrackPtPrior[i]  && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
1427     else if(type == kFlagZ  && fh1FFZPrior[i]   && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0)       ) hPrior = fh1FFZPrior[i];
1428     else if(type == kFlagXi && fh1FFXiPrior[i]  && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0)      ) hPrior = fh1FFXiPrior[i];
1429
1430
1431     TString histNameTHn = hist->GetName();
1432     histNameTHn.ReplaceAll("TH1","THn");
1433
1434     TString priorNameTHn; 
1435     if(hPrior){
1436       priorNameTHn = hPrior->GetName();
1437       priorNameTHn.ReplaceAll("TH1","THn");
1438     }
1439
1440     TString histNameBackFolded = hist->GetName();
1441     histNameBackFolded.Append("_backfold");
1442
1443     TString histNameRatioFolded = hist->GetName();
1444     histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
1445     histNameRatioFolded.Append("_unfold");
1446
1447     TString histNameRatioBackFolded = hist->GetName();
1448     histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
1449     histNameRatioBackFolded.Append("_backfold");
1450  
1451     THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
1452     THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff 
1453     THnSparse* hnPrior          = 0;
1454     if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
1455
1456     THnSparse* hnUnfolded 
1457       = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);  
1458      
1459     TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
1460     hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
1461     
1462     if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
1463     if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
1464     if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);
1465
1466     // backfolding: apply response matrix to unfolded spectrum
1467     TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); 
1468     hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
1469
1470     if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
1471     if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
1472     if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
1473     
1474     // ratio unfolded to original histo 
1475     TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
1476     hRatioUnfolded->Reset();
1477     hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
1478
1479     if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
1480     if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
1481     if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;
1482
1483
1484     // ratio backfolded to original histo
1485     TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
1486     hRatioBackFolded->Reset();
1487     hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
1488
1489     if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
1490     if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
1491     if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
1492     
1493     delete hnHist;
1494     delete hnFlatEfficiency;
1495
1496  }
1497 }
1498
1499 //_____________________________________________________________________________________________________
1500 void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
1501 {
1502   
1503   Int_t type = kFlagPt;
1504   UnfoldHistos(nIter, useCorrelatedErrors, type);
1505 }
1506
1507 //_____________________________________________________________________________________________________
1508 void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
1509 {
1510   
1511   Int_t type = kFlagZ;
1512   UnfoldHistos(nIter, useCorrelatedErrors, type);
1513 }
1514
1515 //_____________________________________________________________________________________________________
1516 void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
1517 {
1518   
1519   Int_t type = kFlagXi;
1520   UnfoldHistos(nIter, useCorrelatedErrors, type);
1521 }
1522
1523 //______________________________________________________________________________________________
1524 TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
1525 {
1526   // apply (multiply) response matrix to hist 
1527
1528   TH1F* hOut = new TH1F(*hist);
1529   hOut->Reset();
1530
1531   const Int_t axisM = 0; 
1532   const Int_t axisT = 1;
1533  
1534   Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
1535   Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();
1536
1537   // response matrix normalization
1538   // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
1539   // so for unfolding proper we leave it to CORRFW ...
1540
1541   Double_t normFacResponse[nBinsT];
1542
1543   for(Int_t iT=1; iT<=nBinsT; iT++){
1544
1545     Double_t sumResp = 0;
1546     
1547     for(Int_t iM=1; iM<=nBinsM; iM++){
1548       
1549       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1550       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1551       
1552       Double_t binCoord[] = {coordM,coordT};
1553       
1554       Long64_t binIndex = hnResponse->GetBin(binCoord);
1555       
1556       Double_t resp = hnResponse->GetBinContent(binIndex); 
1557       
1558       sumResp += resp;
1559     }
1560     
1561     normFacResponse[iT] = 0;
1562     if(sumResp) normFacResponse[iT] = 1/sumResp;
1563   }
1564   
1565   
1566   
1567   for(Int_t iM=1; iM<=nBinsM; iM++){
1568     
1569     Double_t contM   = 0;
1570
1571     for(Int_t iT=1; iT<=nBinsT; iT++){
1572
1573       Double_t contT = hist->GetBinContent(iT);
1574       
1575       Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
1576       Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
1577
1578       Double_t binCoord[] = {coordM,coordT};
1579       
1580       Long64_t binIndex = hnResponse->GetBin(binCoord);
1581       
1582       Double_t resp = hnResponse->GetBinContent(binIndex); 
1583       
1584       contM   += resp*normFacResponse[iT]*contT; 
1585     }
1586
1587     hOut->SetBinContent(iM,contM);
1588   }
1589
1590   return hOut;
1591 }
1592
1593 //_______________________________________________________________________________________________________
1594 void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
1595 {
1596   // read reconstruction efficiency from file
1597   // argument strlist optional - read from directory strdir if not specified
1598
1599   // temporary histos to hold histos from file
1600   TH1F* hEffPt[fNJetPtSlices]; 
1601   TH1F* hEffZ[fNJetPtSlices];
1602   TH1F* hEffXi[fNJetPtSlices];
1603   
1604   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
1605   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
1606   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;
1607
1608   TFile f(strfile,"READ");
1609
1610   if(!f.IsOpen()){
1611     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1612     return;
1613   }
1614
1615   if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1616  
1617   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1618
1619   TList* list = 0;
1620
1621   if(strlist && strlist.Length()){
1622    
1623     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1624       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1625       return;
1626     }
1627   }  
1628
1629   for(Int_t i=0; i<fNJetPtSlices; i++){
1630     
1631     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1632     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1633     
1634     TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
1635     TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
1636     TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
1637      
1638    
1639     if(list){
1640       hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt); 
1641       hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ); 
1642       hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi); 
1643     }
1644     else{
1645       hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt); 
1646       hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ); 
1647       hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi); 
1648     }
1649     
1650     if(!hEffPt[i]){
1651       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
1652     }
1653   
1654     if(!hEffZ[i]){
1655       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
1656     }    
1657
1658     if(!hEffXi[i]){
1659       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
1660     }
1661
1662
1663     if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
1664     if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
1665     if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());
1666
1667     if(hEffPt[i]) hEffPt[i]->SetDirectory(0); 
1668     if(hEffZ[i])  hEffZ[i]->SetDirectory(0); 
1669     if(hEffXi[i]) hEffXi[i]->SetDirectory(0); 
1670
1671   } // jet slices loop
1672
1673   f.Close();
1674
1675   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1676     if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
1677     if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
1678     if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
1679   }
1680 }
1681
1682 //___________________________________________________________________________________________________________
1683 void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
1684 {
1685   // read bgr eff from file
1686   // argument strlist optional - read from directory strdir if not specified
1687  
1688   TH1F* hEffPtBgr[fNJetPtSlices]; 
1689   TH1F* hEffZBgr [fNJetPtSlices];
1690   TH1F* hEffXiBgr[fNJetPtSlices];
1691
1692   for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
1693   for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
1694   for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;
1695
1696
1697   TFile f(strfile,"READ");
1698
1699   if(!f.IsOpen()){
1700     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
1701     return;
1702   }
1703
1704   if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
1705  
1706   if(strdir && strdir.Length()) gDirectory->cd(strdir);
1707
1708   TList* list = 0;
1709
1710   if(strlist && strlist.Length()){
1711        
1712     if(!(list = (TList*) gDirectory->Get(strlist))){ 
1713       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
1714       return;
1715     }
1716   }  
1717
1718   for(Int_t i=0; i<fNJetPtSlices; i++){
1719     
1720     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
1721     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
1722     
1723     TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1724     TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1725     TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
1726   
1727    
1728     if(list){
1729       hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr); 
1730       hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr); 
1731       hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr); 
1732     }
1733     else{
1734       hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr); 
1735       hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr); 
1736       hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr); 
1737     }
1738     
1739     if(!hEffPtBgr[i]){
1740       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
1741     }
1742   
1743     if(!hEffZBgr[i]){
1744       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
1745     }    
1746
1747     if(!hEffXiBgr[i]){
1748       Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
1749     }
1750
1751
1752     if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
1753     if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
1754     if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());
1755
1756     if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0); 
1757     if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0); 
1758     if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0); 
1759
1760   } // jet slices loop
1761
1762   f.Close();
1763
1764   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
1765     if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
1766     if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
1767     if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
1768   }
1769 }
1770
1771 // ________________________________________________
1772 void AliFragmentationFunctionCorrections::EffCorr()
1773 {
1774   // apply efficiency correction
1775
1776   AddCorrectionLevel("eff");
1777
1778   for(Int_t i=0; i<fNJetPtSlices; i++){
1779
1780     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1781     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1782     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1783
1784     TString histNamePt = histPt->GetName();
1785     TString histNameZ  = histZ->GetName();
1786     TString histNameXi = histXi->GetName();
1787
1788     
1789     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1790     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1791     
1792     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1793     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1794     
1795     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1796     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1797     
1798     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1799   }
1800 }
1801
1802 //___________________________________________________
1803 void AliFragmentationFunctionCorrections::EffCorrBgr()
1804 {
1805   // apply efficiency correction to bgr distributions
1806
1807   AddCorrectionLevelBgr("eff");
1808
1809   Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
1810
1811   for(Int_t i=0; i<fNJetPtSlices; i++){
1812
1813     TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
1814     TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
1815     TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
1816     
1817     TString histNamePt = histPt->GetName();
1818     TString histNameZ  = histZ->GetName();
1819     TString histNameXi = histXi->GetName();
1820
1821     
1822     TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
1823     hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
1824     
1825     TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
1826     hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
1827     
1828     TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
1829     hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
1830     
1831     fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
1832   }
1833 }
1834
1835 //______________________________________________________________________
1836 void AliFragmentationFunctionCorrections::XiShift(const Int_t corrLevel)
1837 {
1838   // re-evaluate jet energy after FF corrections from dN/dpt distribution
1839   // 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'
1840   // argument corrlevel: which level of correction to be corrected/shifted to 
1841
1842   if(corrLevel>=fNCorrectionLevels){ 
1843     Printf(" calc xi shift: corrLevel exceeded - do nothing");
1844     return;
1845   }
1846
1847   Double_t* jetPtUncorr = new Double_t[fNJetPtSlices];
1848   Double_t* jetPtCorr   = new Double_t[fNJetPtSlices];
1849   Double_t* deltaXi     = new Double_t[fNJetPtSlices];
1850
1851   for(Int_t i=0; i<fNJetPtSlices; i++){
1852     
1853     TH1F* histPtRaw = fCorrFF[0]->GetTrackPt(i);
1854     TH1F* histPt    = fCorrFF[corrLevel]->GetTrackPt(i);
1855
1856     Double_t ptUncorr = 0;
1857     Double_t ptCorr   = 0;
1858
1859     for(Int_t bin = 1; bin<=histPtRaw->GetNbinsX(); bin++){
1860
1861       Double_t cont   = histPtRaw->GetBinContent(bin);
1862       Double_t width  = histPtRaw->GetBinWidth(bin);
1863       Double_t meanPt = histPtRaw->GetBinCenter(bin);
1864
1865       ptUncorr += meanPt*cont*width;
1866     }
1867     
1868     for(Int_t bin = 1; bin<=histPt->GetNbinsX(); bin++){
1869       
1870       Double_t cont   = histPt->GetBinContent(bin);
1871       Double_t width  = histPt->GetBinWidth(bin);
1872       Double_t meanPt = histPt->GetBinCenter(bin);
1873       
1874       ptCorr += meanPt*cont*width;
1875     }
1876
1877     jetPtUncorr[i] = ptUncorr; 
1878     jetPtCorr[i]   = ptCorr;   
1879   }
1880
1881   // calc dXi from dN/dpt distribution : 
1882   // sum over track pt for raw and corrected FF is equivalent to raw/corrected jet pt 
1883
1884   for(Int_t i=0; i<fNJetPtSlices; i++){
1885
1886     Float_t jetPtLoLim = fJetPtSlices->At(i);
1887     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
1888
1889     Double_t meanJetPt = 0.5*(jetPtUpLim+jetPtLoLim);
1890     
1891     Double_t ptUncorr = jetPtUncorr[i]; 
1892     Double_t ptCorr   = jetPtCorr[i]; 
1893
1894     Double_t dXi = TMath::Log(ptCorr/ptUncorr);
1895     
1896     Printf(" calc xi shift: jet pt slice %d, mean jet pt %f, ptUncorr %f, ptCorr %f, ratio corr/uncorr %f, dXi %f "
1897            ,i,meanJetPt,ptUncorr,ptCorr,ptCorr/ptUncorr,dXi);
1898     
1899     deltaXi[i] = dXi;
1900   }
1901   
1902   // book & fill new dN/dxi histos
1903
1904   for(Int_t i=0; i<fNJetPtSlices; i++){
1905
1906     TH1F* histXi = fCorrFF[corrLevel]->GetXi(i);
1907             
1908     Double_t dXi = deltaXi[i];
1909
1910     Int_t nBins  = histXi->GetNbinsX();
1911     const Double_t* binsVec = histXi->GetXaxis()->GetXbins()->GetArray();
1912     Float_t binsVecNew[nBins+1];
1913     
1914     TString strName = histXi->GetName(); 
1915     strName.Append("_shift");
1916     TString strTit  = histXi->GetTitle();  
1917
1918     TH1F* hXiCorr; 
1919
1920     // create shifted histo ...
1921
1922     if(binsVec){ // binsVec only neq NULL if histo was rebinned before
1923
1924       for(Int_t bin=0; bin<nBins+1; bin++) binsVecNew[bin] = binsVec[bin] + dXi;    
1925       hXiCorr = new TH1F(strName,strTit,nBins,binsVecNew);
1926     }
1927     else{ // uniform bin size
1928       
1929       Double_t xMin  = histXi->GetXaxis()->GetXmin();
1930       Double_t xMax  = histXi->GetXaxis()->GetXmax();
1931
1932       xMin += dXi;
1933       xMax += dXi;
1934       
1935       hXiCorr = new TH1F(strName,strTit,nBins,xMin,xMax);
1936     }
1937
1938     // ... and fill
1939
1940     for(Int_t bin=1; bin<nBins+1; bin++){
1941       Double_t cont = histXi->GetBinContent(bin);
1942       Double_t err  = histXi->GetBinError(bin);
1943       
1944       hXiCorr->SetBinContent(bin,cont);
1945       hXiCorr->SetBinError(bin,err);
1946     }
1947     
1948     new(fh1FFXiShift[i]) TH1F(*hXiCorr);
1949     delete hXiCorr;
1950   }
1951
1952   delete[] jetPtUncorr;
1953   delete[] jetPtCorr;
1954   delete[] deltaXi;
1955 }
1956
1957 //_____________________________________________________
1958 void AliFragmentationFunctionCorrections::SubtractBgr()
1959 {
1960   // subtract bgr distribution from FF
1961   // the current corr level is used for both 
1962   
1963   AddCorrectionLevel("bgrSub");
1964
1965   for(Int_t i=0; i<fNJetPtSlices; i++){
1966
1967     TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
1968     TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
1969     TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
1970     
1971     TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
1972     TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
1973     TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);
1974
1975     TString histNamePt = histPt->GetName();
1976     TString histNameZ  = histZ->GetName();
1977     TString histNameXi = histXi->GetName();
1978     
1979
1980     TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
1981     hFFTrackPtBgrSub->Add(histPtBgr,-1);
1982     
1983     TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
1984     hFFZBgrSub->Add(histZBgr,-1);
1985     
1986     TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
1987     hFFXiBgrSub->Add(histXiBgr,-1);
1988
1989     fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);
1990   }
1991 }
1992
1993 //________________________________________________________________________________________________________________
1994 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile, 
1995                                                               Bool_t updateOutfile, TString strOutDir,TString strPostfix)
1996
1997   // read task ouput from MC and write single track eff - standard dir/list 
1998      
1999   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2000   TString strlist = "fracfunc_" + strID;
2001     
2002   WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
2003 }
2004
2005 //___________________________________________________________________________________________________________________________________
2006 void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, 
2007                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
2008 {
2009   // read task output from MC and write single track eff as function of pt, eta, phi
2010   // argument strlist optional - read from directory strdir if not specified
2011   // write eff to file stroutfile - by default only 'update' file (don't overwrite)
2012
2013
2014   TH1D* hdNdptTracksMCPrimGen;
2015   TH2D* hdNdetadphiTracksMCPrimGen;
2016   
2017   TH1D* hdNdptTracksMCPrimRec;
2018   TH2D* hdNdetadphiTracksMCPrimRec;
2019     
2020
2021   TFile f(strInfile,"READ");
2022
2023   if(!f.IsOpen()){
2024     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2025     return;
2026   }
2027   
2028   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2029  
2030   if(strdir && strdir.Length()){
2031     if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
2032     gDirectory->cd(strdir);
2033   }
2034
2035   TList* list = 0;
2036
2037   if(strlist && strlist.Length()){
2038
2039     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2040       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2041       return;
2042     }
2043   }
2044
2045
2046   TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
2047   if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data()); 
2048
2049   TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
2050   if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());
2051
2052   TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
2053   if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
2054   
2055   TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
2056   if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
2057   
2058
2059   if(list){
2060     hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
2061     hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
2062     
2063     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
2064     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
2065   }
2066   else{
2067     hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
2068     hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);
2069
2070     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
2071     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
2072   }
2073
2074   hdNdptTracksMCPrimGen->SetDirectory(0);
2075   hdNdetadphiTracksMCPrimGen->SetDirectory(0);
2076   hdNdptTracksMCPrimRec->SetDirectory(0);
2077   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2078   
2079   f.Close();
2080
2081   // projections: dN/deta, dN/dphi 
2082
2083   TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
2084   TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
2085  
2086   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2087   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2088
2089   // rebin
2090
2091   TString strNamePtGen = "hTrackPtGenPrim";
2092   TString strNamePtRec = "hTrackPtRecPrim";
2093
2094   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
2095   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
2096  
2097     // efficiency: divide 
2098
2099   TString hNameTrackEffPt = "hSingleTrackEffPt";
2100   if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
2101                                                
2102   TString hNameTrackEffEta = "hSingleTrackEffEta";
2103   if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());
2104
2105   TString hNameTrackEffPhi = "hSingleTrackEffPhi";
2106   if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());
2107
2108
2109   TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
2110   hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors
2111
2112   TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
2113   hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
2114   
2115   TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
2116   hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
2117   
2118   
2119   TString outfileOption = "RECREATE";
2120   if(updateOutfile)  outfileOption = "UPDATE";
2121
2122   TFile out(strOutfile,outfileOption);
2123
2124   if(!out.IsOpen()){
2125     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2126     return;
2127   }
2128
2129   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2130
2131   if(strOutDir && strOutDir.Length()){
2132     
2133     TDirectory* dir;
2134     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2135     else{
2136       dir = out.mkdir(strOutDir);
2137       dir->cd(); 
2138     } 
2139   }
2140
2141   hSingleTrackEffPt->Write();
2142   hSingleTrackEffEta->Write();
2143   hSingleTrackEffPhi->Write();
2144   
2145   out.Close();
2146 }
2147
2148 //________________________________________________________________________________________________________________
2149 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile, 
2150                                                                   Bool_t updateOutfile, TString strOutDir)
2151
2152   // read task ouput from MC and write single track eff - standard dir/list 
2153      
2154   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2155   TString strlist = "fracfunc_" + strID;
2156     
2157   WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2158 }
2159
2160 //___________________________________________________________________________________________________________________________________
2161 void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, 
2162                                                                   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2163 {
2164   // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
2165   // argument strlist optional - read from directory strdir if not specified
2166   // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)
2167
2168   TH1D* hdNdptTracksMCPrimRec;
2169   TH2D* hdNdetadphiTracksMCPrimRec;
2170   
2171   TH1D* hdNdptTracksMCSecRec;
2172   TH2D* hdNdetadphiTracksMCSecRec;
2173
2174   TFile f(strInfile,"READ");
2175
2176   if(!f.IsOpen()){
2177     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2178     return;
2179   }
2180   
2181   if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2182  
2183   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2184
2185   TList* list = 0;
2186
2187   if(strlist && strlist.Length()){
2188
2189     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2190       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2191       return;
2192     }
2193   }
2194
2195
2196   if(list){
2197     hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
2198     hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
2199     
2200     hdNdptTracksMCSecRec       = (TH1D*) list->FindObject("fh1TrackQAPtSecRec");
2201     hdNdetadphiTracksMCSecRec  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRec");
2202   }
2203   else{
2204     hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
2205     hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");
2206
2207     hdNdptTracksMCSecRec       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRec");
2208     hdNdetadphiTracksMCSecRec  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRec");
2209   }
2210   
2211   hdNdptTracksMCPrimRec->SetDirectory(0);
2212   hdNdetadphiTracksMCPrimRec->SetDirectory(0);
2213   hdNdptTracksMCSecRec->SetDirectory(0);
2214   hdNdetadphiTracksMCSecRec->SetDirectory(0);
2215   
2216   f.Close();
2217
2218   // projections: dN/deta, dN/dphi 
2219
2220   TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
2221   TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
2222  
2223   TH1D* hdNdetaTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionX("hdNdetaTracksMcSecRec");
2224   TH1D* hdNdphiTracksMCSecRec = (TH1D*) hdNdetadphiTracksMCSecRec->ProjectionY("hdNdphiTracksMcSecRec");
2225
2226  
2227   // rebin
2228
2229   TString strNamePtPrim = "hTrackPtPrim";
2230   TString strNamePtSec  = "hTrackPtSec";
2231
2232   if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
2233   if(fNHistoBinsSinglePt) hdNdptTracksMCSecRec  = (TH1D*) hdNdptTracksMCSecRec->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
2234
2235
2236   // secondary correction factor: divide prim/(prim+sec)
2237
2238   TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSingleTrackSecCorrPt");
2239   TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec->Clone("hSumPrimSecPt");
2240   hSumPrimSecPt->Add(hdNdptTracksMCPrimRec);
2241   hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors
2242
2243   TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSingleTrackSecCorrEta");
2244   TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCSecRec->Clone("hSumPrimSecEta");
2245   hSumPrimSecEta->Add(hdNdetaTracksMCPrimRec);
2246   hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors
2247
2248   TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSingleTrackSecCorrPhi");
2249   TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCSecRec->Clone("hSumPrimSecPhi");
2250   hSumPrimSecPhi->Add(hdNdphiTracksMCPrimRec);
2251   hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors
2252   
2253   
2254   TString outfileOption = "RECREATE";
2255   if(updateOutfile)  outfileOption = "UPDATE";
2256
2257   TFile out(strOutfile,outfileOption);
2258
2259   if(!out.IsOpen()){
2260     Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2261     return;
2262   }
2263
2264   if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2265
2266   if(strOutDir && strOutDir.Length()){  
2267
2268     TDirectory* dir;
2269     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
2270     else{
2271       dir = out.mkdir(strOutDir);
2272       dir->cd(); 
2273     } 
2274   } 
2275
2276   hdNdptTracksMCSecRec->Write();
2277   hdNdetaTracksMCSecRec->Write();
2278   hdNdphiTracksMCSecRec->Write();
2279
2280   hSingleTrackSecCorrPt->Write();
2281   hSingleTrackSecCorrEta->Write();
2282   hSingleTrackSecCorrPhi->Write();
2283   
2284   out.Close();
2285 }
2286
2287 //________________________________________________________________________________________________________________
2288 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile, 
2289                                                               Bool_t updateOutfile, TString strOutDir)
2290
2291   // read task ouput from MC and write single track eff - standard dir/list 
2292      
2293   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2294   TString strlist = "fracfunc_" + strID;
2295     
2296   WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
2297 }
2298
2299
2300 //_____________________________________________________________________________________________________________________________________
2301 void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
2302                                                               TString strOutfile, Bool_t updateOutfile, TString strOutDir)
2303 {
2304   // read 2d THnSparse response matrices in pt from file
2305   // project TH2 
2306   // write to strOutfile 
2307
2308   THnSparse* hnResponseSinglePt; 
2309   TH2F*      h2ResponseSinglePt; 
2310  
2311   TFile f(strInfile,"READ");
2312
2313   if(!f.IsOpen()){
2314     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2315     return;
2316   }
2317
2318   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2319  
2320   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2321
2322   TList* list = 0;
2323
2324   if(strlist && strlist.Length()){
2325     
2326     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2327       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2328       return;
2329     }
2330   }  
2331   
2332   if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
2333   else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
2334   
2335
2336   if(!hnResponseSinglePt){
2337     Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
2338     return;
2339   }
2340
2341   f.Close();
2342
2343
2344   // project 2d histo 
2345   h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
2346   h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt",""); 
2347     
2348   
2349   // write 
2350
2351   TString outfileOption = "RECREATE";
2352   if(updateOutfile)  outfileOption = "UPDATE";
2353   
2354   TFile out(strOutfile,outfileOption);
2355   
2356   if(!out.IsOpen()){
2357     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2358     return;
2359   }
2360
2361   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2362   
2363   if(strOutDir && strOutDir.Length()){  
2364
2365     TDirectory* dir;
2366     if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
2367     else{
2368       dir = out.mkdir(strOutDir);
2369       dir->cd(); 
2370     } 
2371   }
2372   
2373   hnResponseSinglePt->Write();
2374   h2ResponseSinglePt->Write();
2375   
2376   out.Close();  
2377 }
2378
2379 //________________________________________________________________________________________________________________
2380 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile, 
2381                                                            Bool_t updateOutfile)
2382
2383   // read task ouput from MC and write single track eff - standard dir/list 
2384      
2385   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2386   TString strlist = "fracfunc_" + strID;
2387     
2388   WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
2389 }
2390
2391 //___________________________________________________________________________________________________________________________________
2392 void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, 
2393                                                            TString strOutfile, Bool_t updateOutfile)
2394 {
2395   // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
2396   // argument strlist optional - read from directory strdir if not specified
2397   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2398
2399   TH1F* hEffPt[fNJetPtSlices];
2400   TH1F* hEffXi[fNJetPtSlices];
2401   TH1F* hEffZ[fNJetPtSlices];
2402
2403   TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
2404   TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
2405   TH1F* hdNdzMCPrimGen[fNJetPtSlices];
2406     
2407   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2408   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2409   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2410
2411
2412   TH1F* fh1FFJetPtRecEffGen;
2413
2414   TH2F* fh2FFTrackPtRecEffGen;
2415   TH2F* fh2FFZRecEffGen;
2416   TH2F* fh2FFXiRecEffGen;
2417   
2418   TH2F* fh2FFTrackPtRecEffRec;
2419   TH2F* fh2FFZRecEffRec;
2420   TH2F* fh2FFXiRecEffRec;
2421  
2422
2423   TFile f(strInfile,"READ");
2424
2425   if(!f.IsOpen()){
2426     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2427     return;
2428   }
2429   
2430   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2431  
2432   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2433
2434   TList* list = 0;
2435
2436   if(strlist && strlist.Length()){
2437
2438     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2439       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2440       return;
2441     }
2442   }
2443
2444   if(list){
2445     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2446
2447     fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
2448     fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
2449     fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
2450     
2451     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2452     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2453     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2454   }
2455   else{
2456     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2457
2458     fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffGen");
2459     fh2FFZRecEffGen        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffGen");
2460     fh2FFXiRecEffGen       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffGen");
2461     
2462     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2463     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2464     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2465   }
2466   
2467   fh1FFJetPtRecEffGen->SetDirectory(0); 
2468
2469   fh2FFTrackPtRecEffGen->SetDirectory(0);
2470   fh2FFZRecEffGen->SetDirectory(0);
2471   fh2FFXiRecEffGen->SetDirectory(0);
2472   
2473   fh2FFTrackPtRecEffRec->SetDirectory(0);
2474   fh2FFZRecEffRec->SetDirectory(0);
2475   fh2FFXiRecEffRec->SetDirectory(0);
2476   
2477   f.Close();
2478
2479
2480   // projections: FF for generated and reconstructed primaries 
2481   
2482   for(Int_t i=0; i<fNJetPtSlices; i++){
2483     
2484     Float_t jetPtLoLim = fJetPtSlices->At(i);
2485     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2486
2487     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
2488     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;
2489
2490     if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
2491       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2492       return; 
2493     }
2494     
2495     TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2496     TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2497     TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2498     
2499     TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2500     TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2501     TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2502     
2503     // project 
2504     // appendix 'unbinned' to avoid histos with same name after rebinning
2505
2506     hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2507     hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
2508     hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
2509     
2510     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2511     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
2512     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
2513     
2514     // rebin
2515
2516     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
2517     if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
2518     if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());
2519
2520     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
2521     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
2522     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());
2523
2524     hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
2525     hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
2526     hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
2527     
2528     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
2529     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
2530     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
2531  
2532     // normalize
2533     
2534     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2535
2536     NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin); 
2537     NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin); 
2538     NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin); 
2539
2540     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2541     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2542     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2543     
2544     // divide rec/gen : efficiency
2545
2546     TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2547     TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2548     TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2549  
2550     hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
2551     hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
2552     
2553     hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
2554     hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
2555     
2556     hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
2557     hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
2558   } 
2559   
2560   // write 
2561
2562   TString outfileOption = "RECREATE";
2563   if(updateOutfile)  outfileOption = "UPDATE";
2564
2565   TFile out(strOutfile,outfileOption);
2566   
2567   if(!out.IsOpen()){
2568     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2569     return;
2570   }
2571
2572   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2573
2574 //   if(strdir && strdir.Length()){ 
2575 //     TDirectory* dir = out.mkdir(strdir);
2576 //     dir->cd(); 
2577 //   }
2578
2579   for(Int_t i=0; i<fNJetPtSlices; i++){
2580
2581     hdNdptTracksMCPrimGen[i]->Write();
2582     hdNdxiMCPrimGen[i]->Write();
2583     hdNdzMCPrimGen[i]->Write();
2584     
2585     hdNdptTracksMCPrimRec[i]->Write();
2586     hdNdxiMCPrimRec[i]->Write();
2587     hdNdzMCPrimRec[i]->Write();
2588   
2589     hEffPt[i]->Write();
2590     hEffXi[i]->Write();
2591     hEffZ[i]->Write();
2592   }
2593
2594   out.Close();
2595
2596 }
2597
2598 //________________________________________________________________________________________________________________
2599 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile, 
2600                                                            Bool_t updateOutfile)
2601
2602   // read task ouput from MC and write secondary correction - standard dir/list 
2603      
2604   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2605   TString strlist = "fracfunc_" + strID;
2606     
2607   WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile);
2608 }
2609
2610 //___________________________________________________________________________________________________________________________________
2611 void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, 
2612                                                            TString strOutfile, Bool_t updateOutfile)
2613 {
2614   // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
2615   // argument strlist optional - read from directory strdir if not specified
2616   // write eff to file strOutfile - by default only 'update' file (don't overwrite)
2617
2618   TH1F* hSecCorrPt[fNJetPtSlices];
2619   TH1F* hSecCorrXi[fNJetPtSlices];
2620   TH1F* hSecCorrZ[fNJetPtSlices];
2621
2622   TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
2623   TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
2624   TH1F* hdNdzMCPrimRec[fNJetPtSlices];
2625     
2626   TH1F* hdNdptTracksMCSecRec[fNJetPtSlices];
2627   TH1F* hdNdxiMCSecRec[fNJetPtSlices];
2628   TH1F* hdNdzMCSecRec[fNJetPtSlices];
2629
2630   TH1F* fh1FFJetPtRecEffGen;
2631
2632   TH2F* fh2FFTrackPtRecEffRec;
2633   TH2F* fh2FFZRecEffRec;
2634   TH2F* fh2FFXiRecEffRec;
2635   
2636   TH2F* fh2FFTrackPtSecRec;
2637   TH2F* fh2FFZSecRec;
2638   TH2F* fh2FFXiSecRec;
2639
2640   TFile f(strInfile,"READ");
2641
2642   if(!f.IsOpen()){
2643     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2644     return;
2645   }
2646   
2647   if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
2648  
2649   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2650
2651   TList* list = 0;
2652
2653   if(strlist && strlist.Length()){
2654
2655     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2656       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2657       return;
2658     }
2659   }
2660
2661   if(list){
2662     fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");
2663
2664     fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
2665     fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
2666     fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
2667     
2668     fh2FFTrackPtSecRec  = (TH2F*) list->FindObject("fh2FFTrackPtSecRec");
2669     fh2FFZSecRec        = (TH2F*) list->FindObject("fh2FFZSecRec");
2670     fh2FFXiSecRec       = (TH2F*) list->FindObject("fh2FFXiSecRec");
2671   }
2672   else{
2673     fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->FindObject("fh1FFJetPtRecEffGen");
2674
2675     fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtRecEffRec");
2676     fh2FFZRecEffRec        = (TH2F*) gDirectory->FindObject("fh2FFZRecEffRec");
2677     fh2FFXiRecEffRec       = (TH2F*) gDirectory->FindObject("fh2FFXiRecEffRec");
2678     
2679     fh2FFTrackPtSecRec  = (TH2F*) gDirectory->FindObject("fh2FFTrackPtSecRec");
2680     fh2FFZSecRec        = (TH2F*) gDirectory->FindObject("fh2FFZSecRec");
2681     fh2FFXiSecRec       = (TH2F*) gDirectory->FindObject("fh2FFXiSecRec");
2682   }
2683   
2684   fh1FFJetPtRecEffGen->SetDirectory(0); 
2685
2686   fh2FFTrackPtRecEffRec->SetDirectory(0);
2687   fh2FFZRecEffRec->SetDirectory(0);
2688   fh2FFXiRecEffRec->SetDirectory(0);
2689   
2690   fh2FFTrackPtSecRec->SetDirectory(0);
2691   fh2FFZSecRec->SetDirectory(0);
2692   fh2FFXiSecRec->SetDirectory(0);
2693   
2694   f.Close();
2695
2696
2697   // projections: FF for generated and reconstructed primaries 
2698   
2699   for(Int_t i=0; i<fNJetPtSlices; i++){
2700     
2701     Float_t jetPtLoLim = fJetPtSlices->At(i);
2702     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2703
2704     Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
2705     Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;
2706
2707     if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
2708       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2709       return; 
2710     }
2711     
2712     TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2713     TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2714     TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2715     
2716     TString strNameFFPtSecRec(Form("fh1FFTrackPtRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2717     TString strNameFFZSecRec(Form("fh1FFZRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2718     TString strNameFFXiSecRec(Form("fh1FFXiRecSec_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2719     
2720     // project 
2721     // appendix 'unbinned' to avoid histos with same name after rebinning
2722
2723     hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2724     hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
2725     hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
2726     
2727     hdNdptTracksMCSecRec[i]  = (TH1F*) fh2FFTrackPtSecRec->ProjectionY(strNameFFPtSecRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
2728     hdNdzMCSecRec[i]         = (TH1F*) fh2FFZSecRec->ProjectionY(strNameFFZSecRec+"_unBinned",binLo,binUp,"o");
2729     hdNdxiMCSecRec[i]        = (TH1F*) fh2FFXiSecRec->ProjectionY(strNameFFXiSecRec+"_unBinned",binLo,binUp,"o");
2730     
2731     // rebin
2732
2733     if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
2734     if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
2735     if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());
2736
2737     if(fNHistoBinsPt[i]) hdNdptTracksMCSecRec[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRec,fHistoBinsPt[i]->GetArray());
2738     if(fNHistoBinsZ[i])  hdNdzMCSecRec[i]        = (TH1F*) hdNdzMCSecRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRec,fHistoBinsZ[i]->GetArray());
2739     if(fNHistoBinsXi[i]) hdNdxiMCSecRec[i]       = (TH1F*) hdNdxiMCSecRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRec,fHistoBinsXi[i]->GetArray());
2740
2741     hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
2742     hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
2743     hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
2744     
2745     hdNdptTracksMCSecRec[i]->SetNameTitle(strNameFFPtSecRec,"");
2746     hdNdzMCSecRec[i]->SetNameTitle(strNameFFZSecRec,"");
2747     hdNdxiMCSecRec[i]->SetNameTitle(strNameFFXiSecRec,"");
2748     
2749     // normalize
2750     
2751     Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);
2752
2753     NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
2754     NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
2755     NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
2756
2757     NormalizeTH1(hdNdptTracksMCSecRec[i],nJetsBin); 
2758     NormalizeTH1(hdNdzMCSecRec[i],nJetsBin); 
2759     NormalizeTH1(hdNdxiMCSecRec[i],nJetsBin); 
2760     
2761     // divide rec/gen : efficiency
2762
2763     TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2764     TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2765     TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2766  
2767     hSecCorrPt[i] = (TH1F*) hdNdptTracksMCSecRec[i]->Clone(strNameSecCorrPt);
2768     TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCSecRec[i]->Clone("hSumPrimSecPt");
2769     hSumPrimSecPt->Add(hdNdptTracksMCPrimRec[i]);
2770     hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors
2771
2772     hSecCorrXi[i] = (TH1F*) hdNdxiMCSecRec[i]->Clone(strNameSecCorrXi);
2773     TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCSecRec[i]->Clone("hSumPrimSecXi");
2774     hSumPrimSecXi->Add(hdNdxiMCPrimRec[i]);
2775     hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors
2776
2777     hSecCorrZ[i] = (TH1F*) hdNdzMCSecRec[i]->Clone(strNameSecCorrZ);
2778     TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCSecRec[i]->Clone("hSumPrimSecZ");
2779     hSumPrimSecZ->Add(hdNdzMCPrimRec[i]);
2780     hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors
2781   } 
2782   
2783   // write 
2784
2785   TString outfileOption = "RECREATE";
2786   if(updateOutfile)  outfileOption = "UPDATE";
2787
2788   TFile out(strOutfile,outfileOption);
2789   
2790   if(!out.IsOpen()){
2791     Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
2792     return;
2793   }
2794   
2795   if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
2796
2797   for(Int_t i=0; i<fNJetPtSlices; i++){
2798
2799     //  hdNdptTracksMCSecRec[i]->Write();
2800     //  hdNdxiMCSecRec[i]->Write();
2801     //  hdNdzMCSecRec[i]->Write();
2802   
2803     hSecCorrPt[i]->Write();
2804     hSecCorrXi[i]->Write();
2805     hSecCorrZ[i]->Write();
2806   }
2807
2808   out.Close();
2809 }
2810
2811 //________________________________________________________________________________________________________________
2812 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile, 
2813                                                            Bool_t updateOutfile)
2814
2815   // read task ouput from MC and write single track eff - standard dir/list 
2816      
2817   TString strdir  = "PWG4_FragmentationFunction_" + strID;
2818   TString strlist = "fracfunc_" + strID;
2819     
2820   WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile);
2821 }
2822
2823 //_____________________________________________________________________________________________________________________________________
2824 void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
2825                                                            TString strOutfile, Bool_t updateOutfile)
2826 {
2827   // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
2828   // project THnSparse + TH2 in jet pt slices 
2829   // write to strOutfile 
2830
2831   THnSparse* hn3ResponseJetPt;
2832   THnSparse* hn3ResponseJetZ;
2833   THnSparse* hn3ResponseJetXi;
2834
2835   // 2D response matrices 
2836
2837   THnSparse* hnResponsePt[fNJetPtSlices];
2838   THnSparse* hnResponseZ[fNJetPtSlices];
2839   THnSparse* hnResponseXi[fNJetPtSlices];
2840
2841   TH2F* h2ResponsePt[fNJetPtSlices];
2842   TH2F* h2ResponseZ[fNJetPtSlices];
2843   TH2F* h2ResponseXi[fNJetPtSlices];
2844
2845   // 1D projections on gen pt / rec pt axes
2846
2847   TH1F* h1FFPtRec[fNJetPtSlices]; 
2848   TH1F* h1FFZRec[fNJetPtSlices];
2849   TH1F* h1FFXiRec[fNJetPtSlices];
2850
2851   TH1F* h1FFPtGen[fNJetPtSlices]; 
2852   TH1F* h1FFZGen[fNJetPtSlices];
2853   TH1F* h1FFXiGen[fNJetPtSlices];
2854
2855   TH1F* h1RatioPt[fNJetPtSlices]; 
2856   TH1F* h1RatioZ[fNJetPtSlices];
2857   TH1F* h1RatioXi[fNJetPtSlices];
2858
2859
2860
2861   TFile f(strInfile,"READ");
2862
2863   if(!f.IsOpen()){
2864     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
2865     return;
2866   }
2867
2868   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
2869  
2870   if(strdir && strdir.Length()) gDirectory->cd(strdir);
2871
2872   TList* list = 0;
2873
2874   if(strlist && strlist.Length()){
2875     
2876     if(!(list = (TList*) gDirectory->Get(strlist))){ 
2877       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
2878       return;
2879     }
2880   }  
2881   
2882   if(list){
2883     hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
2884     hn3ResponseJetZ  = (THnSparse*) list->FindObject("fhnResponseJetZ");
2885     hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
2886   }
2887   else{
2888     hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
2889     hn3ResponseJetZ  = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
2890     hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
2891   }
2892
2893   
2894   if(!hn3ResponseJetPt){
2895     Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
2896     return;
2897   }
2898
2899   if(!hn3ResponseJetZ){
2900     Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
2901     return;
2902   }
2903
2904   if(!hn3ResponseJetXi){
2905     Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
2906     return;
2907   }
2908
2909   f.Close();  
2910
2911   // axes 
2912
2913   Int_t axisJetPtTHn3 = -1;
2914   Int_t axisGenPtTHn3 = -1;
2915   Int_t axisRecPtTHn3 = -1;
2916
2917   for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
2918     
2919     TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle(); 
2920
2921     if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i; 
2922     if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i; 
2923     if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i; 
2924   }
2925
2926   if(axisJetPtTHn3 == -1){
2927     Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
2928     return;
2929   }
2930
2931   if(axisGenPtTHn3 == -1){
2932     Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
2933     return;
2934   }
2935
2936   if(axisRecPtTHn3 == -1){
2937     Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
2938     return;
2939   }
2940
2941   for(Int_t i=0; i<fNJetPtSlices; i++){
2942     
2943     Float_t jetPtLoLim = fJetPtSlices->At(i);
2944     Float_t jetPtUpLim = fJetPtSlices->At(i+1);
2945
2946     Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
2947     Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;
2948
2949     if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
2950       Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
2951       return; 
2952     }
2953     
2954     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2955     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2956     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2957
2958     TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2959     TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2960     TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
2961          
2962     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2963     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2964     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2965  
2966     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2967     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2968     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2969  
2970     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2971     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2972     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
2973  
2974     
2975     // 2D projections in jet pt range
2976   
2977     hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
2978     hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
2979     hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
2980
2981     Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};  
2982
2983     hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
2984     hnResponseZ[i]  = hn3ResponseJetZ->Projection(2,axesProj);
2985     hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
2986   
2987     hnResponsePt[i]->SetNameTitle(strNameRespPt,""); 
2988     hnResponseZ[i]->SetNameTitle(strNameRespZ,""); 
2989     hnResponseXi[i]->SetNameTitle(strNameRespXi,""); 
2990
2991     h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
2992     h2ResponseZ[i]  = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
2993     h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
2994  
2995     h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,""); 
2996     h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,""); 
2997     h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,""); 
2998
2999
3000     // 1D projections
3001
3002     Int_t axisGenPtTHn2 = -1;
3003     Int_t axisRecPtTHn2 = -1;
3004
3005     for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
3006     
3007       TString title = hnResponsePt[i]->GetAxis(d)->GetTitle(); 
3008
3009       if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d; 
3010       if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d; 
3011     }
3012
3013     
3014     if(axisGenPtTHn2 == -1){
3015       Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
3016       return;
3017     }
3018     
3019     if(axisRecPtTHn2 == -1){
3020       Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
3021       return;
3022     }
3023     
3024
3025     h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3026     h1FFZRec[i]  = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3027     h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
3028     
3029     h1FFPtRec[i]->SetNameTitle(strNameRecPt,""); 
3030     h1FFZRec[i]->SetNameTitle(strNameRecZ,""); 
3031     h1FFXiRec[i]->SetNameTitle(strNameRecXi,""); 
3032
3033     h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3034     h1FFZGen[i]  = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3035     h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
3036     
3037     h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
3038     h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
3039     h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");
3040
3041     // normalize 1D projections
3042
3043     if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
3044     if(fNHistoBinsZ[i])  h1FFZRec[i]  = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3045     if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3046     
3047     if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3048     if(fNHistoBinsZ[i])  h1FFZGen[i]  = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3049     if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3050        
3051     NormalizeTH1(h1FFPtRec[i],fNJets->At(i)); 
3052     NormalizeTH1(h1FFZRec[i],fNJets->At(i));
3053     NormalizeTH1(h1FFXiRec[i],fNJets->At(i));
3054
3055     NormalizeTH1(h1FFPtGen[i],fNJets->At(i)); 
3056     NormalizeTH1(h1FFZGen[i],fNJets->At(i));
3057     NormalizeTH1(h1FFXiGen[i],fNJets->At(i));
3058
3059     // ratios 1D projections
3060
3061     h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
3062     h1RatioPt[i]->Reset();
3063     h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");
3064
3065     h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
3066     h1RatioZ[i]->Reset();
3067     h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
3068     
3069     h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
3070     h1RatioXi[i]->Reset();
3071     h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
3072   }
3073   
3074   
3075   // write 
3076
3077   TString outfileOption = "RECREATE";
3078   if(updateOutfile)  outfileOption = "UPDATE";
3079   
3080   TFile out(strOutfile,outfileOption);
3081   
3082   if(!out.IsOpen()){
3083     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3084     return;
3085   }
3086
3087   if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
3088
3089   //if(strdir && strdir.Length()){ 
3090   //  TDirectory* dir = out.mkdir(strdir);
3091   //  dir->cd(); 
3092   //}
3093
3094   for(Int_t i=0; i<fNJetPtSlices; i++){
3095
3096     hnResponsePt[i]->Write();
3097     hnResponseXi[i]->Write();
3098     hnResponseZ[i]->Write();
3099
3100     h2ResponsePt[i]->Write();
3101     h2ResponseXi[i]->Write();
3102     h2ResponseZ[i]->Write();
3103
3104     h1FFPtRec[i]->Write(); 
3105     h1FFZRec[i]->Write();
3106     h1FFXiRec[i]->Write();
3107     
3108     h1FFPtGen[i]->Write(); 
3109     h1FFZGen[i]->Write();
3110     h1FFXiGen[i]->Write();
3111     
3112     h1RatioPt[i]->Write(); 
3113     h1RatioZ[i]->Write();
3114     h1RatioXi[i]->Write();
3115
3116   }
3117
3118   out.Close();  
3119 }
3120
3121 //______________________________________________________________________________________________________
3122 void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
3123 {
3124   // read response matrices from file
3125   // argument strlist optional - read from directory strdir if not specified
3126   // note: THnSparse are not rebinned 
3127
3128   THnSparse* hRespPt[fNJetPtSlices]; 
3129   THnSparse* hRespZ[fNJetPtSlices];
3130   THnSparse* hRespXi[fNJetPtSlices];
3131   
3132   for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
3133   for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i]  = 0;
3134   for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;
3135
3136   TFile f(strfile,"READ");
3137
3138   if(!f.IsOpen()){
3139     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3140     return;
3141   }
3142
3143   if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3144  
3145   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3146
3147   TList* list = 0;
3148
3149   if(strlist && strlist.Length()){
3150    
3151     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3152       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3153       return;
3154     }
3155   }  
3156
3157   for(Int_t i=0; i<fNJetPtSlices; i++){
3158     
3159     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3160     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3161     
3162     TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3163     TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
3164     TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
3165         
3166     if(list){
3167       hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt); 
3168       hRespZ[i]  = (THnSparse*) list->FindObject(strNameRespZ); 
3169       hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi); 
3170     }
3171     else{
3172       hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt); 
3173       hRespZ[i]  = (THnSparse*) gDirectory->Get(strNameRespZ); 
3174       hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi); 
3175     }
3176     
3177     if(!hRespPt[i]){
3178       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
3179     }
3180   
3181     if(!hRespZ[i]){
3182       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
3183     }    
3184
3185     if(!hRespXi[i]){
3186       Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
3187     }
3188    
3189     //    if(0){ // can't rebin THnSparse ...
3190     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
3191     //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
3192     
3193     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
3194     //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
3195     
3196     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
3197     //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
3198     //     }
3199    
3200  
3201   } // jet slices loop
3202
3203   f.Close();  // THnSparse pointers still valid even if file closed
3204
3205 //   for(Int_t i=0; i<fNJetPtSlices; i++){  // no copy c'tor ...
3206 //     if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
3207 //     if(hRespZ[i])  new(fhnResponseZ[i])  THnSparseF(*hRespZ[i]);
3208 //     if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
3209 //   }
3210
3211   for(Int_t i=0; i<fNJetPtSlices; i++){ 
3212     fhnResponsePt[i] = hRespPt[i];
3213     fhnResponseZ[i]  = hRespZ[i];
3214     fhnResponseXi[i] = hRespXi[i];
3215   }
3216 }
3217
3218 //______________________________________________________________________________________________________________________
3219 void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
3220 {
3221   // read priors from file: rec primaries, gen pt dist
3222
3223   if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3224   
3225   // temporary histos to store pointers from file
3226   TH1F* hist[fNJetPtSlices]; 
3227
3228   for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;
3229
3230   TFile f(strfile,"READ");
3231
3232   if(!f.IsOpen()){
3233     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3234     return;
3235   }
3236
3237   for(Int_t i=0; i<fNJetPtSlices; i++){
3238     
3239     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3240     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3241     
3242     TString strName;
3243     
3244     if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); 
3245     if(type == kFlagZ)  strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);       
3246     if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);      
3247     
3248     hist[i] = (TH1F*) gDirectory->Get(strName); 
3249     
3250     if(!hist[i]){
3251       Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
3252     }
3253     
3254
3255     //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());
3256
3257     if(hist[i]) hist[i]->SetDirectory(0); 
3258     
3259   } // jet slices loop
3260   
3261   f.Close();
3262
3263   
3264   for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
3265     if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]); 
3266     if(hist[i] && type == kFlagZ)  new(fh1FFZPrior[i])       TH1F(*hist[i]); 
3267     if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i])      TH1F(*hist[i]); 
3268   }
3269 }
3270
3271 //_____________________________________________________
3272 // void AliFragmentationFunctionCorrections::RatioRecGen()
3273 // {
3274 //   // create ratio reconstructed over generated FF
3275 //   // use current highest corrLevel 
3276   
3277 //   Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);
3278
3279 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3280  
3281 //     TH1F* histPtRec =  fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
3282 //     TH1F* histZRec  =  fCorrFF[fNCorrectionLevels-1]->GetZ(i);       // levels -1: latest corr level
3283 //     TH1F* histXiRec =  fCorrFF[fNCorrectionLevels-1]->GetXi(i);      // levels -1: latest corr level
3284     
3285 //     TH1F* histPtGen = fh1FFTrackPtGenPrim[i]; 
3286 //     TH1F* histZGen  = fh1FFZGenPrim[i];
3287 //     TH1F* histXiGen = fh1FFXiGenPrim[i];
3288
3289 //     TString histNamePt = histPtRec->GetName();
3290 //     TString histNameZ  = histZRec->GetName();
3291 //     TString histNameXi = histXiRec->GetName();
3292
3293 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3294 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3295 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");
3296
3297 //     // ratio 
3298 //     TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
3299 //     hRatioRecGenPt->Reset();
3300 //     hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");
3301
3302 //     TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
3303 //     hRatioRecGenZ->Reset();
3304 //     hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
3305     
3306 //     TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
3307 //     hRatioRecGenXi->Reset();
3308 //     hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");
3309
3310 //     new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
3311 //     new(fh1FFRatioRecGenZ[i])  TH1F(*hRatioRecGenZ);
3312 //     new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
3313 //   }
3314 // }
3315
3316 // //___________________________________________________________
3317 // void AliFragmentationFunctionCorrections::RatioRecPrimaries()
3318 // {
3319 //   // create ratio reconstructed tracks over reconstructed primaries 
3320 //   // use raw FF (corrLevel 0)
3321   
3322 //   Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);
3323
3324 //   for(Int_t i=0; i<fNJetPtSlices; i++){
3325  
3326 //     const Int_t corrLevel = 0; 
3327
3328 //     TH1F* histPtRec =  fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
3329 //     TH1F* histZRec  =  fCorrFF[corrLevel]->GetZ(i);       // levels -1: latest corr level
3330 //     TH1F* histXiRec =  fCorrFF[corrLevel]->GetXi(i);      // levels -1: latest corr level
3331     
3332 //     TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i]; 
3333 //     TH1F* histZRecPrim  = fh1FFZRecPrim[i];
3334 //     TH1F* histXiRecPrim = fh1FFXiRecPrim[i];
3335
3336 //     TString histNamePt = histPtRec->GetName();
3337 //     TString histNameZ  = histZRec->GetName();
3338 //     TString histNameXi = histXiRec->GetName();
3339
3340 //     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3341 //     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3342 //     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
3343
3344 //     // ratio 
3345 //     TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
3346 //     hRatioRecPrimPt->Reset();
3347 //     hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");
3348
3349 //     TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
3350 //     hRatioRecPrimZ->Reset();
3351 //     hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
3352     
3353 //     TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
3354 //     hRatioRecPrimXi->Reset();
3355 //     hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");
3356
3357
3358 //     new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
3359 //     new(fh1FFRatioRecPrimZ[i])  TH1F(*hRatioRecPrimZ);
3360 //     new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
3361 //   }
3362 // }
3363
3364 // __________________________________________________________________________________
3365 void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
3366 {
3367
3368   // project response matrices on both axes: 
3369   // FF for rec primaries, in terms of generated and reconstructed momentum
3370   // write FF and ratios to outFile
3371
3372   Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3373
3374   TH1F* hFFPtRec[fNJetPtSlices]; 
3375   TH1F* hFFZRec[fNJetPtSlices];
3376   TH1F* hFFXiRec[fNJetPtSlices];
3377
3378   TH1F* hFFPtGen[fNJetPtSlices]; 
3379   TH1F* hFFZGen[fNJetPtSlices];
3380   TH1F* hFFXiGen[fNJetPtSlices];
3381
3382   TH1F* hRatioPt[fNJetPtSlices]; 
3383   TH1F* hRatioZ[fNJetPtSlices];
3384   TH1F* hRatioXi[fNJetPtSlices];
3385
3386
3387   Int_t axisGenPt = 1;
3388   Int_t axisRecPt = 0;
3389
3390   for(Int_t i=0; i<fNJetPtSlices; i++){ 
3391
3392     Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
3393     Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
3394  
3395     TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3396     TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3397     TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3398  
3399     TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3400     TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3401     TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
3402  
3403     TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3404     TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3405     TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
3406  
3407
3408     hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
3409     hFFZRec[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
3410     hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
3411     
3412     hFFPtRec[i]->SetNameTitle(strNameRecPt,""); 
3413     hFFZRec[i]->SetNameTitle(strNameRecZ,""); 
3414     hFFXiRec[i]->SetNameTitle(strNameRecXi,""); 
3415
3416
3417     hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
3418     hFFZGen[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
3419     hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
3420     
3421     hFFPtGen[i]->SetNameTitle(strNameGenPt,""); 
3422     hFFZGen[i]->SetNameTitle(strNameGenZ,""); 
3423     hFFXiGen[i]->SetNameTitle(strNameGenXi,""); 
3424    
3425
3426     if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
3427     if(fNHistoBinsZ[i])  hFFZRec[i]   = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3428     if(fNHistoBinsXi[i]) hFFXiRec[i]  = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3429     
3430     if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
3431     if(fNHistoBinsZ[i])  hFFZGen[i]  = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
3432     if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
3433     
3434     NormalizeTH1(hFFPtGen[i],fNJets->At(i)); 
3435     NormalizeTH1(hFFZGen[i],fNJets->At(i));
3436     NormalizeTH1(hFFXiGen[i],fNJets->At(i));
3437    
3438     NormalizeTH1(hFFPtRec[i],fNJets->At(i)); 
3439     NormalizeTH1(hFFZRec[i],fNJets->At(i));
3440     NormalizeTH1(hFFXiRec[i],fNJets->At(i));
3441
3442
3443     hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
3444     hRatioPt[i]->Reset();
3445     hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");
3446
3447     hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
3448     hRatioZ[i]->Reset();
3449     hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
3450     
3451     hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
3452     hRatioXi[i]->Reset();
3453     hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
3454   }
3455     
3456
3457
3458   // write 
3459   
3460   TFile out(strOutfile,"RECREATE");
3461   
3462   if(!out.IsOpen()){
3463     Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3464     return;
3465   }
3466   
3467   for(Int_t i=0; i<fNJetPtSlices; i++){
3468     
3469     hFFPtRec[i]->Write(); 
3470     hFFZRec[i]->Write();
3471     hFFXiRec[i]->Write();
3472     
3473     hFFPtGen[i]->Write(); 
3474     hFFZGen[i]->Write();
3475     hFFXiGen[i]->Write();
3476     
3477     hRatioPt[i]->Write(); 
3478     hRatioZ[i]->Write();
3479     hRatioXi[i]->Write();
3480   }
3481
3482   out.Close();  
3483 }
3484
3485 // ____________________________________________________________________________________________________________________________
3486 void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir )
3487
3488   // project response matrix on both axes: 
3489   // pt spec for rec primaries, in terms of generated and reconstructed momentum
3490   // write spec and ratios to outFile
3491   
3492   Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
3493   
3494   TH1F* hSpecPtRec;  
3495   TH1F* hSpecPtGen; 
3496   TH1F* hRatioPt; 
3497   
3498   Int_t axisGenPt = 1;
3499   Int_t axisRecPt = 0;
3500   
3501   TString strNameRecPt   = "h1SpecTrackPtRecPrim_recPt";
3502   TString strNameGenPt   = "h1SpecTrackPtRecPrim_genPt";
3503   TString strNameRatioPt = "h1RatioTrackPtRecPrim";
3504     
3505   hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
3506   hSpecPtRec->SetNameTitle(strNameRecPt,""); 
3507   
3508   hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
3509   hSpecPtGen->SetNameTitle(strNameGenPt,""); 
3510   
3511   hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
3512   hRatioPt->Reset();
3513   hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
3514   
3515   TString outfileOption = "RECREATE";
3516   if(updateOutfile) outfileOption = "UPDATE";
3517
3518   TFile out(strOutfile,outfileOption);
3519   
3520   if(!out.IsOpen()){
3521     Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
3522     return;
3523   }
3524       
3525
3526   if(strOutDir && strOutDir.Length()){  
3527     
3528     TDirectory* dir;
3529     if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
3530     else{
3531       dir = out.mkdir(strOutDir);
3532       dir->cd(); 
3533     } 
3534   }
3535   
3536   hSpecPtRec->Write(); 
3537   hSpecPtGen->Write(); 
3538   hRatioPt->Write(); 
3539   
3540   out.Close();  
3541 }
3542
3543
3544 //__________________________________________________________________________________________________________________________________________________________________
3545 void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
3546
3547   // rebin histo, rescale bins according to new width
3548   // only correct for input histos with equal bin size
3549
3550   // args: jetPtSlice, type, use current corr level
3551  
3552   // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
3553   // array size of binsLimits: nBinsLimits 
3554   // array size of binsWidth: nBinsLimits-1 
3555   // binsLimits have to be in increasing order
3556   // if binning undefined for any slice, original binning will be kept
3557
3558   if(!fNJetPtSlices){
3559     Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
3560     return;
3561   }
3562
3563   if(jetPtSlice>=fNJetPtSlices){
3564     Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
3565     return;
3566   }
3567
3568
3569   Double_t binLimitMin = binsLimits[0];
3570   Double_t binLimitMax = binsLimits[nBinsLimits-1];
3571
3572   Double_t binLimit = binLimitMin; // start value 
3573   
3574   Int_t sizeUpperLim = 1000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1; -  only if first bin has smallest width, but not case for dN/dxi ...
3575   TArrayD binsArray(sizeUpperLim);
3576   Int_t nBins = 0; 
3577   binsArray.SetAt(binLimitMin,nBins++);
3578
3579   while(binLimit<binLimitMax && nBins<sizeUpperLim){
3580
3581     Int_t currentSlice = -1;
3582     for(Int_t i=0; i<nBinsLimits; i++){
3583       if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
3584     }
3585     
3586     Double_t currentBinWidth = binsWidth[currentSlice];
3587     binLimit += currentBinWidth;
3588
3589     binsArray.SetAt(binLimit,nBins++);
3590   }
3591   
3592
3593   TH1F* hist = 0;
3594   if(type == kFlagPt)            hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice); 
3595   else if(type == kFlagZ)        hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);       
3596   else if(type == kFlagXi)       hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);      
3597   else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0); 
3598   else{ 
3599     Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
3600     return;
3601   }
3602   
3603   Double_t binWidthNoRebin = hist->GetBinWidth(1);
3604
3605   Double_t* bins = binsArray.GetArray();
3606
3607   hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
3608  
3609   for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){
3610
3611     Double_t binWidthRebin = hist->GetBinWidth(bin);
3612     Double_t scaleF = binWidthNoRebin / binWidthRebin;
3613
3614     Double_t binCont  = hist->GetBinContent(bin);
3615     Double_t binErr   = hist->GetBinError(bin);
3616     
3617     binCont *= scaleF;
3618     binErr  *= scaleF;
3619
3620     hist->SetBinContent(bin,binCont);
3621     hist->SetBinError(bin,binErr);
3622   }
3623
3624
3625
3626   TH1F* temp = new TH1F(*hist);
3627
3628   if(type == kFlagPt)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
3629   if(type == kFlagZ)        fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
3630   if(type == kFlagXi)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
3631   if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);
3632
3633
3634   delete temp;
3635 }
3636 //__________________________________________________________________________________________________________________________________________________________________
3637 void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
3638
3639  
3640  if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
3641  
3642   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3643
3644   TList* list = 0;
3645
3646   if(strlist && strlist.Length()){
3647     
3648     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3649       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3650       return;
3651     }
3652   }  
3653   
3654   THnSparse* hn6ResponseJetPt  = (THnSparse*) list->FindObject("fhnCorrelation");
3655  
3656   Int_t axis6RecJetPt = 0;
3657   Int_t axis6GenJetPt = 3;
3658
3659   hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)"); 
3660   hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)"); 
3661   
3662   Int_t nBinsRecPt    = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins(); 
3663   Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
3664   Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);
3665
3666   Int_t nBinsGenPt    = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins(); 
3667   Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
3668   Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);
3669
3670   Int_t nBinsTrackPt = 200;
3671   Int_t loLimTrackPt = 0;
3672   Int_t upLimTrackPt = 200;
3673   
3674
3675   Int_t    nBinsResponse[4]  = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
3676   Double_t binMinResponse[4] = {loLimRecPt,loLimTrackPt,loLimGenPt,loLimTrackPt};
3677   Double_t binMaxResponse[4] = {upLimRecPt,upLimTrackPt,upLimGenPt,upLimTrackPt};
3678
3679   const char* labelsResponseSinglePt[4] = {"rec jet p_{T} (GeV/c)", "rec track p_{T} (GeV/c)", "gen jet p_{T} (GeV/c)", "gen track p_{T} (GeV/c)"};
3680   
3681   THnSparseD* hn4ResponseTrackPtJetPt =  new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
3682   
3683   for(Int_t i=0; i<4; i++){
3684     hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]); 
3685   }
3686   
3687
3688   // fill 
3689   
3690
3691 }
3692
3693 //_____________________________________________________________________________________________________________________________________
3694 void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
3695 {
3696
3697   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);
3698
3699 }
3700
3701 //_____________________________________________________________________________________________________________________________________
3702 void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
3703 {
3704
3705   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);
3706
3707 }
3708
3709 //_____________________________________________________________________________________________________________________________________
3710 void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
3711 {
3712
3713   ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);
3714
3715 }
3716
3717 //______________________________________________________________________________________________________________________________________________________
3718 void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
3719 {
3720   // read single track correction (pt) from file
3721   // type: efficiency / response / secondaries correction 
3722
3723   if(!((type == kFlagEfficiency)  || (type == kFlagResponse) || (type == kFlagSecondaries))){
3724     Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
3725     return;
3726   }
3727
3728   TFile f(strfile,"READ");
3729
3730   if(!f.IsOpen()){
3731     Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
3732     return;
3733   }
3734
3735   if(fDebug>0 && type==kFlagEfficiency)  Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3736   if(fDebug>0 && type==kFlagResponse)    Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3737   if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3738  
3739   if(strdir && strdir.Length()) gDirectory->cd(strdir);
3740
3741   TList* list = 0;
3742
3743   if(strlist && strlist.Length()){
3744     
3745     if(!(list = (TList*) gDirectory->Get(strlist))){ 
3746       Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3747       return;
3748     }
3749   }  
3750   
3751   TH1F*      h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
3752   THnSparse* hnCorrHist = 0;
3753   
3754   if(type == kFlagEfficiency || type == kFlagSecondaries){
3755
3756     if(list) h1CorrHist = (TH1F*) list->FindObject(strname); 
3757     else     h1CorrHist = (TH1F*) gDirectory->Get(strname);
3758
3759     if(!h1CorrHist){
3760       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3761       return;
3762     }
3763
3764   }
3765   else if(type == kFlagResponse){
3766     
3767     if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
3768     else     hnCorrHist = (THnSparse*) gDirectory->Get(strname);
3769
3770     if(!hnCorrHist){
3771       Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
3772       return;
3773     }
3774  
3775   }
3776
3777   if(h1CorrHist) h1CorrHist->SetDirectory(0);
3778   //if(hnCorrHist) hnCorrHist->SetDirectory(0);
3779   
3780   f.Close();  
3781
3782   if(type == kFlagEfficiency)       fh1EffSinglePt      = h1CorrHist; 
3783   else if(type == kFlagResponse)    fhnResponseSinglePt = hnCorrHist;
3784   else if(type == kFlagSecondaries) fh1SecCorrSinglePt  = h1CorrHist; 
3785  
3786 }
3787
3788 //________________________________________________________________________________________________________________
3789 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
3790
3791   // read track pt spec from task ouput - standard dir/list 
3792      
3793   TString strdir  = "PWG4_FragmentationFunction_" + strID;
3794   TString strlist = "fracfunc_" + strID;
3795     
3796   ReadRawPtSpec(strInfile,strdir,strlist);
3797 }
3798
3799 //_______________________________________________________________________________________________________
3800 void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
3801 {
3802   // get raw pt spectra from file 
3803
3804
3805   // book histos
3806   fNCorrectionLevelsSinglePt = 0; 
3807   fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
3808   AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum
3809
3810   // get raw pt spec from input file, normalize
3811
3812   TFile f(strfile,"READ");
3813
3814   if(!f.IsOpen()){
3815     Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
3816     return;
3817   }
3818
3819   if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
3820
3821   gDirectory->cd(strdir);
3822
3823   TList* list = 0;
3824   
3825   if(!(list = (TList*) gDirectory->Get(strlist))){ 
3826     Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
3827     return;
3828   }
3829
3830   TString hnameTrackPt("fh1TrackQAPtRecCuts");
3831   TString hnameEvtSel("fh1EvtSelection");
3832
3833   TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
3834   TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);  
3835
3836   if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
3837   if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
3838  
3839   
3840   //Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
3841
3842   // evts after physics selection
3843   Float_t nEvents = fh1EvtSel->GetEntries() 
3844     - fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(1))   // ev