More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.cxx
index 34464d80aecc8863a4675862edfb2e05909be8fb..76575473ca8fc07970b510add1284273e584c505 100644 (file)
 // Class to do the energy correction of FMD ESD data
 //
 #include "AliFMDEnergyFitter.h"
+#include "AliForwardCorrectionManager.h"
 #include <AliESDFMD.h>
 #include <TAxis.h>
 #include <TList.h>
 #include <TH1.h>
 #include <TF1.h>
 #include <TMath.h>
-#include "AliFMDAnaParameters.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 
-
-#define N_A(N)  (4+N-2)
-#define N2_A(N) (4+(N-2)*3)
-#define N2_D(N) (4+(N-2)*3+1)
-#define N2_X(N) (4+(N-2)*3+2)
-
-//____________________________________________________________________
 namespace {
-  Double_t 
-  NLandau(Double_t* xp, Double_t* pp) 
-  {
-    Double_t  e        = xp[0];
-    Double_t  constant = pp[0];
-    Double_t  mpv      = pp[1];
-    Double_t  fwhm     = pp[2];
-    Int_t     n        = Int_t(pp[3]);
-    Double_t  result   = 0;
-    for (Int_t i = 1; i <= n; i++) {
-      Double_t mpvi  =  i*(mpv+fwhm*TMath::Log(i));
-      Double_t fwhmi =  i*fwhm;
-      Double_t ai    =  (i == 1 ? 1 : pp[N_A(i)]);
-      result         += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE);
-    }
-    result *= constant;
-    return result;
-  }
-
-  Double_t 
-  NLandau2(Double_t* xp, Double_t* pp) 
-  {
-    Double_t  e        = xp[0];
-    Double_t  constant = pp[0];
-    Double_t  mpv      = pp[1];
-    Double_t  fwhm     = pp[2];
-    Int_t     n        = Int_t(pp[3]);
-    Double_t  result   = TMath::Landau(e,mpv,fwhm,kTRUE);
-    for (Int_t i = 2; i <= n; i++) {
-      Double_t ai    =  pp[N2_A(i)];
-      Double_t mpvi  =  pp[N2_D(i)];
-      Double_t fwhmi =  pp[N2_X(i)];
-      result         += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE);
-    }
-    result *= constant;
-    return result;
-  }
-
-  Double_t 
-  TripleLandau(Double_t *x, Double_t *par) 
-  {
-    Double_t energy   = x[0];
-    Double_t constant = par[0];
-    Double_t mpv      = par[1];
-    Double_t fwhm     = par[2];
-    Double_t alpha    = par[3];
-    Double_t beta     = par[4];
-    Double_t mpv2     = 2*mpv+2*fwhm*TMath::Log(2);
-    Double_t mpv3     = 3*mpv+3*fwhm*TMath::Log(3);
-
-    Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+
-                            alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE)+
-                            beta  * TMath::Landau(energy,mpv3,3*fwhm,kTRUE));
-  
-    return f;
-  }
-  Double_t 
-  DoubleLandau(Double_t *x, Double_t *par) 
-  {
-    Double_t energy   = x[0];
-    Double_t constant = par[0];
-    Double_t mpv      = par[1];
-    Double_t fwhm     = par[2];
-    Double_t alpha    = par[3];
-    Double_t mpv2     = 2*mpv+2*fwhm*TMath::Log(2);
-    
-    Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+
-                            alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE));
-  
-    return f;
-  }
-  Double_t 
-  SingleLandau(Double_t *x, Double_t *par) 
-  {
-    Double_t energy   = x[0];
-    Double_t constant = par[0];
-    Double_t mpv      = par[1];
-    Double_t fwhm     = par[2];
-    
-    Double_t f = constant * TMath::Landau(energy,mpv,fwhm,kTRUE);
-  
-    return f;
-  }
+  const char* fgkEDistFormat = "%s_etabin%03d";
 }
 
+
 //____________________________________________________________________
 AliFMDEnergyFitter::AliFMDEnergyFitter()
   : TNamed(), 
     fRingHistos(),
     fLowCut(0.3),
-    fNLandau(3),
+    fNParticles(3),
     fMinEntries(100),
-    fBinsToSubtract(4),
+    fFitRangeBinWidth(4),
     fDoFits(false),
+    fDoMakeObject(false),
     fEtaAxis(),
+    fMaxE(10),
+    fNEbins(300), 
+    fUseIncreasingBins(true),
+    fMaxRelParError(.25),
+    fMaxChi2PerNDF(20), 
+    fMinWeight(1e-7),
     fDebug(0)
 {}
 
@@ -129,11 +51,18 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
   : TNamed("fmdEnergyFitter", title), 
     fRingHistos(), 
     fLowCut(0.3),
-    fNLandau(3),
+    fNParticles(3),
     fMinEntries(100),
