]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Made member functions virtual
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDCorrELossFit.cxx
CommitLineData
0bd4b00f 1#include "AliFMDCorrELossFit.h"
2#include "AliForwardUtil.h"
3#include <TF1.h>
4#include <TBrowser.h>
5#include <TVirtualPad.h>
6#include <THStack.h>
7#include <TH1D.h>
8#include <AliLog.h>
9#include <TList.h>
10#include <iostream>
11#include <iomanip>
12
13Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .2;
14Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
15Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 20;
16
17//____________________________________________________________________
18AliFMDCorrELossFit::ELossFit::ELossFit()
19 : fN(0),
20 fNu(0),
21 fChi2(0),
22 fC(0),
23 fDelta(0),
24 fXi(0),
25 fSigma(0),
26 fSigmaN(0),
27 fA(0),
28 fEC(0),
29 fEDelta(0),
30 fEXi(0),
31 fESigma(0),
32 fESigmaN(0),
33 fEA(0),
34 fQuality(0),
35 fDet(0),
36 fRing('\0'),
37 fBin(0)
38{}
39//____________________________________________________________________
40AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
41 : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ?
42 f.GetParameter(AliForwardUtil::ELossFitter::kN) :
43 1),
44 fNu(f.GetNDF()),
45 fChi2(f.GetChisquare()),
46 fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
47 fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
48 fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
49 fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
50 fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
51 fA(0),
52 fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
53 fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
54 fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
55 fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
56 fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
57 fEA(0),
58 fQuality(quality),
59 fDet(0),
60 fRing('\0'),
61 fBin(0)
62{
63 if (fN <= 0) return;
64 fA = new Double_t[fN];
65 fEA = new Double_t[fN];
66 for (Int_t i = 0; i < fN-1; i++) {
67 fA[i] = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
68 fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
69 }
70 fA[fN-1] = -9999;
71 fEA[fN-1] = -9999;
72}
73
74//____________________________________________________________________
75AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality,UShort_t n,
76 Double_t chi2, UShort_t nu,
77 Double_t c, Double_t ec,
78 Double_t delta, Double_t edelta,
79 Double_t xi, Double_t exi,
80 Double_t sigma, Double_t esigma,
81 Double_t sigman, Double_t esigman,
82 Double_t* a, Double_t* ea)
83 : fN(n),
84 fNu(nu),
85 fChi2(chi2),
86 fC(c),
87 fDelta(delta),
88 fXi(xi),
89 fSigma(sigma),
90 fSigmaN(sigman),
91 fA(0),
92 fEC(ec),
93 fEDelta(edelta),
94 fEXi(exi),
95 fESigma(esigma),
96 fESigmaN(esigman),
97 fEA(0),
98 fQuality(quality),
99 fDet(0),
100 fRing('\0'),
101 fBin(0)
102{
103 if (fN <= 0) return;
104 fA = new Double_t[fN];
105 fEA = new Double_t[fN];
106 for (Int_t i = 0; i < fN-1; i++) {
107 fA[i] = a[i];
108 fEA[i] = ea[i];
109 }
110 fA[fN-1] = -9999;
111 fEA[fN-1] = -9999;
112}
113//____________________________________________________________________
114AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
115 : TObject(o),
116 fN(o.fN),
117 fNu(o.fNu),
118 fChi2(o.fChi2),
119 fC(o.fC),
120 fDelta(o.fDelta),
121 fXi(o.fXi),
122 fSigma(o.fSigma),
123 fSigmaN(o.fSigmaN),
124 fA(0),
125 fEC(o.fEC),
126 fEDelta(o.fEDelta),
127 fEXi(o.fEXi),
128 fESigma(o.fESigma),
129 fESigmaN(o.fESigmaN),
130 fEA(0),
131 fQuality(o.fQuality),
132 fDet(o.fDet),
133 fRing(o.fRing),
134 fBin(o.fBin)
135{
136 if (fN <= 0) return;
137 fA = new Double_t[fN];
138 fEA = new Double_t[fN];
139 for (Int_t i = 0; i < fN-1; i++) {
140 fA[i] = o.fA[i];
141 fEA[i] = o.fEA[i];
142 }
143 fA[fN-1] = -9999;
144 fEA[fN-1] = -9999;
145}
146
147//____________________________________________________________________
148AliFMDCorrELossFit::ELossFit&
149AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
150{
151 fN = o.fN;
152 fNu = o.fNu;
153 fChi2 = o.fChi2;
154 fC = o.fC;
155 fDelta = o.fDelta;
156 fXi = o.fXi;
157 fSigma = o.fSigma;
158 fSigmaN = o.fSigmaN;
159 fEC = o.fEC;
160 fEDelta = o.fEDelta;
161 fEXi = o.fEXi;
162 fESigma = o.fESigma;
163 fESigmaN = o.fESigmaN;
164 fQuality = o.fQuality;
165 fDet = o.fDet;
166 fRing = o.fRing;
167 fBin = o.fBin;
168 if (fA) delete [] fA;
169 if (fEA) delete [] fEA;
170 fA = 0;
171 fEA = 0;
172
173 if (fN <= 0) return *this;
174 fA = new Double_t[fN];
175 fEA = new Double_t[fN];
176 for (Int_t i = 0; i < fN; i++) {
177 fA[i] = o.fA[i];
178 fEA[i] = o.fEA[i];
179 }
180
181 return *this;
182}
183
184//____________________________________________________________________
185AliFMDCorrELossFit::ELossFit::~ELossFit()
186{
187 if (fA) delete[] fA;
188 if (fEA) delete[] fEA;
189}
190
191
192//____________________________________________________________________
193Int_t
194AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError,
195 Double_t leastWeight,
196 UShort_t maxN) const
197{
198 Int_t n = TMath::Min(maxN, UShort_t(fN-1));
199 Int_t m = 1;
200 // fN is one larger than we have data
201 for (Int_t i = 0; i < n-1; i++, m++) {
202 if (fA[i] < leastWeight) break;
203 if (fEA[i] / fA[i] > maxRelError) break;
204 }
205 return m;
206}
207
208//____________________________________________________________________
209Double_t
210AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x,
211 UShort_t maxN) const
212{
213 return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN,
214 TMath::Min(maxN, UShort_t(fN)), fA);
215}
216
217//____________________________________________________________________
218Double_t
219AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x,
220 UShort_t maxN) const
221{
222 UShort_t n = TMath::Min(maxN, UShort_t(fN-1));
223 Double_t num = 0;
224 Double_t den = 0;
225 for (Int_t i = 1; i <= n; i++) {
226 Double_t a = (i == 1 ? 1 : fA[i-1]);
227 if (fA[i-1] < 0) break;
228 Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
229 num += i * a * f;
230 den += a * f;
231 }
232 if (den <= 0) return 1;
233 return num / den;
234}
235
236
237#define OUTPAR(N,V,E) \
238 std::setprecision(9) << \
239 std::setw(12) << N << ": " << \
240 std::setw(14) << V << " +/- " << \
241 std::setw(14) << E << " (" << \
242 std::setprecision(-1) << \
243 std::setw(5) << 100*(V>0?E/V:1) << "%)\n"
244
245
246//____________________________________________________________________
247Int_t
248AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
249{
250 const ELossFit* other = static_cast<const ELossFit*>(o);
251 if (this->fQuality < other->fQuality) return -1;
252 if (this->fQuality > other->fQuality) return +1;
253 Double_t chi2nu = (fNu == 0 ? 1000 : fChi2 / fNu);
254 Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
255 if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
256 if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
257 return 0;
258}
259
260//____________________________________________________________________
261void
262AliFMDCorrELossFit::ELossFit::Print(Option_t*) const
263{
264 std::cout << GetName() << ":\n"
265 << " chi^2/nu = " << fChi2 << "/" << fNu << " = "
266 << (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
267 << " Quality: " << fQuality << "\n"
268 << " NParticles: " << fN << " (" << FindMaxWeight() << ")\n"
269 << OUTPAR("Delta", fDelta, fEDelta)
270 << OUTPAR("xi", fXi, fEXi)
271 << OUTPAR("sigma", fSigma, fESigma)
272 << OUTPAR("sigma_n", fSigmaN, fESigmaN);
273 for (Int_t i = 0; i < fN-1; i++)
274 std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
275 std::cout << std::flush;
276}
277
278//____________________________________________________________________
279const Char_t*
280AliFMDCorrELossFit::ELossFit::GetName() const
281{
282 return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
283}
284
285//____________________________________________________________________
286void
287AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
288{
289 // Draw this one
290 Draw(b ? b->GetDrawOption() : "comp");
291 gPad->SetLogy();
292 gPad->Update();
293}
294
295//____________________________________________________________________
296void
297AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
298{
299 TString opt(option);
300 opt.ToUpper();
301 bool comp = false;
302 if (opt.Contains("COMP")) {
303 opt.ReplaceAll("COMP","");
304 comp = true;
305 }
306 if (!opt.Contains("SAME")) {
307 gPad->Clear();
308 }
309
310 TObjArray cleanup;
311 TF1* tot = AliForwardUtil::MakeNLandauGaus(1,
312 fDelta, fXi,
313 fSigma, fSigmaN,
314 fN, fA,
315 0.01, 10);
316 tot->SetLineColor(kBlack);
317 tot->SetLineWidth(2);
318 tot->SetLineStyle(1);
319 tot->SetTitle(GetName());
320 Double_t max = tot->GetMaximum();
321
322 if (!opt.Contains("SAME")) {
323 TH1* frame = new TH1F(GetName(),
324 Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
325 100, 0, 10);
326 frame->SetMinimum(max/10000);
327 frame->SetMaximum(max*1.4);
328 frame->SetXTitle("#Delta / #Delta_{mip}");
329 frame->Draw();
330 opt.Append(" SAME");
331 }
332 tot->DrawCopy(opt.Data());
333 cleanup.Add(tot);
334
335 if (!comp) {
336 gPad->cd();
337 return;
338 }
339
340 Double_t min = max;
341 opt.Append(" same");
342 Int_t maxW = FindMaxWeight();
343 for (Int_t i=1; i <= fN; i++) {
344 TF1* f = AliForwardUtil::MakeILandauGaus((i == 1 ? 1 : fA[i-2]),
345 fDelta, fXi,
346 fSigma, fSigmaN,
347 i, 0.01, 10);
348 f->SetLineWidth(2);
349 f->SetLineStyle(i > maxW ? 2 : 1);
350 min = TMath::Min(f->GetMaximum(), min);
351 f->DrawCopy(opt.Data());
352 cleanup.Add(f);
353 }
354 min /= 100;
355 tot->GetHistogram()->SetMaximum(max);
356 tot->GetHistogram()->SetMinimum(min);
357 tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
358
359 gPad->cd();
360}
361
362
363//____________________________________________________________________
364#define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
365
366//____________________________________________________________________
367void
368AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
369 Double_t maxRelError,
370 Double_t leastWeight)
371{
372 Int_t qual = 0;
373 if (fNu > 0 && fChi2 / fNu < maxChi2nu) qual += 4;;
374 if (CHECKPAR(fDelta, fEDelta, maxRelError)) qual++;
375 if (CHECKPAR(fXi, fEXi, maxRelError)) qual++;
376 if (CHECKPAR(fSigma, fESigma, maxRelError)) qual++;
377 if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
378 qual += FindMaxWeight(1.5*maxRelError, leastWeight, fN);
379 fQuality = qual;
380}
381
382//____________________________________________________________________
383AliFMDCorrELossFit::AliFMDCorrELossFit()
384 : TObject(),
385 fRings(),
386 fEtaAxis(0,0,0),
387 fLowCut(0)
388{
389 fRings.SetOwner(kTRUE);
390 fEtaAxis.SetTitle("#eta");
391 fEtaAxis.SetName("etaAxis");
392 fRings.SetName("rings");
393}
394
395//____________________________________________________________________
396AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
397 : TObject(o),
398 fRings(o.fRings),
399 fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
400 fLowCut(0)
401{
402 fEtaAxis.SetTitle("#eta");
403 fEtaAxis.SetName("etaAxis");
404}
405
406//____________________________________________________________________
407AliFMDCorrELossFit::~AliFMDCorrELossFit()
408{
409 fRings.Clear();
410}
411
412//____________________________________________________________________
413AliFMDCorrELossFit&
414AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
415{
416 fRings = o.fRings;
417 fLowCut = o.fLowCut;
418 SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
419
420 return *this;
421}
422//____________________________________________________________________
423Int_t
424AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
425{
426 if (fEtaAxis.GetXmin() == fEtaAxis.GetXmax() || fEtaAxis.GetNbins() == 0) {
427 AliWarning("No eta axis defined");
428 return -1;
429 }
430 Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
431 if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
432 return bin;
433}
434
435//____________________________________________________________________
436Bool_t
437AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
438{
439 TObjArray* ringArray = GetOrMakeRingArray(d, r);
440 if (!ringArray) {
441 AliError(Form("Failed to make ring array for FMD%d%c", d, r));
442 return kFALSE;
443 }
444 if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
445 AliError(Form("bin=%d is out of range [%d,%d]",
446 etaBin, 1, fEtaAxis.GetNbins()));
447 return kFALSE;
448 }
449 // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
450 ringArray->AddAtAndExpand(fit, etaBin);
451 return kTRUE;
452}
453
454//____________________________________________________________________
455Bool_t
456AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
457{
458 Int_t bin = FindEtaBin(eta);
459 if (bin <= 0) {
460 AliError(Form("eta=%f is out of range [%f,%f]",
461 eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
462 return kFALSE;
463 }
464
465 return SetFit(d, r, bin, fit);
466}
467//____________________________________________________________________
468Bool_t
469AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r,
470 Double_t eta,
471 Int_t quality,UShort_t n,
472 Double_t chi2, UShort_t nu,
473 Double_t c, Double_t ec,
474 Double_t delta, Double_t edelta,
475 Double_t xi, Double_t exi,
476 Double_t sigma, Double_t esigma,
477 Double_t sigman, Double_t esigman,
478 Double_t* a, Double_t* ea)
479{
480 ELossFit* e = new ELossFit(quality, n,
481 chi2, nu,
482 c, ec,
483 delta, edelta,
484 xi, exi,
485 sigma, esigma,
486 sigman, esigman,
487 a, ea);
488 if (!SetFit(d, r, eta, e)) {
489 delete e;
490 return kFALSE;
491 }
492 return kTRUE;
493}
494//____________________________________________________________________
495Bool_t
496AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta,
497 Int_t quality, const TF1& f)
498{
499 ELossFit* e = new ELossFit(quality, f);
500 if (!SetFit(d, r, eta, e)) {
501 delete e;
502 return kFALSE;
503 }
504 return kTRUE;
505}
506//____________________________________________________________________
507AliFMDCorrELossFit::ELossFit*
508AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Int_t etabin) const
509{
510 TObjArray* ringArray = GetRingArray(d, r);
511 if (!ringArray) {
512 AliError(Form("Failed to make ring array for FMD%d%c", d, r));
513 return 0;
514 }
515 if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
516 AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
517 etabin, 1, fEtaAxis.GetNbins(), d, r));
518 return 0;
519 }
520 if (etabin > ringArray->GetEntriesFast()) {
521 AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
522 etabin, 1, ringArray->GetEntriesFast(), d, r));
523 return 0;
524 }
525 else if (etabin >= ringArray->GetEntriesFast()) {
526 AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, "
527 "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
528 etabin-1));
529 etabin--;
530 }
531 else if (!ringArray->At(etabin)) {
532 AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d",
533 etabin, d, r, etabin+1));
534 etabin++;
535 }
536 return static_cast<ELossFit*>(ringArray->At(etabin));
537}
538//____________________________________________________________________
539AliFMDCorrELossFit::ELossFit*
540AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Double_t eta) const
541{
542 Int_t etabin = FindEtaBin(eta);
543 return FindFit(d, r, etabin);
544}
545//____________________________________________________________________
546TObjArray*
547AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
548{
549 Int_t idx = -1;
550 switch (d) {
551 case 1: idx = 0; break;
552 case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
553 case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
554 }
555 if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
556 return static_cast<TObjArray*>(fRings.At(idx));
557}
558//____________________________________________________________________
559TObjArray*
560AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
561{
562 Int_t idx = -1;
563 switch (d) {
564 case 1: idx = 0; break;
565 case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
566 case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
567 }
568 if (idx < 0) return 0;
569 if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
570 TObjArray* a = new TObjArray(0);
571 // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
572 a->SetName(Form("FMD%d%c", d, r));
573 a->SetOwner();
574 fRings.AddAtAndExpand(a, idx);
575 }
576 return static_cast<TObjArray*>(fRings.At(idx));
577}
578
579namespace {
580 TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
581 Int_t color)
582 {
583 TH1D* h = new TH1D(Form("%s_%s", name, title),
584 Form("%s %s", name, title),
585 axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
586 h->SetDirectory(0);
587 h->SetMarkerStyle(20);
588 h->SetMarkerColor(color);
589 h->SetMarkerSize(0.5);
590 h->SetFillColor(color);
591 h->SetFillStyle(3001);
592 h->SetLineColor(color);
593 return h;
594 }
595}
596
597//____________________________________________________________________
598void
599AliFMDCorrELossFit::Draw(Option_t* option)
600{
601 TString opt(Form("nostack %s", option));
602 opt.ToLower();
603 Bool_t rel = (opt.Contains("rel"));
604 Bool_t err = (opt.Contains("err"));
605 if (rel) opt.ReplaceAll("rel","");
606 if (err) opt.ReplaceAll("err","");
607 Int_t nRings = fRings.GetEntriesFast();
608 UShort_t maxN = 0;
609 for (Int_t i = 0; i < nRings; i++) {
610 if (!fRings.At(i)) continue;
611 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
612 Int_t nFits = a->GetEntriesFast();
613
614 for (Int_t j = 0; j < nFits; j++) {
615 ELossFit* fit = static_cast<ELossFit*>(a->At(j));
616 if (!fit) continue;
617 maxN = TMath::Max(maxN, UShort_t(fit->fN));
618 }
619 }
620 // AliInfo(Form("Maximum N is %d", maxN));
621 Int_t nPad = 7+maxN-1; // 7 regular params, and maxN-1 weights
622 TVirtualPad* pad = gPad;
623 pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
624
625 THStack* chi2nu;
626 THStack* c;
627 THStack* delta;
628 THStack* xi;
629 THStack* sigma;
630 THStack* sigman;
631 THStack* n;
632 TList stacks;
633 stacks.AddAt(chi2nu= new THStack("chi2", "#chi^{2}/#nu"), 0);
634 stacks.AddAt(c = new THStack("c", "C"), 1);
635 stacks.AddAt(delta = new THStack("delta", "#Delta_{mp}"), 2);
636 stacks.AddAt(xi = new THStack("xi", "#xi"), 3);
637 stacks.AddAt(sigma = new THStack("sigma", "#sigma"), 4);
638 stacks.AddAt(sigman= new THStack("sigman", "#sigma_{n}"), 5);
639 stacks.AddAt(n = new THStack("n", "N"), 6);
640 for (Int_t i = 1; i <= maxN; i++) {
641 stacks.AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), 6+i);
642 }
643
644 for (Int_t i = 0; i < nRings; i++) {
645 if (!fRings.At(i)) continue;
646 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
647 Int_t nFits = a->GetEntriesFast();
648 Int_t color = 0;
649 switch (i) {
650 case 0: color = kRed+2; break;
651 case 1: color = kGreen+2; break;
652 case 2: color = kGreen-2; break;
653 case 3: color = kBlue+2; break;
654 case 4: color = kBlue-2; break;
655 }
656
657 TH1D* hChi = MakeHist(fEtaAxis,a->GetName(), "chi2", color);
658 TH1D* hC = MakeHist(fEtaAxis,a->GetName(), "c", color);
659 TH1D* hDelta = MakeHist(fEtaAxis,a->GetName(), "delta", color);
660 TH1D* hXi = MakeHist(fEtaAxis,a->GetName(), "xi", color);
661 TH1D* hSigma = MakeHist(fEtaAxis,a->GetName(), "sigma", color);
662 TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
663 TH1D* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
664 TH1D* hA[maxN];
665 const char* ho = (rel || !err ? "hist" : "e");
666 chi2nu->Add(hChi, "hist"); // 0
667 c ->Add(hC, ho); // 1
668 delta ->Add(hDelta, ho); // 2
669 xi ->Add(hXi, ho); // 3
670 sigma ->Add(hSigma, ho); // 4
671 sigman->Add(hSigmaN,ho); // 5
672 n ->Add(hN, "hist"); // 6
673 hChi->SetFillColor(color);
674 hChi->SetFillStyle(3001);
675 hN->SetFillColor(color);
676 hN->SetFillStyle(3001);
677
678 for (Int_t k = 1; k <= maxN; k++) {
679 hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
680 static_cast<THStack*>(stacks.At(6+k))->Add(hA[k-1], ho);
681 }
682
683 for (Int_t j = 0; j < nFits; j++) {
684 ELossFit* f = static_cast<ELossFit*>(a->At(j));
685 if (!f) continue;
686
687 Int_t b = f->fBin;
688 Int_t nW = f->FindMaxWeight();
689 hChi ->SetBinContent(b, (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
690 hN ->SetBinContent(b, nW);
691
692 if (rel) {
693 hC ->SetBinContent(b, f->fC >0 ?f->fEC /f->fC :0);
694 hDelta ->SetBinContent(b, f->fDelta >0 ?f->fEDelta /f->fDelta :0);
695 hXi ->SetBinContent(b, f->fXi >0 ?f->fEXi /f->fXi :0);
696 hSigma ->SetBinContent(b, f->fSigma >0 ?f->fESigma /f->fSigma :0);
697 hSigmaN->SetBinContent(b, f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0);
698 }
699 else {
700 hC ->SetBinContent(b, f->fC);
701 hDelta ->SetBinContent(b, f->fDelta);
702 hXi ->SetBinContent(b, f->fXi);
703 hSigma ->SetBinContent(b, f->fSigma);
704 hSigmaN->SetBinContent(b, f->fSigmaN);
705 hC ->SetBinError(b, f->fEC);
706 hDelta ->SetBinError(b, f->fEDelta);
707 hXi ->SetBinError(b, f->fEXi);
708 hSigma ->SetBinError(b, f->fESigma);
709 hSigmaN->SetBinError(b, f->fESigmaN);
710 }
711 for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
712 if (rel)
713 hA[k]->SetBinContent(b, f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0);
714 else {
715 hA[k]->SetBinContent(b, f->fA[k]);
716 hA[k]->SetBinError(b, f->fEA[k]);
717 }
718 }
719 }
720 }
721 Int_t nPad2 = (nPad+1) / 2;
722 for (Int_t i = 0; i < nPad; i++) {
723 Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
724 TVirtualPad* p = pad->cd(iPad);
725 p->SetLeftMargin(.15);
726 p->SetFillColor(0);
727 p->SetFillStyle(0);
728 p->SetGridx();
729 p->SetGridy();
730 if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
731
732 THStack* stack = static_cast<THStack*>(stacks.At(i));
733 // AliInfo(Form("Drawing %s (%d) in pad %d", stack->GetName(), i, iPad));
734 stack->Draw(opt.Data());
735
736 TString tit(stack->GetTitle());
737 if (rel && i != 0 && i != 6)
738 tit = Form("#delta %s/%s", tit.Data(), tit.Data());
739 TH1* hist = stack->GetHistogram();
740 TAxis* yaxis = hist->GetYaxis();
741 yaxis->SetTitle(tit.Data());
742 yaxis->SetTitleSize(0.15);
743 yaxis->SetLabelSize(0.08);
744 yaxis->SetTitleOffset(0.35);
745 yaxis->SetNdivisions(5);
746
747 TAxis* xaxis = stack->GetHistogram()->GetXaxis();
748 xaxis->SetTitle("#eta");
749 xaxis->SetTitleSize(0.15);
750 xaxis->SetLabelSize(0.08);
751 xaxis->SetTitleOffset(0.35);
752 xaxis->SetNdivisions(10);
753
754
755 if (!rel) {
756 switch (i) {
757 case 0: break; // chi^2/nu
758 case 1: break; // C
759 case 2: stack->SetMinimum(0.4); break; // Delta
760 case 3: stack->SetMinimum(0.02);break; // xi
761 case 4: stack->SetMinimum(0.05);break; // sigma
762 case 5: break; // sigmaN
763 case 6:
764 stack->SetMinimum(-.1);
765 stack->SetMaximum(maxN+.5);
766 break; // N
767 default: break; // Weights
768 }
769 }
770 stack->DrawClone(opt.Data());
771 }
772 pad->cd();
773}
774
775//____________________________________________________________________
776void
777AliFMDCorrELossFit::Print(Option_t* option) const
778{
779 TString opt(option);
780 opt.ToUpper();
781 Int_t nRings = fRings.GetEntriesFast();
782 bool recurse = opt.Contains("R");
783
784 std::cout << "Low cut in fit range: " << fLowCut << "\n"
785 << "Eta axis: " << fEtaAxis.GetNbins()
786 << " bins, range [" << fEtaAxis.GetXmin() << ","
787 << fEtaAxis.GetXmax() << "]" << std::endl;
788
789 for (Int_t i = 0; i < nRings; i++) {
790 if (!fRings.At(i)) continue;
791 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
792 Int_t nFits = a->GetEntriesFast();
793
794 std::cout << a->GetName() << " [" << nFits << " entries]"
795 << (recurse ? ":\n" : "\t");
796 Int_t min = fEtaAxis.GetNbins()+1;
797 Int_t max = 0;
798 for (Int_t j = 0; j < nFits; j++) {
799 if (!a->At(j)) continue;
800
801 min = TMath::Min(j, min);
802 max = TMath::Max(j, max);
803
804 if (recurse) {
805 std::cout << "Bin # " << j << "\t";
806 ELossFit* fit = static_cast<ELossFit*>(a->At(j));
807 fit->Print(option);
808 }
809 }
810 if (!recurse)
811 std::cout << " bin range: " << std::setw(3) << min
812 << "-" << std::setw(3) << max << " " << std::setw(3)
813 << (max-min+1) << " bins" << std::endl;
814 }
815}
816
817//____________________________________________________________________
818void
819AliFMDCorrELossFit::Browse(TBrowser* b)
820{
821 b->Add(&fRings);
822 b->Add(&fEtaAxis);
823}
824
825
826
827//____________________________________________________________________
828//
829// EOF
830//