PWGLF/FORWARD/analysis2/AliFMDCorrELossFit.cxx

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