//
// Class to do the energy correction of FMD ESD data
//
#include "AliFMDEnergyFitter.h"
#include "AliForwardUtil.h"
#include "AliLandauGausFitter.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TMath.h>
#include <AliLog.h>
#include <TClonesArray.h>
#include <TFitResult.h>
#include <THStack.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>

ClassImp(AliFMDEnergyFitter)
#if 0
; // This is for Emacs
#endif 
namespace {
  const char* fgkEDistFormat = "%s_etabin%03d";
}


//____________________________________________________________________
AliFMDEnergyFitter::AliFMDEnergyFitter()
  : TNamed(), 
    fRingHistos(),
    fLowCut(0.4),
    fNParticles(5),
    fMinEntries(1000),
    fFitRangeBinWidth(4),
    fDoFits(false),
    fDoMakeObject(false),
    fEtaAxis(),
    fCentralityAxis(),
    fMaxE(10),
    fNEbins(300), 
    fUseIncreasingBins(true),
    fMaxRelParError(AliFMDCorrELossFit::ELossFit::fgMaxRelError),
    fMaxChi2PerNDF(AliFMDCorrELossFit::ELossFit::fgMaxChi2nu), 
    fMinWeight(AliFMDCorrELossFit::ELossFit::fgLeastWeight),
    fDebug(0),
    fResidualMethod(kNoResiduals),
    fSkips(0),
    fRegularizationCut(3e6)
{
  // 
  // Default Constructor - do not use 
  //
  // DGUARD(fDebug, 1, "Default CTOR of AliFMDEnergyFitter");
  fRingHistos.SetOwner();
}

//____________________________________________________________________
AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
  : TNamed("fmdEnergyFitter", title), 
    fRingHistos(), 
    fLowCut(0.4),
    fNParticles(5),
    fMinEntries(1000),
    fFitRangeBinWidth(4),
    fDoFits(false),
    fDoMakeObject(false),
    fEtaAxis(0,0,0),
    fCentralityAxis(0,0,0),
    fMaxE(10),
    fNEbins(300), 
    fUseIncreasingBins(true),
    fMaxRelParError(AliFMDCorrELossFit::ELossFit::fgMaxRelError),
    fMaxChi2PerNDF(AliFMDCorrELossFit::ELossFit::fgMaxChi2nu), 
    fMinWeight(AliFMDCorrELossFit::ELossFit::fgLeastWeight),
    fDebug(3),
    fResidualMethod(kNoResiduals),
    fSkips(0),
    fRegularizationCut(3e6)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    title Title of object  - not significant 
  //
  // DGUARD(fDebug, 1, "Named CTOR of AliFMDEnergyFitter: %s", title);
  fEtaAxis.SetName("etaAxis");
  fEtaAxis.SetTitle("#eta");
  fCentralityAxis.SetName("centralityAxis");
  fCentralityAxis.SetTitle("Centrality [%]");
}

//____________________________________________________________________
AliFMDEnergyFitter::~AliFMDEnergyFitter()
{
  // 
  // Destructor
  //
  DGUARD(fDebug, 1, "DTOR of AliFMDEnergyFitter");
  // fRingHistos.Delete();
}

//____________________________________________________________________
AliFMDEnergyFitter::RingHistos*
AliFMDEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const
{
  return new RingHistos(d,r);
}

//____________________________________________________________________
AliFMDEnergyFitter::RingHistos*
AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
{
  // 
  // Get the ring histogram container 
  // 
  // Parameters:
  //    d Detector
  //    r Ring 
  // 
  // Return:
  //    Ring histogram container 
  //
  Int_t idx = -1;
  switch (d) { 
  case 1: idx = 0; break;
  case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
  case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
  }
  if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
  
  return static_cast<RingHistos*>(fRingHistos.At(idx));
}

//____________________________________________________________________
void
AliFMDEnergyFitter::Init()
{
  // Create the ring histograms.  
  // 
  // Should be called from task initialization. 
  DGUARD(1,fDebug, "Creating histogram caches for each ring");
  fRingHistos.Add(CreateRingHistos(1, 'I'));
  fRingHistos.Add(CreateRingHistos(2, 'I'));
  fRingHistos.Add(CreateRingHistos(2, 'O'));
  fRingHistos.Add(CreateRingHistos(3, 'I'));
  fRingHistos.Add(CreateRingHistos(3, 'O'));
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->fDebug = fDebug;
  }
}

//____________________________________________________________________
void
AliFMDEnergyFitter::CreateOutputObjects(TList* dir)
{
  // 
  // Define the output histograms.  These are put in a sub list of the
  // passed list.   The histograms are merged before the parent task calls 
  // AliAnalysisTaskSE::Terminate.  Called on task initialization on slaves.  
  // 
  // Parameters:
  //    dir Directory to add to 
  //
  DGUARD(fDebug, 1, "Define output in AliFMDEnergyFitter");
  TList* d = new TList;
  d->SetName(GetName());
  d->SetOwner(true);
  dir->Add(d);

  // Store eta axis as a histogram, since that can be merged!
  TH1* hEta = 0;
  if (fEtaAxis.GetXbins()->GetArray()) 
    hEta = new TH1I(fEtaAxis.GetName(), fEtaAxis.GetTitle(), 
		    fEtaAxis.GetNbins(), fEtaAxis.GetXbins()->GetArray());
  else 
    hEta = new TH1I(fEtaAxis.GetName(), fEtaAxis.GetTitle(), 
		    fEtaAxis.GetNbins(),fEtaAxis.GetXmin(),fEtaAxis.GetXmax());
  hEta->SetXTitle("#eta");
  hEta->SetYTitle("Nothing");
  hEta->SetDirectory(0);

  d->Add(hEta);
  d->Add(AliForwardUtil::MakeParameter("lowCut",        fLowCut));
  d->Add(AliForwardUtil::MakeParameter("nParticles",    fNParticles));
  d->Add(AliForwardUtil::MakeParameter("minEntries",    fMinEntries));
  d->Add(AliForwardUtil::MakeParameter("subtractBins",  fFitRangeBinWidth));
  d->Add(AliForwardUtil::MakeParameter("doFits",        fDoFits));
  d->Add(AliForwardUtil::MakeParameter("doObject",      fDoMakeObject));
  d->Add(AliForwardUtil::MakeParameter("maxE",          fMaxE));
  d->Add(AliForwardUtil::MakeParameter("nEbins",        fNEbins));
  d->Add(AliForwardUtil::MakeParameter("increasingBins",fUseIncreasingBins));
  d->Add(AliForwardUtil::MakeParameter("maxRelPerError",fMaxRelParError));
  d->Add(AliForwardUtil::MakeParameter("maxChi2PerNDF", fMaxChi2PerNDF));
  d->Add(AliForwardUtil::MakeParameter("minWeight",     fMinWeight));
  d->Add(AliForwardUtil::MakeParameter("regCut",        fRegularizationCut));
  d->Add(AliForwardUtil::MakeParameter("deltaShift", 
				       AliLandauGaus::EnableSigmaShift()));

  if (fRingHistos.GetEntries() <= 0) { 
    AliFatal("No ring histograms where defined - giving up!");
    return;
  }
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->fDebug = fDebug;
    o->CreateOutputObjects(d);
  }
}

