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