]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
small fix to the merging efficiency on or off
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDCorrELossFit.cxx
CommitLineData
7984e5f7 1// Object holding the Energy loss fit 'correction'
2//
3// These are generated from Monte-Carlo or real ESDs.
0bd4b00f 4#include "AliFMDCorrELossFit.h"
5#include "AliForwardUtil.h"
6#include <TF1.h>
7#include <TBrowser.h>
8#include <TVirtualPad.h>
9#include <THStack.h>
10#include <TH1D.h>
11#include <AliLog.h>
12#include <TList.h>
13#include <iostream>
14#include <iomanip>
15
16Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .2;
17Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
18Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu = 20;
19
20//____________________________________________________________________
21AliFMDCorrELossFit::ELossFit::ELossFit()
22 : fN(0),
23 fNu(0),
24 fChi2(0),
25 fC(0),
26 fDelta(0),
27 fXi(0),
28 fSigma(0),
29 fSigmaN(0),
30 fA(0),
31 fEC(0),
32 fEDelta(0),
33 fEXi(0),
34 fESigma(0),
35 fESigmaN(0),
36 fEA(0),
37 fQuality(0),
38 fDet(0),
39 fRing('\0'),
40 fBin(0)
7984e5f7 41{
42 //
43 // Default constructor
44 //
45 //
46}
0bd4b00f 47//____________________________________________________________________
48AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
49 : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ?
50 f.GetParameter(AliForwardUtil::ELossFitter::kN) :
51 1),
52 fNu(f.GetNDF()),
53 fChi2(f.GetChisquare()),
54 fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
55 fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
56 fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
57 fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
58 fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
59 fA(0),
60 fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
61 fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
62 fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
63 fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
64 fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
65 fEA(0),
66 fQuality(quality),
67 fDet(0),
68 fRing('\0'),
69 fBin(0)
70{
7984e5f7 71 //
72 // Construct from a function
73 //
74 // Parameters:
75 // quality Quality flag
76 // f Function
77 //
0bd4b00f 78 if (fN <= 0) return;
79 fA = new Double_t[fN];
80 fEA = new Double_t[fN];
81 for (Int_t i = 0; i < fN-1; i++) {
82 fA[i] = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
83 fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
84 }
85 fA[fN-1] = -9999;
86 fEA[fN-1] = -9999;
87}
88
89//____________________________________________________________________
90AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality,UShort_t n,
91 Double_t chi2, UShort_t nu,
92 Double_t c, Double_t ec,
93 Double_t delta, Double_t edelta,
94 Double_t xi, Double_t exi,
95 Double_t sigma, Double_t esigma,
96 Double_t sigman, Double_t esigman,
97 Double_t* a, Double_t* ea)
98 : fN(n),
99 fNu(nu),
100 fChi2(chi2),
101 fC(c),
102 fDelta(delta),
103 fXi(xi),
104 fSigma(sigma),
105 fSigmaN(sigman),
106 fA(0),
107 fEC(ec),
108 fEDelta(edelta),
109 fEXi(exi),
110 fESigma(esigma),
111 fESigmaN(esigman),
112 fEA(0),
113 fQuality(quality),
114 fDet(0),
115 fRing('\0'),
116 fBin(0)
117{
7984e5f7 118 //
119 // Constructor with full parameter set
120 //
121 // Parameters:
122 // quality Quality flag
123 // n @f$ N@f$ - Number of fitted peaks
124 // chi2 @f$ \chi^2 @f$
125 // nu @f$ \nu @f$ - number degrees of freedom
126 // c @f$ C@f$ - scale constant
127 // ec @f$ \delta C@f$ - error on @f$ C@f$
128 // delta @f$ \Delta@f$ - Most probable value
129 // edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
130 // xi @f$ \xi@f$ - width
131 // exi @f$ \delta\xi@f$ - error on @f$\xi@f$
132 // sigma @f$ \sigma@f$ - Width of Gaussian
133 // esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
134 // sigman @f$ \sigma_n@f$ - Noise width
135 // esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
136 // a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
137 // @f$ i=2,\ldots@f$
138 // ea Array of @f$ N-1@f$ error on the weights @f$ a_i@f$ for
139 // @f$ i=2,\ldots@f$
140 //
0bd4b00f 141 if (fN <= 0) return;
142 fA = new Double_t[fN];
143 fEA = new Double_t[fN];
144 for (Int_t i = 0; i < fN-1; i++) {
145 fA[i] = a[i];
146 fEA[i] = ea[i];
147 }
148 fA[fN-1] = -9999;
149 fEA[fN-1] = -9999;
150}
151//____________________________________________________________________
152AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
153 : TObject(o),
154 fN(o.fN),
155 fNu(o.fNu),
156 fChi2(o.fChi2),
157 fC(o.fC),
158 fDelta(o.fDelta),
159 fXi(o.fXi),
160 fSigma(o.fSigma),
161 fSigmaN(o.fSigmaN),
162 fA(0),
163 fEC(o.fEC),
164 fEDelta(o.fEDelta),
165 fEXi(o.fEXi),
166 fESigma(o.fESigma),
167 fESigmaN(o.fESigmaN),
168 fEA(0),
169 fQuality(o.fQuality),
170 fDet(o.fDet),
171 fRing(o.fRing),
172 fBin(o.fBin)
173{
7984e5f7 174 //
175 // Copy constructor
176 //
177 // Parameters:
178 // o Object to copy from
179 //
0bd4b00f 180 if (fN <= 0) return;
181 fA = new Double_t[fN];
182 fEA = new Double_t[fN];
183 for (Int_t i = 0; i < fN-1; i++) {
184 fA[i] = o.fA[i];
185 fEA[i] = o.fEA[i];
186 }
187 fA[fN-1] = -9999;
188 fEA[fN-1] = -9999;
189}
190
191//____________________________________________________________________
192AliFMDCorrELossFit::ELossFit&
193AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
194{
7984e5f7 195 //
196 // Assignment operator
197 //
198 // Parameters:
199 // o Object to assign from
200 //
201 // Return:
202 // Reference to this object
203 //
0bd4b00f 204 fN = o.fN;
205 fNu = o.fNu;
206 fChi2 = o.fChi2;
207 fC = o.fC;
208 fDelta = o.fDelta;
209 fXi = o.fXi;
210 fSigma = o.fSigma;
211 fSigmaN = o.fSigmaN;
212 fEC = o.fEC;
213 fEDelta = o.fEDelta;
214 fEXi = o.fEXi;
215 fESigma = o.fESigma;
216 fESigmaN = o.fESigmaN;
217 fQuality = o.fQuality;
218 fDet = o.fDet;
219 fRing = o.fRing;
220 fBin = o.fBin;
221 if (fA) delete [] fA;
222 if (fEA) delete [] fEA;
223 fA = 0;
224 fEA = 0;
225
226 if (fN <= 0) return *this;
227 fA = new Double_t[fN];
228 fEA = new Double_t[fN];
229 for (Int_t i = 0; i < fN; i++) {
230 fA[i] = o.fA[i];
231 fEA[i] = o.fEA[i];
232 }
233
234 return *this;
235}
236
237//____________________________________________________________________
238AliFMDCorrELossFit::ELossFit::~ELossFit()
239{
240 if (fA) delete[] fA;
241 if (fEA) delete[] fEA;
242}
243
244
245//____________________________________________________________________
246Int_t
247AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError,
248 Double_t leastWeight,
249 UShort_t maxN) const
250{
7984e5f7 251 //
252 // Find the maximum weight to use. The maximum weight is the
253 // largest i for which
254 //
255 // - @f$ i \leq \max{N}@f$
256 // - @f$ a_i > \min{a}@f$
257 // - @f$ \delta a_i/a_i > \delta_{max}@f$
258 //
259 // Parameters:
260 // maxRelError @f$ \min{a}@f$
261 // leastWeight @f$ \delta_{max}@f$
262 // maxN @f$ \max{N}@f$
263 //
264 // Return:
265 // The largest index @f$ i@f$ for which the above
266 // conditions hold. Will never return less than 1.
267 //
0bd4b00f 268 Int_t n = TMath::Min(maxN, UShort_t(fN-1));
269 Int_t m = 1;
270 // fN is one larger than we have data
271 for (Int_t i = 0; i < n-1; i++, m++) {
272 if (fA[i] < leastWeight) break;
273 if (fEA[i] / fA[i] > maxRelError) break;
274 }
275 return m;
276}
277
278//____________________________________________________________________
279Double_t
280AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x,
281 UShort_t maxN) const
282{
7984e5f7 283 //
284 // Evaluate
285 // @f[
286 // f_N(x;\Delta,\xi,\sigma') =
287 // \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
288 // @f]
289 //
290 // (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
291 // that fulfills the requirements
292 //
293 // Parameters:
294 // x Where to evaluate
295 // maxN @f$ \max{N}@f$
296 //
297 // Return:
298 // @f$ f_N(x;\Delta,\xi,\sigma')@f$
299 //
0bd4b00f 300 return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN,
301 TMath::Min(maxN, UShort_t(fN)), fA);
302}
303
304//____________________________________________________________________
305Double_t
306AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x,
307 UShort_t maxN) const
308{
7984e5f7 309 //
310 // Evaluate
311 // @f[
312 // f_W(x;\Delta,\xi,\sigma') =
313 // \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
314 // f_N(x;\Delta,\xi,\sigma')} =
315 // \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
316 // \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
317 // @f]
318 // where @f$ n@f$ fulfills the requirements (see FindMaxWeight).
319 //
320 // If the denominator is zero, then 1 is returned.
321 //
322 // See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
323 // for more information on the evaluated functions.
324 //
325 // Parameters:
326 // x Where to evaluate
327 // maxN @f$ \max{N}@f$
328 //
329 // Return:
330 // @f$ f_W(x;\Delta,\xi,\sigma')@f$.
331 //
0bd4b00f 332 UShort_t n = TMath::Min(maxN, UShort_t(fN-1));
333 Double_t num = 0;
334 Double_t den = 0;
335 for (Int_t i = 1; i <= n; i++) {
336 Double_t a = (i == 1 ? 1 : fA[i-1]);
337 if (fA[i-1] < 0) break;
338 Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
339 num += i * a * f;
340 den += a * f;
341 }
342 if (den <= 0) return 1;
343 return num / den;
344}
345
346
347#define OUTPAR(N,V,E) \
348 std::setprecision(9) << \
349 std::setw(12) << N << ": " << \
350 std::setw(14) << V << " +/- " << \
351 std::setw(14) << E << " (" << \
352 std::setprecision(-1) << \
353 std::setw(5) << 100*(V>0?E/V:1) << "%)\n"
354
355
356//____________________________________________________________________
357Int_t
358AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
359{
7984e5f7 360 //
361 // Compare to another ELossFit object.
362 //
363 // - +1, if this quality is better (larger) than other objects quality
364 // - -1, if this quality is worse (smaller) than other objects quality
365 // - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
366 // - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
367 // - 0 otherwise
368 //
369 // Parameters:
370 // o Other object to compare to
371 //
0bd4b00f 372 const ELossFit* other = static_cast<const ELossFit*>(o);
373 if (this->fQuality < other->fQuality) return -1;
374 if (this->fQuality > other->fQuality) return +1;
375 Double_t chi2nu = (fNu == 0 ? 1000 : fChi2 / fNu);
376 Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
377 if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
378 if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
379 return 0;
380}
381
382//____________________________________________________________________
383void
384AliFMDCorrELossFit::ELossFit::Print(Option_t*) const
385{
7984e5f7 386 //
387 // Information to standard output
388 //
389 // Parameters:
390 // option Not used
391 //
0bd4b00f 392 std::cout << GetName() << ":\n"
393 << " chi^2/nu = " << fChi2 << "/" << fNu << " = "
394 << (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
395 << " Quality: " << fQuality << "\n"
396 << " NParticles: " << fN << " (" << FindMaxWeight() << ")\n"
397 << OUTPAR("Delta", fDelta, fEDelta)
398 << OUTPAR("xi", fXi, fEXi)
399 << OUTPAR("sigma", fSigma, fESigma)
400 << OUTPAR("sigma_n", fSigmaN, fESigmaN);
401 for (Int_t i = 0; i < fN-1; i++)
402 std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
403 std::cout << std::flush;
404}
405
406//____________________________________________________________________
407const Char_t*
408AliFMDCorrELossFit::ELossFit::GetName() const
409{
7984e5f7 410 //
411 // Get the name of this object
412 //
413 //
414 // Return:
415 //
416 //
0bd4b00f 417 return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
418}
419
420//____________________________________________________________________
421void
422AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
423{
7984e5f7 424 //
425 // Browse this object
426 //
427 // Parameters:
428 // b Browser
429 //
0bd4b00f 430 Draw(b ? b->GetDrawOption() : "comp");
431 gPad->SetLogy();
432 gPad->Update();
433}
434
435//____________________________________________________________________
436void
437AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
438{
7984e5f7 439 //
440 // Draw this fit
441 //
442 // Parameters:
443 // option Options
444 // - COMP Draw components too
445 //
0bd4b00f 446 TString opt(option);
447 opt.ToUpper();
448 bool comp = false;
449 if (opt.Contains("COMP")) {
450 opt.ReplaceAll("COMP","");
451 comp = true;
452 }
453 if (!opt.Contains("SAME")) {
454 gPad->Clear();
455 }
456
457 TObjArray cleanup;
458 TF1* tot = AliForwardUtil::MakeNLandauGaus(1,
459 fDelta, fXi,
460 fSigma, fSigmaN,
461 fN, fA,
462 0.01, 10);
463 tot->SetLineColor(kBlack);
464 tot->SetLineWidth(2);
465 tot->SetLineStyle(1);
466 tot->SetTitle(GetName());
467 Double_t max = tot->GetMaximum();
468
469 if (!opt.Contains("SAME")) {
470 TH1* frame = new TH1F(GetName(),
471 Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
472 100, 0, 10);
473 frame->SetMinimum(max/10000);
474 frame->SetMaximum(max*1.4);
475 frame->SetXTitle("#Delta / #Delta_{mip}");
476 frame->Draw();
477 opt.Append(" SAME");
478 }
479 tot->DrawCopy(opt.Data());
480 cleanup.Add(tot);
481
482 if (!comp) {
483 gPad->cd();
484 return;
485 }
486
487 Double_t min = max;
488 opt.Append(" same");
489 Int_t maxW = FindMaxWeight();
490 for (Int_t i=1; i <= fN; i++) {
491 TF1* f = AliForwardUtil::MakeILandauGaus((i == 1 ? 1 : fA[i-2]),
492 fDelta, fXi,
493 fSigma, fSigmaN,
494 i, 0.01, 10);
495 f->SetLineWidth(2);
496 f->SetLineStyle(i > maxW ? 2 : 1);
497 min = TMath::Min(f->GetMaximum(), min);
498 f->DrawCopy(opt.Data());
499 cleanup.Add(f);
500 }
501 min /= 100;
502 tot->GetHistogram()->SetMaximum(max);
503 tot->GetHistogram()->SetMinimum(min);
504 tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
505
506 gPad->cd();
507}
508
509
510//____________________________________________________________________
511#define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
512
513//____________________________________________________________________
514void
515AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
516 Double_t maxRelError,
517 Double_t leastWeight)
518{
7984e5f7 519 //
520 // Calculate the quality
521 //
0bd4b00f 522 Int_t qual = 0;
523 if (fNu > 0 && fChi2 / fNu < maxChi2nu) qual += 4;;
524 if (CHECKPAR(fDelta, fEDelta, maxRelError)) qual++;
525 if (CHECKPAR(fXi, fEXi, maxRelError)) qual++;
526 if (CHECKPAR(fSigma, fESigma, maxRelError)) qual++;
527 if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
528 qual += FindMaxWeight(1.5*maxRelError, leastWeight, fN);
529 fQuality = qual;
530}
531
532//____________________________________________________________________
533AliFMDCorrELossFit::AliFMDCorrELossFit()
534 : TObject(),
535 fRings(),
536 fEtaAxis(0,0,0),
537 fLowCut(0)
538{
7984e5f7 539 //
540 // Default constructor
541 //
0bd4b00f 542 fRings.SetOwner(kTRUE);
543 fEtaAxis.SetTitle("#eta");
544 fEtaAxis.SetName("etaAxis");
545 fRings.SetName("rings");
546}
547
548//____________________________________________________________________
549AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
550 : TObject(o),
551 fRings(o.fRings),
552 fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
553 fLowCut(0)
554{
7984e5f7 555 //
556 // Copy constructor
557 //
558 // Parameters:
559 // o Object to copy from
560 //
0bd4b00f 561 fEtaAxis.SetTitle("#eta");
562 fEtaAxis.SetName("etaAxis");
563}
564
565//____________________________________________________________________
566AliFMDCorrELossFit::~AliFMDCorrELossFit()
567{
7984e5f7 568 //
569 // Destructor
570 //
0bd4b00f 571 fRings.Clear();
572}
573
574//____________________________________________________________________
575AliFMDCorrELossFit&
576AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
577{
7984e5f7 578 //
579 // Assignment operator
580 //
581 // Parameters:
582 // o Object to assign from
583 //
584 // Return:
585 // Reference to this object
586 //
0bd4b00f 587 fRings = o.fRings;
588 fLowCut = o.fLowCut;
589 SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
590
591 return *this;
592}
593//____________________________________________________________________
594Int_t
595AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
596{
7984e5f7 597 //
598 // Find the eta bin corresponding to the given eta
599 //
600 // Parameters:
601 // eta Eta value
602 //
603 // Return:
604 // Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
605 // eta, or 0 if out of range.
606 //
0bd4b00f 607 if (fEtaAxis.GetXmin() == fEtaAxis.GetXmax() || fEtaAxis.GetNbins() == 0) {
608 AliWarning("No eta axis defined");
609 return -1;
610 }
611 Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
612 if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
613 return bin;
614}
615
616//____________________________________________________________________
617Bool_t
618AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
619{
7984e5f7 620 //
621 // Set the fit parameters from a function
622 //
623 // Parameters:
624 // d Detector
625 // r Ring
626 // etaBin Eta (bin number, 1->nBins)
627 // f ELoss fit result - note, the object will take ownership
628 //
0bd4b00f 629 TObjArray* ringArray = GetOrMakeRingArray(d, r);
630 if (!ringArray) {
631 AliError(Form("Failed to make ring array for FMD%d%c", d, r));
632 return kFALSE;
633 }
634 if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
635 AliError(Form("bin=%d is out of range [%d,%d]",
636 etaBin, 1, fEtaAxis.GetNbins()));
637 return kFALSE;
638 }
639 // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
640 ringArray->AddAtAndExpand(fit, etaBin);
641 return kTRUE;
642}
643
644//____________________________________________________________________
645Bool_t
646AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
647{
7984e5f7 648 //
649 // Set the fit parameters from a function
650 //
651 // Parameters:
652 // d Detector
653 // r Ring
654 // eta Eta
655 // f ELoss fit result - note, the object will take ownership
656 //
0bd4b00f 657 Int_t bin = FindEtaBin(eta);
658 if (bin <= 0) {
659 AliError(Form("eta=%f is out of range [%f,%f]",
660 eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
661 return kFALSE;
662 }
663
664 return SetFit(d, r, bin, fit);
665}
666//____________________________________________________________________
667Bool_t
668AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r,
669 Double_t eta,
670 Int_t quality,UShort_t n,
671 Double_t chi2, UShort_t nu,
672 Double_t c, Double_t ec,
673 Double_t delta, Double_t edelta,
674 Double_t xi, Double_t exi,
675 Double_t sigma, Double_t esigma,
676 Double_t sigman, Double_t esigman,
677 Double_t* a, Double_t* ea)
678{
7984e5f7 679 //
680 // Set the fit parameters from a function
681 //
682 // Parameters:
683 // d Detector number
684 // r Ring identifier
685 // eta Eta value
686 // quality Quality flag
687 // n @f$ N@f$ - Number of fitted peaks
688 // chi2 @f$ \chi^2 @f$
689 // nu @f$ \nu @f$ - number degrees of freedom
690 // c @f$ C@f$ - scale constant
691 // ec @f$ \delta C@f$ - error on @f$ C@f$
692 // delta @f$ \Delta@f$ - most probable value
693 // edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
694 // xi @f$ \xi@f$ - Landau width
695 // exi @f$ \delta\xi@f$ - error on @f$\xi@f$
696 // sigma @f$ \sigma@f$ - Gaussian width
697 // esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
698 // sigman @f$ \sigma_n@f$ - Noise width
699 // esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
700 // a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
701 // @f$ i=2,\ldots@f$
702 // ea Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for
703 // @f$ i=2,\ldots@f$
704 //
0bd4b00f 705 ELossFit* e = new ELossFit(quality, n,
706 chi2, nu,
707 c, ec,
708 delta, edelta,
709 xi, exi,
710 sigma, esigma,
711 sigman, esigman,
712 a, ea);
713 if (!SetFit(d, r, eta, e)) {
714 delete e;
715 return kFALSE;
716 }
717 return kTRUE;
718}
719//____________________________________________________________________
720Bool_t
721AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta,
722 Int_t quality, const TF1& f)
723{
7984e5f7 724 //
725 // Set the fit parameters from a function
726 //
727 // Parameters:
728 // d Detector
729 // r Ring
730 // eta Eta
731 // quality Quality flag
732 // f Function from fit
733 //
0bd4b00f 734 ELossFit* e = new ELossFit(quality, f);
735 if (!SetFit(d, r, eta, e)) {
736 delete e;
737 return kFALSE;
738 }
739 return kTRUE;
740}
741//____________________________________________________________________
742AliFMDCorrELossFit::ELossFit*
743AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Int_t etabin) const
744{
7984e5f7 745 //
746 // Find the fit corresponding to the specified parameters
747 //
748 // Parameters:
749 // d Detector
750 // r Ring
751 // etabin Eta bin (1 based)
752 //
753 // Return:
754 // Fit parameters or null in case of problems
755 //
0bd4b00f 756 TObjArray* ringArray = GetRingArray(d, r);
757 if (!ringArray) {
758 AliError(Form("Failed to make ring array for FMD%d%c", d, r));
759 return 0;
760 }
761 if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
762 AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
763 etabin, 1, fEtaAxis.GetNbins(), d, r));
764 return 0;
765 }
766 if (etabin > ringArray->GetEntriesFast()) {
767 AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
768 etabin, 1, ringArray->GetEntriesFast(), d, r));
769 return 0;
770 }
771 else if (etabin >= ringArray->GetEntriesFast()) {
772 AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, "
773 "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
774 etabin-1));
775 etabin--;
776 }
777 else if (!ringArray->At(etabin)) {
778 AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d",
779 etabin, d, r, etabin+1));
780 etabin++;
781 }
782 return static_cast<ELossFit*>(ringArray->At(etabin));
783}
784//____________________________________________________________________
785AliFMDCorrELossFit::ELossFit*
786AliFMDCorrELossFit::FindFit(UShort_t d, Char_t r, Double_t eta) const
787{
7984e5f7 788 //
789 // Find the fit corresponding to the specified parameters
790 //
791 // Parameters:
792 // d Detector
793 // r Ring
794 // eta Eta value
795 //
796 // Return:
797 // Fit parameters or null in case of problems
798 //
0bd4b00f 799 Int_t etabin = FindEtaBin(eta);
800 return FindFit(d, r, etabin);
801}
802//____________________________________________________________________
803TObjArray*
804AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
805{
7984e5f7 806 //
807 // Get the ring array corresponding to the specified ring
808 //
809 // Parameters:
810 // d Detector
811 // r Ring
812 //
813 // Return:
814 // Pointer to ring array, or null in case of problems
815 //
0bd4b00f 816 Int_t idx = -1;
817 switch (d) {
818 case 1: idx = 0; break;
819 case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
820 case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
821 }
822 if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
823 return static_cast<TObjArray*>(fRings.At(idx));
824}
825//____________________________________________________________________
826TObjArray*
827AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
828{
7984e5f7 829 //
830 // Get the ring array corresponding to the specified ring
831 //
832 // Parameters:
833 // d Detector
834 // r Ring
835 //
836 // Return:
837 // Pointer to ring array, or newly created container
838 //
0bd4b00f 839 Int_t idx = -1;
840 switch (d) {
841 case 1: idx = 0; break;
842 case 2: idx = (r == 'i' || r == 'I') ? 1 : 2; break;
843 case 3: idx = (r == 'o' || r == 'I') ? 3 : 4; break;
844 }
845 if (idx < 0) return 0;
846 if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
847 TObjArray* a = new TObjArray(0);
848 // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
849 a->SetName(Form("FMD%d%c", d, r));
850 a->SetOwner();
851 fRings.AddAtAndExpand(a, idx);
852 }
853 return static_cast<TObjArray*>(fRings.At(idx));
854}
855
856namespace {
857 TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
858 Int_t color)
859 {
860 TH1D* h = new TH1D(Form("%s_%s", name, title),
861 Form("%s %s", name, title),
862 axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
863 h->SetDirectory(0);
864 h->SetMarkerStyle(20);
865 h->SetMarkerColor(color);
866 h->SetMarkerSize(0.5);
867 h->SetFillColor(color);
868 h->SetFillStyle(3001);
869 h->SetLineColor(color);
870 return h;
871 }
872}
873
874//____________________________________________________________________
875void
876AliFMDCorrELossFit::Draw(Option_t* option)
877{
7984e5f7 878 //
879 // Draw this object
880 //
881 // Parameters:
882 // option Options. Possible values are
883 // - err Plot error bars
884 //
0bd4b00f 885 TString opt(Form("nostack %s", option));
886 opt.ToLower();
887 Bool_t rel = (opt.Contains("rel"));
888 Bool_t err = (opt.Contains("err"));
889 if (rel) opt.ReplaceAll("rel","");
890 if (err) opt.ReplaceAll("err","");
891 Int_t nRings = fRings.GetEntriesFast();
892 UShort_t maxN = 0;
893 for (Int_t i = 0; i < nRings; i++) {
894 if (!fRings.At(i)) continue;
895 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
896 Int_t nFits = a->GetEntriesFast();
897
898 for (Int_t j = 0; j < nFits; j++) {
899 ELossFit* fit = static_cast<ELossFit*>(a->At(j));
900 if (!fit) continue;
901 maxN = TMath::Max(maxN, UShort_t(fit->fN));
902 }
903 }
904 // AliInfo(Form("Maximum N is %d", maxN));
905 Int_t nPad = 7+maxN-1; // 7 regular params, and maxN-1 weights
906 TVirtualPad* pad = gPad;
907 pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
908
909 THStack* chi2nu;
910 THStack* c;
911 THStack* delta;
912 THStack* xi;
913 THStack* sigma;
914 THStack* sigman;
915 THStack* n;
916 TList stacks;
917 stacks.AddAt(chi2nu= new THStack("chi2", "#chi^{2}/#nu"), 0);
918 stacks.AddAt(c = new THStack("c", "C"), 1);
919 stacks.AddAt(delta = new THStack("delta", "#Delta_{mp}"), 2);
920 stacks.AddAt(xi = new THStack("xi", "#xi"), 3);
921 stacks.AddAt(sigma = new THStack("sigma", "#sigma"), 4);
922 stacks.AddAt(sigman= new THStack("sigman", "#sigma_{n}"), 5);
923 stacks.AddAt(n = new THStack("n", "N"), 6);
924 for (Int_t i = 1; i <= maxN; i++) {
925 stacks.AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), 6+i);
926 }
927
928 for (Int_t i = 0; i < nRings; i++) {
929 if (!fRings.At(i)) continue;
930 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
931 Int_t nFits = a->GetEntriesFast();
932 Int_t color = 0;
933 switch (i) {
934 case 0: color = kRed+2; break;
935 case 1: color = kGreen+2; break;
936 case 2: color = kGreen-2; break;
937 case 3: color = kBlue+2; break;
938 case 4: color = kBlue-2; break;
939 }
940
941 TH1D* hChi = MakeHist(fEtaAxis,a->GetName(), "chi2", color);
942 TH1D* hC = MakeHist(fEtaAxis,a->GetName(), "c", color);
943 TH1D* hDelta = MakeHist(fEtaAxis,a->GetName(), "delta", color);
944 TH1D* hXi = MakeHist(fEtaAxis,a->GetName(), "xi", color);
945 TH1D* hSigma = MakeHist(fEtaAxis,a->GetName(), "sigma", color);
946 TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
947 TH1D* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
948 TH1D* hA[maxN];
949 const char* ho = (rel || !err ? "hist" : "e");
950 chi2nu->Add(hChi, "hist"); // 0
951 c ->Add(hC, ho); // 1
952 delta ->Add(hDelta, ho); // 2
953 xi ->Add(hXi, ho); // 3
954 sigma ->Add(hSigma, ho); // 4
955 sigman->Add(hSigmaN,ho); // 5
956 n ->Add(hN, "hist"); // 6
957 hChi->SetFillColor(color);
958 hChi->SetFillStyle(3001);
959 hN->SetFillColor(color);
960 hN->SetFillStyle(3001);
961
962 for (Int_t k = 1; k <= maxN; k++) {
963 hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
964 static_cast<THStack*>(stacks.At(6+k))->Add(hA[k-1], ho);
965 }
966
967 for (Int_t j = 0; j < nFits; j++) {
968 ELossFit* f = static_cast<ELossFit*>(a->At(j));
969 if (!f) continue;
970
971 Int_t b = f->fBin;
972 Int_t nW = f->FindMaxWeight();
973 hChi ->SetBinContent(b, (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
974 hN ->SetBinContent(b, nW);
975
976 if (rel) {
977 hC ->SetBinContent(b, f->fC >0 ?f->fEC /f->fC :0);
978 hDelta ->SetBinContent(b, f->fDelta >0 ?f->fEDelta /f->fDelta :0);
979 hXi ->SetBinContent(b, f->fXi >0 ?f->fEXi /f->fXi :0);
980 hSigma ->SetBinContent(b, f->fSigma >0 ?f->fESigma /f->fSigma :0);
981 hSigmaN->SetBinContent(b, f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0);
982 }
983 else {
984 hC ->SetBinContent(b, f->fC);
985 hDelta ->SetBinContent(b, f->fDelta);
986 hXi ->SetBinContent(b, f->fXi);
987 hSigma ->SetBinContent(b, f->fSigma);
988 hSigmaN->SetBinContent(b, f->fSigmaN);
989 hC ->SetBinError(b, f->fEC);
990 hDelta ->SetBinError(b, f->fEDelta);
991 hXi ->SetBinError(b, f->fEXi);
992 hSigma ->SetBinError(b, f->fESigma);
993 hSigmaN->SetBinError(b, f->fESigmaN);
994 }
995 for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
996 if (rel)
997 hA[k]->SetBinContent(b, f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0);
998 else {
999 hA[k]->SetBinContent(b, f->fA[k]);
1000 hA[k]->SetBinError(b, f->fEA[k]);
1001 }
1002 }
1003 }
1004 }
1005 Int_t nPad2 = (nPad+1) / 2;
1006 for (Int_t i = 0; i < nPad; i++) {
1007 Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
1008 TVirtualPad* p = pad->cd(iPad);
1009 p->SetLeftMargin(.15);
1010 p->SetFillColor(0);
1011 p->SetFillStyle(0);
1012 p->SetGridx();
1013 p->SetGridy();
1014 if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
1015
1016 THStack* stack = static_cast<THStack*>(stacks.At(i));
1017 // AliInfo(Form("Drawing %s (%d) in pad %d", stack->GetName(), i, iPad));
1018 stack->Draw(opt.Data());
1019
1020 TString tit(stack->GetTitle());
1021 if (rel && i != 0 && i != 6)
1022 tit = Form("#delta %s/%s", tit.Data(), tit.Data());
1023 TH1* hist = stack->GetHistogram();
1024 TAxis* yaxis = hist->GetYaxis();
1025 yaxis->SetTitle(tit.Data());
1026 yaxis->SetTitleSize(0.15);
1027 yaxis->SetLabelSize(0.08);
1028 yaxis->SetTitleOffset(0.35);
1029 yaxis->SetNdivisions(5);
1030
1031 TAxis* xaxis = stack->GetHistogram()->GetXaxis();
1032 xaxis->SetTitle("#eta");
1033 xaxis->SetTitleSize(0.15);
1034 xaxis->SetLabelSize(0.08);
1035 xaxis->SetTitleOffset(0.35);
1036 xaxis->SetNdivisions(10);
1037
1038
1039 if (!rel) {
1040 switch (i) {
1041 case 0: break; // chi^2/nu
1042 case 1: break; // C
1043 case 2: stack->SetMinimum(0.4); break; // Delta
1044 case 3: stack->SetMinimum(0.02);break; // xi
1045 case 4: stack->SetMinimum(0.05);break; // sigma
1046 case 5: break; // sigmaN
1047 case 6:
1048 stack->SetMinimum(-.1);
1049 stack->SetMaximum(maxN+.5);
1050 break; // N
1051 default: break; // Weights
1052 }
1053 }
1054 stack->DrawClone(opt.Data());
1055 }
1056 pad->cd();
1057}
1058
1059//____________________________________________________________________
1060void
1061AliFMDCorrELossFit::Print(Option_t* option) const
1062{
1063 TString opt(option);
1064 opt.ToUpper();
1065 Int_t nRings = fRings.GetEntriesFast();
1066 bool recurse = opt.Contains("R");
1067
1068 std::cout << "Low cut in fit range: " << fLowCut << "\n"
1069 << "Eta axis: " << fEtaAxis.GetNbins()
1070 << " bins, range [" << fEtaAxis.GetXmin() << ","
1071 << fEtaAxis.GetXmax() << "]" << std::endl;
1072
1073 for (Int_t i = 0; i < nRings; i++) {
1074 if (!fRings.At(i)) continue;
1075 TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1076 Int_t nFits = a->GetEntriesFast();
1077
1078 std::cout << a->GetName() << " [" << nFits << " entries]"
1079 << (recurse ? ":\n" : "\t");
1080 Int_t min = fEtaAxis.GetNbins()+1;
1081 Int_t max = 0;
1082 for (Int_t j = 0; j < nFits; j++) {
1083 if (!a->At(j)) continue;
1084
1085 min = TMath::Min(j, min);
1086 max = TMath::Max(j, max);
1087
1088 if (recurse) {
1089 std::cout << "Bin # " << j << "\t";
1090 ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1091 fit->Print(option);
1092 }
1093 }
1094 if (!recurse)
1095 std::cout << " bin range: " << std::setw(3) << min
1096 << "-" << std::setw(3) << max << " " << std::setw(3)
1097 << (max-min+1) << " bins" << std::endl;
1098 }
1099}
1100
1101//____________________________________________________________________
1102void
1103AliFMDCorrELossFit::Browse(TBrowser* b)
1104{
7984e5f7 1105 //
1106 // Browse this object
1107 //
1108 // Parameters:
1109 // b
1110 //
0bd4b00f 1111 b->Add(&fRings);
1112 b->Add(&fEtaAxis);
1113}
1114
1115
1116
1117//____________________________________________________________________
1118//
1119// EOF
1120//