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