More code clean up.
[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 "AliForwardCorrectionManager.h"
6 #include <AliESDFMD.h>
7 #include <TAxis.h>
8 #include <TList.h>
9 #include <TH1.h>
10 #include <TF1.h>
11 #include <TMath.h>
12 #include <AliLog.h>
13 #include <TClonesArray.h>
14 #include <TFitResult.h>
15 #include <THStack.h>
16 #include <TROOT.h>
17 #include <iostream>
18 #include <iomanip>
19
20 ClassImp(AliFMDEnergyFitter)
21 #if 0
22 ; // This is for Emacs
23 #endif 
24 namespace {
25   const char* fgkEDistFormat = "%s_etabin%03d";
26 }
27
28
29 //____________________________________________________________________
30 AliFMDEnergyFitter::AliFMDEnergyFitter()
31   : TNamed(), 
32     fRingHistos(),
33     fLowCut(0.3),
34     fNParticles(3),
35     fMinEntries(100),
36     fFitRangeBinWidth(4),
37     fDoFits(false),
38     fDoMakeObject(false),
39     fEtaAxis(),
40     fMaxE(10),
41     fNEbins(300), 
42     fUseIncreasingBins(true),
43     fMaxRelParError(.25),
44     fMaxChi2PerNDF(20), 
45     fMinWeight(1e-7),
46     fDebug(0)
47 {}
48
49 //____________________________________________________________________
50 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
51   : TNamed("fmdEnergyFitter", title), 
52     fRingHistos(), 
53     fLowCut(0.3),
54     fNParticles(3),
55     fMinEntries(100),
56     fFitRangeBinWidth(4),
57     fDoFits(false),
58     fDoMakeObject(false),
59     fEtaAxis(0,0,0),
60     fMaxE(10),
61     fNEbins(300), 
62     fUseIncreasingBins(true),
63     fMaxRelParError(.25),
64     fMaxChi2PerNDF(20),  
65     fMinWeight(1e-7),
66     fDebug(3)
67 {
68   fEtaAxis.SetName("etaAxis");
69   fEtaAxis.SetTitle("#eta");
70   fRingHistos.Add(new RingHistos(1, 'I'));
71   fRingHistos.Add(new RingHistos(2, 'I'));
72   fRingHistos.Add(new RingHistos(2, 'O'));
73   fRingHistos.Add(new RingHistos(3, 'I'));
74   fRingHistos.Add(new RingHistos(3, 'O'));
75 }
76
77 //____________________________________________________________________
78 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
79   : TNamed(o), 
80     fRingHistos(), 
81     fLowCut(o.fLowCut),
82     fNParticles(o.fNParticles),
83     fMinEntries(o.fMinEntries),
84     fFitRangeBinWidth(o.fFitRangeBinWidth),
85     fDoFits(o.fDoFits),
86     fDoMakeObject(o.fDoMakeObject),
87     fEtaAxis(o.fEtaAxis),
88     fMaxE(o.fMaxE),
89     fNEbins(o.fNEbins), 
90     fUseIncreasingBins(o.fUseIncreasingBins),
91     fMaxRelParError(o.fMaxRelParError),
92     fMaxChi2PerNDF(o.fMaxChi2PerNDF), 
93     fMinWeight(o.fMinWeight),
94     fDebug(o.fDebug)
95 {
96   TIter    next(&o.fRingHistos);
97   TObject* obj = 0;
98   while ((obj = next())) fRingHistos.Add(obj);
99 }
100
101 //____________________________________________________________________
102 AliFMDEnergyFitter::~AliFMDEnergyFitter()
103 {
104   fRingHistos.Delete();
105 }
106
107 //____________________________________________________________________
108 AliFMDEnergyFitter&
109 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
110 {
111   TNamed::operator=(o);
112
113   fLowCut        = o.fLowCut;
114   fNParticles       = o.fNParticles;
115   fMinEntries    = o.fMinEntries;
116   fFitRangeBinWidth= o.fFitRangeBinWidth;
117   fDoFits        = o.fDoFits;
118   fDoMakeObject  = o.fDoMakeObject;
119   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
120   fDebug         = o.fDebug;
121   fMaxRelParError= o.fMaxRelParError;
122   fMaxChi2PerNDF = o.fMaxChi2PerNDF;
123   fMinWeight     = o.fMinWeight;
124
125   fRingHistos.Delete();
126   TIter    next(&o.fRingHistos);
127   TObject* obj = 0;
128   while ((obj = next())) fRingHistos.Add(obj);
129   
130   return *this;
131 }
132
133 //____________________________________________________________________
134 AliFMDEnergyFitter::RingHistos*
135 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
136 {
137   Int_t idx = -1;
138   switch (d) { 
139   case 1: idx = 0; break;
140   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
141   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
142   }
143   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
144   
145   return static_cast<RingHistos*>(fRingHistos.At(idx));
146 }
147
148 //____________________________________________________________________
149 void
150 AliFMDEnergyFitter::Init(const TAxis& eAxis)
151 {
152   if (fEtaAxis.GetNbins() == 0 || 
153       fEtaAxis.GetXmin() == fEtaAxis.GetXmax()) 
154     SetEtaAxis(eAxis);
155   TIter    next(&fRingHistos);
156   RingHistos* o = 0;
157   while ((o = static_cast<RingHistos*>(next())))
158     o->Init(fEtaAxis, fMaxE, fNEbins, fUseIncreasingBins);
159 }  
160 //____________________________________________________________________
161 void
162 AliFMDEnergyFitter::SetEtaAxis(const TAxis& eAxis)
163 {
164   SetEtaAxis(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
165 }
166 //____________________________________________________________________
167 void
168 AliFMDEnergyFitter::SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax)
169 {
170   fEtaAxis.Set(nBins,etaMin,etaMax);
171 }
172
173 //____________________________________________________________________
174 Bool_t
175 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, 
176                                Bool_t           empty)
177 {
178   for(UShort_t d = 1; d <= 3; d++) {
179     Int_t nRings = (d == 1 ? 1 : 2);
180     for (UShort_t q = 0; q < nRings; q++) {
181       Char_t      r    = (q == 0 ? 'I' : 'O');
182       UShort_t    nsec = (q == 0 ?  20 :  40);
183       UShort_t    nstr = (q == 0 ? 512 : 256);
184
185       RingHistos* histos = GetRingHistos(d, r);
186       
187       for(UShort_t s = 0; s < nsec;  s++) {
188         for(UShort_t t = 0; t < nstr; t++) {
189           Float_t mult = input.Multiplicity(d,r,s,t);
190           
191           // Keep dead-channel information. 
192           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0) 
193             continue;
194
195           // Get the pseudo-rapidity 
196           Double_t eta1 = input.Eta(d,r,s,t);
197           Int_t ieta = fEtaAxis.FindBin(eta1);
198           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 
199
200           histos->Fill(empty, ieta-1, mult);
201
202         } // for strip
203       } // for sector
204     } // for ring 
205   } // for detector
206
207   return kTRUE;
208 }
209
210 //____________________________________________________________________
211 void
212 AliFMDEnergyFitter::Fit(TList* dir)
213 {
214   if (!fDoFits) return;
215
216   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
217   if (!d) return;
218
219   // +1          for chi^2
220   // +3          for 1 landau 
221   // +1          for N 
222   // +fNParticles-1 for weights 
223   Int_t nStack = kN+fNParticles;
224   THStack* stack[nStack]; 
225   stack[0]         = new THStack("chi2",   "#chi^{2}/#nu");
226   stack[kC     +1] = new THStack("c",      "Constant");
227   stack[kDelta +1] = new THStack("delta",  "#Delta_{p}");
228   stack[kXi    +1] = new THStack("xi",     "#xi");
229   stack[kSigma +1] = new THStack("sigma",  "#sigma");
230   stack[kSigmaN+1] = new THStack("sigman", "#sigma_{n}");
231   stack[kN     +1] = new THStack("n",      "# Particles");
232   for (Int_t i = 2; i <= fNParticles; i++) 
233     stack[kN+i] = new THStack(Form("a%d", i), Form("a_{%d}", i));
234   for (Int_t i = 0; i < nStack; i++) 
235     d->Add(stack[i]);
236
237   TIter    next(&fRingHistos);
238   RingHistos* o = 0;
239   while ((o = static_cast<RingHistos*>(next()))) {
240     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNParticles,
241                           fMinEntries, fFitRangeBinWidth,
242                           fMaxRelParError, fMaxChi2PerNDF);
243     if (!l) continue;
244     for (Int_t i = 0; i < l->GetEntriesFast(); i++) { 
245       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i))); 
246     }
247   }
248
249   if (!fDoMakeObject) return;
250
251   MakeCorrectionsObject(d);
252 }
253
254 //____________________________________________________________________
255 void
256 AliFMDEnergyFitter::MakeCorrectionsObject(TList* d)
257 {
258   AliForwardCorrectionManager& mgr = AliForwardCorrectionManager::Instance();
259     
260   AliFMDCorrELossFit* obj = new AliFMDCorrELossFit;
261   obj->SetEtaAxis(fEtaAxis);
262   obj->SetLowCut(fLowCut);
263
264   TIter    next(&fRingHistos);
265   RingHistos* o = 0;
266   while ((o = static_cast<RingHistos*>(next()))) {
267     o->FindBestFits(d, *obj, fEtaAxis, fMaxRelParError, 
268                     fMaxChi2PerNDF, fMinWeight);
269   }
270   
271   TString oName(mgr.GetObjectName(AliForwardCorrectionManager::kELossFits));
272   TString fileName(mgr.GetFilePath(AliForwardCorrectionManager::kELossFits));
273   AliInfo(Form("Object %s created in output - should be extracted and copied "
274                "to %s", oName.Data(), fileName.Data()));
275   d->Add(obj, oName.Data());
276 }
277
278 //____________________________________________________________________
279 void
280 AliFMDEnergyFitter::DefineOutput(TList* dir)
281 {
282   TList* d = new TList;
283   d->SetName(GetName());
284   dir->Add(d);
285   d->Add(&fEtaAxis);
286   TIter    next(&fRingHistos);
287   RingHistos* o = 0;
288   while ((o = static_cast<RingHistos*>(next()))) {
289     o->Output(d);
290   }
291 }
292 //____________________________________________________________________
293 void
294 AliFMDEnergyFitter::SetDebug(Int_t dbg)
295 {
296   fDebug = dbg;
297   TIter    next(&fRingHistos);
298   RingHistos* o = 0;
299   while ((o = static_cast<RingHistos*>(next())))
300     o->fDebug = dbg;
301 }  
302
303 //____________________________________________________________________
304 void
305 AliFMDEnergyFitter::Print(Option_t*) const
306 {
307   char ind[gROOT->GetDirLevel()+1];
308   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
309   ind[gROOT->GetDirLevel()] = '\0';
310   std::cout << ind << "AliFMDEnergyFitter: " << GetName() << '\n'
311             << ind << " Low cut:                " << fLowCut << " E/E_mip\n"
312             << ind << " Max(particles):         " << fNParticles << '\n'
313             << ind << " Min(entries):           " << fMinEntries << '\n'
314             << ind << " Fit range:              " 
315             << fFitRangeBinWidth << " bins\n"
316             << ind << " Make fits:              " 
317             << (fDoFits ? "yes\n" : "no\n")
318             << ind << " Make object:            " 
319             << (fDoMakeObject ? "yes\n" : "no\n")
320             << ind << " Max E/E_mip:            " << fMaxE << '\n'
321             << ind << " N bins:                 " << fNEbins << '\n'
322             << ind << " Increasing bins:        " 
323             << (fUseIncreasingBins ?"yes\n" : "no\n")
324             << ind << " max(delta p/p):         " << fMaxRelParError << '\n'
325             << ind << " max(chi^2/nu):          " << fMaxChi2PerNDF << '\n'
326             << ind << " min(a_i):               " << fMinWeight << std::endl;
327 }
328   
329 //====================================================================
330 AliFMDEnergyFitter::RingHistos::RingHistos()
331   : AliForwardUtil::RingHistos(), 
332     fEDist(0), 
333     fEmpty(0),
334     fEtaEDists(), 
335     fList(0),
336     fFits("AliFMDCorrELossFit::ELossFit"),
337     fDebug(0)
338 {}
339
340 //____________________________________________________________________
341 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
342   : AliForwardUtil::RingHistos(d,r), 
343     fEDist(0), 
344     fEmpty(0),
345     fEtaEDists(), 
346     fList(0),
347     fFits("AliFMDCorrELossFit::ELossFit"),
348     fDebug(0)
349 {
350   fEtaEDists.SetName("EDists");
351 }
352 //____________________________________________________________________
353 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
354   : AliForwardUtil::RingHistos(o), 
355     fEDist(o.fEDist), 
356     fEmpty(o.fEmpty),
357     fEtaEDists(), 
358     fList(0),
359     fFits("AliFMDCorrELossFit::ELossFit"),
360     fDebug(0)
361 {
362   fFits.Clear();
363   TIter next(&o.fEtaEDists);
364   TObject* obj = 0;
365   while ((obj = next())) fEtaEDists.Add(obj->Clone());
366   if (o.fList) {
367     fList = new TList;
368     fList->SetName(fName);
369     TIter next2(o.fList);
370     while ((obj = next2())) fList->Add(obj->Clone());
371   }
372 }
373
374 //____________________________________________________________________
375 AliFMDEnergyFitter::RingHistos&
376 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
377 {
378   AliForwardUtil::RingHistos::operator=(o);
379   
380   if (fEDist) delete fEDist;
381   if (fEmpty) delete fEmpty;
382   fEtaEDists.Delete();
383   if (fList) fList->Delete();
384   fFits.Clear();
385
386   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
387   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
388   
389   TIter next(&o.fEtaEDists);
390   TObject* obj = 0;
391   while ((obj = next())) fEtaEDists.Add(obj->Clone());
392
393   if (o.fList) {
394     fList = new TList;
395     fList->SetName(fName);
396     TIter next2(o.fList);
397     while ((obj = next2())) fList->Add(obj->Clone());
398   }
399
400   return *this;
401 }
402 //____________________________________________________________________
403 AliFMDEnergyFitter::RingHistos::~RingHistos()
404 {
405   if (fEDist) delete fEDist;
406   fEtaEDists.Delete();
407 }
408
409 //____________________________________________________________________
410 void
411 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
412 {
413   if (empty) { 
414     fEmpty->Fill(mult);
415     return;
416   }
417   fEDist->Fill(mult);
418   
419   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
420   
421   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
422   if (!h) return;
423
424   h->Fill(mult);
425 }
426
427 //__________________________________________________________________
428 TArrayD
429 AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min, 
430                                                    Double_t max) const
431 {
432   // Service function to define a logarithmic axis. 
433   // Parameters: 
434   //   n    Number of bins 
435   //   min  Minimum of axis 
436   //   max  Maximum of axis 
437   TArrayD bins(n+1);
438   Double_t dx = (max-min) / n;
439   bins[0]     = min;
440   Int_t    i  = 1;
441   for (i = 1; i < n+1; i++) {
442     Double_t dI   = float(i)/n;
443     Double_t next = bins[i-1] + dx + dI * dI * dx;
444     bins[i]       = next;
445     if (next > max) break;
446   }
447   bins.Set(i+1);
448   return bins;
449 }
450
451 //____________________________________________________________________
452 void
453 AliFMDEnergyFitter::RingHistos::Make(Int_t    ieta, 
454                                      Double_t emin, 
455                                      Double_t emax,
456                                      Double_t deMax, 
457                                      Int_t    nDeBins, 
458                                      Bool_t   incr)
459 {
460   TH1D* h = 0;
461   TString name  = Form(fgkEDistFormat, fName.Data(), ieta);
462   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
463                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
464   if (!incr) 
465     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
466   else { 
467     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
468     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
469   }
470     
471   h->SetDirectory(0);
472   h->SetXTitle("#DeltaE/#DeltaE_{mip}");        
473   h->SetFillColor(Color());
474   h->SetMarkerColor(Color());
475   h->SetLineColor(Color());
476   h->SetFillStyle(3001);
477   h->Sumw2();
478   
479   fEtaEDists.AddAt(h, ieta-1);
480 }
481 //____________________________________________________________________
482 TH1D* AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
483                                               const char* title, 
484                                               const TAxis& eta) const
485 {
486   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
487                      Form("%s for %s", title, fName.Data()), 
488                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
489   h->SetXTitle("#eta");
490   h->SetYTitle(title);
491   h->SetDirectory(0);
492   h->SetFillColor(Color());
493   h->SetMarkerColor(Color());
494   h->SetLineColor(Color());
495   h->SetFillStyle(3001);
496   h->Sumw2();
497
498   return h;
499 }
500 //____________________________________________________________________
501 TH1D*
502 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, 
503                                           const char* title, 
504                                           const TAxis& eta, 
505                                           Int_t low, 
506                                           Int_t high, 
507                                           Double_t val, 
508                                           Double_t err) const
509 {
510   Double_t xlow  = eta.GetBinLowEdge(low);
511   Double_t xhigh = eta.GetBinUpEdge(high);
512   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
513                      Form("%s for %s", title, fName.Data()), 
514                      1, xlow, xhigh);
515   h->SetBinContent(1, val);
516   h->SetBinError(1, err);
517   h->SetXTitle("#eta");
518   h->SetYTitle(title);
519   h->SetDirectory(0);
520   h->SetFillColor(0);
521   h->SetMarkerColor(Color());
522   h->SetLineColor(Color());
523   h->SetFillStyle(0);
524   h->SetLineStyle(2);
525   h->SetLineWidth(2);
526
527   return h;
528 }
529                      
530 //____________________________________________________________________
531 void
532 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, 
533                                      Double_t maxDE, Int_t nDEbins, 
534                                      Bool_t useIncrBin)
535 {
536   fEDist = new TH1D(Form("%s_edist", fName.Data()), 
537                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 
538                     200, 0, 6);
539   fEmpty = new TH1D(Form("%s_empty", fName.Data()), 
540                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)", 
541                          fName.Data()), 200, 0, 6);
542   fList->Add(fEDist);
543   fList->Add(fEmpty);
544   // fEtaEDists.Expand(eAxis.GetNbins());
545   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { 
546     Double_t min = eAxis.GetBinLowEdge(i);
547     Double_t max = eAxis.GetBinUpEdge(i);
548     Make(i, min, max, maxDE, nDEbins, useIncrBin);
549   }
550   fList->Add(&fEtaEDists);
551 }
552 //____________________________________________________________________
553 TObjArray*
554 AliFMDEnergyFitter::RingHistos::Fit(TList*           dir, 
555                                     const TAxis&     eta,
556                                     Double_t         lowCut, 
557                                     UShort_t         nParticles,
558                                     UShort_t         minEntries,
559                                     UShort_t         minusBins, 
560                                     Double_t         relErrorCut, 
561                                     Double_t         chi2nuCut) const
562 {
563   // Get our ring list from the output container 
564   TList* l = GetOutputList(dir);
565   if (!l) return 0; 
566
567   // Get the energy distributions from the output container 
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   // Container of the fit results histograms 
577   TObjArray* pars  = new TObjArray(3+nParticles+1);
578   pars->SetName("FitResults");
579   l->Add(pars);
580
581   // Result objects stored in separate list on the output 
582   TH1* hChi2   = 0;
583   TH1* hC      = 0;
584   TH1* hDelta  = 0;
585   TH1* hXi     = 0;
586   TH1* hSigma  = 0;
587   TH1* hSigmaN = 0;
588   TH1* hN      = 0;
589   TH1* hA[nParticles-1];
590   pars->Add(hChi2   = MakePar("chi2",   "#chi^{2}/#nu", eta));
591   pars->Add(hC      = MakePar("c",      "Constant",     eta));
592   pars->Add(hDelta  = MakePar("delta",  "#Delta_{p}",   eta));
593   pars->Add(hXi     = MakePar("xi",     "#xi",          eta));
594   pars->Add(hSigma  = MakePar("sigma",  "#sigma",       eta));
595   pars->Add(hSigmaN = MakePar("sigman", "#sigma_{n}",   eta));
596   pars->Add(hN      = MakePar("n",      "N", eta));
597   for (UShort_t i = 1; i < nParticles; i++) 
598     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
599
600   
601   Int_t nDists = dists->GetEntries();
602   Int_t low    = nDists;
603   Int_t high   = 0;
604   for (Int_t i = 0; i < nDists; i++) { 
605     TH1D* dist = static_cast<TH1D*>(dists->At(i));
606     // Ignore empty histograms altoghether 
607     if (!dist || dist->GetEntries() <= 0) continue; 
608
609     // Scale to the bin-width
610     dist->Scale(1., "width");
611     
612     // Normalize peak to 1 
613     Double_t max = dist->GetMaximum(); 
614     if (max <= 0) continue;
615     dist->Scale(1/max);
616     
617     // Check that we have enough entries 
618     if (dist->GetEntries() <= minEntries) { 
619       AliWarning(Form("Histogram at %3d (%s) has too few entries (%d <= %d)",
620                       i, dist->GetName(), Int_t(dist->GetEntries()), 
621                       minEntries));
622       continue;
623     }
624
625     // Now fit 
626     TF1* res = FitHist(dist,lowCut,nParticles,minusBins,
627                        relErrorCut,chi2nuCut);
628     if (!res) continue;
629     // dist->GetListOfFunctions()->Add(res);
630
631     // Store eta limits 
632     low   = TMath::Min(low,i+1);
633     high  = TMath::Max(high,i+1);
634
635     // Get the reduced chi square
636     Double_t chi2 = res->GetChisquare();
637     Int_t    ndf  = res->GetNDF();
638     
639     // Store results of best fit in output histograms 
640     res->SetLineWidth(4);
641     hChi2   ->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
642     hC      ->SetBinContent(i+1, res->GetParameter(kC));   
643     hDelta  ->SetBinContent(i+1, res->GetParameter(kDelta)); 
644     hXi     ->SetBinContent(i+1, res->GetParameter(kXi));   
645     hSigma  ->SetBinContent(i+1, res->GetParameter(kSigma));   
646     hSigmaN ->SetBinContent(i+1, res->GetParameter(kSigmaN));   
647     hN      ->SetBinContent(i+1, res->GetParameter(kN));   
648
649     hC     ->SetBinError(i+1, res->GetParError(kC));
650     hDelta ->SetBinError(i+1, res->GetParError(kDelta));
651     hXi    ->SetBinError(i+1, res->GetParError(kXi));
652     hSigma ->SetBinError(i+1, res->GetParError(kSigma));
653     hSigmaN->SetBinError(i+1, res->GetParError(kSigmaN));
654     hN     ->SetBinError(i+1, res->GetParError(kN));
655
656     for (UShort_t j = 0; j < nParticles-1; j++) {
657       hA[j]->SetBinContent(i+1, res->GetParameter(kA+j));
658       hA[j]->SetBinError(i+1, res->GetParError(kA+j));
659     }
660   }
661
662   // Fit the full-ring histogram 
663   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
664   if (total && total->GetEntries() >= minEntries) { 
665
666     // Scale to the bin-width
667     total->Scale(1., "width");
668     
669     // Normalize peak to 1 
670     Double_t max = total->GetMaximum(); 
671     if (max > 0) total->Scale(1/max);
672
673     TF1* res = FitHist(total,lowCut,nParticles,minusBins,
674                        relErrorCut,chi2nuCut);
675     if (res) { 
676       // Make histograms for the result of this fit 
677       Double_t chi2 = res->GetChisquare();
678       Int_t    ndf  = res->GetNDF();
679       pars->Add(MakeTotal("t_chi2",   "#chi^{2}/#nu", eta, low, high,
680                           ndf > 0 ? chi2/ndf : 0, 0));
681       pars->Add(MakeTotal("t_c",      "Constant",     eta, low, high,
682                           res->GetParameter(kC),
683                           res->GetParError(kC)));
684       pars->Add(MakeTotal("t_delta",  "#Delta_{p}",   eta, low, high,
685                           res->GetParameter(kDelta),
686                           res->GetParError(kDelta)));
687       pars->Add(MakeTotal("t_xi",     "#xi",          eta, low, high,
688                           res->GetParameter(kXi),
689                           res->GetParError(kXi)));
690       pars->Add(MakeTotal("t_sigma",  "#sigma",       eta, low, high,
691                           res->GetParameter(kSigma),
692                           res->GetParError(kSigma)));
693       pars->Add(MakeTotal("t_sigman", "#sigma_{n}",   eta, low, high,
694                           res->GetParameter(kSigmaN),
695                           res->GetParError(kSigmaN)));
696       pars->Add(MakeTotal("t_n",    "N",              eta, low, high,
697                           res->GetParameter(kN),0));
698       for (UShort_t j = 0; j < nParticles-1; j++) 
699         pars->Add(MakeTotal(Form("t_a%d",j+2), 
700                             Form("a_{%d}",j+2), eta, low, high,
701                             res->GetParameter(kA+j), 
702                             res->GetParError(kA+j)));
703     }
704   }
705     
706   // Clean up list of histogram.  Histograms with no entries or 
707   // no functions are deleted.  We have to do this using the TObjLink 
708   // objects stored in the list since ROOT cannot guaranty the validity 
709   // of iterators when removing from a list - tsck.  Should just implement
710   // TIter::Remove(). 
711   TObjLink* lnk = dists->FirstLink();
712   while (lnk) {
713     TH1* h = static_cast<TH1*>(lnk->GetObject());
714     bool remove = false;
715     if (h->GetEntries() <= 0) { 
716       // AliWarning(Form("No entries in %s - removing", h->GetName()));
717       remove = true;
718     }
719     else if (h->GetListOfFunctions()->GetEntries() <= 0) {
720       // AliWarning(Form("No fuctions associated with %s - removing",
721       //            h->GetName()));
722       remove = true;
723     }
724     if (remove) {
725       TObjLink* keep = lnk->Next();
726       dists->Remove(lnk);
727       lnk = keep;
728       continue;
729     }
730     lnk = lnk->Next();
731   }
732
733   return pars;
734 }
735
736 //____________________________________________________________________
737 TF1*
738 AliFMDEnergyFitter::RingHistos::FitHist(TH1*     dist,
739                                         Double_t lowCut, 
740                                         UShort_t nParticles, 
741                                         UShort_t minusBins, 
742                                         Double_t relErrorCut, 
743                                         Double_t chi2nuCut) const
744 {
745   Double_t maxRange = 10;
746
747   // Create a fitter object 
748   AliForwardUtil::ELossFitter f(lowCut, maxRange, minusBins); 
749   f.Clear();
750   
751   
752   // If we are only asked to fit a single particle, return this fit, 
753   // no matter what. 
754   if (nParticles == 1) {
755     TF1* r = f.Fit1Particle(dist, 0);
756     if (!r) return 0;
757     return new TF1(*r);
758   }
759
760   // Fit from 2 upto n particles  
761   for (Int_t i = 2; i <= nParticles; i++) f.FitNParticle(dist, i, 0);
762
763
764   // Now, we need to select the best fit 
765   Int_t nFits = f.fFitResults.GetEntriesFast();
766   TF1*  good[nFits];
767   for (Int_t i = nFits-1; i >= 0; i--) { 
768     good[i] = 0;
769     TF1* ff = static_cast<TF1*>(f.fFunctions.At(i));
770     // ff->SetLineColor(Color());
771     ff->SetRange(0, maxRange);
772     dist->GetListOfFunctions()->Add(new TF1(*ff));
773     if (CheckResult(static_cast<TFitResult*>(f.fFitResults.At(i)),
774                     relErrorCut, chi2nuCut)) {
775       good[i] = ff;
776       ff->SetLineWidth(2);
777       // f.fFitResults.At(i)->Print("V");
778     }
779   }
780   // If all else fails, use the 1 particle fit 
781   TF1* ret = static_cast<TF1*>(f.fFunctions.At(0));
782
783   // Find the fit with the most valid particles 
784   for (Int_t i = nFits-1; i > 0; i--) {
785     if (!good[i]) continue;
786     if (fDebug > 1) {
787       AliInfo(Form("Choosing fit with n=%d", i+1));
788       f.fFitResults.At(i)->Print();
789     }
790     ret = good[i];
791     break;
792   }
793   // Give a warning if we're using fall-back 
794   if (ret == f.fFunctions.At(0)) {
795     AliWarning("Choosing fall-back 1 particle fit");
796   }
797   // Copy our result and return (the functions are owned by the fitter)
798   TF1* fret = new TF1(*ret);
799   return fret;
800 }
801
802 //____________________________________________________________________
803 Bool_t
804 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult* r,
805                                             Double_t relErrorCut, 
806                                             Double_t chi2nuCut) const
807 {
808   if (fDebug > 10) r->Print();
809   TString  n    = r->GetName();
810   n.ReplaceAll("TFitResult-", "");
811   Double_t chi2 = r->Chi2();
812   Int_t    ndf  = r->Ndf();
813   // Double_t prob = r.Prob();
814   Bool_t ret = kTRUE;
815   
816   // Check that the reduced chi square isn't larger than 5
817   if (ndf <= 0 || chi2 / ndf > chi2nuCut) { 
818     if (fDebug > 2) {
819       AliWarning(Form("%s: chi^2/ndf=%12.5f/%3d=%12.5f>%12.5f", 
820                       n.Data(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf),
821                       chi2nuCut)); }
822     ret = kFALSE;
823   }
824     
825   // Check each parameter 
826   UShort_t nPar = r->NPar();
827   for (UShort_t i = 0; i < nPar; i++) { 
828     if (i == kN) continue;  // Skip the number parameter 
829     
830     // Get value and error and check value 
831     Double_t v  = r->Parameter(i);
832     Double_t e  = r->ParError(i);
833     if (v == 0) continue;
834
835     // Calculate the relative error and check it 
836     Double_t re  = e / v;
837     Double_t cut = relErrorCut * (i < kN ? 1 : 2);
838     if (re > cut) { 
839       if (fDebug > 2) {
840         AliWarning(Form("%s: Delta %s/%s=%9.5f/%9.5f=%5.1f%%>%5.1f%%",
841                         n.Data(), r->ParName(i).c_str(), 
842                         r->ParName(i).c_str(), e, v, 
843                         100*(v == 0 ? 0 : e/v),
844                         100*(cut))); }
845       ret = kFALSE;
846     }
847   }
848
849   // Check if we have scale parameters 
850   if (nPar > kN) { 
851     
852     // Check that the last particle has a significant contribution 
853     Double_t lastScale = r->Parameter(nPar-1);
854     if (lastScale <= 1e-7) { 
855       if (fDebug) {
856         AliWarning(Form("%s: %s=%9.6f<1e-7", 
857                         n.Data(), r->ParName(nPar-1).c_str(), lastScale)); }
858       ret = kFALSE;
859     }
860   }
861   return ret;
862 }
863
864
865 //__________________________________________________________________
866 void
867 AliFMDEnergyFitter::RingHistos::FindBestFits(TList*              d, 
868                                              AliFMDCorrELossFit& obj,
869                                              const TAxis&        eta,     
870                                              Double_t            relErrorCut, 
871                                              Double_t            chi2nuCut,
872                                              Double_t            minWeightCut)
873 {
874   // Get our ring list from the output container 
875   TList* l = GetOutputList(d);
876   if (!l) return; 
877
878   // Get the energy distributions from the output container 
879   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
880   if (!dists) { 
881     AliWarning(Form("Didn't find %s_EtaEDists in %s", 
882                     fName.Data(), l->GetName()));
883     l->ls();
884     return;
885   }
886   Int_t nBin = eta.GetNbins();
887
888   for (Int_t b = 1; b <= nBin; b++) { 
889     TString n(Form(fgkEDistFormat, fName.Data(), b));
890     TH1D*   dist = static_cast<TH1D*>(dists->FindObject(n));
891     // Ignore empty histograms altoghether 
892     if (!dist || dist->GetEntries() <= 0) continue; 
893     
894     AliFMDCorrELossFit::ELossFit* best = FindBestFit(dist,
895                                                      relErrorCut,
896                                                      chi2nuCut,  
897                                                      minWeightCut);
898     best->fDet  = fDet; 
899     best->fRing = fRing;
900     best->fBin  = b; // 
901     // Double_t eta = fAxis->GetBinCenter(b);
902     obj.SetFit(fDet, fRing, b, new AliFMDCorrELossFit::ELossFit(*best));
903   }
904 }
905
906 //__________________________________________________________________
907 AliFMDCorrELossFit::ELossFit* 
908 AliFMDEnergyFitter::RingHistos::FindBestFit(TH1* dist,
909                                             Double_t relErrorCut, 
910                                             Double_t chi2nuCut,
911                                             Double_t minWeightCut) 
912 {
913   TList* funcs = dist->GetListOfFunctions();
914   TIter  next(funcs);
915   TF1*   func = 0;
916   fFits.Clear();
917   Int_t  i = 0;
918   // Info("FindBestFit", "%s", dist->GetName());
919   while ((func = static_cast<TF1*>(next()))) { 
920     AliFMDCorrELossFit::ELossFit* fit = 
921       new(fFits[i++]) AliFMDCorrELossFit::ELossFit(0,*func);
922     fit->CalculateQuality(chi2nuCut, relErrorCut, minWeightCut);
923   }
924   fFits.Sort(false);
925   // fFits.Print();
926   return static_cast<AliFMDCorrELossFit::ELossFit*>(fFits.At(0));
927 }
928
929
930 //____________________________________________________________________
931 void
932 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
933 {
934   fList = DefineOutputList(dir);
935 }
936
937 //____________________________________________________________________
938 //
939 // EOF
940 //