// *************************************************************************
// *                                                                       *
// *               corrections to Fragmentation Functions                  *
// *                                                                       *
// *************************************************************************


/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id: */

#include "TMath.h"
#include "TH2F.h"
#include "THnSparse.h"
#include "TFile.h"
#include "TDirectory.h"
#include "AliCFUnfolding.h"
#include "AliFragmentationFunctionCorrections.h"
#include <iostream> // OB TEST!!!
#include <fstream>
using std::cout;
using std::endl;
using std::cerr;
using std::fstream;

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  correction class for ouput of AliAnalysisTaskFragmentationFunction       //     
//                                                                           //
//  corrections for: reconstruction efficiency, momentum resolution,         //
//                   secondaries, UE / HI background, fluctuations           // 
//                   back-correction on jet energy on dN/dxi                 //
//                                                                           //
//  read MC ouput and write out efficiency histos, response matrices etc.    //
//  read measured distributions and apply efficiency, response matrices etc. //
//                                                                           //
//  contact: o.busch@gsi.de                                                  //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////


ClassImp(AliFragmentationFunctionCorrections)

//________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections()
   : TObject()
   ,fDebug(0)
   ,fNJetPtSlices(0)
   ,fJetPtSlices(0)   
   ,fNJets(0)
   ,fNJetsBgr(0)
   ,fNHistoBinsSinglePt(0)
   ,fHistoBinsSinglePt(0)
   ,fNHistoBinsPt(0)
   ,fNHistoBinsZ(0)
   ,fNHistoBinsXi(0)
   ,fHistoBinsPt(0)
   ,fHistoBinsZ(0)
   ,fHistoBinsXi(0)
   ,fNCorrectionLevels(0)
   ,fCorrFF(0)
   ,fNCorrectionLevelsBgr(0)
   ,fCorrBgr(0)
   ,fNCorrectionLevelsSinglePt(0)
   ,fCorrSinglePt(0)
   ,fh1EffSinglePt(0)
   ,fh1EffPt(0) 
   ,fh1EffZ(0) 
   ,fh1EffXi(0) 
   ,fh1EffBgrPt(0)
   ,fh1EffBgrZ(0)
   ,fh1EffBgrXi(0)
   ,fh1BbBCorrSinglePt(0) 
   ,fh1BbBPt(0)
   ,fh1BbBZ(0)
   ,fh1BbBXi(0)
   ,fh1BbBBgrPt(0)
   ,fh1BbBBgrZ(0)
   ,fh1BbBBgrXi(0)
   ,fh1FoldingCorrPt(0)
   ,fh1FoldingCorrZ(0)
   ,fh1FoldingCorrXi(0)
   ,fh1SecCorrPt(0)
   ,fh1SecCorrZ(0)
   ,fh1SecCorrXi(0)
   ,fh1SecCorrBgrPt(0)
   ,fh1SecCorrBgrZ(0)
   ,fh1SecCorrBgrXi(0)
   ,fh1FFTrackPtBackFolded(0)   
   ,fh1FFZBackFolded(0)
   ,fh1FFXiBackFolded(0)        
   ,fh1FFRatioTrackPtFolded(0)
   ,fh1FFRatioZFolded(0)
   ,fh1FFRatioXiFolded(0)
   ,fh1FFRatioTrackPtBackFolded(0)
   ,fh1FFRatioZBackFolded(0)
   ,fh1FFRatioXiBackFolded(0)
   ,fh1SingleTrackPtBackFolded(0)   
   ,fh1RatioSingleTrackPtFolded(0) 
   ,fh1RatioSingleTrackPtBackFolded(0)
   ,fhnResponseSinglePt(0)
   ,fhnResponsePt(0)
   ,fhnResponseZ(0) 
   ,fhnResponseXi(0)       
   ,fh1FFTrackPtPrior(0)
   ,fh1FFZPrior(0)
   ,fh1FFXiPrior(0)
   ,fh1SecCorrSinglePt(0)
{
   // default constructor
}

//________________________________________________________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragmentationFunctionCorrections(const  AliFragmentationFunctionCorrections &copy)
  : TObject()
  ,fDebug(copy.fDebug)                        
  ,fNJetPtSlices(copy.fNJetPtSlices)                 
  ,fJetPtSlices(copy.fJetPtSlices)                  
  ,fNJets(copy.fNJets)                        
  ,fNJetsBgr(copy.fNJetsBgr)
  ,fNHistoBinsSinglePt(copy.fNHistoBinsSinglePt)  
  ,fHistoBinsSinglePt(copy.fHistoBinsSinglePt) 
  ,fNHistoBinsPt(copy.fNHistoBinsPt)                 
  ,fNHistoBinsZ(copy.fNHistoBinsZ)                  
  ,fNHistoBinsXi(copy.fNHistoBinsXi)                 
  ,fHistoBinsPt(copy.fHistoBinsPt)                  
  ,fHistoBinsZ(copy.fHistoBinsZ)                   
  ,fHistoBinsXi(copy.fHistoBinsXi)                  
  ,fNCorrectionLevels(copy.fNCorrectionLevels)            
  ,fCorrFF(copy.fCorrFF)                       
  ,fNCorrectionLevelsBgr(copy.fNCorrectionLevelsBgr)         
  ,fCorrBgr(copy.fCorrBgr)    
  ,fNCorrectionLevelsSinglePt(copy.fNCorrectionLevelsSinglePt)
  ,fCorrSinglePt(copy.fCorrSinglePt)   
  ,fh1EffSinglePt(copy.fh1EffSinglePt)
  ,fh1EffPt(copy.fh1EffPt)                      
  ,fh1EffZ(copy.fh1EffZ)                       
  ,fh1EffXi(copy.fh1EffXi)                      
  ,fh1EffBgrPt(copy.fh1EffBgrPt)                   
  ,fh1EffBgrZ(copy.fh1EffBgrZ)                    
  ,fh1EffBgrXi(copy.fh1EffBgrXi)
  ,fh1BbBCorrSinglePt(copy.fh1BbBCorrSinglePt) 
  ,fh1BbBPt(copy.fh1BbBPt)
  ,fh1BbBZ(copy.fh1BbBZ)
  ,fh1BbBXi(copy.fh1BbBXi)
  ,fh1BbBBgrPt(copy.fh1BbBBgrPt)
  ,fh1BbBBgrZ(copy.fh1BbBBgrZ)
  ,fh1BbBBgrXi(copy.fh1BbBBgrXi)
  ,fh1FoldingCorrPt(copy.fh1FoldingCorrPt)
  ,fh1FoldingCorrZ(copy.fh1FoldingCorrZ)
  ,fh1FoldingCorrXi(copy.fh1FoldingCorrXi)
  ,fh1SecCorrPt(copy.fh1SecCorrPt)
  ,fh1SecCorrZ(copy.fh1SecCorrZ)
  ,fh1SecCorrXi(copy.fh1SecCorrXi)
  ,fh1SecCorrBgrPt(copy.fh1SecCorrBgrPt)
  ,fh1SecCorrBgrZ(copy.fh1SecCorrBgrZ) 
  ,fh1SecCorrBgrXi(copy.fh1SecCorrBgrXi)
  ,fh1FFTrackPtBackFolded(copy.fh1FFTrackPtBackFolded)        
  ,fh1FFZBackFolded(copy.fh1FFZBackFolded)              
  ,fh1FFXiBackFolded(copy.fh1FFXiBackFolded)             
  ,fh1FFRatioTrackPtFolded(copy.fh1FFRatioTrackPtFolded)       
  ,fh1FFRatioZFolded(copy.fh1FFRatioZFolded)             
  ,fh1FFRatioXiFolded(copy.fh1FFRatioXiFolded)            
  ,fh1FFRatioTrackPtBackFolded(copy.fh1FFRatioTrackPtBackFolded)   
  ,fh1FFRatioZBackFolded(copy.fh1FFRatioZBackFolded)         
  ,fh1FFRatioXiBackFolded(copy.fh1FFRatioXiBackFolded)        
  ,fh1SingleTrackPtBackFolded(copy.fh1SingleTrackPtBackFolded)     
  ,fh1RatioSingleTrackPtFolded(copy.fh1RatioSingleTrackPtFolded)    
  ,fh1RatioSingleTrackPtBackFolded(copy.fh1RatioSingleTrackPtBackFolded)
  ,fhnResponseSinglePt(copy.fhnResponseSinglePt)                 
  ,fhnResponsePt(copy.fhnResponsePt)                 
  ,fhnResponseZ(copy.fhnResponseZ)                  
  ,fhnResponseXi(copy.fhnResponseXi)                 
  ,fh1FFTrackPtPrior(copy.fh1FFTrackPtPrior)   
  ,fh1FFZPrior(copy.fh1FFZPrior)
  ,fh1FFXiPrior(copy.fh1FFXiPrior)
  ,fh1SecCorrSinglePt(copy.fh1SecCorrSinglePt)
{
  // copy constructor
  
}

// ______________________________________________________________________________________________________________________________
AliFragmentationFunctionCorrections& AliFragmentationFunctionCorrections::operator=(const AliFragmentationFunctionCorrections& o)
{
  // assignment
  
  if(this!=&o){
    TObject::operator=(o);
    fDebug                          = o.fDebug;                        
    fNJetPtSlices                   = o.fNJetPtSlices;                 
    fJetPtSlices                    = o.fJetPtSlices;                  
    fNJets                          = o.fNJets;                        
    fNJetsBgr                       = o.fNJetsBgr;                     
    fNHistoBinsSinglePt             = o.fNHistoBinsSinglePt;  
    fHistoBinsSinglePt              = o.fHistoBinsSinglePt; 
    fNHistoBinsPt                   = o.fNHistoBinsPt;                 
    fNHistoBinsZ                    = o.fNHistoBinsZ;                  
    fNHistoBinsXi                   = o.fNHistoBinsXi;                 
    fHistoBinsPt                    = o.fHistoBinsPt;                  
    fHistoBinsZ                     = o.fHistoBinsZ;                   
    fHistoBinsXi                    = o.fHistoBinsXi;                  
    fNCorrectionLevels              = o.fNCorrectionLevels;            
    fCorrFF                         = o.fCorrFF;                       
    fNCorrectionLevelsBgr           = o.fNCorrectionLevelsBgr;         
    fCorrBgr                        = o.fCorrBgr;                      
    fNCorrectionLevelsSinglePt      = o.fNCorrectionLevelsSinglePt;
    fCorrSinglePt                   = o.fCorrSinglePt;
    fh1EffSinglePt                  = o.fh1EffSinglePt;
    fh1EffPt                        = o.fh1EffPt;                      
    fh1EffZ                         = o.fh1EffZ;                       
    fh1EffXi                        = o.fh1EffXi;                      
    fh1EffBgrPt                     = o.fh1EffBgrPt;                   
    fh1EffBgrZ                      = o.fh1EffBgrZ;                    
    fh1EffBgrXi                     = o.fh1EffBgrXi;                   
    fh1BbBCorrSinglePt              = o.fh1BbBCorrSinglePt; 
    fh1BbBPt                        = o.fh1BbBPt;
    fh1BbBZ                         = o.fh1BbBZ;
    fh1BbBXi                        = o.fh1BbBXi;
    fh1BbBBgrPt                     = o.fh1BbBBgrPt;
    fh1BbBBgrZ                      = o.fh1BbBBgrZ;
    fh1BbBBgrXi                     = o.fh1BbBBgrXi;
    fh1FoldingCorrPt                = o.fh1FoldingCorrPt;
    fh1FoldingCorrZ                 = o.fh1FoldingCorrZ;
    fh1FoldingCorrXi                = o.fh1FoldingCorrXi;
    fh1SecCorrPt                    = o.fh1SecCorrPt;
    fh1SecCorrZ                     = o.fh1SecCorrZ;
    fh1SecCorrXi                    = o.fh1SecCorrXi;
    fh1SecCorrBgrPt                 = o.fh1SecCorrBgrPt;
    fh1SecCorrBgrZ                  = o.fh1SecCorrBgrZ;
    fh1SecCorrBgrXi                 = o.fh1SecCorrBgrXi;
    fh1FFTrackPtBackFolded          = o.fh1FFTrackPtBackFolded;        
    fh1FFZBackFolded                = o.fh1FFZBackFolded;              
    fh1FFXiBackFolded               = o.fh1FFXiBackFolded;             
    fh1FFRatioTrackPtFolded         = o.fh1FFRatioTrackPtFolded;       
    fh1FFRatioZFolded               = o.fh1FFRatioZFolded;             
    fh1FFRatioXiFolded              = o.fh1FFRatioXiFolded;            
    fh1FFRatioTrackPtBackFolded     = o.fh1FFRatioTrackPtBackFolded;   
    fh1FFRatioZBackFolded           = o.fh1FFRatioZBackFolded;         
    fh1FFRatioXiBackFolded          = o.fh1FFRatioXiBackFolded;        
    fh1SingleTrackPtBackFolded      = o.fh1SingleTrackPtBackFolded;     
    fh1RatioSingleTrackPtFolded     = o.fh1RatioSingleTrackPtFolded;    
    fh1RatioSingleTrackPtBackFolded = o.fh1RatioSingleTrackPtBackFolded;
    fhnResponseSinglePt             = o.fhnResponseSinglePt;                 
    fhnResponsePt                   = o.fhnResponsePt;                 
    fhnResponseZ                    = o.fhnResponseZ;                  
    fhnResponseXi                   = o.fhnResponseXi;                 
    fh1FFTrackPtPrior               = o.fh1FFTrackPtPrior;
    fh1FFZPrior                     = o.fh1FFZPrior;
    fh1FFXiPrior                    = o.fh1FFXiPrior;
    fh1SecCorrSinglePt              = o.fh1SecCorrSinglePt;
  }
  
  return *this;
}

//_________________________________________________________________________
AliFragmentationFunctionCorrections::~AliFragmentationFunctionCorrections()
{
  // destructor  

  if(fJetPtSlices) delete fJetPtSlices;
  if(fNJets)       delete fNJets;
  if(fNJetsBgr)    delete fNJetsBgr;

  DeleteHistoArray(fh1EffPt);
  DeleteHistoArray(fh1EffXi);
  DeleteHistoArray(fh1EffZ );

  DeleteHistoArray(fh1EffBgrPt);
  DeleteHistoArray(fh1EffBgrXi);
  DeleteHistoArray(fh1EffBgrZ);

  // BbB
  DeleteHistoArray(fh1BbBPt);
  DeleteHistoArray(fh1BbBXi);
  DeleteHistoArray(fh1BbBZ);

  DeleteHistoArray(fh1BbBBgrPt);
  DeleteHistoArray(fh1BbBBgrXi);
  DeleteHistoArray(fh1BbBBgrZ);

  DeleteHistoArray(fh1FoldingCorrPt);
  DeleteHistoArray(fh1FoldingCorrXi);
  DeleteHistoArray(fh1FoldingCorrZ);

  // sec corr
  DeleteHistoArray(fh1SecCorrPt);
  DeleteHistoArray(fh1SecCorrXi);
  DeleteHistoArray(fh1SecCorrZ);

  DeleteHistoArray(fh1SecCorrBgrPt);
  DeleteHistoArray(fh1SecCorrBgrXi);
  DeleteHistoArray(fh1SecCorrBgrZ);

  // unfolding

  DeleteHistoArray(fh1FFTrackPtBackFolded);
  DeleteHistoArray(fh1FFZBackFolded);
  DeleteHistoArray(fh1FFXiBackFolded);

  DeleteHistoArray(fh1FFRatioTrackPtFolded);
  DeleteHistoArray(fh1FFRatioZFolded);
  DeleteHistoArray(fh1FFRatioXiFolded);

  DeleteHistoArray(fh1FFRatioTrackPtBackFolded);
  DeleteHistoArray(fh1FFRatioZBackFolded);
  DeleteHistoArray(fh1FFRatioXiBackFolded);

  DeleteTHnSparseArray(fhnResponsePt);
  DeleteTHnSparseArray(fhnResponseZ);
  DeleteTHnSparseArray(fhnResponseXi);

  DeleteHistoArray(fh1FFTrackPtPrior);
  DeleteHistoArray(fh1FFZPrior);
  DeleteHistoArray(fh1FFXiPrior);

  // clean up corrected FF 
  
  for(Int_t c=0; c<fNCorrectionLevels; c++) delete  fCorrFF[c];
  delete[] fCorrFF;

  // clean up bgr

  for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) delete fCorrBgr[c];
  delete[] fCorrBgr;

  // clean up inclusive pt 
  for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) delete fCorrSinglePt[c];
  delete[] fCorrSinglePt;

  delete[] fNHistoBinsPt;
  delete[] fNHistoBinsZ;
  delete[] fNHistoBinsXi;

  // clean up histo bins

  if(fHistoBinsSinglePt) delete fHistoBinsSinglePt; 

  for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsPt[i];
  for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsZ[i];
  for(Int_t i=0; i<fNJetPtSlices; i++) delete fHistoBinsXi[i];

  delete[] fHistoBinsPt;
  delete[] fHistoBinsZ;
  delete[] fHistoBinsXi;
}

//_________________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos()
  : TObject()
  ,fArraySize(0)
  ,fh1CorrFFTrackPt(0)
  ,fh1CorrFFZ(0)
  ,fh1CorrFFXi(0)
  ,fCorrLabel(0)
{
  // default constructor
  
}

//__________________________________________________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& copy)
  : TObject()
  ,fArraySize(copy.fArraySize)
  ,fh1CorrFFTrackPt(0)
  ,fh1CorrFFZ(0)
  ,fh1CorrFFXi(0)
  ,fCorrLabel(copy.fCorrLabel)
{
  // copy constructor

  fh1CorrFFTrackPt = new TH1F*[copy.fArraySize];
  fh1CorrFFZ       = new TH1F*[copy.fArraySize];
  fh1CorrFFXi      = new TH1F*[copy.fArraySize];

  for(Int_t i=0; i<copy.fArraySize; i++){
    fh1CorrFFTrackPt[i] = new TH1F(*(copy.fh1CorrFFTrackPt[i]));
    fh1CorrFFZ[i]       = new TH1F(*(copy.fh1CorrFFZ[i]));
    fh1CorrFFXi[i]      = new TH1F(*(copy.fh1CorrFFXi[i]));
  }
}

//_______________________________________________________________________________________________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::operator=(const AliFragmentationFunctionCorrections::AliFragFuncCorrHistos& o)
{
  // assignment
  
  if(this!=&o){
    TObject::operator=(o);
    Int_t fArraySize_new = o.fArraySize;
    fCorrLabel       = o.fCorrLabel;
 
    TH1F** fh1CorrFFTrackPt_new = new TH1F*[fArraySize_new];
    TH1F** fh1CorrFFZ_new       = new TH1F*[fArraySize_new];
    TH1F** fh1CorrFFXi_new      = new TH1F*[fArraySize_new];

    for(Int_t i=0; i<fArraySize_new; i++){
      fh1CorrFFTrackPt_new[i] = new TH1F(*(o.fh1CorrFFTrackPt[i]));
      fh1CorrFFZ_new[i]       = new TH1F(*(o.fh1CorrFFZ[i]));
      fh1CorrFFXi_new[i]      = new TH1F(*(o.fh1CorrFFXi[i]));
    }
    
    // ---

    if(fArraySize){
      for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
      for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
      for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];          
    }

    if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
    if(fh1CorrFFZ)       delete[] fh1CorrFFZ;
    if(fh1CorrFFXi)      delete[] fh1CorrFFXi;
        
    // ---

    fArraySize = fArraySize_new;

    fh1CorrFFTrackPt = fh1CorrFFTrackPt_new;
    fh1CorrFFZ       = fh1CorrFFZ_new;
    fh1CorrFFXi      = fh1CorrFFXi_new;
    
    for(Int_t i=0; i<fArraySize; i++){
      fh1CorrFFTrackPt[i] = fh1CorrFFTrackPt_new[i];
      fh1CorrFFZ[i]       = fh1CorrFFZ_new[i];
      fh1CorrFFXi[i]      = fh1CorrFFXi_new[i];
    }
  }
  
  return *this;
}

//__________________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::~AliFragFuncCorrHistos()
{
  // destructor 

  if(fArraySize){
    for(int i=0; i<fArraySize; i++) delete fh1CorrFFTrackPt[i];
    for(int i=0; i<fArraySize; i++) delete fh1CorrFFZ[i];
    for(int i=0; i<fArraySize; i++) delete fh1CorrFFXi[i];          
  }

  if(fh1CorrFFTrackPt) delete[] fh1CorrFFTrackPt;
  if(fh1CorrFFZ)       delete[] fh1CorrFFZ;
  if(fh1CorrFFXi)      delete[] fh1CorrFFXi;

}

//___________________________________________________________________________________________________________________
AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AliFragFuncCorrHistos(const char* label, Int_t arraySize)
  : TObject()
  ,fArraySize(0)
  ,fh1CorrFFTrackPt(0)
  ,fh1CorrFFZ(0)
  ,fh1CorrFFXi(0)
  ,fCorrLabel(label)
{
  // constructor

  fArraySize = arraySize;
  fh1CorrFFTrackPt = new TH1F*[arraySize];
  fh1CorrFFZ       = new TH1F*[arraySize];
  fh1CorrFFXi      = new TH1F*[arraySize];

  for(int i=0; i<arraySize; i++) fh1CorrFFTrackPt[i] = new TH1F;
  for(int i=0; i<arraySize; i++) fh1CorrFFZ[i]       = new TH1F;
  for(int i=0; i<arraySize; i++) fh1CorrFFXi[i]      = new TH1F;
}

//_______________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::AddCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
{
  // place histo in array, add corrLabel to histo name

  if(slice>=fArraySize){
    Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
    return;
  }

  if(histPt){
    TString name = histPt->GetName();
    if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
    histPt->SetName(name);
    
    if(!(histPt->GetSumw2())) histPt->Sumw2(); 

    new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
  }
  
  if(histZ){
    TString name = histZ->GetName();
    if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
    histZ->SetName(name);

    if(!(histZ->GetSumw2())) histZ->Sumw2(); 

    new(fh1CorrFFZ[slice]) TH1F(*histZ);
  }

  if(histXi){
    TString name = histXi->GetName();
    if(fCorrLabel.Length()>0) name += "_"+fCorrLabel;
    histXi->SetName(name);

    if(!(histXi->GetSumw2())) histXi->Sumw2(); 

    new(fh1CorrFFXi[slice]) TH1F(*histXi);
  }
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::ReplaceCorrHistos(Int_t slice,TH1F* histPt,TH1F* histZ,TH1F* histXi)
{
  // replace histo in array 

  if(slice>=fArraySize){
    Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
    return;
  }

  if(histPt){
    if(!(histPt->GetSumw2())) histPt->Sumw2(); 

    TString name = histPt->GetName();
    histPt->SetName(name);

    delete fh1CorrFFTrackPt[slice];
    fh1CorrFFTrackPt[slice] =  new TH1F;
    new(fh1CorrFFTrackPt[slice]) TH1F(*histPt);
  }
  
  if(histZ){
    if(!(histZ->GetSumw2())) histZ->Sumw2(); 

    TString name = histZ->GetName();
    histZ->SetName(name);

    delete fh1CorrFFZ[slice];
    fh1CorrFFZ[slice] = new TH1F;
    new(fh1CorrFFZ[slice]) TH1F(*histZ);
  }

  if(histXi){
    if(!(histXi->GetSumw2())) histXi->Sumw2(); 
    
    TString name = histXi->GetName();
    histXi->SetName(name);

    delete fh1CorrFFXi[slice];
    fh1CorrFFXi[slice] = new TH1F;
    new(fh1CorrFFXi[slice]) TH1F(*histXi);
  }
}

// ___________________________________________________________________________________________
TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetTrackPt(const Int_t slice)
{ 
  // return pt histo 
  
  if(slice>=fArraySize){
    Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
    return 0;
  }

  return fh1CorrFFTrackPt[slice]; 

}

// ______________________________________________________________________________________
TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetZ(const Int_t slice)
{ 
  // return z histo 
  
  if(slice>=fArraySize){
    Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
    return 0;
  }

  return fh1CorrFFZ[slice]; 
}

// ________________________________________________________________________________________
TH1F* AliFragmentationFunctionCorrections::AliFragFuncCorrHistos::GetXi(const Int_t slice)
{ 
  // return xi histo

  if(slice>=fArraySize){
    Printf("%s:%d -- slice > array size", (char*)__FILE__,__LINE__);
    return 0;
  }

  return fh1CorrFFXi[slice]; 
}

// __________________________________________________________________________
void AliFragmentationFunctionCorrections::DeleteHistoArray(TH1F** hist) const
{
  // delete array of TH1 
 
  if(hist){
    for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
    delete[] hist;
  }
}

// ____________________________________________________________________________________
void AliFragmentationFunctionCorrections::DeleteTHnSparseArray(THnSparse** hist) const
{
  // delete array of THnSparse 
 
  if(hist){
    for(Int_t i=0; i<fNJetPtSlices; i++) delete hist[i];
    delete[] hist;
  }
}

// _________________________________________________________
TH1F** AliFragmentationFunctionCorrections::BookHistoArray()
{
  // book array of TH1

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return 0;
  }

  TH1F** hist = new TH1F*[fNJetPtSlices];
  for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new TH1F();
  
  return hist;
}

//__________________________________________________________________
THnSparse** AliFragmentationFunctionCorrections::BookTHnSparseArray()
{
  // book array of THnSparse

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return 0;
  }

  THnSparse** hist = new THnSparse*[fNJetPtSlices];
  for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = new THnSparseF(); 
  
  return hist;
}

//_____________________________________________________________________________
void AliFragmentationFunctionCorrections::AddCorrectionLevel(const char* label)
{
  // increase corr level 

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return;
  }

  if(fNCorrectionLevels >= fgMaxNCorrectionLevels){
    Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
    return;
  }

  fCorrFF[fNCorrectionLevels] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
  fNCorrectionLevels++;
}