//____________________________________________________________________
void
AliFMDEnergyFitter::SetupForData(const TAxis& eAxis)
{
  // 
  // Initialise the task - called at first event 
  // 
  // Parameters:
  //    etaAxis The eta axis to use.  Note, that if the eta axis
  // has already been set (using SetEtaAxis), then this parameter will be 
  // ignored
  //
  DGUARD(fDebug, 1, "Initialize of AliFMDEnergyFitter");
  if (fEtaAxis.GetNbins() == 0 || 
      TMath::Abs(fEtaAxis.GetXmax() - fEtaAxis.GetXmin()) < 1e-7) 
    SetEtaAxis(eAxis);
  if (fCentralityAxis.GetNbins() == 0) {
    UShort_t n = 12;
    Double_t bins[] = {  0.,  5., 10., 15., 20., 30., 
			 40., 50., 60., 70., 80., 100. };
    SetCentralityAxis(n, bins);
  }
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next())))
    o->SetupForData(fEtaAxis, fCentralityAxis, fMaxE, 
		    fNEbins, fUseIncreasingBins);
}  
//____________________________________________________________________
void
AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
{
  // 
  // Set the eta axis to use.  This will force the code to use this
  // eta axis definition - irrespective of whatever axis is passed to
  // the Init member function.  Therefore, this member function can be
  // used to force another eta axis than one found in the correction
  // objects. 
  // 
  // Parameters:
  //    etaAxis Eta axis to use 
  //
  SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
}
//____________________________________________________________________
void
AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
{
  // 
  // Set the eta axis to use.  This will force the code to use this
  // eta axis definition - irrespective of whatever axis is passed to
  // the Init member function.  Therefore, this member function can be
  // used to force another eta axis than one found in the correction
  // objects. 
  // 
  // Parameters:
  //    nBins  Number of bins 
  //    etaMin Minimum of the eta axis 
  //    etaMax Maximum of the eta axis 
  //
  fEtaAxis.Set(nBins,etaMin,etaMax);
}
//____________________________________________________________________
void
AliFMDEnergyFitter::SetCentralityAxis(UShort_t n, Double_t* bins)
{
  fCentralityAxis.Set(n-1, bins);
}
//____________________________________________________________________
void
AliFMDEnergyFitter::SetEnableDeltaShift(Bool_t use) 
{
  AliLandauGaus::EnableSigmaShift(use ? 1 : 0);
}

//____________________________________________________________________
Bool_t
AliFMDEnergyFitter::Accumulate(const AliESDFMD& input,
			       Double_t         cent,
			       Bool_t           empty)
{
  // 
  // Fitter the input AliESDFMD object
  // 
  // Parameters:
  //    input     Input 
  //    cent      Centrality 
  //    empty     Whether the event is 'empty'
  // 
  // Return:
  //    True on success, false otherwise 
  DGUARD(fDebug, 5, "Accumulate statistics in AliFMDEnergyFitter - cholm");
  Int_t icent = fCentralityAxis.FindBin(cent);
  if (icent < 1 || icent > fCentralityAxis.GetNbins()) icent = 0;

  UShort_t nFills = 0;
  for(UShort_t d = 1; d <= 3; d++) {
    Int_t nRings = (d == 1 ? 1 : 2);
    for (UShort_t q = 0; q < nRings; q++) {
      Char_t      r    = (q == 0 ? 'I' : 'O');
      UShort_t    nsec = (q == 0 ?  20 :  40);
      UShort_t    nstr = (q == 0 ? 512 : 256);

      RingHistos* histos = GetRingHistos(d, r);
      if (!histos) continue;
      
      for(UShort_t s = 0; s < nsec;  s++) {
	for(UShort_t t = 0; t < nstr; t++) {
	  Float_t mult = input.Multiplicity(d,r,s,t);
	  
	  // Keep dead-channel information. 
	  if (mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0) 
	    continue;

	  // Get the pseudo-rapidity 
	  Double_t eta1 = input.Eta(d,r,s,t);

	  // Int_t ieta = fEtaAxis.FindBin(eta1);
	  // if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 

	  // histos->Fill(empty, ieta-1, icent, mult);
	  histos->Fill(empty, eta1, icent, mult);
	  nFills++;
	} // for strip
      } // for sector
    } // for ring 
  } // for detector

  DMSG(fDebug, 3, "Found a total of %d signals for c=%f, and %sempty event", 
       nFills, cent, (empty ? "" : "non-"));
  return kTRUE;
}

//____________________________________________________________________
void
AliFMDEnergyFitter::Fit(const TList* dir)
{
  // 
  // Scale the histograms to the total number of events 
  // 
  // Parameters:
  //    dir Where the histograms are  
  //
  DGUARD(fDebug, 1, "Fit distributions in AliFMDEnergyFitter");
  if (!fDoFits) {
    AliInfo("Not asked to do fits, returning");
    return;
  }

  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
  if (!d) {
    AliWarningF("No list named %s found in %s", GetName(), dir->GetName());
    return;
  }    

  // +1          for chi^2
  // +3          for 1 landau 
  // +1          for N 
  // +fNParticles-1 for weights 
  Int_t nStack = kN+fNParticles;
  THStack* stack[nStack]; 
  stack[0]         = new THStack("chi2",   "#chi^{2}/#nu");
  stack[kC     +1] = new THStack("c",      "Constant");
  stack[kDelta +1] = new THStack("delta",  "#Delta_{p}");
  stack[kXi    +1] = new THStack("xi",     "#xi");
  stack[kSigma +1] = new THStack("sigma",  "#sigma");
  stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}");
  stack[kN     +1] = new THStack("n",      "# Particles");
  for (Int_t i = 2; i <= fNParticles; i++) 
    stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i));
  for (Int_t i = 0; i < nStack; i++) 
    d->Add(stack[i]);

  // If we have no ring histograms, re-init. 
  if (fRingHistos.GetEntries() <= 0) Init();

  AliInfoF("Will do fits for %d rings", fRingHistos.GetEntries());
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    AliInfoF("Fill fit for FMD%d%c", o->fDet, o->fRing);
    if (CheckSkip(o->fDet, o->fRing, fSkips)) {
      AliWarningF("Skipping FMD%d%c for fitting", o->fDet, o->fRing);
      continue;
    }
    
    TObjArray* l = o->Fit(d, fLowCut, fNParticles,
			  fMinEntries, fFitRangeBinWidth,
			  fMaxRelParError, fMaxChi2PerNDF,
			  fMinWeight, fRegularizationCut,
			  fResidualMethod);
    if (!l) continue;
    for (Int_t i = 0; i < l->GetEntriesFast()-1; i++) { // Last is status 
      stack[i % nStack]->Add(static_cast<TH1*>(l->At(i))); 
    }
  }

  if (!fDoMakeObject) return;

  MakeCorrectionsObject(d);
}

