More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDCorrELossFit.cxx
1 #include "AliFMDCorrELossFit.h"
2 #include "AliForwardUtil.h"
3 #include <TF1.h>
4 #include <TBrowser.h>
5 #include <TVirtualPad.h>
6 #include <THStack.h>
7 #include <TH1D.h>
8 #include <AliLog.h>
9 #include <TList.h>
10 #include <iostream>
11 #include <iomanip>
12
13 Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .2;
14 Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
15 Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu   = 20;
16
17 //____________________________________________________________________
18 AliFMDCorrELossFit::ELossFit::ELossFit()
19   : fN(0),
20     fNu(0),
21     fChi2(0),
22     fC(0),
23     fDelta(0),
24     fXi(0),
25     fSigma(0),
26     fSigmaN(0),
27     fA(0),
28     fEC(0),
29     fEDelta(0),
30     fEXi(0),
31     fESigma(0),
32     fESigmaN(0),
33     fEA(0),
34     fQuality(0), 
35     fDet(0), 
36     fRing('\0'),
37     fBin(0)
38 {}
39 //____________________________________________________________________
40 AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
41   : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ? 
42        f.GetParameter(AliForwardUtil::ELossFitter::kN) : 
43        1),
44     fNu(f.GetNDF()),
45     fChi2(f.GetChisquare()),
46     fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
47     fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
48     fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
49     fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
50     fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
51     fA(0),
52     fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
53     fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
54     fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
55     fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
56     fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
57     fEA(0),
58     fQuality(quality),
59     fDet(0), 
60     fRing('\0'),
61     fBin(0)
62 {
63   if (fN <= 0) return;
64   fA  = new Double_t[fN];
65   fEA = new Double_t[fN];
66   for (Int_t i = 0; i < fN-1; i++) { 
67     fA[i]  = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
68     fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
69   }
70   fA[fN-1]  = -9999;
71   fEA[fN-1] = -9999;
72 }
73
74 //____________________________________________________________________
75 AliFMDCorrELossFit::ELossFit::ELossFit(Int_t     quality,UShort_t  n, 
76                                        Double_t  chi2,   UShort_t  nu, 
77                                        Double_t  c,      Double_t  ec, 
78                                        Double_t  delta,  Double_t  edelta, 
79                                        Double_t  xi,     Double_t  exi,
80                                        Double_t  sigma,  Double_t  esigma, 
81                                        Double_t  sigman, Double_t  esigman, 
82                                        Double_t* a,      Double_t* ea)
83   : fN(n),
84     fNu(nu),
85     fChi2(chi2),
86     fC(c),
87     fDelta(delta),
88     fXi(xi),
89     fSigma(sigma),
90     fSigmaN(sigman),
91     fA(0),
92     fEC(ec),
93     fEDelta(edelta),
94     fEXi(exi),
95     fESigma(esigma),
96     fESigmaN(esigman),
97     fEA(0),
98     fQuality(quality),
99     fDet(0), 
100     fRing('\0'),
101     fBin(0)
102 {
103   if (fN <= 0) return;
104   fA  = new Double_t[fN];
105   fEA = new Double_t[fN];
106   for (Int_t i = 0; i < fN-1; i++) { 
107     fA[i]  = a[i];
108     fEA[i] = ea[i];
109   }
110   fA[fN-1]  = -9999;
111   fEA[fN-1] = -9999;
112 }
113 //____________________________________________________________________
114 AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
115   : TObject(o), 
116     fN(o.fN),
117     fNu(o.fNu),
118     fChi2(o.fChi2),
119     fC(o.fC),
120     fDelta(o.fDelta),
121     fXi(o.fXi),
122     fSigma(o.fSigma),
123     fSigmaN(o.fSigmaN),
124     fA(0),
125     fEC(o.fEC),
126     fEDelta(o.fEDelta),
127     fEXi(o.fEXi),
128     fESigma(o.fESigma),
129     fESigmaN(o.fESigmaN),
130     fEA(0),
131     fQuality(o.fQuality),
132     fDet(o.fDet), 
133     fRing(o.fRing),
134     fBin(o.fBin)
135 {
136   if (fN <= 0) return;
137   fA  = new Double_t[fN];
138   fEA = new Double_t[fN];
139   for (Int_t i = 0; i < fN-1; i++) { 
140     fA[i]  = o.fA[i];
141     fEA[i] = o.fEA[i];
142   }
143   fA[fN-1]  = -9999;
144   fEA[fN-1] = -9999;
145 }
146
147 //____________________________________________________________________
148 AliFMDCorrELossFit::ELossFit&
149 AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
150 {
151   fN       = o.fN;
152   fNu      = o.fNu;
153   fChi2    = o.fChi2;
154   fC       = o.fC;
155   fDelta   = o.fDelta;
156   fXi      = o.fXi;
157   fSigma   = o.fSigma;
158   fSigmaN  = o.fSigmaN;
159   fEC      = o.fEC;
160   fEDelta  = o.fEDelta;
161   fEXi     = o.fEXi;
162   fESigma  = o.fESigma;
163   fESigmaN = o.fESigmaN;
164   fQuality = o.fQuality;
165   fDet     = o.fDet; 
166   fRing    = o.fRing;
167   fBin     = o.fBin;
168   if (fA)  delete [] fA;
169   if (fEA) delete [] fEA; 
170   fA  = 0;
171   fEA = 0;
172
173   if (fN <= 0) return *this;
174   fA  = new Double_t[fN];
175   fEA = new Double_t[fN];
176   for (Int_t i = 0; i < fN; i++) { 
177     fA[i]  = o.fA[i];
178     fEA[i] = o.fEA[i];
179   }
180
181   return *this;
182 }
183
184 //____________________________________________________________________
185 AliFMDCorrELossFit::ELossFit::~ELossFit()
186 {
187   if (fA)  delete[] fA;
188   if (fEA) delete[] fEA;
189 }
190
191
192 //____________________________________________________________________
193 Int_t
194 AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError, 
195                                             Double_t leastWeight, 
196                                             UShort_t  maxN) const
197 {
198   Int_t n = TMath::Min(maxN, UShort_t(fN-1));
199   Int_t m = 1;
200   // fN is one larger than we have data 
201   for (Int_t i = 0; i < n-1; i++, m++) { 
202     if (fA[i] < leastWeight)  break;
203     if (fEA[i] / fA[i] > maxRelError) break;
204   }
205   return m;
206 }
207
208 //____________________________________________________________________
209 Double_t 
210 AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x, 
211                                        UShort_t maxN) const
212 {
213   return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN, 
214                                      TMath::Min(maxN, UShort_t(fN)), fA);
215 }
216
217 //____________________________________________________________________
218 Double_t 
219 AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x, 
220                                                UShort_t maxN) const
221 {                                                                       
222   UShort_t n   = TMath::Min(maxN, UShort_t(fN-1));
223   Double_t num = 0;
224   Double_t den = 0;
225   for (Int_t i = 1; i <= n; i++) {
226     Double_t a = (i == 1 ? 1 : fA[i-1]);
227     if (fA[i-1] < 0) break;
228     Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
229     num += i * a * f;
230     den += a * f;
231   }
232   if (den <= 0) return 1;
233   return num / den;
234 }
235
236
237 #define OUTPAR(N,V,E)                   \
238   std::setprecision(9) <<               \
239   std::setw(12) << N << ": " <<         \
240   std::setw(14) << V << " +/- " <<      \
241   std::setw(14) << E << " (" <<         \
242   std::setprecision(-1) <<              \
243   std::setw(5)  << 100*(V>0?E/V:1) << "%)\n"
244   
245
246 //____________________________________________________________________
247 Int_t
248 AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
249 {
250   const ELossFit* other = static_cast<const ELossFit*>(o);
251   if (this->fQuality < other->fQuality) return -1;
252   if (this->fQuality > other->fQuality) return +1;
253   Double_t chi2nu  = (fNu == 0 ? 1000 : fChi2 / fNu);
254   Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
255   if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
256   if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
257   return 0;
258 }
259
260 //____________________________________________________________________
261 void
262 AliFMDCorrELossFit::ELossFit::Print(Option_t*) const
263 {
264   std::cout << GetName() << ":\n"
265             << " chi^2/nu = " << fChi2 << "/" << fNu << " = " 
266             << (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
267             << "     Quality:   " << fQuality << "\n" 
268             << "  NParticles:   " << fN << "  (" << FindMaxWeight() << ")\n"
269             << OUTPAR("Delta", fDelta, fEDelta) 
270             << OUTPAR("xi", fXi, fEXi)
271             << OUTPAR("sigma", fSigma, fESigma)
272             << OUTPAR("sigma_n", fSigmaN, fESigmaN);
273   for (Int_t i = 0; i < fN-1; i++) 
274     std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
275   std::cout << std::flush;
276 }
277
278 //____________________________________________________________________
279 const Char_t*
280 AliFMDCorrELossFit::ELossFit::GetName() const 
281 {
282   return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
283 }
284
285 //____________________________________________________________________
286 void
287 AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
288 {
289   // Draw this one 
290   Draw(b ? b->GetDrawOption() : "comp");
291   gPad->SetLogy();
292   gPad->Update();
293 }
294      
295 //____________________________________________________________________
296 void
297 AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
298 {
299   TString opt(option);
300   opt.ToUpper();
301   bool comp = false;
302   if (opt.Contains("COMP")) { 
303     opt.ReplaceAll("COMP","");
304     comp = true;
305   }
306   if (!opt.Contains("SAME")) { 
307     gPad->Clear();
308   }
309
310   TObjArray cleanup;
311   TF1* tot = AliForwardUtil::MakeNLandauGaus(1, 
312                                              fDelta, fXi, 
313                                              fSigma, fSigmaN, 
314                                              fN,     fA, 
315                                              0.01,   10);
316   tot->SetLineColor(kBlack);
317   tot->SetLineWidth(2);
318   tot->SetLineStyle(1);
319   tot->SetTitle(GetName());
320   Double_t max = tot->GetMaximum();
321
322   if (!opt.Contains("SAME")) {
323     TH1* frame = new TH1F(GetName(), 
324                           Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
325                           100, 0, 10);
326     frame->SetMinimum(max/10000);
327     frame->SetMaximum(max*1.4);
328     frame->SetXTitle("#Delta / #Delta_{mip}");
329     frame->Draw();
330     opt.Append(" SAME");
331   }
332   tot->DrawCopy(opt.Data());
333   cleanup.Add(tot);
334
335   if (!comp) { 
336     gPad->cd();
337     return;
338   }
339
340   Double_t min = max;
341   opt.Append(" same");
342   Int_t maxW = FindMaxWeight();
343   for (Int_t i=1; i <= fN; i++) { 
344     TF1* f = AliForwardUtil::MakeILandauGaus((i == 1 ? 1 : fA[i-2]), 
345                                              fDelta, fXi, 
346                                              fSigma, fSigmaN, 
347                                              i,      0.01, 10);
348     f->SetLineWidth(2);
349     f->SetLineStyle(i > maxW ? 2 : 1);
350     min = TMath::Min(f->GetMaximum(), min);
351     f->DrawCopy(opt.Data());
352     cleanup.Add(f);
353   }
354   min /= 100;
355   tot->GetHistogram()->SetMaximum(max);
356   tot->GetHistogram()->SetMinimum(min);
357   tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
358
359   gPad->cd();
360 }
361
362                                                  
363 //____________________________________________________________________
364 #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
365
366 //____________________________________________________________________
367 void 
368 AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu, 
369                                                Double_t maxRelError, 
370                                                Double_t leastWeight)
371 {
372   Int_t qual = 0;
373   if (fNu > 0 && fChi2 / fNu < maxChi2nu) qual += 4;;
374   if (CHECKPAR(fDelta,  fEDelta,  maxRelError)) qual++;
375   if (CHECKPAR(fXi,     fEXi,     maxRelError)) qual++;
376   if (CHECKPAR(fSigma,  fESigma,  maxRelError)) qual++;
377   if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
378   qual += FindMaxWeight(1.5*maxRelError, leastWeight, fN);
379   fQuality = qual;
380 }
381
382 //____________________________________________________________________
383 AliFMDCorrELossFit::AliFMDCorrELossFit()
384   : TObject(), 
385     fRings(), 
386     fEtaAxis(0,0,0), 
387     fLowCut(0)
388 {
389   fRings.SetOwner(kTRUE);
390   fEtaAxis.SetTitle("#eta");
391   fEtaAxis.SetName("etaAxis");
392   fRings.SetName("rings");
393 }
394
395 //____________________________________________________________________
396 AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
397   : TObject(o), 
398     fRings(o.fRings),
399     fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()), 
400     fLowCut(0)
401 {
402   fEtaAxis.SetTitle("#eta");
403   fEtaAxis.SetName("etaAxis");
404 }
405
406 //____________________________________________________________________
407 AliFMDCorrELossFit::~AliFMDCorrELossFit()
408 {
409   fRings.Clear();
410 }
411
412 //____________________________________________________________________
413 AliFMDCorrELossFit&
414 AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
415 {
416   fRings = o.fRings;
417   fLowCut = o.fLowCut;
418   SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
419
420   return *this;
421 }
422 //____________________________________________________________________
423 Int_t
424 AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
425 {
426   if (fEtaAxis.GetXmin() == fEtaAxis.GetXmax() || fEtaAxis.GetNbins() == 0) {
427     AliWarning("No eta axis defined");
428     return -1;
429   }
430   Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
431   if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
432   return bin;
433 }
434
435 //____________________________________________________________________
436 Bool_t
437 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
438 {
439   TObjArray* ringArray = GetOrMakeRingArray(d, r);
440   if (!ringArray) { 
441     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
442     return kFALSE;
443   }
444   if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) { 
445     AliError(Form("bin=%d is out of range [%d,%d]", 
446                   etaBin, 1, fEtaAxis.GetNbins()));
447     return kFALSE;
448   }
449   // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
450   ringArray->AddAtAndExpand(fit, etaBin);
451   return kTRUE;
452 }
453
454 //____________________________________________________________________
455 Bool_t
456 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
457 {
458   Int_t bin = FindEtaBin(eta);
459   if (bin <= 0) { 
460     AliError(Form("eta=%f is out of range [%f,%f]", 
461                   eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
462     return kFALSE;
463   }
464
465   return SetFit(d, r, bin, fit);
466 }
467 //____________________________________________________________________
468 Bool_t
469 AliFMDCorrELossFit::SetFit(UShort_t  d,      Char_t    r, 
470                            Double_t  eta, 
471                            Int_t     quality,UShort_t  n, 
472                            Double_t  chi2,   UShort_t  nu, 
473                            Double_t  c,      Double_t  ec, 
474                            Double_t  delta,  Double_t  edelta, 
475                            Double_t  xi,     Double_t  exi,
476                            Double_t  sigma,  Double_t  esigma, 
477                            Double_t  sigman, Double_t  esigman, 
478                            Double_t* a,      Double_t* ea)
479 {
480   ELossFit* e = new ELossFit(quality, n, 
481                              chi2,    nu,
482                              c,       ec,
483                              delta,   edelta,
484                              xi,      exi,
485                              sigma,   esigma,
486                              sigman,  esigman,
487                              a,       ea);
488   if (!SetFit(d, r, eta, e)) { 
489     delete e;
490     return kFALSE;
491   }
492   return kTRUE;
493 }
494 //____________________________________________________________________
495 Bool_t
496 AliFMDCorrELossFit::SetFit(UShort_t  d, Char_t r, Double_t eta, 
497                            Int_t quality, const TF1& f)
498 {
499   ELossFit* e = new ELossFit(quality, f);
500   if (!SetFit(d, r, eta, e)) { 
501     delete e;
502     return kFALSE;
503   }
504   return kTRUE;
505 }
506 //____________________________________________________________________
507 AliFMDCorrELossFit::ELossFit*
508 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Int_t etabin) const
509 {
510   TObjArray* ringArray = GetRingArray(d, r);
511   if (!ringArray) { 
512     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
513     return 0;
514   }
515   if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) { 
516     AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
517                   etabin, 1, fEtaAxis.GetNbins(), d, r));
518     return 0;
519   }
520   if (etabin > ringArray->GetEntriesFast()) { 
521     AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
522                   etabin, 1, ringArray->GetEntriesFast(), d, r));
523     return 0;
524   }
525   else if (etabin >= ringArray->GetEntriesFast()) { 
526     AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, " 
527                     "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
528                     etabin-1));
529     etabin--;
530   }
531   else if (!ringArray->At(etabin)) { 
532     AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d", 
533                     etabin, d, r, etabin+1));
534     etabin++;
535   }
536   return static_cast<ELossFit*>(ringArray->At(etabin));
537 }
538 //____________________________________________________________________
539 AliFMDCorrELossFit::ELossFit*
540 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Double_t eta) const
541 {
542   Int_t etabin = FindEtaBin(eta);
543   return FindFit(d, r, etabin);
544 }
545 //____________________________________________________________________
546 TObjArray*
547 AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
548 {
549   Int_t idx = -1;
550   switch (d) { 
551   case 1:   idx = 0; break;
552   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
553   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
554   }
555   if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
556   return static_cast<TObjArray*>(fRings.At(idx));
557 }
558 //____________________________________________________________________
559 TObjArray*
560 AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
561 {
562   Int_t idx = -1;
563   switch (d) { 
564   case 1:   idx = 0; break;
565   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
566   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
567   }
568   if (idx < 0) return 0;
569   if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
570     TObjArray* a = new TObjArray(0);
571     // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
572     a->SetName(Form("FMD%d%c", d, r));
573     a->SetOwner();
574     fRings.AddAtAndExpand(a, idx);
575   }
576   return static_cast<TObjArray*>(fRings.At(idx));
577 }
578
579 namespace { 
580   TH1D* MakeHist(const TAxis& axis, const char* name, const char* title, 
581                  Int_t color)
582   {
583     TH1D* h = new TH1D(Form("%s_%s", name, title), 
584                        Form("%s %s", name, title), 
585                        axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
586     h->SetDirectory(0);
587     h->SetMarkerStyle(20);
588     h->SetMarkerColor(color);
589     h->SetMarkerSize(0.5);
590     h->SetFillColor(color);
591     h->SetFillStyle(3001);
592     h->SetLineColor(color);
593     return h;
594   }
595 }
596
597 //____________________________________________________________________
598 void
599 AliFMDCorrELossFit::Draw(Option_t* option)
600 {
601   TString opt(Form("nostack %s", option));
602   opt.ToLower();
603   Bool_t  rel = (opt.Contains("rel"));
604   Bool_t  err = (opt.Contains("err"));
605   if (rel) opt.ReplaceAll("rel","");
606   if (err) opt.ReplaceAll("err","");
607   Int_t    nRings = fRings.GetEntriesFast();
608   UShort_t maxN   = 0;
609   for (Int_t i = 0; i < nRings; i++) { 
610     if (!fRings.At(i)) continue;
611     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
612     Int_t      nFits = a->GetEntriesFast();
613
614     for (Int_t j = 0; j < nFits; j++) {
615       ELossFit* fit = static_cast<ELossFit*>(a->At(j));
616       if (!fit) continue;
617       maxN          = TMath::Max(maxN, UShort_t(fit->fN));
618     }
619   }
620   // AliInfo(Form("Maximum N is %d", maxN));
621   Int_t nPad = 7+maxN-1; // 7 regular params, and maxN-1 weights
622   TVirtualPad* pad = gPad;
623   pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
624   
625   THStack* chi2nu;
626   THStack* c;
627   THStack* delta;
628   THStack* xi;
629   THStack* sigma;
630   THStack* sigman;
631   THStack* n;
632   TList stacks;
633   stacks.AddAt(chi2nu= new THStack("chi2",   "#chi^{2}/#nu"), 0);
634   stacks.AddAt(c     = new THStack("c",       "C"),           1);
635   stacks.AddAt(delta = new THStack("delta",  "#Delta_{mp}"),  2);
636   stacks.AddAt(xi    = new THStack("xi",     "#xi"),          3);
637   stacks.AddAt(sigma = new THStack("sigma",  "#sigma"),       4);
638   stacks.AddAt(sigman= new THStack("sigman", "#sigma_{n}"),   5);
639   stacks.AddAt(n     = new THStack("n",      "N"),            6);
640   for (Int_t i = 1; i <= maxN; i++) {
641     stacks.AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), 6+i);
642   }
643
644   for (Int_t i = 0; i < nRings; i++) { 
645     if (!fRings.At(i)) continue;
646     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
647     Int_t      nFits = a->GetEntriesFast();
648     Int_t      color = 0;
649     switch (i) { 
650     case 0: color = kRed+2;   break;
651     case 1: color = kGreen+2; break;
652     case 2: color = kGreen-2; break;
653     case 3: color = kBlue+2;  break;
654     case 4: color = kBlue-2;  break;
655     }
656
657     TH1D* hChi    = MakeHist(fEtaAxis,a->GetName(), "chi2",   color);
658     TH1D* hC      = MakeHist(fEtaAxis,a->GetName(), "c",      color);
659     TH1D* hDelta  = MakeHist(fEtaAxis,a->GetName(), "delta",  color);
660     TH1D* hXi     = MakeHist(fEtaAxis,a->GetName(), "xi",     color);
661     TH1D* hSigma  = MakeHist(fEtaAxis,a->GetName(), "sigma",  color);
662     TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
663     TH1D* hN      = MakeHist(fEtaAxis,a->GetName(), "n",      color);
664     TH1D* hA[maxN];
665     const char* ho = (rel || !err ? "hist" : "e");
666     chi2nu->Add(hChi,   "hist"); // 0
667     c     ->Add(hC,     ho);     // 1
668     delta ->Add(hDelta, ho);     // 2
669     xi    ->Add(hXi,    ho);     // 3
670     sigma ->Add(hSigma, ho);     // 4
671     sigman->Add(hSigmaN,ho);     // 5
672     n     ->Add(hN,     "hist"); // 6
673     hChi->SetFillColor(color);
674     hChi->SetFillStyle(3001);
675     hN->SetFillColor(color);
676     hN->SetFillStyle(3001);
677
678     for (Int_t k = 1; k <= maxN; k++) { 
679       hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
680       static_cast<THStack*>(stacks.At(6+k))->Add(hA[k-1], ho);
681     }
682                           
683     for (Int_t j = 0; j < nFits; j++) {
684       ELossFit* f = static_cast<ELossFit*>(a->At(j));
685       if (!f) continue;
686
687       Int_t     b = f->fBin;
688       Int_t     nW = f->FindMaxWeight();
689       hChi   ->SetBinContent(b, (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
690       hN     ->SetBinContent(b, nW);
691
692       if (rel) { 
693         hC     ->SetBinContent(b, f->fC     >0 ?f->fEC     /f->fC      :0);
694         hDelta ->SetBinContent(b, f->fDelta >0 ?f->fEDelta /f->fDelta  :0);
695         hXi    ->SetBinContent(b, f->fXi    >0 ?f->fEXi    /f->fXi     :0);
696         hSigma ->SetBinContent(b, f->fSigma >0 ?f->fESigma /f->fSigma  :0);
697         hSigmaN->SetBinContent(b, f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0);
698       }
699       else {
700         hC     ->SetBinContent(b, f->fC);
701         hDelta ->SetBinContent(b, f->fDelta);
702         hXi    ->SetBinContent(b, f->fXi);
703         hSigma ->SetBinContent(b, f->fSigma);
704         hSigmaN->SetBinContent(b, f->fSigmaN);
705         hC     ->SetBinError(b, f->fEC);
706         hDelta ->SetBinError(b, f->fEDelta);
707         hXi    ->SetBinError(b, f->fEXi);
708         hSigma ->SetBinError(b, f->fESigma);
709         hSigmaN->SetBinError(b, f->fESigmaN);
710       }
711       for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) { 
712         if (rel) 
713           hA[k]->SetBinContent(b, f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0);
714         else {
715           hA[k]->SetBinContent(b, f->fA[k]);
716           hA[k]->SetBinError(b, f->fEA[k]);
717         }
718       }
719     }
720   }
721   Int_t nPad2 = (nPad+1) / 2;
722   for (Int_t i = 0; i < nPad; i++) {
723     Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
724     TVirtualPad* p = pad->cd(iPad);
725     p->SetLeftMargin(.15);
726     p->SetFillColor(0);
727     p->SetFillStyle(0);
728     p->SetGridx();
729     p->SetGridy();
730     if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
731
732     THStack* stack = static_cast<THStack*>(stacks.At(i));
733     // AliInfo(Form("Drawing %s (%d) in pad %d", stack->GetName(), i, iPad));
734     stack->Draw(opt.Data());
735
736     TString tit(stack->GetTitle());
737     if (rel && i != 0 && i != 6)
738       tit = Form("#delta %s/%s", tit.Data(), tit.Data());
739     TH1*   hist  = stack->GetHistogram();
740     TAxis* yaxis = hist->GetYaxis();
741     yaxis->SetTitle(tit.Data());
742     yaxis->SetTitleSize(0.15);
743     yaxis->SetLabelSize(0.08);
744     yaxis->SetTitleOffset(0.35);
745     yaxis->SetNdivisions(5);
746
747     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
748     xaxis->SetTitle("#eta");
749     xaxis->SetTitleSize(0.15);
750     xaxis->SetLabelSize(0.08);
751     xaxis->SetTitleOffset(0.35);
752     xaxis->SetNdivisions(10);
753
754
755     if (!rel) {
756       switch (i) { 
757       case 0:                         break; // chi^2/nu 
758       case 1:                         break; // C 
759       case 2: stack->SetMinimum(0.4); break; // Delta 
760       case 3: stack->SetMinimum(0.02);break; // xi
761       case 4: stack->SetMinimum(0.05);break; // sigma
762       case 5:                         break; // sigmaN
763       case 6: 
764         stack->SetMinimum(-.1); 
765         stack->SetMaximum(maxN+.5); 
766         break;                              // N
767       default:                       break; // Weights
768       }
769     }
770     stack->DrawClone(opt.Data());
771   }
772   pad->cd();      
773 }
774
775 //____________________________________________________________________
776 void
777 AliFMDCorrELossFit::Print(Option_t* option) const
778 {
779   TString opt(option);
780   opt.ToUpper();
781   Int_t nRings = fRings.GetEntriesFast();
782   bool recurse = opt.Contains("R");
783
784   std::cout << "Low cut in fit range: " << fLowCut << "\n"
785             << "Eta axis:             " << fEtaAxis.GetNbins() 
786             << " bins, range [" << fEtaAxis.GetXmin() << "," 
787             << fEtaAxis.GetXmax() << "]" << std::endl;
788   
789   for (Int_t i = 0; i < nRings; i++) { 
790     if (!fRings.At(i)) continue;
791     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
792     Int_t      nFits = a->GetEntriesFast();
793
794     std::cout << a->GetName() << " [" << nFits << " entries]" 
795               << (recurse ? ":\n" : "\t");
796     Int_t min = fEtaAxis.GetNbins()+1;
797     Int_t max = 0;
798     for (Int_t j = 0; j < nFits; j++) {
799       if (!a->At(j)) continue;
800       
801       min = TMath::Min(j, min);
802       max = TMath::Max(j, max);
803
804       if (recurse) {
805         std::cout << "Bin # " << j << "\t";
806         ELossFit* fit = static_cast<ELossFit*>(a->At(j));
807         fit->Print(option);
808       }
809     }
810     if (!recurse) 
811       std::cout << " bin range: " << std::setw(3) << min 
812                 << "-" << std::setw(3) << max << " " << std::setw(3) 
813                 << (max-min+1) << " bins" << std::endl;
814   }
815 }
816
817 //____________________________________________________________________
818 void
819 AliFMDCorrELossFit::Browse(TBrowser* b)
820 {
821   b->Add(&fRings);
822   b->Add(&fEtaAxis);
823 }
824
825
826
827 //____________________________________________________________________
828 //
829 // EOF
830 //