-    fBinsToSubtract(4),
+    fFitRangeBinWidth(4),
     fDoFits(false),
-    fEtaAxis(100,-4,6),
+    fDoMakeObject(false),
+    fEtaAxis(0,0,0),
+    fMaxE(10),
+    fNEbins(300), 
+    fUseIncreasingBins(true),
+    fMaxRelParError(.25),
+    fMaxChi2PerNDF(20),  
+    fMinWeight(1e-7),
     fDebug(3)
 {
   fEtaAxis.SetName("etaAxis");
@@ -150,11 +79,18 @@ AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
   : TNamed(o), 
     fRingHistos(), 
     fLowCut(o.fLowCut),
-    fNLandau(o.fNLandau),
+    fNParticles(o.fNParticles),
     fMinEntries(o.fMinEntries),
-    fBinsToSubtract(o.fBinsToSubtract),
+    fFitRangeBinWidth(o.fFitRangeBinWidth),
     fDoFits(o.fDoFits),
+    fDoMakeObject(o.fDoMakeObject),
     fEtaAxis(o.fEtaAxis),
+    fMaxE(o.fMaxE),
+    fNEbins(o.fNEbins), 
+    fUseIncreasingBins(o.fUseIncreasingBins),
+    fMaxRelParError(o.fMaxRelParError),
+    fMaxChi2PerNDF(o.fMaxChi2PerNDF), 
+    fMinWeight(o.fMinWeight),
     fDebug(o.fDebug)
 {
   TIter    next(&o.fRingHistos);
@@ -175,12 +111,16 @@ AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
   TNamed::operator=(o);
 
   fLowCut        = o.fLowCut;
-  fNLandau       = o.fNLandau;
+  fNParticles       = o.fNParticles;
   fMinEntries    = o.fMinEntries;
-  fBinsToSubtract= o.fBinsToSubtract;
+  fFitRangeBinWidth= o.fFitRangeBinWidth;
   fDoFits        = o.fDoFits;
+  fDoMakeObject  = o.fDoMakeObject;
   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
   fDebug         = o.fDebug;
+  fMaxRelParError= o.fMaxRelParError;
+  fMaxChi2PerNDF = o.fMaxChi2PerNDF;
+  fMinWeight     = o.fMinWeight;
 
   fRingHistos.Delete();
   TIter    next(&o.fRingHistos);
@@ -209,12 +149,27 @@ AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
 void
 AliFMDEnergyFitter::Init(const TAxis& eAxis)
 {
-  fEtaAxis.Set(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
+  if (fEtaAxis.GetNbins() == 0 || 
+      fEtaAxis.GetXmin() == fEtaAxis.GetXmax()) 
+    SetEtaAxis(eAxis);
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   while ((o = static_cast<RingHistos*>(next())))
-    o->Init(eAxis);
+    o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
 }  
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
+{
+  SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
+}
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
+{
+  fEtaAxis.Set(nBins,etaMin,etaMax);
+}
+
 //____________________________________________________________________
 Bool_t
 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, 
@@ -264,30 +219,60 @@ AliFMDEnergyFitter::Fit(TList* dir)
   // +1          for chi^2
   // +3          for 1 landau 
   // +1          for N 
-  // +fNLandau-1 for weights 
-  Int_t nStack = 1+3+1+fNLandau-1;
+  // +fNParticles-1 for weights 
+  Int_t nStack = kN+fNParticles;
   THStack* stack[nStack]; 
-  stack[0] = new THStack("chi2", "#chi^{2}/#nu");
-  stack[1] = new THStack("c",    "constant");
-  stack[2] = new THStack("mpv",  "#Delta_{p}");
-  stack[3] = new THStack("w",    "w");
-  stack[4] = new THStack("n",    "# of Landaus");
-  for (Int_t i = 2; i <= fNLandau; i++) 
-    stack[i-1+4] = new THStack(Form("a%d", i), 
-                              Form("Weight of %d signal", i));
+  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]);
 
   TIter    next(&fRingHistos);
   RingHistos* o = 0;
   while ((o = static_cast<RingHistos*>(next()))) {
-    TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNLandau,
-                         fMinEntries, fBinsToSubtract);
+    TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
+                         fMinEntries, fFitRangeBinWidth,
+                         fMaxRelParError, fMaxChi2PerNDF);
     if (!l) continue;
     for (Int_t i = 0; i < l->GetEntriesFast(); i++) { 
       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i))); 
     }
   }
