New energy fitter that uses Landaus convolved with gaussians.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.cxx
1 //
2 // Class to do the energy correction of FMD ESD data
3 //
4 #include "AliFMDEnergyFitter.h"
5 #include <AliESDFMD.h>
6 #include <TAxis.h>
7 #include <TList.h>
8 #include <TH1.h>
9 #include <TF1.h>
10 #include <TMath.h>
11 #include "AliFMDAnaParameters.h"
12 #include <AliLog.h>
13 #include <TClonesArray.h>
14 #include <TFitResult.h>
15 #include <THStack.h>
16
17 ClassImp(AliFMDEnergyFitter)
18 #if 0
19 ; // This is for Emacs
20 #endif 
21
22 #define N_A(N)  (4+N-2)
23 #define N2_A(N) (4+(N-2)*3)
24 #define N2_D(N) (4+(N-2)*3+1)
25 #define N2_X(N) (4+(N-2)*3+2)
26
27 //____________________________________________________________________
28 namespace {
29   Double_t 
30   NLandau(Double_t* xp, Double_t* pp) 
31   {
32     Double_t  e        = xp[0];
33     Double_t  constant = pp[0];
34     Double_t  mpv      = pp[1];
35     Double_t  fwhm     = pp[2];
36     Int_t     n        = Int_t(pp[3]);
37     Double_t  result   = 0;
38     for (Int_t i = 1; i <= n; i++) {
39       Double_t mpvi  =  i*(mpv+fwhm*TMath::Log(i));
40       Double_t fwhmi =  i*fwhm;
41       Double_t ai    =  (i == 1 ? 1 : pp[N_A(i)]);
42       result         += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE);
43     }
44     result *= constant;
45     return result;
46   }
47
48   Double_t 
49   NLandau2(Double_t* xp, Double_t* pp) 
50   {
51     Double_t  e        = xp[0];
52     Double_t  constant = pp[0];
53     Double_t  mpv      = pp[1];
54     Double_t  fwhm     = pp[2];
55     Int_t     n        = Int_t(pp[3]);
56     Double_t  result   = TMath::Landau(e,mpv,fwhm,kTRUE);
57     for (Int_t i = 2; i <= n; i++) {
58       Double_t ai    =  pp[N2_A(i)];
59       Double_t mpvi  =  pp[N2_D(i)];
60       Double_t fwhmi =  pp[N2_X(i)];
61       result         += ai * TMath::Landau(e,mpvi,fwhmi,kTRUE);
62     }
63     result *= constant;
64     return result;
65   }
66
67   Double_t 
68   TripleLandau(Double_t *x, Double_t *par) 
69   {
70     Double_t energy   = x[0];
71     Double_t constant = par[0];
72     Double_t mpv      = par[1];
73     Double_t fwhm     = par[2];
74     Double_t alpha    = par[3];
75     Double_t beta     = par[4];
76     Double_t mpv2     = 2*mpv+2*fwhm*TMath::Log(2);
77     Double_t mpv3     = 3*mpv+3*fwhm*TMath::Log(3);
78
79     Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+
80                              alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE)+
81                              beta  * TMath::Landau(energy,mpv3,3*fwhm,kTRUE));
82   
83     return f;
84   }
85   Double_t 
86   DoubleLandau(Double_t *x, Double_t *par) 
87   {
88     Double_t energy   = x[0];
89     Double_t constant = par[0];
90     Double_t mpv      = par[1];
91     Double_t fwhm     = par[2];
92     Double_t alpha    = par[3];
93     Double_t mpv2     = 2*mpv+2*fwhm*TMath::Log(2);
94     
95     Double_t f = constant * (TMath::Landau(energy,mpv,fwhm,kTRUE)+
96                              alpha * TMath::Landau(energy,mpv2,2*fwhm,kTRUE));
97   
98     return f;
99   }
100   Double_t 
101   SingleLandau(Double_t *x, Double_t *par) 
102   {
103     Double_t energy   = x[0];
104     Double_t constant = par[0];
105     Double_t mpv      = par[1];
106     Double_t fwhm     = par[2];
107     
108     Double_t f = constant * TMath::Landau(energy,mpv,fwhm,kTRUE);
109   
110     return f;
111   }
112 }
113
114 //____________________________________________________________________
115 AliFMDEnergyFitter::AliFMDEnergyFitter()
116   : TNamed(), 
117     fRingHistos(),
118     fLowCut(0.3),
119     fNLandau(3),
120     fMinEntries(100),
121     fBinsToSubtract(4),
122     fDoFits(false),
123     fEtaAxis(),
124     fMaxE(10),
125     fNEbins(300), 
126     fUseIncreasingBins(true),
127     fDebug(0)
128 {}
129
130 //____________________________________________________________________
131 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
132   : TNamed("fmdEnergyFitter", title), 
133     fRingHistos(), 
134     fLowCut(0.3),
135     fNLandau(3),
136     fMinEntries(100),
137     fBinsToSubtract(4),
138     fDoFits(false),
139     fEtaAxis(0,0,0),
140     fMaxE(10),
141     fNEbins(300), 
142     fUseIncreasingBins(true),
143     fDebug(3)
144 {
145   fEtaAxis.SetName("etaAxis");
146   fEtaAxis.SetTitle("#eta");
147   fRingHistos.Add(new RingHistos(1, 'I'));
148   fRingHistos.Add(new RingHistos(2, 'I'));
149   fRingHistos.Add(new RingHistos(2, 'O'));
150   fRingHistos.Add(new RingHistos(3, 'I'));
151   fRingHistos.Add(new RingHistos(3, 'O'));
152 }
153
154 //____________________________________________________________________
155 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
156   : TNamed(o), 
157     fRingHistos(), 
158     fLowCut(o.fLowCut),
159     fNLandau(o.fNLandau),
160     fMinEntries(o.fMinEntries),
161     fBinsToSubtract(o.fBinsToSubtract),
162     fDoFits(o.fDoFits),
163     fEtaAxis(o.fEtaAxis),
164     fMaxE(o.fMaxE),
165     fNEbins(o.fNEbins), 
166     fUseIncreasingBins(o.fUseIncreasingBins),
167     fDebug(o.fDebug)
168 {
169   TIter    next(&o.fRingHistos);
170   TObject* obj = 0;
171   while ((obj = next())) fRingHistos.Add(obj);
172 }
173
174 //____________________________________________________________________
175 AliFMDEnergyFitter::~AliFMDEnergyFitter()
176 {
177   fRingHistos.Delete();
178 }
179
180 //____________________________________________________________________
181 AliFMDEnergyFitter&
182 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
183 {
184   TNamed::operator=(o);
185
186   fLowCut        = o.fLowCut;
187   fNLandau       = o.fNLandau;
188   fMinEntries    = o.fMinEntries;
189   fBinsToSubtract= o.fBinsToSubtract;
190   fDoFits        = o.fDoFits;
191   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
192   fDebug         = o.fDebug;
193
194   fRingHistos.Delete();
195   TIter    next(&o.fRingHistos);
196   TObject* obj = 0;
197   while ((obj = next())) fRingHistos.Add(obj);
198   
199   return *this;
200 }
201
202 //____________________________________________________________________
203 AliFMDEnergyFitter::RingHistos*
204 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
205 {
206   Int_t idx = -1;
207   switch (d) { 
208   case 1: idx = 0; break;
209   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
210   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
211   }
212   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
213   
214   return static_cast<RingHistos*>(fRingHistos.At(idx));
215 }
216
217 //____________________________________________________________________
218 void
219 AliFMDEnergyFitter::Init(const TAxis& eAxis)
220 {
221   if (fEtaAxis.GetNbins() == 0 || 
222       fEtaAxis.GetXmin() == fEtaAxis.GetXmax()) 
223     SetEtaAxis(eAxis);
224   TIter    next(&fRingHistos);
225   RingHistos* o = 0;
226   while ((o = static_cast<RingHistos*>(next())))
227     o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
228 }  
229 //____________________________________________________________________
230 void
231 AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
232 {
233   SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
234 }
235 //____________________________________________________________________
236 void
237 AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
238 {
239   fEtaAxis.Set(nBins,etaMin,etaMax);
240 }
241
242 //____________________________________________________________________
243 Bool_t
244 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, 
245                                Bool_t           empty)
246 {
247   for(UShort_t d = 1; d <= 3; d++) {
248     Int_t nRings = (d == 1 ? 1 : 2);
249     for (UShort_t q = 0; q < nRings; q++) {
250       Char_t      r    = (q == 0 ? 'I' : 'O');
251       UShort_t    nsec = (q == 0 ?  20 :  40);
252       UShort_t    nstr = (q == 0 ? 512 : 256);
253
254       RingHistos* histos = GetRingHistos(d, r);
255       
256       for(UShort_t s = 0; s < nsec;  s++) {
257         for(UShort_t t = 0; t < nstr; t++) {
258           Float_t mult = input.Multiplicity(d,r,s,t);
259           
260           // Keep dead-channel information. 
261           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0) 
262             continue;
263
264           // Get the pseudo-rapidity 
265           Double_t eta1 = input.Eta(d,r,s,t);
266           Int_t ieta = fEtaAxis.FindBin(eta1);
267           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 
268
269           histos->Fill(empty, ieta-1, mult);
270
271         } // for strip
272       } // for sector
273     } // for ring 
274   } // for detector
275
276   return kTRUE;
277 }
278
279 //____________________________________________________________________
280 void
281 AliFMDEnergyFitter::Fit(TList* dir)
282 {
283   if (!fDoFits) return;
284
285   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
286   if (!d) return;
287
288   // +1          for chi^2
289   // +3          for 1 landau 
290   // +1          for N 
291   // +fNLandau-1 for weights 
292   Int_t nStack = 1+3+1+fNLandau-1;
293   THStack* stack[nStack]; 
294   stack[0] = new THStack("chi2", "#chi^{2}/#nu");
295   stack[1] = new THStack("c",    "constant");
296   stack[2] = new THStack("mpv",  "#Delta_{p}");
297   stack[3] = new THStack("w",    "w");
298   stack[4] = new THStack("n",    "# of Landaus");
299   for (Int_t i = 2; i <= fNLandau; i++) 
300     stack[i-1+4] = new THStack(Form("a%d", i), 
301                                Form("Weight of %d signal", i));
302   for (Int_t i = 0; i < nStack; i++) 
303     d->Add(stack[i]);
304
305   TIter    next(&fRingHistos);
306   RingHistos* o = 0;
307   while ((o = static_cast<RingHistos*>(next()))) {
308     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNLandau,
309                           fMinEntries, fBinsToSubtract);
310     if (!l) continue;
311     for (Int_t i = 0; i < l->GetEntriesFast(); i++) { 
312       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i))); 
313     }
314   }
315 }
316
317 //____________________________________________________________________
318 void
319 AliFMDEnergyFitter::DefineOutput(TList* dir)
320 {
321   TList* d = new TList;
322   d->SetName(GetName());
323   dir->Add(d);
324   d->Add(&fEtaAxis);
325   TIter    next(&fRingHistos);
326   RingHistos* o = 0;
327   while ((o = static_cast<RingHistos*>(next()))) {
328     o->Output(d);
329   }
330 }
331 //____________________________________________________________________
332 void
333 AliFMDEnergyFitter::SetDebug(Int_t dbg)
334 {
335   fDebug = dbg;
336   TIter    next(&fRingHistos);
337   RingHistos* o = 0;
338   while ((o = static_cast<RingHistos*>(next())))
339     o->fDebug = dbg;
340 }  
341   
342 //====================================================================
343 AliFMDEnergyFitter::RingHistos::RingHistos()
344   : AliForwardUtil::RingHistos(), 
345     fEDist(0), 
346     fEmpty(0),
347     fEtaEDists(), 
348     fList(0),
349     fDebug(0)
350 {}
351
352 //____________________________________________________________________
353 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
354   : AliForwardUtil::RingHistos(d,r), 
355     fEDist(0), 
356     fEmpty(0),
357     fEtaEDists(), 
358     fList(0),
359     fDebug(0)
360 {
361   fEtaEDists.SetName("EDists");
362 }
363 //____________________________________________________________________
364 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
365   : AliForwardUtil::RingHistos(o), 
366     fEDist(o.fEDist), 
367     fEmpty(o.fEmpty),
368     fEtaEDists(), 
369     fList(0),
370     fDebug(0)
371 {
372   TIter next(&o.fEtaEDists);
373   TObject* obj = 0;
374   while ((obj = next())) fEtaEDists.Add(obj->Clone());
375   if (o.fList) {
376     fList = new TList;
377     fList->SetName(fName);
378     TIter next2(o.fList);
379     while ((obj = next2())) fList->Add(obj->Clone());
380   }
381 }
382
383 //____________________________________________________________________
384 AliFMDEnergyFitter::RingHistos&
385 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
386 {
387   AliForwardUtil::RingHistos::operator=(o);
388   
389   if (fEDist) delete fEDist;
390   if (fEmpty) delete fEmpty;
391   fEtaEDists.Delete();
392   if (fList) fList->Delete();
393
394   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
395   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
396   
397   TIter next(&o.fEtaEDists);
398   TObject* obj = 0;
399   while ((obj = next())) fEtaEDists.Add(obj->Clone());
400
401   if (o.fList) {
402     fList = new TList;
403     fList->SetName(fName);
404     TIter next2(o.fList);
405     while ((obj = next2())) fList->Add(obj->Clone());
406   }
407
408   return *this;
409 }
410 //____________________________________________________________________
411 AliFMDEnergyFitter::RingHistos::~RingHistos()
412 {
413   if (fEDist) delete fEDist;
414   fEtaEDists.Delete();
415 }
416
417 //____________________________________________________________________
418 void
419 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
420 {
421   if (empty) { 
422     fEmpty->Fill(mult);
423     return;
424   }
425   fEDist->Fill(mult);
426   
427   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
428   
429   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
430   if (!h) return;
431
432   h->Fill(mult);
433 }
434
435 //__________________________________________________________________
436 TArrayD
437 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, 
438                                                    Double_t max) const
439 {
440   // Service function to define a logarithmic axis. 
441   // Parameters: 
442   //   n    Number of bins 
443   //   min  Minimum of axis 
444   //   max  Maximum of axis 
445   TArrayD bins(n+1);
446   Double_t dx = (max-min) / n;
447   bins[0]     = min;
448   Int_t    i  = 1;
449   for (i = 1; i < n+1; i++) {
450     Double_t dI   = float(i)/n;
451     Double_t next = bins[i-1] + dx + dI * dI * dx;
452     bins[i]       = next;
453     if (next > max) break;
454   }
455   bins.Set(i+1);
456   return bins;
457 }
458
459 //____________________________________________________________________
460 void
461 AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax,
462                                      Double_t deMax, Int_t nDeBins, 
463                                      Bool_t incr)
464 {
465   TH1D* h = 0;
466   TString name  = Form("%s_etabin%03d", fName.Data(), ieta);
467   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
468                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
469   if (!incr) 
470     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
471   else { 
472     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
473     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
474   }
475     
476   h->SetDirectory(0);
477   h->SetXTitle("#DeltaE/#DeltaE_{mip}");        
478   h->SetFillColor(Color());
479   h->SetMarkerColor(Color());
480   h->SetLineColor(Color());
481   h->SetFillStyle(3001);
482   h->Sumw2();
483   
484   fEtaEDists.AddAt(h, ieta-1);
485 }
486 //____________________________________________________________________
487 TH1D*
488 AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
489                                         const char* title, 
490                                         const TAxis& eta) const
491 {
492   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
493                      Form("%s for %s", title, fName.Data()), 
494                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
495   h->SetXTitle("#eta");
496   h->SetYTitle(title);
497   h->SetDirectory(0);
498   h->SetFillColor(Color());
499   h->SetMarkerColor(Color());
500   h->SetLineColor(Color());
501   h->SetFillStyle(3001);
502   h->Sumw2();
503
504   return h;
505 }
506 //____________________________________________________________________
507 TH1D*
508 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, 
509                                           const char* title, 
510                                           const TAxis& eta, 
511                                           Int_t low, 
512                                           Int_t high, 
513                                           Double_t val, 
514                                           Double_t err) const
515 {
516   Double_t xlow  = eta.GetBinLowEdge(low);
517   Double_t xhigh = eta.GetBinUpEdge(high);
518   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
519                      Form("%s for %s", title, fName.Data()), 
520                      1, xlow, xhigh);
521   h->SetBinContent(1, val);
522   h->SetBinError(1, err);
523   h->SetXTitle("#eta");
524   h->SetYTitle(title);
525   h->SetDirectory(0);
526   h->SetFillColor(0);
527   h->SetMarkerColor(Color());
528   h->SetLineColor(Color());
529   h->SetFillStyle(0);
530   h->SetLineStyle(2);
531   h->SetLineWidth(2);
532
533   return h;
534 }
535                      
536 //____________________________________________________________________
537 void
538 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, 
539                                      Double_t maxDE, Int_t nDEbins, 
540                                      Bool_t useIncrBin)
541 {
542   fEDist = new TH1D(Form("%s_edist", fName.Data()), 
543                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 
544                     200, 0, 6);
545   fEmpty = new TH1D(Form("%s_empty", fName.Data()), 
546                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)", 
547                          fName.Data()), 200, 0, 6);
548   fList->Add(fEDist);
549   fList->Add(fEmpty);
550   // fEtaEDists.Expand(eAxis.GetNbins());
551   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { 
552     Double_t min = eAxis.GetBinLowEdge(i);
553     Double_t max = eAxis.GetBinUpEdge(i);
554     Make(i, min, max, maxDE, nDEbins, useIncrBin);
555   }
556   fList->Add(&fEtaEDists);
557 }
558 //____________________________________________________________________
559 TObjArray*
560 AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
561                                     Double_t lowCut, UShort_t nLandau,
562                                     UShort_t minEntries,
563                                     UShort_t minusBins) const
564 {
565   TList* l = GetOutputList(dir);
566   if (!l) return 0; 
567
568   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
569   if (!dists) { 
570     AliWarning(Form("Didn't find %s_EtaEDists in %s", 
571                     fName.Data(), l->GetName()));
572     l->ls();
573     return 0;
574   }
575
576   TObjArray* pars  = new TObjArray(3+nLandau+1);
577   pars->SetName("FitResults");
578   l->Add(pars);
579
580   TH1* hChi2 = 0;
581   TH1* hC    = 0;
582   TH1* hMpv  = 0;
583   TH1* hW    = 0;
584   TH1* hS    = 0;
585   TH1* hN    = 0;
586   TH1* hA[nLandau-1];
587   pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta));
588   pars->Add(hC    = MakePar("c",    "Constant", eta));
589   pars->Add(hMpv  = MakePar("mpv",  "#Delta_{p}", eta));
590   pars->Add(hW    = MakePar("w",    "#xi", eta));
591   pars->Add(hS    = MakePar("s",    "#sigma", eta));
592   pars->Add(hN    = MakePar("n",    "N", eta));
593   for (UShort_t i = 1; i < nLandau; i++) 
594     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
595
596   
597   Int_t nDists = dists->GetEntries();
598   Int_t low    = nDists;
599   Int_t high   = 0;
600   for (Int_t i = 0; i < nDists; i++) { 
601     TH1D* dist = static_cast<TH1D*>(dists->At(i));
602     if (!dist || dist->GetEntries() <= minEntries) continue;
603
604
605     TF1* res = FitHist(dist,lowCut,nLandau,minusBins);
606     if (!res) continue;
607     
608     low   = TMath::Min(low,i+1);
609     high  = TMath::Max(high,i+1);
610
611     Double_t chi2 = res->GetChisquare();
612     Int_t    ndf  = res->GetNDF();
613     hChi2->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
614     hC  ->SetBinContent(i+1, res->GetParameter(0));   
615     hMpv->SetBinContent(i+1, res->GetParameter(1)); 
616     hW  ->SetBinContent(i+1, res->GetParameter(2));   
617     hN  ->SetBinContent(i+1, res->GetParameter(3));   
618
619     hC  ->SetBinError(i+1, res->GetParError(1));
620     hMpv->SetBinError(i+1, res->GetParError(2));
621     hW  ->SetBinError(i+1, res->GetParError(2));
622     // hN  ->SetBinError(i, res->GetParError(3));
623
624     for (UShort_t j = 0; j < nLandau-1; j++) {
625       hA[j]->SetBinContent(i+1, res->GetParameter(4+j));
626       hA[j]->SetBinError(i+1, res->GetParError(4+j));
627     }
628   }
629
630   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
631   if (total && total->GetEntries() >= minEntries) { 
632     TF1* res = FitHist(total,lowCut,nLandau,minusBins);
633     if (res) { 
634       Double_t chi2 = res->GetChisquare();
635       Int_t    ndf  = res->GetNDF();
636       pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
637                           ndf > 0 ? chi2/ndf : 0, 0));
638       pars->Add(MakeTotal("t_c",    "Constant",     eta, low, high,
639                           res->GetParameter(0),res->GetParError(0)));
640       pars->Add(MakeTotal("t_mpv",  "#Delta_{p}",   eta, low, high,
641                           res->GetParameter(1),res->GetParError(1)));
642       pars->Add(MakeTotal("t_w",    "#xi",          eta, low, high,
643                           res->GetParameter(2),res->GetParError(2)));
644       pars->Add(MakeTotal("t_n",    "N",            eta, low, high,
645                           res->GetParameter(3),0));
646       for (UShort_t j = 1; j < nLandau; j++) 
647         pars->Add(MakeTotal(Form("a%d_t",j+1), 
648                             Form("a_{%d}",j+1), eta, low, high,
649                             res->GetParameter(3+j), res->GetParError(3+j)));
650     }
651   }
652     
653   TObjLink* lnk = dists->FirstLink();
654   while (lnk) {
655     TH1* h = static_cast<TH1*>(lnk->GetObject());
656     if (h->GetEntries() <= 0 || 
657         h->GetListOfFunctions()->GetEntries() <= 0) {
658       TObjLink* keep = lnk->Next();
659       dists->Remove(lnk);
660       lnk = keep;
661       continue;
662     }
663     lnk = lnk->Next();
664   }
665
666   return pars;
667 }
668
669 //____________________________________________________________________
670 TF1*
671 AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
672                                         Double_t lowCut, 
673                                         UShort_t nLandau, 
674                                         UShort_t minusBins) const
675 {
676   Double_t maxRange = 10;
677
678   AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins); 
679   f.Clear();
680   
681   // If we are only asked to fit a single particle, return this fit, 
682   // no matter what. 
683   if (nLandau == 1) {
684     TF1* r = f.Fit1Particle(dist, 0);
685     if (!r) return 0;
686     return new TF1(*r);
687   }
688
689   // Fit from 2 upto n particles  
690   for (Int_t i = 2; i <= nLandau; i++) f.FitNParticle(dist, i, 0);
691
692
693   // Now, we need to select the best fit 
694   Int_t nFits = f.fFitResults.GetEntriesFast();
695   TF1*  good[nFits];
696   for (Int_t i = nFits-1; i >= 0; i--) { 
697     if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i))))
698       good[i] = static_cast<TF1*>(f.fFunctions.At(i));
699   }
700   // If all else fails, use the 1 particle fit 
701   TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
702   for (Int_t i = nFits-1; i > 0; i--) {
703     if (!good[i]) continue;
704     ret = good[i];
705     break;
706   }
707   return new TF1(*ret);
708 }
709
710 //____________________________________________________________________
711 Bool_t
712 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r) const
713 {
714   Double_t chi2 = r->Chi2();
715   Int_t    ndf  = r->Ndf();
716   // Double_t prob = r.Prob();
717   if (ndf <= 0 || chi2 / ndf > 5) { 
718     if (fDebug > 2)
719       AliWarning(Form("Fit %s has chi^2/ndf=%f/%d=%f", 
720                       r->GetName(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf)));
721     return kFALSE;
722   }
723     
724   UShort_t nPar = r->NPar();
725   for (UShort_t i = 0; i < nPar; i++) { 
726     if (i == 3) continue; 
727     
728     Double_t v = r->Parameter(i);
729     Double_t e = r->ParError(i);
730     if (v == 0) continue;
731     if (v == 0 || e / v > 0.2) { 
732       if (fDebug > 2)
733         AliWarning(Form("Fit %s has Delta %s/%s=%f/%f=%f%%",
734                         r->GetName(), r->ParName(i).c_str(), 
735                         r->ParName(i).c_str(), e, v, 100*(v == 0 ? 0 : e/v)));
736       return kFALSE;
737     }
738   }
739   if (nPar > 5) { 
740     Double_t lastScale = r->Parameter(nPar-1);
741     if (lastScale <= 1e-7) { 
742       if (fDebug)
743         AliWarning(Form("Last scale %s is too small %f<1e-7", 
744                         r->ParName(nPar-1).c_str(), lastScale));
745       return kFALSE;
746     }
747   }
748   return kTRUE;
749 }
750
751
752 //____________________________________________________________________
753 void
754 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
755 {
756   fList = DefineOutputList(dir);
757 }
758
759 //____________________________________________________________________
760 //
761 // EOF
762 //