//________________________________________________________________________________
void AliFragmentationFunctionCorrections::AddCorrectionLevelBgr(const char* label)
{
  // increase corr level for bgr FF

  if(!fNJetPtSlices){
    if(fDebug>0)  Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return;
  }
  
  if(fNCorrectionLevelsBgr >= fgMaxNCorrectionLevels){
    Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
    return;
  }

  fCorrBgr[fNCorrectionLevelsBgr] = new AliFragFuncCorrHistos(label,fNJetPtSlices);
  fNCorrectionLevelsBgr++;
}

//_____________________________________________________________________________________
void AliFragmentationFunctionCorrections::AddCorrectionLevelSinglePt(const char* label)
{
  // increase corr level single pt spec

  Int_t nJetPtSlicesSingle = 1; // pro forma 

  if(fNCorrectionLevelsSinglePt >= fgMaxNCorrectionLevels){
    Printf("%s:%d -- max correction level exceeded", (char*)__FILE__,__LINE__);
    return;
  }

  fCorrSinglePt[fNCorrectionLevelsSinglePt] = new AliFragFuncCorrHistos(label,nJetPtSlicesSingle);
  fNCorrectionLevelsSinglePt++;
}

//_____________________________________________________________________________________________
void AliFragmentationFunctionCorrections::SetJetPtSlices(Float_t* bins, const Int_t nJetPtSlices)
{ 
  // set jet pt slices
  // once slices are known, can book all other histos 

  fNJetPtSlices = nJetPtSlices;

  const Int_t size = nJetPtSlices+1;
  fJetPtSlices  = new TArrayF(size,bins);
 
  // nJets array

  fNJets = new TArrayF(size);
  fNJets->Reset(0);

  fNJetsBgr = new TArrayF(size);
  fNJetsBgr->Reset(0);

  // histos 

  fNCorrectionLevels = 0; 
  fCorrFF = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
  AddCorrectionLevel(); // first 'correction' level = raw FF

  fNCorrectionLevelsBgr = 0; 
  fCorrBgr = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
  AddCorrectionLevelBgr(); // first 'correction' level = raw bgr dist

  // eff histos

  fh1EffPt = BookHistoArray();
  fh1EffXi = BookHistoArray();
  fh1EffZ  = BookHistoArray();

  fh1EffBgrPt = BookHistoArray();
  fh1EffBgrXi = BookHistoArray();
  fh1EffBgrZ  = BookHistoArray();


  // BbB

  fh1BbBPt = BookHistoArray();
  fh1BbBXi = BookHistoArray();
  fh1BbBZ  = BookHistoArray();

  fh1BbBBgrPt = BookHistoArray();
  fh1BbBBgrXi = BookHistoArray();
  fh1BbBBgrZ  = BookHistoArray();

  fh1FoldingCorrPt = BookHistoArray();
  fh1FoldingCorrXi = BookHistoArray();
  fh1FoldingCorrZ  = BookHistoArray();

  // SecCorr

  fh1SecCorrPt = BookHistoArray();
  fh1SecCorrXi = BookHistoArray();
  fh1SecCorrZ  = BookHistoArray();

  fh1SecCorrBgrPt = BookHistoArray();
  fh1SecCorrBgrXi = BookHistoArray();
  fh1SecCorrBgrZ  = BookHistoArray();

  // unfolding

  fh1FFTrackPtBackFolded = BookHistoArray();
  fh1FFXiBackFolded      = BookHistoArray();
  fh1FFZBackFolded       = BookHistoArray();

  fh1FFRatioTrackPtFolded = BookHistoArray();
  fh1FFRatioXiFolded      = BookHistoArray();
  fh1FFRatioZFolded       = BookHistoArray();
  
  fh1FFRatioTrackPtBackFolded = BookHistoArray();
  fh1FFRatioXiBackFolded      = BookHistoArray();
  fh1FFRatioZBackFolded       = BookHistoArray();

  //

  fhnResponsePt      = BookTHnSparseArray(); 
  fhnResponseZ       = BookTHnSparseArray();
  fhnResponseXi      = BookTHnSparseArray();

  fh1FFTrackPtPrior =  BookHistoArray();
  fh1FFZPrior       =  BookHistoArray();
  fh1FFXiPrior      =  BookHistoArray();

  // histos bins 

  fNHistoBinsPt = new Int_t[fNJetPtSlices];
  fNHistoBinsXi = new Int_t[fNJetPtSlices];
  fNHistoBinsZ  = new Int_t[fNJetPtSlices];

  for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsPt[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsXi[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) fNHistoBinsZ[i]  = 0;
  
  fHistoBinsPt = new TArrayD*[fNJetPtSlices];
  fHistoBinsXi = new TArrayD*[fNJetPtSlices];
  fHistoBinsZ  = new TArrayD*[fNJetPtSlices];

  for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsPt[i] = new TArrayD(0);
  for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsXi[i] = new TArrayD(0);
  for(Int_t i=0; i<fNJetPtSlices; i++) fHistoBinsZ[i]  = new TArrayD(0);
}

//_____________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t sizeBins, Double_t* bins, const Int_t type)
{ 
  // set histo bins for jet pt slice
  // if binning undefined for any slice, original binning will be used

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return;
  }
  
  if(jetPtSlice>=fNJetPtSlices){
    Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
    return;
  }
  
  if(type == kFlagPt){
    fNHistoBinsPt[jetPtSlice] = sizeBins-1;
    fHistoBinsPt[jetPtSlice]->Set(sizeBins,bins); 
  }
  else if(type==kFlagZ){
    fNHistoBinsZ[jetPtSlice] = sizeBins-1;
    fHistoBinsZ[jetPtSlice]->Set(sizeBins,bins); 
  }

  else if(type==kFlagXi){
    fNHistoBinsXi[jetPtSlice] = sizeBins-1;
    fHistoBinsXi[jetPtSlice]->Set(sizeBins,bins); 
  }
}

//__________________________________________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::SetHistoBins(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
{ 
  // set histo bins for jet pt slice
  // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
  // array size of binsLimits: nBinsLimits 
  // array size of binsWidth: nBinsLimits-1 
  // binsLimits have to be in increasing order
  // if binning undefined for any slice, original binning will be used

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return;
  }

  if(jetPtSlice>=fNJetPtSlices){
    Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
    return;
  }


  Double_t binLimitMin = binsLimits[0];
  Double_t binLimitMax = binsLimits[nBinsLimits-1];

  Double_t binLimit = binLimitMin; // start value 
  
  Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
  TArrayD binsArray(sizeUpperLim);
  Int_t nBins = 0; 
  binsArray.SetAt(binLimitMin,nBins++);

  while(binLimit<0.999999*binLimitMax && nBins<sizeUpperLim){ // 0.999 numerical safety factor

    Int_t currentSlice = -1;
    for(Int_t i=0; i<nBinsLimits; i++){
      if(binLimit >= 0.999999*binsLimits[i]) currentSlice = i;  
    }
    
    Double_t currentBinWidth = binsWidth[currentSlice];
    binLimit += currentBinWidth;

    // if(type == kFlagZ) std::cout<<" SetHistoBins: nBins "<<nBins<<" binLimit "<<binLimit<<" binLimitMax "<<binLimitMax
    // 				<<" currentSlice "<<currentSlice<<std::endl;

    binsArray.SetAt(binLimit,nBins++);
  }
  
  Double_t* bins = binsArray.GetArray();

  SetHistoBins(jetPtSlice,nBins,bins,type); 
}

//_____________________________________________________________________________________________________________________________________
TArrayD* AliFragmentationFunctionCorrections::GetHistoBins(const Int_t& jetPtSlice, const Int_t& type)
{ 
  // set histo bins for jet pt slice
  // if binning undefined for any slice, original binning will be used

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return 0;
  }
  
  if(jetPtSlice>=fNJetPtSlices){
    Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
    return 0;
  }
  
  if(type == kFlagPt){
    return fHistoBinsPt[jetPtSlice]; 
  }
  else if(type==kFlagZ){
    return fHistoBinsZ[jetPtSlice]; 
  }
  
  else if(type==kFlagXi){
    return fHistoBinsXi[jetPtSlice]; 
  }
  else{ 
    Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
    return 0;
  }
}

//__________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t sizeBins, Double_t* bins)
{ 
  // set histo bins for inclusive pt spectra
  // if binning undefined, original binning will be used

  fNHistoBinsSinglePt = sizeBins-1;
  
  fHistoBinsSinglePt = new TArrayD(sizeBins);
  fHistoBinsSinglePt->Set(sizeBins,bins);  
}

//__________________________________________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::SetHistoBinsSinglePt(const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth)
{ 
  // set histo bins for inclusive pt spectra 
  // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
  // array size of binsLimits: nBinsLimits 
  // array size of binsWidth: nBinsLimits-1 
  // binsLimits have to be in increasing order
  // if binning undefined for any slice, original binning will be used


  Double_t binLimitMin = binsLimits[0];
  Double_t binLimitMax = binsLimits[nBinsLimits-1];

  Double_t binLimit = binLimitMin; // start value 
  
  Int_t sizeUpperLim = 10000; //static_cast<Int_t>(binLimitMax/binsWidth[0])+1;
  TArrayD binsArray(sizeUpperLim);
  Int_t nBins = 0; 
  binsArray.SetAt(binLimitMin,nBins++);

  while(binLimit<binLimitMax && nBins<sizeUpperLim){

    Int_t currentSlice = -1;
    for(Int_t i=0; i<nBinsLimits; i++){
      if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
    }
    
    Double_t currentBinWidth = binsWidth[currentSlice];
    binLimit += currentBinWidth;

    binsArray.SetAt(binLimit,nBins++);
  }
  
  Double_t* bins = binsArray.GetArray();
  
  SetHistoBinsSinglePt(nBins,bins); 
}

//____________________________________________________________________________________
void AliFragmentationFunctionCorrections::NormalizeTH1(TH1* hist, const Float_t nJets)
{
  // FF normalization: divide by bin width and normalize to entries/jets
  // should also work for TH2, but to be tested !!!

  if(nJets == 0){ // nothing to do
    if(fDebug>0)  Printf("%s:%d -- normalize: nJets = 0, do nothing", (char*)__FILE__,__LINE__);
    return; 
  }
  
  Int_t nBins = hist->GetNbinsX()*hist->GetNbinsY()*hist->GetNbinsZ();

  for(Int_t bin=0; bin<=nBins; bin++){ // include under-/overflow (?)

    Double_t binWidth = hist->GetBinWidth(bin);
    Double_t binCont  = hist->GetBinContent(bin);
    Double_t binErr   = hist->GetBinError(bin);
    
    binCont /= binWidth;
    binErr  /= binWidth;

    hist->SetBinContent(bin,binCont);
    hist->SetBinError(bin,binErr);
  }

  hist->Scale(1/nJets);
}

//_____________________________________________________
void AliFragmentationFunctionCorrections::NormalizeFF()
{
  // normalize FF

  if(fNCorrectionLevels>1){
    Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
  }

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    if(fDebug>0) Printf(" normalizeFF: i %d, nJets %f",i,fNJets->At(i));

    NormalizeTH1(fCorrFF[0]->GetTrackPt(i),fNJets->At(i)); // always normalize corr level 0 = raw FF
    NormalizeTH1(fCorrFF[0]->GetZ(i),fNJets->At(i));
    NormalizeTH1(fCorrFF[0]->GetXi(i),fNJets->At(i));
  } 
}

//______________________________________________________
void AliFragmentationFunctionCorrections::NormalizeBgr()
{
  // normalize bgr/UE FF

  if(fNCorrectionLevelsBgr>1){
    Printf("%s:%d -- FF normalization should be done BEFORE any correction !!!", (char*)__FILE__,__LINE__);
  }

  for(Int_t i=0; i<fNJetPtSlices; i++){
    NormalizeTH1(fCorrBgr[0]->GetTrackPt(i), fNJetsBgr->At(i)); // always normalize corr level 0 = raw FF
    NormalizeTH1(fCorrBgr[0]->GetZ(i), fNJetsBgr->At(i));
    NormalizeTH1(fCorrBgr[0]->GetXi(i),fNJetsBgr->At(i));
  } 

}

//__________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strID, TString strFFID)
{ 
  // read raw FF - standard dir/list name
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  ReadRawFF(strfile,strdir,strlist,strFFID);
}

//____________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawFF(TString strfile, TString strdir, TString strlist, TString strFFID)
{
  // get raw FF from input file, project in jet pt slice
  // normalization done separately 

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read FF from file %s, dir %s ",(char*)__FILE__,__LINE__,strfile.Data(), strdir.Data());

  gDirectory->cd(strdir);

  TList* list = 0;
  
  if(strlist && strlist.Length()){
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }

  TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data()));
  TString hnameTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
  TString hnameZ(Form("fh2FFZ%s",strFFID.Data()));
  TString hnameXi(Form("fh2FFXi%s",strFFID.Data()));

  TH1F* fh1FFJetPt   = 0;
  TH2F* fh2FFTrackPt = 0;
  TH2F* fh2FFZ       = 0;
  TH2F* fh2FFXi      = 0;
  
  if(list){
    fh1FFJetPt   = (TH1F*) list->FindObject(hnameJetPt);
    fh2FFTrackPt = (TH2F*) list->FindObject(hnameTrackPt);
    fh2FFZ       = (TH2F*) list->FindObject(hnameZ);  
    fh2FFXi      = (TH2F*) list->FindObject(hnameXi); 
  }
  else{
    fh1FFJetPt   = (TH1F*) gDirectory->Get(hnameJetPt);
    fh2FFTrackPt = (TH2F*) gDirectory->Get(hnameTrackPt);
    fh2FFZ       = (TH2F*) gDirectory->Get(hnameZ);  
    fh2FFXi      = (TH2F*) gDirectory->Get(hnameXi);
  }


  if(!fh1FFJetPt)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());   return; }
  if(!fh2FFTrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
  if(!fh2FFZ)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameZ.Data());       return; }
  if(!fh2FFXi)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameXi.Data());      return; }

  fh1FFJetPt->SetDirectory(0);
  fh2FFTrackPt->SetDirectory(0);
  fh2FFZ->SetDirectory(0);  
  fh2FFXi->SetDirectory(0); 

  f.Close();  


  // nJets per bin

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);
    
    Int_t binLo = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh1FFJetPt->FindBin(jetPtUpLim)) - 1;
				     
    Float_t nJetsBin = fh1FFJetPt->Integral(binLo,binUp);
    
    fNJets->SetAt(nJetsBin,i); 

    if(fDebug>0) Printf("jet pt %d to %d: nJets %f",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),fNJets->At(i));
  }
  
  // projections: FF 
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);

    Int_t binLo = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh2FFTrackPt->GetXaxis()->FindBin(jetPtUpLim))-1;

    if(binUp > fh2FFTrackPt->GetNbinsX()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameFFPt(Form("fh1FFTrackPt%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZ(Form("fh1FFZ%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXi(Form("fh1FFXi%s_%02d_%02d",strFFID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    // appendix 'unbinned' to avoid histos with same name after rebinning
    TH1F* projPt = (TH1F*) fh2FFTrackPt->ProjectionY(strNameFFPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    TH1F* projZ  = (TH1F*) fh2FFZ->ProjectionY(strNameFFZ+"_unBinned",binLo,binUp,"o");
    TH1F* projXi = (TH1F*) fh2FFXi->ProjectionY(strNameFFXi+"_unBinned",binLo,binUp,"o");
    
    if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameFFPt,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameFFZ,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameFFXi,fHistoBinsXi[i]->GetArray());

    projPt->SetNameTitle(strNameFFPt,"");
    projZ->SetNameTitle(strNameFFZ,"");
    projXi->SetNameTitle(strNameFFXi,"");

    // raw FF = corr level 0
    fCorrFF[0]->AddCorrHistos(i,projPt,projZ,projXi);
  }


  delete fh1FFJetPt;
  delete fh2FFTrackPt;
  delete fh2FFZ;  
  delete fh2FFXi; 
}

//_____________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strID, TString strBgrID, TString strFFID)
{ 
  // read raw FF - standard dir/list name
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  ReadRawBgr(strfile,strdir,strlist,strBgrID,strFFID);
}

//_______________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawBgr(TString strfile, TString strdir, TString strlist, TString strBgrID, TString strFFID)
{
  // get raw FF from input file, project in jet pt slice
  // use jet dN/dpt corresponding to strFFID, bgr FF to strBgrID+strFFID
  // e.g. "fh1FFJetPtRecCuts", "fh2FFXiBgrPerpRecCuts"
  // normalization done separately 

  TString strID = strBgrID + strFFID;
  
  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s, dir %s ",(char*)__FILE__,__LINE__,strBgrID.Data(),strfile.Data(),strdir.Data());

  gDirectory->cd(strdir);

  TList* list = 0;
  
  if(strlist && strlist.Length()){
    if(!(list = (TList*) gDirectory->Get(strlist))){
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }

  TString hnameNJets = "fh1nRecJetsCuts"; 
  TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); // not: strID.Data() !!! would not be proper normalization
  TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strID.Data()));
  TString hnameBgrZ(Form("fh2FFZ%s",strID.Data()));
  TString hnameBgrXi(Form("fh2FFXi%s",strID.Data()));

  TH1F* fh1NJets;
  TH1F* fh1FFJetPtBgr;  
  TH2F* fh2FFTrackPtBgr;
  TH2F* fh2FFZBgr;
  TH2F* fh2FFXiBgr;     
  
  if(list){
    fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
    fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
    fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
    fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);  
    fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi); 
  }
  else{
    fh1NJets        = (TH1F*) gDirectory->Get(hnameNJets); // needed for normalization of bgr out of 2 jets
    fh1FFJetPtBgr   = (TH1F*) gDirectory->Get(hnameJetPt);
    fh2FFTrackPtBgr = (TH2F*) gDirectory->Get(hnameBgrTrackPt);
    fh2FFZBgr       = (TH2F*) gDirectory->Get(hnameBgrZ);  
    fh2FFXiBgr      = (TH2F*) gDirectory->Get(hnameBgrXi);  
  }


  if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
  if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());              }
  if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
  if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
  if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }

  fh1FFJetPtBgr->SetDirectory(0);
  fh2FFTrackPtBgr->SetDirectory(0);
  fh2FFZBgr->SetDirectory(0);  
  fh2FFXiBgr->SetDirectory(0); 
  if(fh1NJets) fh1NJets->SetDirectory(0);

  f.Close();  

  // nJets per bin

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);
    
    Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
				     
    Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
    Double_t scaleF = 1;

    //if(strBgrID.Contains("Out2Jets")){  // scale by ratio 2 jets events / all events
    //  scaleF = fh1NJets->Integral(fh1NJets->FindBin(2),fh1NJets->GetNbinsX()) 
    //  / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());}


    if(strBgrID.Contains("OutAllJets") ){  // scale by ratio >3 jets events / all events

      if(fh1NJets){
	scaleF = fh1NJets->Integral(fh1NJets->FindBin(4),fh1NJets->GetNbinsX()) 
	  / fh1NJets->Integral(fh1NJets->FindBin(1),fh1NJets->GetNbinsX());
      }
      else{
	Printf("%s:%d -- use bgr OutAllJets, but histo fhn1NJets not found",(char*)__FILE__,__LINE__);    
      }
    }
    
 
    fNJetsBgr->SetAt(nJetsBin*scaleF,i); 

    if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
			static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);

  }
  
  // projections: FF 
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);

    Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;

    if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    // appendix 'unbinned' to avoid histos with same name after rebinning
    TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
    TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
    
    if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());

    projPt->SetNameTitle(strNameBgrPt,"");
    projZ->SetNameTitle(strNameBgrZ,"");
    projXi->SetNameTitle(strNameBgrXi,"");
    
    // raw bgr = corr level 0
    fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
  }

  delete fh1FFJetPtBgr;
  if(fh1NJets) delete fh1NJets;
  delete fh2FFTrackPtBgr;
  delete fh2FFZBgr;  
  delete fh2FFXiBgr; 
}

//_____________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strID, TString strFFID)
{ 
  // read raw FF - standard dir/list name
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  ReadRawBgrEmbedding(strfile,strdir,strlist,strFFID);
}

//_______________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawBgrEmbedding(TString strfile, TString strdir, TString strlist, TString strFFID)
{
  // get raw FF from input file, project in jet pt slice
  // for embedding, the bgr FF are taken from histos "fh1FFJetPtRecCuts", "fh2FFXiRecCuts"
  // normalization done separately 

  TString strBgrID = "BckgEmbed";
  TString strID = strBgrID + strFFID;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read Bgr %s from file %s ",(char*)__FILE__,__LINE__,strFFID.Data(),strfile.Data());

  gDirectory->cd(strdir);

  TList* list = 0;
  
  if(!(list = (TList*) gDirectory->Get(strlist))){
    Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
    return;
  }

  TString hnameNJets = "fh1nRecJetsCuts"; 
  TString hnameJetPt(Form("fh1FFJetPt%s",strFFID.Data())); 
  TString hnameBgrTrackPt(Form("fh2FFTrackPt%s",strFFID.Data()));
  TString hnameBgrZ(Form("fh2FFZ%s",strFFID.Data()));
  TString hnameBgrXi(Form("fh2FFXi%s",strFFID.Data()));

  TH1F* fh1NJets        = (TH1F*) list->FindObject(hnameNJets); // needed for normalization of bgr out of 2 jets
  TH1F* fh1FFJetPtBgr   = (TH1F*) list->FindObject(hnameJetPt);
  TH2F* fh2FFTrackPtBgr = (TH2F*) list->FindObject(hnameBgrTrackPt);
  TH2F* fh2FFZBgr       = (TH2F*) list->FindObject(hnameBgrZ);  
  TH2F* fh2FFXiBgr      = (TH2F*) list->FindObject(hnameBgrXi); 

  if(!fh1FFJetPtBgr)  { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameJetPt.Data());      return; }
  if(!fh1NJets)       { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameNJets.Data());      return; }
  if(!fh2FFTrackPtBgr){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrTrackPt.Data()); return; }
  if(!fh2FFZBgr)      { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrZ.Data());       return; }
  if(!fh2FFXiBgr)     { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameBgrXi.Data());      return; }

  fh1FFJetPtBgr->SetDirectory(0);
  fh1NJets->SetDirectory(0);
  fh2FFTrackPtBgr->SetDirectory(0);
  fh2FFZBgr->SetDirectory(0);  
  fh2FFXiBgr->SetDirectory(0); 

  f.Close();  

  // nJets per bin

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);
    
    Int_t binLo = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh1FFJetPtBgr->FindBin(jetPtUpLim)) - 1;
				     
    Float_t nJetsBin = fh1FFJetPtBgr->Integral(binLo,binUp);
    Double_t scaleF = 1;

    fNJetsBgr->SetAt(nJetsBin*scaleF,i); 

    if(fDebug>0) Printf("bgr jet pt %d to %d: nJets %f, scaleF %.2f",
			static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),nJetsBin,scaleF);

  }
  
  // projections: FF 
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);

    Int_t binLo = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh2FFTrackPtBgr->GetXaxis()->FindBin(jetPtUpLim))-1;

    if(binUp > fh2FFTrackPtBgr->GetNbinsX()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameBgrPt(Form("fh1BgrTrackPt%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameBgrZ(Form("fh1BgrZ%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameBgrXi(Form("fh1BgrXi%s_%02d_%02d",strID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    // appendix 'unbinned' to avoid histos with same name after rebinning
    TH1F* projPt = (TH1F*) fh2FFTrackPtBgr->ProjectionY(strNameBgrPt+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    TH1F* projZ  = (TH1F*) fh2FFZBgr->ProjectionY(strNameBgrZ+"_unBinned",binLo,binUp,"o");
    TH1F* projXi = (TH1F*) fh2FFXiBgr->ProjectionY(strNameBgrXi+"_unBinned",binLo,binUp,"o");
    
    if(fNHistoBinsPt[i]) projPt = (TH1F*) projPt->Rebin(fNHistoBinsPt[i],strNameBgrPt,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  projZ  = (TH1F*) projZ->Rebin(fNHistoBinsZ[i],strNameBgrZ,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) projXi = (TH1F*) projXi->Rebin(fNHistoBinsXi[i],strNameBgrXi,fHistoBinsXi[i]->GetArray());

    projPt->SetNameTitle(strNameBgrPt,"");
    projZ->SetNameTitle(strNameBgrZ,"");
    projXi->SetNameTitle(strNameBgrXi,"");
    
    // raw bgr = corr level 0
    fCorrBgr[0]->AddCorrHistos(i,projPt,projZ,projXi);
  }

  delete fh1FFJetPtBgr;
  delete fh1NJets;
  delete fh2FFTrackPtBgr;
  delete fh2FFZBgr;  
  delete fh2FFXiBgr; 
}


//__________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteOutput(TString strfile, TString strdir, Bool_t updateOutfile)
{
  // write histos to file
  // skip histos with 0 entries

  TString outfileOption = "RECREATE";
  if(updateOutfile) outfileOption = "UPDATE";

  TFile f(strfile,outfileOption); 

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write FF to file %s ",(char*)__FILE__,__LINE__,strfile.Data());

  if(strdir && strdir.Length()){
    TDirectory* dir = f.mkdir(strdir);
    dir->cd(); 
  } 

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetTrackPt(i)->GetEntries()) fCorrFF[c]->GetTrackPt(i)->Write();
    for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetZ(i)->GetEntries())       fCorrFF[c]->GetZ(i)->Write();
    for(Int_t c=0; c<fNCorrectionLevels; c++) if(fCorrFF[c]->GetXi(i)->GetEntries())      fCorrFF[c]->GetXi(i)->Write();

    for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetTrackPt(i)->GetEntries()) fCorrBgr[c]->GetTrackPt(i)->Write();
    for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetZ(i)->GetEntries())       fCorrBgr[c]->GetZ(i)->Write();
    for(Int_t c=0; c<fNCorrectionLevelsBgr; c++) if(fCorrBgr[c]->GetXi(i)->GetEntries())      fCorrBgr[c]->GetXi(i)->Write();

    
    if(fh1FFTrackPtBackFolded[i] && fh1FFTrackPtBackFolded[i]->GetEntries()) fh1FFTrackPtBackFolded[i]->Write();
    if(fh1FFZBackFolded[i]       && fh1FFZBackFolded[i]->GetEntries())       fh1FFZBackFolded[i]->Write();
    if(fh1FFXiBackFolded[i]      && fh1FFXiBackFolded[i]->GetEntries())      fh1FFXiBackFolded[i]->Write();


    if(fh1FFRatioTrackPtFolded[i] && fh1FFRatioTrackPtFolded[i]->GetEntries()) fh1FFRatioTrackPtFolded[i]->Write();
    if(fh1FFRatioZFolded[i]       && fh1FFRatioZFolded[i]->GetEntries())       fh1FFRatioZFolded[i]->Write();
    if(fh1FFRatioXiFolded[i]      && fh1FFRatioXiFolded[i]->GetEntries())      fh1FFRatioXiFolded[i]->Write();

    if(fh1FFRatioTrackPtBackFolded[i] && fh1FFRatioTrackPtBackFolded[i]->GetEntries()) fh1FFRatioTrackPtBackFolded[i]->Write();
    if(fh1FFRatioZBackFolded[i] && fh1FFRatioZBackFolded[i]->GetEntries())             fh1FFRatioZBackFolded[i]->Write();
    if(fh1FFRatioXiBackFolded[i] &&fh1FFRatioXiBackFolded[i]->GetEntries())            fh1FFRatioXiBackFolded[i]->Write();

    if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetEntries()) fh1FoldingCorrPt[i]->Write(); // TEST!!!
    if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetEntries())     fh1FoldingCorrZ[i]->Write();  // TEST!!!
    if(fh1FoldingCorrXi[i] && fh1FoldingCorrXi[i]->GetEntries()) fh1FoldingCorrXi[i]->Write(); // TEST!!!
  }
  
  // inclusive track pt

  for(Int_t c=0; c<fNCorrectionLevelsSinglePt; c++) if(fCorrSinglePt[c]->GetTrackPt(0)->GetEntries()) fCorrSinglePt[c]->GetTrackPt(0)->Write();
  if(fh1SingleTrackPtBackFolded)      fh1SingleTrackPtBackFolded->Write();  
  if(fh1RatioSingleTrackPtFolded)     fh1RatioSingleTrackPtFolded->Write();  
  if(fh1RatioSingleTrackPtBackFolded) fh1RatioSingleTrackPtBackFolded->Write();  

  


  f.Close();  
}