+
+  if (!fDoMakeObject) return;
+
+  MakeCorrectionsObject(d);
+}
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
+{
+  AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
+    
+  AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
+  obj->SetEtaAxis(fEtaAxis);
+  obj->SetLowCut(fLowCut);
+
+  TIter    next(&fRingHistos);
+  RingHistos* o = 0;
+  while ((o = static_cast<RingHistos*>(next()))) {
+    o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError, 
+                   fMaxChi2PerNDF, fMinWeight);
+  }
+  
+  TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
+  TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
+  AliInfo(Form("Object %s created in output - should be extracted and copied "
+              "to %s", oName.Data(), fileName.Data()));
+  d->Add(obj, oName.Data());
 }
 
 //____________________________________________________________________
@@ -314,6 +299,32 @@ AliFMDEnergyFitter::SetDebug(Int_t dbg)
   while ((o = static_cast<RingHistos*>(next())))
     o->fDebug = dbg;
 }  
+
+//____________________________________________________________________
+void
+AliFMDEnergyFitter::Print(Option_t*) const
+{
+  char ind[gROOT->GetDirLevel()+1];
+  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
+  ind[gROOT->GetDirLevel()] = '\0';
+  std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n'
+           << ind << " Low cut:                " << fLowCut << " E/E_mip\n"
+           << ind << " Max(particles):         " << fNParticles << '\n'
+           << ind << " Min(entries):           " << fMinEntries << '\n'
+           << ind << " Fit range:              " 
+           << fFitRangeBinWidth << " bins\n"
+           << ind << " Make fits:              " 
+           << (fDoFits ? "yes\n" : "no\n")
+           << ind << " Make object:            " 
+           << (fDoMakeObject ? "yes\n" : "no\n")
+           << ind << " Max E/E_mip:            " << fMaxE << '\n'
+           << ind << " N bins:                 " << fNEbins << '\n'
+           << ind << " Increasing bins:        " 
+           << (fUseIncreasingBins ?"yes\n" : "no\n")
+           << ind << " max(delta p/p):         " << fMaxRelParError << '\n'
+           << ind << " max(chi^2/nu):          " << fMaxChi2PerNDF << '\n'
+           << ind << " min(a_i):               " << fMinWeight << std::endl;
+}
   
 //====================================================================
 AliFMDEnergyFitter::RingHistos::RingHistos()
@@ -322,6 +333,7 @@ AliFMDEnergyFitter::RingHistos::RingHistos()
     fEmpty(0),
     fEtaEDists(), 
     fList(0),
+    fFits("AliFMDCorrELossFit::ELossFit"),
     fDebug(0)
 {}
 
@@ -332,6 +344,7 @@ AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
     fEmpty(0),
     fEtaEDists(), 
     fList(0),
+    fFits("AliFMDCorrELossFit::ELossFit"),
     fDebug(0)
 {
   fEtaEDists.SetName("EDists");
@@ -343,8 +356,10 @@ AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
     fEmpty(o.fEmpty),
     fEtaEDists(), 
     fList(0),
+    fFits("AliFMDCorrELossFit::ELossFit"),
     fDebug(0)
 {
+  fFits.Clear();
   TIter next(&o.fEtaEDists);
   TObject* obj = 0;
   while ((obj = next())) fEtaEDists.Add(obj->Clone());
@@ -366,6 +381,7 @@ AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
   if (fEmpty) delete fEmpty;
   fEtaEDists.Delete();
   if (fList) fList->Delete();
+  fFits.Clear();
 
   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
@@ -434,12 +450,15 @@ AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
 
 //____________________________________________________________________
 void
-AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax,
-                                    Double_t deMax, Int_t nDeBins, 
-                                    Bool_t incr)
+AliFMDEnergyFitter::RingHistos::Make(Int_t    ieta, 
+                                    Double_t emin, 
+                                    Double_t emax,
+                                    Double_t deMax, 
+                                    Int_t    nDeBins, 
+                                    Bool_t   incr)
 {
   TH1D* h = 0;
-  TString name  = Form("%s_etabin%03d", fName.Data(), ieta);
+  TString name  = Form(fgkEDistFormat, fName.Data(), ieta);
   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
                       "(#eta bin %d)", fName.Data(), emin, emax, ieta);
   if (!incr) 
@@ -460,10 +479,9 @@ AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax,
   fEtaEDists.AddAt(h, ieta-1);
 }
 //____________________________________________________________________
-TH1D*
-AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
-                                       const char* title, 
-                                       const TAxis& eta) const
+TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
+                                             const char* title, 
+                                             const TAxis& eta) const
 {
   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
                     Form("%s for %s", title, fName.Data()), 
@@ -533,14 +551,20 @@ AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis,
 }
 //____________________________________________________________________
 TObjArray*
-AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
-                                   Double_t lowCut, UShort_t nLandau,
-                                   UShort_t minEntries,
-                                   UShort_t minusBins) const
+AliFMDEnergyFitter::RingHistos::Fit(TList*           dir, 
+                                   const TAxis&     eta,
+                                   Double_t         lowCut, 
+                                   UShort_t         nParticles,
+                                   UShort_t         minEntries,
+                                   UShort_t         minusBins, 
+                                   Double_t         relErrorCut, 
+                                   Double_t         chi2nuCut) const
 {
+  // Get our ring list from the output container 
   TList* l = GetOutputList(dir);
   if (!l) return 0; 
 
+  // Get the energy distributions from the output container 
   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
   if (!dists) { 
     AliWarning(Form("Didn't find %s_EtaEDists in %s", 
@@ -549,24 +573,28 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
     return 0;
   }
 
-  TObjArray* pars  = new TObjArray(3+nLandau+1);
+  // Container of the fit results histograms 
+  TObjArray* pars  = new TObjArray(3+nParticles+1);
   pars->SetName("FitResults");
   l->Add(pars);
 
-  TH1* hChi2 = 0;
-  TH1* hC    = 0;
-  TH1* hMpv  = 0;
-  TH1* hW    = 0;
-  TH1* hS    = 0;
-  TH1* hN    = 0;
-  TH1* hA[nLandau-1];
-  pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta));
-  pars->Add(hC    = MakePar("c",    "Constant", eta));
-  pars->Add(hMpv  = MakePar("mpv",  "#Delta_{p}", eta));
-  pars->Add(hW    = MakePar("w",    "#xi", eta));
-  pars->Add(hS    = MakePar("s",    "#sigma", eta));
-  pars->Add(hN    = MakePar("n",    "N", eta));
-  for (UShort_t i = 1; i < nLandau; i++) 
+  // 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));
 
   
@@ -575,62 +603,125 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
   Int_t high   = 0;
   for (Int_t i = 0; i < nDists; i++) { 
     TH1D* dist = static_cast<TH1D*>(dists->At(i));
+    // Ignore empty histograms altoghether 
+    if (!dist || dist->GetEntries() <= 0) continue; 
 
-    if (!dist) continue;
+    // Scale to the bin-width
+    dist->Scale(1., "width");
+    
+    // Normalize peak to 1 
+    Double_t max = dist->GetMaximum(); 
+    if (max <= 0) continue;
+    dist->Scale(1/max);
+    
+    // Check that we have enough entries 
+    if (dist->GetEntries() <= minEntries) { 
+      AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
+                     i, dist->GetName(), Int_t(dist->GetEntries()), 
+                     minEntries));
+      continue;
+    }
 
-    TF1* res = FitHist(dist,lowCut,nLandau,minEntries,minusBins);
+    // Now fit 
+    TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
+                      relErrorCut,chi2nuCut);
     if (!res) continue;
-    
+    // dist->GetListOfFunctions()->Add(res);
+
+    // Store eta limits 
     low   = TMath::Min(low,i+1);
     high  = TMath::Max(high,i+1);
 
+    // Get the reduced chi square
     Double_t chi2 = res->GetChisquare();
     Int_t    ndf  = res->GetNDF();
-    hChi2->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
-    hC  ->SetBinContent(i+1, res->GetParameter(0));   
-    hMpv->SetBinContent(i+1, res->GetParameter(1)); 
-    hW  ->SetBinContent(i+1, res->GetParameter(2));   
-    hN  ->SetBinContent(i+1, res->GetParameter(3));   
-
-    hC  ->SetBinError(i+1, res->GetParError(1));
-    hMpv->SetBinError(i+1, res->GetParError(2));
-    hW  ->SetBinError(i+1, res->GetParError(2));
-    // hN  ->SetBinError(i, res->GetParError(3));
-
-    for (UShort_t j = 0; j < nLandau-1; j++) {
-      hA[j]->SetBinContent(i+1, res->GetParameter(4+j));
-      hA[j]->SetBinError(i+1, res->GetParError(4+j));
+    
+    // Store results of best fit in output histograms 
+    res->SetLineWidth(4);
+    hChi2   ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
+    hC      ->SetBinContent(i+1, res->GetParameter(kC));   
+    hDelta  ->SetBinContent(i+1, res->GetParameter(kDelta)); 
+    hXi     ->SetBinContent(i+1, res->GetParameter(kXi));   
+    hSigma  ->SetBinContent(i+1, res->GetParameter(kSigma));   
+    hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));   
+    hN      ->SetBinContent(i+1, res->GetParameter(kN));   
+
+    hC     ->SetBinError(i+1, res->GetParError(kC));
+    hDelta ->SetBinError(i+1, res->GetParError(kDelta));
+    hXi    ->SetBinError(i+1, res->GetParError(kXi));
+    hSigma ->SetBinError(i+1, res->GetParError(kSigma));
+    hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
+    hN     ->SetBinError(i+1, res->GetParError(kN));
+
+    for (UShort_t j = 0; j < nParticles-1; j++) {
+      hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
+      hA[j]->SetBinError(i+1, res->GetParError(kA+j));
     }
   }
 