//____________________________________________________________________
void
AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
{
  // 
  // Generate the corrections object 
  // 
  // Parameters:
  //    dir List to analyse 
  //
  DGUARD(fDebug, 1, "Make the correction objec in AliFMDEnergyFitter");
    
  AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
  obj->SetEtaAxis(fEtaAxis);
  obj->SetLowCut(fLowCut);
  if (AliLandauGaus::EnableSigmaShift()) 
    obj->SetBit(AliFMDCorrELossFit::kHasShift);

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    if (CheckSkip(o->fDet, o->fRing, fSkips)) {
      AliWarningF("Skipping FMD%d%c for correction object", o->fDet, o->fRing);
      continue;
    }
    
    o->FindBestFits(d, *obj, fEtaAxis);
  }
  d->Add(obj, "elossFits");
}

//____________________________________________________________________
void
AliFMDEnergyFitter::SetDebug(Int_t dbg)
{
  // 
  // Set the debug level.  The higher the value the more output 
  // 
  // Parameters:
  //    dbg Debug level 
  //
  fDebug = dbg;
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next())))
    o->fDebug = dbg;
}  

//____________________________________________________________________
namespace {
  template <typename T>
  void GetParam(Bool_t& ret, const TCollection* col, 
		const TString& name, T& var)
  {
    TObject* o = col->FindObject(name);	      
    if (o) AliForwardUtil::GetParameter(o,var); 
    else   ret = false;				
  }
}

//____________________________________________________________________
Bool_t
AliFMDEnergyFitter::ReadParameters(const TCollection* col)
{
  // Read parameters of this object from a collection
  //
  // Parameters:
  //    col   Collection to read parameters from 
  // 
  // Return value:
  //   true on success, false otherwise 
  //
  if (!col) return false;
  Bool_t ret  = true;
  TH1*   hist = static_cast<TH1*>(col->FindObject("etaAxis"));
  if (!hist) ret = false;
  else {
    TAxis* axis = hist->GetXaxis();
    if (axis->GetXbins()->GetArray()) 
      fEtaAxis.Set(axis->GetNbins(), axis->GetXbins()->GetArray());
   else 
      fEtaAxis.Set(axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
  } 
  GetParam(ret,col,"lowCut",        fLowCut);
  GetParam(ret,col,"nParticles",    fNParticles);
  GetParam(ret,col,"minEntries",    fMinEntries);
  GetParam(ret,col,"subtractBins",  fFitRangeBinWidth);
  GetParam(ret,col,"doFits",        fDoFits);
  GetParam(ret,col,"doObject",      fDoMakeObject);
  GetParam(ret,col,"maxE",          fMaxE);
  GetParam(ret,col,"nEbins",        fNEbins);
  GetParam(ret,col,"increasingBins",fUseIncreasingBins);
  GetParam(ret,col,"maxRelPerError",fMaxRelParError);
  GetParam(ret,col,"maxChi2PerNDF", fMaxChi2PerNDF);
  GetParam(ret,col,"minWeight",     fMinWeight);
  Bool_t dummy;
  GetParam(dummy,col,"regCut",      fRegularizationCut);

  return ret;
}

//____________________________________________________________________
Bool_t
AliFMDEnergyFitter::CheckSkip(UShort_t d, Char_t r, UShort_t skips) 
{
  UShort_t q  = (r == 'I' || r == 'i' ? 0 : 1);
  UShort_t c = 1 << (d-1);
  UShort_t t = 1 << (c+q-1);

  return (t & skips) == t;
}

#define PF(N,V,...)					\
  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
#define PFB(N,FLAG)				\
  do {									\
    AliForwardUtil::PrintName(N);					\
    std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
  } while(false)
#define PFV(N,VALUE)					\
  do {							\
    AliForwardUtil::PrintName(N);			\
    std::cout << (VALUE) << std::endl; } while(false)

//____________________________________________________________________
void
AliFMDEnergyFitter::Print(Option_t*) const
{
  // 
  // Print information
  // 
  // Parameters:
  //    option Not used 
  //
  AliForwardUtil::PrintTask(*this);

  gROOT->IncreaseDirLevel();
  PFV("Low cut [E/E_mip]",	fLowCut);
  PFV("Max(particles)",	        fNParticles);
  PFV("Min(entries)",	        fMinEntries);
  PFV("Fit range [bins]",       fFitRangeBinWidth);
  PFB("Make fits",              fDoFits);
  PFB("Make object",            fDoMakeObject);
  PFV("Max E/E_mip",	        fMaxE);
  PFV("N bins",	                fNEbins);
  PFB("Increasing bins",        fUseIncreasingBins);
  PFV("max(delta p/p)",    	fMaxRelParError);
  PFV("max(chi^2/nu)",	        fMaxChi2PerNDF);
  PFV("min(a_i)",	        fMinWeight);
  PFV("Regularization cut",     fRegularizationCut);
  TString r = "";
  switch (fResidualMethod) { 
  case kNoResiduals:              r = "None";       break;
  case kResidualDifference:       r = "Difference"; break;
  case kResidualScaledDifference: r = "Scaled difference"; break;
  case kResidualSquareDifference: r = "Square difference"; break;
  }
  PFV("Residuals", r);
  gROOT->DecreaseDirLevel();
}
  
//====================================================================
AliFMDEnergyFitter::RingHistos::RingHistos()
  : AliForwardUtil::RingHistos(), 
    fEDist(0), 
    fEmpty(0),
    // fEtaEDists(0), 
    fHist(0),
    fList(0),
    fBest(0),
    fFits("AliFMDCorrELossFit::ELossFit", 200),
    fDebug(0)
{
  // 
  // Default CTOR
  //
  DGUARD(fDebug, 3, "Default CTOR AliFMDEnergyFitter::RingHistos");
  fBest.Expand(0);
}

//____________________________________________________________________
AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliForwardUtil::RingHistos(d,r), 
    fEDist(0), 
    fEmpty(0),
    // fEtaEDists(0), 
    fHist(0),
    fList(0),
    fBest(0),
    fFits("AliFMDCorrELossFit::ELossFit", 200),
    fDebug(0)
{
  // 
  // Constructor
  // 
  // Parameters:
  //    d detector
  //    r ring 
  //
  DGUARD(fDebug, 3, "Named CTOR AliFMDEnergyFitter::RingHistos: FMD%d%c",
	 d, r);
  fBest.Expand(0);
}
//____________________________________________________________________
AliFMDEnergyFitter::RingHistos::~RingHistos()
{
  // 
  // Destructor 
  //
  DGUARD(fDebug, 3, "DTOR of AliFMDEnergyFitter::RingHistos");
  // if (fEDist) delete fEDist;
  // fEtaEDists.Delete();
}
//__________________________________________________________________
TArrayD
AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, 
						   Double_t max) const
{
  // 
  // Make an axis with increasing bins 
  // 
  // Parameters:
  //    n    Number of bins 
  //    min  Minimum 
  //    max  Maximum
  // 
  // Return:
  //    An axis with quadratically increasing bin size 
  //

  TArrayD bins(n+1);
  Double_t dx = (max-min) / n;
  bins[0]     = min;
  Int_t    i  = 1;
  for (i = 1; i < n+1; i++) {
    Double_t dI   = float(i)/n;
    Double_t next = bins[i-1] + dx + dI * dI * dx;
    bins[i]       = next;
    if (next > max) break;
  }
  bins.Set(i+1);
  return bins;
}

