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