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