//____________________________________________________________________
TH2*
AliFMDEnergyFitter::RingHistos::Make(const char*  name, 
				     const char*  title, 
				     const TAxis& eAxis, 
				     Double_t     deMax, 
				     Int_t        nDeBins, 
				     Bool_t       incr)
{
  // 
  // Make E/E_mip histogram 
  // 
  // Parameters:
  //    deMax   Maximum energy loss 
  //    nDeBins Number energy loss bins 
  //    incr    Whether to make bins of increasing size
  //
  TH2* h = 0;
  TAxis mAxis;
  if (incr) {
    TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
    mAxis.Set(deAxis.GetSize()-1, deAxis.GetArray());
  }
  else 
    mAxis.Set(nDeBins, 0, deMax);
  
  if (mAxis.GetXbins()->GetArray()) {
    // Variable bin since in Delta 
    if (eAxis.GetXbins()->GetArray()) {
      // Variadic bin size in eta 
      h = new TH2D(name, title, 
		   eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
		   mAxis.GetNbins(), mAxis.GetXbins()->GetArray());
    }
    else { 
      // Evenly spaced bins in eta
      h = new TH2D(name, title, 
		   eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), 
		   mAxis.GetNbins(), mAxis.GetXbins()->GetArray());
    }
  }
  else { 
    // Make increasing bins axis 
    if (eAxis.GetXbins()->GetArray()) {
      // Variable size eta bins 
      h = new TH2D(name, title, 
		   eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
		   mAxis.GetNbins(), mAxis.GetXmin(), mAxis.GetXmax());
    }
    else {
      // Fixed size eta bins 
      h = new TH2D(name, title, 
		   eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(), 
		   mAxis.GetNbins(), mAxis.GetXmin(), mAxis.GetXmax());
    }
  }
  h->SetDirectory(0);
  h->SetYTitle("#sum#DeltaE/#DeltaE_{mip}");	
  h->SetXTitle("#eta");	
  h->SetFillColor(Color());
  h->SetMarkerColor(Color());
  h->SetLineColor(Color());
  h->SetFillStyle(3001);
  h->SetMarkerStyle(20);
  h->Sumw2();

  return h;
}
//____________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::CreateOutputObjects(TList* dir)
{
  // 
  // Define outputs
  // 
  // Parameters:
  //    dir 
  //
  DGUARD(fDebug, 2, "Define output in AliFMDEnergyFitter::RingHistos");
  fList = DefineOutputList(dir);
  
  // fEtaEDists = new TList;
  // fEtaEDists->SetOwner();
  // fEtaEDists->SetName("EDists");
  // 
  // fList->Add(fEtaEDists);
}
//____________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis, 
					     const TAxis& /* cAxis */,
					     Double_t     maxDE, 
					     Int_t        nDEbins, 
					     Bool_t       useIncrBin)
{
  // 
  // Initialise object 
  // 
  // Parameters:
  //    eAxis      Eta axis
  //    maxDE      Max energy loss to consider 
  //    nDEbins    Number of bins 
  //    useIncrBin Whether to use an increasing bin size 
  //
  DGUARD(fDebug, 2, "Initialize in AliFMDEnergyFitter::RingHistos");
  fEDist = new TH1D(Form("%s_edist", fName.Data()), 
		    Form("#sum#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 
		    200, 0, 6);
  fEDist->SetXTitle("#sum#Delta/#Delta_{mip}");
  fEDist->SetFillColor(Color());
  fEDist->SetLineColor(Color());
  fEDist->SetMarkerColor(Color());
  fEDist->SetFillStyle(3001);
  fEDist->SetMarkerStyle(20);
  fEDist->Sumw2();
  fEDist->SetDirectory(0);

  fEmpty = static_cast<TH1D*>(fEDist->Clone(Form("%s_empty", fName.Data())));
  fEmpty->SetTitle(Form("#sum#DeltaE/#DeltaE_{mip} for %s (empty events)", 
			fName.Data()));
  fEmpty->SetDirectory(0);

  fList->Add(fEDist);
  fList->Add(fEmpty);
  fHist = Make("eloss", "#sum#Delta/#Delta_{mip}", eAxis, 
	       maxDE, nDEbins, useIncrBin);
  fList->Add(fHist);
  // fList->ls();
  // fEtaEDists.ls();
}

//____________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::Fill(Bool_t   empty, 
				     Double_t eta, 
				     Int_t /* icent */,  
				     Double_t mult)
{
  // 
  // Fill histogram 
  // 
  // Parameters:
  //    empty  True if event is empty
  //    ieta   Eta bin (0-based)
  //    icent  Centrality bin (1-based)
  //    mult   Signal 
  //
  DGUARD(fDebug, 10, "Filling in AliFMDEnergyFitter::RingHistos");
  if (empty) { 
    fEmpty->Fill(mult);
    return;
  }
  fEDist->Fill(mult);

  // if (!fEtaEDists) { 
  if (!fHist) {
    Warning("Fill", "No list of E dists defined");
    return;
  }
  fHist->Fill(eta, mult);
  // Here, we should first look up in a centrality array, and then in
  // that array look up the eta bin - or perhaps we should do it the
  // other way around?
  // if (ieta < 0 || ieta >= fEtaEDists->GetEntries()) return;
  
  // TH1D* h = static_cast<TH1D*>(fEtaEDists->At(ieta));
  // if (!h) return;

  // h->Fill(mult);
}

