AliPhysics  9c66e61 (9c66e61)
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),
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;
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;
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 //____________________________________________________________________
459 const Char_t*
461 {
462  //
463  // Get the name of this object
464  //
465  //
466  // Return:
467  //
468  //
469  return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
470 }
471 
472 //____________________________________________________________________
473 void
475 {
476  //
477  // Browse this object
478  //
479  // Parameters:
480  // b Browser
481  //
482  Draw(b ? b->GetDrawOption() : "comp values");
483  gPad->SetLogy();
484  gPad->Update();
485 }
486 
487 //____________________________________________________________________
488 void
490 {
491  //
492  // Draw this fit
493  //
494  // Parameters:
495  // option Options
496  // - COMP Draw components too
497  //
498  TString opt(option);
499  opt.ToUpper();
500  bool comp = false;
501  bool good = false;
502  bool vals = false;
503  bool legd = false;
504  bool peak = false;
505  if (opt.Contains("COMP")) {
506  opt.ReplaceAll("COMP","");
507  comp = true;
508  }
509  if (opt.Contains("GOOD")) {
510  opt.ReplaceAll("GOOD","");
511  good = true;
512  }
513  if (opt.Contains("VALUES")) {
514  opt.ReplaceAll("VALUES","");
515  vals = true;
516  }
517  if (opt.Contains("LEGEND")) {
518  opt.ReplaceAll("LEGEND","");
519  legd = comp;
520  }
521  if (!opt.Contains("SAME")) {
522  gPad->Clear();
523  }
524  if (opt.Contains("PEAK")) {
525  peak = true;
526  }
527  TLegend* l = 0;
528  if (legd) {
529  l = new TLegend(.3, .5, .59, .94);
530  l->SetBorderSize(0);
531  l->SetFillColor(0);
532  l->SetFillStyle(0);
533  }
534  TObjArray cleanup;
535  Int_t maxW = FindMaxWeight();
536  TF1* tot = AliLandauGaus::MakeFn(fC * 1, fDelta, fXi, fSigma, fSigmaN,
537  maxW/*fN*/, fA, 0.01, 10);
538  tot->SetLineColor(kBlack);
539  tot->SetLineWidth(2);
540  tot->SetLineStyle(1);
541  tot->SetTitle(GetName());
542  if (l) l->AddEntry(tot, "Total", "l");
543  Double_t max = tot->GetMaximum();
544 
545 
546  if (!opt.Contains("SAME")) {
547  TH1* frame = new TH1F(GetName(),
548  Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
549  100, 0, 10);
550  frame->SetMinimum(max/10000);
551  frame->SetMaximum(max*1.4);
552  frame->SetXTitle("#Delta / #Delta_{mip}");
553  frame->Draw();
554  opt.Append(" SAME");
555  }
556  TF1* cpy = tot->DrawCopy(opt.Data());
557  cleanup.Add(tot);
558 
559  if (vals) {
560  Double_t x1 = .72;
561  Double_t x2 = .73;
562  Double_t y = .90;
563  Double_t dy = .05;
564  TLatex* ltx1 = new TLatex(x1, y, "");
565  TLatex* ltx2 = new TLatex(x2, y, "");
566  ltx1->SetNDC();
567  ltx1->SetTextAlign(33);
568  ltx1->SetTextFont(132);
569  ltx1->SetTextSize(dy-.01);
570  ltx2->SetNDC();
571  ltx2->SetTextAlign(13);
572  ltx2->SetTextFont(132);
573  ltx2->SetTextSize(dy-.01);
574 
575  ltx1->DrawLatex(x1, y, "Quality");
576  ltx2->DrawLatex(x2, y, Form("%d", fQuality));
577  y -= dy;
578 
579  ltx1->DrawLatex(x1, y, "#chi^{2}/#nu");
580  ltx2->DrawLatex(x2, y, Form("%7.3f", (fNu > 0 ? fChi2 / fNu : -1)));
581  y -= dy;
582 
583  const Char_t* pn[] = { "C", "#Delta", "#xi", "#sigma" };
584  Double_t pv[] = { fC, fDelta, fXi, fSigma };
585  Double_t pe[] = { fEC, fEDelta, fEXi, fESigma };
586  for (Int_t i = 0; i < 4; i++) {
587  ltx1->DrawLatex(x1, y, pn[i]);
588  ltx2->DrawLatex(x2, y, Form("%6.4f#pm%6.4f", pv[i], pe[i]));
589  y -= dy;
590  }
591  for (Int_t i=2; i <= fN; i++) {
592  if (i > maxW) {
593  ltx1->SetTextColor(kRed+3);
594  ltx2->SetTextColor(kRed+3);
595  }
596  ltx1->DrawLatex(x1, y, Form("a_{%d}", i));
597  ltx2->DrawLatex(x2, y, Form("%6.4f#pm%6.4f", fA[i-2], fEA[i-2]));
598  y -= dy;
599  }
600 
601  }
602 
603  if (peak) {
604  TLine* pl = new TLine(fDelta, 0.01*max, fDelta, 1.5*max);
605  pl->SetLineStyle(2);
606  pl->SetLineColor(kBlack);
607  pl->Draw();
608  }
609  if (!comp) {
610  gPad->cd();
611  return;
612  }
613 
614  Double_t min = max;
615  opt.Append(" same");
616  for (Int_t i=1; i <= fN; i++) {
617  if (good && i > maxW) break;
618  TF1* f = AliLandauGaus::MakeFi(fC*(i == 1 ? 1 : fA[i-2]),
619  fDelta, fXi,
620  fSigma, fSigmaN,
621  i, 0.01, 10);
622  f->SetLineWidth(2);
623  f->SetLineStyle(i > maxW ? 2 : 1);
624  min = TMath::Min(f->GetMaximum(), min);
625  f->DrawCopy(opt.Data());
626  if (l) l->AddEntry(f, Form("%d MIP%s", i, (i>1 ? "s" : "")), "l");
627 
628  cleanup.Add(f);
629  }
630  min /= 100;
631  if (max <= 0) max = 0.1;
632  if (min <= 0) min = 1e-4;
633  cpy->SetMaximum(max);
634  cpy->SetMinimum(min);
635  cpy->GetHistogram()->SetMaximum(max);
636  cpy->GetHistogram()->SetMinimum(min);
637  cpy->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
638  if (l) l->Draw();
639 
640  gPad->cd();
641 }
642 
643 
644 //____________________________________________________________________
645 #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
646 
647 //____________________________________________________________________
648 Double_t
650 {
651  return f * fDelta;
652 }
653 //____________________________________________________________________
654 Double_t
656 {
657  return fDelta - f * fXi;
658 }
659 //____________________________________________________________________
660 Double_t
662 {
663  return fDelta - f * (fXi+fSigma);
664 }
665 
666 //____________________________________________________________________
667 Double_t
669 {
670  Double_t sum = fXi / (fEXi*fEXi) + fSigma / (fESigma+fESigma);
671  Double_t nor = 1 / (fEXi*fEXi) + 1 / (fESigma+fESigma);
672  return fDelta - f * sum / nor;
673 }
674 
675 //____________________________________________________________________
676 Double_t
678 {
679  Double_t ret = 1000;
680  TF1* f = 0;
681  TGraph* g = 0;
682  try {
683  if (!(f = GetF1(1,5))) // First landau up to Delta/Delta_{mip}=5
684  throw TString("Didn't TF1 object");
685  if (!(g = new TGraph(f, "i")))
686  throw TString("Failed to integrate function");
687 
688  Int_t n = g->GetN();
689  Double_t total = g->GetY()[n-1];
690  if (total <= 0)
691  throw TString::Format("Invalid integral: %lf", total);
692 
693  for (Int_t i = 0; i < n; i++) {
694  Double_t prob = g->GetY()[i] / total;
695  if (prob > low) {
696  ret = g->GetX()[i];
697  break;
698  }
699  }
700  if (ret >= 1000)
701  throw TString::Format("Couldn't find x value for cut %lf", low);
702  }
703  catch (const TString& str) {
704  AliWarningF("%s: %s", GetName(), str.Data());
705  }
706  if (f) delete f;
707  if (g) delete g;
708  return ret;
709 }
710 
711 //____________________________________________________________________
712 void
714  Double_t maxRelError,
715  Double_t leastWeight)
716 {
717  //
718  // Calculate the quality
719  //
720  Double_t decline = maxChi2nu;
721  Int_t qual = 0;
722  if (fNu > 0) {
723  Double_t red = fChi2 / fNu;
724  if (red < maxChi2nu) qual += 4;
725  else {
726  Int_t q = Int_t((maxChi2nu+decline - red) / decline * 4);
727  if (q > 0) qual += q;
728  }
729  }
730  if (CHECKPAR(fDelta, fEDelta, maxRelError)) qual++;
731  if (CHECKPAR(fXi, fEXi, maxRelError)) qual++;
732  if (CHECKPAR(fSigma, fESigma, maxRelError)) qual++;
733  if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
734  // Large penalty for large sigma - often a bad fit. (factor was 10)
735  if (fSigma > 4*fXi) qual -= 4;
736  if (fXi < 0.01) qual -= 2;
737  if (fSigma < 0.01) qual -= 2;
738  //Info("CalculateQuality", "FMD%d%c [%3d] sigma=%f xi=%f sigma/xi=%f qual=%d",
739  // fDet, fRing, fBin, fSigma, fXi, fSigma/fXi, qual);
740  qual += FindMaxWeight(2*maxRelError, leastWeight, fN);
741  fQuality = qual;
742 }
743 
744 //____________________________________________________________________
746  : TObject(),
747  fRings(),
748  fEtaAxis(0,0,0),
749  fLowCut(0),
750  fCache(0)
751 {
752  //
753  // Default constructor
754  //
755  fRings.SetOwner(kTRUE);
756  fEtaAxis.SetTitle("#eta");
757  fEtaAxis.SetName("etaAxis");
758  fRings.SetName("rings");
759  SetBit(kIsGoodAsserted,false);
760 }
761 
762 //____________________________________________________________________
764  : TObject(o),
765  fRings(o.fRings),
766  fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
767  fLowCut(0),
768  fCache(0)
769 {
770  //
771  // Copy constructor
772  //
773  // Parameters:
774  // o Object to copy from
775  //
776  fEtaAxis.SetTitle("#eta");
777  fEtaAxis.SetName("etaAxis");
778  SetBit(kIsGoodAsserted,false);
779 }
780 
781 //____________________________________________________________________
783 {
784  //
785  // Destructor
786  //
787  fRings.Clear();
788 }
789 
790 //____________________________________________________________________
793 {
794  //
795  // Assignment operator
796  //
797  // Parameters:
798  // o Object to assign from
799  //
800  // Return:
801  // Reference to this object
802  //
803  if (&o == this) return *this;
804  fRings = o.fRings;
805  fLowCut = o.fLowCut;
806  fCache = o.fCache;
807  SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
808 
809  return *this;
810 }
811 #define CACHE(BIN,IDX,OFFSET) fCache[IDX*OFFSET+BIN-1]
812 #define CACHEDR(BIN,D,R,OFFSET) \
813  CACHE(BIN,(D == 1 ? 0 : (D - 2) * 2 + 1 + (R=='I' || R=='i' ? 0 : 1)),OFFSET)
814 
815 //____________________________________________________________________
816 void
818 {
820  if (fCache.GetSize() > 0) return;
821 
822  Int_t nRings = fRings.GetEntriesFast();
823  Int_t offset = fEtaAxis.GetNbins();
824 
825  fCache.Set(nRings*offset);
826  fCache.Reset(-1);
827 
828  for (Int_t i = 0; i < nRings; i++) {
829  TObjArray* ringArray = static_cast<TObjArray*>(fRings.At(i));
830 
831  // First loop to find where we actually have fits
832  Int_t nFits = 0;
833  Int_t nGood = 0;
834  Int_t minBin = offset+1;
835  Int_t maxBin = -1;
836  Int_t realMinBin = offset+1;
837  Int_t realMaxBin = -1;
838  for (Int_t j = 1; j < ringArray->GetEntriesFast(); j++) {
839  ELossFit* fit = static_cast<ELossFit*>(ringArray->At(j));
840  if (!fit) continue;
841  nFits++;
842 
843  // Update our range off possible fits
844  realMinBin = TMath::Min(j, realMinBin);
845  realMaxBin = TMath::Max(j, realMaxBin);
846 
847  // Check the quality of the fit
851  if (minQuality > 0 && fit->fQuality < minQuality) continue;
852  nGood++;
853 
854  // Check this bin
855  CACHE(j,i,offset) = j;
856  minBin = TMath::Min(j, minBin);
857  maxBin = TMath::Max(j, maxBin);
858  }
859  AliInfoF("Out of %d bins, %d had fits, of which %d are good (%5.1f%%)",
860  offset, nFits, nGood, (nFits > 0 ? 100*float(nGood)/nFits : 0));
861 
862  // Now loop and set neighbors
863  realMinBin = TMath::Max(1, realMinBin-1); // Include one more
864  realMaxBin = TMath::Min(offset, realMaxBin+1); // Include one more
865  for (Int_t j = realMinBin; j <= realMaxBin; j++) {
866  if (CACHE(j,i,offset) > 0) continue;
867 
868  Int_t nK = TMath::Max(realMaxBin - j, j - realMinBin);
869  Int_t found = -1;
870  for (Int_t k = 1; k <= nK; k++) {
871  Int_t left = j - k;
872  Int_t right = j + k;
873  if (left > realMinBin &&
874  CACHE(left,i,offset) == left) found = left;
875  else if (right < realMaxBin &&
876  CACHE(right,i,offset) == right) found = right;
877  if (found > 0) break;
878  }
879  // Now check that we found something
880  if (found) CACHE(j,i,offset) = CACHE(found,i,offset);
881  else AliWarningF("No fit found for etabin=%d in ring=%d", j, i);
882  }
883  }
884 }
885 
886 
887 //____________________________________________________________________
888 Int_t
890 {
891  //
892  // Find the eta bin corresponding to the given eta
893  //
894  // Parameters:
895  // eta Eta value
896  //
897  // Return:
898  // Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
899  // eta, or 0 if out of range.
900  //
901  if (TMath::Abs(fEtaAxis.GetXmin() - fEtaAxis.GetXmax()) < 1e-6
902  || fEtaAxis.GetNbins() == 0) {
903  AliWarning("No eta axis defined");
904  return -1;
905  }
906  Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
907  if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
908  return bin;
909 }
910 
911 //____________________________________________________________________
912 Bool_t
914 {
915  //
916  // Set the fit parameters from a function
917  //
918  // Parameters:
919  // d Detector
920  // r Ring
921  // etaBin Eta (bin number, 1->nBins)
922  // f ELoss fit result - note, the object will take ownership
923  //
924  TObjArray* ringArray = GetOrMakeRingArray(d, r);
925  if (!ringArray) {
926  AliError(Form("Failed to make ring array for FMD%d%c", d, r));
927  return kFALSE;
928  }
929  if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
930  AliError(Form("bin=%d is out of range [%d,%d]",
931  etaBin, 1, fEtaAxis.GetNbins()));
932  return kFALSE;
933  }
934  // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
935  ringArray->AddAtAndExpand(fit, etaBin);
936  return kTRUE;
937 }
938 
939 //____________________________________________________________________
940 Bool_t
942 {
943  //
944  // Set the fit parameters from a function
945  //
946  // Parameters:
947  // d Detector
948  // r Ring
949  // eta Eta
950  // f ELoss fit result - note, the object will take ownership
951  //
952  Int_t bin = FindEtaBin(eta);
953  if (bin <= 0) {
954  AliError(Form("eta=%f is out of range [%f,%f]",
955  eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
956  return kFALSE;
957  }
958 
959  return SetFit(d, r, bin, fit);
960 }
961 //____________________________________________________________________
962 Bool_t
964  Double_t eta,
965  Int_t quality,UShort_t n,
966  Double_t chi2, UShort_t nu,
967  Double_t c, Double_t ec,
968  Double_t delta, Double_t edelta,
969  Double_t xi, Double_t exi,
970  Double_t sigma, Double_t esigma,
971  Double_t sigman, Double_t esigman,
972  Double_t* a, Double_t* ea)
973 {
974  //
975  // Set the fit parameters from a function
976  //
977  // Parameters:
978  // d Detector number
979  // r Ring identifier
980  // eta Eta value
981  // quality Quality flag
982  // n @f$ N@f$ - Number of fitted peaks
983  // chi2 @f$ \chi^2 @f$
984  // nu @f$ \nu @f$ - number degrees of freedom
985  // c @f$ C@f$ - scale constant
986  // ec @f$ \delta C@f$ - error on @f$ C@f$
987  // delta @f$ \Delta@f$ - most probable value
988  // edelta @f$ \delta\Delta@f$ - error on @f$\Delta@f$
989  // xi @f$ \xi@f$ - Landau width
990  // exi @f$ \delta\xi@f$ - error on @f$\xi@f$
991  // sigma @f$ \sigma@f$ - Gaussian width
992  // esigma @f$ \delta\sigma@f$ - error on @f$\sigma@f$
993  // sigman @f$ \sigma_n@f$ - Noise width
994  // esigman @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
995  // a Array of @f$ N-1@f$ weights @f$ a_i@f$ for
996  // @f$ i=2,\ldots@f$
997  // ea Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for
998  // @f$ i=2,\ldots@f$
999  //
1000  ELossFit* e = new ELossFit(quality, n,
1001  chi2, nu,
1002  c, ec,
1003  delta, edelta,
1004  xi, exi,
1005  sigma, esigma,
1006  sigman, esigman,
1007  a, ea);
1008  if (!SetFit(d, r, eta, e)) {
1009  delete e;
1010  return kFALSE;
1011  }
1012  return kTRUE;
1013 }
1014 //____________________________________________________________________
1015 Bool_t
1017  Int_t quality, const TF1& f)
1018 {
1019  //
1020  // Set the fit parameters from a function
1021  //
1022  // Parameters:
1023  // d Detector
1024  // r Ring
1025  // eta Eta
1026  // quality Quality flag
1027  // f Function from fit
1028  //
1029  ELossFit* e = new ELossFit(quality, f);
1030  if (!SetFit(d, r, eta, e)) {
1031  delete e;
1032  return kFALSE;
1033  }
1034  return kTRUE;
1035 }
1036 //____________________________________________________________________
1039 {
1040  //
1041  // Get the fit corresponding to the specified parameters
1042  //
1043  // Parameters:
1044  // d Detector
1045  // r Ring
1046  // etabin Eta bin (1 based)
1047  //
1048  // Return:
1049  // Fit parameters or null in case of problems
1050  //
1051  TObjArray* ringArray = GetRingArray(d, r);
1052  if (!ringArray) return 0;
1053  if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) return 0;
1054  if (etabin > ringArray->GetEntriesFast()) return 0;
1055  else if (etabin >= ringArray->GetEntriesFast()) etabin--;
1056  else if (!ringArray->At(etabin)) etabin++;
1057  return static_cast<ELossFit*>(ringArray->At(etabin));
1058 }
1059 //____________________________________________________________________
1062 {
1063  //
1064  // Find the fit corresponding to the specified parameters
1065  //
1066  // Parameters:
1067  // d Detector
1068  // r Ring
1069  // eta Eta value
1070  //
1071  // Return:
1072  // Fit parameters or null in case of problems
1073  //
1074  Int_t etabin = FindEtaBin(eta);
1075  return GetFit(d, r, etabin);
1076 }
1077 
1078 //____________________________________________________________________
1081  UShort_t minQ) const
1082 {
1083  //
1084  // Find the fit corresponding to the specified parameters
1085  //
1086  // Parameters:
1087  // d Detector
1088  // r Ring
1089  // etabin Eta bin (1 based)
1090  //
1091  // Return:
1092  // Fit parameters or null in case of problems
1093  //
1094  if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
1095  // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
1096  // etabin, 1, fEtaAxis.GetNbins(), d, r));
1097  return 0;
1098  }
1099 
1100  TObjArray* ringArray = GetRingArray(d, r);
1101  if (!ringArray) {
1102  AliError(Form("Failed to make ring array for FMD%d%c", d, r));
1103  return 0;
1104  }
1105  DMSG(fDebug, 10, "Got ringArray %s for FMD%d%c", ringArray->GetName(), d, r);
1106  if (fCache.GetSize() <= 0) CacheBins(minQ);
1107 #if 0
1108  Int_t idx = (d == 1 ? 0 :
1109  (d - 2) * 2 + 1 + (r=='I' || r=='i' ? 0 : 1));
1110 #endif
1111  Int_t bin = CACHEDR(etabin, d, r, fEtaAxis.GetNbins());
1112 
1113  if (bin < 0 || bin > ringArray->GetEntriesFast()) return 0;
1114 
1115  return static_cast<ELossFit*>(ringArray->At(bin));
1116 }
1117 //____________________________________________________________________
1120  UShort_t minQ) const
1121 {
1122  //
1123  // Find the fit corresponding to the specified parameters
1124  //
1125  // Parameters:
1126  // d Detector
1127  // r Ring
1128  // eta Eta value
1129  //
1130  // Return:
1131  // Fit parameters or null in case of problems
1132  //
1133  Int_t etabin = FindEtaBin(eta);
1134  return FindFit(d, r, etabin, minQ);
1135 }
1136 //____________________________________________________________________
1137 Int_t
1139 {
1140  switch (d) {
1141  case 1: return 0;
1142  case 2: return (r == 'i' || r == 'I') ? 1 : 2;
1143  case 3: return (r == 'i' || r == 'I') ? 3 : 4;
1144  }
1145  return -1;
1146 }
1147 
1148 //____________________________________________________________________
1149 TObjArray*
1151 {
1152  //
1153  // Get the ring array corresponding to the specified ring
1154  //
1155  // Parameters:
1156  // d Detector
1157  // r Ring
1158  //
1159  // Return:
1160  // Pointer to ring array, or null in case of problems
1161  //
1162  Int_t idx = GetRingIndex(d, r);
1163  if (idx < 0) return 0;
1164  return static_cast<TObjArray*>(fRings.At(idx));
1165 }
1166 //____________________________________________________________________
1167 TObjArray*
1169 {
1170  //
1171  // Get the ring array corresponding to the specified ring
1172  //
1173  // Parameters:
1174  // d Detector
1175  // r Ring
1176  //
1177  // Return:
1178  // Pointer to ring array, or newly created container
1179  //
1180  Int_t idx = GetRingIndex(d, r);
1181  if (idx < 0) return 0;
1182  if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
1183  TObjArray* a = new TObjArray(0);
1184  // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
1185  a->SetName(Form("FMD%d%c", d, r));
1186  a->SetOwner();
1187  fRings.AddAtAndExpand(a, idx);
1188  }
1189  return static_cast<TObjArray*>(fRings.At(idx));
1190 }
1191 
1192 //____________________________________________________________________
1193 Double_t
1195  Char_t r,
1196  Int_t etabin,
1197  Double_t f) const
1198 {
1199  ELossFit* fit = FindFit(d, r, etabin, 20);
1200  if (!fit) return -1024;
1201  return fit->GetMpvCut(f);
1202 }
1203 
1204 //____________________________________________________________________
1205 Double_t
1207  Char_t r,
1208  Int_t etabin,
1209  Double_t f) const
1210 {
1211  ELossFit* fit = FindFit(d, r, etabin, 20);
1212  if (!fit) return -1024;
1213  return fit->GetXiCut(f);
1214 }
1215 //____________________________________________________________________
1216 Double_t
1218  Char_t r,
1219  Int_t etabin,
1220  Double_t f) const
1221 {
1222  ELossFit* fit = FindFit(d, r, etabin, 20);
1223  if (!fit) return -1024;
1224  return fit->GetXiSigmaCut(f);
1225 }
1226 //____________________________________________________________________
1227 Double_t
1229  Char_t r,
1230  Int_t etabin,
1231  Double_t f) const
1232 {
1233  ELossFit* fit = FindFit(d, r, etabin, 20);
1234  if (!fit) return -1024;
1235  return fit->GetAvgXiSigmaCut(f);
1236 }
1237 //____________________________________________________________________
1238 Double_t
1240  Char_t r,
1241  Int_t etabin,
1242  Double_t f) const
1243 {
1244  ELossFit* fit = FindFit(d, r, etabin, 20);
1245  if (!fit) return -1024;
1246  return fit->FindProbabilityCut(f);
1247 }
1248 
1249 //____________________________________________________________________
1250 Bool_t
1252  Double_t minRate,
1253  Int_t maxGap,
1254  Int_t minInner,
1255  Int_t minOuter,
1256  Int_t minQuality)
1257 {
1258  if (TestBit(kIsGoodAsserted)) return TestBit(kIsGood);
1259 
1260  Bool_t isGood = true; // Assume success
1261 
1262  if (fRings.GetEntries() <= 0) isGood = false;
1263  TIter nextRing(&fRings);
1264  TObjArray* ringArray = 0;
1265  while ((ringArray = static_cast<TObjArray*>(nextRing()))) {
1266  Char_t r = ringArray->GetName()[4];
1267 
1268  Int_t gap = 0;
1269  Int_t good = 0;
1270  Int_t max = 0;
1271  Int_t len = 0;
1272  Bool_t first = true;
1273  for (Int_t iEta = 1; iEta <= fEtaAxis.GetNbins(); iEta++) {
1274  if (iEta > ringArray->GetEntriesFast()) {
1275  max = TMath::Max(gap, max);
1276  break;
1277  }
1278  TObject* o = ringArray->At(iEta);
1279  if (!o) {
1280  if (!first) { gap++; len++; }
1281  continue;
1282  }
1283 
1284  len++;
1285  first = false;
1286  ELossFit* fit = static_cast<ELossFit*>(o);
1287  if (fit->GetQuality() < minQuality) {
1288  gap++;
1289  continue;
1290  }
1291  good++;
1292  max = TMath::Max(max, gap);
1293  gap = 0;
1294  }
1295  Bool_t thisGood = true; // Local result, assume success
1296  Int_t min = (r == 'I' ? minInner : minOuter);
1297  Double_t rate = len > 1 ? Double_t(good)/(len-1) : 0;
1298  if (rate < minRate) thisGood = false;
1299  if (len < min) thisGood = false;
1300  if (max > maxGap) thisGood = false;
1301  if (!thisGood) isGood = false; // Global result
1302  if (verbose)
1303  Printf("%s: %2d/%2d=%5.1f%%%-2s%5.1f%% good (%d) fits (%-2s %2d) "
1304  "max gap %2d (%-2s %2d) -> %s",
1305  ringArray->GetName(), good, len, 100*rate,
1306  (rate < minRate ? "<" : ">="), 100*minRate,
1307  minQuality,
1308  (len < min ? "<" : ">="), min, max,
1309  (max > maxGap ? ">" : "<="), maxGap,
1310  (thisGood ? "good" : "bad"));
1311 
1312  }
1313  SetBit(kIsGoodAsserted,true);
1314  SetBit(kIsGood, isGood);
1315 
1316  return isGood;
1317 }
1318 
1319 
1320 //____________________________________________________________________
1321 namespace {
1322  TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
1323  Int_t color)
1324  {
1325  TH1D* h = new TH1D(Form("%s_%s", name, title),
1326  Form("%s %s", name, title),
1327  axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
1328  h->SetDirectory(0);
1329  h->SetMarkerStyle(20);
1330  h->SetMarkerColor(color);
1331  h->SetMarkerSize(0.5);
1332  h->SetFillColor(color);
1333  h->SetFillStyle(3001);
1334  h->SetLineColor(color);
1335  return h;
1336  }
1337 }
1338 
1339 
1340 
1341 #define IDX2RING(I) (i == 0 || i == 1 || i == 3 ? 'I' : 'O')
1342 #define IDX2DET(I) (i == 0 ? 1 : (i == 1 || i == 2 ? 2 : 3))
1343 //____________________________________________________________________
1344 TList*
1346  Bool_t rel,
1347  Bool_t good,
1348  UShort_t maxN) const
1349 {
1350  // Get a list of THStacks
1351  Int_t nRings = fRings.GetEntriesFast();
1352  // Int_t nPad = 6+maxN-1; // 7 regular params, and maxN-1 weights
1353 
1354  enum {
1355  kChi2nu = 0,
1356  kC = 1,
1357  kDelta = 2,
1358  kXi = 3,
1359  kSigma = 4,
1360  kN = 5
1361  };
1362 
1363  THStack* sChi2nu;
1364  THStack* sC;
1365  THStack* sDelta;
1366  THStack* sXi;
1367  THStack* sSigma;
1368  // THStack* sigman;
1369  THStack* n;
1370  TList* stacks = new TList;
1371  stacks->AddAt(sChi2nu= new THStack("chi2", "#chi^{2}/#nu"), kChi2nu);
1372  stacks->AddAt(sC = new THStack("c", "C"), kC);
1373  stacks->AddAt(sDelta = new THStack("delta", "#Delta_{mp}"), kDelta);
1374  stacks->AddAt(sXi = new THStack("xi", "#xi"), kXi);
1375  stacks->AddAt(sSigma = new THStack("sigma", "#sigma"), kSigma);
1376  //stacks->AddAt(sigman= new THStack("sigman", "#sigma_{n}"), 5);
1377  stacks->AddAt(n = new THStack("n", "N"), kN);
1378  if (rel) {
1379  sChi2nu->SetName("qual");
1380  sChi2nu->SetTitle("Quality");
1381  n->SetName("good");
1382  n->SetTitle("Bin map");
1383  }
1384  for (Int_t i = 1; i <= maxN; i++) {
1385  stacks->AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), kN+i);
1386  }
1387 
1388  // TArrayD min(nPad);
1389  // TArrayD max(nPad);
1390  // min.Reset(100000);
1391  // max.Reset(-100000);
1392 
1393  for (Int_t i = 0; i < nRings; i++) {
1394  if (!fRings.At(i)) continue;
1395  TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1396  Int_t nFits = a->GetEntriesFast();
1397  UShort_t d = IDX2DET(i);
1398  Char_t r = IDX2RING(i);
1400 
1401  TH1* hChi = MakeHist(fEtaAxis,a->GetName(), "chi2", color);
1402  TH1* hC = MakeHist(fEtaAxis,a->GetName(), "c", color);
1403  TH1* hDelta = MakeHist(fEtaAxis,a->GetName(), "delta", color);
1404  TH1* hXi = MakeHist(fEtaAxis,a->GetName(), "xi", color);
1405  TH1* hSigma = MakeHist(fEtaAxis,a->GetName(), "sigma", color);
1406  // TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
1407  TH1* hN = MakeHist(fEtaAxis,a->GetName(), "n", color);
1408  TH1* hA[maxN];
1409  for (Int_t j = 0; j < maxN; j++) hA[j] = 0;
1410  const char* ho = (rel || !err ? "hist" : "e");
1411  sChi2nu->Add(hChi, "hist"); // 0
1412  sC ->Add(hC, ho); // 1
1413  sDelta ->Add(hDelta, ho); // 2
1414  sXi ->Add(hXi, ho); // 3
1415  sSigma ->Add(hSigma, ho); // 4
1416  // sigman->Add(hSigmaN,ho); // 5
1417  n ->Add(hN, "hist"); // 5
1418  hChi->SetFillColor(color);
1419  hChi->SetFillStyle(3001);
1420  hN->SetFillColor(color);
1421  hN->SetFillStyle(3001);
1422 
1423  for (Int_t k = 1; k <= maxN; k++) {
1424  hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
1425  static_cast<THStack*>(stacks->At(kN+k))->Add(hA[k-1], ho);
1426  }
1427 
1428  if (good) {
1429  for (Int_t j = 1; j <= fEtaAxis.GetNbins(); j++) {
1430  ELossFit* f = FindFit(d, r, j, 20);
1431  if (!f) continue;
1432 
1433  UpdateStackHist(f, rel, j, hChi, hN, hC, hDelta, hXi, hSigma, maxN, hA);
1434  }
1435  }
1436  else {
1437  for (Int_t j = 0; j < nFits; j++) {
1438  ELossFit* f = static_cast<ELossFit*>(a->At(j));
1439  if (!f) continue;
1440 
1441  UpdateStackHist(f, rel, CACHE(j,i,fEtaAxis.GetNbins()),
1442  hChi, hN, hC, hDelta, hXi, hSigma, maxN, hA);
1443  }
1444  }
1445  }
1446  return stacks;
1447 }
1448 //____________________________________________________________________
1449 void
1451  Int_t used,
1452  TH1* hChi, TH1* hN,
1453  TH1* hC, TH1* hDelta,
1454  TH1* hXi, TH1* hSigma,
1455  Int_t maxN, TH1** hA) const
1456 {
1457  Int_t b =f->fBin;
1458  Double_t chi2nu =(rel ? f->fQuality : (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
1459  Double_t c =(rel ? (f->fC >0 ?f->fEC /f->fC :0) : f->fC);
1460  Double_t delta =(rel ? (f->fDelta >0 ?f->fEDelta /f->fDelta :0) : f->fDelta);
1461  Double_t xi =(rel ? (f->fXi >0 ?f->fEXi /f->fXi :0) : f->fXi);
1462  Double_t sigma =(rel ? (f->fSigma >0 ?f->fESigma /f->fSigma :0) : f->fSigma);
1463  Int_t w =(rel ? used : f->FindMaxWeight());
1464  // Double_t sigman = (rel ? (f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0)
1465  // : f->SigmaN);
1466  hChi ->SetBinContent(b, chi2nu);
1467  hN ->SetBinContent(b, w);
1468  hC ->SetBinContent(b, c);
1469  hDelta ->SetBinContent(b, delta);
1470  hXi ->SetBinContent(b, xi);
1471  hSigma ->SetBinContent(b, sigma);
1472 
1473  if (!rel) {
1474  hC ->SetBinError(b, f->fEC);
1475  hDelta ->SetBinError(b, f->fEDelta);
1476  hXi ->SetBinError(b, f->fEXi);
1477  hSigma ->SetBinError(b, f->fESigma);
1478  // hSigmaN->SetBinError(b, f->fESigmaN);
1479  }
1480  for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
1481  Double_t a = (rel ? (f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0) : f->fA[k]);
1482  hA[k]->SetBinContent(b, a);
1483  if (!rel) hA[k]->SetBinError(b, f->fEA[k]);
1484  }
1485 }
1486 
1487 //____________________________________________________________________
1488 void
1490 {
1491  //
1492  // Draw this object
1493  //
1494  // Parameters:
1495  // option Options. Possible values are
1496  // - err Plot error bars
1497  //
1498  TString opt(Form("nostack %s", option));
1499  opt.ToLower();
1500  Bool_t rel = (opt.Contains("relative"));
1501  Bool_t err = (opt.Contains("error"));
1502  Bool_t clr = (opt.Contains("clear"));
1503  Bool_t gdd = (opt.Contains("good"));
1504  if (rel) opt.ReplaceAll("relative","");
1505  if (err) opt.ReplaceAll("error","");
1506  if (clr) opt.ReplaceAll("clear", "");
1507  if (gdd) opt.ReplaceAll("good", "");
1508 
1509  UShort_t maxN = 0;
1510  Int_t nRings = fRings.GetEntriesFast();
1511  for (Int_t i = 0; i < nRings; i++) {
1512  if (!fRings.At(i)) continue;
1513  TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1514  Int_t nFits = a->GetEntriesFast();
1515 
1516  for (Int_t j = 0; j < nFits; j++) {
1517  ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1518  if (!fit) continue;
1519  maxN = TMath::Max(maxN, UShort_t(fit->fN));
1520  }
1521  }
1522  // AliInfo(Form("Maximum N is %d", maxN));
1523  Int_t nPad = 6+maxN-1; // 7 regular params, and maxN-1 weights
1524  TVirtualPad* pad = gPad;
1525  if (clr) {
1526  pad->Clear();
1527  pad->SetTopMargin(0.02);
1528  pad->SetRightMargin(0.02);
1529  pad->SetBottomMargin(0.15);
1530  pad->SetLeftMargin(0.10);
1531  }
1532  pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
1533 
1534  TList* stacks = GetStacks(err, rel, gdd, maxN);
1535 
1536  Int_t nPad2 = (nPad+1) / 2;
1537  for (Int_t i = 0; i < nPad; i++) {
1538  Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
1539  TVirtualPad* p = pad->cd(iPad);
1540  p->SetLeftMargin(.15);
1541  p->SetFillColor(0);
1542  p->SetFillStyle(0);
1543  p->SetGridx();
1544  p->SetGridy();
1545  if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
1546 
1547 
1548  THStack* stack = static_cast<THStack*>(stacks->At(i));
1549  if (!stack->GetHists() || stack->GetHists()->GetEntries() <= 0) {
1550  AliWarningF("No histograms in %s", stack->GetName());
1551  continue;
1552  }
1553  // Double_t powMax = TMath::Log10(max[i]);
1554  // Double_t powMin = min[i] <= 0 ? powMax : TMath::Log10(min[i]);
1555  // if (powMax-powMin > 2. && min[i] != 0) p->SetLogy();
1556 
1557  // stack->SetMinimum(min[i]);
1558  // stack->SetMaximum(max[i]);
1559  stack->Draw(opt.Data());
1560 
1561  TString tit(stack->GetTitle());
1562  if (rel && i != 0 && i != 5)
1563  tit = Form("#delta %s/%s", tit.Data(), tit.Data());
1564  TH1* hist = stack->GetHistogram();
1565  TAxis* yaxis = hist->GetYaxis();
1566  yaxis->SetTitle(tit.Data());
1567  yaxis->SetTitleSize(0.15);
1568  yaxis->SetLabelSize(0.08);
1569  yaxis->SetTitleOffset(0.35);
1570  yaxis->SetTitleFont(132);
1571  yaxis->SetLabelFont(132);
1572  yaxis->SetNdivisions(5);
1573 
1574 
1575  TAxis* xaxis = stack->GetHistogram()->GetXaxis();
1576  xaxis->SetTitle("#eta");
1577  xaxis->SetTitleSize(0.15);
1578  xaxis->SetLabelSize(0.08);
1579  xaxis->SetTitleOffset(0.35);
1580  xaxis->SetTitleFont(132);
1581  xaxis->SetLabelFont(132);
1582  xaxis->SetNdivisions(10);
1583 
1584  stack->Draw(opt.Data());
1585  }
1586  pad->cd();
1587 }
1588 
1589 //____________________________________________________________________
1590 void
1592 {
1593  //
1594  // Print this object.
1595  //
1596  // Parameters:
1597  // option Options
1598  // - R Print recursive
1599  //
1600  //
1601  TString opt(option);
1602  opt.ToUpper();
1603  Int_t nRings = fRings.GetEntriesFast();
1604  bool recurse = opt.Contains("R");
1605  bool cache = opt.Contains("C") && fCache.GetSize() > 0;
1606  Int_t nBins = fEtaAxis.GetNbins();
1607 
1608  std::cout << "Low cut in fit range: " << fLowCut << "\n"
1609  << "Eta axis: " << nBins
1610  << " bins, range [" << fEtaAxis.GetXmin() << ","
1611  << fEtaAxis.GetXmax() << "]" << std::endl;
1612 
1613  for (Int_t i = 0; i < nRings; i++) {
1614  if (!fRings.At(i)) continue;
1615  TObjArray* a = static_cast<TObjArray*>(fRings.At(i));
1616  Int_t nFits = a->GetEntriesFast();
1617 
1618  std::cout << a->GetName() << " [" << nFits << " entries]"
1619  << (recurse ? ":\n" : "\t");
1620  Int_t min = fEtaAxis.GetNbins()+1;
1621  Int_t max = 0;
1622  for (Int_t j = 0; j < nFits; j++) {
1623  if (!a->At(j)) continue;
1624 
1625  min = TMath::Min(j, min);
1626  max = TMath::Max(j, max);
1627 
1628  if (recurse) {
1629  std::cout << "Bin # " << j << "\t";
1630  ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1631  fit->Print(option);
1632  }
1633  }
1634  if (!recurse)
1635  std::cout << " bin range: " << std::setw(3) << min
1636  << "-" << std::setw(3) << max << " " << std::setw(3)
1637  << (max-min+1) << " bins" << std::endl;
1638  }
1639 
1640  if (!cache) return;
1641 
1642  std::cout << " eta bin | Fit bin \n"
1643  << " # range | FMD1i FMD2i FMD2o FMD3i FMD3o"
1644  // << "----+-----+++------+-----------------------------------"
1645  << std::endl;
1646  size_t oldPrec = std::cout.precision();
1647  std::cout.precision(3);
1648  for (Int_t i = 1; i <= nBins; i++) {
1649  std::cout << std::setw(4) << i << " "
1650  << std::setw(5) << std::showpos << fEtaAxis.GetBinLowEdge(i)
1651  << " - " << std::setw(5) << fEtaAxis.GetBinUpEdge(i)
1652  << std::noshowpos << " | ";
1653  for (Int_t j = 0; j < 5; j++) {
1654  Int_t bin = CACHE(i,j,nBins);
1655  if (bin <= 0) std::cout << " ";
1656  else std::cout << std::setw(5) << bin
1657  << (bin == i ? ' ' : '*') << ' ';
1658  }
1659  std::cout << std::endl;
1660  }
1661  std::cout.precision(oldPrec);
1662 }
1663 
1664 //____________________________________________________________________
1665 void
1667 {
1668  //
1669  // Browse this object
1670  //
1671  // Parameters:
1672  // b
1673  //
1674  b->Add(&fRings);
1675  b->Add(&fEtaAxis);
1676 }
1677 
1678 
1679 
1680 //____________________________________________________________________
1681 //
1682 // EOF
1683 //
Int_t color[]
print message on plot with ok/not ok
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:27
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)
Double_t GetXiCut(UShort_t d, Char_t r, Int_t etaBin, Double_t f) const
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
AliStack * stack
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
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.
Double_t GetAvgXiSigmaCut(UShort_t d, Char_t r, Int_t etaBin, Double_t f) const
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)
Double_t GetXiSigmaCut(UShort_t d, Char_t r, Int_t etaBin, Double_t f) const
Definition: External.C:212
Double_t GetProbabilityCut(UShort_t d, Char_t r, Int_t etaBin, Double_t f) const
Double_t GetXiSigmaCut(Double_t f) const
void SetEtaAxis(const TAxis &axis)
Double_t GetMpvCut(Double_t f) const
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
Double_t GetAvgXiSigmaCut(Double_t f) 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
Double_t GetMpvCut(UShort_t d, Char_t r, Int_t etaBin, Double_t f) const
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
Double_t GetXiCut(Double_t f) const
Definition: External.C:196
Int_t FindMaxWeight(Double_t maxRelError=2 *fgMaxRelError, Double_t leastWeight=fgLeastWeight, UShort_t maxN=999) const