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   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   tot->GetHistogram()->SetMaximum(max);
 592   tot->GetHistogram()->SetMinimum(min);
 593   tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
 594   if (l) l->Draw();
 595
 596   gPad->cd();
 597 }
 598
 599
 600 //____________________________________________________________________
 601 #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
 602
 603 //____________________________________________________________________
 604 Double_t
 605 AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f) const
 606 {
 607   //
 608   // Return
 609   //    Delta * f
 610   return f * fDelta;
 611 }
 612 //____________________________________________________________________
 613 Double_t
 614 AliFMDCorrELossFit::ELossFit::GetLowerBound(Double_t f,
 615                                             Bool_t includeSigma) const
 616 {
 617   //
 618   // Return
 619   //    Delta - f * (xi + sigma)
 620   return fDelta - f * (fXi + (includeSigma ? fSigma : 0));
 621 }
 622
 623 //____________________________________________________________________
 624 void
 625 AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
 626                                                Double_t maxRelError,
 627                                                Double_t leastWeight)
 628 {
 629   //
 630   // Calculate the quality
 631   //
 632   Double_t decline = maxChi2nu;
 633   Int_t qual = 0;
 634   if (fNu > 0) {
 635     Double_t red = fChi2 / fNu;
 636     if (red < maxChi2nu) qual += 4;
 637     else {
 638       Int_t q = (maxChi2nu+decline - red) / decline * 4;
 639       if (q > 0) qual += q;
 640     }
 641   }
 642   if (CHECKPAR(fDelta,  fEDelta,  maxRelError)) qual++;
 643   if (CHECKPAR(fXi,     fEXi,     maxRelError)) qual++;
 644   if (CHECKPAR(fSigma,  fESigma,  maxRelError)) qual++;
 645   if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
 646   qual += FindMaxWeight(2*maxRelError, leastWeight, fN);
 647   fQuality = qual;
 648 }
 649
 650 //____________________________________________________________________
 651 AliFMDCorrELossFit::AliFMDCorrELossFit()
 652   : TObject(),
 653     fRings(),
 654     fEtaAxis(0,0,0),
 655     fLowCut(0),
 656     fCache(0)
 657 {
 658   //
 659   // Default constructor
 660   //
 661   fRings.SetOwner(kTRUE);
 662   fEtaAxis.SetTitle("#eta");
 663   fEtaAxis.SetName("etaAxis");
 664   fRings.SetName("rings");
 665 }
 666
 667 //____________________________________________________________________
 668 AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
 669   : TObject(o),
 670     fRings(o.fRings),
 671     fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
 672     fLowCut(0),
 673     fCache(0)
 674 {
 675   //
 676   // Copy constructor
 677   //
 678   // Parameters:
 679   //    o Object to copy from
 680   //
 681   fEtaAxis.SetTitle("#eta");
 682   fEtaAxis.SetName("etaAxis");
 683 }
 684
 685 //____________________________________________________________________
 686 AliFMDCorrELossFit::~AliFMDCorrELossFit()
 687 {
 688   //
 689   // Destructor
 690   //
 691   fRings.Clear();
 692 }
 693
 694 //____________________________________________________________________
 695 AliFMDCorrELossFit&
 696 AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
 697 {
 698   //
 699   // Assignment operator
 700   //
 701   // Parameters:
 702   //    o Object to assign from
 703   //
 704   // Return:
 705   //    Reference to this object
 706   //
 707   if (&o == this) return *this;
 708   fRings = o.fRings;
 709   fLowCut = o.fLowCut;
 710   fCache  = o.fCache;
 711   SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
 712
 713   return *this;
 714 }
 715 #define CACHE(BIN,IDX,OFFSET) fCache[IDX*OFFSET+BIN-1]
 716
 717 //____________________________________________________________________
 718 void
 719 AliFMDCorrELossFit::CacheBins(UShort_t minQuality) const
 720 {
 721   if (fCache.GetSize() > 0) return;
 722
 723   Int_t nRings = fRings.GetEntriesFast();
 724   Int_t offset = fEtaAxis.GetNbins();
 725
 726   fCache.Set(nRings*offset);
 727   fCache.Reset(-1);
 728
 729   for (Int_t i = 0; i < nRings; i++) {
 730     TObjArray* ringArray  = static_cast<TObjArray*>(fRings.At(i));
 731
 732     // First loop to find where we actually have fits
 733     Int_t nFits      = 0;
 734     Int_t nGood      = 0;
 735     Int_t minBin     = offset+1;
 736     Int_t maxBin     = -1;
 737     Int_t realMinBin = offset+1;
 738     Int_t realMaxBin = -1;
 739     for (Int_t j = 1; j < ringArray->GetEntriesFast(); j++) {
 740       ELossFit* fit = static_cast<ELossFit*>(ringArray->At(j));
 741       if (!fit) continue;
 742       nFits++;
 743
 744       // Update our range off possible fits
 745       realMinBin = TMath::Min(j, realMinBin);
 746       realMaxBin = TMath::Max(j, realMaxBin);
 747
 748       // Check the quality of the fit
 749       fit->CalculateQuality(AliFMDCorrELossFit::ELossFit::fgMaxChi2nu,
 750                             AliFMDCorrELossFit::ELossFit::fgMaxRelError,
 751                             AliFMDCorrELossFit::ELossFit::fgLeastWeight);
 752       if (minQuality > 0 && fit->fQuality < minQuality) continue;
 753       nGood++;
 754
 755       // Check this bin
 756       CACHE(j,i,offset) = j;
 757       minBin            = TMath::Min(j, minBin);
 758       maxBin            = TMath::Max(j, maxBin);
 759     }
 760     AliInfoF("Out of %d bins, %d had fits, of which %d are good (%5.1f%%)",
 761              offset, nFits, nGood, (nFits > 0 ? 100*float(nGood)/nFits : 0));
 762
 763     // Now loop and set neighbors
 764     realMinBin = TMath::Max(1,      realMinBin-1); // Include one more
 765     realMaxBin = TMath::Min(offset, realMaxBin+1); // Include one more
 766     for (Int_t j = realMinBin; j <= realMaxBin; j++) {
 767       if (CACHE(j,i,offset) > 0) continue;
 768
 769       Int_t nK    = TMath::Max(realMaxBin - j, j - realMinBin);
 770       Int_t found = -1;
 771       for (Int_t k = 1; k <= nK; k++) {
 772         Int_t left  = j - k;
 773         Int_t right = j + k;
 774         if      (left  > realMinBin &&
 775                  CACHE(left,i,offset)  == left) found = left;
 776         else if (right < realMaxBin &&
 777                  CACHE(right,i,offset) == right) found = right;
 778         if (found > 0) break;
 779       }
 780       // Now check that we found something
 781       if (found) CACHE(j,i,offset) = CACHE(found,i,offset);
 782       else AliWarningF("No fit found for etabin=%d in ring=%d", j, i);
 783     }
 784   }
 785 }
 786
 787
 788 //____________________________________________________________________
 789 Int_t
 790 AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
 791 {
 792   //
 793   // Find the eta bin corresponding to the given eta
 794   //
 795   // Parameters:
 796   //    eta  Eta value
 797   //
 798   // Return:
 799   //    Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
 800   // eta, or 0 if out of range.
 801   //
 802   if (TMath::Abs(fEtaAxis.GetXmin() - fEtaAxis.GetXmax()) < 1e-6
 803       || fEtaAxis.GetNbins() == 0) {
 804     AliWarning("No eta axis defined");
 805     return -1;
 806   }
 807   Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
 808   if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
 809   return bin;
 810 }
 811
 812 //____________________________________________________________________
 813 Bool_t
 814 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
 815 {
 816   //
 817   // Set the fit parameters from a function
 818   //
 819   // Parameters:
 820   //    d       Detector
 821   //    r       Ring
 822   //    etaBin  Eta (bin number, 1->nBins)
 823   //    f       ELoss fit result - note, the object will take ownership
 824   //
 825   TObjArray* ringArray = GetOrMakeRingArray(d, r);
 826   if (!ringArray) {
 827     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
 828     return kFALSE;
 829   }
 830   if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
 831     AliError(Form("bin=%d is out of range [%d,%d]",
 832                   etaBin, 1, fEtaAxis.GetNbins()));
 833     return kFALSE;
 834   }
 835   // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
 836   ringArray->AddAtAndExpand(fit, etaBin);
 837   return kTRUE;
 838 }
 839
 840 //____________________________________________________________________
 841 Bool_t
 842 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
 843 {
 844   //
 845   // Set the fit parameters from a function
 846   //
 847   // Parameters:
 848   //    d    Detector
 849   //    r    Ring
 850   //    eta  Eta
 851   //    f    ELoss fit result - note, the object will take ownership
 852   //
 853   Int_t bin = FindEtaBin(eta);
 854   if (bin <= 0) {
 855     AliError(Form("eta=%f is out of range [%f,%f]",
 856                   eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
 857     return kFALSE;
 858   }
 859
 860   return SetFit(d, r, bin, fit);
 861 }
 862 //____________________________________________________________________
 863 Bool_t
 864 AliFMDCorrELossFit::SetFit(UShort_t  d,      Char_t    r,
 865                            Double_t  eta,
 866                            Int_t     quality,UShort_t  n,
 867                            Double_t  chi2,   UShort_t  nu,
 868                            Double_t  c,      Double_t  ec,
 869                            Double_t  delta,  Double_t  edelta,
 870                            Double_t  xi,     Double_t  exi,
 871                            Double_t  sigma,  Double_t  esigma,
 872                            Double_t  sigman, Double_t  esigman,
 873                            Double_t* a,      Double_t* ea)
 874 {
 875   //
 876   // Set the fit parameters from a function
 877   //
 878   // Parameters:
 879   //    d         Detector number
 880   //    r         Ring identifier
 881   //    eta       Eta value
 882   //    quality   Quality flag
 883   //    n         @f$ N@f$ - Number of fitted peaks
 884   //    chi2      @f$ \chi^2 @f$
 885   //    nu        @f$ \nu @f$ - number degrees of freedom
 886   //    c         @f$ C@f$ - scale constant
 887   //    ec        @f$ \delta C@f$ - error on @f$ C@f$
 888   //    delta     @f$ \Delta@f$ - most probable value
 889   //    edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$
 890   //    xi        @f$ \xi@f$ - Landau width
 891   //    exi       @f$ \delta\xi@f$ - error on @f$\xi@f$
 892   //    sigma     @f$ \sigma@f$ - Gaussian width
 893   //    esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$
 894   //    sigman    @f$ \sigma_n@f$ - Noise width
 895   //    esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$
 896   //    a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for
 897   //                  @f$ i=2,\ldots@f$
 898   //    ea        Array of @f$ N-1@f$ errors on weights @f$ a_i@f$ for
 899   //                  @f$ i=2,\ldots@f$
 900   //
 901   ELossFit* e = new ELossFit(quality, n,
 902                              chi2,    nu,
 903                              c,       ec,
 904                              delta,   edelta,
 905                              xi,      exi,
 906                              sigma,   esigma,
 907                              sigman,  esigman,
 908                              a,       ea);
 909   if (!SetFit(d, r, eta, e)) {
 910     delete e;
 911     return kFALSE;
 912   }
 913   return kTRUE;
 914 }
 915 //____________________________________________________________________
 916 Bool_t
 917 AliFMDCorrELossFit::SetFit(UShort_t  d, Char_t r, Double_t eta,
 918                            Int_t quality, const TF1& f)
 919 {
 920   //
 921   // Set the fit parameters from a function
 922   //
 923   // Parameters:
 924   //    d        Detector
 925   //    r        Ring
 926   //    eta      Eta
 927   //    quality  Quality flag
 928   //    f        Function from fit
 929   //
 930   ELossFit* e = new ELossFit(quality, f);
 931   if (!SetFit(d, r, eta, e)) {
 932     delete e;
 933     return kFALSE;
 934   }
 935   return kTRUE;
 936 }
 937 //____________________________________________________________________
 938 AliFMDCorrELossFit::ELossFit*
 939 AliFMDCorrELossFit::GetFit(UShort_t  d, Char_t r, Int_t etabin) const
 940 {
 941   //
 942   // Get the fit corresponding to the specified parameters
 943   //
 944   // Parameters:
 945   //    d      Detector
 946   //    r      Ring
 947   //    etabin Eta bin (1 based)
 948   //
 949   // Return:
 950   //    Fit parameters or null in case of problems
 951   //
 952   TObjArray* ringArray = GetRingArray(d, r);
 953   if (!ringArray)                                   return 0;
 954   if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) return 0;
 955   if      (etabin > ringArray->GetEntriesFast())    return 0;
 956   else if (etabin >= ringArray->GetEntriesFast())   etabin--;
 957   else if (!ringArray->At(etabin))                  etabin++;
 958   return static_cast<ELossFit*>(ringArray->At(etabin));
 959 }
 960 //____________________________________________________________________
 961 AliFMDCorrELossFit::ELossFit*
 962 AliFMDCorrELossFit::GetFit(UShort_t  d, Char_t r, Double_t eta) const
 963 {
 964   //
 965   // Find the fit corresponding to the specified parameters
 966   //
 967   // Parameters:
 968   //    d   Detector
 969   //    r   Ring
 970   //    eta Eta value
 971   //
 972   // Return:
 973   //    Fit parameters or null in case of problems
 974   //
 975   Int_t etabin = FindEtaBin(eta);
 976   return GetFit(d, r, etabin);
 977 }
 978
 979 //____________________________________________________________________
 980 AliFMDCorrELossFit::ELossFit*
 981 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Int_t etabin,
 982                             UShort_t minQ) const
 983 {
 984   //
 985   // Find the fit corresponding to the specified parameters
 986   //
 987   // Parameters:
 988   //    d      Detector
 989   //    r      Ring
 990   //    etabin Eta bin (1 based)
 991   //
 992   // Return:
 993   //    Fit parameters or null in case of problems
 994   //
 995   if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
 996     // AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
 997     //               etabin, 1, fEtaAxis.GetNbins(), d, r));
 998     return 0;
 999   }