+  // Fit the full-ring histogram 
   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
-  if (total) { 
-    TF1* res = FitHist(total,lowCut,nLandau,minEntries,minusBins);
+  if (total && total->GetEntries() >= minEntries) { 
+
+    // Scale to the bin-width
+    total->Scale(1., "width");
+    
+    // Normalize peak to 1 
+    Double_t max = total->GetMaximum(); 
+    if (max > 0) total->Scale(1/max);
+
+    TF1* res = FitHist(total,lowCut,nParticles,minusBins,
+                      relErrorCut,chi2nuCut);
     if (res) { 
+      // Make histograms for the result of this fit 
       Double_t chi2 = res->GetChisquare();
       Int_t    ndf  = res->GetNDF();
-      pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
+      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,
-                         res->GetParameter(0),res->GetParError(0)));
-      pars->Add(MakeTotal("t_mpv",  "#Delta_{p}",   eta, low, high,
-                         res->GetParameter(1),res->GetParError(1)));
-      pars->Add(MakeTotal("t_w",    "#xi",          eta, low, high,
-                         res->GetParameter(2),res->GetParError(2)));
-      pars->Add(MakeTotal("t_n",    "N",            eta, low, high,
-                         res->GetParameter(3),0));
-      for (UShort_t j = 1; j < nLandau; j++) 
-       pars->Add(MakeTotal(Form("a%d_t",j+1), 
-                           Form("a_{%d}",j+1), eta, low, high,
-                           res->GetParameter(3+j), res->GetParError(3+j)));
+      pars->Add(MakeTotal("t_c",      "Constant",     eta, low, high,
+                         res->GetParameter(kC),
+                         res->GetParError(kC)));
+      pars->Add(MakeTotal("t_delta",  "#Delta_{p}",   eta, low, high,
+                         res->GetParameter(kDelta),
+                         res->GetParError(kDelta)));
+      pars->Add(MakeTotal("t_xi",     "#xi",          eta, low, high,
+                         res->GetParameter(kXi),
+                         res->GetParError(kXi)));
+      pars->Add(MakeTotal("t_sigma",  "#sigma",       eta, low, high,
+                         res->GetParameter(kSigma),
+                         res->GetParError(kSigma)));
+      pars->Add(MakeTotal("t_sigman", "#sigma_{n}",   eta, low, high,
+                         res->GetParameter(kSigmaN),
+                         res->GetParError(kSigmaN)));
+      pars->Add(MakeTotal("t_n",    "N",              eta, low, high,
+                         res->GetParameter(kN),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,
+                           res->GetParameter(kA+j), 
+                           res->GetParError(kA+j)));
     }
   }
     
