]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliFMDEnergyFitter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added event inspector and energy distribution fitter.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.cxx
CommitLineData
f8715167 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
17ClassImp(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//____________________________________________________________________
28namespace {
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//____________________________________________________________________
115AliFMDEnergyFitter::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//____________________________________________________________________
128AliFMDEnergyFitter::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//____________________________________________________________________
149AliFMDEnergyFitter::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//____________________________________________________________________
166AliFMDEnergyFitter::~AliFMDEnergyFitter()
167{
168 fRingHistos.Delete();
169}
170
171//____________________________________________________________________
172AliFMDEnergyFitter&
173AliFMDEnergyFitter::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//____________________________________________________________________
194AliFMDEnergyFitter::RingHistos*
195AliFMDEnergyFitter::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//____________________________________________________________________
209void
210AliFMDEnergyFitter::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//____________________________________________________________________
219Bool_t
220AliFMDEnergyFitter::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//____________________________________________________________________
256void
257AliFMDEnergyFitter::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//____________________________________________________________________
294void
295AliFMDEnergyFitter::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//____________________________________________________________________
308void
309AliFMDEnergyFitter::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//====================================================================
319AliFMDEnergyFitter::RingHistos::RingHistos()
320 : AliForwardUtil::RingHistos(),
321 fEDist(0),
322 fEmpty(0),
323 fEtaEDists(),
324 fList(0),
325 fDebug(0)
326{}
327
328//____________________________________________________________________
329AliFMDEnergyFitter::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//____________________________________________________________________
340AliFMDEnergyFitter::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//____________________________________________________________________
360AliFMDEnergyFitter::RingHistos&
361AliFMDEnergyFitter::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//____________________________________________________________________
387AliFMDEnergyFitter::RingHistos::~RingHistos()
388{
389 if (fEDist) delete fEDist;
390 fEtaEDists.Delete();
391}
392
393//____________________________________________________________________
394void
395AliFMDEnergyFitter::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//__________________________________________________________________
412TArrayD
66b34429 413AliFMDEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t n, Double_t min,
414 Double_t max) const
f8715167 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;
66b34429 427 Double_t next = bins[i-1] + dx + dI * dI * dx;
f8715167 428 bins[i] = next;
429 if (next > max) break;
430 }
431 bins.Set(i+1);
432 return bins;
433}
434
435//____________________________________________________________________
436void
437AliFMDEnergyFitter::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//____________________________________________________________________
463TH1D*
464AliFMDEnergyFitter::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//____________________________________________________________________
483TH1D*
484AliFMDEnergyFitter::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//____________________________________________________________________
513void
514AliFMDEnergyFitter::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//____________________________________________________________________
535TObjArray*
536AliFMDEnergyFitter::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;
66b34429 560 TH1* hS = 0;
f8715167 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//____________________________________________________________________
646TF1*
647AliFMDEnergyFitter::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//____________________________________________________________________
690TF1*
691AliFMDEnergyFitter::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//____________________________________________________________________
762TF1*
763AliFMDEnergyFitter::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//____________________________________________________________________
847Bool_t
848AliFMDEnergyFitter::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//____________________________________________________________________
888void
889AliFMDEnergyFitter::RingHistos::Output(TList* dir)
890{
891 fList = DefineOutputList(dir);
892}
893
894//____________________________________________________________________
895//
896// EOF
897//