Added ROOT THtml docs to .cxx and top .h
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDCorrELossFit.cxx
1 // Object holding the Energy loss fit 'correction'
2 // 
3 // These are generated from Monte-Carlo or real ESDs. 
4 #include "AliFMDCorrELossFit.h"
5 #include "AliForwardUtil.h"
6 #include <TF1.h>
7 #include <TBrowser.h>
8 #include <TVirtualPad.h>
9 #include <THStack.h>
10 #include <TH1D.h>
11 #include <AliLog.h>
12 #include <TList.h>
13 #include <iostream>
14 #include <iomanip>
15
16 Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .2;
17 Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
18 Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu   = 20;
19
20 //____________________________________________________________________
21 AliFMDCorrELossFit::ELossFit::ELossFit()
22   : fN(0),
23     fNu(0),
24     fChi2(0),
25     fC(0),
26     fDelta(0),
27     fXi(0),
28     fSigma(0),
29     fSigmaN(0),
30     fA(0),
31     fEC(0),
32     fEDelta(0),
33     fEXi(0),
34     fESigma(0),
35     fESigmaN(0),
36     fEA(0),
37     fQuality(0), 
38     fDet(0), 
39     fRing('\0'),
40     fBin(0)
41 {
42   //
43   // Default constructor 
44   // 
45   //
46 }
47 //____________________________________________________________________
48 AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
49   : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ? 
50        f.GetParameter(AliForwardUtil::ELossFitter::kN) : 
51        1),
52     fNu(f.GetNDF()),
53     fChi2(f.GetChisquare()),
54     fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
55     fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
56     fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
57     fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
58     fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
59     fA(0),
60     fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
61     fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
62     fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
63     fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
64     fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
65     fEA(0),
66     fQuality(quality),
67     fDet(0), 
68     fRing('\0'),
69     fBin(0)
70 {
71   // 
72   // Construct from a function
73   // 
74   // Parameters:
75   //    quality Quality flag
76   //    f       Function
77   //
78   if (fN <= 0) return;
79   fA  = new Double_t[fN];
80   fEA = new Double_t[fN];
81   for (Int_t i = 0; i < fN-1; i++) { 
82     fA[i]  = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
83     fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
84   }
85   fA[fN-1]  = -9999;
86   fEA[fN-1] = -9999;
87 }
88
89 //____________________________________________________________________
90 AliFMDCorrELossFit::ELossFit::ELossFit(Int_t     quality,UShort_t  n, 
91                                        Double_t  chi2,   UShort_t  nu, 
92                                        Double_t  c,      Double_t  ec, 
93                                        Double_t  delta,  Double_t  edelta, 
94                                        Double_t  xi,     Double_t  exi,
95                                        Double_t  sigma,  Double_t  esigma, 
96                                        Double_t  sigman, Double_t  esigman, 
97                                        Double_t* a,      Double_t* ea)
98   : fN(n),
99     fNu(nu),
100     fChi2(chi2),
101     fC(c),
102     fDelta(delta),
103     fXi(xi),
104     fSigma(sigma),
105     fSigmaN(sigman),
106     fA(0),
107     fEC(ec),
108     fEDelta(edelta),
109     fEXi(exi),
110     fESigma(esigma),
111     fESigmaN(esigman),
112     fEA(0),
113     fQuality(quality),
114     fDet(0), 
115     fRing('\0'),
116     fBin(0)
117 {
118   // 
119   // Constructor with full parameter set
120   // 
121   // Parameters:
122   //    quality   Quality flag
123   //    n         @f$ N@f$ - Number of fitted peaks
124   //    chi2      @f$ \chi^2 @f$
125   //    nu        @f$ \nu @f$ - number degrees of freedom
126   //    c         @f$ C@f$ - scale constant
127   //    ec        @f$ \delta C@f$ - error on @f$ C@f$ 
128   //    delta     @f$ \Delta@f$ - Most probable value             
129   //    edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$ 
130   //    xi        @f$ \xi@f$ - width  
131   //    exi       @f$ \delta\xi@f$ - error on @f$\xi@f$ 
132   //    sigma     @f$ \sigma@f$ - Width of Gaussian                
133   //    esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$ 
134   //    sigman    @f$ \sigma_n@f$ - Noise width           
135   //    esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$ 
136   //    a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for 
137   //                  @f$ i=2,\ldots@f$ 
138   //    ea        Array of @f$ N-1@f$ error on the weights @f$ a_i@f$ for 
139   //                  @f$ i=2,\ldots@f$ 
140   //
141   if (fN <= 0) return;
142   fA  = new Double_t[fN];
143   fEA = new Double_t[fN];
144   for (Int_t i = 0; i < fN-1; i++) { 
145     fA[i]  = a[i];
146     fEA[i] = ea[i];
147   }
148   fA[fN-1]  = -9999;
149   fEA[fN-1] = -9999;
150 }
151 //____________________________________________________________________
152 AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
153   : TObject(o), 
154     fN(o.fN),
155     fNu(o.fNu),
156     fChi2(o.fChi2),
157     fC(o.fC),
158     fDelta(o.fDelta),
159     fXi(o.fXi),
160     fSigma(o.fSigma),
161     fSigmaN(o.fSigmaN),
162     fA(0),
163     fEC(o.fEC),
164     fEDelta(o.fEDelta),
165     fEXi(o.fEXi),
166     fESigma(o.fESigma),
167     fESigmaN(o.fESigmaN),
168     fEA(0),
169     fQuality(o.fQuality),
170     fDet(o.fDet), 
171     fRing(o.fRing),
172     fBin(o.fBin)
173 {
174   // 
175   // Copy constructor 
176   // 
177   // Parameters:
178   //    o Object to copy from 
179   //
180   if (fN <= 0) return;
181   fA  = new Double_t[fN];
182   fEA = new Double_t[fN];
183   for (Int_t i = 0; i < fN-1; i++) { 
184     fA[i]  = o.fA[i];
185     fEA[i] = o.fEA[i];
186   }
187   fA[fN-1]  = -9999;
188   fEA[fN-1] = -9999;
189 }
190
191 //____________________________________________________________________
192 AliFMDCorrELossFit::ELossFit&
193 AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
194 {
195   // 
196   // Assignment operator 
197   // 
198   // Parameters:
199   //    o Object to assign from 
200   // 
201   // Return:
202   //    Reference to this object 
203   //
204   fN       = o.fN;
205   fNu      = o.fNu;
206   fChi2    = o.fChi2;
207   fC       = o.fC;
208   fDelta   = o.fDelta;
209   fXi      = o.fXi;
210   fSigma   = o.fSigma;
211   fSigmaN  = o.fSigmaN;
212   fEC      = o.fEC;
213   fEDelta  = o.fEDelta;
214   fEXi     = o.fEXi;
215   fESigma  = o.fESigma;
216   fESigmaN = o.fESigmaN;
217   fQuality = o.fQuality;
218   fDet     = o.fDet; 
219   fRing    = o.fRing;
220   fBin     = o.fBin;
221   if (fA)  delete [] fA;
222   if (fEA) delete [] fEA; 
223   fA  = 0;
224   fEA = 0;
225
226   if (fN <= 0) return *this;
227   fA  = new Double_t[fN];
228   fEA = new Double_t[fN];
229   for (Int_t i = 0; i < fN; i++) { 
230     fA[i]  = o.fA[i];
231     fEA[i] = o.fEA[i];
232   }
233
234   return *this;
235 }
236
237 //____________________________________________________________________
238 AliFMDCorrELossFit::ELossFit::~ELossFit()
239 {
240   if (fA)  delete[] fA;
241   if (fEA) delete[] fEA;
242 }
243
244
245 //____________________________________________________________________
246 Int_t
247 AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError, 
248                                             Double_t leastWeight, 
249                                             UShort_t  maxN) const
250 {
251   // 
252   // Find the maximum weight to use.  The maximum weight is the
253   // largest i for which 
254   // 
255   // - @f$ i \leq \max{N}@f$ 
256   // - @f$ a_i > \min{a}@f$ 
257   // - @f$ \delta a_i/a_i > \delta_{max}@f$ 
258   // 
259   // Parameters:
260   //    maxRelError @f$ \min{a}@f$ 
261   //    leastWeight @f$ \delta_{max}@f$ 
262   //    maxN        @f$ \max{N}@f$      
263   // 
264   // Return:
265   //    The largest index @f$ i@f$ for which the above
266   // conditions hold.  Will never return less than 1. 
267   //
268   Int_t n = TMath::Min(maxN, UShort_t(fN-1));
269   Int_t m = 1;
270   // fN is one larger than we have data 
271   for (Int_t i = 0; i < n-1; i++, m++) { 
272     if (fA[i] < leastWeight)  break;
273     if (fEA[i] / fA[i] > maxRelError) break;
274   }
275   return m;
276 }
277
278 //____________________________________________________________________
279 Double_t 
280 AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x, 
281                                        UShort_t maxN) const
282 {
283   // 
284   // Evaluate 
285   // @f[ 
286   //  f_N(x;\Delta,\xi,\sigma') = 
287   //     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
288   // @f] 
289   //
290   // (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
291   // that fulfills the requirements 
292   // 
293   // Parameters:
294   //    x           Where to evaluate 
295   //    maxN      @f$ \max{N}@f$    
296   // 
297   // Return:
298   //    @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
299   //
300   return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN, 
301                                      TMath::Min(maxN, UShort_t(fN)), fA);
302 }
303
304 //____________________________________________________________________
305 Double_t 
306 AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x, 
307                                                UShort_t maxN) const
308 {                                                                       
309   // 
310   // Evaluate 
311   // @f[ 
312   //   f_W(x;\Delta,\xi,\sigma') = 
313   //   \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
314   //     f_N(x;\Delta,\xi,\sigma')} = 
315   //   \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
316   //     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
317   // @f] 
318   // where @f$ n@f$ fulfills the requirements (see FindMaxWeight). 
319   //
320   // If the denominator is zero, then 1 is returned. 
321   //
322   // See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
323   // for more information on the evaluated functions. 
324   // 
325   // Parameters:
326   //    x           Where to evaluate 
327   //    maxN      @f$ \max{N}@f$      
328   // 
329   // Return:
330   //    @f$ f_W(x;\Delta,\xi,\sigma')@f$.  
331   //
332   UShort_t n   = TMath::Min(maxN, UShort_t(fN-1));
333   Double_t num = 0;
334   Double_t den = 0;
335   for (Int_t i = 1; i <= n; i++) {
336     Double_t a = (i == 1 ? 1 : fA[i-1]);
337     if (fA[i-1] < 0) break;
338     Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
339     num += i * a * f;
340     den += a * f;
341   }
342   if (den <= 0) return 1;
343   return num / den;
344 }
345
346
347 #define OUTPAR(N,V,E)                   \
348   std::setprecision(9) <<               \
349   std::setw(12) << N << ": " <<         \
350   std::setw(14) << V << " +/- " <<      \
351   std::setw(14) << E << " (" <<         \
352   std::setprecision(-1) <<              \
353   std::setw(5)  << 100*(V>0?E/V:1) << "%)\n"
354   
355
356 //____________________________________________________________________
357 Int_t
358 AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
359 {
360   // 
361   // Compare to another ELossFit object. 
362   // 
363   // - +1, if this quality is better (larger) than other objects quality
364   // - -1, if this quality is worse (smaller) than other objects quality
365   // - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
366   // - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
367   // - 0 otherwise 
368   // 
369   // Parameters:
370   //    o Other object to compare to 
371   //
372   const ELossFit* other = static_cast<const ELossFit*>(o);
373   if (this->fQuality < other->fQuality) return -1;
374   if (this->fQuality > other->fQuality) return +1;
375   Double_t chi2nu  = (fNu == 0 ? 1000 : fChi2 / fNu);
376   Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
377   if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
378   if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
379   return 0;
380 }
381
382 //____________________________________________________________________
383 void
384 AliFMDCorrELossFit::ELossFit::Print(Option_t*) const
385 {
386   // 
387   // Information to standard output 
388   // 
389   // Parameters:
390   //    option Not used 
391   //
392   std::cout << GetName() << ":\n"
393             << " chi^2/nu = " << fChi2 << "/" << fNu << " = " 
394             << (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
395             << "     Quality:   " << fQuality << "\n" 
396             << "  NParticles:   " << fN << "  (" << FindMaxWeight() << ")\n"
397             << OUTPAR("Delta", fDelta, fEDelta) 
398             << OUTPAR("xi", fXi, fEXi)
399             << OUTPAR("sigma", fSigma, fESigma)
400             << OUTPAR("sigma_n", fSigmaN, fESigmaN);
401   for (Int_t i = 0; i < fN-1; i++) 
402     std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
403   std::cout << std::flush;
404 }
405
406 //____________________________________________________________________
407 const Char_t*
408 AliFMDCorrELossFit::ELossFit::GetName() const 
409 {
410   // 
411   // Get the name of this object 
412   // 
413   // 
414   // Return:
415   //    
416   //
417   return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
418 }
419
420 //____________________________________________________________________
421 void
422 AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
423 {
424   // 
425   // Browse this object 
426   // 
427   // Parameters:
428   //    b Browser
429   //
430   Draw(b ? b->GetDrawOption() : "comp");
431   gPad->SetLogy();
432   gPad->Update();
433 }
434      
435 //____________________________________________________________________
436 void
437 AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
438 {
439   // 
440   // Draw this fit 
441   // 
442   // Parameters:
443   //    option Options 
444   //  - COMP  Draw components too 
445   //
446   TString opt(option);
447   opt.ToUpper();
448   bool comp = false;
449   if (opt.Contains("COMP")) { 
450     opt.ReplaceAll("COMP","");
451     comp = true;
452   }
453   if (!opt.Contains("SAME")) { 
454     gPad->Clear();
455   }
456
457   TObjArray cleanup;
458   TF1* tot = AliForwardUtil::MakeNLandauGaus(1, 
459                                              fDelta, fXi, 
460                                              fSigma, fSigmaN, 
461                                              fN,     fA, 
462                                              0.01,   10);
463   tot->SetLineColor(kBlack);
464   tot->SetLineWidth(2);
465   tot->SetLineStyle(1);
466   tot->SetTitle(GetName());
467   Double_t max = tot->GetMaximum();
468
469   if (!opt.Contains("SAME")) {
470     TH1* frame = new TH1F(GetName(), 
471                           Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
472                           100, 0, 10);
473     frame->SetMinimum(max/10000);
474     frame->SetMaximum(max*1.4);
475     frame->SetXTitle("#Delta / #Delta_{mip}");
476     frame->Draw();
477     opt.Append(" SAME");
478   }
479   tot->DrawCopy(opt.Data());
480   cleanup.Add(tot);
481
482   if (!comp) { 
483     gPad->cd();
484     return;
485   }
486
487   Double_t min = max;
488   opt.Append(" same");
489   Int_t maxW = FindMaxWeight();
490   for (Int_t i=1; i <= fN; i++) { 
491     TF1* f = AliForwardUtil::MakeILandauGaus((i == 1 ? 1 : fA[i-2]), 
492                                              fDelta, fXi, 
493                                              fSigma, fSigmaN, 
494                                              i,      0.01, 10);
495     f->SetLineWidth(2);
496     f->SetLineStyle(i > maxW ? 2 : 1);
497     min = TMath::Min(f->GetMaximum(), min);
498     f->DrawCopy(opt.Data());
499     cleanup.Add(f);
500   }
501   min /= 100;
502   tot->GetHistogram()->SetMaximum(max);
503   tot->GetHistogram()->SetMinimum(min);
504   tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
505
506   gPad->cd();
507 }
508
509                                                  
510 //____________________________________________________________________
511 #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
512
513 //____________________________________________________________________
514 void 
515 AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu, 
516                                                Double_t maxRelError, 
517                                                Double_t leastWeight)
518 {
519   // 
520   // Calculate the quality 
521   //
522   Int_t qual = 0;
523   if (fNu > 0 && fChi2 / fNu < maxChi2nu) qual += 4;;
524   if (CHECKPAR(fDelta,  fEDelta,  maxRelError)) qual++;
525   if (CHECKPAR(fXi,     fEXi,     maxRelError)) qual++;
526   if (CHECKPAR(fSigma,  fESigma,  maxRelError)) qual++;
527   if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
528   qual += FindMaxWeight(1.5*maxRelError, leastWeight, fN);
529   fQuality = qual;
530 }
531
532 //____________________________________________________________________
533 AliFMDCorrELossFit::AliFMDCorrELossFit()
534   : TObject(), 
535     fRings(), 
536     fEtaAxis(0,0,0), 
537     fLowCut(0)
538 {
539   // 
540   // Default constructor 
541   //
542   fRings.SetOwner(kTRUE);
543   fEtaAxis.SetTitle("#eta");
544   fEtaAxis.SetName("etaAxis");
545   fRings.SetName("rings");
546 }
547
548 //____________________________________________________________________
549 AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
550   : TObject(o), 
551     fRings(o.fRings),
552     fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()), 
553     fLowCut(0)
554 {
555   // 
556   // Copy constructor 
557   // 
558   // Parameters:
559   //    o Object to copy from 
560   //
561   fEtaAxis.SetTitle("#eta");
562   fEtaAxis.SetName("etaAxis");
563 }
564
565 //____________________________________________________________________
566 AliFMDCorrELossFit::~AliFMDCorrELossFit()
567 {
568   // 
569   // Destructor 
570   //
571   fRings.Clear();
572 }
573
574 //____________________________________________________________________
575 AliFMDCorrELossFit&
576 AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
577 {
578   // 
579   // Assignment operator 
580   // 
581   // Parameters:
582   //    o Object to assign from 
583   // 
584   // Return:
585   //    Reference to this object 
586   //
587   fRings = o.fRings;
588   fLowCut = o.fLowCut;
589   SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
590
591   return *this;
592 }
593 //____________________________________________________________________
594 Int_t
595 AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
596 {
597   // 
598   // Find the eta bin corresponding to the given eta 
599   // 
600   // Parameters:
601   //    eta  Eta value 
602   // 
603   // Return:
604   //    Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
605   // eta, or 0 if out of range.
606   //
607   if (fEtaAxis.GetXmin() == fEtaAxis.GetXmax() || fEtaAxis.GetNbins() == 0) {
608     AliWarning("No eta axis defined");
609     return -1;
610   }
611   Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
612   if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
613   return bin;
614 }
615
616 //____________________________________________________________________
617 Bool_t
618 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
619 {
620   // 
621   // Set the fit parameters from a function 
622   // 
623   // Parameters:
624   //    d       Detector
625   //    r       Ring 
626   //    etaBin  Eta (bin number, 1->nBins)
627   //    f       ELoss fit result - note, the object will take ownership
628   //  
629   TObjArray* ringArray = GetOrMakeRingArray(d, r);
630   if (!ringArray) { 
631     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
632     return kFALSE;
633   }
634   if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) { 
635     AliError(Form("bin=%d is out of range [%d,%d]", 
636                   etaBin, 1, fEtaAxis.GetNbins()));
637     return kFALSE;
638   }
639   // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
640   ringArray->AddAtAndExpand(fit, etaBin);
641   return kTRUE;
642 }
643
644 //____________________________________________________________________
645 Bool_t
646 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
647 {
648   // 
649   // Set the fit parameters from a function 
650   // 
651   // Parameters:
652   //    d    Detector
653   //    r    Ring 
654   //    eta  Eta 
655   //    f    ELoss fit result - note, the object will take ownership
656   //  
657   Int_t bin = FindEtaBin(eta);
658   if (bin <= 0) { 
659     AliError(Form("eta=%f is out of range [%f,%f]", 
660                   eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
661     return kFALSE;
662   }
663
664   return SetFit(d, r, bin, fit);
665 }
666 //____________________________________________________________________
667 Bool_t
668 AliFMDCorrELossFit::SetFit(UShort_t  d,      Char_t    r, 
669                            Double_t  eta, 
670                            Int_t     quality,UShort_t  n, 
671                            Double_t  chi2,   UShort_t  nu, 
672                            Double_t  c,      Double_t  ec, 
673                            Double_t  delta,  Double_t  edelta, 
674                            Double_t  xi,     Double_t  exi,
675                            Double_t  sigma,  Double_t  esigma, 
676                            Double_t  sigman, Double_t  esigman, 
677                            Double_t* a,      Double_t* ea)
678 {
679   // 
680   // Set the fit parameters from a function 
681   // 
682   // Parameters:
683   //    d         Detector number
684   //    r         Ring identifier 
685   //    eta       Eta value
686   //    quality   Quality flag
687   //    n         @f$ N@f$ - Number of fitted peaks
688   //    chi2      @f$ \chi^2 @f$
689   //    nu        @f$ \nu @f$ - number degrees of freedom
690   //    c         @f$ C@f$ - scale constant
691   //    ec        @f$ \delta C@f$ - error on @f$ C@f$ 
692   //    delta     @f$ \Delta@f$ - most probable value
693   //    edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$ 
694   //    xi        @f$ \xi@f$ - Landau width               
695   //    exi       @f$ \delta\xi@f$ - error on @f$\xi@f$ 
696   //    sigma     @f$ \sigma@f$ - Gaussian width
697   //    esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$ 
698   //    sigman    @f$ \sigma_n@f$ - Noise width           
699   //    esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$ 
700   //    a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for 
701   //                  @f$ i=2,\ldots@f$ 
702   //    ea        Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for 
703   //                  @f$ i=2,\ldots@f$ 
704   //
705   ELossFit* e = new ELossFit(quality, n, 
706                              chi2,    nu,
707                              c,       ec,
708                              delta,   edelta,
709                              xi,      exi,
710                              sigma,   esigma,
711                              sigman,  esigman,
712                              a,       ea);
713   if (!SetFit(d, r, eta, e)) { 
714     delete e;
715     return kFALSE;
716   }
717   return kTRUE;
718 }
719 //____________________________________________________________________
720 Bool_t
721 AliFMDCorrELossFit::SetFit(UShort_t  d, Char_t r, Double_t eta, 
722                            Int_t quality, const TF1& f)
723 {
724   // 
725   // Set the fit parameters from a function 
726   // 
727   // Parameters:
728   //    d        Detector
729   //    r        Ring 
730   //    eta      Eta 
731   //    quality  Quality flag
732   //    f        Function from fit 
733   //  
734   ELossFit* e = new ELossFit(quality, f);
735   if (!SetFit(d, r, eta, e)) { 
736     delete e;
737     return kFALSE;
738   }
739   return kTRUE;
740 }
741 //____________________________________________________________________
742 AliFMDCorrELossFit::ELossFit*
743 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Int_t etabin) const
744 {
745   // 
746   // Find the fit corresponding to the specified parameters 
747   // 
748   // Parameters:
749   //    d      Detector 
750   //    r      Ring 
751   //    etabin Eta bin (1 based)
752   // 
753   // Return:
754   //    Fit parameters or null in case of problems 
755   //
756   TObjArray* ringArray = GetRingArray(d, r);
757   if (!ringArray) { 
758     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
759     return 0;
760   }
761   if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) { 
762     AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
763                   etabin, 1, fEtaAxis.GetNbins(), d, r));
764     return 0;
765   }
766   if (etabin > ringArray->GetEntriesFast()) { 
767     AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c", 
768                   etabin, 1, ringArray->GetEntriesFast(), d, r));
769     return 0;
770   }
771   else if (etabin >= ringArray->GetEntriesFast()) { 
772     AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, " 
773                     "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
774                     etabin-1));
775     etabin--;
776   }
777   else if (!ringArray->At(etabin)) { 
778     AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d", 
779                     etabin, d, r, etabin+1));
780     etabin++;
781   }
782   return static_cast<ELossFit*>(ringArray->At(etabin));
783 }
784 //____________________________________________________________________
785 AliFMDCorrELossFit::ELossFit*
786 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Double_t eta) const
787 {
788   // 
789   // Find the fit corresponding to the specified parameters 
790   // 
791   // Parameters:
792   //    d   Detector 
793   //    r   Ring 
794   //    eta Eta value 
795   // 
796   // Return:
797   //    Fit parameters or null in case of problems 
798   //
799   Int_t etabin = FindEtaBin(eta);
800   return FindFit(d, r, etabin);
801 }
802 //____________________________________________________________________
803 TObjArray*
804 AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
805 {
806   // 
807   // Get the ring array corresponding to the specified ring
808   // 
809   // Parameters:
810   //    d Detector 
811   //    r Ring 
812   // 
813   // Return:
814   //    Pointer to ring array, or null in case of problems
815   //
816   Int_t idx = -1;
817   switch (d) { 
818   case 1:   idx = 0; break;
819   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
820   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
821   }
822   if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
823   return static_cast<TObjArray*>(fRings.At(idx));
824 }
825 //____________________________________________________________________
826 TObjArray*
827 AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
828 {
829   // 
830   // Get the ring array corresponding to the specified ring
831   // 
832   // Parameters:
833   //    d Detector 
834   //    r Ring 
835   // 
836   // Return:
837   //    Pointer to ring array, or newly created container 
838   //
839   Int_t idx = -1;
840   switch (d) { 
841   case 1:   idx = 0; break;
842   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
843   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
844   }
845   if (idx < 0) return 0;
846   if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
847     TObjArray* a = new TObjArray(0);
848     // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
849     a->SetName(Form("FMD%d%c", d, r));
850     a->SetOwner();
851     fRings.AddAtAndExpand(a, idx);
852   }
853   return static_cast<TObjArray*>(fRings.At(idx));
854 }
855
856 namespace { 
857   TH1D* MakeHist(const TAxis& axis, const char* name, const char* title, 
858                  Int_t color)
859   {
860     TH1D* h = new TH1D(Form("%s_%s", name, title), 
861                        Form("%s %s", name, title), 
862                        axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
863     h->SetDirectory(0);
864     h->SetMarkerStyle(20);
865     h->SetMarkerColor(color);
866     h->SetMarkerSize(0.5);
867     h->SetFillColor(color);
868     h->SetFillStyle(3001);
869     h->SetLineColor(color);
870     return h;
871   }
872 }
873
874 //____________________________________________________________________
875 void
876 AliFMDCorrELossFit::Draw(Option_t* option)
877 {
878   // 
879   // Draw this object 
880   // 
881   // Parameters:
882   //    option Options.  Possible values are 
883   //  - err Plot error bars 
884   //
885   TString opt(Form("nostack %s", option));
886   opt.ToLower();
887   Bool_t  rel = (opt.Contains("rel"));
888   Bool_t  err = (opt.Contains("err"));
889   if (rel) opt.ReplaceAll("rel","");
890   if (err) opt.ReplaceAll("err","");
891   Int_t    nRings = fRings.GetEntriesFast();
892   UShort_t maxN   = 0;
893   for (Int_t i = 0; i < nRings; i++) { 
894     if (!fRings.At(i)) continue;
895     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
896     Int_t      nFits = a->GetEntriesFast();
897
898     for (Int_t j = 0; j < nFits; j++) {
899       ELossFit* fit = static_cast<ELossFit*>(a->At(j));
900       if (!fit) continue;
901       maxN          = TMath::Max(maxN, UShort_t(fit->fN));
902     }
903   }
904   // AliInfo(Form("Maximum N is %d", maxN));
905   Int_t nPad = 7+maxN-1; // 7 regular params, and maxN-1 weights
906   TVirtualPad* pad = gPad;
907   pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
908   
909   THStack* chi2nu;
910   THStack* c;
911   THStack* delta;
912   THStack* xi;
913   THStack* sigma;
914   THStack* sigman;
915   THStack* n;
916   TList stacks;
917   stacks.AddAt(chi2nu= new THStack("chi2",   "#chi^{2}/#nu"), 0);
918   stacks.AddAt(c     = new THStack("c",       "C"),           1);
919   stacks.AddAt(delta = new THStack("delta",  "#Delta_{mp}"),  2);
920   stacks.AddAt(xi    = new THStack("xi",     "#xi"),          3);
921   stacks.AddAt(sigma = new THStack("sigma",  "#sigma"),       4);
922   stacks.AddAt(sigman= new THStack("sigman", "#sigma_{n}"),   5);
923   stacks.AddAt(n     = new THStack("n",      "N"),            6);
924   for (Int_t i = 1; i <= maxN; i++) {
925     stacks.AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), 6+i);
926   }
927
928   for (Int_t i = 0; i < nRings; i++) { 
929     if (!fRings.At(i)) continue;
930     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
931     Int_t      nFits = a->GetEntriesFast();
932     Int_t      color = 0;
933     switch (i) { 
934     case 0: color = kRed+2;   break;
935     case 1: color = kGreen+2; break;
936     case 2: color = kGreen-2; break;
937     case 3: color = kBlue+2;  break;
938     case 4: color = kBlue-2;  break;
939     }
940
941     TH1D* hChi    = MakeHist(fEtaAxis,a->GetName(), "chi2",   color);
942     TH1D* hC      = MakeHist(fEtaAxis,a->GetName(), "c",      color);
943     TH1D* hDelta  = MakeHist(fEtaAxis,a->GetName(), "delta",  color);
944     TH1D* hXi     = MakeHist(fEtaAxis,a->GetName(), "xi",     color);
945     TH1D* hSigma  = MakeHist(fEtaAxis,a->GetName(), "sigma",  color);
946     TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
947     TH1D* hN      = MakeHist(fEtaAxis,a->GetName(), "n",      color);
948     TH1D* hA[maxN];
949     const char* ho = (rel || !err ? "hist" : "e");
950     chi2nu->Add(hChi,   "hist"); // 0
951     c     ->Add(hC,     ho);     // 1
952     delta ->Add(hDelta, ho);     // 2
953     xi    ->Add(hXi,    ho);     // 3
954     sigma ->Add(hSigma, ho);     // 4
955     sigman->Add(hSigmaN,ho);     // 5
956     n     ->Add(hN,     "hist"); // 6
957     hChi->SetFillColor(color);
958     hChi->SetFillStyle(3001);
959     hN->SetFillColor(color);
960     hN->SetFillStyle(3001);
961
962     for (Int_t k = 1; k <= maxN; k++) { 
963       hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
964       static_cast<THStack*>(stacks.At(6+k))->Add(hA[k-1], ho);
965     }
966                           
967     for (Int_t j = 0; j < nFits; j++) {
968       ELossFit* f = static_cast<ELossFit*>(a->At(j));
969       if (!f) continue;
970
971       Int_t     b = f->fBin;
972       Int_t     nW = f->FindMaxWeight();
973       hChi   ->SetBinContent(b, (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
974       hN     ->SetBinContent(b, nW);
975
976       if (rel) { 
977         hC     ->SetBinContent(b, f->fC     >0 ?f->fEC     /f->fC      :0);
978         hDelta ->SetBinContent(b, f->fDelta >0 ?f->fEDelta /f->fDelta  :0);
979         hXi    ->SetBinContent(b, f->fXi    >0 ?f->fEXi    /f->fXi     :0);
980         hSigma ->SetBinContent(b, f->fSigma >0 ?f->fESigma /f->fSigma  :0);
981         hSigmaN->SetBinContent(b, f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0);
982       }
983       else {
984         hC     ->SetBinContent(b, f->fC);
985         hDelta ->SetBinContent(b, f->fDelta);
986         hXi    ->SetBinContent(b, f->fXi);
987         hSigma ->SetBinContent(b, f->fSigma);
988         hSigmaN->SetBinContent(b, f->fSigmaN);
989         hC     ->SetBinError(b, f->fEC);
990         hDelta ->SetBinError(b, f->fEDelta);
991         hXi    ->SetBinError(b, f->fEXi);
992         hSigma ->SetBinError(b, f->fESigma);
993         hSigmaN->SetBinError(b, f->fESigmaN);
994       }
995       for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) { 
996         if (rel) 
997           hA[k]->SetBinContent(b, f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0);
998         else {
999           hA[k]->SetBinContent(b, f->fA[k]);
1000           hA[k]->SetBinError(b, f->fEA[k]);
1001         }
1002       }
1003     }
1004   }
1005   Int_t nPad2 = (nPad+1) / 2;
1006   for (Int_t i = 0; i < nPad; i++) {
1007     Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
1008     TVirtualPad* p = pad->cd(iPad);
1009     p->SetLeftMargin(.15);
1010     p->SetFillColor(0);
1011     p->SetFillStyle(0);
1012     p->SetGridx();
1013     p->SetGridy();
1014     if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
1015
1016     THStack* stack = static_cast<THStack*>(stacks.At(i));
1017     // AliInfo(Form("Drawing %s (%d) in pad %d", stack->GetName(), i, iPad));
1018     stack->Draw(opt.Data());
1019
1020     TString tit(stack->GetTitle());
1021     if (rel && i != 0 && i != 6)
1022       tit = Form("#delta %s/%s", tit.Data(), tit.Data());
1023     TH1*   hist  = stack->GetHistogram();
1024     TAxis* yaxis = hist->GetYaxis();
1025     yaxis->SetTitle(tit.Data());
1026     yaxis->SetTitleSize(0.15);
1027     yaxis->SetLabelSize(0.08);
1028     yaxis->SetTitleOffset(0.35);
1029     yaxis->SetNdivisions(5);
1030
1031     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
1032     xaxis->SetTitle("#eta");
1033     xaxis->SetTitleSize(0.15);
1034     xaxis->SetLabelSize(0.08);
1035     xaxis->SetTitleOffset(0.35);
1036     xaxis->SetNdivisions(10);
1037
1038
1039     if (!rel) {
1040       switch (i) { 
1041       case 0:                         break; // chi^2/nu 
1042       case 1:                         break; // C 
1043       case 2: stack->SetMinimum(0.4); break; // Delta 
1044       case 3: stack->SetMinimum(0.02);break; // xi
1045       case 4: stack->SetMinimum(0.05);break; // sigma
1046       case 5:                         break; // sigmaN
1047       case 6: 
1048         stack->SetMinimum(-.1); 
1049         stack->SetMaximum(maxN+.5); 
1050         break;                              // N
1051       default:                       break; // Weights
1052       }
1053     }
1054     stack->DrawClone(opt.Data());
1055   }
1056   pad->cd();      
1057 }
1058
1059 //____________________________________________________________________
1060 void
1061 AliFMDCorrELossFit::Print(Option_t* option) const
1062 {
1063   TString opt(option);
1064   opt.ToUpper();
1065   Int_t nRings = fRings.GetEntriesFast();
1066   bool recurse = opt.Contains("R");
1067
1068   std::cout << "Low cut in fit range: " << fLowCut << "\n"
1069             << "Eta axis:             " << fEtaAxis.GetNbins() 
1070             << " bins, range [" << fEtaAxis.GetXmin() << "," 
1071             << fEtaAxis.GetXmax() << "]" << std::endl;
1072   
1073   for (Int_t i = 0; i < nRings; i++) { 
1074     if (!fRings.At(i)) continue;
1075     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
1076     Int_t      nFits = a->GetEntriesFast();
1077
1078     std::cout << a->GetName() << " [" << nFits << " entries]" 
1079               << (recurse ? ":\n" : "\t");
1080     Int_t min = fEtaAxis.GetNbins()+1;
1081     Int_t max = 0;
1082     for (Int_t j = 0; j < nFits; j++) {
1083       if (!a->At(j)) continue;
1084       
1085       min = TMath::Min(j, min);
1086       max = TMath::Max(j, max);
1087
1088       if (recurse) {
1089         std::cout << "Bin # " << j << "\t";
1090         ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1091         fit->Print(option);
1092       }
1093     }
1094     if (!recurse) 
1095       std::cout << " bin range: " << std::setw(3) << min 
1096                 << "-" << std::setw(3) << max << " " << std::setw(3) 
1097                 << (max-min+1) << " bins" << std::endl;
1098   }
1099 }
1100
1101 //____________________________________________________________________
1102 void
1103 AliFMDCorrELossFit::Browse(TBrowser* b)
1104 {
1105   // 
1106   // Browse this object 
1107   // 
1108   // Parameters:
1109   //    b 
1110   //
1111   b->Add(&fRings);
1112   b->Add(&fEtaAxis);
1113 }
1114
1115
1116
1117 //____________________________________________________________________
1118 //
1119 // EOF
1120 //