//____________________________________________________________________________________________________________________________________
THnSparse* AliFragmentationFunctionCorrections::TH1toSparse(const TH1F* hist, TString strName, TString strTit, const Bool_t fillConst)
{
  // copy 1-dimensional histo to THnSparse 
  // if fillConst TRUE, create THnSparse with same binning as hist but all entries = 1
  // histos with variable bin size are supported

  // note: function returns pointer to 'new' THnSparse on heap, object needs to be deleted by user

  THnSparseF* fhnHist;

  Int_t nBins       = hist->GetXaxis()->GetNbins();
  Int_t nBinsVec[1] = { nBins };

  const Double_t* binsVec = hist->GetXaxis()->GetXbins()->GetArray();
  
  if(binsVec){ // binsVec only neq NULL if histo was rebinned before

    fhnHist = new THnSparseF(strName,strTit,1,nBinsVec/*,binMinVec,binMaxVec*/);
    fhnHist->SetBinEdges(0,binsVec);
  }
  else{ // uniform bin size
    
    Double_t xMin  = hist->GetXaxis()->GetXmin();
    Double_t xMax  = hist->GetXaxis()->GetXmax();
    
    Double_t binMinVec[1] = { xMin };
    Double_t binMaxVec[1] = { xMax };
    
    fhnHist = new THnSparseF(strName,strTit,1,nBinsVec,binMinVec,binMaxVec);
  }
 

  for(Int_t bin=1; bin<=nBins; bin++){ 

    Double_t binCenter = fhnHist->GetAxis(0)->GetBinCenter(bin);

    Double_t binCoord[] = {binCenter}; 
    fhnHist->Fill(binCoord,1); // initially need to create the bin 

    Long64_t binIndex = fhnHist->GetBin(binCoord);

    Double_t cont =  hist->GetBinContent(bin);
    Double_t err  =  hist->GetBinError(bin);
    
    if(fillConst){
      cont = 1;
      err  = 0;
    }

    fhnHist->SetBinContent(binIndex,cont);
    fhnHist->SetBinError(binIndex,err);
 }

 return fhnHist;
}

//______________________________________________________________________________________________________________________________________________
TH1F* AliFragmentationFunctionCorrections::Unfold(THnSparse* hnHist, const THnSparse* hnResponse, const THnSparse* hnEff, const Int_t nIter,
						  const Bool_t useCorrelatedErrors, const THnSparse* hnPrior)
{
  // unfold input histo 

  AliCFUnfolding unfolding("unfolding","",1,hnResponse,hnEff,hnHist,hnPrior);  // arg3: nVar; hnEff required, hnPrior not (defaults to 0x0)
  unfolding.SetMaxNumberOfIterations(nIter);
  // unfolding.SetMaxChi2PerDOF(1.e-07); // OBSOLETE !!!
  // if(useSmoothing) unfolding.UseSmoothing(); 

  if(useCorrelatedErrors) unfolding.SetUseCorrelatedErrors();
  unfolding.Unfold();

  THnSparse* unfolded = unfolding.GetUnfolded();

  TString hnameUnf = hnHist->GetName();
  hnameUnf.Append("_unf");
  unfolded->SetNameTitle(hnameUnf,hnHist->GetTitle());
  
  TH1F* h1Unfolded = (TH1F*) unfolded->Projection(0); 
  h1Unfolded->SetNameTitle(hnHist->GetName(),hnHist->GetTitle());
  
  return h1Unfolded;
}

//___________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::UnfoldHistos(const Int_t nIter, const Bool_t useCorrelatedErrors, const Int_t type)
{
  // loop over jet pt slices and unfold dN/dpt spectra
  
  TString labelCorr = fCorrFF[fNCorrectionLevels-1]->GetLabel(); 
  if(!labelCorr.Contains("unfold")) AddCorrectionLevel("unfold");
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
 
    TH1F* hist = 0;
    if(type == kFlagPt)      hist = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i); // level -2: before unfolding, level -1: unfolded
    else if(type == kFlagZ)  hist = fCorrFF[fNCorrectionLevels-2]->GetZ(i);       // level -2: before unfolding, level -1: unfolded
    else if(type == kFlagXi) hist = fCorrFF[fNCorrectionLevels-2]->GetXi(i);      // level -2: before unfolding, level -1: unfolded
    else{ 
      Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
      return;
    }
    
    if(!hist){
      Printf("%s%d no histo found ",(char*)__FILE__,__LINE__);
      return;
    }


    THnSparse* hnResponse = 0;
    if(type == kFlagPt)      hnResponse = fhnResponsePt[i];
    else if(type == kFlagZ)  hnResponse = fhnResponseZ[i];
    else if(type == kFlagXi) hnResponse = fhnResponseXi[i];

    TH1F* hPrior = 0;
    if(type == kFlagPt && fh1FFTrackPtPrior[i]  && ((TString(fh1FFTrackPtPrior[i]->GetName())).Length() > 0) ) hPrior = fh1FFTrackPtPrior[i];
    else if(type == kFlagZ  && fh1FFZPrior[i]   && ((TString(fh1FFZPrior[i]->GetName())).Length() > 0)       ) hPrior = fh1FFZPrior[i];
    else if(type == kFlagXi && fh1FFXiPrior[i]  && ((TString(fh1FFXiPrior[i]->GetName())).Length() > 0)      ) hPrior = fh1FFXiPrior[i];


    TString histNameTHn;
    histNameTHn = hist->GetName();
    histNameTHn.ReplaceAll("TH1","THn");
    

    TString priorNameTHn; 
    if(hPrior){
      priorNameTHn = hPrior->GetName();
      priorNameTHn.ReplaceAll("TH1","THn");
    }

    TString histNameBackFolded;
    histNameBackFolded = hist->GetName();
    histNameBackFolded.Append("_backfold");
    
    TString histNameRatioFolded;
    histNameRatioFolded = hist->GetName();
    histNameRatioFolded.ReplaceAll("fh1FF","hRatioFF");
    
    histNameRatioFolded.Append("_unfold");

    TString histNameRatioBackFolded;
    histNameRatioBackFolded = hist->GetName();
    histNameRatioBackFolded.ReplaceAll("fh1FF","hRatioFF");
    histNameRatioBackFolded.Append("_backfold");
 
    THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
    THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff 
    THnSparse* hnPrior          = 0;
    if(hPrior) hnPrior = TH1toSparse(hPrior,priorNameTHn,hPrior->GetTitle());
    
    //THnSparse* hnUnfolded 
    //  = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);  
     
    //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
    //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
    
    TH1F* hUnfolded 
      = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors,hnPrior);  
     
    //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
    //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());

    // errors
    for(Int_t bin=1; bin<=hist->GetNbinsX(); bin++){
      
      Double_t contOrg = hist->GetBinContent(bin);
      Double_t errOrg  = hist->GetBinError(bin);
      
      Double_t contUnf = hUnfolded->GetBinContent(bin);
      //Double_t errUnf  = hUnfolded->GetBinError(bin);

      Double_t relErrOrg = 0;
      if(contOrg) relErrOrg = errOrg/contOrg;
      
      //       std::cout<<" cont original "<<contOrg<<" err original "<<errOrg<<std::endl;
      //       std::cout<<" cont unf "<<contUnf<<" errUnf"<<errUnf<<std::endl;
      //       std::cout<<" err/cont original "<<relErrOrg<<std::endl;
      //       if(contUnf)  std::cout<<" err/conf unf "<<errUnf/contUnf<<std::endl;
      
      hUnfolded->SetBinError(bin,relErrOrg*contUnf);
      
      //       if(hUnfolded->GetBinContent(bin)){
      // 	std::cout<<" modified err unfolded "<<hUnfolded->GetBinError(bin)/hUnfolded->GetBinContent(bin)<<std::endl; 
      //       }
    }
    

    if(type == kFlagPt) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hUnfolded,0,0);
    if(type == kFlagZ)  fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,hUnfolded,0);
    if(type == kFlagXi) fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,0,0,hUnfolded);

    // backfolding: apply response matrix to unfolded spectrum
    TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); 
    hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());

    if(type == kFlagPt) fh1FFTrackPtBackFolded[i] = hBackFolded;
    if(type == kFlagZ)  fh1FFZBackFolded[i]       = hBackFolded;
    if(type == kFlagXi) fh1FFXiBackFolded[i]      = hBackFolded;
    
    // ratio unfolded to original histo 
    TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
    hRatioUnfolded->Reset();
    hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");

    if(type == kFlagPt) fh1FFRatioTrackPtFolded[i] = hRatioUnfolded;
    if(type == kFlagZ)  fh1FFRatioZFolded[i]       = hRatioUnfolded;
    if(type == kFlagXi) fh1FFRatioXiFolded[i]      = hRatioUnfolded;


    // ratio backfolded to original histo
    TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
    hRatioBackFolded->Reset();
    hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");

    if(type == kFlagPt) fh1FFRatioTrackPtBackFolded[i] = hRatioBackFolded;
    if(type == kFlagZ)  fh1FFRatioZBackFolded[i]       = hRatioBackFolded;
    if(type == kFlagXi) fh1FFRatioXiBackFolded[i]      = hRatioBackFolded;
    
    delete hnHist;
    delete hnFlatEfficiency;
 }
}

//_____________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::UnfoldPt(const Int_t nIter, const Bool_t useCorrelatedErrors)
{
  
  Int_t type = kFlagPt;
  UnfoldHistos(nIter, useCorrelatedErrors, type);
}

//_____________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::UnfoldZ(const Int_t nIter, const Bool_t useCorrelatedErrors)
{
  
  Int_t type = kFlagZ;
  UnfoldHistos(nIter, useCorrelatedErrors, type);
}

//_____________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::UnfoldXi(const Int_t nIter, const Bool_t useCorrelatedErrors)
{
  
  Int_t type = kFlagXi;
  UnfoldHistos(nIter, useCorrelatedErrors, type);
}


//______________________________________________________________________________________________
TH1F* AliFragmentationFunctionCorrections::ApplyResponse(const TH1F* hist, THnSparse* hnResponse)
{
  // apply (multiply) response matrix to hist 

  TH1F* hOut = new TH1F(*hist);
  hOut->Reset();

  const Int_t axisM = 0; 
  const Int_t axisT = 1;
 
  Int_t nBinsM = hnResponse->GetAxis(axisM)->GetNbins();
  Int_t nBinsT = hnResponse->GetAxis(axisT)->GetNbins();

  // response matrix normalization
  // do it in this function and not when reading response, since for 'proper' normalization errors are difficult to assign
  // so for unfolding proper we leave it to CORRFW ...

  Double_t normFacResponse[nBinsT+1];

  for(Int_t iT=1; iT<=nBinsT; iT++){

    Double_t sumResp = 0;
    
    for(Int_t iM=1; iM<=nBinsM; iM++){
      
      Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
      Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);
      
      Double_t binCoord[] = {coordM,coordT};
      
      Long64_t binIndex = hnResponse->GetBin(binCoord);
      
      Double_t resp = hnResponse->GetBinContent(binIndex); 
      
      sumResp += resp;
    }
    
    normFacResponse[iT] = 0;
    if(sumResp) normFacResponse[iT] = 1/sumResp;
  }
  
  
  
  for(Int_t iM=1; iM<=nBinsM; iM++){
    
    Double_t contM   = 0;

    for(Int_t iT=1; iT<=nBinsT; iT++){

      Double_t contT = hist->GetBinContent(iT);
      
      Double_t coordM = hnResponse->GetAxis(axisM)->GetBinCenter(iM);
      Double_t coordT = hnResponse->GetAxis(axisT)->GetBinCenter(iT);

      Double_t binCoord[] = {coordM,coordT};
      
      Long64_t binIndex = hnResponse->GetBin(binCoord);
      
      Double_t resp = hnResponse->GetBinContent(binIndex); 
      
      contM   += resp*normFacResponse[iT]*contT; 
    }

    hOut->SetBinContent(iM,contM);
  }

  return hOut;
}

