]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/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
Fixes for energy loss generation
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDCorrELossFit.cxx
CommitLineData
7dacdf84 1// Object holding the Energy loss fit 'correction'
7984e5f7 2//
3// These are generated from Monte-Carlo or real ESDs.
7dacdf84 4//
0bd4b00f 5#include "AliFMDCorrELossFit.h"
6#include "AliForwardUtil.h"
7#include <TF1.h>
8#include <TBrowser.h>
9#include <TVirtualPad.h>
10#include <THStack.h>
8449e3e0 11#include <TLatex.h>
12#include <TLegend.h>
0bd4b00f 13#include <TH1D.h>
14#include <AliLog.h>
7dacdf84 15#include <TMath.h>
0bd4b00f 16#include <TList.h>
17#include <iostream>
18#include <iomanip>
19
81775aba 20Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .25;
21Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-7;
0bd4b00f 22Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 20;
23
24//____________________________________________________________________
25AliFMDCorrELossFit::ELossFit::ELossFit()
26 : fN(0),
27 fNu(0),
28 fChi2(0),
29 fC(0),
30 fDelta(0),
31 fXi(0),
32 fSigma(0),
33 fSigmaN(0),
34 fA(0),
35 fEC(0),
36 fEDelta(0),
37 fEXi(0),
38 fESigma(0),
39 fESigmaN(0),
40 fEA(0),
41 fQuality(0),
42 fDet(0),
43 fRing('\0'),
8449e3e0 44 fBin(0),
45 fMaxWeight(0)
7984e5f7 46{
47 //
48 // Default constructor
49 //
50 //
51}
0bd4b00f 52//____________________________________________________________________
53AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
54 : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ?
6feacd76 55 Int_t(f.GetParameter(AliForwardUtil::ELossFitter::kN)) :
0bd4b00f 56 1),
57 fNu(f.GetNDF()),
58 fChi2(f.GetChisquare()),
59 fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
60 fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
61 fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
62 fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
63 fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
64 fA(0),
65 fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
66 fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
67 fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
68 fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
69 fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
70 fEA(0),
71 fQuality(quality),
72 fDet(0),
73 fRing('\0'),
8449e3e0 74 fBin(0),
75 fMaxWeight(0)
0bd4b00f 76{
7984e5f7 77 //
78 // Construct from a function
79 //
80 // Parameters:
81 // quality Quality flag
82 // f Function
83 //
0bd4b00f 84 if (fN <= 0) return;
85 fA = new Double_t[fN];
86 fEA = new Double_t[fN];
87 for (Int_t i = 0; i < fN-1; i++) {
88 fA[i] = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
89 fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
90 }
91 fA[fN-1] = -9999;
92 fEA[fN-1] = -9999;
93}
94
95//____________________________________________________________________
96AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality,UShort_t n,
97 Double_t chi2, UShort_t nu,
98 Double_t c, Double_t ec,
99 Double_t delta, Double_t edelta,
100 Double_t xi, Double_t exi,
101 Double_t sigma, Double_t esigma,
102 Double_t sigman, Double_t esigman,
fb3430ac 103 const Double_t* a,const Double_t* ea)
0bd4b00f 104 : fN(n),
105 fNu(nu),
106 fChi2(chi2),
107 fC(c),
108 fDelta(delta),
109 fXi(xi),
110 fSigma(sigma),
111 fSigmaN(sigman),
112 fA(0),
113 fEC(ec),
114 fEDelta(edelta),
115 fEXi(exi),
116 fESigma(esigma),
117 fESigmaN(esigman),
118 fEA(0),
119 fQuality(quality),
120 fDet(0),
121 fRing('\0'),
8449e3e0 122 fBin(0),
123 fMaxWeight(0)
0bd4b00f 124{
7984e5f7 125 //
126 // Constructor with full parameter set
127 //
128 // Parameters:
129 // quality Quality flag
130 // n @f$ N@f$ - Number of fitted peaks
131 // chi2 @f$ \chi^2 @f$
132 // nu @f$ \nu @f$ - number degrees of freedom
133 // c @f$ C@f$ - scale constant
134 // ec @f$ \delta C@f$ - error on @f$ C@f$
135 // delta @f$ \Delta@f$ - Most probable value
136 // edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
137 // xi @f$ \xi@f$ - width
138 // exi @f$ \delta\xi@f$ - error on @f$\xi@f$
139 // sigma @f$ \sigma@f$ - Width of Gaussian
140 // esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
141 // sigman @f$ \sigma_n@f$ - Noise width
142 // esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
143 // a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
144 // @f$ i=2,\ldots@f$
145 // ea Array of @f$ N-1@f$ error on the weights @f$ a_i@f$ for
146 // @f$ i=2,\ldots@f$
147 //
0bd4b00f 148 if (fN <= 0) return;
149 fA = new Double_t[fN];
150 fEA = new Double_t[fN];
151 for (Int_t i = 0; i < fN-1; i++) {
152 fA[i] = a[i];
153 fEA[i] = ea[i];
154 }
155 fA[fN-1] = -9999;
156 fEA[fN-1] = -9999;
157}
158//____________________________________________________________________
159AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
160 : TObject(o),
161 fN(o.fN),
162 fNu(o.fNu),
163 fChi2(o.fChi2),
164 fC(o.fC),
165 fDelta(o.fDelta),
166 fXi(o.fXi),
167 fSigma(o.fSigma),
168 fSigmaN(o.fSigmaN),
169 fA(0),
170 fEC(o.fEC),
171 fEDelta(o.fEDelta),
172 fEXi(o.fEXi),
173 fESigma(o.fESigma),
174 fESigmaN(o.fESigmaN),
175 fEA(0),
176 fQuality(o.fQuality),
177 fDet(o.fDet),
178 fRing(o.fRing),
8449e3e0 179 fBin(o.fBin),
180 fMaxWeight(o.fMaxWeight)
0bd4b00f 181{
7984e5f7 182 //
183 // Copy constructor
184 //
185 // Parameters:
186 // o Object to copy from
187 //
0bd4b00f 188 if (fN <= 0) return;
189 fA = new Double_t[fN];
190 fEA = new Double_t[fN];
191 for (Int_t i = 0; i < fN-1; i++) {
192 fA[i] = o.fA[i];
193 fEA[i] = o.fEA[i];
194 }
195 fA[fN-1] = -9999;
196 fEA[fN-1] = -9999;
197}
198
199//____________________________________________________________________
200AliFMDCorrELossFit::ELossFit&
201AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
202{
7984e5f7 203 //
204 // Assignment operator
205 //
206 // Parameters:
207 // o Object to assign from
208 //
209 // Return:
210 // Reference to this object
211 //
d015ecfe 212 if (&o == this) return *this;
8449e3e0 213 fN = o.fN;
214 fNu = o.fNu;
215 fChi2 = o.fChi2;
216 fC = o.fC;
217 fDelta = o.fDelta;
218 fXi = o.fXi;
219 fSigma = o.fSigma;
220 fSigmaN = o.fSigmaN;
221 fEC = o.fEC;
222 fEDelta = o.fEDelta;
223 fEXi = o.fEXi;
224 fESigma = o.fESigma;
225 fESigmaN = o.fESigmaN;
226 fQuality = o.fQuality;
227 fDet = o.fDet;
228 fRing = o.fRing;
229 fBin = o.fBin;
230 fMaxWeight = o.fMaxWeight;
0bd4b00f 231 if (fA) delete [] fA;
232 if (fEA) delete [] fEA;
233 fA = 0;
234 fEA = 0;
235
236 if (fN <= 0) return *this;
237 fA = new Double_t[fN];
238 fEA = new Double_t[fN];
239 for (Int_t i = 0; i < fN; i++) {
240 fA[i] = o.fA[i];
241 fEA[i] = o.fEA[i];
242 }
243
244 return *this;
245}
246
247//____________________________________________________________________
248AliFMDCorrELossFit::ELossFit::~ELossFit()
249{
250 if (fA) delete[] fA;
251 if (fEA) delete[] fEA;
252}
253
254
255//____________________________________________________________________
256Int_t
257AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError,
258 Double_t leastWeight,
259 UShort_t maxN) const
260{
7984e5f7 261 //
262 // Find the maximum weight to use. The maximum weight is the
263 // largest i for which
264 //
265 // - @f$ i \leq \max{N}@f$
266 // - @f$ a_i > \min{a}@f$
267 // - @f$ \delta a_i/a_i > \delta_{max}@f$
268 //
269 // Parameters:
270 // maxRelError @f$ \min{a}@f$
271 // leastWeight @f$ \delta_{max}@f$
272 // maxN @f$ \max{N}@f$
273 //
274 // Return:
275 // The largest index @f$ i@f$ for which the above
276 // conditions hold. Will never return less than 1.
277 //
8449e3e0 278 if (fMaxWeight > 0) return fMaxWeight;
0bd4b00f 279 Int_t n = TMath::Min(maxN, UShort_t(fN-1));
280 Int_t m = 1;
281 // fN is one larger than we have data
282 for (Int_t i = 0; i < n-1; i++, m++) {
283 if (fA[i] < leastWeight) break;
284 if (fEA[i] / fA[i] > maxRelError) break;
285 }
8449e3e0 286 return fMaxWeight = m;
0bd4b00f 287}
288
289//____________________________________________________________________
290Double_t
291AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x,
292 UShort_t maxN) const
293{
7984e5f7 294 //
295 // Evaluate
296 // @f[
297 // f_N(x;\Delta,\xi,\sigma') =
298 // \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
299 // @f]
300 //
301 // (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
302 // that fulfills the requirements
303 //
304 // Parameters:
305 // x Where to evaluate
306 // maxN @f$ \max{N}@f$
307 //
308 // Return:
309 // @f$ f_N(x;\Delta,\xi,\sigma')@f$
310 //
0bd4b00f 311 return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN,
312 TMath::Min(maxN, UShort_t(fN)), fA);
313}
314
315//____________________________________________________________________
316Double_t
317AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x,
318 UShort_t maxN) const
319{
7984e5f7 320 //
321 // Evaluate
322 // @f[
323 // f_W(x;\Delta,\xi,\sigma') =
324 // \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
325 // f_N(x;\Delta,\xi,\sigma')} =
326 // \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
327 // \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
328 // @f]
329 // where @f$ n@f$ fulfills the requirements (see FindMaxWeight).
330 //
331 // If the denominator is zero, then 1 is returned.
332 //
333 // See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
334 // for more information on the evaluated functions.
335 //
336 // Parameters:
337 // x Where to evaluate
338 // maxN @f$ \max{N}@f$
339 //
340 // Return:
341 // @f$ f_W(x;\Delta,\xi,\sigma')@f$.
342 //
0bd4b00f 343 UShort_t n = TMath::Min(maxN, UShort_t(fN-1));
344 Double_t num = 0;
345 Double_t den = 0;
346 for (Int_t i = 1; i <= n; i++) {
347 Double_t a = (i == 1 ? 1 : fA[i-1]);
348 if (fA[i-1] < 0) break;
349 Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
350 num += i * a * f;
351 den += a * f;
352 }
353 if (den <= 0) return 1;
354 return num / den;
355}
356
357
358#define OUTPAR(N,V,E) \
359 std::setprecision(9) << \
360 std::setw(12) << N << ": " << \
361 std::setw(14) << V << " +/- " << \
362 std::setw(14) << E << " (" << \
363 std::setprecision(-1) << \
364 std::setw(5) << 100*(V>0?E/V:1) << "%)\n"
365
366
367//____________________________________________________________________
368Int_t
369AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
370{
7984e5f7 371 //
372 // Compare to another ELossFit object.
373 //
374 // - +1, if this quality is better (larger) than other objects quality
375 // - -1, if this quality is worse (smaller) than other objects quality
376 // - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
377 // - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
378 // - 0 otherwise
379 //
380 // Parameters:
381 // o Other object to compare to
382 //
0bd4b00f 383 const ELossFit* other = static_cast<const ELossFit*>(o);
384 if (this->fQuality < other->fQuality) return -1;
385 if (this->fQuality > other->fQuality) return +1;
386 Double_t chi2nu = (fNu == 0 ? 1000 : fChi2 / fNu);
387 Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
388 if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
389 if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
390 return 0;
391}
392
393//____________________________________________________________________
394void
81775aba 395AliFMDCorrELossFit::ELossFit::Print(Option_t* option) const
0bd4b00f 396{
7984e5f7 397 //
398 // Information to standard output
399 //
400 // Parameters:
401 // option Not used
402 //
81775aba 403 TString o(option);
404 if (o.Contains("S", TString::kIgnoreCase)) {
405 Printf("%15s: q=%2d n=%1d chi2/nu=%6.3f",
406 GetName(), fQuality, fN, (fNu <= 0 ? 999 : fChi2 / fNu));
407 return;
408 }
409
0bd4b00f 410 std::cout << GetName() << ":\n"
411 << " chi^2/nu = " << fChi2 << "/" << fNu << " = "
412 << (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
413 << " Quality: " << fQuality << "\n"
414 << " NParticles: " << fN << " (" << FindMaxWeight() << ")\n"
415 << OUTPAR("Delta", fDelta, fEDelta)
416 << OUTPAR("xi", fXi, fEXi)
417 << OUTPAR("sigma", fSigma, fESigma)
418 << OUTPAR("sigma_n", fSigmaN, fESigmaN);
419 for (Int_t i = 0; i < fN-1; i++)
420 std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
421 std::cout << std::flush;
422}
423
424//____________________________________________________________________
425const Char_t*
426AliFMDCorrELossFit::ELossFit::GetName() const
427{
7984e5f7 428 //
429 // Get the name of this object
430 //
431 //
432 // Return:
433 //
434 //
0bd4b00f 435 return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
436}
437
438//____________________________________________________________________
439void
440AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
441{
7984e5f7 442 //
443 // Browse this object
444 //
445 // Parameters:
446 // b Browser
447 //
8449e3e0 448 Draw(b ? b->GetDrawOption() : "comp values");
0bd4b00f 449 gPad->SetLogy();
450 gPad->Update();
451}
452
453//____________________________________________________________________
454void
455AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
456{
7984e5f7 457 //
458 // Draw this fit
459 //
460 // Parameters:
461 // option Options
462 // - COMP Draw components too
463 //
0bd4b00f 464 TString opt(option);
465 opt.ToUpper();
466 bool comp = false;
8449e3e0 467 bool good = false;
468 bool vals = false;
469 bool legd = false;
0bd4b00f 470 if (opt.Contains("COMP")) {
471 opt.ReplaceAll("COMP","");
472 comp = true;
473 }
8449e3e0 474 if (opt.Contains("GOOD")) {
475 opt.ReplaceAll("GOOD","");
476 good = true;
477 }
478 if (opt.Contains("VALUES")) {
479 opt.ReplaceAll("VALUES","");
480 vals = true;
481 }
482 if (opt.Contains("LEGEND")) {
483 opt.ReplaceAll("LEGEND","");
484 legd = comp;
485 }
0bd4b00f 486 if (!opt.Contains("SAME")) {
487 gPad->Clear();
488 }
489
8449e3e0 490 TLegend* l = 0;
491 if (legd) {
492 l = new TLegend(.3, .5, .59, .94);
493 l->SetBorderSize(0);
494 l->SetFillColor(0);
495 l->SetFillStyle(0);
496 }
0bd4b00f 497 TObjArray cleanup;
8449e3e0 498 Int_t maxW = FindMaxWeight();
499 TF1* tot = AliForwardUtil::MakeNLandauGaus(fC * 1,
0bd4b00f 500 fDelta, fXi,
501 fSigma, fSigmaN,
8449e3e0 502 maxW/*fN*/, fA,
0bd4b00f 503 0.01, 10);
504 tot->SetLineColor(kBlack);
505 tot->SetLineWidth(2);
506 tot->SetLineStyle(1);
507 tot->SetTitle(GetName());
8449e3e0 508 if (l) l->AddEntry(tot, "Total", "l");
0bd4b00f 509 Double_t max = tot->GetMaximum();
510
8449e3e0 511
0bd4b00f 512 if (!opt.Contains("SAME")) {
513 TH1* frame = new TH1F(GetName(),
514 Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
515 100, 0, 10);
516 frame->SetMinimum(max/10000);
517 frame->SetMaximum(max*1.4);
518 frame->SetXTitle("#Delta / #Delta_{mip}");
519 frame->Draw();
520 opt.Append(" SAME");
521 }
522 tot->DrawCopy(opt.Data());
523 cleanup.Add(tot);
524
8449e3e0 525 if (vals) {
526 Double_t x1 = .72;
527 Double_t x2 = .73;
528 Double_t y = .90;
529 Double_t dy = .05;
530 TLatex* ltx1 = new TLatex(x1, y, "");
531 TLatex* ltx2 = new TLatex(x2, y, "");
532 ltx1->SetNDC();
533 ltx1->SetTextAlign(33);
534 ltx1->SetTextFont(132);
535 ltx1->SetTextSize(dy-.01);
536 ltx2->SetNDC();
537 ltx2->SetTextAlign(13);
538 ltx2->SetTextFont(132);
539 ltx2->SetTextSize(dy-.01);
540
541 ltx1->DrawLatex(x1, y, "Quality");
542 ltx2->DrawLatex(x2, y, Form("%d", fQuality));
543 y -= dy;
544
545 ltx1->DrawLatex(x1, y, "#chi^{2}/#nu");
546 ltx2->DrawLatex(x2, y, Form("%7.3f", (fNu > 0 ? fChi2 / fNu : -1)));
547 y -= dy;
548
549 const Char_t* pn[] = { "C", "#Delta", "#xi", "#sigma" };
550 Double_t pv[] = { fC, fDelta, fXi, fSigma };
551 Double_t pe[] = { fEC, fEDelta, fEXi, fESigma };
552 for (Int_t i = 0; i < 4; i++) {
553 ltx1->DrawLatex(x1, y, pn[i]);
554 ltx2->DrawLatex(x2, y, Form("%6.4f#pm%6.4f", pv[i], pe[i]));
555 y -= dy;
556 }
557 for (Int_t i=2; i <= fN; i++) {
558 if (i > maxW) {
559 ltx1->SetTextColor(kRed+3);
560 ltx2->SetTextColor(kRed+3);
561 }
562 ltx1->DrawLatex(x1, y, Form("a_{%d}", i));
563 ltx2->DrawLatex(x2, y, Form("%6.4f#pm%6.4f", fA[i-2], fEA[i-2]));
564 y -= dy;
565 }
566
567 }
568
0bd4b00f 569 if (!comp) {
570 gPad->cd();
571 return;
572 }
573
574 Double_t min = max;
575 opt.Append(" same");
0bd4b00f 576 for (Int_t i=1; i <= fN; i++) {
8449e3e0 577 if (good && i > maxW) break;
578 TF1* f = AliForwardUtil::MakeILandauGaus(fC*(i == 1 ? 1 : fA[i-2]),
0bd4b00f 579 fDelta, fXi,
580 fSigma, fSigmaN,
581 i, 0.01, 10);
582 f->SetLineWidth(2);
583 f->SetLineStyle(i > maxW ? 2 : 1);
584 min = TMath::Min(f->GetMaximum(), min);
585 f->DrawCopy(opt.Data());
8449e3e0 586 if (l) l->AddEntry(f, Form("%d MIP%s", i, (i>1 ? "s" : "")), "l");
587
0bd4b00f 588 cleanup.Add(f);
589 }
590 min /= 100;
591 tot->GetHistogram()->SetMaximum(max);
592 tot->GetHistogram()->SetMinimum(min);
593 tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
8449e3e0 594 if (l) l->Draw();
0bd4b00f 595
596 gPad->cd();
597}
598
599
600//____________________________________________________________________
601#define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
602
16ac05a3 603//____________________________________________________________________
604Double_t
605AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f) const
606{
607 //
608 // Return
609 // Delta * f
610 return f * fDelta;
611}
5bb5d1f6 612//____________________________________________________________________
613Double_t
614AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f,
615 Bool_t includeSigma) const
616{
617 //
618 // Return
619 // Delta - f * (xi + sigma)
620 return fDelta - f * (fXi + (includeSigma ? fSigma : 0));
621}
622
0bd4b00f 623//____________________________________________________________________
624void
625AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
626 Double_t maxRelError,
627 Double_t leastWeight)
628{
7984e5f7 629 //
630 // Calculate the quality
631 //
81775aba 632 Double_t decline = maxChi2nu;
0bd4b00f 633 Int_t qual = 0;
81775aba 634 if (fNu > 0) {
635 Double_t red = fChi2 / fNu;
636 if (red < maxChi2nu) qual += 4;
637 else {
638 Int_t q = (maxChi2nu+decline - red) / decline * 4;
639 if (q > 0) qual += q;
640 }
641 }
0bd4b00f 642 if (CHECKPAR(fDelta, fEDelta, maxRelError)) qual++;
643 if (CHECKPAR(fXi, fEXi, maxRelError)) qual++;
644 if (CHECKPAR(fSigma, fESigma, maxRelError)) qual++;
645 if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
81775aba 646 qual += FindMaxWeight(2*maxRelError, leastWeight, fN);
0bd4b00f 647 fQuality = qual;
648}
649
650//____________________________________________________________________
651AliFMDCorrELossFit::AliFMDCorrELossFit()
652 : TObject(),
653 fRings(),
654 fEtaAxis(0,0,0),
8449e3e0 655 fLowCut(0),
656 fCache(0)
0bd4b00f 657{
7984e5f7 658 //
659 // Default constructor
660 //
0bd4b00f 661 fRings.SetOwner(kTRUE);
662 fEtaAxis.SetTitle("#eta");
663 fEtaAxis.SetName("etaAxis");
664 fRings.SetName("rings");
665}
666
667//____________________________________________________________________
668AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
669 : TObject(o),
670 fRings(o.fRings),
671 fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
8449e3e0 672 fLowCut(0),
673 fCache(0)
0bd4b00f 674{
7984e5f7 675 //
676 // Copy constructor
677 //
678 // Parameters:
679 // o Object to copy from
680 //
0bd4b00f 681 fEtaAxis.SetTitle("#eta");
682 fEtaAxis.SetName("etaAxis");
683}
684
685//____________________________________________________________________
686AliFMDCorrELossFit::~AliFMDCorrELossFit()
687{
7984e5f7 688 //
689 // Destructor
690 //
0bd4b00f 691 fRings.Clear();
692}
693
694//____________________________________________________________________
695AliFMDCorrELossFit&
696AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
697{
7984e5f7 698 //
699 // Assignment operator
700 //
701 // Parameters:
702 // o Object to assign from
703 //
704 // Return:
705 // Reference to this object
706 //
d015ecfe 707 if (&o == this) return *this;
0bd4b00f 708 fRings = o.fRings;
709 fLowCut = o.fLowCut;
8449e3e0 710 fCache = o.fCache;
0bd4b00f 711 SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
712
713 return *this;
714}
8449e3e0 715#define CACHE(BIN,IDX,OFFSET) fCache[IDX*OFFSET+BIN-1]
716
717//____________________________________________________________________
718void
719AliFMDCorrELossFit::CacheBins(UShort_t minQuality) const
720{
721 if (fCache.GetSize() > 0) return;
722
723 Int_t nRings = fRings.GetEntriesFast();
724 Int_t offset = fEtaAxis.GetNbins();
725
726 fCache.Set(nRings*offset);
727 fCache.Reset(-1);
728
729 for (Int_t i = 0; i < nRings; i++) {
730 TObjArray* ringArray = static_cast<TObjArray*>(fRings.At(i));
731
732 // First loop to find where we actually have fits
733 Int_t nFits = 0;
734 Int_t nGood = 0;
735 Int_t minBin = offset+1;
736 Int_t maxBin = -1;
737 Int_t realMinBin = offset+1;
738 Int_t realMaxBin = -1;
739 for (Int_t j = 1; j < ringArray->GetEntriesFast(); j++) {
740 ELossFit* fit = static_cast<ELossFit*>(ringArray->At(j));
741 if (!fit) continue;
742 nFits++;
743
744 // Update our range off possible fits
745 realMinBin = TMath::Min(j, realMinBin);
746 realMaxBin = TMath::Max(j, realMaxBin);
747
748 // Check the quality of the fit
81775aba 749 fit->CalculateQuality(AliFMDCorrELossFit::ELossFit::fgMaxChi2nu,
750 AliFMDCorrELossFit::ELossFit::fgMaxRelError,
751 AliFMDCorrELossFit::ELossFit::fgLeastWeight);
8449e3e0 752 if (minQuality > 0 && fit->fQuality < minQuality) continue;
753 nGood++;
754
755 // Check this bin
756 CACHE(j,i,offset) = j;
757 minBin = TMath::Min(j, minBin);
758 maxBin = TMath::Max(j, maxBin);
759 }
760 AliInfoF("Out of %d bins, %d had fits, of which %d are good (%5.1f%%)",
997ba0f4 761 offset, nFits, nGood, (nFits > 0 ? 100*float(nGood)/nFits : 0));
762
8449e3e0 763 // Now loop and set neighbors
764 realMinBin = TMath::Max(1, realMinBin-1); // Include one more
765 realMaxBin = TMath::Min(offset, realMaxBin+1); // Include one more
766 for (Int_t j = realMinBin; j <= realMaxBin; j++) {
767 if (CACHE(j,i,offset) > 0) continue;
768
769 Int_t nK = TMath::Max(realMaxBin - j, j - realMinBin);
770 Int_t found = -1;
771 for (Int_t k = 1; k <= nK; k++) {
772 Int_t left = j - k;
773 Int_t right = j + k;
774 if (left > realMinBin &&
775 CACHE(left,i,offset) == left) found = left;
776 else if (right < realMaxBin &&
777 CACHE(right,i,offset) == right) found = right;
778 if (found > 0) break;
779 }
780 // Now check that we found something
781 if (found) CACHE(j,i,offset) = CACHE(found,i,offset);
782 else AliWarningF("No fit found for etabin=%d in ring=%d", j, i);
783 }
784 }
785}
786
787
0bd4b00f 788//____________________________________________________________________
789Int_t
790AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
791{
7984e5f7 792 //
793 // Find the eta bin corresponding to the given eta
794 //
795 // Parameters:
796 // eta Eta value
797 //
798 // Return:
799 // Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
800 // eta, or 0 if out of range.
801 //
fb3430ac 802 if (TMath::Abs(fEtaAxis.GetXmin() - fEtaAxis.GetXmax()) < 1e-6
803 || fEtaAxis.GetNbins() == 0) {
0bd4b00f 804 AliWarning("No eta axis defined");
805 return -1;
806 }
807 Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
808 if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
809 return bin;
810}
811
812//____________________________________________________________________
813Bool_t
814AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
815{
7984e5f7 816 //
817 // Set the fit parameters from a function
818 //
819 // Parameters:
820 // d Detector
821 // r Ring
822 // etaBin Eta (bin number, 1->nBins)
823 // f ELoss fit result - note, the object will take ownership
824 //
0bd4b00f 825 TObjArray* ringArray = GetOrMakeRingArray(d, r);
826 if (!ringArray) {
827 AliError(Form("Failed to make ring array for FMD%d%c", d, r));
828 return kFALSE;
829 }
830 if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
831 AliError(Form("bin=%d is out of range [%d,%d]",
832 etaBin, 1, fEtaAxis.GetNbins()));
833 return kFALSE;
834 }
835 // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
836 ringArray->AddAtAndExpand(fit, etaBin);
837 return kTRUE;
838}
839
840//____________________________________________________________________
841Bool_t
842AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
843{
7984e5f7 844 //
845 // Set the fit parameters from a function
846 //
847 // Parameters:
848 // d Detector
849 // r Ring
850 // eta Eta
851 // f ELoss fit result - note, the object will take ownership
852 //
0bd4b00f 853 Int_t bin = FindEtaBin(eta);
854 if (bin <= 0) {
855 AliError(Form("eta=%f is out of range [%f,%f]",
856 eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
857 return kFALSE;
858 }
859
860 return SetFit(d, r, bin, fit);
861}
862//____________________________________________________________________
863Bool_t
864AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r,
865 Double_t eta,
866 Int_t quality,UShort_t n,
867 Double_t chi2, UShort_t nu,
868 Double_t c, Double_t ec,
869 Double_t delta, Double_t edelta,
870 Double_t xi, Double_t exi,
871 Double_t sigma, Double_t esigma,
872 Double_t sigman, Double_t esigman,
873 Double_t* a, Double_t* ea)
874{
7984e5f7 875 //
876 // Set the fit parameters from a function
877 //
878 // Parameters:
879 // d Detector number
880 // r Ring identifier
881 // eta Eta value
882 // quality Quality flag
883 // n @f$ N@f$ - Number of fitted peaks
884 // chi2 @f$ \chi^2 @f$
885 // nu @f$ \nu @f$ - number degrees of freedom
886 // c @f$ C@f$ - scale constant
887 // ec @f$ \delta C@f$ - error on @f$ C@f$
888 // delta @f$ \Delta@f$ - most probable value
889 // edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
890 // xi @f$ \xi@f$ - Landau width
891 // exi @f$ \delta\xi@f$ - error on @f$\xi@f$
892 // sigma @f$ \sigma@f$ - Gaussian width
893 // esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
894 // sigman @f$ \sigma_n@f$ - Noise width
895 // esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
896 // a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
897 // @f$ i=2,\ldots@f$
898 // ea Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for
899 // @f$ i=2,\ldots@f$
900 //
0bd4b00f 901 ELossFit* e = new ELossFit(quality, n,
902 chi2, nu,
903 c, ec,
904 delta, edelta,
905 xi, exi,
906 sigma, esigma,
907 sigman, esigman,
908 a, ea);
909 if (!SetFit(d, r, eta, e)) {
910 delete e;
911 return kFALSE;
912 }
913 return kTRUE;
914}
915//____________________________________________________________________
916Bool_t
917AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta,
918 Int_t quality, const TF1& f)
919{
7984e5f7 920 //
921 // Set the fit parameters from a function
922 //
923 // Parameters:
924 // d Detector
925 // r Ring
926 // eta Eta
927 // quality Quality flag
928 // f Function from fit
929 //
0bd4b00f 930 ELossFit* e = new ELossFit(quality, f);
931 if (!SetFit(d, r, eta, e)) {
932 delete e;
933 return kFALSE;
934 }
935 return kTRUE;
936}
1174780f 937//____________________________________________________________________
938AliFMDCorrELossFit::ELossFit*
939AliFMDCorrELossFit::GetFit(UShort_t d, Char_t r, Int_t etabin) const
940{
941 //
942 // Get the fit corresponding to the specified parameters
943 //
944 // Parameters:
945 // d Detector
946 // r Ring
947 // etabin Eta bin (1 based)
948 //
949 // Return:
950 // Fit parameters or null in case of problems
951 //
952 TObjArray* ringArray = GetRingArray(d, r);
953 if (!ringArray) return 0;
954 if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) return 0;
955 if (etabin > ringArray->GetEntriesFast()) return 0;
956 else if (etabin >= ringArray->GetEntriesFast()) etabin--;
957 else if (!ringArray->At(etabin)) etabin++;
958 return static_cast<ELossFit*>(ringArray->At(etabin));
959}
960//____________________________________________________________________
961AliFMDCorrELossFit::ELossFit*
962AliFMDCorrELossFit::GetFit(UShort_t d, Char_t r, Double_t eta) const
963{
964 //
965 // Find the fit corresponding to the specified parameters
966 //
967 // Parameters:
968 // d Detector
969 // r Ring
970 // eta Eta value
971 //
972 // Return:
973 // Fit parameters or null in case of problems
974 //
975 Int_t etabin = FindEtaBin(eta);
976 return GetFit(d, r, etabin);
977}
978
0bd4b00f 979//____________________________________________________________________
980AliFMDCorrELossFit::ELossFit*
8449e3e0 981AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Int_t etabin,
982 UShort_t minQ) const
0bd4b00f 983{
7984e5f7 984 //
985 // Find the fit corresponding to the specified parameters
986 //
987 // Parameters:
988 // d Detector
989 // r Ring
990 // etabin Eta bin (1 based)
991 //
992 // Return:
993 // Fit parameters or null in case of problems
994 //
0bd4b00f 995 if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
5bb5d1f6 996 // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
997 // etabin, 1, fEtaAxis.GetNbins(), d, r));
0bd4b00f 998 return 0;
999 }
8449e3e0 1000
1001 TObjArray* ringArray = GetRingArray(d, r);
1002 if (!ringArray) {
1003 AliError(Form("Failed to make ring array for FMD%d%c", d, r));
0bd4b00f 1004 return 0;
1005 }
8449e3e0 1006 if (fCache.GetSize() <= 0) CacheBins(minQ);
1007 Int_t idx = (d == 1 ? 0 :
1008 (d - 2) * 2 + 1 + (r=='I' || r=='i' ? 0 : 1));
1009 Int_t bin = CACHE(etabin, idx, fEtaAxis.GetNbins());
1010
1011 if (bin < 0 || bin > ringArray->GetEntriesFast()) return 0;
1012
1013 return static_cast<ELossFit*>(ringArray->At(bin));
0bd4b00f 1014}
1015//____________________________________________________________________
1016AliFMDCorrELossFit::ELossFit*
8449e3e0 1017AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Double_t eta,
1018 UShort_t minQ) const
0bd4b00f 1019{
7984e5f7 1020 //
1021 // Find the fit corresponding to the specified parameters
1022 //
1023 // Parameters:
1024 // d Detector
1025 // r Ring
1026 // eta Eta value
1027 //
1028 // Return:
1029 // Fit parameters or null in case of problems
1030 //
0bd4b00f 1031 Int_t etabin = FindEtaBin(eta);
8449e3e0 1032 return FindFit(d, r, etabin, minQ);
0bd4b00f 1033}
1034//____________________________________________________________________
1035TObjArray*
1036AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
1037{
7984e5f7 1038 //
1039 // Get the ring array corresponding to the specified ring
1040 //
1041 // Parameters:
1042 // d Detector
1043 // r Ring
1044 //
1045 // Return:
1046 // Pointer to ring array, or null in case of problems
1047 //
0bd4b00f 1048 Int_t idx = -1;
1049 switch (d) {
1050 case 1: idx = 0; break;
1051 case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
1052 case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
1053 }
1054 if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
1055 return static_cast<TObjArray*>(fRings.At(idx));
1056}
1057//____________________________________________________________________
1058TObjArray*
1059AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
1060{
7984e5f7 1061 //
1062 // Get the ring array corresponding to the specified ring
1063 //
1064 // Parameters:
1065 // d Detector
1066 // r Ring
1067 //
1068 // Return:
1069 // Pointer to ring array, or newly created container
1070 //
0bd4b00f 1071 Int_t idx = -1;
1072 switch (d) {
1073 case 1: idx = 0; break;
1074 case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
1075 case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
1076 }
1077 if (idx < 0) return 0;
1078 if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
1079 TObjArray* a = new TObjArray(0);
1080 // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
1081 a->SetName(Form("FMD%d%c", d, r));
1082 a->SetOwner();
1083 fRings.AddAtAndExpand(a, idx);
1084 }
1085 return static_cast<TObjArray*>(fRings.At(idx));
1086}
1087
16ac05a3 1088//____________________________________________________________________
1089Double_t
1090AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Int_t etabin,
1091 Double_t f) const
1092{
8449e3e0 1093 ELossFit* fit = FindFit(d, r, etabin, 20);
16ac05a3 1094 if (!fit) return -1024;
1095 return fit->GetLowerBound(f);
1096}
1097//____________________________________________________________________
1098Double_t
1099AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Double_t eta,
1100 Double_t f) const
1101{
1102 Int_t bin = FindEtaBin(eta);
1103 if (bin <= 0) return -1024;
1104 return GetLowerBound(d, r, Int_t(bin), f);
1105}
5bb5d1f6 1106//____________________________________________________________________
1107Double_t
1108AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Int_t etabin,
1109 Double_t f, Bool_t showErrors,
1110 Bool_t includeSigma) const
1111{
8449e3e0 1112 ELossFit* fit = FindFit(d, r, etabin, 20);
5bb5d1f6 1113 if (!fit) {
1114 if (showErrors) {
1115 AliWarning(Form("No fit for FMD%d%c @ etabin=%d", d, r, etabin));
1116 }
1117 return -1024;
1118 }
1119 return fit->GetLowerBound(f, includeSigma);
1120}
1121
1122//____________________________________________________________________
1123Double_t
1124AliFMDCorrELossFit::GetLowerBound(UShort_t d, Char_t r, Double_t eta,
1125 Double_t f, Bool_t showErrors,
1126 Bool_t includeSigma) const
1127{
1128 Int_t bin = FindEtaBin(eta);
1129 if (bin <= 0) {
1130 if (showErrors)
1131 AliError(Form("eta=%f out of bounds for FMD%d%c", eta, d, r));
1132 return -1024;
1133 }
1134 return GetLowerBound(d, r, bin, f, showErrors, includeSigma);
1135}
1136
1137//____________________________________________________________________
0bd4b00f 1138namespace {
1139 TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
1140 Int_t color)
1141 {
1142 TH1D* h = new TH1D(Form("%s_%s", name, title),
1143 Form("%s %s", name, title),
1144 axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
1145 h->SetDirectory(0);
1146 h->SetMarkerStyle(20);
1147 h->SetMarkerColor(color);
1148 h->SetMarkerSize(0.5);
1149 h->SetFillColor(color);
1150 h->SetFillStyle(3001);
1151 h->SetLineColor(color);
1152 return h;
1153 }
1154}
1155
5bb5d1f6 1156
1157
7dacdf84 1158#define IDX2RING(I) (i == 0 || i == 1 || i == 3 ? 'I' : 'O')
1159#define IDX2DET(I) (i == 0 ? 1 : (i == 1 || i == 2 ? 2 : 3))
0bd4b00f 1160//____________________________________________________________________
f7cfc454 1161TList*
8449e3e0 1162AliFMDCorrELossFit::GetStacks(Bool_t err,
1163 Bool_t rel,
1164 Bool_t /*good*/,
1165 UShort_t maxN) const
0bd4b00f 1166{
f7cfc454 1167 // Get a list of THStacks
1168 Int_t nRings = fRings.GetEntriesFast();
8449e3e0 1169 // Int_t nPad = 6+maxN-1; // 7 regular params, and maxN-1 weights
7dacdf84 1170
1171 enum {
1172 kChi2nu = 0,
1173 kC = 1,
1174 kDelta = 2,
1175 kXi = 3,
1176 kSigma = 4,
1177 kN = 5
1178 };
0bd4b00f 1179
7dacdf84 1180 THStack* sChi2nu;
1181 THStack* sC;
1182 THStack* sDelta;
1183 THStack* sXi;
1184 THStack* sSigma;
1185 // THStack* sigman;
0bd4b00f 1186 THStack* n;
f7cfc454 1187 TList* stacks = new TList;
1188 stacks->AddAt(sChi2nu= new THStack("chi2", "#chi^{2}/#nu"), kChi2nu);
1189 stacks->AddAt(sC = new THStack("c", "C"), kC);
1190 stacks->AddAt(sDelta = new THStack("delta", "#Delta_{mp}"), kDelta);
1191 stacks->AddAt(sXi = new THStack("xi", "#xi"), kXi);
1192 stacks->AddAt(sSigma = new THStack("sigma", "#sigma"), kSigma);
1193 //stacks->AddAt(sigman= new THStack("sigman", "#sigma_{n}"), 5);
8449e3e0 1194 stacks->AddAt(n = new THStack("n", "N"), kN);
1195 if (rel) {
1196 sChi2nu->SetName("qual");
1197 sChi2nu->SetTitle("Quality");
1198 n->SetName("good");
1199 n->SetTitle("Bin map");
1200 }
0bd4b00f 1201 for (Int_t i = 1; i <= maxN; i++) {
f7cfc454 1202 stacks->AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), kN+i);
0bd4b00f 1203 }
f7cfc454 1204
8449e3e0 1205 // TArrayD min(nPad);
1206 // TArrayD max(nPad);
1207 // min.Reset(100000);
1208 // max.Reset(-100000);
7dacdf84 1209
0bd4b00f 1210 for (Int_t i = 0; i < nRings; i++) {
1211 if (!fRings.At(i)) continue;
1212 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1213 Int_t nFits = a->GetEntriesFast();
7dacdf84 1214 Int_t color = AliForwardUtil::RingColor(IDX2DET(i), IDX2RING(i));
0bd4b00f 1215
1216 TH1D* hChi = MakeHist(fEtaAxis,a->GetName(), "chi2", color);
1217 TH1D* hC = MakeHist(fEtaAxis,a->GetName(), "c", color);
1218 TH1D* hDelta = MakeHist(fEtaAxis,a->GetName(), "delta", color);
1219 TH1D* hXi = MakeHist(fEtaAxis,a->GetName(), "xi", color);
1220 TH1D* hSigma = MakeHist(fEtaAxis,a->GetName(), "sigma", color);
7dacdf84 1221 // TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
0bd4b00f 1222 TH1D* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
1223 TH1D* hA[maxN];
1224 const char* ho = (rel || !err ? "hist" : "e");
7dacdf84 1225 sChi2nu->Add(hChi, "hist"); // 0
1226 sC ->Add(hC, ho); // 1
1227 sDelta ->Add(hDelta, ho); // 2
1228 sXi ->Add(hXi, ho); // 3
1229 sSigma ->Add(hSigma, ho); // 4
1230 // sigman->Add(hSigmaN,ho); // 5
1231 n ->Add(hN, "hist"); // 5
0bd4b00f 1232 hChi->SetFillColor(color);
1233 hChi->SetFillStyle(3001);
1234 hN->SetFillColor(color);
1235 hN->SetFillStyle(3001);
1236
1237 for (Int_t k = 1; k <= maxN; k++) {
1238 hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
f7cfc454 1239 static_cast<THStack*>(stacks->At(kN+k))->Add(hA[k-1], ho);
0bd4b00f 1240 }
1241
1242 for (Int_t j = 0; j < nFits; j++) {
1243 ELossFit* f = static_cast<ELossFit*>(a->At(j));
1244 if (!f) continue;
1245
8449e3e0 1246 Int_t b = f->fBin;
1247 Double_t vChi2nu = (rel ? f->fQuality
1248 : (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
7dacdf84 1249 Double_t vC = (rel ? (f->fC >0 ?f->fEC /f->fC :0)
1250 : f->fC);
1251 Double_t vDelta = (rel ? (f->fDelta >0 ?f->fEDelta /f->fDelta :0)
1252 : f->fDelta);
1253 Double_t vXi = (rel ? (f->fXi >0 ?f->fEXi /f->fXi :0)
1254 : f->fXi);
1255 Double_t vSigma = (rel ? (f->fSigma >0 ?f->fESigma /f->fSigma :0)
1256 : f->fSigma);
8449e3e0 1257 Int_t nW = (rel ? CACHE(j,i,fEtaAxis.GetNbins()) :
1258 f->FindMaxWeight());
7dacdf84 1259 // Double_t sigman = (rel ? (f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0)
1260 // : f->SigmaN);
1261 hChi ->SetBinContent(b, vChi2nu);
0bd4b00f 1262 hN ->SetBinContent(b, nW);
7dacdf84 1263 hC ->SetBinContent(b, vC);
1264 hDelta ->SetBinContent(b, vDelta);
1265 hXi ->SetBinContent(b, vXi);
1266 hSigma ->SetBinContent(b, vSigma);
8449e3e0 1267 // if (vChi2nu > 1e-12) {
1268 // min[kChi2nu] = TMath::Min(min[kChi2nu], vChi2nu);
1269 // max[kChi2nu] = TMath::Max(max[kChi2nu], vChi2nu);
1270 // }
1271 // if (vC > 1e-12) {
1272 // min[kC] = TMath::Min(min[kC], vC);
1273 // max[kC] = TMath::Max(max[kC], vC);
1274 // }
1275 // if (vDelta > 1e-12) {
1276 // min[kDelta] = TMath::Min(min[kDelta], vDelta);
1277 // max[kDelta] = TMath::Max(max[kDelta], vDelta);
1278 // }
1279 // if (vXi > 1e-12) {
1280 // min[kXi] = TMath::Min(min[kXi], vXi);
1281 // max[kXi] = TMath::Max(max[kXi], vXi);
1282 // }
1283 // if (vSigma > 1e-12) {
1284 // min[kSigma] = TMath::Min(min[kSigma], vSigma);
1285 // max[kSigma] = TMath::Max(max[kSigma], vSigma);
1286 // }
1287 // if (nW > 1e-12) {
1288 // min[kN] = TMath::Min(min[kN], Double_t(nW));
1289 // max[kN] = TMath::Max(max[kN], Double_t(nW));
1290 // }
7dacdf84 1291 // hSigmaN->SetBinContent(b, sigman);
1292 if (!rel) {
0bd4b00f 1293 hC ->SetBinError(b, f->fEC);
1294 hDelta ->SetBinError(b, f->fEDelta);
1295 hXi ->SetBinError(b, f->fEXi);
1296 hSigma ->SetBinError(b, f->fESigma);
7dacdf84 1297 // hSigmaN->SetBinError(b, f->fESigmaN);
0bd4b00f 1298 }
1299 for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
7dacdf84 1300 Double_t vA = (rel ? (f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0)
1301 : f->fA[k]);
1302 hA[k]->SetBinContent(b, vA);
1303 if (!rel) hA[k]->SetBinError(b, f->fEA[k]);
8449e3e0 1304 // if (vA > 1e-12) {
1305 // min[kN+1+k] = TMath::Min(min[kN+1+k], vA);
1306 // max[kN+1+k] = TMath::Max(max[kN+1+k], vA);
1307 // }
0bd4b00f 1308 }
1309 }
1310 }
f7cfc454 1311 return stacks;
1312}
1313
1314//____________________________________________________________________
1315void
1316AliFMDCorrELossFit::Draw(Option_t* option)
1317{
1318 //
1319 // Draw this object
1320 //
1321 // Parameters:
1322 // option Options. Possible values are
1323 // - err Plot error bars
1324 //
1325 TString opt(Form("nostack %s", option));
1326 opt.ToLower();
1327 Bool_t rel = (opt.Contains("relative"));
1328 Bool_t err = (opt.Contains("error"));
8449e3e0 1329 Bool_t clr = (opt.Contains("clear"));
1330 Bool_t gdd = (opt.Contains("good"));
f7cfc454 1331 if (rel) opt.ReplaceAll("relative","");
1332 if (err) opt.ReplaceAll("error","");
8449e3e0 1333 if (clr) opt.ReplaceAll("clear", "");
1334 if (gdd) opt.ReplaceAll("good", "");
f7cfc454 1335
1336 UShort_t maxN = 0;
1337 Int_t nRings = fRings.GetEntriesFast();
1338 for (Int_t i = 0; i < nRings; i++) {
1339 if (!fRings.At(i)) continue;
1340 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1341 Int_t nFits = a->GetEntriesFast();
1342
1343 for (Int_t j = 0; j < nFits; j++) {
1344 ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1345 if (!fit) continue;
1346 maxN = TMath::Max(maxN, UShort_t(fit->fN));
1347 }
1348 }
1349 // AliInfo(Form("Maximum N is %d", maxN));
1350 Int_t nPad = 6+maxN-1; // 7 regular params, and maxN-1 weights
1351 TVirtualPad* pad = gPad;
8449e3e0 1352 if (clr) {
1353 pad->Clear();
1354 pad->SetTopMargin(0.02);
1355 pad->SetRightMargin(0.02);
1356 pad->SetBottomMargin(0.15);
1357 pad->SetLeftMargin(0.10);
1358 }
f7cfc454 1359 pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
1360
8449e3e0 1361 TList* stacks = GetStacks(err, rel, gdd, maxN);
f7cfc454 1362
0bd4b00f 1363 Int_t nPad2 = (nPad+1) / 2;
1364 for (Int_t i = 0; i < nPad; i++) {
1365 Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
1366 TVirtualPad* p = pad->cd(iPad);
1367 p->SetLeftMargin(.15);
1368 p->SetFillColor(0);
1369 p->SetFillStyle(0);
1370 p->SetGridx();
1371 p->SetGridy();
1372 if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
1373
7dacdf84 1374
f7cfc454 1375 THStack* stack = static_cast<THStack*>(stacks->At(i));
1376
1377 // Double_t powMax = TMath::Log10(max[i]);
1378 // Double_t powMin = min[i] <= 0 ? powMax : TMath::Log10(min[i]);
1379 // if (powMax-powMin > 2. && min[i] != 0) p->SetLogy();
1380
7dacdf84 1381 // stack->SetMinimum(min[i]);
1382 // stack->SetMaximum(max[i]);
0bd4b00f 1383 stack->Draw(opt.Data());
1384
1385 TString tit(stack->GetTitle());
8449e3e0 1386 if (rel && i != 0 && i != 5)
0bd4b00f 1387 tit = Form("#delta %s/%s", tit.Data(), tit.Data());
1388 TH1* hist = stack->GetHistogram();
1389 TAxis* yaxis = hist->GetYaxis();
1390 yaxis->SetTitle(tit.Data());
1391 yaxis->SetTitleSize(0.15);
1392 yaxis->SetLabelSize(0.08);
1393 yaxis->SetTitleOffset(0.35);
7dacdf84 1394 yaxis->SetTitleFont(132);
1395 yaxis->SetLabelFont(132);
0bd4b00f 1396 yaxis->SetNdivisions(5);
1397
7dacdf84 1398
0bd4b00f 1399 TAxis* xaxis = stack->GetHistogram()->GetXaxis();
1400 xaxis->SetTitle("#eta");
1401 xaxis->SetTitleSize(0.15);
1402 xaxis->SetLabelSize(0.08);
1403 xaxis->SetTitleOffset(0.35);
7dacdf84 1404 xaxis->SetTitleFont(132);
1405 xaxis->SetLabelFont(132);
0bd4b00f 1406 xaxis->SetNdivisions(10);
1407
7dacdf84 1408 stack->Draw(opt.Data());
0bd4b00f 1409 }
1410 pad->cd();
1411}
1412
1413//____________________________________________________________________
1414void
1415AliFMDCorrELossFit::Print(Option_t* option) const
1416{
fb3430ac 1417 //
1418 // Print this object.
1419 //
1420 // Parameters:
1421 // option Options
1422 // - R Print recursive
1423 //
1424 //
0bd4b00f 1425 TString opt(option);
1426 opt.ToUpper();
8449e3e0 1427 Int_t nRings = fRings.GetEntriesFast();
1428 bool recurse = opt.Contains("R");
1429 bool cache = opt.Contains("C") && fCache.GetSize() > 0;
1430 Int_t nBins = fEtaAxis.GetNbins();
0bd4b00f 1431
1432 std::cout << "Low cut in fit range: " << fLowCut << "\n"
8449e3e0 1433 << "Eta axis: " << nBins
0bd4b00f 1434 << " bins, range [" << fEtaAxis.GetXmin() << ","
1435 << fEtaAxis.GetXmax() << "]" << std::endl;
1436
1437 for (Int_t i = 0; i < nRings; i++) {
1438 if (!fRings.At(i)) continue;
1439 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1440 Int_t nFits = a->GetEntriesFast();
1441
1442 std::cout << a->GetName() << " [" << nFits << " entries]"
1443 << (recurse ? ":\n" : "\t");
1444 Int_t min = fEtaAxis.GetNbins()+1;
1445 Int_t max = 0;
1446 for (Int_t j = 0; j < nFits; j++) {
1447 if (!a->At(j)) continue;
1448
1449 min = TMath::Min(j, min);
1450 max = TMath::Max(j, max);
1451
1452 if (recurse) {
1453 std::cout << "Bin # " << j << "\t";
1454 ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1455 fit->Print(option);
1456 }
1457 }
1458 if (!recurse)
1459 std::cout << " bin range: " << std::setw(3) << min
1460 << "-" << std::setw(3) << max << " " << std::setw(3)
1461 << (max-min+1) << " bins" << std::endl;
1462 }
8449e3e0 1463
1464 if (!cache) return;
1465
1466 std::cout << " eta bin | Fit bin \n"
1467 << " # range | FMD1i FMD2i FMD2o FMD3i FMD3o"
1468 // << "----+-----+++------+-----------------------------------"
1469 << std::endl;
1470 size_t oldPrec = std::cout.precision();
1471 std::cout.precision(3);
1472 for (Int_t i = 1; i <= nBins; i++) {
1473 std::cout << std::setw(4) << i << " "
1474 << std::setw(5) << std::showpos << fEtaAxis.GetBinLowEdge(i)
1475 << " - " << std::setw(5) << fEtaAxis.GetBinUpEdge(i)
1476 << std::noshowpos << " | ";
1477 for (Int_t j = 0; j < 5; j++) {
1478 Int_t bin = CACHE(i,j,nBins);
1479 if (bin <= 0) std::cout << " ";
1480 else std::cout << std::setw(5) << bin
1481 << (bin == i ? ' ' : '*') << ' ';
1482 }
1483 std::cout << std::endl;
1484 }
1485 std::cout.precision(oldPrec);
0bd4b00f 1486}
1487
1488//____________________________________________________________________
1489void
1490AliFMDCorrELossFit::Browse(TBrowser* b)
1491{
7984e5f7 1492 //
1493 // Browse this object
1494 //
1495 // Parameters:
1496 // b
1497 //
0bd4b00f 1498 b->Add(&fRings);
1499 b->Add(&fEtaAxis);
1500}
1501
1502
1503
1504//____________________________________________________________________
1505//
1506// EOF
1507//