PWG2/FORWARD/analysis2/AliFMDCorrELossFit.cxx

   1 #include "AliFMDCorrELossFit.h"
   2 #include "AliForwardUtil.h"
   3 #include <TF1.h>
   4 #include <TBrowser.h>
   5 #include <TVirtualPad.h>
   6 #include <THStack.h>
   7 #include <TH1D.h>
   8 #include <AliLog.h>
   9 #include <TList.h>
  10 #include <iostream>
  11 #include <iomanip>
  12
  13 Double_t AliFMDCorrELossFit::ELossFit::fgMaxRelError = .2;
  14 Double_t AliFMDCorrELossFit::ELossFit::fgLeastWeight = 1e-5;
  15 Double_t AliFMDCorrELossFit::ELossFit::fgMaxChi2nu   = 20;
  16
  17 //____________________________________________________________________
  18 AliFMDCorrELossFit::ELossFit::ELossFit()
  19   : fN(0),
  20     fNu(0),
  21     fChi2(0),
  22     fC(0),
  23     fDelta(0),
  24     fXi(0),
  25     fSigma(0),
  26     fSigmaN(0),
  27     fA(0),
  28     fEC(0),
  29     fEDelta(0),
  30     fEXi(0),
  31     fESigma(0),
  32     fESigmaN(0),
  33     fEA(0),
  34     fQuality(0),
  35     fDet(0),
  36     fRing('\0'),
  37     fBin(0)
  38 {}
  39 //____________________________________________________________________
  40 AliFMDCorrELossFit::ELossFit::ELossFit(Int_t quality, const TF1& f)
  41   : fN(f.GetNpar() > AliForwardUtil::ELossFitter::kN ?
  42        f.GetParameter(AliForwardUtil::ELossFitter::kN) :
  43        1),
  44     fNu(f.GetNDF()),
  45     fChi2(f.GetChisquare()),
  46     fC(f.GetParameter(AliForwardUtil::ELossFitter::kC)),
  47     fDelta(f.GetParameter(AliForwardUtil::ELossFitter::kDelta)),
  48     fXi(f.GetParameter(AliForwardUtil::ELossFitter::kXi)),
  49     fSigma(f.GetParameter(AliForwardUtil::ELossFitter::kSigma)),
  50     fSigmaN(f.GetParameter(AliForwardUtil::ELossFitter::kSigmaN)),
  51     fA(0),
  52     fEC(f.GetParError(AliForwardUtil::ELossFitter::kC)),
  53     fEDelta(f.GetParError(AliForwardUtil::ELossFitter::kDelta)),
  54     fEXi(f.GetParError(AliForwardUtil::ELossFitter::kXi)),
  55     fESigma(f.GetParError(AliForwardUtil::ELossFitter::kSigma)),
  56     fESigmaN(f.GetParError(AliForwardUtil::ELossFitter::kSigmaN)),
  57     fEA(0),
  58     fQuality(quality),
  59     fDet(0),
  60     fRing('\0'),
  61     fBin(0)
  62 {
  63   if (fN <= 0) return;
  64   fA  = new Double_t[fN];
  65   fEA = new Double_t[fN];
  66   for (Int_t i = 0; i < fN-1; i++) {
  67     fA[i]  = f.GetParameter(AliForwardUtil::ELossFitter::kA+i);
  68     fEA[i] = f.GetParError(AliForwardUtil::ELossFitter::kA+i);
  69   }
  70   fA[fN-1]  = -9999;
  71   fEA[fN-1] = -9999;
  72 }
  73
  74 //____________________________________________________________________
  75 AliFMDCorrELossFit::ELossFit::ELossFit(Int_t     quality,UShort_t  n,
  76                                        Double_t  chi2,   UShort_t  nu,
  77                                        Double_t  c,      Double_t  ec,
  78                                        Double_t  delta,  Double_t  edelta,
  79                                        Double_t  xi,     Double_t  exi,
  80                                        Double_t  sigma,  Double_t  esigma,
  81                                        Double_t  sigman, Double_t  esigman,
  82                                        Double_t* a,      Double_t* ea)
  83   : fN(n),
  84     fNu(nu),
  85     fChi2(chi2),
  86     fC(c),
  87     fDelta(delta),
  88     fXi(xi),
  89     fSigma(sigma),
  90     fSigmaN(sigman),
  91     fA(0),
  92     fEC(ec),
  93     fEDelta(edelta),
  94     fEXi(exi),
  95     fESigma(esigma),
  96     fESigmaN(esigman),
  97     fEA(0),
  98     fQuality(quality),
  99     fDet(0),
 100     fRing('\0'),
 101     fBin(0)
 102 {
 103   if (fN <= 0) return;
 104   fA  = new Double_t[fN];
 105   fEA = new Double_t[fN];
 106   for (Int_t i = 0; i < fN-1; i++) {
 107     fA[i]  = a[i];
 108     fEA[i] = ea[i];
 109   }
 110   fA[fN-1]  = -9999;
 111   fEA[fN-1] = -9999;
 112 }
 113 //____________________________________________________________________
 114 AliFMDCorrELossFit::ELossFit::ELossFit(const ELossFit& o)
 115   : TObject(o),
 116     fN(o.fN),
 117     fNu(o.fNu),
 118     fChi2(o.fChi2),
 119     fC(o.fC),
 120     fDelta(o.fDelta),
 121     fXi(o.fXi),
 122     fSigma(o.fSigma),
 123     fSigmaN(o.fSigmaN),
 124     fA(0),
 125     fEC(o.fEC),
 126     fEDelta(o.fEDelta),
 127     fEXi(o.fEXi),
 128     fESigma(o.fESigma),
 129     fESigmaN(o.fESigmaN),
 130     fEA(0),
 131     fQuality(o.fQuality),
 132     fDet(o.fDet),
 133     fRing(o.fRing),
 134     fBin(o.fBin)
 135 {
 136   if (fN <= 0) return;
 137   fA  = new Double_t[fN];
 138   fEA = new Double_t[fN];
 139   for (Int_t i = 0; i < fN-1; i++) {
 140     fA[i]  = o.fA[i];
 141     fEA[i] = o.fEA[i];
 142   }
 143   fA[fN-1]  = -9999;
 144   fEA[fN-1] = -9999;
 145 }
 146
 147 //____________________________________________________________________
 148 AliFMDCorrELossFit::ELossFit&
 149 AliFMDCorrELossFit::ELossFit::operator=(const ELossFit& o)
 150 {
 151   fN       = o.fN;
 152   fNu      = o.fNu;
 153   fChi2    = o.fChi2;
 154   fC       = o.fC;
 155   fDelta   = o.fDelta;
 156   fXi      = o.fXi;
 157   fSigma   = o.fSigma;
 158   fSigmaN  = o.fSigmaN;
 159   fEC      = o.fEC;
 160   fEDelta  = o.fEDelta;
 161   fEXi     = o.fEXi;
 162   fESigma  = o.fESigma;
 163   fESigmaN = o.fESigmaN;
 164   fQuality = o.fQuality;
 165   fDet     = o.fDet;
 166   fRing    = o.fRing;
 167   fBin     = o.fBin;
 168   if (fA)  delete [] fA;
 169   if (fEA) delete [] fEA;
 170   fA  = 0;
 171   fEA = 0;
 172
 173   if (fN <= 0) return *this;
 174   fA  = new Double_t[fN];
 175   fEA = new Double_t[fN];
 176   for (Int_t i = 0; i < fN; i++) {
 177     fA[i]  = o.fA[i];
 178     fEA[i] = o.fEA[i];
 179   }
 180
 181   return *this;
 182 }
 183
 184 //____________________________________________________________________
 185 AliFMDCorrELossFit::ELossFit::~ELossFit()
 186 {
 187   if (fA)  delete[] fA;
 188   if (fEA) delete[] fEA;
 189 }
 190
 191
 192 //____________________________________________________________________
 193 Int_t
 194 AliFMDCorrELossFit::ELossFit::FindMaxWeight(Double_t maxRelError,
 195                                             Double_t leastWeight,
 196                                             UShort_t  maxN) const
 197 {
 198   Int_t n = TMath::Min(maxN, UShort_t(fN-1));
 199   Int_t m = 1;
 200   // fN is one larger than we have data
 201   for (Int_t i = 0; i < n-1; i++, m++) {
 202     if (fA[i] < leastWeight)  break;
 203     if (fEA[i] / fA[i] > maxRelError) break;
 204   }
 205   return m;
 206 }
 207
 208 //____________________________________________________________________
 209 Double_t
 210 AliFMDCorrELossFit::ELossFit::Evaluate(Double_t x,
 211                                        UShort_t maxN) const
 212 {
 213   return AliForwardUtil::NLandauGaus(x, fDelta, fXi, fSigma, fSigmaN,
 214                                      TMath::Min(maxN, UShort_t(fN)), fA);
 215 }
 216
 217 //____________________________________________________________________
 218 Double_t
 219 AliFMDCorrELossFit::ELossFit::EvaluateWeighted(Double_t x,
 220                                                UShort_t maxN) const
 221 {
 222   UShort_t n   = TMath::Min(maxN, UShort_t(fN-1));
 223   Double_t num = 0;
 224   Double_t den = 0;
 225   for (Int_t i = 1; i <= n; i++) {
 226     Double_t a = (i == 1 ? 1 : fA[i-1]);
 227     if (fA[i-1] < 0) break;
 228     Double_t f = AliForwardUtil::ILandauGaus(x,fDelta,fXi,fSigma,fSigmaN,i);
 229     num += i * a * f;
 230     den += a * f;
 231   }
 232   if (den <= 0) return 1;
 233   return num / den;
 234 }
 235
 236
 237 #define OUTPAR(N,V,E)                   \
 238   std::setprecision(9) <<               \
 239   std::setw(12) << N << ": " <<         \
 240   std::setw(14) << V << " +/- " <<      \
 241   std::setw(14) << E << " (" <<         \
 242   std::setprecision(-1) <<              \
 243   std::setw(5)  << 100*(V>0?E/V:1) << "%)\n"
 244
 245
 246 //____________________________________________________________________
 247 Int_t
 248 AliFMDCorrELossFit::ELossFit::Compare(const TObject* o) const
 249 {
 250   const ELossFit* other = static_cast<const ELossFit*>(o);
 251   if (this->fQuality < other->fQuality) return -1;
 252   if (this->fQuality > other->fQuality) return +1;
 253   Double_t chi2nu  = (fNu == 0 ? 1000 : fChi2 / fNu);
 254   Double_t oChi2nu = (other->fNu == 0 ? 1000 : other->fChi2 / other->fNu);
 255   if (TMath::Abs(chi2nu-1) < TMath::Abs(oChi2nu-1)) return +1;
 256   if (TMath::Abs(chi2nu-1) > TMath::Abs(oChi2nu-1)) return -1;
 257   return 0;
 258 }
 259
 260 //____________________________________________________________________
 261 void
 262 AliFMDCorrELossFit::ELossFit::Print(Option_t*) const
 263 {
 264   std::cout << GetName() << ":\n"
 265             << " chi^2/nu = " << fChi2 << "/" << fNu << " = "
 266             << (fNu == 0 ? 999 : fChi2 / fNu) << "\n"
 267             << "     Quality:   " << fQuality << "\n"
 268             << "  NParticles:   " << fN << "  (" << FindMaxWeight() << ")\n"
 269             << OUTPAR("Delta", fDelta, fEDelta)
 270             << OUTPAR("xi", fXi, fEXi)
 271             << OUTPAR("sigma", fSigma, fESigma)
 272             << OUTPAR("sigma_n", fSigmaN, fESigmaN);
 273   for (Int_t i = 0; i < fN-1; i++)
 274     std::cout << OUTPAR(Form("a%d", i+2), fA[i], fEA[i]);
 275   std::cout << std::flush;
 276 }
 277
 278 //____________________________________________________________________
 279 const Char_t*
 280 AliFMDCorrELossFit::ELossFit::GetName() const
 281 {
 282   return Form("FMD%d%c_etabin%03d", fDet, fRing, fBin);
 283 }
 284
 285 //____________________________________________________________________
 286 void
 287 AliFMDCorrELossFit::ELossFit::Browse(TBrowser* b)
 288 {
 289   // Draw this one
 290   Draw(b ? b->GetDrawOption() : "comp");
 291   gPad->SetLogy();
 292   gPad->Update();
 293 }
 294
 295 //____________________________________________________________________
 296 void
 297 AliFMDCorrELossFit::ELossFit::Draw(Option_t* option)
 298 {
 299   TString opt(option);
 300   opt.ToUpper();
 301   bool comp = false;
 302   if (opt.Contains("COMP")) {
 303     opt.ReplaceAll("COMP","");
 304     comp = true;
 305   }
 306   if (!opt.Contains("SAME")) {
 307     gPad->Clear();
 308   }
 309
 310   TObjArray cleanup;
 311   TF1* tot = AliForwardUtil::MakeNLandauGaus(1,
 312                                              fDelta, fXi,
 313                                              fSigma, fSigmaN,
 314                                              fN,     fA,
 315                                              0.01,   10);
 316   tot->SetLineColor(kBlack);
 317   tot->SetLineWidth(2);
 318   tot->SetLineStyle(1);
 319   tot->SetTitle(GetName());
 320   Double_t max = tot->GetMaximum();
 321
 322   if (!opt.Contains("SAME")) {
 323     TH1* frame = new TH1F(GetName(),
 324                           Form("FMD%d%c, eta bin %d",fDet,fRing,fBin),
 325                           100, 0, 10);
 326     frame->SetMinimum(max/10000);
 327     frame->SetMaximum(max*1.4);
 328     frame->SetXTitle("#Delta / #Delta_{mip}");
 329     frame->Draw();
 330     opt.Append(" SAME");
 331   }
 332   tot->DrawCopy(opt.Data());
 333   cleanup.Add(tot);
 334
 335   if (!comp) {
 336     gPad->cd();
 337     return;
 338   }
 339
 340   Double_t min = max;
 341   opt.Append(" same");
 342   Int_t maxW = FindMaxWeight();
 343   for (Int_t i=1; i <= fN; i++) {
 344     TF1* f = AliForwardUtil::MakeILandauGaus((i == 1 ? 1 : fA[i-2]),
 345                                              fDelta, fXi,
 346                                              fSigma, fSigmaN,
 347                                              i,      0.01, 10);
 348     f->SetLineWidth(2);
 349     f->SetLineStyle(i > maxW ? 2 : 1);
 350     min = TMath::Min(f->GetMaximum(), min);
 351     f->DrawCopy(opt.Data());
 352     cleanup.Add(f);
 353   }
 354   min /= 100;
 355   tot->GetHistogram()->SetMaximum(max);
 356   tot->GetHistogram()->SetMinimum(min);
 357   tot->GetHistogram()->GetYaxis()->SetRangeUser(min, max);
 358
 359   gPad->cd();
 360 }
 361
 362
 363 //____________________________________________________________________
 364 #define CHECKPAR(V,E,T) ((V > 0) && (E / V < T))
 365
 366 //____________________________________________________________________
 367 void
 368 AliFMDCorrELossFit::ELossFit::CalculateQuality(Double_t maxChi2nu,
 369                                                Double_t maxRelError,
 370                                                Double_t leastWeight)
 371 {
 372   Int_t qual = 0;
 373   if (fNu > 0 && fChi2 / fNu < maxChi2nu) qual += 4;;
 374   if (CHECKPAR(fDelta,  fEDelta,  maxRelError)) qual++;
 375   if (CHECKPAR(fXi,     fEXi,     maxRelError)) qual++;
 376   if (CHECKPAR(fSigma,  fESigma,  maxRelError)) qual++;
 377   if (CHECKPAR(fSigmaN, fESigmaN, maxRelError)) qual++;
 378   qual += FindMaxWeight(1.5*maxRelError, leastWeight, fN);
 379   fQuality = qual;
 380 }
 381
 382 //____________________________________________________________________
 383 AliFMDCorrELossFit::AliFMDCorrELossFit()
 384   : TObject(),
 385     fRings(),
 386     fEtaAxis(0,0,0),
 387     fLowCut(0)
 388 {
 389   fRings.SetOwner(kTRUE);
 390   fEtaAxis.SetTitle("#eta");
 391   fEtaAxis.SetName("etaAxis");
 392   fRings.SetName("rings");
 393 }
 394
 395 //____________________________________________________________________
 396 AliFMDCorrELossFit::AliFMDCorrELossFit(const AliFMDCorrELossFit& o)
 397   : TObject(o),
 398     fRings(o.fRings),
 399     fEtaAxis(o.fEtaAxis.GetNbins(),o.fEtaAxis.GetXmin(),o.fEtaAxis.GetXmax()),
 400     fLowCut(0)
 401 {
 402   fEtaAxis.SetTitle("#eta");
 403   fEtaAxis.SetName("etaAxis");
 404 }
 405
 406 //____________________________________________________________________
 407 AliFMDCorrELossFit::~AliFMDCorrELossFit()
 408 {
 409   fRings.Clear();
 410 }
 411
 412 //____________________________________________________________________
 413 AliFMDCorrELossFit&
 414 AliFMDCorrELossFit::operator=(const AliFMDCorrELossFit& o)
 415 {
 416   fRings = o.fRings;
 417   fLowCut = o.fLowCut;
 418   SetEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax());
 419
 420   return *this;
 421 }
 422 //____________________________________________________________________
 423 Int_t
 424 AliFMDCorrELossFit::FindEtaBin(Double_t eta) const
 425 {
 426   if (fEtaAxis.GetXmin() == fEtaAxis.GetXmax() || fEtaAxis.GetNbins() == 0) {
 427     AliWarning("No eta axis defined");
 428     return -1;
 429   }
 430   Int_t bin = const_cast<TAxis&>(fEtaAxis).FindBin(eta);
 431   if (bin <= 0 || bin > fEtaAxis.GetNbins()) return 0;
 432   return bin;
 433 }
 434
 435 //____________________________________________________________________
 436 Bool_t
 437 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* fit)
 438 {
 439   TObjArray* ringArray = GetOrMakeRingArray(d, r);
 440   if (!ringArray) {
 441     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
 442     return kFALSE;
 443   }
 444   if (etaBin <= 0 || etaBin >= fEtaAxis.GetNbins()+1) {
 445     AliError(Form("bin=%d is out of range [%d,%d]",
 446                   etaBin, 1, fEtaAxis.GetNbins()));
 447     return kFALSE;
 448   }
 449   // AliInfo(Form("Adding fit %p at %3d", fit, etaBin));
 450   ringArray->AddAtAndExpand(fit, etaBin);
 451   return kTRUE;
 452 }
 453
 454 //____________________________________________________________________
 455 Bool_t
 456 AliFMDCorrELossFit::SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* fit)
 457 {
 458   Int_t bin = FindEtaBin(eta);
 459   if (bin <= 0) {
 460     AliError(Form("eta=%f is out of range [%f,%f]",
 461                   eta, fEtaAxis.GetXmin(), fEtaAxis.GetXmax()));
 462     return kFALSE;
 463   }
 464
 465   return SetFit(d, r, bin, fit);
 466 }
 467 //____________________________________________________________________
 468 Bool_t
 469 AliFMDCorrELossFit::SetFit(UShort_t  d,      Char_t    r,
 470                            Double_t  eta,
 471                            Int_t     quality,UShort_t  n,
 472                            Double_t  chi2,   UShort_t  nu,
 473                            Double_t  c,      Double_t  ec,
 474                            Double_t  delta,  Double_t  edelta,
 475                            Double_t  xi,     Double_t  exi,
 476                            Double_t  sigma,  Double_t  esigma,
 477                            Double_t  sigman, Double_t  esigman,
 478                            Double_t* a,      Double_t* ea)
 479 {
 480   ELossFit* e = new ELossFit(quality, n,
 481                              chi2,    nu,
 482                              c,       ec,
 483                              delta,   edelta,
 484                              xi,      exi,
 485                              sigma,   esigma,
 486                              sigman,  esigman,
 487                              a,       ea);
 488   if (!SetFit(d, r, eta, e)) {
 489     delete e;
 490     return kFALSE;
 491   }
 492   return kTRUE;
 493 }
 494 //____________________________________________________________________
 495 Bool_t
 496 AliFMDCorrELossFit::SetFit(UShort_t  d, Char_t r, Double_t eta,
 497                            Int_t quality, const TF1& f)
 498 {
 499   ELossFit* e = new ELossFit(quality, f);
 500   if (!SetFit(d, r, eta, e)) {
 501     delete e;
 502     return kFALSE;
 503   }
 504   return kTRUE;
 505 }
 506 //____________________________________________________________________
 507 AliFMDCorrELossFit::ELossFit*
 508 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Int_t etabin) const
 509 {
 510   TObjArray* ringArray = GetRingArray(d, r);
 511   if (!ringArray) {
 512     AliError(Form("Failed to make ring array for FMD%d%c", d, r));
 513     return 0;
 514   }
 515   if (etabin <= 0 || etabin >= fEtaAxis.GetNbins()) {
 516     AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
 517                   etabin, 1, fEtaAxis.GetNbins(), d, r));
 518     return 0;
 519   }
 520   if (etabin > ringArray->GetEntriesFast()) {
 521     AliError(Form("Eta bin=%3d out of bounds [%d,%d] for FMD%d%c",
 522                   etabin, 1, ringArray->GetEntriesFast(), d, r));
 523     return 0;
 524   }
 525   else if (etabin >= ringArray->GetEntriesFast()) {
 526     AliWarning(Form("Eta bin=%3d out of bounds by +1 [%d,%d] for FMD%d%c, "
 527                     "trying %3d", etabin, 1, ringArray->GetEntriesFast(), d, r,
 528                     etabin-1));
 529     etabin--;
 530   }
 531   else if (!ringArray->At(etabin)) {
 532     AliWarning(Form("Eta bin=%d has no fit for FMD%d%c, trying %03d",
 533                     etabin, d, r, etabin+1));
 534     etabin++;
 535   }
 536   return static_cast<ELossFit*>(ringArray->At(etabin));
 537 }
 538 //____________________________________________________________________
 539 AliFMDCorrELossFit::ELossFit*
 540 AliFMDCorrELossFit::FindFit(UShort_t  d, Char_t r, Double_t eta) const
 541 {
 542   Int_t etabin = FindEtaBin(eta);
 543   return FindFit(d, r, etabin);
 544 }
 545 //____________________________________________________________________
 546 TObjArray*
 547 AliFMDCorrELossFit::GetRingArray(UShort_t d, Char_t r) const
 548 {
 549   Int_t idx = -1;
 550   switch (d) {
 551   case 1:   idx = 0; break;
 552   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
 553   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
 554   }
 555   if (idx < 0 || idx >= fRings.GetEntriesFast()) return 0;
 556   return static_cast<TObjArray*>(fRings.At(idx));
 557 }
 558 //____________________________________________________________________
 559 TObjArray*
 560 AliFMDCorrELossFit::GetOrMakeRingArray(UShort_t d, Char_t r)
 561 {
 562   Int_t idx = -1;
 563   switch (d) {
 564   case 1:   idx = 0; break;
 565   case 2:   idx = (r == 'i' || r == 'I') ? 1 : 2; break;
 566   case 3:   idx = (r == 'o' || r == 'I') ? 3 : 4; break;
 567   }
 568   if (idx < 0) return 0;
 569   if (idx >= fRings.GetEntriesFast() || !fRings.At(idx)) {
 570     TObjArray* a = new TObjArray(0);
 571     // TOrdCollection* a = new TOrdCollection(fEtaAxis.GetNbins());
 572     a->SetName(Form("FMD%d%c", d, r));
 573     a->SetOwner();
 574     fRings.AddAtAndExpand(a, idx);
 575   }
 576   return static_cast<TObjArray*>(fRings.At(idx));
 577 }
 578
 579 namespace {
 580   TH1D* MakeHist(const TAxis& axis, const char* name, const char* title,
 581                  Int_t color)
 582   {
 583     TH1D* h = new TH1D(Form("%s_%s", name, title),
 584                        Form("%s %s", name, title),
 585                        axis.GetNbins(), axis.GetXmin(), axis.GetXmax());
 586     h->SetDirectory(0);
 587     h->SetMarkerStyle(20);
 588     h->SetMarkerColor(color);
 589     h->SetMarkerSize(0.5);
 590     h->SetFillColor(color);
 591     h->SetFillStyle(3001);
 592     h->SetLineColor(color);
 593     return h;
 594   }
 595 }
 596
 597 //____________________________________________________________________
 598 void
 599 AliFMDCorrELossFit::Draw(Option_t* option)
 600 {
 601   TString opt(Form("nostack %s", option));
 602   opt.ToLower();
 603   Bool_t  rel = (opt.Contains("rel"));
 604   Bool_t  err = (opt.Contains("err"));
 605   if (rel) opt.ReplaceAll("rel","");
 606   if (err) opt.ReplaceAll("err","");
 607   Int_t    nRings = fRings.GetEntriesFast();
 608   UShort_t maxN   = 0;
 609   for (Int_t i = 0; i < nRings; i++) {
 610     if (!fRings.At(i)) continue;
 611     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
 612     Int_t      nFits = a->GetEntriesFast();
 613
 614     for (Int_t j = 0; j < nFits; j++) {
 615       ELossFit* fit = static_cast<ELossFit*>(a->At(j));
 616       if (!fit) continue;
 617       maxN          = TMath::Max(maxN, UShort_t(fit->fN));
 618     }
 619   }
 620   // AliInfo(Form("Maximum N is %d", maxN));
 621   Int_t nPad = 7+maxN-1; // 7 regular params, and maxN-1 weights
 622   TVirtualPad* pad = gPad;
 623   pad->Divide(2, (nPad+1)/2, 0.1, 0, 0);
 624
 625   THStack* chi2nu;
 626   THStack* c;
 627   THStack* delta;
 628   THStack* xi;
 629   THStack* sigma;
 630   THStack* sigman;
 631   THStack* n;
 632   TList stacks;
 633   stacks.AddAt(chi2nu= new THStack("chi2",   "#chi^{2}/#nu"), 0);
 634   stacks.AddAt(c     = new THStack("c",       "C"),           1);
 635   stacks.AddAt(delta = new THStack("delta",  "#Delta_{mp}"),  2);
 636   stacks.AddAt(xi    = new THStack("xi",     "#xi"),          3);
 637   stacks.AddAt(sigma = new THStack("sigma",  "#sigma"),       4);
 638   stacks.AddAt(sigman= new THStack("sigman", "#sigma_{n}"),   5);
 639   stacks.AddAt(n     = new THStack("n",      "N"),            6);
 640   for (Int_t i = 1; i <= maxN; i++) {
 641     stacks.AddAt(new THStack(Form("a_%02d", i+1), Form("a_{%d}", i+1)), 6+i);
 642   }
 643
 644   for (Int_t i = 0; i < nRings; i++) {
 645     if (!fRings.At(i)) continue;
 646     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
 647     Int_t      nFits = a->GetEntriesFast();
 648     Int_t      color = 0;
 649     switch (i) {
 650     case 0: color = kRed+2;   break;
 651     case 1: color = kGreen+2; break;
 652     case 2: color = kGreen-2; break;
 653     case 3: color = kBlue+2;  break;
 654     case 4: color = kBlue-2;  break;
 655     }
 656
 657     TH1D* hChi    = MakeHist(fEtaAxis,a->GetName(), "chi2",   color);
 658     TH1D* hC      = MakeHist(fEtaAxis,a->GetName(), "c",      color);
 659     TH1D* hDelta  = MakeHist(fEtaAxis,a->GetName(), "delta",  color);
 660     TH1D* hXi     = MakeHist(fEtaAxis,a->GetName(), "xi",     color);
 661     TH1D* hSigma  = MakeHist(fEtaAxis,a->GetName(), "sigma",  color);
 662     TH1D* hSigmaN = MakeHist(fEtaAxis,a->GetName(), "sigman", color);
 663     TH1D* hN      = MakeHist(fEtaAxis,a->GetName(), "n",      color);
 664     TH1D* hA[maxN];
 665     const char* ho = (rel || !err ? "hist" : "e");
 666     chi2nu->Add(hChi,   "hist"); // 0
 667     c     ->Add(hC,     ho);     // 1
 668     delta ->Add(hDelta, ho);     // 2
 669     xi    ->Add(hXi,    ho);     // 3
 670     sigma ->Add(hSigma, ho);     // 4
 671     sigman->Add(hSigmaN,ho);     // 5
 672     n     ->Add(hN,     "hist"); // 6
 673     hChi->SetFillColor(color);
 674     hChi->SetFillStyle(3001);
 675     hN->SetFillColor(color);
 676     hN->SetFillStyle(3001);
 677
 678     for (Int_t k = 1; k <= maxN; k++) {
 679       hA[k-1] = MakeHist(fEtaAxis,a->GetName(), Form("a%02d",k+1), color);
 680       static_cast<THStack*>(stacks.At(6+k))->Add(hA[k-1], ho);
 681     }
 682
 683     for (Int_t j = 0; j < nFits; j++) {
 684       ELossFit* f = static_cast<ELossFit*>(a->At(j));
 685       if (!f) continue;
 686
 687       Int_t     b = f->fBin;
 688       Int_t     nW = f->FindMaxWeight();
 689       hChi   ->SetBinContent(b, (f->fNu <= 0 ? 0 : f->fChi2 / f->fNu));
 690       hN     ->SetBinContent(b, nW);
 691
 692       if (rel) {
 693         hC     ->SetBinContent(b, f->fC     >0 ?f->fEC     /f->fC      :0);
 694         hDelta ->SetBinContent(b, f->fDelta >0 ?f->fEDelta /f->fDelta  :0);
 695         hXi    ->SetBinContent(b, f->fXi    >0 ?f->fEXi    /f->fXi     :0);
 696         hSigma ->SetBinContent(b, f->fSigma >0 ?f->fESigma /f->fSigma  :0);
 697         hSigmaN->SetBinContent(b, f->fSigmaN>0 ?f->fESigmaN/f->fSigmaN :0);
 698       }
 699       else {
 700         hC     ->SetBinContent(b, f->fC);
 701         hDelta ->SetBinContent(b, f->fDelta);
 702         hXi    ->SetBinContent(b, f->fXi);
 703         hSigma ->SetBinContent(b, f->fSigma);
 704         hSigmaN->SetBinContent(b, f->fSigmaN);
 705         hC     ->SetBinError(b, f->fEC);
 706         hDelta ->SetBinError(b, f->fEDelta);
 707         hXi    ->SetBinError(b, f->fEXi);
 708         hSigma ->SetBinError(b, f->fESigma);
 709         hSigmaN->SetBinError(b, f->fESigmaN);
 710       }
 711       for (Int_t k = 0; k < f->fN-1 && k < maxN; k++) {
 712         if (rel)
 713           hA[k]->SetBinContent(b, f->fA[k] > 0 ? f->fEA[k] / f->fA[k] : 0);
 714         else {
 715           hA[k]->SetBinContent(b, f->fA[k]);
 716           hA[k]->SetBinError(b, f->fEA[k]);
 717         }
 718       }
 719     }
 720   }
 721   Int_t nPad2 = (nPad+1) / 2;
 722   for (Int_t i = 0; i < nPad; i++) {
 723     Int_t iPad = 1 + i/nPad2 + 2 * (i % nPad2);
 724     TVirtualPad* p = pad->cd(iPad);
 725     p->SetLeftMargin(.15);
 726     p->SetFillColor(0);
 727     p->SetFillStyle(0);
 728     p->SetGridx();
 729     p->SetGridy();
 730     if (rel && i != 0 && i != 6 && i != 5 && i != 4) p->SetLogy();
 731
 732     THStack* stack = static_cast<THStack*>(stacks.At(i));
 733     // AliInfo(Form("Drawing %s (%d) in pad %d", stack->GetName(), i, iPad));
 734     stack->Draw(opt.Data());
 735
 736     TString tit(stack->GetTitle());
 737     if (rel && i != 0 && i != 6)
 738       tit = Form("#delta %s/%s", tit.Data(), tit.Data());
 739     TH1*   hist  = stack->GetHistogram();
 740     TAxis* yaxis = hist->GetYaxis();
 741     yaxis->SetTitle(tit.Data());
 742     yaxis->SetTitleSize(0.15);
 743     yaxis->SetLabelSize(0.08);
 744     yaxis->SetTitleOffset(0.35);
 745     yaxis->SetNdivisions(5);
 746
 747     TAxis* xaxis = stack->GetHistogram()->GetXaxis();
 748     xaxis->SetTitle("#eta");
 749     xaxis->SetTitleSize(0.15);
 750     xaxis->SetLabelSize(0.08);
 751     xaxis->SetTitleOffset(0.35);
 752     xaxis->SetNdivisions(10);
 753
 754
 755     if (!rel) {
 756       switch (i) {
 757       case 0:                         break; // chi^2/nu
 758       case 1:                         break; // C
 759       case 2: stack->SetMinimum(0.4); break; // Delta
 760       case 3: stack->SetMinimum(0.02);break; // xi
 761       case 4: stack->SetMinimum(0.05);break; // sigma
 762       case 5:                         break; // sigmaN
 763       case 6:
 764         stack->SetMinimum(-.1);
 765         stack->SetMaximum(maxN+.5);
 766         break;                              // N
 767       default:                       break; // Weights
 768       }
 769     }
 770     stack->DrawClone(opt.Data());
 771   }
 772   pad->cd();
 773 }
 774
 775 //____________________________________________________________________
 776 void
 777 AliFMDCorrELossFit::Print(Option_t* option) const
 778 {
 779   TString opt(option);
 780   opt.ToUpper();
 781   Int_t nRings = fRings.GetEntriesFast();
 782   bool recurse = opt.Contains("R");
 783
 784   std::cout << "Low cut in fit range: " << fLowCut << "\n"
 785             << "Eta axis:             " << fEtaAxis.GetNbins()
 786             << " bins, range [" << fEtaAxis.GetXmin() << ","
 787             << fEtaAxis.GetXmax() << "]" << std::endl;
 788
 789   for (Int_t i = 0; i < nRings; i++) {
 790     if (!fRings.At(i)) continue;
 791     TObjArray* a     = static_cast<TObjArray*>(fRings.At(i));
 792     Int_t      nFits = a->GetEntriesFast();
 793
 794     std::cout << a->GetName() << " [" << nFits << " entries]"
 795               << (recurse ? ":\n" : "\t");
 796     Int_t min = fEtaAxis.GetNbins()+1;
 797     Int_t max = 0;
 798     for (Int_t j = 0; j < nFits; j++) {
 799       if (!a->At(j)) continue;
 800
 801       min = TMath::Min(j, min);
 802       max = TMath::Max(j, max);
 803
 804       if (recurse) {
 805         std::cout << "Bin # " << j << "\t";
 806         ELossFit* fit = static_cast<ELossFit*>(a->At(j));
 807         fit->Print(option);
 808       }
 809     }
 810     if (!recurse)
 811       std::cout << " bin range: " << std::setw(3) << min
 812                 << "-" << std::setw(3) << max << " " << std::setw(3)
 813                 << (max-min+1) << " bins" << std::endl;
 814   }
 815 }
 816
 817 //____________________________________________________________________
 818 void
 819 AliFMDCorrELossFit::Browse(TBrowser* b)
 820 {
 821   b->Add(&fRings);
 822   b->Add(&fEtaAxis);
 823 }
 824
 825
 826
 827 //____________________________________________________________________
 828 //
 829 // EOF
 830 //