//_______________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadEfficiency(TString strfile, TString strdir, TString strlist)
{
  // read reconstruction efficiency from file
  // argument strlist optional - read from directory strdir if not specified

  // temporary histos to hold histos from file
  TH1F* hEffPt[fNJetPtSlices]; 
  TH1F* hEffZ[fNJetPtSlices];
  TH1F* hEffXi[fNJetPtSlices];
  
  for(Int_t i=0; i<fNJetPtSlices; i++) hEffPt[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hEffZ[i]  = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hEffXi[i] = 0;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
   
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
    
    TString strNameEffPt(Form("hEffPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameEffZ(Form("hEffZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameEffXi(Form("hEffXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
     
   
    if(list){
      hEffPt[i] = (TH1F*) list->FindObject(strNameEffPt); 
      hEffZ[i]  = (TH1F*) list->FindObject(strNameEffZ); 
      hEffXi[i] = (TH1F*) list->FindObject(strNameEffXi); 
    }
    else{
      hEffPt[i] = (TH1F*) gDirectory->Get(strNameEffPt); 
      hEffZ[i]  = (TH1F*) gDirectory->Get(strNameEffZ); 
      hEffXi[i] = (TH1F*) gDirectory->Get(strNameEffXi); 
    }
    
    if(!hEffPt[i]){
      Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPt.Data());
    }
  
    if(!hEffZ[i]){
      Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZ.Data());
    }    

    if(!hEffXi[i]){
      Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXi.Data());
    }


    if(fNHistoBinsPt[i]) hEffPt[i] = (TH1F*) hEffPt[i]->Rebin(fNHistoBinsPt[i],strNameEffPt+"_rebin",fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hEffZ[i]  = (TH1F*) hEffZ[i]->Rebin(fNHistoBinsZ[i],strNameEffZ+"_rebin",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hEffXi[i] = (TH1F*) hEffXi[i]->Rebin(fNHistoBinsXi[i],strNameEffXi+"_rebin",fHistoBinsXi[i]->GetArray());

    if(hEffPt[i]) hEffPt[i]->SetDirectory(0); 
    if(hEffZ[i])  hEffZ[i]->SetDirectory(0); 
    if(hEffXi[i]) hEffXi[i]->SetDirectory(0); 

  } // jet slices loop

  f.Close();

  for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
    if(hEffPt[i]) new(fh1EffPt[i]) TH1F(*hEffPt[i]);
    if(hEffZ[i])  new(fh1EffZ[i])  TH1F(*hEffZ[i]);
    if(hEffXi[i]) new(fh1EffXi[i]) TH1F(*hEffXi[i]);
  }
}

//___________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadBgrEfficiency(TString strfile, TString strdir, TString strlist)
{
  // read bgr eff from file
  // argument strlist optional - read from directory strdir if not specified
 
  TH1F* hEffPtBgr[fNJetPtSlices]; 
  TH1F* hEffZBgr [fNJetPtSlices];
  TH1F* hEffXiBgr[fNJetPtSlices];

  for(Int_t i=0; i<fNJetPtSlices; i++) hEffPtBgr[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hEffZBgr[i]  = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hEffXiBgr[i] = 0;


  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read bgr efficiencies from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
       
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
    
    TString strNameEffPtBgr(Form("hEffPtBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
    TString strNameEffZBgr(Form("hEffZBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
    TString strNameEffXiBgr(Form("hEffXiBgr%02dto%02d",jetPtLoLim,jetPtUpLim));
  
   
    if(list){
      hEffPtBgr[i] = (TH1F*) list->FindObject(strNameEffPtBgr); 
      hEffZBgr[i]  = (TH1F*) list->FindObject(strNameEffZBgr); 
      hEffXiBgr[i] = (TH1F*) list->FindObject(strNameEffXiBgr); 
    }
    else{
      hEffPtBgr[i] = (TH1F*) gDirectory->Get(strNameEffPtBgr); 
      hEffZBgr[i]  = (TH1F*) gDirectory->Get(strNameEffZBgr); 
      hEffXiBgr[i] = (TH1F*) gDirectory->Get(strNameEffXiBgr); 
    }
    
    if(!hEffPtBgr[i]){
      Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffPtBgr.Data());
    }
  
    if(!hEffZBgr[i]){
      Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffZBgr.Data());
    }    

    if(!hEffXiBgr[i]){
      Printf("%s:%d -- error retrieving efficiency %s", (char*)__FILE__,__LINE__,strNameEffXiBgr.Data());
    }


    if(fNHistoBinsPt[i]) hEffPtBgr[i] = (TH1F*) hEffPtBgr[i]->Rebin(fNHistoBinsPt[i],strNameEffPtBgr+"_rebin",fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hEffZBgr[i]  = (TH1F*) hEffZBgr[i]->Rebin(fNHistoBinsZ[i],strNameEffZBgr+"_rebin",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hEffXiBgr[i] = (TH1F*) hEffXiBgr[i]->Rebin(fNHistoBinsXi[i],strNameEffXiBgr+"_rebin",fHistoBinsXi[i]->GetArray());

    if(hEffPtBgr[i]) hEffPtBgr[i]->SetDirectory(0); 
    if(hEffZBgr[i])  hEffZBgr[i]->SetDirectory(0); 
    if(hEffXiBgr[i]) hEffXiBgr[i]->SetDirectory(0); 

  } // jet slices loop

  f.Close();

  for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
    if(hEffPtBgr[i]) new(fh1EffBgrPt[i]) TH1F(*hEffPtBgr[i]);
    if(hEffZBgr[i])  new(fh1EffBgrZ[i])  TH1F(*hEffZBgr[i]);
    if(hEffXiBgr[i]) new(fh1EffBgrXi[i]) TH1F(*hEffXiBgr[i]);
  }
}

// ________________________________________________
void AliFragmentationFunctionCorrections::EffCorr()
{
  // apply efficiency correction

  AddCorrectionLevel("eff");

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
    TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();

    
    TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
    hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
    
    TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
    hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
    
    TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
    hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
    
    fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
  }
}

//___________________________________________________
void AliFragmentationFunctionCorrections::EffCorrBgr()
{
  // apply efficiency correction to bgr distributions

  AddCorrectionLevelBgr("eff");

  Printf("%s:%d -- apply efficiency correction, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
    TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);
    
    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();

    
    TH1F* hFFTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
    hFFTrackPtEffCorr->Divide(histPt,fh1EffPt[i],1,1,"");
    
    TH1F* hFFZEffCorr = (TH1F*) histZ->Clone(histNameZ);
    hFFZEffCorr->Divide(histZ,fh1EffZ[i],1,1,"");
    
    TH1F* hFFXiEffCorr = (TH1F*) histXi->Clone(histNameXi);
    hFFXiEffCorr->Divide(histXi,fh1EffXi[i],1,1,"");
    
    fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtEffCorr,hFFZEffCorr,hFFXiEffCorr);
  }
}

//_____________________________________________________
void AliFragmentationFunctionCorrections::SubtractBgr(Double_t sysErr)
{
  // subtract bgr distribution from FF
  // the current corr level is used for both 
  
  AddCorrectionLevel("bgrSub");

  fstream ascii_out_pt;
  fstream ascii_out_z;
  fstream ascii_out_xi;

  if(sysErr) ascii_out_pt.open("sysErrUE_pt.txt",std::ios::out);
  if(sysErr) ascii_out_z.open("sysErrUE_z.txt",std::ios::out);
  if(sysErr) ascii_out_xi.open("sysErrUE_xi.txt",std::ios::out);

  for(Int_t i=0; i<fNJetPtSlices; i++){ // jet slices

    TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
    TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);
    
    TH1F* histPtBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetTrackPt(i);
    TH1F* histZBgr  = fCorrBgr[fNCorrectionLevelsBgr-1]->GetZ(i);
    TH1F* histXiBgr = fCorrBgr[fNCorrectionLevelsBgr-1]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();
    

    TH1F* hFFTrackPtBgrSub = (TH1F*) histPt->Clone(histNamePt);
    hFFTrackPtBgrSub->Add(histPtBgr,-1);
    
    TH1F* hFFZBgrSub =  (TH1F*) histZ->Clone(histNameZ);
    hFFZBgrSub->Add(histZBgr,-1);
    
    TH1F* hFFXiBgrSub = (TH1F*) histXi->Clone(histNameXi);
    hFFXiBgrSub->Add(histXiBgr,-1);

    fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBgrSub,hFFZBgrSub,hFFXiBgrSub);

    if(sysErr){

      for(Int_t bin=1; bin<=histPt->GetNbinsX(); bin++){

	//Double_t sigPlBgr = histPt->GetBinContent(bin);
	Double_t sig      = hFFTrackPtBgrSub->GetBinContent(bin);
	Double_t bgr      = histPtBgr->GetBinContent(bin);
	
	Double_t relErr = 0; 
	if(sig>0)  relErr = sysErr*bgr/sig;

// 	std::cout<<" pt bin "<<bin<<" mean "<<histPt->GetBinCenter(bin)
// 		 <<" sigPlBgr "<<sigPlBgr<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
	
	ascii_out_pt<<i<<" "<<bin<<" "<<relErr<<std::endl;
      }


      for(Int_t bin=1; bin<=histZ->GetNbinsX(); bin++){
	
	//Double_t sigPlBgr = histZ->GetBinContent(bin);
	Double_t sig      = hFFZBgrSub->GetBinContent(bin);
	Double_t bgr      = histZBgr->GetBinContent(bin);
	
	Double_t relErr = 0; 
	if(sig>0)  relErr = sysErr*bgr/sig;
	
 	std::cout<<" z bin "<<bin<<" mean "<<histZ->GetBinCenter(bin)
 		 <<" sigPlBgr "<<histZ->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
	
	ascii_out_z<<i<<" "<<bin<<" "<<relErr<<std::endl;
      }


      for(Int_t bin=1; bin<=histXi->GetNbinsX(); bin++){
	
	//Double_t sigPlBgr = histXi->GetBinContent(bin);
	Double_t sig      = hFFXiBgrSub->GetBinContent(bin);
	Double_t bgr      = histXiBgr->GetBinContent(bin);
	
	Double_t relErr = 0; 
	if(sig>0)  relErr = sysErr*bgr/sig;
	
 	std::cout<<" xi bin "<<bin<<" mean "<<histXi->GetBinCenter(bin)
 		 <<" sigPlBgr "<<histXi->GetBinContent(bin)<<" sig "<<sig<<" bgr "<<bgr<<" relErr "<<relErr<<std::endl;
	
	ascii_out_xi<<i<<" "<<bin<<" "<<relErr<<std::endl;
      }
    }
  }

  if(sysErr) ascii_out_pt.close(); 
  if(sysErr) ascii_out_xi.close(); 


}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strID, TString strOutfile, 
							      Bool_t updateOutfile, TString strOutDir,TString strPostfix)
{ 
  // read task ouput from MC and write single track eff - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteSingleTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir,strPostfix);
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteSingleTrackEff(TString strInfile, TString strdir, TString strlist, 
							      TString strOutfile, Bool_t updateOutfile, TString strOutDir,TString strPostfix)
{
  // read task output from MC and write single track eff as function of pt, eta, phi
  // argument strlist optional - read from directory strdir if not specified
  // write eff to file stroutfile - by default only 'update' file (don't overwrite)


  TH1D* hdNdptTracksMCPrimGen;
  TH2D* hdNdetadphiTracksMCPrimGen;
  
  TH1D* hdNdptTracksMCPrimRec;
  TH2D* hdNdetadphiTracksMCPrimRec;
    

  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()){
    if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: change dir to %s",(char*)__FILE__,__LINE__,strdir.Data());
    gDirectory->cd(strdir);
  }

  TList* list = 0;

  if(strlist && strlist.Length()){

    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }


  TString hnamePtRecEffGen = "fh1TrackQAPtRecEffGen";
  if(strPostfix.Length()) hnamePtRecEffGen.Form("fh1TrackQAPtRecEffGen_%s",strPostfix.Data()); 

  TString hnamePtRecEffRec = "fh1TrackQAPtRecEffRec";
  if(strPostfix.Length()) hnamePtRecEffRec.Form("fh1TrackQAPtRecEffRec_%s",strPostfix.Data());

  TString hnameEtaPhiRecEffGen = "fh2TrackQAEtaPhiRecEffGen";
  if(strPostfix.Length()) hnameEtaPhiRecEffGen.Form("fh2TrackQAEtaPhiRecEffGen_%s",strPostfix.Data());
  
  TString hnameEtaPhiRecEffRec = "fh2TrackQAEtaPhiRecEffRec";
  if(strPostfix.Length()) hnameEtaPhiRecEffRec.Form("fh2TrackQAEtaPhiRecEffRec_%s",strPostfix.Data());
  

  if(list){
    hdNdptTracksMCPrimGen       = (TH1D*) list->FindObject(hnamePtRecEffGen);
    hdNdetadphiTracksMCPrimGen  = (TH2D*) list->FindObject(hnameEtaPhiRecEffGen);
    
    hdNdptTracksMCPrimRec       = (TH1D*) list->FindObject(hnamePtRecEffRec);
    hdNdetadphiTracksMCPrimRec  = (TH2D*) list->FindObject(hnameEtaPhiRecEffRec);
  }
  else{
    hdNdptTracksMCPrimGen       = (TH1D*) gDirectory->Get(hnamePtRecEffGen);
    hdNdetadphiTracksMCPrimGen  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffGen);

    hdNdptTracksMCPrimRec       = (TH1D*) gDirectory->Get(hnamePtRecEffRec);
    hdNdetadphiTracksMCPrimRec  = (TH2D*) gDirectory->Get(hnameEtaPhiRecEffRec);
  }

  hdNdptTracksMCPrimGen->SetDirectory(0);
  hdNdetadphiTracksMCPrimGen->SetDirectory(0);
  hdNdptTracksMCPrimRec->SetDirectory(0);
  hdNdetadphiTracksMCPrimRec->SetDirectory(0);
  
  f.Close();

  // projections: dN/deta, dN/dphi 

  TH1D* hdNdetaTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionX("hdNdetaTracksMcPrimGen");
  TH1D* hdNdphiTracksMCPrimGen = (TH1D*) hdNdetadphiTracksMCPrimGen->ProjectionY("hdNdphiTracksMcPrimGen");
 
  TH1D* hdNdetaTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
  TH1D* hdNdphiTracksMCPrimRec = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");

  // rebin

  TString strNamePtGen = "hTrackPtGenPrim";
  TString strNamePtRec = "hTrackPtRecPrim";

  if(fNHistoBinsSinglePt) hdNdptTracksMCPrimGen = (TH1D*) hdNdptTracksMCPrimGen->Rebin(fNHistoBinsSinglePt,strNamePtGen,fHistoBinsSinglePt->GetArray());
  if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtRec,fHistoBinsSinglePt->GetArray());
 
    // efficiency: divide 

  TString hNameTrackEffPt = "hSingleTrackEffPt";
  if(strPostfix.Length()) hNameTrackEffPt.Form("hSingleTrackEffPt_%s",strPostfix.Data());
					       
  TString hNameTrackEffEta = "hSingleTrackEffEta";
  if(strPostfix.Length()) hNameTrackEffEta.Form("hSingleTrackEffEta_%s",strPostfix.Data());

  TString hNameTrackEffPhi = "hSingleTrackEffPhi";
  if(strPostfix.Length()) hNameTrackEffPhi.Form("hSingleTrackEffPhi_%s",strPostfix.Data());


  TH1F* hSingleTrackEffPt = (TH1F*) hdNdptTracksMCPrimRec->Clone(hNameTrackEffPt);
  hSingleTrackEffPt->Divide(hdNdptTracksMCPrimRec,hdNdptTracksMCPrimGen,1,1,"B"); // binominal errors

  TH1F* hSingleTrackEffEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone(hNameTrackEffEta);
  hSingleTrackEffEta->Divide(hdNdetaTracksMCPrimRec,hdNdetaTracksMCPrimGen,1,1,"B"); // binominal errors
  
  TH1F* hSingleTrackEffPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone(hNameTrackEffPhi);
  hSingleTrackEffPhi->Divide(hdNdphiTracksMCPrimRec,hdNdphiTracksMCPrimGen,1,1,"B"); // binominal errors
  
  
  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);

  if(!out.IsOpen()){
    Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

  if(strOutDir && strOutDir.Length()){
    
    TDirectory* dir;
    if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }

  hSingleTrackEffPt->Write();
  hSingleTrackEffEta->Write();
  hSingleTrackEffPhi->Write();
  
  out.Close();

  delete hdNdptTracksMCPrimGen;
  delete hdNdetadphiTracksMCPrimGen;
  delete hdNdptTracksMCPrimRec;
  delete hdNdetadphiTracksMCPrimRec;
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strID, TString strOutfile, 
								  Bool_t updateOutfile, TString strOutDir)
{ 
  // read task ouput from MC and write single track eff - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteSingleTrackSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteSingleTrackSecCorr(TString strInfile, TString strdir, TString strlist, 
								  TString strOutfile, Bool_t updateOutfile, TString strOutDir)
{
  // read task output from MC and write single track secondaries contamination correction as function of pt, eta, phi
  // argument strlist optional - read from directory strdir if not specified
  // write corr factor to file stroutfile - by default only 'update' file (don't overwrite)

  TH1D* hdNdptTracksMCPrimRec;
  TH2D* hdNdetadphiTracksMCPrimRec;
  
  TH1D* hdNdptTracksMCSecRecNS;
  TH2D* hdNdetadphiTracksMCSecRecNS;

  TH1D* hdNdptTracksMCSecRecSsc;
  TH2D* hdNdetadphiTracksMCSecRecSsc;

  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeSingleTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){

    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }


  if(list){
    hdNdptTracksMCPrimRec        = (TH1D*) list->FindObject("fh1TrackQAPtRecEffGen");
    hdNdetadphiTracksMCPrimRec   = (TH2D*) list->FindObject("fh2TrackQAEtaPhiRecEffGen");
    
    hdNdptTracksMCSecRecNS       = (TH1D*) list->FindObject("fh1TrackQAPtSecRecNS");
    hdNdetadphiTracksMCSecRecNS  = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecNS");

    hdNdptTracksMCSecRecSsc      = (TH1D*) list->FindObject("fh1TrackQAPtSecRecSsc");
    hdNdetadphiTracksMCSecRecSsc = (TH2D*) list->FindObject("fh2TrackQAEtaPhiSecRecSsc");

  }
  else{
    hdNdptTracksMCPrimRec        = (TH1D*) gDirectory->Get("fh1TrackQAPtRecEffGen");
    hdNdetadphiTracksMCPrimRec   = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiRecEffGen");

    hdNdptTracksMCSecRecNS       = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecNS");
    hdNdetadphiTracksMCSecRecNS  = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecNS");

    hdNdptTracksMCSecRecSsc      = (TH1D*) gDirectory->Get("fh1TrackQAPtSecRecSsc");
    hdNdetadphiTracksMCSecRecSsc = (TH2D*) gDirectory->Get("fh2TrackQAEtaPhiSecRecSsc");
  }
  
  hdNdptTracksMCPrimRec->SetDirectory(0);
  hdNdetadphiTracksMCPrimRec->SetDirectory(0);

  hdNdptTracksMCSecRecNS->SetDirectory(0);
  hdNdetadphiTracksMCSecRecNS->SetDirectory(0);

  hdNdptTracksMCSecRecSsc->SetDirectory(0);
  hdNdetadphiTracksMCSecRecSsc->SetDirectory(0);
  
  f.Close();

  // projections: dN/deta, dN/dphi 

  TH1D* hdNdetaTracksMCPrimRec   = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionX("hdNdetaTracksMcPrimRec");
  TH1D* hdNdphiTracksMCPrimRec   = (TH1D*) hdNdetadphiTracksMCPrimRec->ProjectionY("hdNdphiTracksMcPrimRec");
 
  TH1D* hdNdetaTracksMCSecRecNS  = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionX("hdNdetaTracksMcSecRecNS");
  TH1D* hdNdphiTracksMCSecRecNS  = (TH1D*) hdNdetadphiTracksMCSecRecNS->ProjectionY("hdNdphiTracksMcSecRecNS");

  TH1D* hdNdetaTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionX("hdNdetaTracksMcSecRecSsc");
  TH1D* hdNdphiTracksMCSecRecSsc = (TH1D*) hdNdetadphiTracksMCSecRecSsc->ProjectionY("hdNdphiTracksMcSecRecSsc");
 
  // rebin

  TString strNamePtPrim = "hTrackPtPrim";
  TString strNamePtSec  = "hTrackPtSec";

  if(fNHistoBinsSinglePt) hdNdptTracksMCPrimRec   = (TH1D*) hdNdptTracksMCPrimRec->Rebin(fNHistoBinsSinglePt,strNamePtPrim,fHistoBinsSinglePt->GetArray());
  if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecNS  = (TH1D*) hdNdptTracksMCSecRecNS->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());
  if(fNHistoBinsSinglePt) hdNdptTracksMCSecRecSsc = (TH1D*) hdNdptTracksMCSecRecSsc->Rebin(fNHistoBinsSinglePt,strNamePtSec,fHistoBinsSinglePt->GetArray());


  // secondary correction factor: divide prim/(prim+sec)

  TH1F* hSingleTrackSecCorrPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSingleTrackSecCorrPt");
  TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec->Clone("hSumPrimSecPt");
  hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS);
  hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc);
  hSingleTrackSecCorrPt->Divide(hdNdptTracksMCPrimRec,hSumPrimSecPt,1,1,"B"); // binominal errors


  TH1F* hSingleTrackSecCorrEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSingleTrackSecCorrEta");
  TH1F* hSumPrimSecEta = (TH1F*) hdNdetaTracksMCPrimRec->Clone("hSumPrimSecEta");
  hSumPrimSecEta->Add(hdNdetaTracksMCSecRecNS);
  hSumPrimSecEta->Add(hdNdetaTracksMCSecRecSsc);
  hSingleTrackSecCorrEta->Divide(hdNdetaTracksMCPrimRec,hSumPrimSecEta,1,1,"B"); // binominal errors

  TH1F* hSingleTrackSecCorrPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSingleTrackSecCorrPhi");
  TH1F* hSumPrimSecPhi = (TH1F*) hdNdphiTracksMCPrimRec->Clone("hSumPrimSecPhi");
  hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecNS);
  hSumPrimSecPhi->Add(hdNdphiTracksMCSecRecSsc);
  hSingleTrackSecCorrPhi->Divide(hdNdphiTracksMCPrimRec,hSumPrimSecPhi,1,1,"B"); // binominal errors

  // 
  
  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);

  if(!out.IsOpen()){
    Printf("%s:%d -- error opening secCorr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write secCorr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

  if(strOutDir && strOutDir.Length()){  

    TDirectory* dir;
    if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  } 

  hdNdptTracksMCSecRecNS->Write();
  hdNdetaTracksMCSecRecNS->Write();
  hdNdphiTracksMCSecRecNS->Write();

  hdNdptTracksMCSecRecSsc->Write();
  hdNdetaTracksMCSecRecSsc->Write();
  hdNdphiTracksMCSecRecSsc->Write();

  hSingleTrackSecCorrPt->Write();
  hSingleTrackSecCorrEta->Write();
  hSingleTrackSecCorrPhi->Write();
  
  out.Close();

  delete hdNdptTracksMCPrimRec;       
  delete hdNdetadphiTracksMCPrimRec;  
  delete hdNdptTracksMCSecRecNS;
  delete hdNdetadphiTracksMCSecRecNS;
  delete hdNdptTracksMCSecRecSsc;
  delete hdNdetadphiTracksMCSecRecSsc;
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strID, TString strOutfile, 
							      Bool_t updateOutfile, TString strOutDir)
{ 
  // read task ouput from MC and write single track eff - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteSingleResponse(strInfile,strdir,strlist,strOutfile,updateOutfile,strOutDir);
}


//_____________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteSingleResponse(TString strInfile, TString strdir, TString strlist,
							      TString strOutfile, Bool_t updateOutfile, TString strOutDir)
{
  // read 2d THnSparse response matrices in pt from file
  // project TH2 
  // write to strOutfile 

  THnSparse* hnResponseSinglePt; 
  TH2F*      h2ResponseSinglePt; 
 
  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
    
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  
  
  if(list) hnResponseSinglePt = (THnSparse*) list->FindObject("fhnResponseSinglePt");
  else     hnResponseSinglePt = (THnSparse*) gDirectory->Get("fhnResponseSinglePt");
  

  if(!hnResponseSinglePt){
    Printf("%s:%d -- error retrieving response matrix fhnResponseSinglePt",(char*)__FILE__,__LINE__);
    return;
  }

  f.Close();


  // project 2d histo 
  h2ResponseSinglePt = (TH2F*) hnResponseSinglePt->Projection(1,0);// note convention: yDim,xDim
  h2ResponseSinglePt->SetNameTitle("h2ResponseSinglePt",""); 
    
  
  // write 

  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";
  
  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());
  
  if(strOutDir && strOutDir.Length()){  

    TDirectory* dir;
    if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }
  
  hnResponseSinglePt->Write();
  h2ResponseSinglePt->Write();
  
  out.Close();  
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strID, TString strOutfile, 
							   Bool_t updateOutfile)
{ 
  // read task ouput from MC and write single track eff - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteJetTrackEff(strInfile,strdir,strlist,strOutfile,updateOutfile);
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetTrackEff(TString strInfile, TString strdir, TString strlist, 
							   TString strOutfile, Bool_t updateOutfile)
{
  // read task output from MC and write track eff in jet pt slices as function of pt, z, xi
  // argument strlist optional - read from directory strdir if not specified
  // write eff to file strOutfile - by default only 'update' file (don't overwrite)

  TH1F* hEffPt[fNJetPtSlices];
  TH1F* hEffXi[fNJetPtSlices];
  TH1F* hEffZ[fNJetPtSlices];

  TH1F* hdNdptTracksMCPrimGen[fNJetPtSlices];
  TH1F* hdNdxiMCPrimGen[fNJetPtSlices];
  TH1F* hdNdzMCPrimGen[fNJetPtSlices];
    
  TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
  TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
  TH1F* hdNdzMCPrimRec[fNJetPtSlices];


  TH1F* fh1FFJetPtRecEffGen;

  TH2F* fh2FFTrackPtRecEffGen;
  TH2F* fh2FFZRecEffGen;
  TH2F* fh2FFXiRecEffGen;
  
  TH2F* fh2FFTrackPtRecEffRec;
  TH2F* fh2FFZRecEffRec;
  TH2F* fh2FFXiRecEffRec;
 

  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeJetTrackEff: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){

    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }

  if(list){
    fh1FFJetPtRecEffGen    = (TH1F*) list->FindObject("fh1FFJetPtRecEffGen");

    fh2FFTrackPtRecEffGen  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffGen");
    fh2FFZRecEffGen        = (TH2F*) list->FindObject("fh2FFZRecEffGen");
    fh2FFXiRecEffGen       = (TH2F*) list->FindObject("fh2FFXiRecEffGen");
    
    fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
    fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
    fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
  }
  else{
    fh1FFJetPtRecEffGen    = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffGen");

    fh2FFTrackPtRecEffGen  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffGen");
    fh2FFZRecEffGen        = (TH2F*) gDirectory->Get("fh2FFZRecEffGen");
    fh2FFXiRecEffGen       = (TH2F*) gDirectory->Get("fh2FFXiRecEffGen");
    
    fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
    fh2FFZRecEffRec        = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
    fh2FFXiRecEffRec       = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
  }
  
  fh1FFJetPtRecEffGen->SetDirectory(0); 

  fh2FFTrackPtRecEffGen->SetDirectory(0);
  fh2FFZRecEffGen->SetDirectory(0);
  fh2FFXiRecEffGen->SetDirectory(0);
  
  fh2FFTrackPtRecEffRec->SetDirectory(0);
  fh2FFZRecEffRec->SetDirectory(0);
  fh2FFXiRecEffRec->SetDirectory(0);
  
  f.Close();


  // projections: FF for generated and reconstructed primaries 
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);

    Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffGen->GetXaxis()->FindBin(jetPtUpLim))-1;

    if(binUp > fh2FFTrackPtRecEffGen->GetNbinsX()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameFFPtGen(Form("fh1FFTrackPtGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZGen(Form("fh1FFZGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiGen(Form("fh1FFXiGenPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    TString strNameFFPtRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    // project 
    // appendix 'unbinned' to avoid histos with same name after rebinning

    hdNdptTracksMCPrimGen[i] = (TH1F*) fh2FFTrackPtRecEffGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCPrimGen[i]        = (TH1F*) fh2FFZRecEffGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
    hdNdxiMCPrimGen[i]       = (TH1F*) fh2FFXiRecEffGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
    
    hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
    hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
    
    // rebin

    if(fNHistoBinsPt[i]) hdNdptTracksMCPrimGen[i] = (TH1F*) hdNdptTracksMCPrimGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCPrimGen[i]  = (TH1F*) hdNdzMCPrimGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCPrimGen[i] = (TH1F*) hdNdxiMCPrimGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());

    if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]  = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i] = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());

    hdNdptTracksMCPrimGen[i]->SetNameTitle(strNameFFPtGen,"");
    hdNdzMCPrimGen[i]->SetNameTitle(strNameFFZGen,"");
    hdNdxiMCPrimGen[i]->SetNameTitle(strNameFFXiGen,"");
    
    hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtRec,"");
    hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZRec,"");
    hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiRec,"");
 
    // normalize
    
    Double_t nJetsBin = fh1FFJetPtRecEffGen->Integral(binLo,binUp);

    NormalizeTH1(hdNdptTracksMCPrimGen[i],nJetsBin); 
    NormalizeTH1(hdNdzMCPrimGen[i],nJetsBin); 
    NormalizeTH1(hdNdxiMCPrimGen[i],nJetsBin); 

    NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
    NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
    NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 
    
    // divide rec/gen : efficiency

    TString strNameEffPt(Form("hEffPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameEffZ(Form("hEffZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameEffXi(Form("hEffXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
 
    hEffPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameEffPt);
    hEffPt[i]->Divide(hdNdptTracksMCPrimRec[i],hdNdptTracksMCPrimGen[i],1,1,"B"); // binominal errors
    
    hEffXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameEffXi);
    hEffXi[i]->Divide(hdNdxiMCPrimRec[i],hdNdxiMCPrimGen[i],1,1,"B"); // binominal errors
    
    hEffZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameEffZ);
    hEffZ[i]->Divide(hdNdzMCPrimRec[i],hdNdzMCPrimGen[i],1,1,"B"); // binominal errors
  } 
  
  // write 

  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write efficiency to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