//____________________________________________________________________
TH1* 
AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
					const char* title, 
					const TAxis& eta) const
{
  // 
  // Make a parameter histogram
  // 
  // Parameters:
  //    name   Name of histogram.
  //    title  Title of histogram. 
  //    eta    Eta axis 
  // 
  // Return:
  //    
  //
  TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
		     Form("%s for %s", title, fName.Data()), 
		     eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
  h->SetXTitle("#eta");
  h->SetYTitle(title);
  h->SetDirectory(0);
  h->SetFillColor(Color());
  h->SetMarkerColor(Color());
  h->SetLineColor(Color());
  h->SetFillStyle(3001);
  h->Sumw2();

  return h;
}
//____________________________________________________________________
TH1*
AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, 
					  const char* title, 
					  const TAxis& eta, 
					  Int_t low, 
					  Int_t high, 
					  Double_t val, 
					  Double_t err) const
{
  // 
  // Make a histogram that contains the results of the fit over the full ring 
  // 
  // Parameters:
  //    name  Name 
  //    title Title
  //    eta   Eta axis 
  //    low   Least bin
  //    high  Largest bin
  //    val   Value of parameter 
  //    err   Error on parameter 
  // 
  // Return:
  //    The newly allocated histogram 
  //
  Double_t xlow  = eta.GetBinLowEdge(low);
  Double_t xhigh = eta.GetBinUpEdge(high);
  TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
		     Form("%s for %s", title, fName.Data()), 
		     1, xlow, xhigh);
  h->SetBinContent(1, val);
  h->SetBinError(1, err);
  h->SetXTitle("#eta");
  h->SetYTitle(title);
  h->SetDirectory(0);
  h->SetFillColor(0);
  h->SetMarkerColor(Color());
  h->SetLineColor(Color());
  h->SetFillStyle(0);
  h->SetLineStyle(2);
  h->SetLineWidth(2);

  return h;
}
		     
//____________________________________________________________________
TObjArray*
AliFMDEnergyFitter::RingHistos::Fit(TList*           dir, 
				    Double_t         lowCut, 
				    UShort_t         nParticles,
				    UShort_t         minEntries,
				    UShort_t         minusBins, 
				    Double_t         relErrorCut, 
				    Double_t         chi2nuCut,
				    Double_t         minWeight,
				    Double_t         regCut,
				    EResidualMethod  residuals) const
{
  return FitSlices(dir, "eloss", lowCut, nParticles, minEntries, 
		   minusBins, relErrorCut, chi2nuCut, minWeight, regCut, 
		   residuals, true, &fBest);
}