+  // Clean up list of histogram.  Histograms with no entries or 
+  // no functions are deleted.  We have to do this using the TObjLink 
+  // objects stored in the list since ROOT cannot guaranty the validity 
+  // of iterators when removing from a list - tsck.  Should just implement
+  // TIter::Remove(). 
   TObjLink* lnk = dists->FirstLink();
   while (lnk) {
     TH1* h = static_cast<TH1*>(lnk->GetObject());
-    if (h->GetEntries() <= 0 || 
-       h->GetListOfFunctions()->GetEntries() <= 0) {
+    bool remove = false;
+    if (h->GetEntries() <= 0) { 
+      // AliWarning(Form("No entries in %s - removing", h->GetName()));
+      remove = true;
+    }
+    else if (h->GetListOfFunctions()->GetEntries() <= 0) {
+      // AliWarning(Form("No fuctions associated with %s - removing",
+      //            h->GetName()));
+      remove = true;
+    }
+    if (remove) {
       TObjLink* keep = lnk->Next();
       dists->Remove(lnk);
       lnk = keep;
@@ -644,243 +735,195 @@ AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
 
 //____________________________________________________________________
 TF1*
-AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,Double_t lowCut, 
-                                       UShort_t nLandau, 
-                                       UShort_t minEntries,
-                                       UShort_t minusBins) const
+AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
+                                       Double_t lowCut, 
+                                       UShort_t nParticles, 
+                                       UShort_t minusBins, 
+                                       Double_t relErrorCut, 
+                                       Double_t chi2nuCut) const
 {
   Double_t maxRange = 10;
 
-  if (dist->GetEntries() <= minEntries) return 0;
-
-  // Find the fit range 
-  dist->GetXaxis()->SetRangeUser(lowCut, maxRange);
-  
-  // Normalize peak to 1 
-  Double_t max = dist->GetMaximum(); 
-  dist->Scale(1/max);
+  // Create a fitter object 
+  AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins); 
+  f.Clear();
   
-  // Get the bin with maximum 
-  Int_t    maxBin = dist->GetMaximumBin();
-  Double_t maxE   = dist->GetBinLowEdge(maxBin);
   
-  // Get the low edge 
-  dist->GetXaxis()->SetRangeUser(lowCut, maxE);
-  Int_t    minBin = maxBin - minusBins; // dist->GetMinimumBin();
-  Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),lowCut);
-  Double_t maxEE  = dist->GetBinCenter(maxBin+2*minusBins);
-
-  // Restore the range 
-  dist->GetXaxis()->SetRangeUser(0, maxRange);
-  
-  // First do a single landau fit 
-  TF1*          landau1 = new TF1("landau1", "landau", minE, maxEE);
-  // landau1->SetParameters(1,1,1,1);
-  landau1->SetParNames("C","#Delta_{p}","#xi");
-  landau1->SetParLimits(1,minE,maxRange);
-  landau1->SetParLimits(2,0,maxRange);
-  landau1->SetLineColor(Color());
-  // Tight fit around peak - make sure we get that right. 
-  TFitResultPtr r = dist->Fit(landau1, "RQNS", "", minE, maxEE);
-
-  return FitMore2(dist, nLandau, *r, landau1, minE, maxRange);
-}
-
-//____________________________________________________________________
-TF1*
-AliFMDEnergyFitter::RingHistos::FitMore(TH1*        dist,
-                                       UShort_t    nLandau,
-                                       TFitResult& r, 
-                                       TF1*        landau1,
-                                       Double_t    minE,
-                                       Double_t    maxRange) const
-{
-  static TClonesArray res("TFitResult");
-  static TObjArray funcs;
-  res.Clear();
-  funcs.SetOwner();
-  funcs.Clear();
-  Int_t nRes = 0;
-
-  //r->Print();
-  new(res[nRes++]) TFitResult(r);
-  funcs.AddAtAndExpand(landau1, 0);
-
-  // Now try to fit 
-  for (UShort_t n = 2; n <= nLandau; n++) { 
-    TFitResult* rr = static_cast<TFitResult*>(res.At(nRes-1));
-    if (!rr) { 
-      AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n));
-      return 0;
-    }
-    // Create the function object 
-    Double_t mpvi  = rr->Parameter(1);
-    Double_t wi    = rr->Parameter(2);
-    Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi;
-    TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,3+n);
-    landaui->SetLineStyle((n % 10)+1);
-    landaui->SetLineWidth(2);
-    landaui->SetLineColor(Color());
-    landaui->SetParameter(0, rr->Parameter(0));
-    landaui->SetParameter(1, rr->Parameter(1));
-    landaui->SetParameter(2, rr->Parameter(2));
-    landaui->SetParLimits(1, minE, maxRange);
-    landaui->SetParLimits(2,0,maxRange);
-    landaui->FixParameter(3, n);
-    landaui->SetParNames("C","#Delta_{p}","#xi", "N");
-    for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit 
-      landaui->SetParameter(N_A(i), rr->Parameter(N_A(i)));
-      landaui->SetParLimits(N_A(i), 0,1);
-      landaui->SetParName(i, Form("a_{%d}", i));
-    }
-    landaui->SetParameter(N_A(n), (n == 2 ? 0.05 : rr->Parameter(N_A(n-1))/5));
-    landaui->SetParLimits(N_A(n), 0, 1);
-    landaui->SetParName(N_A(n), Form("a_{%d}", n));
-    // landaui->Print();
-
-    TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi);
-    // tr->Print();
-    if (CheckResult(*tr)) { 
-      new(res[nRes++]) TFitResult(*tr);
-      funcs.AddAtAndExpand(landaui,n-1);
-      continue;
-    }
-    // Stop on bad fit 
-    break;
+  // 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) return 0;
+    return new TF1(*r);
   }
-  TF1* ret = 0;
-  if (funcs.At(nRes-1)) ret = static_cast<TF1*>(funcs.At(nRes-1)->Clone());
-  dist->GetListOfFunctions()->Add(ret);
 
-  res.Clear();
-  funcs.Clear();
-
-  return ret;
-}
-
-//____________________________________________________________________
-TF1*
-AliFMDEnergyFitter::RingHistos::FitMore2(TH1*        dist,
-                                        UShort_t    nLandau,
-                                        TFitResult& r, 
-                                        TF1*        landau1,
-                                        Double_t    minE,
-                                        Double_t    maxRange) const
-{
-  static TClonesArray res("TFitResult");
-  static TObjArray funcs;
-  res.Clear();
-  funcs.SetOwner();
-  funcs.Clear();
-  Int_t nRes = 0;
-
-  //r->Print();
-  new(res[nRes++]) TFitResult(r);
-  funcs.AddAtAndExpand(landau1, 0);
-
-  // Now try to fit 
-  for (UShort_t n = 2; n <= nLandau; n++) { 
-    TFitResult* rr = static_cast<TFitResult*>(res.At(nRes-1));
-    if (!rr) { 
-      AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n));
-      return 0;
+  // 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.fFitResults.GetEntriesFast();
+  TF1*  good[nFits];
+  for (Int_t i = nFits-1; i >= 0; i--) { 
+    good[i] = 0;
+    TF1* ff = static_cast<TF1*>(f.fFunctions.At(i));
+    // ff->SetLineColor(Color());
+    ff->SetRange(0, maxRange);
+    dist->GetListOfFunctions()->Add(new TF1(*ff));
+    if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i)),
+                   relErrorCut, chi2nuCut)) {
+      good[i] = ff;
+      ff->SetLineWidth(2);
+      // f.fFitResults.At(i)->Print("V");
     }