//   if(strdir && strdir.Length()){ 
//     TDirectory* dir = out.mkdir(strdir);
//     dir->cd(); 
//   }

  for(Int_t i=0; i<fNJetPtSlices; i++){

    hdNdptTracksMCPrimGen[i]->Write();
    hdNdxiMCPrimGen[i]->Write();
    hdNdzMCPrimGen[i]->Write();
    
    hdNdptTracksMCPrimRec[i]->Write();
    hdNdxiMCPrimRec[i]->Write();
    hdNdzMCPrimRec[i]->Write();
  
    hEffPt[i]->Write();
    hEffXi[i]->Write();
    hEffZ[i]->Write();
  }

  out.Close();

 
  delete fh1FFJetPtRecEffGen; 

  delete fh2FFTrackPtRecEffGen;
  delete fh2FFZRecEffGen;
  delete fh2FFXiRecEffGen;
  
  delete fh2FFTrackPtRecEffRec;
  delete fh2FFZRecEffRec;
  delete fh2FFXiRecEffRec;
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strID, TString strOutfile, 
							  Bool_t updateOutfile, TString strOutDir)
{ 
  // read task ouput from MC and write secondary correction - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kFALSE,"");
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteBgrJetSecCorr(TString strInfile, TString strBgrID, TString strID, TString strOutfile, 
							     Bool_t updateOutfile, TString strOutDir,Double_t scaleFacBgrRec)
{ 
  // read task ouput from MC and write secondary correction - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteJetSecCorr(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir,kTRUE,strBgrID,scaleFacBgrRec);
}


//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetSecCorr(TString strInfile, TString strdir, TString strlist, 
							  TString strOutfile, Bool_t updateOutfile, TString strOutDir,
							  Bool_t writeBgr, TString strBgrID,Double_t scaleFacBgrRec)
{
  // read task output from MC and write secondary correction in jet pt slices as function of pt, z, xi
  // argument strlist optional - read from directory strdir if not specified
  // write eff to file strOutfile - by default only 'update' file (don't overwrite)

  TH1F* hSecCorrPt[fNJetPtSlices];
  TH1F* hSecCorrXi[fNJetPtSlices];
  TH1F* hSecCorrZ[fNJetPtSlices];

  // corr factors using naive Pythia strangeness - for sys uncertainties
  TH1F* hSecCorrPt_nonSc[fNJetPtSlices];
  TH1F* hSecCorrXi_nonSc[fNJetPtSlices];
  TH1F* hSecCorrZ_nonSc[fNJetPtSlices];

  TH1F* hdNdptTracksMCPrimRec[fNJetPtSlices];
  TH1F* hdNdxiMCPrimRec[fNJetPtSlices];
  TH1F* hdNdzMCPrimRec[fNJetPtSlices];
    
  TH1F* hdNdptTracksMCSecRecNS[fNJetPtSlices];
  TH1F* hdNdxiMCSecRecNS[fNJetPtSlices];
  TH1F* hdNdzMCSecRecNS[fNJetPtSlices];

  TH1F* hdNdptTracksMCSecRecS[fNJetPtSlices];
  TH1F* hdNdxiMCSecRecS[fNJetPtSlices];
  TH1F* hdNdzMCSecRecS[fNJetPtSlices];

  TH1F* hdNdptTracksMCSecRecSsc[fNJetPtSlices];
  TH1F* hdNdxiMCSecRecSsc[fNJetPtSlices];
  TH1F* hdNdzMCSecRecSsc[fNJetPtSlices];
  
  // ---

  TH1F* fh1FFJetPtRecEffRec;

  TH2F* fh2FFTrackPtRecEffRec;
  TH2F* fh2FFZRecEffRec;
  TH2F* fh2FFXiRecEffRec;
  
  TH2F* fh2FFTrackPtSecRecNS;
  TH2F* fh2FFZSecRecNS;
  TH2F* fh2FFXiSecRecNS;

  TH2F* fh2FFTrackPtSecRecS;
  TH2F* fh2FFZSecRecS;
  TH2F* fh2FFXiSecRecS;

  TH2F* fh2FFTrackPtSecRecSsc;
  TH2F* fh2FFZSecRecSsc;
  TH2F* fh2FFXiSecRecSsc;
  
  // ---

  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeJetTrackSecCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){

    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }

  if(list){
    fh1FFJetPtRecEffRec    = (TH1F*) list->FindObject("fh1FFJetPtRecEffRec");

    if(writeBgr){ 
      fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
      fh2FFZRecEffRec        = (TH2F*) list->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
      fh2FFXiRecEffRec       = (TH2F*) list->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data())); 

      fh2FFTrackPtSecRecNS   = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
      fh2FFZSecRecNS         = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
      fh2FFXiSecRecNS        = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));

      fh2FFTrackPtSecRecS    = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
      fh2FFZSecRecS          = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
      fh2FFXiSecRecS         = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
      
      fh2FFTrackPtSecRecSsc  = (TH2F*) list->FindObject(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
      fh2FFZSecRecSsc        = (TH2F*) list->FindObject(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
      fh2FFXiSecRecSsc       = (TH2F*) list->FindObject(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
    }
    else{
      fh2FFTrackPtRecEffRec  = (TH2F*) list->FindObject("fh2FFTrackPtRecEffRec");
      fh2FFZRecEffRec        = (TH2F*) list->FindObject("fh2FFZRecEffRec");
      fh2FFXiRecEffRec       = (TH2F*) list->FindObject("fh2FFXiRecEffRec");
      
      fh2FFTrackPtSecRecNS   = (TH2F*) list->FindObject("fh2FFTrackPtSecRecNS");
      fh2FFZSecRecNS         = (TH2F*) list->FindObject("fh2FFZSecRecNS");
      fh2FFXiSecRecNS        = (TH2F*) list->FindObject("fh2FFXiSecRecNS");
      
      fh2FFTrackPtSecRecS    = (TH2F*) list->FindObject("fh2FFTrackPtSecRecS");
      fh2FFZSecRecS          = (TH2F*) list->FindObject("fh2FFZSecRecS");
      fh2FFXiSecRecS         = (TH2F*) list->FindObject("fh2FFXiSecRecS");
      
      fh2FFTrackPtSecRecSsc  = (TH2F*) list->FindObject("fh2FFTrackPtSecRecSsc");
      fh2FFZSecRecSsc        = (TH2F*) list->FindObject("fh2FFZSecRecSsc");
      fh2FFXiSecRecSsc       = (TH2F*) list->FindObject("fh2FFXiSecRecSsc");
    }
  }
  else{
    fh1FFJetPtRecEffRec    = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");

    if(writeBgr){
      fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
      fh2FFZRecEffRec        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
      fh2FFXiRecEffRec       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data())); 

      fh2FFTrackPtSecRecNS   = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecNS",strBgrID.Data()));
      fh2FFZSecRecNS         = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecNS",strBgrID.Data()));
      fh2FFXiSecRecNS        = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecNS",strBgrID.Data()));

      fh2FFTrackPtSecRecS    = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecS",strBgrID.Data()));
      fh2FFZSecRecS          = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecS",strBgrID.Data()));
      fh2FFXiSecRecS         = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecS",strBgrID.Data()));
      
      fh2FFTrackPtSecRecSsc  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sSecRecSsc",strBgrID.Data()));
      fh2FFZSecRecSsc        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sSecRecSsc",strBgrID.Data()));
      fh2FFXiSecRecSsc       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sSecRecSsc",strBgrID.Data()));
    }
    else{
      fh2FFTrackPtRecEffRec  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
      fh2FFZRecEffRec        = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
      fh2FFXiRecEffRec       = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec");
    
      fh2FFTrackPtSecRecNS   = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecNS");
      fh2FFZSecRecNS         = (TH2F*) gDirectory->Get("fh2FFZSecRecNS");
      fh2FFXiSecRecNS        = (TH2F*) gDirectory->Get("fh2FFXiSecRecNS");
      
      fh2FFTrackPtSecRecS    = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecS");
      fh2FFZSecRecS          = (TH2F*) gDirectory->Get("fh2FFZSecRecS");
      fh2FFXiSecRecS         = (TH2F*) gDirectory->Get("fh2FFXiSecRecS");
      
      fh2FFTrackPtSecRecSsc  = (TH2F*) gDirectory->Get("fh2FFTrackPtSecRecSsc");
      fh2FFZSecRecSsc        = (TH2F*) gDirectory->Get("fh2FFZSecRecSsc");
      fh2FFXiSecRecSsc       = (TH2F*) gDirectory->Get("fh2FFXiSecRecSsc");
    }
  }
  
  fh1FFJetPtRecEffRec->SetDirectory(0); 
   
  fh2FFTrackPtRecEffRec->SetDirectory(0);
  fh2FFZRecEffRec->SetDirectory(0);
  fh2FFXiRecEffRec->SetDirectory(0);
  
  fh2FFTrackPtSecRecNS->SetDirectory(0);
  fh2FFZSecRecNS->SetDirectory(0);
  fh2FFXiSecRecNS->SetDirectory(0);
  
  fh2FFTrackPtSecRecS->SetDirectory(0);
  fh2FFZSecRecS->SetDirectory(0);
  fh2FFXiSecRecS->SetDirectory(0);
  
  fh2FFTrackPtSecRecSsc->SetDirectory(0);
  fh2FFZSecRecSsc->SetDirectory(0);
  fh2FFXiSecRecSsc->SetDirectory(0);

  f.Close();


  // projections: FF for generated and reconstructed primaries 
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);

    Int_t binLo = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh2FFTrackPtRecEffRec->GetXaxis()->FindBin(jetPtUpLim))-1;

    if(binUp > fh2FFTrackPtRecEffRec->GetNbinsX()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameFFPtPrimRec(Form("fh1FFTrackPtRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZPrimRec(Form("fh1FFZRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiPrimRec(Form("fh1FFXiRecPrim_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    TString strNameFFPtSecRecNS(Form("fh1FFTrackPtRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZSecRecNS(Form("fh1FFZRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiSecRecNS(Form("fh1FFXiRecSecNS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
     
    TString strNameFFPtSecRecS(Form("fh1FFTrackPtRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZSecRecS(Form("fh1FFZRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiSecRecS(Form("fh1FFXiRecSecS_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
     
    TString strNameFFPtSecRecSsc(Form("fh1FFTrackPtRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZSecRecSsc(Form("fh1FFZRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiSecRecSsc(Form("fh1FFXiRecSecSsc_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));

    if(writeBgr){
      strNameFFPtPrimRec.Form("fh1BgrTrackPt%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFZPrimRec.Form("fh1BgrZ%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFXiPrimRec.Form("fh1BgrXi%sRecPrim_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      
      strNameFFPtSecRecNS.Form("fh1BgrTrackPt%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFZSecRecNS.Form("fh1BgrZ%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFXiSecRecNS.Form("fh1BgrXi%sRecSecNS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      
      strNameFFPtSecRecS.Form("fh1BgrTrackPt%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFZSecRecS.Form("fh1BgrZ%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFXiSecRecS.Form("fh1BgrXi%sRecSecS_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      
      strNameFFPtSecRecSsc.Form("fh1BgrTrackPt%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFZSecRecSsc.Form("fh1BgrZ%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
      strNameFFXiSecRecSsc.Form("fh1BgrXi%sRecSecSsc_%02d_%02d",strBgrID.Data(),static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim));
    }
    
    // project 
    // appendix 'unbinned' to avoid histos with same name after rebinning

    hdNdptTracksMCPrimRec[i] = (TH1F*) fh2FFTrackPtRecEffRec->ProjectionY(strNameFFPtPrimRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCPrimRec[i]        = (TH1F*) fh2FFZRecEffRec->ProjectionY(strNameFFZPrimRec+"_unBinned",binLo,binUp,"o");
    hdNdxiMCPrimRec[i]       = (TH1F*) fh2FFXiRecEffRec->ProjectionY(strNameFFXiPrimRec+"_unBinned",binLo,binUp,"o");
    
    hdNdptTracksMCSecRecNS[i]  = (TH1F*) fh2FFTrackPtSecRecNS->ProjectionY(strNameFFPtSecRecNS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCSecRecNS[i]         = (TH1F*) fh2FFZSecRecNS->ProjectionY(strNameFFZSecRecNS+"_unBinned",binLo,binUp,"o");
    hdNdxiMCSecRecNS[i]        = (TH1F*) fh2FFXiSecRecNS->ProjectionY(strNameFFXiSecRecNS+"_unBinned",binLo,binUp,"o");
    
    hdNdptTracksMCSecRecS[i]  = (TH1F*) fh2FFTrackPtSecRecS->ProjectionY(strNameFFPtSecRecS+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCSecRecS[i]         = (TH1F*) fh2FFZSecRecS->ProjectionY(strNameFFZSecRecS+"_unBinned",binLo,binUp,"o");
    hdNdxiMCSecRecS[i]        = (TH1F*) fh2FFXiSecRecS->ProjectionY(strNameFFXiSecRecS+"_unBinned",binLo,binUp,"o");

    hdNdptTracksMCSecRecSsc[i]  = (TH1F*) fh2FFTrackPtSecRecSsc->ProjectionY(strNameFFPtSecRecSsc+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCSecRecSsc[i]         = (TH1F*) fh2FFZSecRecSsc->ProjectionY(strNameFFZSecRecSsc+"_unBinned",binLo,binUp,"o");
    hdNdxiMCSecRecSsc[i]        = (TH1F*) fh2FFXiSecRecSsc->ProjectionY(strNameFFXiSecRecSsc+"_unBinned",binLo,binUp,"o");


    // rebin
    if(fNHistoBinsPt[i]) hdNdptTracksMCPrimRec[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtPrimRec,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCPrimRec[i]        = (TH1F*) hdNdzMCPrimRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZPrimRec,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCPrimRec[i]       = (TH1F*) hdNdxiMCPrimRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiPrimRec,fHistoBinsXi[i]->GetArray());

    if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecNS[i] = (TH1F*) hdNdptTracksMCSecRecNS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecNS,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCSecRecNS[i]        = (TH1F*) hdNdzMCSecRecNS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecNS,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCSecRecNS[i]       = (TH1F*) hdNdxiMCSecRecNS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecNS,fHistoBinsXi[i]->GetArray());

    if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecS[i] = (TH1F*) hdNdptTracksMCSecRecS[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecS,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCSecRecS[i]        = (TH1F*) hdNdzMCSecRecS[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecS,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCSecRecS[i]       = (TH1F*) hdNdxiMCSecRecS[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecS,fHistoBinsXi[i]->GetArray());

    if(fNHistoBinsPt[i]) hdNdptTracksMCSecRecSsc[i] = (TH1F*) hdNdptTracksMCSecRecSsc[i]->Rebin(fNHistoBinsPt[i],strNameFFPtSecRecSsc,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCSecRecSsc[i]        = (TH1F*) hdNdzMCSecRecSsc[i]->Rebin(fNHistoBinsZ[i],strNameFFZSecRecSsc,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCSecRecSsc[i]       = (TH1F*) hdNdxiMCSecRecSsc[i]->Rebin(fNHistoBinsXi[i],strNameFFXiSecRecSsc,fHistoBinsXi[i]->GetArray());


    hdNdptTracksMCPrimRec[i]->SetNameTitle(strNameFFPtPrimRec,"");
    hdNdzMCPrimRec[i]->SetNameTitle(strNameFFZPrimRec,"");
    hdNdxiMCPrimRec[i]->SetNameTitle(strNameFFXiPrimRec,"");
    
    hdNdptTracksMCSecRecNS[i]->SetNameTitle(strNameFFPtSecRecNS,"");
    hdNdzMCSecRecNS[i]->SetNameTitle(strNameFFZSecRecNS,"");
    hdNdxiMCSecRecNS[i]->SetNameTitle(strNameFFXiSecRecNS,"");
    
    hdNdptTracksMCSecRecS[i]->SetNameTitle(strNameFFPtSecRecS,"");
    hdNdzMCSecRecS[i]->SetNameTitle(strNameFFZSecRecS,"");
    hdNdxiMCSecRecS[i]->SetNameTitle(strNameFFXiSecRecS,"");

    hdNdptTracksMCSecRecSsc[i]->SetNameTitle(strNameFFPtSecRecSsc,"");
    hdNdzMCSecRecSsc[i]->SetNameTitle(strNameFFZSecRecSsc,"");
    hdNdxiMCSecRecSsc[i]->SetNameTitle(strNameFFXiSecRecSsc,"");
    
    // normalize 
    Double_t nJetsBin = fh1FFJetPtRecEffRec->Integral(binLo,binUp);

    // scale fac for perp2 bgr
    if(writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBin /= scaleFacBgrRec;

    NormalizeTH1(hdNdptTracksMCPrimRec[i],nJetsBin); 
    NormalizeTH1(hdNdzMCPrimRec[i],nJetsBin); 
    NormalizeTH1(hdNdxiMCPrimRec[i],nJetsBin); 

    NormalizeTH1(hdNdptTracksMCSecRecNS[i],nJetsBin); 
    NormalizeTH1(hdNdzMCSecRecNS[i],nJetsBin); 
    NormalizeTH1(hdNdxiMCSecRecNS[i],nJetsBin); 

    NormalizeTH1(hdNdptTracksMCSecRecS[i],nJetsBin); 
    NormalizeTH1(hdNdzMCSecRecS[i],nJetsBin); 
    NormalizeTH1(hdNdxiMCSecRecS[i],nJetsBin); 

    NormalizeTH1(hdNdptTracksMCSecRecSsc[i],nJetsBin); 
    NormalizeTH1(hdNdzMCSecRecSsc[i],nJetsBin); 
    NormalizeTH1(hdNdxiMCSecRecSsc[i],nJetsBin); 

    
    // divide prim / (prim+sec) : corr factor
    TString strNameSecCorrPt(Form("hSecCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameSecCorrZ(Form("hSecCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameSecCorrXi(Form("hSecCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    
    if(writeBgr){
      strNameSecCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameSecCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameSecCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
    }
     
    hSecCorrPt[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt);
    TH1F* hSumPrimSecPt = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt");
    hSumPrimSecPt->Add(hdNdptTracksMCSecRecNS[i]);
    hSumPrimSecPt->Add(hdNdptTracksMCSecRecSsc[i]);
    hSecCorrPt[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt,1,1,"B"); // binominal errors

    hSecCorrXi[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi);
    TH1F* hSumPrimSecXi = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi");
    hSumPrimSecXi->Add(hdNdxiMCSecRecNS[i]);
    hSumPrimSecXi->Add(hdNdxiMCSecRecSsc[i]);
    hSecCorrXi[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi,1,1,"B"); // binominal errors

    hSecCorrZ[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ);
    TH1F* hSumPrimSecZ = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ");
    hSumPrimSecZ->Add(hdNdzMCSecRecNS[i]);
    hSumPrimSecZ->Add(hdNdzMCSecRecSsc[i]);
    hSecCorrZ[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ,1,1,"B"); // binominal errors

    // the same using unscaled strangeness
    TString strNameSecCorrPt_nonSc(Form("hSecCorrPt_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameSecCorrZ_nonSc(Form("hSecCorrZ_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameSecCorrXi_nonSc(Form("hSecCorrXi_nonSc_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));

    if(writeBgr){
      strNameSecCorrPt_nonSc.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameSecCorrZ_nonSc.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameSecCorrXi_nonSc.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
    }

    hSecCorrPt_nonSc[i] = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone(strNameSecCorrPt_nonSc);
    TH1F* hSumPrimSecPt_nonSc = (TH1F*) hdNdptTracksMCPrimRec[i]->Clone("hSumPrimSecPt_nonSc");
    hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecNS[i]);
    hSumPrimSecPt_nonSc->Add(hdNdptTracksMCSecRecS[i]); // non-scaled secondaries from strangeness
    hSecCorrPt_nonSc[i]->Divide(hdNdptTracksMCPrimRec[i],hSumPrimSecPt_nonSc,1,1,"B"); // binominal errors
    
    hSecCorrZ_nonSc[i] = (TH1F*) hdNdzMCPrimRec[i]->Clone(strNameSecCorrZ_nonSc);
    TH1F* hSumPrimSecZ_nonSc = (TH1F*) hdNdzMCPrimRec[i]->Clone("hSumPrimSecZ_nonSc");
    hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecNS[i]);
    hSumPrimSecZ_nonSc->Add(hdNdzMCSecRecS[i]); // non-scaled secondaries from strangeness
    hSecCorrZ_nonSc[i]->Divide(hdNdzMCPrimRec[i],hSumPrimSecZ_nonSc,1,1,"B"); // binominal errors
    
    hSecCorrXi_nonSc[i] = (TH1F*) hdNdxiMCPrimRec[i]->Clone(strNameSecCorrXi_nonSc);
    TH1F* hSumPrimSecXi_nonSc = (TH1F*) hdNdxiMCPrimRec[i]->Clone("hSumPrimSecXi_nonSc");
    hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecNS[i]);
    hSumPrimSecXi_nonSc->Add(hdNdxiMCSecRecS[i]); // non-scaled secondaries from strangeness
    hSecCorrXi_nonSc[i]->Divide(hdNdxiMCPrimRec[i],hSumPrimSecXi_nonSc,1,1,"B"); // binominal errors
  } 
  
  // write 

  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening sec corr output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- write jet sec corr to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

  if(strOutDir && strOutDir.Length()){  

    TDirectory* dir;
    if((dir = ((TDirectory*)  gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }
  

  for(Int_t i=0; i<fNJetPtSlices; i++){
  
    hSecCorrPt[i]->Write();
    hSecCorrXi[i]->Write();
    hSecCorrZ[i]->Write();
  
    hSecCorrPt_nonSc[i]->Write();
    hSecCorrXi_nonSc[i]->Write();
    hSecCorrZ_nonSc[i]->Write();
  }

  TString strSpectraDir      = "spectraSecCorr";
  if(writeBgr) strSpectraDir = "spectraBgrSecCorr";
  
  TDirectory *dOut  = gDirectory->mkdir(strSpectraDir);
  dOut->cd();
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    hdNdptTracksMCPrimRec[i]->Write(); 
    hdNdzMCPrimRec[i]->Write(); 
    hdNdxiMCPrimRec[i]->Write(); 
    
    hdNdptTracksMCSecRecNS[i]->Write(); 
    hdNdzMCSecRecNS[i]->Write(); 
    hdNdxiMCSecRecNS[i]->Write(); 
    
    hdNdptTracksMCSecRecS[i]->Write(); 
    hdNdzMCSecRecS[i]->Write(); 
    hdNdxiMCSecRecS[i]->Write(); 
    
    hdNdptTracksMCSecRecSsc[i]->Write(); 
    hdNdzMCSecRecSsc[i]->Write(); 
    hdNdxiMCSecRecSsc[i]->Write();  
  }

  out.Close();
  
  delete fh1FFJetPtRecEffRec; 

  delete fh2FFTrackPtRecEffRec;
  delete fh2FFZRecEffRec;
  delete fh2FFXiRecEffRec;
  
  delete fh2FFTrackPtSecRecNS;
  delete fh2FFZSecRecNS;
  delete fh2FFXiSecRecNS;
  
  delete fh2FFTrackPtSecRecS;
  delete fh2FFZSecRecS;
  delete fh2FFXiSecRecS;
  
  delete fh2FFTrackPtSecRecSsc;
  delete fh2FFZSecRecSsc;
  delete fh2FFXiSecRecSsc;
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strID, TString strOutfile, 
							   Bool_t updateOutfile, TString strOutDir )
{ 
  // read task ouput from MC and write single track eff - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  WriteJetResponse(strInfile,strdir,strlist,strOutfile,updateOutfile, strOutDir);
}

//_____________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetResponse(TString strInfile, TString strdir, TString strlist,
							   TString strOutfile, Bool_t updateOutfile, TString strOutDir)
{
  // read 3d THnSparse response matrices in pt,z,xi vs jet pt from file
  // project THnSparse + TH2 in jet pt slices 
  // write to strOutfile 

  THnSparse* hn3ResponseJetPt;
  THnSparse* hn3ResponseJetZ;
  THnSparse* hn3ResponseJetXi;

  // 2D response matrices 

  THnSparse* hnResponsePt[fNJetPtSlices];
  THnSparse* hnResponseZ[fNJetPtSlices];
  THnSparse* hnResponseXi[fNJetPtSlices];

  TH2F* h2ResponsePt[fNJetPtSlices];
  TH2F* h2ResponseZ[fNJetPtSlices];
  TH2F* h2ResponseXi[fNJetPtSlices];

  // 1D projections on gen pt / rec pt axes

  TH1F* h1FFPtRec[fNJetPtSlices]; 
  TH1F* h1FFZRec[fNJetPtSlices];
  TH1F* h1FFXiRec[fNJetPtSlices];

  TH1F* h1FFPtGen[fNJetPtSlices]; 
  TH1F* h1FFZGen[fNJetPtSlices];
  TH1F* h1FFXiGen[fNJetPtSlices];

  TH1F* h1RatioPt[fNJetPtSlices]; 
  TH1F* h1RatioZ[fNJetPtSlices];
  TH1F* h1RatioXi[fNJetPtSlices];



  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
    
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  
  
  if(list){
    hn3ResponseJetPt = (THnSparse*) list->FindObject("fhnResponseJetTrackPt");
    hn3ResponseJetZ  = (THnSparse*) list->FindObject("fhnResponseJetZ");
    hn3ResponseJetXi = (THnSparse*) list->FindObject("fhnResponseJetXi");
  }
  else{
    hn3ResponseJetPt = (THnSparse*) gDirectory->Get("fhnResponseJetTrackPt");
    hn3ResponseJetZ  = (THnSparse*) gDirectory->Get("fhnResponseJetZ");
    hn3ResponseJetXi = (THnSparse*) gDirectory->Get("fhnResponseJetXi");
  }

  
  if(!hn3ResponseJetPt){
    Printf("%s:%d -- error retrieving response matrix fhnResponseJetTrackPt",(char*)__FILE__,__LINE__);
    return;
  }

  if(!hn3ResponseJetZ){
    Printf("%s:%d -- error retrieving response matrix fhnResponseJetZ",(char*)__FILE__,__LINE__);
    return;
  }

  if(!hn3ResponseJetXi){
    Printf("%s:%d -- error retrieving response matrix fhnResponseJetXi",(char*)__FILE__,__LINE__);
    return;
  }

  f.Close();  

  // axes 

  Int_t axisJetPtTHn3 = -1;
  Int_t axisGenPtTHn3 = -1;
  Int_t axisRecPtTHn3 = -1;

  for(Int_t i=0; i<hn3ResponseJetPt->GetNdimensions(); i++){
    
    TString title = hn3ResponseJetPt->GetAxis(i)->GetTitle(); 

    if(title.Contains("jet p_{T}")) axisJetPtTHn3 = i; 
    if(title.Contains("gen p_{T}")) axisGenPtTHn3 = i; 
    if(title.Contains("rec p_{T}")) axisRecPtTHn3 = i; 
  }

  if(axisJetPtTHn3 == -1){
    Printf("%s:%d -- error axisJetPtTHn3",(char*)__FILE__,__LINE__);
    return;
  }

  if(axisGenPtTHn3 == -1){
    Printf("%s:%d -- error axisGenPtTHn3",(char*)__FILE__,__LINE__);
    return;
  }

  if(axisRecPtTHn3 == -1){
    Printf("%s:%d -- error axisRecPtTHn3",(char*)__FILE__,__LINE__);
    return;
  }

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);

    Int_t binLo = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->FindBin(jetPtUpLim))-1;

    if(binUp > hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->GetNbins()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameRespPt(Form("hnResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameRespZ(Form("hnResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameRespXi(Form("hnResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));

    TString strNameTH2RespPt(Form("h2ResponsePt_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameTH2RespZ(Form("h2ResponseZ_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameTH2RespXi(Form("h2ResponseXi_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
         
    TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
 
    TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
 
    TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
 
    
    // 2D projections in jet pt range
  
    hn3ResponseJetPt->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
    hn3ResponseJetZ->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 
    hn3ResponseJetXi->GetAxis(axisJetPtTHn3)->SetRange(binLo,binUp); 

    Int_t axesProj[2] = {axisRecPtTHn3,axisGenPtTHn3};  

    hnResponsePt[i] = hn3ResponseJetPt->Projection(2,axesProj);
    hnResponseZ[i]  = hn3ResponseJetZ->Projection(2,axesProj);
    hnResponseXi[i] = hn3ResponseJetXi->Projection(2,axesProj);
  
    hnResponsePt[i]->SetNameTitle(strNameRespPt,""); 
    hnResponseZ[i]->SetNameTitle(strNameRespZ,""); 
    hnResponseXi[i]->SetNameTitle(strNameRespXi,""); 

    h2ResponsePt[i] = (TH2F*) hnResponsePt[i]->Projection(1,0);// note convention: yDim,xDim
    h2ResponseZ[i]  = (TH2F*) hnResponseZ[i]->Projection(1,0); // note convention: yDim,xDim
    h2ResponseXi[i] = (TH2F*) hnResponseXi[i]->Projection(1,0);// note convention: yDim,xDim
 
    h2ResponsePt[i]->SetNameTitle(strNameTH2RespPt,""); 
    h2ResponseZ[i]->SetNameTitle(strNameTH2RespZ,""); 
    h2ResponseXi[i]->SetNameTitle(strNameTH2RespXi,""); 


    // 1D projections

    Int_t axisGenPtTHn2 = -1;
    Int_t axisRecPtTHn2 = -1;

    for(Int_t d=0; d<hnResponsePt[i]->GetNdimensions(); d++){
    
      TString title = hnResponsePt[i]->GetAxis(d)->GetTitle(); 

      if(title.Contains("gen p_{T}")) axisGenPtTHn2 = d; 
      if(title.Contains("rec p_{T}")) axisRecPtTHn2 = d; 
    }

    
    if(axisGenPtTHn2 == -1){
      Printf("%s:%d -- error axisGenPtTHn2",(char*)__FILE__,__LINE__);
      return;
    }
    
    if(axisRecPtTHn2 == -1){
      Printf("%s:%d -- error axisRecPtTHn2",(char*)__FILE__,__LINE__);
      return;
    }
    

    h1FFPtRec[i] = (TH1F*) hnResponsePt[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
    h1FFZRec[i]  = (TH1F*) hnResponseZ[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
    h1FFXiRec[i] = (TH1F*) hnResponseXi[i]->Projection(axisRecPtTHn2);// note convention: yDim,xDim
    
    h1FFPtRec[i]->SetNameTitle(strNameRecPt,""); 
    h1FFZRec[i]->SetNameTitle(strNameRecZ,""); 
    h1FFXiRec[i]->SetNameTitle(strNameRecXi,""); 

    h1FFPtGen[i] = (TH1F*) hnResponsePt[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
    h1FFZGen[i]  = (TH1F*) hnResponseZ[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
    h1FFXiGen[i] = (TH1F*) hnResponseXi[i]->Projection(axisGenPtTHn2);// note convention: yDim,xDim
    
    h1FFPtGen[i]->SetNameTitle(strNameGenPt,"");
    h1FFZGen[i]->SetNameTitle(strNameGenZ,"");
    h1FFXiGen[i]->SetNameTitle(strNameGenXi,"");

    // normalize 1D projections

    if(fNHistoBinsPt[i]) h1FFPtRec[i] = (TH1F*) h1FFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
    if(fNHistoBinsZ[i])  h1FFZRec[i]  = (TH1F*) h1FFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) h1FFXiRec[i] = (TH1F*) h1FFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
    
    if(fNHistoBinsPt[i]) h1FFPtGen[i] = (TH1F*) h1FFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  h1FFZGen[i]  = (TH1F*) h1FFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) h1FFXiGen[i] = (TH1F*) h1FFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
       
    NormalizeTH1(h1FFPtRec[i],fNJets->At(i)); 
    NormalizeTH1(h1FFZRec[i],fNJets->At(i));
    NormalizeTH1(h1FFXiRec[i],fNJets->At(i));

    NormalizeTH1(h1FFPtGen[i],fNJets->At(i)); 
    NormalizeTH1(h1FFZGen[i],fNJets->At(i));
    NormalizeTH1(h1FFXiGen[i],fNJets->At(i));

    // ratios 1D projections

    h1RatioPt[i] = (TH1F*) h1FFPtRec[i]->Clone(strNameRatioPt);
    h1RatioPt[i]->Reset();
    h1RatioPt[i]->Divide(h1FFPtRec[i],h1FFPtGen[i],1,1,"B");

    h1RatioZ[i] = (TH1F*) h1FFZRec[i]->Clone(strNameRatioZ);
    h1RatioZ[i]->Reset();
    h1RatioZ[i]->Divide(h1FFZRec[i],h1FFZGen[i],1,1,"B");
    
    h1RatioXi[i] = (TH1F*) h1FFXiRec[i]->Clone(strNameRatioXi);
    h1RatioXi[i]->Reset();
    h1RatioXi[i]->Divide(h1FFXiRec[i],h1FFXiGen[i],1,1,"B");
  }
   
  
  // write 

  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";
  
  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write response matrices to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

  if(strOutDir && strOutDir.Length()){  
    
    TDirectory* dir;
    if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }


  for(Int_t i=0; i<fNJetPtSlices; i++){

    hnResponsePt[i]->Write();
    hnResponseXi[i]->Write();
    hnResponseZ[i]->Write();

    h2ResponsePt[i]->Write();
    h2ResponseXi[i]->Write();
    h2ResponseZ[i]->Write();

    h1FFPtRec[i]->Write(); 
    h1FFZRec[i]->Write();
    h1FFXiRec[i]->Write();
    
    h1FFPtGen[i]->Write(); 
    h1FFZGen[i]->Write();
    h1FFXiGen[i]->Write();
    
    h1RatioPt[i]->Write(); 
    h1RatioZ[i]->Write();
    h1RatioXi[i]->Write();

  }

  out.Close();  
}

//______________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadResponse(TString strfile, TString strdir, TString strlist)
{
  // read response matrices from file
  // argument strlist optional - read from directory strdir if not specified
  // note: THnSparse are not rebinned 

  THnSparse* hRespPt[fNJetPtSlices]; 
  THnSparse* hRespZ[fNJetPtSlices];
  THnSparse* hRespXi[fNJetPtSlices];
  
  for(Int_t i=0; i<fNJetPtSlices; i++) hRespPt[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hRespZ[i]  = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hRespXi[i] = 0;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read response matrices from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
   
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
    
    TString strNameRespPt(Form("hnResponsePt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameRespZ(Form("hnResponseZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameRespXi(Form("hnResponseXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
        
    if(list){
      hRespPt[i] = (THnSparse*) list->FindObject(strNameRespPt); 
      hRespZ[i]  = (THnSparse*) list->FindObject(strNameRespZ); 
      hRespXi[i] = (THnSparse*) list->FindObject(strNameRespXi); 
    }
    else{
      hRespPt[i] = (THnSparse*) gDirectory->Get(strNameRespPt); 
      hRespZ[i]  = (THnSparse*) gDirectory->Get(strNameRespZ); 
      hRespXi[i] = (THnSparse*) gDirectory->Get(strNameRespXi); 
    }
    
    if(!hRespPt[i]){
      Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespPt.Data());
    }
  
    if(!hRespZ[i]){
      Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespZ.Data());
    }    

    if(!hRespXi[i]){
      Printf("%s:%d -- error retrieving response %s", (char*)__FILE__,__LINE__,strNameRespXi.Data());
    }
   
    //    if(0){ // can't rebin THnSparse ...
    //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(0,fHistoBinsPt[i]->GetArray());
    //       if(fNHistoBinsPt[i]) hRespPt[i]->SetBinEdges(1,fHistoBinsPt[i]->GetArray());
    
    //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(0,fHistoBinsZ[i]->GetArray());
    //       if(fNHistoBinsZ[i])  hRespZ[i]->SetBinEdges(1,fHistoBinsZ[i]->GetArray());
    
    //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(0,fHistoBinsXi[i]->GetArray());
    //       if(fNHistoBinsXi[i]) hRespXi[i]->SetBinEdges(1,fHistoBinsXi[i]->GetArray());
    //     }
   
 
  } // jet slices loop

  f.Close();  // THnSparse pointers still valid even if file closed

//   for(Int_t i=0; i<fNJetPtSlices; i++){  // no copy c'tor ...
//     if(hRespPt[i]) new(fhnResponsePt[i]) THnSparseF(*hRespPt[i]);
//     if(hRespZ[i])  new(fhnResponseZ[i])  THnSparseF(*hRespZ[i]);
//     if(hRespXi[i]) new(fhnResponseXi[i]) THnSparseF(*hRespXi[i]);
//   }

  for(Int_t i=0; i<fNJetPtSlices; i++){ 
    fhnResponsePt[i] = hRespPt[i];
    fhnResponseZ[i]  = hRespZ[i];
    fhnResponseXi[i] = hRespXi[i];
  }
}

//______________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadPriors(TString strfile,const Int_t type)
{
  // read priors from file: rec primaries, gen pt dist

  if(fDebug>0) Printf("%s:%d -- read priors from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
  
  // temporary histos to store pointers from file
  TH1F* hist[fNJetPtSlices]; 

  for(Int_t i=0; i<fNJetPtSlices; i++) hist[i] = 0;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
    
    TString strName;
    
    if(type == kFlagPt) strName.Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim); 
    if(type == kFlagZ)  strName.Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);       
    if(type == kFlagXi) strName.Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim);      
    
    hist[i] = (TH1F*) gDirectory->Get(strName); 
    
    if(!hist[i]){
      Printf("%s:%d -- error retrieving prior %s", (char*)__FILE__,__LINE__,strName.Data());
    }
    

    //if(fNHistoBinsPt[i]) hist[i] = (TH1F*) hist[i]->Rebin(fNHistoBinsPt[i],hist[i]->GetName()+"_rebin",fHistoBinsPt[i]->GetArray());

    if(hist[i]) hist[i]->SetDirectory(0); 
    
  } // jet slices loop
  
  f.Close();

  
  for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
    if(hist[i] && type == kFlagPt) new(fh1FFTrackPtPrior[i]) TH1F(*hist[i]); 
    if(hist[i] && type == kFlagZ)  new(fh1FFZPrior[i])       TH1F(*hist[i]); 
    if(hist[i] && type == kFlagXi) new(fh1FFXiPrior[i])      TH1F(*hist[i]); 
  }
}


//_____________________________________________________
// void AliFragmentationFunctionCorrections::RatioRecGen()
// {
//   // create ratio reconstructed over generated FF
//   // use current highest corrLevel 
  
//   Printf("%s:%d -- build ratio rec.gen, corrLevel %d",(char*)__FILE__,__LINE__,fNCorrectionLevels-1);

//   for(Int_t i=0; i<fNJetPtSlices; i++){
 
//     TH1F* histPtRec =  fCorrFF[fNCorrectionLevels-1]->GetTrackPt(i); // levels -1: latest corr level
//     TH1F* histZRec  =  fCorrFF[fNCorrectionLevels-1]->GetZ(i);       // levels -1: latest corr level
//     TH1F* histXiRec =  fCorrFF[fNCorrectionLevels-1]->GetXi(i);      // levels -1: latest corr level
    
//     TH1F* histPtGen = fh1FFTrackPtGenPrim[i]; 
//     TH1F* histZGen  = fh1FFZGenPrim[i];
//     TH1F* histXiGen = fh1FFXiGenPrim[i];

//     TString histNamePt = histPtRec->GetName();
//     TString histNameZ  = histZRec->GetName();
//     TString histNameXi = histXiRec->GetName();

//     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecGen");
//     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecGen");
//     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecGen");

//     // ratio 
//     TH1F* hRatioRecGenPt = (TH1F*) histPtRec->Clone(histNamePt);
//     hRatioRecGenPt->Reset();
//     hRatioRecGenPt->Divide(histPtRec,histPtGen,1,1,"B");

//     TH1F* hRatioRecGenZ = (TH1F*) histZRec->Clone(histNameZ);
//     hRatioRecGenZ->Reset();
//     hRatioRecGenZ->Divide(histZRec,histZGen,1,1,"B");
    
//     TH1F* hRatioRecGenXi = (TH1F*) histXiRec->Clone(histNameXi);
//     hRatioRecGenXi->Reset();
//     hRatioRecGenXi->Divide(histXiRec,histXiGen,1,1,"B");

//     new(fh1FFRatioRecGenPt[i]) TH1F(*hRatioRecGenPt);
//     new(fh1FFRatioRecGenZ[i])  TH1F(*hRatioRecGenZ);
//     new(fh1FFRatioRecGenXi[i]) TH1F(*hRatioRecGenXi);
//   }
// }

// //___________________________________________________________
// void AliFragmentationFunctionCorrections::RatioRecPrimaries()
// {
//   // create ratio reconstructed tracks over reconstructed primaries 
//   // use raw FF (corrLevel 0)
  
//   Printf("%s:%d -- build ratio rec tracks /rec primaries",(char*)__FILE__,__LINE__);

//   for(Int_t i=0; i<fNJetPtSlices; i++){
 
//     const Int_t corrLevel = 0; 

//     TH1F* histPtRec =  fCorrFF[corrLevel]->GetTrackPt(i); // levels -1: latest corr level
//     TH1F* histZRec  =  fCorrFF[corrLevel]->GetZ(i);       // levels -1: latest corr level
//     TH1F* histXiRec =  fCorrFF[corrLevel]->GetXi(i);      // levels -1: latest corr level
    
//     TH1F* histPtRecPrim = fh1FFTrackPtRecPrim[i]; 
//     TH1F* histZRecPrim  = fh1FFZRecPrim[i];
//     TH1F* histXiRecPrim = fh1FFXiRecPrim[i];

//     TString histNamePt = histPtRec->GetName();
//     TString histNameZ  = histZRec->GetName();
//     TString histNameXi = histXiRec->GetName();

//     histNamePt.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
//     histNameZ.ReplaceAll("fh1FF","fh1FFRatioRecPrim");
//     histNameXi.ReplaceAll("fh1FF","fh1FFRatioRecPrim");

//     // ratio 
//     TH1F* hRatioRecPrimPt = (TH1F*) histPtRec->Clone(histNamePt);
//     hRatioRecPrimPt->Reset();
//     hRatioRecPrimPt->Divide(histPtRec,histPtRecPrim,1,1,"B");

//     TH1F* hRatioRecPrimZ = (TH1F*) histZRec->Clone(histNameZ);
//     hRatioRecPrimZ->Reset();
//     hRatioRecPrimZ->Divide(histZRec,histZRecPrim,1,1,"B");
    
//     TH1F* hRatioRecPrimXi = (TH1F*) histXiRec->Clone(histNameXi);
//     hRatioRecPrimXi->Reset();
//     hRatioRecPrimXi->Divide(histXiRec,histXiRecPrim,1,1,"B");


//     new(fh1FFRatioRecPrimPt[i]) TH1F(*hRatioRecPrimPt);
//     new(fh1FFRatioRecPrimZ[i])  TH1F(*hRatioRecPrimZ);
//     new(fh1FFRatioRecPrimXi[i]) TH1F(*hRatioRecPrimXi);
//   }
// }

// __________________________________________________________________________________
void AliFragmentationFunctionCorrections::ProjectJetResponseMatrices(TString strOutfile)
{

  // project response matrices on both axes: 
  // FF for rec primaries, in terms of generated and reconstructed momentum
  // write FF and ratios to outFile

  Printf("%s:%d -- project response matrices, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());

  TH1F* hFFPtRec[fNJetPtSlices]; 
  TH1F* hFFZRec[fNJetPtSlices];
  TH1F* hFFXiRec[fNJetPtSlices];

  TH1F* hFFPtGen[fNJetPtSlices]; 
  TH1F* hFFZGen[fNJetPtSlices];
  TH1F* hFFXiGen[fNJetPtSlices];

  TH1F* hRatioPt[fNJetPtSlices]; 
  TH1F* hRatioZ[fNJetPtSlices];
  TH1F* hRatioXi[fNJetPtSlices];


  Int_t axisGenPt = 1;
  Int_t axisRecPt = 0;

  for(Int_t i=0; i<fNJetPtSlices; i++){ 

    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
 
    TString strNameRecPt(Form("h1FFTrackPtRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameRecZ(Form("h1FFZRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameRecXi(Form("h1FFXiRecPrim_recPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
 
    TString strNameGenPt(Form("h1FFTrackPtRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameGenZ(Form("h1FFZRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameGenXi(Form("h1FFXiRecPrim_genPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
 
    TString strNameRatioPt(Form("h1RatioTrackPtRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameRatioZ(Form("h1RatioZRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameRatioXi(Form("h1RatioXiRecPrim_%02d_%02d",jetPtLoLim,jetPtUpLim));
 

    hFFPtRec[i] = (TH1F*) fhnResponsePt[i]->Projection(axisRecPt);// note convention: yDim,xDim
    hFFZRec[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisRecPt);// note convention: yDim,xDim
    hFFXiRec[i] = (TH1F*) fhnResponseXi[i]->Projection(axisRecPt);// note convention: yDim,xDim
    
    hFFPtRec[i]->SetNameTitle(strNameRecPt,""); 
    hFFZRec[i]->SetNameTitle(strNameRecZ,""); 
    hFFXiRec[i]->SetNameTitle(strNameRecXi,""); 


    hFFPtGen[i] = (TH1F*) fhnResponsePt[i]->Projection(axisGenPt);// note convention: yDim,xDim
    hFFZGen[i]  = (TH1F*) fhnResponseZ[i]->Projection(axisGenPt);// note convention: yDim,xDim
    hFFXiGen[i] = (TH1F*) fhnResponseXi[i]->Projection(axisGenPt);// note convention: yDim,xDim
    
    hFFPtGen[i]->SetNameTitle(strNameGenPt,""); 
    hFFZGen[i]->SetNameTitle(strNameGenZ,""); 
    hFFXiGen[i]->SetNameTitle(strNameGenXi,""); 
   

    if(fNHistoBinsPt[i]) hFFPtRec[i] = (TH1F*) hFFPtRec[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray()); 
    if(fNHistoBinsZ[i])  hFFZRec[i]   = (TH1F*) hFFZRec[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hFFXiRec[i]  = (TH1F*) hFFXiRec[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
    
    if(fNHistoBinsPt[i]) hFFPtGen[i] = (TH1F*) hFFPtGen[i]->Rebin(fNHistoBinsPt[i],"",fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hFFZGen[i]  = (TH1F*) hFFZGen[i]->Rebin(fNHistoBinsZ[i],"",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hFFXiGen[i] = (TH1F*) hFFXiGen[i]->Rebin(fNHistoBinsXi[i],"",fHistoBinsXi[i]->GetArray());
    
    NormalizeTH1(hFFPtGen[i],fNJets->At(i)); 
    NormalizeTH1(hFFZGen[i],fNJets->At(i));
    NormalizeTH1(hFFXiGen[i],fNJets->At(i));
   
    NormalizeTH1(hFFPtRec[i],fNJets->At(i)); 
    NormalizeTH1(hFFZRec[i],fNJets->At(i));
    NormalizeTH1(hFFXiRec[i],fNJets->At(i));


    hRatioPt[i] = (TH1F*) hFFPtRec[i]->Clone(strNameRatioPt);
    hRatioPt[i]->Reset();
    hRatioPt[i]->Divide(hFFPtRec[i],hFFPtGen[i],1,1,"B");

    hRatioZ[i] = (TH1F*) hFFZRec[i]->Clone(strNameRatioZ);
    hRatioZ[i]->Reset();
    hRatioZ[i]->Divide(hFFZRec[i],hFFZGen[i],1,1,"B");
    
    hRatioXi[i] = (TH1F*) hFFXiRec[i]->Clone(strNameRatioXi);
    hRatioXi[i]->Reset();
    hRatioXi[i]->Divide(hFFXiRec[i],hFFXiGen[i],1,1,"B");
  }
    


  // write 
  
  TFile out(strOutfile,"RECREATE");
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening response matrix output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    hFFPtRec[i]->Write(); 
    hFFZRec[i]->Write();
    hFFXiRec[i]->Write();
    
    hFFPtGen[i]->Write(); 
    hFFZGen[i]->Write();
    hFFXiGen[i]->Write();
    
    hRatioPt[i]->Write(); 
    hRatioZ[i]->Write();
    hRatioXi[i]->Write();
  }

  out.Close();  
}

// ____________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ProjectSingleResponseMatrix(TString strOutfile, Bool_t updateOutfile, TString strOutDir)
{ 
  // project response matrix on both axes: 
  // pt spec for rec primaries, in terms of generated and reconstructed momentum
  // write spec and ratios to outFile
  
  Printf("%s:%d -- project single pt response matrix, write to %s",(char*)__FILE__,__LINE__,strOutfile.Data());
  
  TH1F* hSpecPtRec;  
  TH1F* hSpecPtGen; 
  TH1F* hRatioPt; 
  
  Int_t axisGenPt = 1;
  Int_t axisRecPt = 0;
  
  TString strNameRecPt   = "h1SpecTrackPtRecPrim_recPt";
  TString strNameGenPt   = "h1SpecTrackPtRecPrim_genPt";
  TString strNameRatioPt = "h1RatioTrackPtRecPrim";
    
  hSpecPtRec = (TH1F*) fhnResponseSinglePt->Projection(axisRecPt);// note convention: yDim,xDim
  hSpecPtRec->SetNameTitle(strNameRecPt,""); 
  
  hSpecPtGen = (TH1F*) fhnResponseSinglePt->Projection(axisGenPt);// note convention: yDim,xDim
  hSpecPtGen->SetNameTitle(strNameGenPt,""); 
  
  hRatioPt = (TH1F*) hSpecPtRec->Clone(strNameRatioPt);
  hRatioPt->Reset();
  hRatioPt->Divide(hSpecPtRec,hSpecPtGen,1,1,"B");
  
  TString outfileOption = "RECREATE";
  if(updateOutfile) outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening reponse matrix projections output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }
      

  if(strOutDir && strOutDir.Length()){  
    
    TDirectory* dir;
    if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }
  
  hSpecPtRec->Write(); 
  hSpecPtGen->Write(); 
  hRatioPt->Write(); 
  
  out.Close();  
}


//__________________________________________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::RebinHisto(const Int_t jetPtSlice, const Int_t nBinsLimits, Double_t* binsLimits, Double_t* binsWidth, const Int_t type)
{ 
  // rebin histo, rescale bins according to new width
  // only correct for input histos with equal bin size

  // args: jetPtSlice, type, use current corr level
 
  // function takes array of limits and widths (e.g. 1 GeV bins up to 10 GeV, 2 GeV width up to 20 GeV, ...)  
  // array size of binsLimits: nBinsLimits 
  // array size of binsWidth: nBinsLimits-1 
  // binsLimits have to be in increasing order
  // if binning undefined for any slice, original binning will be kept

  if(!fNJetPtSlices){
    Printf("%s:%d -- jetPtSlices not defined", (char*)__FILE__,__LINE__);
    return;
  }

  if(jetPtSlice>=fNJetPtSlices){
    Printf("%s:%d -- jetPtSlice %d exceeds max",(char*)__FILE__,__LINE__,jetPtSlice);
    return;
  }


  Double_t binLimitMin = binsLimits[0];
  Double_t binLimitMax = binsLimits[nBinsLimits-1];

  Double_t binLimit = binLimitMin; // start value 
  
  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 ...
  TArrayD binsArray(sizeUpperLim);
  Int_t nBins = 0; 
  binsArray.SetAt(binLimitMin,nBins++);

  while(binLimit<binLimitMax && nBins<sizeUpperLim){

    Int_t currentSlice = -1;
    for(Int_t i=0; i<nBinsLimits; i++){
      if(binLimit >= 0.999*binsLimits[i]) currentSlice = i; // 0.999 numerical saftey factor 
    }
    
    Double_t currentBinWidth = binsWidth[currentSlice];
    binLimit += currentBinWidth;

    binsArray.SetAt(binLimit,nBins++);
  }
  

  TH1F* hist = 0;
  if(type == kFlagPt)            hist = fCorrFF[fNCorrectionLevels-1]->GetTrackPt(jetPtSlice); 
  else if(type == kFlagZ)        hist = fCorrFF[fNCorrectionLevels-1]->GetZ(jetPtSlice);       
  else if(type == kFlagXi)       hist = fCorrFF[fNCorrectionLevels-1]->GetXi(jetPtSlice);      
  else if(type == kFlagSinglePt) hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->GetTrackPt(0); 
  else{ 
    Printf("%s%d unknown type",(char*)__FILE__,__LINE__);
    return;
  }
  
  Double_t binWidthNoRebin = hist->GetBinWidth(1);

  Double_t* bins = binsArray.GetArray();

  hist = (TH1F*) hist->Rebin(nBins-1,"",bins);
 
  for(Int_t bin=0; bin <= hist->GetNbinsX(); bin++){

    Double_t binWidthRebin = hist->GetBinWidth(bin);
    Double_t scaleF = binWidthNoRebin / binWidthRebin;

    Double_t binCont  = hist->GetBinContent(bin);
    Double_t binErr   = hist->GetBinError(bin);
    
    binCont *= scaleF;
    binErr  *= scaleF;

    hist->SetBinContent(bin,binCont);
    hist->SetBinError(bin,binErr);
  }



  TH1F* temp = new TH1F(*hist);

  if(type == kFlagPt)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,temp,0,0);
  if(type == kFlagZ)        fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,temp,0);
  if(type == kFlagXi)       fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(jetPtSlice,0,0,temp);
  if(type == kFlagSinglePt) fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->ReplaceCorrHistos(0,temp,0,0);


  delete temp;
}
//__________________________________________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteJetSpecResponse(TString strInfile, TString strdir, TString strlist/*, TString strOutfile*/)
{ 
 
 if(fDebug>0) Printf("%s:%d -- read jet spectrum response matrix from file %s ",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){ 
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  
  if(list == 0)return; // catch strlist.Lenght() == 0;
   
  THnSparse* hn6ResponseJetPt  = (THnSparse*) list->FindObject("fhnCorrelation");
 
  Int_t axis6RecJetPt = 0;
  Int_t axis6GenJetPt = 3;

  hn6ResponseJetPt->GetAxis(axis6RecJetPt)->SetTitle("rec jet p_{T} (GeV/c)"); 
  hn6ResponseJetPt->GetAxis(axis6GenJetPt)->SetTitle("gen jet p_{T} (GeV/c)"); 
  
  Int_t nBinsRecPt    = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetNbins(); 
  Double_t loLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinLowEdge(1);
  Double_t upLimRecPt = hn6ResponseJetPt->GetAxis(axis6RecJetPt)->GetBinUpEdge(nBinsRecPt);

  Int_t nBinsGenPt    = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetNbins(); 
  Double_t loLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinLowEdge(1);
  Double_t upLimGenPt = hn6ResponseJetPt->GetAxis(axis6GenJetPt)->GetBinUpEdge(nBinsGenPt);

  Int_t nBinsTrackPt = 200;
  Int_t loLimTrackPt = 0;
  Int_t upLimTrackPt = 200;
  

  Int_t    nBinsResponse[4]  = {nBinsRecPt,nBinsTrackPt,nBinsGenPt,nBinsTrackPt};
  Double_t binMinResponse[4] = { static_cast<Double_t>(loLimRecPt), static_cast<Double_t>(loLimTrackPt), static_cast<Double_t>(loLimGenPt), static_cast<Double_t>(loLimTrackPt)};
  Double_t binMaxResponse[4] = { static_cast<Double_t>(upLimRecPt), static_cast<Double_t>(upLimTrackPt), static_cast<Double_t>(upLimGenPt), static_cast<Double_t>(upLimTrackPt)};

  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)"};
  
  THnSparseD* hn4ResponseTrackPtJetPt =  new THnSparseD("hn4ResponseTrackPtJetPt","",4,nBinsResponse,binMinResponse,binMaxResponse);
  
  for(Int_t i=0; i<4; i++){
    hn4ResponseTrackPtJetPt->GetAxis(i)->SetTitle(labelsResponseSinglePt[i]); 
  }
  

  // fill 
  

}

//_____________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadSingleTrackEfficiency(TString strfile, TString strdir, TString strlist, TString strname)
{

  ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagEfficiency);

}

//_____________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadSingleTrackResponse(TString strfile, TString strdir, TString strlist, TString strname)
{

  ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagResponse);

}

//_____________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadSingleTrackSecCorr(TString strfile, TString strdir, TString strlist, TString strname)
{

  ReadSingleTrackCorrection(strfile,strdir,strlist,strname,kFlagSecondaries);

}

//______________________________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadSingleTrackCorrection(TString strfile, TString strdir, TString strlist, TString strname, const Int_t type)
{
  // read single track correction (pt) from file
  // type: efficiency / response / secondaries correction 

  if(!((type == kFlagEfficiency)  || (type == kFlagResponse) || (type == kFlagSecondaries))){
    Printf("%s:%d -- no such correction ",(char*)__FILE__,__LINE__);
    return;
  }

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0 && type==kFlagEfficiency)  Printf("%s:%d -- read single track corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
  if(fDebug>0 && type==kFlagResponse)    Printf("%s:%d -- read single track response from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
  if(fDebug>0 && type==kFlagSecondaries) Printf("%s:%d -- read single track secondaries corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
    
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  
  
  TH1F*      h1CorrHist = 0; // common TObject pointer not possible, need SetDirectory() later
  THnSparse* hnCorrHist = 0;
  
  if(type == kFlagEfficiency || type == kFlagSecondaries){

    if(list) h1CorrHist = (TH1F*) list->FindObject(strname); 
    else     h1CorrHist = (TH1F*) gDirectory->Get(strname);

    if(!h1CorrHist){
      Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
      return;
    }

  }
  else if(type == kFlagResponse){
    
    if(list) hnCorrHist = (THnSparse*) list->FindObject(strname);
    else     hnCorrHist = (THnSparse*) gDirectory->Get(strname);

    if(!hnCorrHist){
      Printf("%s:%d -- error retrieving histo %s", (char*)__FILE__,__LINE__,strname.Data());
      return;
    }
 
  }

  if(h1CorrHist) h1CorrHist->SetDirectory(0);
  //if(hnCorrHist) hnCorrHist->SetDirectory(0);
  
  f.Close();  

  if(type == kFlagEfficiency)       fh1EffSinglePt      = h1CorrHist; 
  else if(type == kFlagResponse)    fhnResponseSinglePt = hnCorrHist;
  else if(type == kFlagSecondaries) fh1SecCorrSinglePt  = h1CorrHist; 
 
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strInfile, TString strID)
{ 
  // read track pt spec from task ouput - standard dir/list 
     
  TString strdir  = "PWGJE_FragmentationFunction_" + strID;
  TString strlist = "fracfunc_" + strID;
    
  ReadRawPtSpec(strInfile,strdir,strlist);
}

//_______________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawPtSpec(TString strfile, TString strdir, TString strlist)
{
  // get raw pt spectra from file 

  // book histos
  fNCorrectionLevelsSinglePt = 0; 
  fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
  AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum

  // get raw pt spec from input file, normalize

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read raw spectra from file %s ",(char*)__FILE__,__LINE__,strfile.Data());

  gDirectory->cd(strdir);

  TList* list = 0;
  
  if(!(list = (TList*) gDirectory->Get(strlist))){ 
    Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
    return;
  }

  TString hnameTrackPt("fh1TrackQAPtRecCuts");
  TString hnameEvtSel("fh1EvtSelection");

  TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
  TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);  

  if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
  if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
   
  Float_t nEvents = fh1EvtSel->GetBinContent(fh1EvtSel->FindBin(0));
  
  fh1TrackPt->SetDirectory(0);
  
  f.Close();  
  
  // rebin + normalize
  if(fNHistoBinsSinglePt) fh1TrackPt = (TH1F*) fh1TrackPt->Rebin(fNHistoBinsSinglePt,"",fHistoBinsSinglePt->GetArray());

  NormalizeTH1(fh1TrackPt,nEvents);

  // raw FF = corr level 0
  fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
}


//_______________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadRawPtSpecQATask(TString strfile, TString strdir, TString strlist)
{
  // get raw pt spectra from file 
  // for output from Martas QA task


  // book histos
  fNCorrectionLevelsSinglePt = 0; 
  fCorrSinglePt = new AliFragFuncCorrHistos*[fgMaxNCorrectionLevels];
  AddCorrectionLevelSinglePt(); // first 'correction' level = raw spectrum

  // get raw pt spec from input file, normalize

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read raw pt spec from QA task output file %s ",(char*)__FILE__,__LINE__,strfile.Data());

  gDirectory->cd(strdir);

  TList* list = 0;
  
  if(!(list = (TList*) gDirectory->Get(strlist))){ 
    Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
    return;
  }

  TString hnameTrackPt("fPtSel");
  TString hnameEvtSel("fNEventAll");

  TH1F* fh1TrackPt = (TH1F*) list->FindObject(hnameTrackPt);
  TH1F* fh1EvtSel  = (TH1F*) list->FindObject(hnameEvtSel);  

  if(!fh1TrackPt){ Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameTrackPt.Data()); return; }
  if(!fh1EvtSel) { Printf("%s:%d -- histo %s not found",(char*)__FILE__,__LINE__,hnameEvtSel.Data()); return; }
 

  // evts after physics selection
  Float_t nEvents = fh1EvtSel->GetEntries();
 
  fh1TrackPt->SetDirectory(0);
  
  f.Close();  

  
  NormalizeTH1(fh1TrackPt,nEvents);

  // raw FF = corr level 0
  fCorrSinglePt[0]->AddCorrHistos(0,fh1TrackPt);
}

// ________________________________________________________
void AliFragmentationFunctionCorrections::EffCorrSinglePt()
{
  // apply efficiency correction to inclusive track pt spec

  AddCorrectionLevelSinglePt("eff");

  TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
 
  if(histPt->GetNbinsX() != fh1EffSinglePt->GetNbinsX()) Printf("%s:%d: inconsistency pt spec and eff corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);

  TString histNamePt = histPt->GetName();
  TH1F* hTrackPtEffCorr = (TH1F*) histPt->Clone(histNamePt);
  hTrackPtEffCorr->Divide(histPt,fh1EffSinglePt,1,1,"");

  fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtEffCorr);
}

//___________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::UnfoldSinglePt(const Int_t nIter, const Bool_t useCorrelatedErrors)
{
  // unfolde inclusive dN/dpt spectra
  
  AddCorrectionLevelSinglePt("unfold");
   
  TH1F* hist = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0); // level -2: before unfolding, level -1: unfolded
  THnSparse* hnResponse = fhnResponseSinglePt;
      
  TString histNameTHn = hist->GetName();
  if(histNameTHn.Contains("TH1")) histNameTHn.ReplaceAll("TH1","THn");
  if(histNameTHn.Contains("fPt")) histNameTHn.ReplaceAll("fPt","fhnPt");
  

  TString histNameBackFolded = hist->GetName();
  histNameBackFolded.Append("_backfold");
  
  TString histNameRatioFolded = hist->GetName();
  if(histNameRatioFolded.Contains("fh1")) histNameRatioFolded.ReplaceAll("fh1","hRatio");
  if(histNameRatioFolded.Contains("fPt")) histNameRatioFolded.ReplaceAll("fPt","hRatioPt");
  histNameRatioFolded.Append("_unfold");
  
  TString histNameRatioBackFolded = hist->GetName();
  if(histNameRatioBackFolded.Contains("fh1")) histNameRatioBackFolded.ReplaceAll("fh1","hRatio");
  if(histNameRatioBackFolded.Contains("fPt")) histNameRatioBackFolded.ReplaceAll("fPt","hRatioPt");
  histNameRatioBackFolded.Append("_backfold");

  THnSparse* hnHist           = TH1toSparse(hist,histNameTHn,hist->GetTitle());
  THnSparse* hnFlatEfficiency = TH1toSparse(hist,"fhnEfficiency","eff",kTRUE); // could optionally also use real eff 
  
  TH1F* hUnfolded 
    = Unfold(hnHist,hnResponse,hnFlatEfficiency,nIter,useCorrelatedErrors);  
  
  //TH1F* hUnfolded = (TH1F*) hnUnfolded->Projection(0); 
  //hUnfolded->SetNameTitle(hist->GetName(),hist->GetTitle());
  
  fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hUnfolded);
  
  // backfolding: apply response matrix to unfolded spectrum
  TH1F* hBackFolded = ApplyResponse(hUnfolded,hnResponse); 
  hBackFolded->SetNameTitle(histNameBackFolded,hist->GetTitle());
  
  fh1SingleTrackPtBackFolded = hBackFolded;
  
  
  // ratio unfolded to original histo 
  TH1F* hRatioUnfolded = (TH1F*) hUnfolded->Clone(histNameRatioFolded);
  hRatioUnfolded->Reset();
  hRatioUnfolded->Divide(hUnfolded,hist,1,1,"B");
  
  fh1RatioSingleTrackPtFolded = hRatioUnfolded;
    
  
  // ratio backfolded to original histo
  TH1F* hRatioBackFolded = (TH1F*) hBackFolded->Clone(histNameRatioBackFolded);
  hRatioBackFolded->Reset();
  hRatioBackFolded->Divide(hBackFolded,hist,1,1,"B");
  
  fh1RatioSingleTrackPtBackFolded = hRatioBackFolded;
    
  delete hnHist;
  delete hnFlatEfficiency;
  
}

// ________________________________________________________
void AliFragmentationFunctionCorrections::SecCorrSinglePt()
{
  // apply efficiency correction to inclusive track pt spec

  AddCorrectionLevelSinglePt("secCorr");

  TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
  
  if(histPt->GetNbinsX() != fh1SecCorrSinglePt->GetNbinsX())
    Printf("%s:%d: inconsistency pt spec and secondaries corr bins - rebin effCorr ...", (char*)__FILE__,__LINE__);

  TString histNamePt = histPt->GetName();
  TH1F* hTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
  
  hTrackPtSecCorr->Multiply(histPt,fh1SecCorrSinglePt,1,1,"");
  
  fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtSecCorr);
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::dNdz2dNdxi()
{
  // transform dN/dz distribution into dN/dxi
  // for current corr level, all jet pt slices
  

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histZ = fCorrFF[fNCorrectionLevels-1]->GetZ(i);
    Int_t nBins = histZ->GetNbinsX();

    Double_t* binLims = new Double_t[nBins+1];
    
    for(Int_t bin = 0; bin<nBins; bin++){
      
      Int_t binZ = nBins-bin;

      Double_t zLo = histZ->GetXaxis()->GetBinLowEdge(binZ);
      Double_t zUp = histZ->GetXaxis()->GetBinUpEdge(binZ);
      
      Double_t xiLo = TMath::Log(1/zUp);    
      Double_t xiUp = TMath::Log(1/zLo);

      if(bin == 0) binLims[0] = xiLo;
      binLims[bin+1] = xiUp;   
    }

    // for(Int_t bin = 0; bin<=nBins; bin++) std::cout<<" bin "<<bin<<" binLims "<<binLims[bin]<<std::endl;

    TString strTitle = histZ->GetTitle();
    TString strName  = histZ->GetName();

    strName.ReplaceAll("Z","XiNew");
    strTitle.ReplaceAll("Z","XiNew");
    

    TH1F* histXiNew = new TH1F("histXiNew","",nBins,binLims);
    histXiNew->SetNameTitle(strName,strTitle);
    

    for(Int_t binZ = 1; binZ<=nBins; binZ++){
      
      Double_t meanZ  = histZ->GetBinCenter(binZ); 
      Double_t cont   = histZ->GetBinContent(binZ); 
      Double_t err    = histZ->GetBinError(binZ); 

      Double_t meanXi = TMath::Log(1/meanZ);   
      Int_t binXi = histXiNew->FindBin(meanXi);
      
      histXiNew->SetBinContent(binXi,cont);
      histXiNew->SetBinError(binXi,err);
      
      //std::cout<<" binZ "<<binZ<<" meanZ "<<meanZ<<" binXi "<<binXi<<" meanXi "<<meanXi<<" cont "<<cont<<" err "<<err<<std::endl;
    }
    
    fCorrFF[fNCorrectionLevels-1]->ReplaceCorrHistos(i,0,0,histXiNew);

    delete histXiNew;
    delete[] binLims;
  }
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strIDGen,  TString strIDRec,  
							    TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
							    TString strOutDir,Double_t scaleFacBgrRec)
{ 
  TString strdirGen  = "PWGJE_FragmentationFunction_" + strIDGen;
  TString strlistGen = "fracfunc_" + strIDGen;
  
  TString strdirRec  = "PWGJE_FragmentationFunction_" + strIDRec;
  TString strlistRec = "fracfunc_" + strIDRec;
  
  WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kFALSE, "",scaleFacBgrRec);
}

//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteBgrBinShiftCorr(TString strInfile, TString strBgrID, TString strIDGen,  TString strIDRec,  
							       TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
							       TString strOutDir,Double_t scaleFacBgrRec)
{ 
  TString strdirGen  = "PWGJE_FragmentationFunction_" + strIDGen;
  TString strlistGen = "fracfunc_" + strIDGen;
  
  TString strdirRec  = "PWGJE_FragmentationFunction_" + strIDRec;
  TString strlistRec = "fracfunc_" + strIDRec;
  
  WriteBinShiftCorr(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim, strOutDir, kTRUE, strBgrID,scaleFacBgrRec);
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteBinShiftCorr(TString strInfile, TString strdirGen, TString strlistGen, 
							    TString strdirRec, TString strlistRec, 
							    TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
							    TString strOutDir,  Bool_t writeBgr, TString strBgrID, Double_t scaleFacBgrRec)
{
 
  if((writeBgr && strBgrID.Length() == 0) || (!writeBgr && strBgrID.Length()>0) ){
    Printf("%s:%d -- inconsistent arguments to WriteBinShiftCorr FF/UE", (char*)__FILE__,__LINE__);
    return; 
  }

  TH1F* hCorrPt[fNJetPtSlices];
  TH1F* hCorrXi[fNJetPtSlices];
  TH1F* hCorrZ[fNJetPtSlices];

  TH1F* hdNdptTracksMCGen[fNJetPtSlices];
  TH1F* hdNdxiMCGen[fNJetPtSlices];
  TH1F* hdNdzMCGen[fNJetPtSlices];
    
  TH1F* hdNdptTracksMCRec[fNJetPtSlices];
  TH1F* hdNdxiMCRec[fNJetPtSlices];
  TH1F* hdNdzMCRec[fNJetPtSlices];

  TH1F* fh1FFJetPtMCGen = 0;
  TH1F* fh1FFJetPtMCRec = 0;

  TH2F* fh2FFTrackPtMCGen = 0;
  TH2F* fh2FFZMCGen       = 0;
  TH2F* fh2FFXiMCGen      = 0;
  
  TH2F* fh2FFTrackPtMCRec = 0;
  TH2F* fh2FFZMCRec       = 0;
  TH2F* fh2FFXiMCRec      = 0;
 
  // gen level FF

  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeBinShiftCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdirGen && strdirGen.Length()) gDirectory->cd(strdirGen);

  TList* listGen = 0;

  if(strlistGen && strlistGen.Length()){

    if(!(listGen = (TList*) gDirectory->Get(strlistGen))){ 
      Printf("%s:%d -- error retrieving listGen %s from directory %s", (char*)__FILE__,__LINE__,strlistGen.Data(),strdirGen.Data());
      return;
    }
  }
 
  if(listGen){ 

    fh1FFJetPtMCGen    = (TH1F*) listGen->FindObject("fh1FFJetPtGen");
    
    if(writeBgr){
      fh2FFTrackPtMCGen  = (TH2F*) listGen->FindObject(Form("fh2FFTrackPt%sGen",strBgrID.Data()));
      fh2FFZMCGen        = (TH2F*) listGen->FindObject(Form("fh2FFZ%sGen",strBgrID.Data()));
      fh2FFXiMCGen       = (TH2F*) listGen->FindObject(Form("fh2FFXi%sGen",strBgrID.Data()));  
    }
    else{
      fh2FFTrackPtMCGen  = (TH2F*) listGen->FindObject("fh2FFTrackPtGen");
      fh2FFZMCGen        = (TH2F*) listGen->FindObject("fh2FFZGen");
      fh2FFXiMCGen       = (TH2F*) listGen->FindObject("fh2FFXiGen"); 
    }
  }
  else{
    fh1FFJetPtMCGen    = (TH1F*) gDirectory->Get("fh1FFJetPtGen");
    if(writeBgr){
      fh2FFTrackPtMCGen  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sGen",strBgrID.Data()));
      fh2FFZMCGen        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sGen",strBgrID.Data()));
      fh2FFXiMCGen       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sGen",strBgrID.Data())); 
    }
    else{
      fh2FFTrackPtMCGen  = (TH2F*) gDirectory->Get("fh2FFTrackPtGen");
      fh2FFZMCGen        = (TH2F*) gDirectory->Get("fh2FFZGen");
      fh2FFXiMCGen       = (TH2F*) gDirectory->Get("fh2FFXiGen"); 
    }
  }
  
  fh1FFJetPtMCGen->SetDirectory(0); 
  fh2FFTrackPtMCGen->SetDirectory(0);
  fh2FFZMCGen->SetDirectory(0);
  fh2FFXiMCGen->SetDirectory(0);

  f.Close();

  // rec level FF

  TFile g(strInfile,"READ");

  if(!g.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeBinShiftCorr: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdirRec && strdirRec.Length()) gDirectory->cd(strdirRec);

  TList* listRec = 0;

  if(strlistRec && strlistRec.Length()){

    if(!(listRec = (TList*) gDirectory->Get(strlistRec))){ 
      Printf("%s:%d -- error retrieving listRec %s from directory %s", (char*)__FILE__,__LINE__,strlistRec.Data(),strdirRec.Data());
      return;
    }
  }
 
  if(useRecPrim){
    if(listRec){ 
      fh1FFJetPtMCRec    = (TH1F*) listRec->FindObject("fh1FFJetPtRecEffRec");

      if(writeBgr){
	fh2FFTrackPtMCRec  = (TH2F*) listRec->FindObject(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
	fh2FFZMCRec        = (TH2F*) listRec->FindObject(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
	fh2FFXiMCRec       = (TH2F*) listRec->FindObject(Form("fh2FFXi%sRecEffRec",strBgrID.Data())); 
      }
      else{
	fh2FFTrackPtMCRec  = (TH2F*) listRec->FindObject("fh2FFTrackPtRecEffRec");
	fh2FFZMCRec        = (TH2F*) listRec->FindObject("fh2FFZRecEffRec");
	fh2FFXiMCRec       = (TH2F*) listRec->FindObject("fh2FFXiRecEffRec"); 
      }
    }
    else{
      fh1FFJetPtMCRec    = (TH1F*) gDirectory->Get("fh1FFJetPtRecEffRec");
      if(writeBgr){
	fh2FFTrackPtMCRec  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecEffRec",strBgrID.Data()));
	fh2FFZMCRec        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecEffRec",strBgrID.Data()));
	fh2FFXiMCRec       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecEffRec",strBgrID.Data())); 
      }
      else{
	fh2FFTrackPtMCRec  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecEffRec");
	fh2FFZMCRec        = (TH2F*) gDirectory->Get("fh2FFZRecEffRec");
	fh2FFXiMCRec       = (TH2F*) gDirectory->Get("fh2FFXiRecEffRec"); 
      }
    }
  }
  else{
    if(listRec){ 
      fh1FFJetPtMCRec    = (TH1F*) listRec->FindObject("fh1FFJetPtRecCuts");
      if(writeBgr){
	fh2FFTrackPtMCRec  = (TH2F*) listRec->FindObject(Form("fh2FFTrackPt%sRecCuts",strBgrID.Data()));
	fh2FFZMCRec        = (TH2F*) listRec->FindObject(Form("fh2FFZ%sRecCuts",strBgrID.Data()));
	fh2FFXiMCRec       = (TH2F*) listRec->FindObject(Form("fh2FFXi%sRecCuts",strBgrID.Data())); 
      }
      else{
	fh2FFTrackPtMCRec  = (TH2F*) listRec->FindObject("fh2FFTrackPtRecCuts");
	fh2FFZMCRec        = (TH2F*) listRec->FindObject("fh2FFZRecCuts");
	fh2FFXiMCRec       = (TH2F*) listRec->FindObject("fh2FFXiRecCuts");
      }
    }
    else{
      fh1FFJetPtMCRec    = (TH1F*) gDirectory->Get("fh1FFJetPtRecCuts");
      if(writeBgr){
	fh2FFTrackPtMCRec  = (TH2F*) gDirectory->Get(Form("fh2FFTrackPt%sRecCuts",strBgrID.Data()));
	fh2FFZMCRec        = (TH2F*) gDirectory->Get(Form("fh2FFZ%sRecCuts",strBgrID.Data()));
	fh2FFXiMCRec       = (TH2F*) gDirectory->Get(Form("fh2FFXi%sRecCuts",strBgrID.Data()));  
      }
      else{
	fh2FFTrackPtMCRec  = (TH2F*) gDirectory->Get("fh2FFTrackPtRecCuts");
	fh2FFZMCRec        = (TH2F*) gDirectory->Get("fh2FFZRecCuts");
	fh2FFXiMCRec       = (TH2F*) gDirectory->Get("fh2FFXiRecCuts"); 
      }
    }
  }
  

  fh1FFJetPtMCRec->SetDirectory(0); 
  fh2FFTrackPtMCRec->SetDirectory(0);
  fh2FFZMCRec->SetDirectory(0);
  fh2FFXiMCRec->SetDirectory(0);  

  g.Close();

  // projections: FF for generated and reconstructed  
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Float_t jetPtLoLim = fJetPtSlices->At(i);
    Float_t jetPtUpLim = fJetPtSlices->At(i+1);
    
    Int_t binLo = static_cast<Int_t>(fh2FFTrackPtMCGen->GetXaxis()->FindBin(jetPtLoLim));
    Int_t binUp = static_cast<Int_t>(fh2FFTrackPtMCGen->GetXaxis()->FindBin(jetPtUpLim))-1;

    if(binUp > fh2FFTrackPtMCGen->GetNbinsX()){
      Printf("%s:%d -- jet pt range %0.3f exceeds histo limits",(char*)__FILE__,__LINE__,jetPtUpLim); 
      return; 
    }
    
    TString strNameFFPtGen(Form("fh1FFTrackPtGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZGen(Form("fh1FFZGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiGen(Form("fh1FFXiGenBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    
    TString strNameFFPtRec(Form("fh1FFTrackPtRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFZRec(Form("fh1FFZRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));
    TString strNameFFXiRec(Form("fh1FFXiRecBinShift_%02d_%02d",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim)));

    if(writeBgr){
      strNameFFPtGen.Form("fh1TrackPtGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
      strNameFFZGen.Form("fh1ZGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
      strNameFFXiGen.Form("fh1XiGenBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
      
      strNameFFPtRec.Form("fh1TrackPtRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
      strNameFFZRec.Form("fh1ZRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
      strNameFFXiRec.Form("fh1XiRecBinShiftBgr_%02d_%02d_%s",static_cast<Int_t> (jetPtLoLim),static_cast<Int_t> (jetPtUpLim),strBgrID.Data());
    }
    
    // project 
    // appendix 'unbinned' to avoid histos with same name after rebinning

    hdNdptTracksMCGen[i] = (TH1F*) fh2FFTrackPtMCGen->ProjectionY(strNameFFPtGen+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCGen[i]        = (TH1F*) fh2FFZMCGen->ProjectionY(strNameFFZGen+"_unBinned",binLo,binUp,"o");
    hdNdxiMCGen[i]       = (TH1F*) fh2FFXiMCGen->ProjectionY(strNameFFXiGen+"_unBinned",binLo,binUp,"o");
    
    hdNdptTracksMCRec[i] = (TH1F*) fh2FFTrackPtMCRec->ProjectionY(strNameFFPtRec+"_unBinned",binLo,binUp,"o"); // option "o": original axis range 
    hdNdzMCRec[i]        = (TH1F*) fh2FFZMCRec->ProjectionY(strNameFFZRec+"_unBinned",binLo,binUp,"o");
    hdNdxiMCRec[i]       = (TH1F*) fh2FFXiMCRec->ProjectionY(strNameFFXiRec+"_unBinned",binLo,binUp,"o");
    

    // rebin

    if(fNHistoBinsPt[i]) hdNdptTracksMCGen[i] = (TH1F*) hdNdptTracksMCGen[i]->Rebin(fNHistoBinsPt[i],strNameFFPtGen,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCGen[i]  = (TH1F*) hdNdzMCGen[i]->Rebin(fNHistoBinsZ[i],strNameFFZGen,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCGen[i] = (TH1F*) hdNdxiMCGen[i]->Rebin(fNHistoBinsXi[i],strNameFFXiGen,fHistoBinsXi[i]->GetArray());

    if(fNHistoBinsPt[i]) hdNdptTracksMCRec[i] = (TH1F*) hdNdptTracksMCRec[i]->Rebin(fNHistoBinsPt[i],strNameFFPtRec,fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hdNdzMCRec[i]  = (TH1F*) hdNdzMCRec[i]->Rebin(fNHistoBinsZ[i],strNameFFZRec,fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hdNdxiMCRec[i] = (TH1F*) hdNdxiMCRec[i]->Rebin(fNHistoBinsXi[i],strNameFFXiRec,fHistoBinsXi[i]->GetArray());


    hdNdptTracksMCGen[i]->SetNameTitle(strNameFFPtGen,"");
    hdNdzMCGen[i]->SetNameTitle(strNameFFZGen,"");
    hdNdxiMCGen[i]->SetNameTitle(strNameFFXiGen,"");
    
    hdNdptTracksMCRec[i]->SetNameTitle(strNameFFPtRec,"");
    hdNdzMCRec[i]->SetNameTitle(strNameFFZRec,"");
    hdNdxiMCRec[i]->SetNameTitle(strNameFFXiRec,"");
    
     // normalize
    
    Double_t nJetsBinGen = fh1FFJetPtMCGen->Integral(binLo,binUp);
    Double_t nJetsBinRec = fh1FFJetPtMCRec->Integral(binLo,binUp);

    // scale fac for perp2 bgr
     if(useRecPrim && writeBgr && scaleFacBgrRec && (scaleFacBgrRec != 1)) nJetsBinRec /= scaleFacBgrRec;

    NormalizeTH1(hdNdptTracksMCGen[i],nJetsBinGen); 
    NormalizeTH1(hdNdzMCGen[i],nJetsBinGen); 
    NormalizeTH1(hdNdxiMCGen[i],nJetsBinGen); 

    NormalizeTH1(hdNdptTracksMCRec[i],nJetsBinRec); 
    NormalizeTH1(hdNdzMCRec[i],nJetsBinRec); 
    NormalizeTH1(hdNdxiMCRec[i],nJetsBinRec); 

    // divide gen/rec : corr factor

    TString strNameCorrPt(Form("hBbBCorrPt_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameCorrZ(Form("hBbBCorrZ_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));
    TString strNameCorrXi(Form("hBbBCorrXi_%02d_%02d",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim)));

    if(writeBgr){
      strNameCorrPt.Form("hBbBCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameCorrZ.Form("hBbBCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameCorrXi.Form("hBbBCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
    }

    hCorrPt[i] = (TH1F*) hdNdptTracksMCGen[i]->Clone(strNameCorrPt);
    hCorrPt[i]->Divide(hdNdptTracksMCGen[i],hdNdptTracksMCRec[i],1,1,"B"); // binominal errors
    
    hCorrXi[i] = (TH1F*) hdNdxiMCGen[i]->Clone(strNameCorrXi);
    hCorrXi[i]->Divide(hdNdxiMCGen[i],hdNdxiMCRec[i],1,1,"B"); // binominal errors
    
    hCorrZ[i] = (TH1F*) hdNdzMCGen[i]->Clone(strNameCorrZ);
    hCorrZ[i]->Divide(hdNdzMCGen[i],hdNdzMCRec[i],1,1,"B"); // binominal errors
  } 
  
  // write 

  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write bin shift correction to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

  if(strOutDir && strOutDir.Length()){  
    
    TDirectory* dir;
    if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }

  for(Int_t i=0; i<fNJetPtSlices; i++){

    hdNdptTracksMCGen[i]->Write();
    hdNdxiMCGen[i]->Write();
    hdNdzMCGen[i]->Write();
    
    hdNdptTracksMCRec[i]->Write();
    hdNdxiMCRec[i]->Write();
    hdNdzMCRec[i]->Write();
  
    hCorrPt[i]->Write();
    hCorrXi[i]->Write();
    hCorrZ[i]->Write();
  }

  // out.Close();
 
  delete fh1FFJetPtMCGen; 
  delete fh1FFJetPtMCRec; 

  delete fh2FFTrackPtMCGen;
  delete fh2FFZMCGen;
  delete fh2FFXiMCGen;
  
  delete fh2FFTrackPtMCRec;
  delete fh2FFZMCRec;
  delete fh2FFXiMCRec;
}

//________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadBgrBinShiftCorr(TString strfile, TString strBgrID, TString strdir, TString strlist)
{

  ReadBinShiftCorr(strfile, strdir, strlist, kTRUE, strBgrID);
}

//___________________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadBinShiftCorr(TString strfile, TString strdir, TString strlist, Bool_t readBgr, TString strBgrID)
{

  if((readBgr && strBgrID.Length() == 0) || (!readBgr && strBgrID.Length()>0) ){
    Printf("%s:%d -- inconsistent arguments to ReadBinShiftCorr FF/UE", (char*)__FILE__,__LINE__);
    return; 
  }

  // temporary histos to hold histos from file
  TH1F* hCorrPt[fNJetPtSlices]; 
  TH1F* hCorrZ[fNJetPtSlices];
  TH1F* hCorrXi[fNJetPtSlices];
  
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i]  = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read FF / UE bin shift correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
   
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
    
    TString strNameCorrPt(Form("hBbBCorrPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameCorrZ(Form("hBbBCorrZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameCorrXi(Form("hBbBCorrXi_%02d_%02d",jetPtLoLim,jetPtUpLim));

    if(readBgr){
      strNameCorrPt.Form("hBbBCorrBgrPt_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameCorrZ.Form("hBbBCorrBgrZ_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
      strNameCorrXi.Form("hBbBCorrBgrXi_%02d_%02d_%s",static_cast<Int_t>(jetPtLoLim),static_cast<Int_t>(jetPtUpLim),strBgrID.Data());
    }

   
    if(list){
      hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt); 
      hCorrZ[i]  = (TH1F*) list->FindObject(strNameCorrZ); 
      hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi); 
    }
    else{
      hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt); 
      hCorrZ[i]  = (TH1F*) gDirectory->Get(strNameCorrZ); 
      hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi); 
    }
    
    if(!hCorrPt[i]){
      Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
    }
  
    if(!hCorrZ[i]){
      Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
    }    

    if(!hCorrXi[i]){
      Printf("%s:%d -- error retrieving bin by bin correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
    }


    if(fNHistoBinsPt[i]) hCorrPt[i] = (TH1F*) hCorrPt[i]->Rebin(fNHistoBinsPt[i],strNameCorrPt+"_rebin",fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hCorrZ[i]  = (TH1F*) hCorrZ[i]->Rebin(fNHistoBinsZ[i],strNameCorrZ+"_rebin",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hCorrXi[i] = (TH1F*) hCorrXi[i]->Rebin(fNHistoBinsXi[i],strNameCorrXi+"_rebin",fHistoBinsXi[i]->GetArray());

    if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0); 
    if(hCorrZ[i])  hCorrZ[i]->SetDirectory(0); 
    if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0); 

  } // jet slices loop

  f.Close();

  if(readBgr){
    for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
      if(hCorrPt[i]) new(fh1BbBBgrPt[i]) TH1F(*hCorrPt[i]);
      if(hCorrZ[i])  new(fh1BbBBgrZ[i])  TH1F(*hCorrZ[i]);
      if(hCorrXi[i]) new(fh1BbBBgrXi[i]) TH1F(*hCorrXi[i]);
    }
  }
  else{
    for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
      if(hCorrPt[i]) new(fh1BbBPt[i]) TH1F(*hCorrPt[i]);
      if(hCorrZ[i])  new(fh1BbBZ[i])  TH1F(*hCorrZ[i]);
      if(hCorrXi[i]) new(fh1BbBXi[i]) TH1F(*hCorrXi[i]);
    }
  }
}

// ________________________________________________
void AliFragmentationFunctionCorrections::BbBCorr()
{
  // apply bin-by-bin correction

  AddCorrectionLevel("BbB");

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
    TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();
    
    TH1F* hFFTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
    hFFTrackPtBbBCorr->Multiply(histPt,fh1BbBPt[i],1,1,"");
    
    TH1F* hFFZBbBCorr = (TH1F*) histZ->Clone(histNameZ);
    hFFZBbBCorr->Multiply(histZ,fh1BbBZ[i],1,1,"");
    
    TH1F* hFFXiBbBCorr = (TH1F*) histXi->Clone(histNameXi);
    hFFXiBbBCorr->Multiply(histXi,fh1BbBXi[i],1,1,"");
    
    fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtBbBCorr,hFFZBbBCorr,hFFXiBbBCorr);
  }
}

// ___________________________________________________
void AliFragmentationFunctionCorrections::BbBCorrBgr()
{
  // apply bin-by-bin correction

  AddCorrectionLevelBgr("BbB");

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
    TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();
    
    TH1F* hBgrTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
    hBgrTrackPtBbBCorr->Multiply(histPt,fh1BbBBgrPt[i],1,1,"");
    
    TH1F* hBgrZBbBCorr = (TH1F*) histZ->Clone(histNameZ);
    hBgrZBbBCorr->Multiply(histZ,fh1BbBBgrZ[i],1,1,"");
    
    TH1F* hBgrXiBbBCorr = (TH1F*) histXi->Clone(histNameXi);
    hBgrXiBbBCorr->Multiply(histXi,fh1BbBBgrXi[i],1,1,"");
    
    fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hBgrTrackPtBbBCorr,hBgrZBbBCorr,hBgrXiBbBCorr);
  }
}

//_______________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadFoldingCorr(TString strfile, TString strdir, TString strlist)
{

  // temporary histos to hold histos from file
  TH1F* hCorrPt[fNJetPtSlices]; 
  TH1F* hCorrZ[fNJetPtSlices];
  TH1F* hCorrXi[fNJetPtSlices];
  
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i]  = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read bin shift correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
   
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));
    
    TString strNameCorrPt(Form("hCorrFoldingPt_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameCorrZ(Form("hCorrFoldingZ_%02d_%02d",jetPtLoLim,jetPtUpLim));
    TString strNameCorrXi(Form("hCorrFoldingXi_%02d_%02d",jetPtLoLim,jetPtUpLim));
     
   
    if(list){
      hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt); 
      hCorrZ[i]  = (TH1F*) list->FindObject(strNameCorrZ); 
      hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi); 
    }
    else{
      hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt); 
      hCorrZ[i]  = (TH1F*) gDirectory->Get(strNameCorrZ); 
      hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi); 
    }
    
    if(!hCorrPt[i]){
      //Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
    }
  
    if(!hCorrZ[i]){
      //Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
    }    

    if(!hCorrXi[i]){
      Printf("%s:%d -- error retrieving folding correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
    }

    if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0); 
    if(hCorrZ[i])  hCorrZ[i]->SetDirectory(0); 
    if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0); 

  } // jet slices loop

  f.Close();

  for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos
    if(hCorrPt[i]) new(fh1FoldingCorrPt[i]) TH1F(*hCorrPt[i]);
    if(hCorrZ[i])  new(fh1FoldingCorrZ[i])  TH1F(*hCorrZ[i]);
    if(hCorrXi[i]) new(fh1FoldingCorrXi[i]) TH1F(*hCorrXi[i]);
  }
}

// ___________________________________________________
void AliFragmentationFunctionCorrections::FoldingCorr()
{
  // apply bin-by-bin correction

  AddCorrectionLevel("FoldingCorr");

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
    TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();
        
    std::cout<<" foldingCorr: i "<<i<<" corr pt "<<fh1FoldingCorrPt[i]<<" z "<<fh1FoldingCorrZ[i]<<" xi "
	     <<fh1FoldingCorrXi[i]<<std::endl;

    std::cout<<" foldingCorr: i "<<i<<" mean corr pt "<<fh1FoldingCorrPt[i]->GetMean()<<" z "<<fh1FoldingCorrZ[i]->GetMean()<<" xi "
	     <<fh1FoldingCorrXi[i]->GetMean()<<std::endl;

    

    TH1F* hFFTrackPtFoldingCorr = (TH1F*) histPt->Clone(histNamePt);
    if(fh1FoldingCorrPt[i] && fh1FoldingCorrPt[i]->GetMean()>0) hFFTrackPtFoldingCorr->Multiply(histPt,fh1FoldingCorrPt[i],1,1,"");
    else hFFTrackPtFoldingCorr->Reset();

    TH1F* hFFZFoldingCorr = (TH1F*) histZ->Clone(histNameZ);
    if(fh1FoldingCorrZ[i] && fh1FoldingCorrZ[i]->GetMean()>0) hFFZFoldingCorr->Multiply(histZ,fh1FoldingCorrZ[i],1,1,"");
    else hFFZFoldingCorr->Reset();
    
    TH1F* hFFXiFoldingCorr = (TH1F*) histXi->Clone(histNameXi);
    if(fh1FoldingCorrXi[i]&& fh1FoldingCorrXi[i]->GetMean()>0) hFFXiFoldingCorr->Multiply(histXi,fh1FoldingCorrXi[i],1,1,"");
    else hFFXiFoldingCorr->Reset();
    
    fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtFoldingCorr,hFFZFoldingCorr,hFFXiFoldingCorr);
  }
}

//________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadBgrJetSecCorr(TString strfile, TString strBgrID, TString strdir, TString strlist, 
							    Bool_t useScaledStrangeness)
{

  ReadJetSecCorr(strfile, strdir, strlist, useScaledStrangeness, kTRUE, strBgrID);
}

//_______________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadJetSecCorr(TString strfile, TString strdir, TString strlist, Bool_t useScaledStrangeness, 
							 Bool_t readBgr, TString strBgrID){

  // read reconstruction efficiency from file
  // argument strlist optional - read from directory strdir if not specified

  // temporary histos to hold histos from file
  TH1F* hCorrPt[fNJetPtSlices]; 
  TH1F* hCorrZ[fNJetPtSlices];
  TH1F* hCorrXi[fNJetPtSlices];
  
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrPt[i] = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrZ[i]  = 0;
  for(Int_t i=0; i<fNJetPtSlices; i++) hCorrXi[i] = 0;

  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read secondary correction from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
   
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    Int_t jetPtLoLim = static_cast<Int_t> (fJetPtSlices->At(i));
    Int_t jetPtUpLim = static_cast<Int_t> (fJetPtSlices->At(i+1));

    TString strNameCorrPt("");
    TString strNameCorrZ("");
    TString strNameCorrXi("");
    
    if(readBgr){
      strNameCorrPt.Form("hSecCorrBgrPt_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
      strNameCorrZ.Form("hSecCorrBgrZ_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
      strNameCorrXi.Form("hSecCorrBgrXi_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
      
      if(!useScaledStrangeness){
	Printf("%s:%d -- readJetSecCorr bgr: using naive Pythia strangenss ", (char*)__FILE__,__LINE__);
	strNameCorrPt.Form("hSecCorrBgrPt_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
	strNameCorrZ.Form("hSecCorrBgrZ_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
	strNameCorrXi.Form("hSecCorrBgrXi_nonSc_%02d_%02d_%s",jetPtLoLim,jetPtUpLim,strBgrID.Data());
      }
    }
    else{
      strNameCorrPt.Form("hSecCorrPt_%02d_%02d",jetPtLoLim,jetPtUpLim);
      strNameCorrZ.Form("hSecCorrZ_%02d_%02d",jetPtLoLim,jetPtUpLim);
      strNameCorrXi.Form("hSecCorrXi_%02d_%02d",jetPtLoLim,jetPtUpLim);
      
      if(!useScaledStrangeness){
	Printf("%s:%d -- readJetSecCorr: using naive Pythia strangenss ", (char*)__FILE__,__LINE__);
	strNameCorrPt.Form("hSecCorrPt_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
	strNameCorrZ.Form("hSecCorrZ_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
	strNameCorrXi.Form("hSecCorrXi_nonSc_%02d_%02d",jetPtLoLim,jetPtUpLim);
      }
    }
   
    if(list){
      hCorrPt[i] = (TH1F*) list->FindObject(strNameCorrPt); 
      hCorrZ[i]  = (TH1F*) list->FindObject(strNameCorrZ); 
      hCorrXi[i] = (TH1F*) list->FindObject(strNameCorrXi); 
    }
    else{
      hCorrPt[i] = (TH1F*) gDirectory->Get(strNameCorrPt); 
      hCorrZ[i]  = (TH1F*) gDirectory->Get(strNameCorrZ); 
      hCorrXi[i] = (TH1F*) gDirectory->Get(strNameCorrXi); 
    }
    
    if(!hCorrPt[i]){
      Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrPt.Data());
    }
  
    if(!hCorrZ[i]){
      Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrZ.Data());
    }    

    if(!hCorrXi[i]){
      Printf("%s:%d -- error retrieving secondaries correction %s", (char*)__FILE__,__LINE__,strNameCorrXi.Data());
    }

    if(fNHistoBinsPt[i]) hCorrPt[i] = (TH1F*) hCorrPt[i]->Rebin(fNHistoBinsPt[i],strNameCorrPt+"_rebin",fHistoBinsPt[i]->GetArray());
    if(fNHistoBinsZ[i])  hCorrZ[i]  = (TH1F*) hCorrZ[i]->Rebin(fNHistoBinsZ[i],strNameCorrZ+"_rebin",fHistoBinsZ[i]->GetArray());
    if(fNHistoBinsXi[i]) hCorrXi[i] = (TH1F*) hCorrXi[i]->Rebin(fNHistoBinsXi[i],strNameCorrXi+"_rebin",fHistoBinsXi[i]->GetArray());

    if(hCorrPt[i]) hCorrPt[i]->SetDirectory(0); 
    if(hCorrZ[i])  hCorrZ[i]->SetDirectory(0); 
    if(hCorrXi[i]) hCorrXi[i]->SetDirectory(0); 

  } // jet slices loop

  f.Close();

  for(Int_t i=0; i<fNJetPtSlices; i++){ // 2nd loop: need to close input file before placing histos

    if(readBgr){
      if(hCorrPt[i]) new(fh1SecCorrBgrPt[i]) TH1F(*hCorrPt[i]);
      if(hCorrZ[i])  new(fh1SecCorrBgrZ[i])  TH1F(*hCorrZ[i]);
      if(hCorrXi[i]) new(fh1SecCorrBgrXi[i]) TH1F(*hCorrXi[i]);
    }
    else{
      if(hCorrPt[i]) new(fh1SecCorrPt[i]) TH1F(*hCorrPt[i]);
      if(hCorrZ[i])  new(fh1SecCorrZ[i])  TH1F(*hCorrZ[i]);
      if(hCorrXi[i]) new(fh1SecCorrXi[i]) TH1F(*hCorrXi[i]);
    }
  }
}

// ___________________________________________________
void AliFragmentationFunctionCorrections::JetSecCorr()
{
  // apply secondaries correction

  AddCorrectionLevel("SecCorr");

  Printf("%s:%d -- apply jet secondaries correction", (char*)__FILE__,__LINE__);

  for(Int_t i=0; i<fNJetPtSlices; i++){

    TH1F* histPt = fCorrFF[fNCorrectionLevels-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrFF[fNCorrectionLevels-2]->GetZ(i);
    TH1F* histXi = fCorrFF[fNCorrectionLevels-2]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();
    
    TH1F* hFFTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
    hFFTrackPtSecCorr->Multiply(histPt,fh1SecCorrPt[i],1,1,"");
    
    TH1F* hFFZSecCorr = (TH1F*) histZ->Clone(histNameZ);
    hFFZSecCorr->Multiply(histZ,fh1SecCorrZ[i],1,1,"");
    
    TH1F* hFFXiSecCorr = (TH1F*) histXi->Clone(histNameXi);
    hFFXiSecCorr->Multiply(histXi,fh1SecCorrXi[i],1,1,"");
    
    fCorrFF[fNCorrectionLevels-1]->AddCorrHistos(i,hFFTrackPtSecCorr,hFFZSecCorr,hFFXiSecCorr);
  }
}



// ___________________________________________________
void AliFragmentationFunctionCorrections::JetSecCorrBgr()
{
  // apply secondaries correction to UE
  
  AddCorrectionLevelBgr("SecCorr"); 

  Printf("%s:%d -- apply jet secondaries correction", (char*)__FILE__,__LINE__);
  
  for(Int_t i=0; i<fNJetPtSlices; i++){
    
    TH1F* histPt = fCorrBgr[fNCorrectionLevelsBgr-2]->GetTrackPt(i);
    TH1F* histZ  = fCorrBgr[fNCorrectionLevelsBgr-2]->GetZ(i);
    TH1F* histXi = fCorrBgr[fNCorrectionLevelsBgr-2]->GetXi(i);

    TString histNamePt = histPt->GetName();
    TString histNameZ  = histZ->GetName();
    TString histNameXi = histXi->GetName();
    
    TH1F* hFFTrackPtSecCorr = (TH1F*) histPt->Clone(histNamePt);
    hFFTrackPtSecCorr->Multiply(histPt,fh1SecCorrBgrPt[i],1,1,"");
    
    TH1F* hFFZSecCorr = (TH1F*) histZ->Clone(histNameZ);
    hFFZSecCorr->Multiply(histZ,fh1SecCorrBgrZ[i],1,1,"");
    
    TH1F* hFFXiSecCorr = (TH1F*) histXi->Clone(histNameXi);
    hFFXiSecCorr->Multiply(histXi,fh1SecCorrBgrXi[i],1,1,"");
    
    fCorrBgr[fNCorrectionLevelsBgr-1]->AddCorrHistos(i,hFFTrackPtSecCorr,hFFZSecCorr,hFFXiSecCorr);
  }
}


//________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteBinShiftCorrSinglePt(TString strInfile, TString strIDGen,  TString strIDRec,  
								    TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
								    TString strOutDir)
{ 
  TString strdirGen  = "PWGJE_FragmentationFunction_" + strIDGen;
  TString strlistGen = "fracfunc_" + strIDGen;
  
  TString strdirRec  = "PWGJE_FragmentationFunction_" + strIDRec;
  TString strlistRec = "fracfunc_" + strIDRec;
  
  WriteBinShiftCorrSinglePt(strInfile,strdirGen,strlistGen,strdirRec,strlistRec,strOutfile,updateOutfile,useRecPrim,strOutDir);
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::WriteBinShiftCorrSinglePt(TString strInfile, TString strdirGen, TString strlistGen, 
								    TString strdirRec, TString strlistRec, 
								    TString strOutfile, Bool_t updateOutfile, Bool_t useRecPrim, 
								    TString strOutDir){
 
  
  TH1F* hCorrPt           = 0;
  TH1F* hdNdptTracksMCGen = 0;
  TH1F* hdNdptTracksMCRec = 0;
 
  // gen level FF

  TFile f(strInfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeBinShiftCorrSinglePt: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdirGen && strdirGen.Length()) gDirectory->cd(strdirGen);

  TList* listGen = 0;

  if(strlistGen && strlistGen.Length()){

    if(!(listGen = (TList*) gDirectory->Get(strlistGen))){ 
      Printf("%s:%d -- error retrieving listGen %s from directory %s", (char*)__FILE__,__LINE__,strlistGen.Data(),strdirGen.Data());
      return;
    }
  }
 
  if(listGen){
    hdNdptTracksMCGen  = (TH1F*) listGen->FindObject("fh1TrackQAPtGen");
  }
  else{
    hdNdptTracksMCGen  = (TH1F*) gDirectory->Get("fh1TrackQAPtGen");
  }

  hdNdptTracksMCGen->SetDirectory(0);
  
  f.Close();

  // rec level FF

  TFile g(strInfile,"READ");

  if(!g.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strInfile.Data());
    return;
  }
  
  if(fDebug>0) Printf("%s:%d -- writeBinShiftCorrSinglePt: open task ouput file %s",(char*)__FILE__,__LINE__,strInfile.Data());
 
  if(strdirRec && strdirRec.Length()) gDirectory->cd(strdirRec);

  TList* listRec = 0;

  if(strlistRec && strlistRec.Length()){

    if(!(listRec = (TList*) gDirectory->Get(strlistRec))){ 
      Printf("%s:%d -- error retrieving listRec %s from directory %s", (char*)__FILE__,__LINE__,strlistRec.Data(),strdirRec.Data());
      return;
    }
  }
 

  if(useRecPrim){
    if(listRec){ 
      hdNdptTracksMCRec = (TH1F*) listRec->FindObject("fh1TrackQAPtRecEffRec");
    }
    else{
      hdNdptTracksMCRec = (TH1F*) gDirectory->Get("fh1TrackQAPtRecEffRec");
    }
  }
  else{
    if(listRec){
      hdNdptTracksMCRec  = (TH1F*) listRec->FindObject("fh1TrackQAPtRecCuts");
    } 
    else{
      hdNdptTracksMCRec  = (TH1F*) gDirectory->Get("fh1TrackQAPtRecCuts");
    }
  }
  
  hdNdptTracksMCRec->SetDirectory(0);
  
  g.Close();
    
  TString strNamePtGen = "fh1SinglePtGenBbB";
  TString strNamePtRec = "fh1SinglePtRecBbB";

  // rebin
  if(fNHistoBinsSinglePt) hdNdptTracksMCGen = (TH1F*) hdNdptTracksMCGen->Rebin(fNHistoBinsSinglePt,strNamePtGen+"_rebin",fHistoBinsSinglePt->GetArray());
  if(fNHistoBinsSinglePt) hdNdptTracksMCRec = (TH1F*) hdNdptTracksMCRec->Rebin(fNHistoBinsSinglePt,strNamePtRec+"_rebin",fHistoBinsSinglePt->GetArray());
  
  hdNdptTracksMCGen->SetNameTitle(strNamePtGen,"");
  hdNdptTracksMCRec->SetNameTitle(strNamePtRec,"");
  
  // corr fac
  TString strTitCorr        = "hBbBCorrSinglePt";
  if(useRecPrim) strTitCorr = "hBbBCorrRecPrimSinglePt";

  hCorrPt = (TH1F*) hdNdptTracksMCGen->Clone(strTitCorr);
  hCorrPt->Divide(hdNdptTracksMCGen,hdNdptTracksMCRec,1,1,"B"); // binominal errors
  
  // write 

  TString outfileOption = "RECREATE";
  if(updateOutfile)  outfileOption = "UPDATE";

  TFile out(strOutfile,outfileOption);
  
  if(!out.IsOpen()){
    Printf("%s:%d -- error opening efficiency output file %s", (char*)__FILE__,__LINE__,strOutfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- write bin shift correction to file %s ",(char*)__FILE__,__LINE__,strOutfile.Data());

  if(strOutDir && strOutDir.Length()){  
    
    TDirectory* dir;
    if((dir = ((TDirectory*) gDirectory->Get(strOutDir)))) dir->cd(); 
    else{
      dir = out.mkdir(strOutDir);
      dir->cd(); 
    } 
  }
  
  
  hdNdptTracksMCGen->Write();
  hdNdptTracksMCRec->Write();
  hCorrPt->Write();
  
  out.Close();
 
  delete hdNdptTracksMCGen; 
  delete hdNdptTracksMCRec;  
  delete hCorrPt; 
}

//___________________________________________________________________________________________________________________________________
void AliFragmentationFunctionCorrections::ReadBinShiftCorrSinglePt(TString strfile, TString strdir, TString strlist, Bool_t useRecPrim)
{
  // read reconstruction efficiency from file
  // argument strlist optional - read from directory strdir if not specified

  // temporary histos to hold histos from file
  TH1F* hBbBCorrPt = 0; 
  
  TFile f(strfile,"READ");

  if(!f.IsOpen()){
    Printf("%s:%d -- error opening raw data file %s", (char*)__FILE__,__LINE__,strfile.Data());
    return;
  }

  if(fDebug>0) Printf("%s:%d -- read BbB corr from file %s ",(char*)__FILE__,__LINE__,strfile.Data());
 
  if(strdir && strdir.Length()) gDirectory->cd(strdir);

  TList* list = 0;

  if(strlist && strlist.Length()){
   
    if(!(list = (TList*) gDirectory->Get(strlist))){ 
      Printf("%s:%d -- error retrieving list %s from directory %s", (char*)__FILE__,__LINE__,strlist.Data(),strdir.Data());
      return;
    }
  }  

    
  TString strNameBbBCorrPt = "hBbBCorrSinglePt";
  if(useRecPrim) strNameBbBCorrPt = "hBbBCorrRecPrimSinglePt";     

  if(list){
    hBbBCorrPt = (TH1F*)  list->FindObject(strNameBbBCorrPt); 
  }
  else{
    hBbBCorrPt = (TH1F*)  gDirectory->Get(strNameBbBCorrPt); 
  }
    
  if(!hBbBCorrPt){
    Printf("%s:%d -- error retrieving BbB corr single pt %s", (char*)__FILE__,__LINE__,strNameBbBCorrPt.Data());
  }
 

  if(fNHistoBinsPt && hBbBCorrPt) 
    hBbBCorrPt = (TH1F*) hBbBCorrPt->Rebin(fNHistoBinsSinglePt,strNameBbBCorrPt+"_rebin",fHistoBinsSinglePt->GetArray());

  if(hBbBCorrPt) hBbBCorrPt->SetDirectory(0); 
 
  f.Close();

  fh1BbBCorrSinglePt = hBbBCorrPt;
}

// ------------------------------------------------------------------

void AliFragmentationFunctionCorrections::BbBCorrSinglePt()
{
  // apply efficiency correction to inclusive track pt spec

  AddCorrectionLevelSinglePt("BbB");

  TH1F* histPt = fCorrSinglePt[fNCorrectionLevelsSinglePt-2]->GetTrackPt(0);
  
  if(histPt->GetNbinsX() != fh1BbBCorrSinglePt->GetNbinsX()){
    Printf("%s:%d: inconsistency pt spec and BbB corr binning ", (char*)__FILE__,__LINE__);
    return; 
  }

  TString histNamePt    = histPt->GetName();
  TH1F* hTrackPtBbBCorr = (TH1F*) histPt->Clone(histNamePt);
  
  hTrackPtBbBCorr->Multiply(histPt,fh1BbBCorrSinglePt,1,1,"");
  
  fCorrSinglePt[fNCorrectionLevelsSinglePt-1]->AddCorrHistos(0,hTrackPtBbBCorr);
}