PWGLF/FORWARD/analysis2/AliFMDCorrELossFit.cxx

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