//____________________________________________________________________
TObjArray*
AliFMDEnergyFitter::RingHistos::FitSlices(TList*           dir, 
					  const char*      name, 
					  Double_t         lowCut, 
					  UShort_t         nParticles,
					  UShort_t         minEntries,
					  UShort_t         minusBins, 
					  Double_t         relErrorCut, 
					  Double_t         chi2nuCut,
					  Double_t         minWeight,
					  Double_t         regCut,
					  EResidualMethod  residuals,
					  Bool_t           scaleToPeak,
					  TObjArray*       best) const
{
  // 
  // Fit each histogram to up to @a nParticles particle responses.
  // 
  // Parameters:
  //    dir         Output list 
  //    eta         Eta axis 
  //    lowCut      Lower cut 
  //    nParticles  Max number of convolved landaus to fit
  //    minEntries  Minimum number of entries 
  //    minusBins   Number of bins from peak to subtract to 
  //                    get the fit range 
  //    relErrorCut Cut applied to relative error of parameter. 
  //                    Note, for multi-particle weights, the cut 
  //                    is loosend by a factor of 2 
  //    chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
  //                    the reduced @f$\chi^2@f$ 
  //
  DGUARD(fDebug, 2, "Fit in AliFMDEnergyFitter::RingHistos");

  // Get our ring list from the output container 
  TList* l = GetOutputList(dir);
  if (!l) return 0; 

  TList* dists = 0;
  // Get the 2D histogram 
  TH2* h = static_cast<TH2*>(l->FindObject(name));
  if (!h) { 
    AliWarningF("Didn't find 2D histogram '%s' in %s", name, l->GetName());
    // Get the energy distributions from the output container 
    dists = static_cast<TList*>(l->FindObject("EDists"));
    if (!dists) { 
      AliWarningF("Didn't find EtaEDists (%s) in %s", 
		  fName.Data(), l->GetName());
      l->ls();
      return 0;
    }
  }
  if (!h && !dists) return 0;

  const TAxis* pEta = (h ? h->GetXaxis() : 
			    static_cast<TAxis*>(dir->FindObject("etaAxis")));
  if (!pEta) { 
    AliWarningF("Didn't find the eta axis - either from histogram %p or "
		"list %p (%s)", h, dir, (dir ? dir->GetName() : "-"));
    return 0;
  }
  const TAxis& eta = *pEta;

  // Create an output list for the fitted distributions 
  TList* out = new TList;
  out->SetOwner();
  out->SetName(Form("%sDists", name));
  l->Add(out);

  // Optional container for residuals 
  TList* resi = 0;
  if (residuals != kNoResiduals) {
    resi = new TList();
    resi->SetName(Form("%sResiduals", name));
    resi->SetOwner();
    l->Add(resi);
  }

  // Container of the fit results histograms 
  TObjArray* pars  = new TObjArray(3+nParticles+1);
  pars->SetName(Form("%sResults", name));
  l->Add(pars);

  // Result objects stored in separate list on the output 
  TH1* hChi2   = 0;
  TH1* hC      = 0;
  TH1* hDelta  = 0;
  TH1* hXi     = 0;
  TH1* hSigma  = 0;
  TH1* hSigmaN = 0;
  TH1* hN      = 0;
  TH1* hA[nParticles-1];
  pars->Add(hChi2   = MakePar("chi2",   "#chi^{2}/#nu", eta));
  pars->Add(hC      = MakePar("c",      "Constant",     eta));
  pars->Add(hDelta  = MakePar("delta",  "#Delta_{p}",   eta));
  pars->Add(hXi     = MakePar("xi",     "#xi",          eta));
  pars->Add(hSigma  = MakePar("sigma",  "#sigma",       eta));
  pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}",   eta));
  pars->Add(hN      = MakePar("n",      "N", eta));
  for (UShort_t i = 1; i < nParticles; i++) 
    pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));

  
  Int_t nDists = h ? h->GetNbinsX() : dists->GetEntries();
  Int_t low    = nDists;
  Int_t high   = 0;
  Int_t nEmpty = 0;
  Int_t nLow   = 0;
  Int_t nFitted= 0;
  if (best) {
    best->Expand(nDists+1);
    best->Clear();
    best->SetOwner(false);
  }
  for (Int_t i = 0; i < nDists; i++) { 
    // Ignore empty histograms altoghether 
    Int_t b    = i+1;
    TH1D* dist = (h ? h->ProjectionY(Form(fgkEDistFormat,GetName(),b),b,b,"e") 
		  : static_cast<TH1D*>(dists->At(i)));
    if (!dist) { 
      // If we got the null pointer, return 0
      nEmpty++;
      continue;
    }
    // Then releasing the histogram from the it's directory
    dist->SetDirectory(0);
    // Set a meaningful title
    dist->SetTitle(Form("#Delta/#Delta_{mip} for %s in %6.2f<#eta<%6.2f",
			GetName(), eta.GetBinLowEdge(b),
			eta.GetBinUpEdge(b)));

    // Now fit 
    UShort_t    status1 = 0;
    ELossFit_t* res     = FitHist(dist,
				  lowCut, 
				  nParticles,
				  minEntries,
				  minusBins,   
				  relErrorCut,
				  chi2nuCut,
				  minWeight,
				  regCut,
				  scaleToPeak,
				  status1);
    if (!res) {
      switch (status1) { 
      case 1: nEmpty++; break;
      case 2: nLow++;   break;
      }
      // Only clean up if we have no input list 
      if (h) delete dist;
      continue;
    }
      
    // Now count as fitted, store as best fits, and add distribution
    // to the dedicated output list
    nFitted++;
    res->fBin = b; // We only have the bin information here 
    if (best) best->AddAt(res, b);
    out->Add(dist);
    // dist->GetListOfFunctions()->Add(res);
    
    // If asked to calculate residuals, do so, and store result on the
    // dedicated output list
    if (residuals != kNoResiduals && resi) 
      CalculateResiduals(residuals, lowCut, dist, res, resi);

    // Store eta limits 
    low   = TMath::Min(low,b);
    high  = TMath::Max(high,b);

    // Get the reduced chi square
    Double_t chi2 = res->fChi2; // GetChisquare();
    Int_t    ndf  = res->fNu;   // GetNDF();
    
    // Store results of best fit in output histograms 
    // res->SetLineWidth(4);
    hChi2   ->SetBinContent(b, ndf > 0 ? chi2 / ndf : 0);
    hC      ->SetBinContent(b, res->fC); 
    hDelta  ->SetBinContent(b, res->fDelta); 
    hXi     ->SetBinContent(b, res->fXi); 
    hSigma  ->SetBinContent(b, res->fSigma); 
    hSigmaN ->SetBinContent(b, res->fSigmaN); 
    hN      ->SetBinContent(b, res->fN); 

    hC     ->SetBinError(b, res->fEC); 
    hDelta ->SetBinError(b, res->fEDelta);
    hXi    ->SetBinError(b, res->fEXi); 
    hSigma ->SetBinError(b, res->fESigma); 
    hSigmaN->SetBinError(b, res->fESigmaN); 
    // hN     ->SetBinError(b, res->fGetParError(kN));

    for (UShort_t j = 0; j < nParticles-1; j++) {
      hA[j]->SetBinContent(b, res->fA[j]); 
      hA[j]->SetBinError(b, res->fEA[j]); 
    }
  }
  printf("%s: Out of %d histograms, %d where empty, %d had too little data,"
	 "leaving %d to be fitted, of which %d succeeded\n",  
	 GetName(), nDists, nEmpty, nLow, nDists-nEmpty-nLow, nFitted);

  // Fit the full-ring histogram 
  TH1*        total   = GetOutputHist(l, Form("%s_edist", fName.Data()));
  if (total) {
    UShort_t    statusT = 0;
    ELossFit_t* resT    = FitHist(total,
				  lowCut, 
				  nParticles,
				  minEntries, 
				  minusBins,
				  relErrorCut,
				  chi2nuCut,
				  minWeight,
				  regCut,
				  scaleToPeak,
				  statusT);
    if (resT) { 
      // Make histograms for the result of this fit 
      Double_t chi2 = resT->GetChi2();
      Int_t    ndf  = resT->GetNu();
      pars->Add(MakeTotal("t_chi2",   "#chi^{2}/#nu", eta, low, high,
			  ndf > 0 ? chi2/ndf : 0, 0));
      pars->Add(MakeTotal("t_c",      "Constant",     eta, low, high,
			  resT->GetC(), resT->GetEC()));
      pars->Add(MakeTotal("t_delta",  "#Delta_{p}",   eta, low, high,
			  resT->GetDelta(), resT->GetEDelta()));
      pars->Add(MakeTotal("t_xi",     "#xi",          eta, low, high,
			  resT->GetXi(), resT->GetEXi()));
      pars->Add(MakeTotal("t_sigma",  "#sigma",       eta, low, high,
			  resT->GetSigma(), resT->GetESigma()));
      pars->Add(MakeTotal("t_sigman", "#sigma_{n}",   eta, low, high,
			  resT->GetSigmaN(), resT->GetESigmaN()));
      pars->Add(MakeTotal("t_n",    "N",              eta, low, high,
			  resT->GetN(), 0));
      for (UShort_t j = 0; j < nParticles-1; j++) 
	pars->Add(MakeTotal(Form("t_a%d",j+2), 
			    Form("a_{%d}",j+2), eta, low, high,
			    resT->GetA(j), resT->GetEA(j)));
    }
  }

  TH1* status = new TH1I(Form("%sStatus",name), "Status of Fits", 5, 0, 5);
  status->GetXaxis()->SetBinLabel(1, "Total");
  status->GetXaxis()->SetBinLabel(2, "Empty");
  status->GetXaxis()->SetBinLabel(3, Form("<%d", minEntries));
  status->GetXaxis()->SetBinLabel(4, "Candidates");
  status->GetXaxis()->SetBinLabel(5, "Fitted");
  status->SetXTitle("Status");
  status->SetYTitle("# of #Delta distributions");
  status->SetBinContent(1, nDists);
  status->SetBinContent(2, nEmpty);
  status->SetBinContent(3, nLow);
  status->SetBinContent(4, nDists-nLow-nEmpty);
  status->SetBinContent(5, nFitted);
  status->SetFillColor(Color());
  status->SetFillStyle(3001);
  status->SetLineColor(Color());
  status->SetDirectory(0);
  status->SetStats(0);
  pars->Add(status);
  return pars;
}


