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