PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx

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