1000
1001   TObjArray* ringArray = GetRingArray(d, r);
1002   if (!ringArray) {
1003     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
1004     return 0;
1005   }
1006   if (fCache.GetSize() <= 0) CacheBins(minQ);
1007   Int_t idx = (d == 1 ? 0 :
1008                (d - 2) * 2 + 1 + (r=='I' || r=='i' ? 0 : 1));
1009   Int_t bin = CACHE(etabin, idx, fEtaAxis.GetNbins());
1010
1011   if (bin < 0 || bin > ringArray->GetEntriesFast()) return 0;
1012
1013   return static_cast<ELossFit*>(ringArray->At(bin));
1014 }
1015 //____________________________________________________________________
1016 AliFMDCorrELossFit::ELossFit*
1017 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Double_t eta,
1018                             UShort_t minQ) const
1019 {
1020   //
1021   // Find the fit corresponding to the specified parameters
1022   //
1023   // Parameters:
1024   //    d   Detector
1025   //    r   Ring
1026   //    eta Eta value
1027   //
1028   // Return:
1029   //    Fit parameters or null in case of problems
1030   //
1031   Int_t etabin = FindEtaBin(eta);
1032   return FindFit(d, r, etabin, minQ);
1033 }
1034 //____________________________________________________________________
1035 TObjArray*
1036 AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
1037 {
1038   //
1039   // Get the ring array corresponding to the specified ring
1040   //
1041   // Parameters:
1042   //    d Detector
1043   //    r Ring
1044   //
1045   // Return:
1046   //    Pointer to ring array, or null in case of problems
1047   //
1048   Int_t idx = -1;
1049   switch (d) {
1050   case 1:   idx = 0; break;
1051   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
1052   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
1053   }
1054   if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
1055   return static_cast<TObjArray*>(fRings.At(idx));
1056 }
1057 //____________________________________________________________________
1058 TObjArray*
1059 AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
1060 {
1061   //
1062   // Get the ring array corresponding to the specified ring
1063   //
1064   // Parameters:
1065   //    d Detector
1066   //    r Ring
1067   //
1068   // Return:
1069   //    Pointer to ring array, or newly created container
1070   //
1071   Int_t idx = -1;
1072   switch (d) {
1073   case 1:   idx = 0; break;
1074   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
1075   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
1076   }
1077   if (idx < 0) return 0;
1078   if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
1079     TObjArray* a = new TObjArray(0);
1080     // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
1081     a->SetName(Form("FMD%d%c", d, r));
1082     a->SetOwner();
1083     fRings.AddAtAndExpand(a, idx);
1084   }
1085   return static_cast<TObjArray*>(fRings.At(idx));
1086 }
1087
1088 //____________________________________________________________________
1089 Double_t
1090 AliFMDCorrELossFit::GetLowerBound(UShort_t  d, Char_t r, Int_t etabin,
1091                                   Double_t f) const
1092 {
1093   ELossFit* fit = FindFit(d, r, etabin, 20);
1094   if (!fit) return -1024;
1095   return fit->GetLowerBound(f);
1096 }
1097 //____________________________________________________________________
1098 Double_t
1099 AliFMDCorrELossFit::GetLowerBound(UShort_t  d, Char_t r, Double_t eta,
1100                                   Double_t f) const
1101 {
1102   Int_t bin = FindEtaBin(eta);
1103   if (bin <= 0) return -1024;
1104   return GetLowerBound(d, r, Int_t(bin), f);
1105 }
1106 //____________________________________________________________________
1107 Double_t
1108 AliFMDCorrELossFit::GetLowerBound(UShort_t  d, Char_t r, Int_t etabin,
1109                                   Double_t f, Bool_t showErrors,
1110                                   Bool_t includeSigma) const
1111 {
1112   ELossFit* fit = FindFit(d, r, etabin, 20);
1113   if (!fit) {
1114     if (showErrors) {
1115       AliWarning(Form("No fit for FMD%d%c @ etabin=%d", d, r, etabin));
1116     }
1117     return -1024;
1118   }
1119   return fit->GetLowerBound(f, includeSigma);
1120 }
1121
1122 //____________________________________________________________________
1123 Double_t
1124 AliFMDCorrELossFit::GetLowerBound(UShort_t  d, Char_t r, Double_t eta,
1125                                   Double_t f, Bool_t showErrors,
1126                                   Bool_t includeSigma) const
1127 {
1128   Int_t bin = FindEtaBin(eta);
1129   if (bin <= 0) {
1130     if (showErrors)
1131       AliError(Form("eta=%f out of bounds for FMD%d%c", eta, d, r));
1132     return -1024;
1133   }
1134   return GetLowerBound(d, r, bin, f, showErrors, includeSigma);
1135 }
1136
1137 //____________________________________________________________________
1138 namespace {
1139   TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
1140                  Int_t color)
1141   {
1142     TH1D* h = new TH1D(Form("%s_%s", name, title),
1143                        Form("%s %s", name, title),
1144                        axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
1145     h->SetDirectory(0);
1146     h->SetMarkerStyle(20);
1147     h->SetMarkerColor(color);
1148     h->SetMarkerSize(0.5);
1149     h->SetFillColor(color);
1150     h->SetFillStyle(3001);
1151     h->SetLineColor(color);
1152     return h;
1153   }
1154 }
1155
1156
1157
1158 #define IDX2RING(I) (i == 0 || i == 1 || i == 3 ? 'I' : 'O')
1159 #define IDX2DET(I)  (i == 0 ? 1 : (i == 1 || i == 2 ? 2 : 3))
1160 //____________________________________________________________________
1161 TList*
1162 AliFMDCorrELossFit::GetStacks(Bool_t   err,
1163                               Bool_t   rel,
1164                               Bool_t   /*good*/,
1165                               UShort_t maxN) const
1166 {
1167   // Get a list of THStacks
1168   Int_t nRings = fRings.GetEntriesFast();
1169   // Int_t nPad   = 6+maxN-1; // 7 regular params, and maxN-1 weights
1170
1171   enum {
1172     kChi2nu = 0,
1173     kC      = 1,
1174     kDelta  = 2,
1175     kXi     = 3,
1176     kSigma  = 4,
1177     kN      = 5
1178   };
1179
1180   THStack* sChi2nu;
1181   THStack* sC;
1182   THStack* sDelta;
1183   THStack* sXi;
1184   THStack* sSigma;
1185   // THStack* sigman;
1186   THStack* n;
1187   TList* stacks = new TList;
1188   stacks->AddAt(sChi2nu= new THStack("chi2",   "#chi^{2}/#nu"), kChi2nu);
1189   stacks->AddAt(sC     = new THStack("c",       "C"),           kC);
1190   stacks->AddAt(sDelta = new THStack("delta",  "#Delta_{mp}"),  kDelta);
1191   stacks->AddAt(sXi    = new THStack("xi",     "#xi"),          kXi);
1192   stacks->AddAt(sSigma = new THStack("sigma",  "#sigma"),       kSigma);
1193   //stacks->AddAt(sigman= new THStack("sigman", "#sigma_{n}"),   5);
1194   stacks->AddAt(n      = new THStack("n",      "N"),            kN);
1195   if (rel) {
1196     sChi2nu->SetName("qual");
1197     sChi2nu->SetTitle("Quality");
1198     n->SetName("good");
1199     n->SetTitle("Bin map");
1200   }
1201   for (Int_t i = 1; i <= maxN; i++) {
1202     stacks->AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), kN+i);
1203   }
1204
1205   // TArrayD min(nPad);
1206   // TArrayD max(nPad);
1207   // min.Reset(100000);
1208   // max.Reset(-100000);
1209
1210   for (Int_t i = 0; i < nRings; i++) {
1211     if (!fRings.At(i)) continue;
1212     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
1213     Int_t      nFits = a->GetEntriesFast();
1214     Int_t      color = AliForwardUtil::RingColor(IDX2DET(i), IDX2RING(i));
1215
1216     TH1D* hChi    = MakeHist(fEtaAxis,a->GetName(), "chi2",   color);
1217     TH1D* hC      = MakeHist(fEtaAxis,a->GetName(), "c",      color);
1218     TH1D* hDelta  = MakeHist(fEtaAxis,a->GetName(), "delta",  color);
1219     TH1D* hXi     = MakeHist(fEtaAxis,a->GetName(), "xi",     color);
1220     TH1D* hSigma  = MakeHist(fEtaAxis,a->GetName(), "sigma",  color);
1221     // TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
1222     TH1D* hN      = MakeHist(fEtaAxis,a->GetName(), "n",      color);
1223     TH1D* hA[maxN];
1224     const char* ho = (rel || !err ? "hist" : "e");
1225     sChi2nu->Add(hChi,   "hist"); // 0
1226     sC     ->Add(hC,     ho);     // 1
1227     sDelta ->Add(hDelta, ho);     // 2
1228     sXi    ->Add(hXi,    ho);     // 3
1229     sSigma ->Add(hSigma, ho);     // 4
1230     // sigman->Add(hSigmaN,ho);     // 5
1231     n     ->Add(hN,     "hist"); // 5
1232     hChi->SetFillColor(color);
1233     hChi->SetFillStyle(3001);
1234     hN->SetFillColor(color);
1235     hN->SetFillStyle(3001);
1236
1237     for (Int_t k = 1; k <= maxN; k++) {
1238       hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
1239       static_cast<THStack*>(stacks->At(kN+k))->Add(hA[k-1], ho);
1240     }
1241
1242     for (Int_t j = 0; j < nFits; j++) {
1243       ELossFit* f = static_cast<ELossFit*>(a->At(j));
1244       if (!f) continue;
1245
1246       Int_t     b       = f->fBin;
1247       Double_t  vChi2nu = (rel ? f->fQuality
1248                            : (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
1249       Double_t  vC      = (rel ? (f->fC     >0 ?f->fEC     /f->fC      :0)
1250                            : f->fC);
1251       Double_t  vDelta  = (rel ? (f->fDelta >0 ?f->fEDelta /f->fDelta  :0)
1252                            : f->fDelta);
1253       Double_t  vXi     = (rel ? (f->fXi    >0 ?f->fEXi    /f->fXi     :0)
1254                            : f->fXi);
1255       Double_t  vSigma  = (rel ? (f->fSigma >0 ?f->fESigma /f->fSigma  :0)
1256                            : f->fSigma);
1257       Int_t     nW      = (rel ? CACHE(j,i,fEtaAxis.GetNbins()) :
1258                            f->FindMaxWeight());
1259       // Double_t  sigman = (rel ? (f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0)
1260       //                     : f->SigmaN);
1261       hChi   ->SetBinContent(b, vChi2nu);
1262       hN     ->SetBinContent(b, nW);
1263       hC     ->SetBinContent(b, vC);
1264       hDelta ->SetBinContent(b, vDelta);
1265       hXi    ->SetBinContent(b, vXi);
1266       hSigma ->SetBinContent(b, vSigma);
1267       // if (vChi2nu > 1e-12) {
1268       //        min[kChi2nu] = TMath::Min(min[kChi2nu], vChi2nu);
1269       //        max[kChi2nu] = TMath::Max(max[kChi2nu], vChi2nu);
1270       // }
1271       // if (vC > 1e-12) {
1272       //        min[kC] = TMath::Min(min[kC], vC);
1273       //        max[kC] = TMath::Max(max[kC], vC);
1274       // }
1275       // if (vDelta > 1e-12) {
1276       //        min[kDelta] = TMath::Min(min[kDelta], vDelta);
1277       //        max[kDelta] = TMath::Max(max[kDelta], vDelta);
1278       // }
1279       // if (vXi > 1e-12) {
1280       //        min[kXi] = TMath::Min(min[kXi], vXi);
1281       //        max[kXi] = TMath::Max(max[kXi], vXi);
1282       // }
1283       // if (vSigma > 1e-12) {
1284       //        min[kSigma] = TMath::Min(min[kSigma], vSigma);
1285       //        max[kSigma] = TMath::Max(max[kSigma], vSigma);
1286       // }
1287       // if (nW > 1e-12) {
1288       //        min[kN] = TMath::Min(min[kN], Double_t(nW));
1289       //        max[kN] = TMath::Max(max[kN], Double_t(nW));
1290       // }
1291       // hSigmaN->SetBinContent(b,  sigman);
1292       if (!rel) {
1293         hC     ->SetBinError(b, f->fEC);
1294         hDelta ->SetBinError(b, f->fEDelta);
1295         hXi    ->SetBinError(b, f->fEXi);
1296         hSigma ->SetBinError(b, f->fESigma);
1297         // hSigmaN->SetBinError(b, f->fESigmaN);
1298       }
1299       for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
1300         Double_t vA = (rel ? (f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0)
1301                        : f->fA[k]);
1302         hA[k]->SetBinContent(b, vA);
1303         if (!rel) hA[k]->SetBinError(b, f->fEA[k]);
1304         // if (vA > 1e-12) {
1305         //   min[kN+1+k] = TMath::Min(min[kN+1+k], vA);
1306         //   max[kN+1+k] = TMath::Max(max[kN+1+k], vA);
1307         // }
1308       }
1309     }
1310   }
1311   return stacks;
1312 }
1313
1314 //____________________________________________________________________
1315 void
1316 AliFMDCorrELossFit::Draw(Option_t* option)
1317 {
1318   //
1319   // Draw this object
1320   //
1321   // Parameters:
1322   //    option Options.  Possible values are
1323   //  - err Plot error bars
1324   //
1325   TString opt(Form("nostack %s", option));
1326   opt.ToLower();
1327   Bool_t  rel = (opt.Contains("relative"));
1328   Bool_t  err = (opt.Contains("error"));
1329   Bool_t  clr = (opt.Contains("clear"));
1330   Bool_t  gdd = (opt.Contains("good"));
1331   if (rel) opt.ReplaceAll("relative","");
1332   if (err) opt.ReplaceAll("error","");
1333   if (clr) opt.ReplaceAll("clear", "");
1334   if (gdd) opt.ReplaceAll("good", "");
1335
1336   UShort_t maxN   = 0;
1337   Int_t nRings = fRings.GetEntriesFast();
1338   for (Int_t i = 0; i < nRings; i++) {
1339     if (!fRings.At(i)) continue;
1340     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
1341     Int_t      nFits = a->GetEntriesFast();
1342
1343     for (Int_t j = 0; j < nFits; j++) {
1344       ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1345       if (!fit) continue;
1346       maxN          = TMath::Max(maxN, UShort_t(fit->fN));
1347     }
1348   }
1349   // AliInfo(Form("Maximum N is %d", maxN));
1350   Int_t nPad = 6+maxN-1; // 7 regular params, and maxN-1 weights
1351   TVirtualPad* pad = gPad;
1352   if (clr) {
1353     pad->Clear();
1354     pad->SetTopMargin(0.02);
1355     pad->SetRightMargin(0.02);
1356     pad->SetBottomMargin(0.15);
1357     pad->SetLeftMargin(0.10);
1358   }
1359   pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
1360
1361   TList* stacks = GetStacks(err, rel, gdd, maxN);
1362
1363   Int_t nPad2 = (nPad+1) / 2;
1364   for (Int_t i = 0; i < nPad; i++) {
1365     Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
1366     TVirtualPad* p = pad->cd(iPad);
1367     p->SetLeftMargin(.15);
1368     p->SetFillColor(0);
1369     p->SetFillStyle(0);
1370     p->SetGridx();
1371     p->SetGridy();
1372     if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
1373
1374
1375     THStack* stack = static_cast<THStack*>(stacks->At(i));
1376
1377     // Double_t powMax = TMath::Log10(max[i]);
1378     // Double_t powMin = min[i] <= 0 ? powMax : TMath::Log10(min[i]);
1379     // if (powMax-powMin > 2. && min[i] != 0) p->SetLogy();
1380
1381     // stack->SetMinimum(min[i]);
1382     // stack->SetMaximum(max[i]);
1383     stack->Draw(opt.Data());
1384
1385     TString tit(stack->GetTitle());
1386     if (rel && i != 0 && i != 5)
1387       tit = Form("#delta %s/%s", tit.Data(), tit.Data());
1388     TH1*   hist  = stack->GetHistogram();
1389     TAxis* yaxis = hist->GetYaxis();
1390     yaxis->SetTitle(tit.Data());
1391     yaxis->SetTitleSize(0.15);
1392     yaxis->SetLabelSize(0.08);
1393     yaxis->SetTitleOffset(0.35);
1394     yaxis->SetTitleFont(132);
1395     yaxis->SetLabelFont(132);
1396     yaxis->SetNdivisions(5);
1397
1398
1399     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
1400     xaxis->SetTitle("#eta");
1401     xaxis->SetTitleSize(0.15);
1402     xaxis->SetLabelSize(0.08);
1403     xaxis->SetTitleOffset(0.35);
1404     xaxis->SetTitleFont(132);
1405     xaxis->SetLabelFont(132);
1406     xaxis->SetNdivisions(10);
1407
1408     stack->Draw(opt.Data());
1409   }
1410   pad->cd();
1411 }
1412
1413 //____________________________________________________________________
1414 void
1415 AliFMDCorrELossFit::Print(Option_t* option) const
1416 {
1417   //
1418   // Print this object.
1419   //
1420   // Parameters:
1421   //    option Options
1422   //   - R   Print recursive
1423   //
1424   //
1425   TString opt(option);
1426   opt.ToUpper();
1427   Int_t nRings  = fRings.GetEntriesFast();
1428   bool  recurse = opt.Contains("R");
1429   bool  cache   = opt.Contains("C") && fCache.GetSize() > 0;
1430   Int_t nBins   = fEtaAxis.GetNbins();
1431
1432   std::cout << "Low cut in fit range: " << fLowCut << "\n"
1433             << "Eta axis:             " << nBins
1434             << " bins, range [" << fEtaAxis.GetXmin() << ","
1435             << fEtaAxis.GetXmax() << "]" << std::endl;
1436
1437   for (Int_t i = 0; i < nRings; i++) {
1438     if (!fRings.At(i)) continue;
1439     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
1440     Int_t      nFits = a->GetEntriesFast();
1441
1442     std::cout << a->GetName() << " [" << nFits << " entries]"
1443               << (recurse ? ":\n" : "\t");
1444     Int_t min = fEtaAxis.GetNbins()+1;
1445     Int_t max = 0;
1446     for (Int_t j = 0; j < nFits; j++) {
1447       if (!a->At(j)) continue;
1448
1449       min = TMath::Min(j, min);
1450       max = TMath::Max(j, max);
1451
1452       if (recurse) {
1453         std::cout << "Bin # " << j << "\t";
1454         ELossFit* fit = static_cast<ELossFit*>(a->At(j));
1455         fit->Print(option);
1456       }
1457     }
1458     if (!recurse)
1459       std::cout << " bin range: " << std::setw(3) << min
1460                 << "-" << std::setw(3) << max << " " << std::setw(3)
1461                 << (max-min+1) << " bins" << std::endl;
1462   }
1463
1464   if (!cache) return;
1465
1466   std::cout << " eta bin           |           Fit bin              \n"
1467             << " #       range     | FMD1i  FMD2i  FMD2o  FMD3i  FMD3o"
1468     // << "----+-----+++------+-----------------------------------"
1469             << std::endl;
1470   size_t oldPrec = std::cout.precision();
1471   std::cout.precision(3);
1472   for (Int_t i = 1; i <= nBins; i++) {
1473     std::cout << std::setw(4) << i << " "
1474               << std::setw(5) << std::showpos << fEtaAxis.GetBinLowEdge(i)
1475               << " - " << std::setw(5) << fEtaAxis.GetBinUpEdge(i)
1476               << std::noshowpos << " | ";
1477     for (Int_t j = 0; j < 5; j++) {
1478       Int_t bin = CACHE(i,j,nBins);
1479       if (bin <= 0) std::cout << "       ";
1480       else          std::cout << std::setw(5) << bin
1481                               << (bin == i ? ' ' : '*') << ' ';
1482     }
1483     std::cout << std::endl;
1484   }
1485   std::cout.precision(oldPrec);
1486 }
1487
1488 //____________________________________________________________________
1489 void
1490 AliFMDCorrELossFit::Browse(TBrowser* b)
1491 {
1492   //
1493   // Browse this object
1494   //
1495   // Parameters:
1496   //    b
1497   //
1498   b->Add(&fRings);
1499   b->Add(&fEtaAxis);
1500 }
1501
1502
1503
1504 //____________________________________________________________________
1505 //
1506 // EOF
1507 //