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