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