First go of energy loss spectrum algorithm.
[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     fDebug(0)
125 {}
126
127 //____________________________________________________________________
128 AliFMDEnergyFitter::AliFMDEnergyFitter(const char* title)
129   : TNamed("fmdEnergyFitter", title), 
130     fRingHistos(), 
131     fLowCut(0.3),
132     fNLandau(3),
133     fMinEntries(100),
134     fBinsToSubtract(4),
135     fDoFits(false),
136     fEtaAxis(100,-4,6),
137     fDebug(3)
138 {
139   fEtaAxis.SetName("etaAxis");
140   fEtaAxis.SetTitle("#eta");
141   fRingHistos.Add(new RingHistos(1, 'I'));
142   fRingHistos.Add(new RingHistos(2, 'I'));
143   fRingHistos.Add(new RingHistos(2, 'O'));
144   fRingHistos.Add(new RingHistos(3, 'I'));
145   fRingHistos.Add(new RingHistos(3, 'O'));
146 }
147
148 //____________________________________________________________________
149 AliFMDEnergyFitter::AliFMDEnergyFitter(const AliFMDEnergyFitter& o)
150   : TNamed(o), 
151     fRingHistos(), 
152     fLowCut(o.fLowCut),
153     fNLandau(o.fNLandau),
154     fMinEntries(o.fMinEntries),
155     fBinsToSubtract(o.fBinsToSubtract),
156     fDoFits(o.fDoFits),
157     fEtaAxis(o.fEtaAxis),
158     fDebug(o.fDebug)
159 {
160   TIter    next(&o.fRingHistos);
161   TObject* obj = 0;
162   while ((obj = next())) fRingHistos.Add(obj);
163 }
164
165 //____________________________________________________________________
166 AliFMDEnergyFitter::~AliFMDEnergyFitter()
167 {
168   fRingHistos.Delete();
169 }
170
171 //____________________________________________________________________
172 AliFMDEnergyFitter&
173 AliFMDEnergyFitter::operator=(const AliFMDEnergyFitter& o)
174 {
175   TNamed::operator=(o);
176
177   fLowCut        = o.fLowCut;
178   fNLandau       = o.fNLandau;
179   fMinEntries    = o.fMinEntries;
180   fBinsToSubtract= o.fBinsToSubtract;
181   fDoFits        = o.fDoFits;
182   fEtaAxis.Set(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax());
183   fDebug         = o.fDebug;
184
185   fRingHistos.Delete();
186   TIter    next(&o.fRingHistos);
187   TObject* obj = 0;
188   while ((obj = next())) fRingHistos.Add(obj);
189   
190   return *this;
191 }
192
193 //____________________________________________________________________
194 AliFMDEnergyFitter::RingHistos*
195 AliFMDEnergyFitter::GetRingHistos(UShort_t d, Char_t r) const
196 {
197   Int_t idx = -1;
198   switch (d) { 
199   case 1: idx = 0; break;
200   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
201   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
202   }
203   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
204   
205   return static_cast<RingHistos*>(fRingHistos.At(idx));
206 }
207
208 //____________________________________________________________________
209 void
210 AliFMDEnergyFitter::Init(const TAxis& eAxis)
211 {
212   fEtaAxis.Set(eAxis.GetNbins(),eAxis.GetXmin(),eAxis.GetXmax());
213   TIter    next(&fRingHistos);
214   RingHistos* o = 0;
215   while ((o = static_cast<RingHistos*>(next())))
216     o->Init(eAxis);
217 }  
218 //____________________________________________________________________
219 Bool_t
220 AliFMDEnergyFitter::Accumulate(const AliESDFMD& input, 
221                                Bool_t           empty)
222 {
223   for(UShort_t d = 1; d <= 3; d++) {
224     Int_t nRings = (d == 1 ? 1 : 2);
225     for (UShort_t q = 0; q < nRings; q++) {
226       Char_t      r    = (q == 0 ? 'I' : 'O');
227       UShort_t    nsec = (q == 0 ?  20 :  40);
228       UShort_t    nstr = (q == 0 ? 512 : 256);
229
230       RingHistos* histos = GetRingHistos(d, r);
231       
232       for(UShort_t s = 0; s < nsec;  s++) {
233         for(UShort_t t = 0; t < nstr; t++) {
234           Float_t mult = input.Multiplicity(d,r,s,t);
235           
236           // Keep dead-channel information. 
237           if(mult == AliESDFMD::kInvalidMult || mult > 10 || mult <= 0) 
238             continue;
239
240           // Get the pseudo-rapidity 
241           Double_t eta1 = input.Eta(d,r,s,t);
242           Int_t ieta = fEtaAxis.FindBin(eta1);
243           if (ieta < 1 || ieta >  fEtaAxis.GetNbins()) continue; 
244
245           histos->Fill(empty, ieta-1, mult);
246
247         } // for strip
248       } // for sector
249     } // for ring 
250   } // for detector
251
252   return kTRUE;
253 }
254
255 //____________________________________________________________________
256 void
257 AliFMDEnergyFitter::Fit(TList* dir)
258 {
259   if (!fDoFits) return;
260
261   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
262   if (!d) return;
263
264   // +1          for chi^2
265   // +3          for 1 landau 
266   // +1          for N 
267   // +fNLandau-1 for weights 
268   Int_t nStack = 1+3+1+fNLandau-1;
269   THStack* stack[nStack]; 
270   stack[0] = new THStack("chi2", "#chi^{2}/#nu");
271   stack[1] = new THStack("c",    "constant");
272   stack[2] = new THStack("mpv",  "#Delta_{p}");
273   stack[3] = new THStack("w",    "w");
274   stack[4] = new THStack("n",    "# of Landaus");
275   for (Int_t i = 2; i <= fNLandau; i++) 
276     stack[i-1+4] = new THStack(Form("a%d", i), 
277                                Form("Weight of %d signal", i));
278   for (Int_t i = 0; i < nStack; i++) 
279     d->Add(stack[i]);
280
281   TIter    next(&fRingHistos);
282   RingHistos* o = 0;
283   while ((o = static_cast<RingHistos*>(next()))) {
284     TObjArray* l = o->Fit(d, fEtaAxis, fLowCut, fNLandau,
285                           fMinEntries, fBinsToSubtract);
286     if (!l) continue;
287     for (Int_t i = 0; i < l->GetEntriesFast(); i++) { 
288       stack[i % nStack]->Add(static_cast<TH1*>(l->At(i))); 
289     }
290   }
291 }
292
293 //____________________________________________________________________
294 void
295 AliFMDEnergyFitter::DefineOutput(TList* dir)
296 {
297   TList* d = new TList;
298   d->SetName(GetName());
299   dir->Add(d);
300   d->Add(&fEtaAxis);
301   TIter    next(&fRingHistos);
302   RingHistos* o = 0;
303   while ((o = static_cast<RingHistos*>(next()))) {
304     o->Output(d);
305   }
306 }
307 //____________________________________________________________________
308 void
309 AliFMDEnergyFitter::SetDebug(Int_t dbg)
310 {
311   fDebug = dbg;
312   TIter    next(&fRingHistos);
313   RingHistos* o = 0;
314   while ((o = static_cast<RingHistos*>(next())))
315     o->fDebug = dbg;
316 }  
317   
318 //====================================================================
319 AliFMDEnergyFitter::RingHistos::RingHistos()
320   : AliForwardUtil::RingHistos(), 
321     fEDist(0), 
322     fEmpty(0),
323     fEtaEDists(), 
324     fList(0),
325     fDebug(0)
326 {}
327
328 //____________________________________________________________________
329 AliFMDEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
330   : AliForwardUtil::RingHistos(d,r), 
331     fEDist(0), 
332     fEmpty(0),
333     fEtaEDists(), 
334     fList(0),
335     fDebug(0)
336 {
337   fEtaEDists.SetName("EDists");
338 }
339 //____________________________________________________________________
340 AliFMDEnergyFitter::RingHistos::RingHistos(const RingHistos& o)
341   : AliForwardUtil::RingHistos(o), 
342     fEDist(o.fEDist), 
343     fEmpty(o.fEmpty),
344     fEtaEDists(), 
345     fList(0),
346     fDebug(0)
347 {
348   TIter next(&o.fEtaEDists);
349   TObject* obj = 0;
350   while ((obj = next())) fEtaEDists.Add(obj->Clone());
351   if (o.fList) {
352     fList = new TList;
353     fList->SetName(fName);
354     TIter next2(o.fList);
355     while ((obj = next2())) fList->Add(obj->Clone());
356   }
357 }
358
359 //____________________________________________________________________
360 AliFMDEnergyFitter::RingHistos&
361 AliFMDEnergyFitter::RingHistos::operator=(const RingHistos& o)
362 {
363   AliForwardUtil::RingHistos::operator=(o);
364   
365   if (fEDist) delete fEDist;
366   if (fEmpty) delete fEmpty;
367   fEtaEDists.Delete();
368   if (fList) fList->Delete();
369
370   fEDist = static_cast<TH1D*>(o.fEDist->Clone());
371   fEmpty = static_cast<TH1D*>(o.fEmpty->Clone());
372   
373   TIter next(&o.fEtaEDists);
374   TObject* obj = 0;
375   while ((obj = next())) fEtaEDists.Add(obj->Clone());
376
377   if (o.fList) {
378     fList = new TList;
379     fList->SetName(fName);
380     TIter next2(o.fList);
381     while ((obj = next2())) fList->Add(obj->Clone());
382   }
383
384   return *this;
385 }
386 //____________________________________________________________________
387 AliFMDEnergyFitter::RingHistos::~RingHistos()
388 {
389   if (fEDist) delete fEDist;
390   fEtaEDists.Delete();
391 }
392
393 //____________________________________________________________________
394 void
395 AliFMDEnergyFitter::RingHistos::Fill(Bool_t empty, Int_t ieta, Double_t mult)
396 {
397   if (empty) { 
398     fEmpty->Fill(mult);
399     return;
400   }
401   fEDist->Fill(mult);
402   
403   if (ieta < 0 || ieta >= fEtaEDists.GetEntries()) return;
404   
405   TH1D* h = static_cast<TH1D*>(fEtaEDists.At(ieta));
406   if (!h) return;
407
408   h->Fill(mult);
409 }
410
411 //__________________________________________________________________
412 TArrayD
413 AliFMDEnergyFitter::RingHistos::::MakeIncreasingAxis(Int_t n, Double_t min, 
414                                                      Double_t max) 
415 {
416   // Service function to define a logarithmic axis. 
417   // Parameters: 
418   //   n    Number of bins 
419   //   min  Minimum of axis 
420   //   max  Maximum of axis 
421   TArrayD bins(n+1);
422   Double_t dx = (max-min) / n;
423   bins[0]     = min;
424   Int_t    i  = 1;
425   for (i = 1; i < n+1; i++) {
426     Double_t dI   = float(i)/n;
427     Double_t next = bins[i-1] + dX1 + dI * dI * dX1;
428     bins[i]       = next;
429     if (next > max) break;
430   }
431   bins.Set(i+1);
432   return bins;
433 }
434
435 //____________________________________________________________________
436 void
437 AliFMDEnergyFitter::RingHistos::Make(Int_t ieta, Double_t emin, Double_t emax,
438                                      Double_t deMax, Int_t nDeBins, 
439                                      Bool_t incr)
440 {
441   TH1D* h = 0;
442   TString name  = Form("%s_etabin%03d", fName.Data(), ieta);
443   TString title = Form("#DeltaE/#DeltaE_{mip} for %s %+5.3f#leq#eta<%+5.3f "
444                        "(#eta bin %d)", fName.Data(), emin, emax, ieta);
445   if (!incr) 
446     h = new TH1D(name.Data(), title.Data(), nDeBins, 0, deMax);
447   else { 
448     TArrayD deAxis = MakeIncreasingAxis(nDeBins, 0, deMax);
449     h = new TH1D(name.Data(), title.Data(), deAxis.fN-1, deAxis.fArray);
450   }
451     
452   h->SetDirectory(0);
453   h->SetXTitle("#DeltaE/#DeltaE_{mip}");        
454   h->SetFillColor(Color());
455   h->SetMarkerColor(Color());
456   h->SetLineColor(Color());
457   h->SetFillStyle(3001);
458   h->Sumw2();
459   
460   fEtaEDists.AddAt(h, ieta-1);
461 }
462 //____________________________________________________________________
463 TH1D*
464 AliFMDEnergyFitter::RingHistos::MakePar(const char* name, 
465                                         const char* title, 
466                                         const TAxis& eta) const
467 {
468   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
469                      Form("%s for %s", title, fName.Data()), 
470                      eta.GetNbins(), eta.GetXmin(), eta.GetXmax());
471   h->SetXTitle("#eta");
472   h->SetYTitle(title);
473   h->SetDirectory(0);
474   h->SetFillColor(Color());
475   h->SetMarkerColor(Color());
476   h->SetLineColor(Color());
477   h->SetFillStyle(3001);
478   h->Sumw2();
479
480   return h;
481 }
482 //____________________________________________________________________
483 TH1D*
484 AliFMDEnergyFitter::RingHistos::MakeTotal(const char* name, 
485                                           const char* title, 
486                                           const TAxis& eta, 
487                                           Int_t low, 
488                                           Int_t high, 
489                                           Double_t val, 
490                                           Double_t err) const
491 {
492   Double_t xlow  = eta.GetBinLowEdge(low);
493   Double_t xhigh = eta.GetBinUpEdge(high);
494   TH1D* h = new TH1D(Form("%s_%s", fName.Data(), name), 
495                      Form("%s for %s", title, fName.Data()), 
496                      1, xlow, xhigh);
497   h->SetBinContent(1, val);
498   h->SetBinError(1, err);
499   h->SetXTitle("#eta");
500   h->SetYTitle(title);
501   h->SetDirectory(0);
502   h->SetFillColor(0);
503   h->SetMarkerColor(Color());
504   h->SetLineColor(Color());
505   h->SetFillStyle(0);
506   h->SetLineStyle(2);
507   h->SetLineWidth(2);
508
509   return h;
510 }
511                      
512 //____________________________________________________________________
513 void
514 AliFMDEnergyFitter::RingHistos::Init(const TAxis& eAxis, 
515                                      Double_t maxDE, Int_t nDEbins, 
516                                      Bool_t useIncrBin)
517 {
518   fEDist = new TH1D(Form("%s_edist", fName.Data()), 
519                     Form("#DeltaE/#DeltaE_{mip} for %s", fName.Data()), 
520                     200, 0, 6);
521   fEmpty = new TH1D(Form("%s_empty", fName.Data()), 
522                     Form("#DeltaE/#DeltaE_{mip} for %s (empty events)", 
523                          fName.Data()), 200, 0, 6);
524   fList->Add(fEDist);
525   fList->Add(fEmpty);
526   // fEtaEDists.Expand(eAxis.GetNbins());
527   for (Int_t i = 1; i <= eAxis.GetNbins(); i++) { 
528     Double_t min = eAxis.GetBinLowEdge(i);
529     Double_t max = eAxis.GetBinUpEdge(i);
530     Make(i, min, max, maxDE, nDEbins, useIncrBin);
531   }
532   fList->Add(&fEtaEDists);
533 }
534 //____________________________________________________________________
535 TObjArray*
536 AliFMDEnergyFitter::RingHistos::Fit(TList* dir, const TAxis& eta,
537                                     Double_t lowCut, UShort_t nLandau,
538                                     UShort_t minEntries,
539                                     UShort_t minusBins) const
540 {
541   TList* l = GetOutputList(dir);
542   if (!l) return 0; 
543
544   TList* dists = static_cast<TList*>(l->FindObject("EDists"));
545   if (!dists) { 
546     AliWarning(Form("Didn't find %s_EtaEDists in %s", 
547                     fName.Data(), l->GetName()));
548     l->ls();
549     return 0;
550   }
551
552   TObjArray* pars  = new TObjArray(3+nLandau+1);
553   pars->SetName("FitResults");
554   l->Add(pars);
555
556   TH1* hChi2 = 0;
557   TH1* hC    = 0;
558   TH1* hMpv  = 0;
559   TH1* hW    = 0;
560   TH1* nS    = 0;
561   TH1* hN    = 0;
562   TH1* hA[nLandau-1];
563   pars->Add(hChi2 = MakePar("chi2", "#chi^{2}/#nu", eta));
564   pars->Add(hC    = MakePar("c",    "Constant", eta));
565   pars->Add(hMpv  = MakePar("mpv",  "#Delta_{p}", eta));
566   pars->Add(hW    = MakePar("w",    "#xi", eta));
567   pars->Add(hS    = MakePar("s",    "#sigma", eta));
568   pars->Add(hN    = MakePar("n",    "N", eta));
569   for (UShort_t i = 1; i < nLandau; i++) 
570     pars->Add(hA[i-1] = MakePar(Form("a%d",i+1), Form("a_{%d}",i+1), eta));
571
572   
573   Int_t nDists = dists->GetEntries();
574   Int_t low    = nDists;
575   Int_t high   = 0;
576   for (Int_t i = 0; i < nDists; i++) { 
577     TH1D* dist = static_cast<TH1D*>(dists->At(i));
578
579     if (!dist) continue;
580
581     TF1* res = FitHist(dist,lowCut,nLandau,minEntries,minusBins);
582     if (!res) continue;
583     
584     low   = TMath::Min(low,i+1);
585     high  = TMath::Max(high,i+1);
586
587     Double_t chi2 = res->GetChisquare();
588     Int_t    ndf  = res->GetNDF();
589     hChi2->SetBinContent(i+1, ndf > 0 ? chi2 / ndf : 0);
590     hC  ->SetBinContent(i+1, res->GetParameter(0));   
591     hMpv->SetBinContent(i+1, res->GetParameter(1)); 
592     hW  ->SetBinContent(i+1, res->GetParameter(2));   
593     hN  ->SetBinContent(i+1, res->GetParameter(3));   
594
595     hC  ->SetBinError(i+1, res->GetParError(1));
596     hMpv->SetBinError(i+1, res->GetParError(2));
597     hW  ->SetBinError(i+1, res->GetParError(2));
598     // hN  ->SetBinError(i, res->GetParError(3));
599
600     for (UShort_t j = 0; j < nLandau-1; j++) {
601       hA[j]->SetBinContent(i+1, res->GetParameter(4+j));
602       hA[j]->SetBinError(i+1, res->GetParError(4+j));
603     }
604   }
605
606   TH1* total = GetOutputHist(l, Form("%s_edist", fName.Data()));
607   if (total) { 
608     TF1* res = FitHist(total,lowCut,nLandau,minEntries,minusBins);
609     if (res) { 
610       Double_t chi2 = res->GetChisquare();
611       Int_t    ndf  = res->GetNDF();
612       pars->Add(MakeTotal("t_chi2", "#chi^{2}/#nu", eta, low, high,
613                           ndf > 0 ? chi2/ndf : 0, 0));
614       pars->Add(MakeTotal("t_c",    "Constant",     eta, low, high,
615                           res->GetParameter(0),res->GetParError(0)));
616       pars->Add(MakeTotal("t_mpv",  "#Delta_{p}",   eta, low, high,
617                           res->GetParameter(1),res->GetParError(1)));
618       pars->Add(MakeTotal("t_w",    "#xi",          eta, low, high,
619                           res->GetParameter(2),res->GetParError(2)));
620       pars->Add(MakeTotal("t_n",    "N",            eta, low, high,
621                           res->GetParameter(3),0));
622       for (UShort_t j = 1; j < nLandau; j++) 
623         pars->Add(MakeTotal(Form("a%d_t",j+1), 
624                             Form("a_{%d}",j+1), eta, low, high,
625                             res->GetParameter(3+j), res->GetParError(3+j)));
626     }
627   }
628     
629   TObjLink* lnk = dists->FirstLink();
630   while (lnk) {
631     TH1* h = static_cast<TH1*>(lnk->GetObject());
632     if (h->GetEntries() <= 0 || 
633         h->GetListOfFunctions()->GetEntries() <= 0) {
634       TObjLink* keep = lnk->Next();
635       dists->Remove(lnk);
636       lnk = keep;
637       continue;
638     }
639     lnk = lnk->Next();
640   }
641
642   return pars;
643 }
644
645 //____________________________________________________________________
646 TF1*
647 AliFMDEnergyFitter::RingHistos::FitHist(TH1* dist,Double_t lowCut, 
648                                         UShort_t nLandau, 
649                                         UShort_t minEntries,
650                                         UShort_t minusBins) const
651 {
652   Double_t maxRange = 10;
653
654   if (dist->GetEntries() <= minEntries) return 0;
655
656   // Find the fit range 
657   dist->GetXaxis()->SetRangeUser(lowCut, maxRange);
658   
659   // Normalize peak to 1 
660   Double_t max = dist->GetMaximum(); 
661   dist->Scale(1/max);
662   
663   // Get the bin with maximum 
664   Int_t    maxBin = dist->GetMaximumBin();
665   Double_t maxE   = dist->GetBinLowEdge(maxBin);
666   
667   // Get the low edge 
668   dist->GetXaxis()->SetRangeUser(lowCut, maxE);
669   Int_t    minBin = maxBin - minusBins; // dist->GetMinimumBin();
670   Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),lowCut);
671   Double_t maxEE  = dist->GetBinCenter(maxBin+2*minusBins);
672
673   // Restore the range 
674   dist->GetXaxis()->SetRangeUser(0, maxRange);
675   
676   // First do a single landau fit 
677   TF1*          landau1 = new TF1("landau1", "landau", minE, maxEE);
678   // landau1->SetParameters(1,1,1,1);
679   landau1->SetParNames("C","#Delta_{p}","#xi");
680   landau1->SetParLimits(1,minE,maxRange);
681   landau1->SetParLimits(2,0,maxRange);
682   landau1->SetLineColor(Color());
683   // Tight fit around peak - make sure we get that right. 
684   TFitResultPtr r = dist->Fit(landau1, "RQNS", "", minE, maxEE);
685
686   return FitMore2(dist, nLandau, *r, landau1, minE, maxRange);
687 }
688
689 //____________________________________________________________________
690 TF1*
691 AliFMDEnergyFitter::RingHistos::FitMore(TH1*        dist,
692                                         UShort_t    nLandau,
693                                         TFitResult& r, 
694                                         TF1*        landau1,
695                                         Double_t    minE,
696                                         Double_t    maxRange) const
697 {
698   static TClonesArray res("TFitResult");
699   static TObjArray funcs;
700   res.Clear();
701   funcs.SetOwner();
702   funcs.Clear();
703   Int_t nRes = 0;
704
705   //r->Print();
706   new(res[nRes++]) TFitResult(r);
707   funcs.AddAtAndExpand(landau1, 0);
708
709   // Now try to fit 
710   for (UShort_t n = 2; n <= nLandau; n++) { 
711     TFitResult* rr = static_cast<TFitResult*>(res.At(nRes-1));
712     if (!rr) { 
713       AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n));
714       return 0;
715     }
716     // Create the function object 
717     Double_t mpvi  = rr->Parameter(1);
718     Double_t wi    = rr->Parameter(2);
719     Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi;
720     TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,3+n);
721     landaui->SetLineStyle((n % 10)+1);
722     landaui->SetLineWidth(2);
723     landaui->SetLineColor(Color());
724     landaui->SetParameter(0, rr->Parameter(0));
725     landaui->SetParameter(1, rr->Parameter(1));
726     landaui->SetParameter(2, rr->Parameter(2));
727     landaui->SetParLimits(1, minE, maxRange);
728     landaui->SetParLimits(2,0,maxRange);
729     landaui->FixParameter(3, n);
730     landaui->SetParNames("C","#Delta_{p}","#xi", "N");
731     for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit 
732       landaui->SetParameter(N_A(i), rr->Parameter(N_A(i)));
733       landaui->SetParLimits(N_A(i), 0,1);
734       landaui->SetParName(i, Form("a_{%d}", i));
735     }
736     landaui->SetParameter(N_A(n), (n == 2 ? 0.05 : rr->Parameter(N_A(n-1))/5));
737     landaui->SetParLimits(N_A(n), 0, 1);
738     landaui->SetParName(N_A(n), Form("a_{%d}", n));
739     // landaui->Print();
740
741     TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi);
742     // tr->Print();
743     if (CheckResult(*tr)) { 
744       new(res[nRes++]) TFitResult(*tr);
745       funcs.AddAtAndExpand(landaui,n-1);
746       continue;
747     }
748     // Stop on bad fit 
749     break;
750   }
751   TF1* ret = 0;
752   if (funcs.At(nRes-1)) ret = static_cast<TF1*>(funcs.At(nRes-1)->Clone());
753   dist->GetListOfFunctions()->Add(ret);
754
755   res.Clear();
756   funcs.Clear();
757
758   return ret;
759 }
760
761 //____________________________________________________________________
762 TF1*
763 AliFMDEnergyFitter::RingHistos::FitMore2(TH1*        dist,
764                                          UShort_t    nLandau,
765                                          TFitResult& r, 
766                                          TF1*        landau1,
767                                          Double_t    minE,
768                                          Double_t    maxRange) const
769 {
770   static TClonesArray res("TFitResult");
771   static TObjArray funcs;
772   res.Clear();
773   funcs.SetOwner();
774   funcs.Clear();
775   Int_t nRes = 0;
776
777   //r->Print();
778   new(res[nRes++]) TFitResult(r);
779   funcs.AddAtAndExpand(landau1, 0);
780
781   // Now try to fit 
782   for (UShort_t n = 2; n <= nLandau; n++) { 
783     TFitResult* rr = static_cast<TFitResult*>(res.At(nRes-1));
784     if (!rr) { 
785       AliError(Form("No previous result (%p, %d) for n=%d", rr, nRes, n));
786       return 0;
787     }
788     // Create the function object 
789     Double_t mpvi  = rr->Parameter(1);
790     Double_t wi    = rr->Parameter(2);
791     Double_t maxEi = n*(mpvi+wi*TMath::Log(n))+2*n*wi;
792     TF1* landaui = new TF1(Form("landau%d", n), &NLandau, minE, maxEi,
793                            n*3+1);
794     landaui->SetLineStyle((n % 10)+1);
795     landaui->SetLineWidth(2);
796     landaui->SetLineColor(Color());
797     landaui->SetParameter(0, rr->Parameter(0));
798     landaui->SetParameter(1, rr->Parameter(1));
799     landaui->SetParameter(2, rr->Parameter(2));
800     landaui->SetParLimits(1, minE, maxRange);
801     landaui->SetParLimits(2,0,maxRange);
802     landaui->FixParameter(3, n);
803     landaui->SetParNames("C","#Delta_{p}","#xi", "N");
804     for (UShort_t i = 2; i < n; i++) {// Take parameters from last fit 
805       landaui->SetParameter(N2_A(i), rr->Parameter(N2_A(i)));
806       landaui->SetParameter(N2_D(i), rr->Parameter(N2_D(i)));
807       landaui->SetParameter(N2_X(i), rr->Parameter(N2_X(i)));
808       landaui->SetParLimits(N2_A(i), 0,1);
809       landaui->SetParLimits(N2_D(i), minE,maxEi);
810       landaui->SetParLimits(N2_X(i), 0,maxRange);
811       landaui->SetParName(N2_A(i), Form("a_{%d}", i));
812       landaui->SetParName(N2_D(i), Form("#Delta_{p,%d}", i));
813       landaui->SetParName(N2_X(i), Form("#xi_{%d}", i));
814     }
815     landaui->SetParameter(N2_A(n), n == 2 ? 0.05 : rr->Parameter(N2_A(n-1))/5);
816     landaui->SetParameter(N2_D(n), mpvi);
817     landaui->SetParameter(N2_X(n), wi);
818     landaui->SetParLimits(N2_A(n), 0,1);
819     landaui->SetParLimits(N2_D(n), minE,maxEi);
820     landaui->SetParLimits(N2_X(n), 0,maxRange);
821     landaui->SetParName(N2_A(n), Form("a_{%d}", n));
822     landaui->SetParName(N2_D(n), Form("#Delta_{p,%d}", n));
823     landaui->SetParName(N2_X(n), Form("#xi_{%d}", n));
824     if (fDebug > 2) landaui->Print();
825
826     TFitResultPtr tr = dist->Fit(landaui, "RNSQ", "", minE, maxEi);
827     if (fDebug > 2)  tr->Print();
828     if (CheckResult(*tr)) { 
829       new(res[nRes++]) TFitResult(*tr);
830       funcs.AddAtAndExpand(landaui,n-1);
831       continue;
832     }
833     // Stop on bad fit 
834     break;
835   }
836   TF1* ret = 0;
837   if (funcs.At(nRes-1)) ret = static_cast<TF1*>(funcs.At(nRes-1)->Clone());
838   dist->GetListOfFunctions()->Add(ret);
839
840   res.Clear();
841   funcs.Clear();
842
843   return ret;
844 }
845
846 //____________________________________________________________________
847 Bool_t
848 AliFMDEnergyFitter::RingHistos::CheckResult(TFitResult& r) const
849 {
850   Double_t chi2 = r.Chi2();
851   Int_t    ndf  = r.Ndf();
852   // Double_t prob = r.Prob();
853   if (ndf <= 0 || chi2 / ndf > 5) { 
854     if (fDebug > 2)
855       AliWarning(Form("Fit %s has chi^2/ndf=%f/%d=%f", 
856                       r.GetName(), chi2, ndf, (ndf<0 ? 0 : chi2/ndf)));
857     return kFALSE;
858   }
859     
860   UShort_t nPar = r.NPar();
861   for (UShort_t i = 0; i < nPar; i++) { 
862     if (i == 3) continue; 
863     
864     Double_t v = r.Parameter(i);
865     Double_t e = r.ParError(i);
866     if (v == 0) continue;
867     if (v == 0 || e / v > 0.2) { 
868       if (fDebug > 2)
869         AliWarning(Form("Fit %s has Delta %s/%s=%f/%f=%f%%",
870                         r.GetName(), r.ParName(i).c_str(), 
871                         r.ParName(i).c_str(), e, v, 100*(v == 0 ? 0 : e/v)));
872       return kFALSE;
873     }
874   }
875   if (nPar > 4) { 
876     if (r.Parameter(nPar-1) <= 1e-10) { 
877       if (fDebug)
878         AliWarning(Form("Last scale %s is too small %f<1e-10", 
879                         r.ParName(nPar-1).c_str(), r.Parameter(nPar-1)));
880       return kFALSE;
881     }
882   }
883   return kTRUE;
884 }
885
886
887 //____________________________________________________________________
888 void
889 AliFMDEnergyFitter::RingHistos::Output(TList* dir)
890 {
891   fList = DefineOutputList(dir);
892 }
893
894 //____________________________________________________________________
895 //
896 // EOF
897 //