-    // Create the function object 
-    Double_t mpvi  = rr->Parameter(1);
-    Double_t wi    = rr->Parameter(2);
-    Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi;
-    TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,
-                          n*3+1);
-    landaui->SetLineStyle((n % 10)+1);
-    landaui->SetLineWidth(2);
-    landaui->SetLineColor(Color());
-    landaui->SetParameter(0, rr->Parameter(0));
-    landaui->SetParameter(1, rr->Parameter(1));
-    landaui->SetParameter(2, rr->Parameter(2));
-    landaui->SetParLimits(1, minE, maxRange);
-    landaui->SetParLimits(2,0,maxRange);
-    landaui->FixParameter(3, n);
-    landaui->SetParNames("C","#Delta_{p}","#xi", "N");
-    for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit 
-      landaui->SetParameter(N2_A(i), rr->Parameter(N2_A(i)));
-      landaui->SetParameter(N2_D(i), rr->Parameter(N2_D(i)));
-      landaui->SetParameter(N2_X(i), rr->Parameter(N2_X(i)));
-      landaui->SetParLimits(N2_A(i), 0,1);
-      landaui->SetParLimits(N2_D(i), minE,maxEi);
-      landaui->SetParLimits(N2_X(i), 0,maxRange);
-      landaui->SetParName(N2_A(i), Form("a_{%d}", i));
-      landaui->SetParName(N2_D(i), Form("#Delta_{p,%d}", i));
-      landaui->SetParName(N2_X(i), Form("#xi_{%d}", i));
-    }
-    landaui->SetParameter(N2_A(n), n == 2 ? 0.05 : rr->Parameter(N2_A(n-1))/5);
-    landaui->SetParameter(N2_D(n), mpvi);
-    landaui->SetParameter(N2_X(n), wi);
-    landaui->SetParLimits(N2_A(n), 0,1);
-    landaui->SetParLimits(N2_D(n), minE,maxEi);
-    landaui->SetParLimits(N2_X(n), 0,maxRange);
-    landaui->SetParName(N2_A(n), Form("a_{%d}", n));
-    landaui->SetParName(N2_D(n), Form("#Delta_{p,%d}", n));
-    landaui->SetParName(N2_X(n), Form("#xi_{%d}", n));
-    if (fDebug > 2) landaui->Print();
-
-    TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi);
-    if (fDebug > 2)  tr->Print();
-    if (CheckResult(*tr)) { 
-      new(res[nRes++]) TFitResult(*tr);
-      funcs.AddAtAndExpand(landaui,n-1);
-      continue;
+  }
+  // If all else fails, use the 1 particle fit 
+  TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
+
+  // Find the fit with the most valid particles 
+  for (Int_t i = nFits-1; i > 0; i--) {
+    if (!good[i]) continue;
+    if (fDebug > 1) {
+      AliInfo(Form("Choosing fit with n=%d", i+1));
+      f.fFitResults.At(i)->Print();
     }
-    // Stop on bad fit 
+    ret = good[i];
     break;
   }
-  TF1* ret = 0;
-  if (funcs.At(nRes-1)) ret = static_cast<TF1*>(funcs.At(nRes-1)->Clone());
-  dist->GetListOfFunctions()->Add(ret);
-
-  res.Clear();
-  funcs.Clear();
-
-  return ret;
+  // Give a warning if we're using fall-back 
+  if (ret == f.fFunctions.At(0)) {
+    AliWarning("Choosing fall-back 1 particle fit");
+  }
+  // Copy our result and return (the functions are owned by the fitter)
+  TF1* fret = new TF1(*ret);
+  return fret;
 }
 
 //____________________________________________________________________
 Bool_t
-AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult& r) const
+AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
+                                           Double_t relErrorCut, 
+                                           Double_t chi2nuCut) const
 {
-  Double_t chi2 = r.Chi2();
-  Int_t    ndf  = r.Ndf();
+  if (fDebug > 10) r->Print();
+  TString  n    = r->GetName();
+  n.ReplaceAll("TFitResult-", "");
+  Double_t chi2 = r->Chi2();
+  Int_t    ndf  = r->Ndf();
   // Double_t prob = r.Prob();
-  if (ndf <= 0 || chi2 / ndf > 5) { 
-    if (fDebug > 2)
-      AliWarning(Form("Fit %s has chi^2/ndf=%f/%d=%f", 
-                     r.GetName(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf)));
-    return kFALSE;
+  Bool_t ret = kTRUE;
+  
+  // Check that the reduced chi square isn't larger than 5
+  if (ndf <= 0 || chi2 / ndf > chi2nuCut) { 
+    if (fDebug > 2) {
+      AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f", 
+                     n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
+                     chi2nuCut)); }
+    ret = kFALSE;
   }
     
-  UShort_t nPar = r.NPar();
+  // Check each parameter 
+  UShort_t nPar = r->NPar();
   for (UShort_t i = 0; i < nPar; i++) { 
-    if (i == 3) continue; 
+    if (i == kN) continue;  // Skip the number parameter 
     
-    Double_t v = r.Parameter(i);
-    Double_t e = r.ParError(i);
+    // Get value and error and check value 
+    Double_t v  = r->Parameter(i);
+    Double_t e  = r->ParError(i);
     if (v == 0) continue;
-    if (v == 0 || e / v > 0.2) { 
-      if (fDebug > 2)
-       AliWarning(Form("Fit %s has Delta %s/%s=%f/%f=%f%%",
-                       r.GetName(), r.ParName(i).c_str(), 
-                       r.ParName(i).c_str(), e, v, 100*(v == 0 ? 0 : e/v)));
-      return kFALSE;
+
+    // Calculate the relative error and check it 
+    Double_t re  = e / v;
+    Double_t cut = relErrorCut * (i < kN ? 1 : 2);
+    if (re > cut) { 
+      if (fDebug > 2) {
+       AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
+                       n.Data(), r->ParName(i).c_str(), 
+                       r->ParName(i).c_str(), e, v, 
+                       100*(v == 0 ? 0 : e/v),
+                       100*(cut))); }
+      ret = kFALSE;
     }
   }
-  if (nPar > 4) { 
-    if (r.Parameter(nPar-1) <= 1e-10) { 
-      if (fDebug)
-       AliWarning(Form("Last scale %s is too small %f<1e-10", 
-                       r.ParName(nPar-1).c_str(), r.Parameter(nPar-1)));
-      return kFALSE;
+
+  // Check if we have scale parameters 
+  if (nPar > kN) { 
+    
+    // Check that the last particle has a significant contribution 
+    Double_t lastScale = r->Parameter(nPar-1);
+    if (lastScale <= 1e-7) { 
+      if (fDebug) {
+       AliWarning(Form("%s: %s=%9.6f<1e-7", 
+                       n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
+      ret = kFALSE;
     }
   }
-  return kTRUE;
+  return ret;
+}
+
+
+//__________________________________________________________________
+void
+AliFMDEnergyFitter::RingHistos::FindBestFits(TList*              d, 
+                                            AliFMDCorrELossFit& obj,
+                                            const TAxis&        eta,     
+                                            Double_t            relErrorCut, 
+                                            Double_t            chi2nuCut,
+                                            Double_t            minWeightCut)
+{
+  // 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("EDists"));
+  if (!dists) { 
+    AliWarning(Form("Didn't find %s_EtaEDists in %s", 
+                   fName.Data(), l->GetName()));
+    l->ls();
+    return;
+  }
+  Int_t nBin = eta.GetNbins();
+
+  for (Int_t b = 1; b <= nBin; b++) { 
+    TString n(Form(fgkEDistFormat, fName.Data(), b));
+    TH1D*   dist = static_cast<TH1D*>(dists->FindObject(n));
+    // Ignore empty histograms altoghether 
+    if (!dist || dist->GetEntries() <= 0) continue; 
+    
+    AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
+                                                    relErrorCut,
+                                                    chi2nuCut,  
+                                                    minWeightCut);
+    best->fDet  = fDet; 
+    best->fRing = fRing;
+    best->fBin  = b; // 
+    // Double_t eta = fAxis->GetBinCenter(b);
+    obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
+  }
+}
+
+//__________________________________________________________________
+AliFMDCorrELossFit::ELossFit* 
+AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist,
+                                           Double_t relErrorCut, 
+                                           Double_t chi2nuCut,
+                                           Double_t minWeightCut) 
+{
+  TList* funcs = dist->GetListOfFunctions();
+  TIter  next(funcs);
+  TF1*   func = 0;
+  fFits.Clear();
+  Int_t  i = 0;
+  // Info("FindBestFit", "%s", dist->GetName());
+  while ((func = static_cast<TF1*>(next()))) { 
+    AliFMDCorrELossFit::ELossFit* fit = 
+      new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
+    fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
+  }
+  fFits.Sort(false);
+  // fFits.Print();
+  return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
 }