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