//____________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::Scale(TH1* dist) const
{
  // Scale to the bin-width
  dist->Scale(1., "width");
}  


//____________________________________________________________________
AliFMDEnergyFitter::RingHistos::ELossFit_t*
AliFMDEnergyFitter::RingHistos::FitHist(TH1*      dist,
					Double_t  lowCut, 
					UShort_t  nParticles, 
					UShort_t  minEntries,
					UShort_t  minusBins, 
					Double_t  relErrorCut, 
					Double_t  chi2nuCut,
					Double_t  minWeight,
					Double_t  regCut,
					Bool_t    scaleToPeak,
					UShort_t& status) const
{
  // 
  // Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
  // maximum bin content in the range @f$ [E_{min},\infty]@f$ is
  // found.  Then the fit range is set to the bin range 
  // @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1 
  // particle signal is fitted to that.  The parameters of that fit 
  // is then used as seeds for a fit of the @f$ N@f$ particle response 
  // to the data in the range 
  // @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
  // 
  // Parameters:
  //    dist        Histogram to fit 
  //    lowCut      Lower cut @f$ E_{min}@f$ on signal 
  //    nParticles  Max number @f$ N@f$ of convolved landaus to fit
  //    minusBins   Number of bins @f$ \Delta b@f$ from peak to 
  //                    subtract to get the fit range 
  //    relErrorCut Cut applied to relative error of parameter. 
  //                    Note, for multi-particle weights, the cut 
  //                    is loosend by a factor of 2 
  //    chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
  //                    the reduced @f$\chi^2@f$ 
  // 
  // Return:
  //    The best fit function 
  //
  DGUARD(fDebug, 2, "Fit histogram in AliFMDEnergyFitter::RingHistos: %s",
	 dist->GetName());
  Double_t maxRange = 10;


  if (dist->GetEntries() <= 0) { 
    status = 1; // `empty'
    return 0;
  }
  Scale(dist);
  
  // Narrow search window for the peak 
  Int_t    cutBin  = TMath::Max(dist->GetXaxis()->FindBin(lowCut),3);
  Int_t    maxBin  = TMath::Min(dist->GetXaxis()->FindBin(10),
				dist->GetNbinsX());
  dist->GetXaxis()->SetRange(cutBin, maxBin);
    
  // Get the bin with maximum 
  Int_t    peakBin = dist->GetMaximumBin();
    
  // Normalize peak to 1 
  // Double_t max = dist->GetMaximum(); 
  Double_t max = dist->GetBinContent(peakBin); // Maximum(); 
  if (max <= 0) {
    status = 1; // `empty'
    return 0;
  }
  if (scaleToPeak) dist->Scale(1/max);
  DMSG(fDebug,5,"max(%s) -> %f", dist->GetName(), max);

  // Check that we have enough entries 
  Double_t nEntries = dist->GetEntries();
  if (nEntries <= minEntries) { 
    AliWarning(Form("Histogram at %s has too few entries (%f <= %d)",
		    dist->GetName(), nEntries, minEntries));
    status = 2;
    return 0;
  }

  // Create a fitter object 
  AliLandauGausFitter f(lowCut, maxRange, minusBins); 
  f.Clear();
  f.SetDebug(fDebug > 2); 

  // regularization cut - should be a parameter of the class 
  if (dist->GetEntries() > regCut) { 
    // We should rescale the errors 
    Double_t s = TMath::Sqrt(dist->GetEntries() / regCut);
    if (fDebug > 0) printf("Error scale: %f ", s);
    for (Int_t i = 1; i <= dist->GetNbinsX(); i++) {
      Double_t e = dist->GetBinError(i);
      dist->SetBinError(i, e * s);
    }
  }
  // If we are only asked to fit a single particle, return this fit, 
  // no matter what. 
  if (nParticles == 1) {
    TF1* r = f.Fit1Particle(dist, 0);
    if (!r) {
      status = 3; // No-fit
      return 0;
    }
    TF1* ff = new TF1(*r);
    dist->GetListOfFunctions()->Add(ff);
    ELossFit_t* ret = new ELossFit_t(0, *ff);
    ret->CalculateQuality(chi2nuCut, relErrorCut, minWeight);
    status = 0; // OK
    return ret;
  }

  // Fit from 2 upto n particles  
  for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);

  // Now, we need to select the best fit 
  Int_t nFits = f.GetFitResults().GetEntriesFast();
  for (Int_t i = nFits-1; i >= 0; i--) { 
    TF1* ff = static_cast<TF1*>(f.GetFunctions().At(i));
    // ff->SetRange(0, maxRange);
    dist->GetListOfFunctions()->Add(new TF1(*ff));
  }
  status = 0; // OK

  // Here, we use the real quality assesor instead of the old
  // `CheckResult' to ensure consitency in all output.
  ELossFit_t* ret = FindBestFit(dist, relErrorCut, chi2nuCut, minWeight);
  if (!ret) status = 3;
  return ret;
}

//__________________________________________________________________
AliFMDEnergyFitter::RingHistos::ELossFit_t* 
AliFMDEnergyFitter::RingHistos::FindBestFit(const TH1* dist,
					    Double_t   relErrorCut, 
					    Double_t   chi2nuCut,
					    Double_t   minWeightCut)  const
{
  // 
  // Find the best fit 
  // 
  // Parameters:
  //    dist           Histogram 
  //    relErrorCut    Cut applied to relative error of parameter. 
  //                       Note, for multi-particle weights, the cut 
  //                       is loosend by a factor of 2 
  //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
  //                       the reduced @f$\chi^2@f$ 
  //    minWeightCut   Least valid @f$ a_i@f$ 
  // 
  // Return:
  //    Best fit or null
  //
  TList* funcs = dist->GetListOfFunctions();
  TF1*   func  = 0;
  Int_t  i     = 0;
  TIter  next(funcs);
  fFits.Clear(); // This is only ever used here

  if (fDebug) printf("Find best fit for %s ... ", dist->GetName());
  if (fDebug > 2) printf("\n");

  // Loop over all functions stored in distribution, 
  // and calculate the quality 
  while ((func = static_cast<TF1*>(next()))) { 
    ELossFit_t* fit = new(fFits[i++]) ELossFit_t(0,*func);
    fit->fDet  = fDet;
    fit->fRing = fRing;
    // fit->fBin  = b;
    fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
    if (fDebug > 2) 
      Printf("%10s: %3d (chi^2/nu: %6.3f)", 
	     func->GetName(), fit->fQuality, 
	     (fit->fNu > 0 ? fit->fChi2 / fit->fNu : 999));
  }

  // Sort all the found fit objects in increasing quality 
  fFits.Sort();
  if (fDebug > 1) fFits.Print("s");

  // Get the top-most fit
  ELossFit_t* ret = static_cast<ELossFit_t*>(fFits.At(i-1));
  if (!ret) {
    AliWarningF("No fit found for %s", GetName());
    return 0;
  }
  if (ret && fDebug > 0) {
    if (fDebug > 1) printf(" %d: ", i-1);
    ret->Print("s");
  }
  // We have to make a copy here, because other wise the clones array
  // will overwrite the address
  return new ELossFit_t(*ret);
}

//____________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::CalculateResiduals(EResidualMethod mode, 
						   Double_t        lowCut,
						   TH1*            dist, 
						   ELossFit_t*     fit, 
						   TCollection*    out) const
{

  // Clone the input, and reset
  TH1*    resi = static_cast<TH1*>(dist->Clone());
  TString tit(resi->GetTitle());
  tit.ReplaceAll("#DeltaE/#DeltaE_{mip}", "Residuals");
  resi->SetTitle(tit);
  resi->SetDirectory(0);

  // Set title on Y axis
  switch (mode) { 
  case kResidualDifference:       
    resi->SetYTitle("h_{i}-f(#Delta_{i}) #pm #delta_{i}");
    break;
  case kResidualScaledDifference:  
    resi->SetYTitle("[h_{i}-f(#Delta_{i})]/#delta_{i}"); break;
  case kResidualSquareDifference:  
    resi->SetYTitle("#chi_{i}^{2}=[h_{i}-f(#Delta_{i})]^{2}/#delta^{2}_{i}");
    break;
  default: 
    resi->SetYTitle("Unknown");
    break;
  }
  out->Add(resi);

  // Try to find the function 
  Double_t highCut = dist->GetXaxis()->GetXmax();
  TString funcName("landau1");
  if (fit->GetN() > 1) 
    funcName = Form("nlandau%d", fit->GetN());
  TF1* func = dist->GetFunction(funcName);
  if (func) func->GetRange(lowCut, highCut);
  resi->Reset("ICES");
  resi->GetListOfFunctions()->Clear();
  resi->SetUniqueID(mode);

    // Reset histogram
  Int_t nX = resi->GetNbinsX();
  for (Int_t i  = 1; i <= nX; i++) { 
    Double_t x  = dist->GetBinCenter(i);
    if (x < lowCut)  continue;
    if (x > highCut) break;

    Double_t h  = dist->GetBinContent(i);
    Double_t e  = dist->GetBinError(i);
    Double_t r  = 0;
    Double_t er = 0;
    if (h > 0 && e > 0) { 
      Double_t f = fit->GetC() * fit->Evaluate(x);
      if (f > 0) { 
	r  = h-f;
	switch (mode) { 
	case kResidualDifference: er = e; break;
	case kResidualScaledDifference:  r /= e; break;
	case kResidualSquareDifference:  r *= r; r /= (e*e); break;
	default: r = 0; break;
	}
      }
    }
    resi->SetBinContent(i, r);
    resi->SetBinError(i, er);
  }  
}

//__________________________________________________________________
void
AliFMDEnergyFitter::RingHistos::FindBestFits(const TList*        d, 
					     AliFMDCorrELossFit& obj,
					     const TAxis&        eta)
{
  // 
  // Find the best fits 
  // 
  // Parameters:
  //    d              Parent list
  //    obj            Object to add fits to
  //    eta            Eta axis 
  //    relErrorCut    Cut applied to relative error of parameter. 
  //                       Note, for multi-particle weights, the cut 
  //                       is loosend by a factor of 2 
  //    chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
  //                       the reduced @f$\chi^2@f$ 
  //    minWeightCut   Least valid @f$ a_i@f$ 
  //

  // Get our ring list from the output container 
  TList* l = GetOutputList(d);
  if (!l) return; 

  // Get the energy distributions from the output container 
  TList* dists = static_cast<TList*>(l->FindObject("elossDists"));
  if (!dists) { 
    AliWarningF("Didn't find elossDists in %s", l->GetName());
    l->ls();
    return;
  }
  Int_t nBin = eta.GetNbins();
  if (fBest.GetEntriesFast() <= 0) { 
    AliWarningF("No fits found for %s", GetName());
    return;
  }
  
  for (Int_t b = 1; b <= nBin; b++) { 
    ELossFit_t* best = static_cast<ELossFit_t*>(fBest.At(b));
    if (!best) { 
      // AliErrorF("No best fit found @ %d for %s", b, GetName());
      continue;
    }
    // FindBestFit(b, dist, relErrorCut, chi2nuCut, minWeightCut);
    best->fDet  = fDet; 
    best->fRing = fRing;
    best->fBin  = b; // 
    if (fDebug > 0) {
      printf("Bin # %3d: ", b); 
      best->Print("s");
    }
    // Double_t eta = fAxis->GetBinCenter(b);
    obj.SetFit(fDet, fRing, b, best); 
    // new ELossFit_t(*best));
  }
}



//____________________________________________________________________
//
// EOF
//