TestF.C

Go to the documentation of this file.

#include "AliLandauGaus.h"
#include <TSystem.h>
#include <TString.h>
#include <TList.h>
#include <TFile.h>
#include <TMath.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TFitResult.h>
#include <TNtuple.h>
#include <TAxis.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TLine.h>
#include <TLatex.h>
#include <TPaveStats.h>

//====================================================================
struct TestF
{
  //__________________________________________________________________
  struct Test { 
    TF1* src;
    Test() : src(0) {}
    Test(const Test& o) 
      : src(static_cast<TF1*>(o.src ? o.src->Clone() : 0)) {}
    virtual ~Test() { delete src; }
    Test& operator=(const Test& o) {
      if (&o == this) return *this;
      if (o.src) src = static_cast<TF1*>(o.src->Clone());
      return *this;
    }
    void Run(const Double_t  c, 
       const Double_t  delta, 
       const Double_t  xi, 
       const Double_t  sigma, 
       const Int_t     n, 
       const Double_t* a,
       const Double_t  xMin, 
       const Double_t  xMax) {
      printf("%8.4f | %8.4f | %8.4f | ", delta, xi, sigma);
      fflush(stdout);
      src = AliLandauGaus::MakeFn(c,delta,xi,sigma,0,n,&(a[1]),xMin,xMax);
      src->SetName("src");
      src->SetLineColor(kBlack);
      src->SetNpx(10000);
      // If the constant is 1, then this is a fake distribtion and we
      // need to scale it. 
      if (c >= 1) { 
  Double_t max = src->GetMaximum();
  src->SetParameter(0, c/max);
      }

      // src->Print();
      DoRun();
    }
    Double_t C()      const { return src->GetParameter(AliLandauGaus::kC); }
    Double_t Delta()  const { return src->GetParameter(AliLandauGaus::kDelta); }
    Double_t Xi()     const { return src->GetParameter(AliLandauGaus::kXi); }
    Double_t Sigma()  const { return src->GetParameter(AliLandauGaus::kSigma); }
    Double_t SigmaN() const { return src->GetParameter(AliLandauGaus::kSigmaN);}
    Double_t N()      const { return src->GetParameter(AliLandauGaus::kN); }
    Double_t A(Int_t i) const { 
      return (i <= 1 ? 1 : src->GetParameter(AliLandauGaus::kN+(i-1))); }
    Double_t XMin() const { return src->GetXmin(); }
    Double_t XMax() const { return src->GetXmax(); }
    TF1* Fi(Double_t c, Double_t delta, Double_t xi, Double_t sigma, 
      Double_t sigmaN, Int_t i) const { 
      return AliLandauGaus::MakeFi(c,delta,xi,sigma,sigmaN,i,XMin(),XMax());
    }
    TF1* Fi(Int_t i) { 
      return Fi(C()*A(i),Delta(),Xi(),Sigma(),SigmaN(),i);
    }
    
    virtual void DoRun() = 0;
    virtual void WriteOut(TDirectory* d) = 0;
    virtual void DrawInPad(TVirtualPad* p) = 0;
    /* @} */
  };
  //__________________________________________________________________
  static TLine*
  VerticalLine(Int_t i, Color_t col, Double_t x, Double_t low, Double_t high)
  {
    TLine* l = new TLine(x,low,x,high);
    l->SetLineStyle(i);
    l->SetLineColor(col);
    return l;
  }
  //__________________________________________________________________
  virtual ~TestF() {}
  virtual Test* MakeTest() = 0;
  Test* DoOne(const Double_t  c, 
        const Double_t  delta, 
        const Double_t  xi, 
        const Double_t  sigma, 
        const Int_t     n, 
        const Double_t* a, 
        const Double_t  xMin, 
        const Double_t  xMax,
        TVirtualPad*    pad=0)
  {
    Test* t = MakeTest();
    t->Run(c, delta, xi, sigma, n, a, xMin, xMax);

    TString name(Form("sigma%04d_xi%04d", Int_t(100*sigma), Int_t(100*xi)));
    TString title(Form("Sigma=%f xi=%f", sigma, xi));
    if (!pad) { 
      TCanvas* can = new TCanvas(name, title, 1200, 900);
      pad = can;
    }
    else 
      pad->Clear();
    t->DrawInPad(pad);

    return t;
  }
  Test* DoUnit(const Double_t sigma, 
         const Double_t xi=1, 
         const Int_t    n=10,
         TVirtualPad*   pad=0) 
  { 
    const Double_t c      = 1;
    const Double_t delta  = 10;
    Double_t       idelta = delta;
    Double_t       ixi    = xi;
    Double_t       isigma = sigma;
    AliLandauGaus::IPars(n+1,idelta,ixi, isigma);
    Double_t       up     = idelta + ixi + isigma; // 8*n
    Double_t       a[n];
    a[0] = 1;
    for (Int_t i = 1; i < n; i++) 
      a[i] = TMath::Exp(0.5 - TMath::Sqrt(2)*(i+1));

    return DoOne(c,delta,xi,sigma,n,a,-3,up, pad);
  }
  Test* DoRealistic(TVirtualPad* pad=0) 
  {
    Double_t c     = 0.4070;
    Double_t delta = 0.5651;
    Double_t xi    = 0.0643;
    Double_t sigma = 0.0611;
    Int_t    n     = 5;
    Double_t a[]   = {1, 0.1179, 0.0333, 0.0068, 0.0012 };

    return DoOne(c, delta, xi, sigma, n, a, 0.02, n, pad);
  }
  virtual const char* Prefix() const = 0;
  virtual TNtuple* MakeNTuple() { return 0; }
  virtual void SetupScan(UShort_t  mode, 
       TFile*&   file, 
       TCanvas*& canvas, 
       TNtuple*& nt) 
  {
    TString  base = Form("%s%s", Prefix(), 
       (mode == 0 ? "SigmaXi" : 
        mode == 1 ? "Sigma" : "Xi"));
    file   = TFile::Open(Form("%s.root", base.Data()), "RECREATE");
    nt     = MakeNTuple();
    canvas = new TCanvas(base, Form("%s.pdf", base.Data()),
             TMath::Sqrt(2)*900, 900);
    TString title("pdf Landscape Title:Scan");
    gSystem->RedirectOutput("/dev/null");
    canvas->Print(Form("%s[", canvas->GetTitle()), title);
    gSystem->RedirectOutput(0);
  }
  void PrintCanvas(TCanvas* c, Double_t xi, Double_t sigma) {
    TString tit = "pdf Landscape Title:";
    Bool_t  cls = false;
    if (xi > 0 && sigma > 0) {
      tit.Append(Form("xi=%8.4f sigma=%8.4f", xi, sigma));
    }
    else { 
      tit.Append("Results");
      cls = true;
    }
    gSystem->RedirectOutput("/dev/null");
    c->Print(c->GetTitle(), tit);
    if (cls) {
      c->Clear();
      c->Print(Form("%s]", c->GetTitle()), tit);
    }
    gSystem->RedirectOutput(0);
    if (cls)
      ::Info("", "Plots saved in %s", c->GetTitle());
  }
  virtual void PreLoop(UShort_t mode, Int_t nVal) = 0;
  virtual void Step(UShort_t mode, 
        Int_t    i, 
        Int_t    j, 
        Int_t    n, 
        Test*    test,
        TNtuple* nt) = 0;
  virtual void PostLoop(UShort_t mode, 
      TCanvas* can, 
      TNtuple* nt, 
      TFile*   out) = 0;
  void ScanOne(Bool_t          scanSigma, 
         Int_t           n, 
         const Double_t* values, 
         Int_t           maxN=10) 
  {
    UShort_t mod = scanSigma ? 1 : 2;
    Scan(mod, n, values, maxN);
  }
  void ScanOne(Bool_t         scanSigma, 
         const TArrayD& values, 
         Int_t          maxN)
  {
    ScanOne(scanSigma, values.GetSize(), values.GetArray(), maxN);
  }
  void ScanTwo(Int_t           n, 
         const Double_t* values, 
         Int_t           maxN=10)
  {
    Scan(0, n, values, maxN);
  }
  void ScanTwo(const TArrayD& values, Int_t maxN)
  {
    ScanTwo(values.GetSize(), values.GetArray(), maxN);
  }
  void Scan(UShort_t        mode, 
      Int_t           n, 
      const Double_t* values, 
      Int_t           maxN)
  {
    TFile*   out = 0;
    TNtuple* nt  = 0;
    TCanvas* can = 0;
    SetupScan(mode, out, can, nt);
    PreLoop(mode, n);

    Test* t  = DoRealistic(can);
    t->WriteOut(out->mkdir("realistic"));
    can->Write();
    PrintCanvas(can, t->Xi(), t->Sigma());
    
    if (mode == 0) LoopTwo(n, values, maxN, can, nt, out);
    else           LoopOne(mode, n, values, maxN, can, nt, out);
    
    PostLoop(mode, can, nt, out);
    out->Write();
    out->Close();

  }
  void LoopOne(UShort_t        mode, 
         Int_t           n, 
         const Double_t* values, 
         Int_t           maxN,
         TCanvas*        can, 
         TNtuple*        nt, 
         TFile*          out) {
    Bool_t      scanSigma = mode == 1;
    const char* var       = (scanSigma ? "sigma" : "xi");
  
    for (Int_t i = 0; i < n; i++) { 
      Double_t    v     = values[i];
      Double_t    sigma = (scanSigma  ? v : 1);
      Double_t    xi    = (!scanSigma ? v : 1);

      StepIt(mode, xi, sigma, maxN, 0, i, n, can, nt, out, var, v);
    }
  }
  void LoopTwo(Int_t           n, 
         const Double_t* values, 
         Int_t           maxN,
         TCanvas*        can, 
         TNtuple*        nt, 
         TFile*          out) 
  { 
    UShort_t mode = 0;
    for (Int_t i = 0; i < n; i++) { 
      Double_t      xi  = values[i];
      TDirectory*   d   = out->mkdir(Form("xi%04d", Int_t(100*xi)));
      for (Int_t j = 0; j < n; j++) { 
  Double_t    sigma = values[j];
  StepIt(mode, xi, sigma, maxN, i, j, n, can, nt, d, "sigma", sigma);
      }
    }
  }
  void StepIt(UShort_t    mode, 
        Double_t    xi, 
        Double_t    sigma, 
        Int_t       maxN, 
        Int_t       i, 
        Int_t       j, 
        Int_t       n, 
        TCanvas*    can, 
        TNtuple*    nt, 
        TDirectory* p, 
        const char* pre, 
        Double_t    v)
  {
    TDirectory* d     = p->mkdir(Form("%s%04d", pre, Int_t(100*v)));
    Test*       t     = DoUnit(sigma, xi, maxN, can);
    t->WriteOut(d);
    can->Write();
    PrintCanvas(can, xi, sigma);
      
    Step(mode, i, j, n, t, nt);
  }
};

//====================================================================
#ifdef TEST_SHIFT
#include "WrappedGraph.C"
#include <TGraphErrors.h>
#include <TGraph2DErrors.h>
#include <TMultiGraph.h>
#include <Math/RootFinder.h>

struct TestShift : public TestF
{
  static Double_t fitFunc(Double_t* xp, Double_t* pp) {
    Double_t x     = *xp;
    Double_t c     = pp[0];
    Double_t p     = pp[1];
    Double_t xi    = pp[2];
    Double_t sigma = pp[3];
    Double_t u     = sigma / xi;
    Double_t a     = c * u * sigma;
    Double_t b     = p * u * TMath::Sqrt(u);
      
    return a / TMath::Power(1+1./x, b);
  }
  static void ScaleGraph(TGraph* g, Double_t scale)
  {
    // Info("", "Scaling graph %p with %f", g, scale);
    if (scale == 1) return;
    for (Int_t i = 0; i < g->GetN(); i++) 
      g->SetPoint(i,g->GetX()[i], g->GetY()[i] * scale);
  }
  //__________________________________________________________________
  struct Test : public TestF::Test 
  {
    TGraph*  zeros;       // Zeros of the derivatives
    TGraph*  diff;        // Difference from calculated Delta to roots
    TList    comps;       // Components 
    TList    derivs;      // Derivatives 
    TF1*     fit;         // Fit to difference 
    Test() : TestF::Test(), zeros(0), diff(0), fit(0) 
    {
      comps.SetName("components");
      derivs.SetName("derivatives");

    }
    void DoRun() 
    { 
      AliLandauGaus::EnableSigmaShift(false); // Disable sigma-shift
      const Int_t n = N();
      zeros = new TGraph(n);
      zeros->SetName("zeros");
      zeros->SetTitle("Zeros of derivatives");
      zeros->SetMarkerStyle(25);
      zeros->SetMarkerSize(1.2);
      zeros->SetMarkerColor(kRed+1);

      diff = new TGraph(n);
      diff->SetName("diff");
      diff->SetTitle("#delta#Delta_{p}");
      diff->SetMarkerStyle(20);
      diff->SetMarkerColor(kBlue+1);

      // Make the components 
      Int_t j = 0;
      for (Int_t i = 1; i <= n; i++, j++) 
  if (!MakeComp(i)) break;
      zeros->Set(j);
      diff->Set(j);

      // Now fit 
      fit = new TF1("fit", &fitFunc, 1, j, 4);
      fit->SetParNames("c", "p", "#xi", "#sigma");
      fit->SetParameter(0, .5);
      fit->SetParameter(1, .5);
      fit->FixParameter(2, Xi());
      fit->FixParameter(3, Sigma());
      
      // fit->SetParNames("c", "#xi", "#sigma");
      // fit->FixParameter(1, xi);


      Int_t ifit = 4;
      do { 
  gSystem->RedirectOutput("/dev/null");
  diff->Fit(fit, Form("MR%s", ifit == 0 ? "Q" : "Q"));
  gSystem->RedirectOutput(0);
  ifit--;
      } while (ifit >= 0);
      Printf("%8.4f | %8.4f %9.5f | %8.4f %9.5f",
       fit->GetChisquare() / fit->GetNDF(),
       fit->GetParameter(0), fit->GetParError(0),
       fit->GetParameter(1), fit->GetParError(1));

    }
    Bool_t MakeComp(const Int_t    i) {
      TF1* fi = Fi(i);
      fi->SetNpx(src->GetNpx());
      // fi->Print();

      TGraph* dF = new TGraph(fi, "d");
      dF->SetName(Form("d%s", fi->GetName()));
      dF->SetTitle(Form("d%s/d#Delta", fi->GetTitle()));
      dF->SetMarkerStyle(1); // A dot 
      // dF->SetMarkerStyle(20+i-1);
      dF->SetMarkerColor(fi->GetLineColor());
      dF->SetLineColor(fi->GetLineColor());

      Double_t max = TMath::MaxElement(dF->GetN(), dF->GetY());
      if (max < 1e-6) return false;

      ScaleGraph(dF, 1./max);
      Double_t delta  = Delta();
      Double_t xi     = Xi();
      Double_t deltaI = i * (delta + xi * TMath::Log(i));
      Double_t xiI    = i * xi;
      Double_t xR     = FindRoot(dF, deltaI-2*xiI, deltaI+2*xiI);
      
      // Printf("Component %2d: Exp: %8.4f  Got: %8.4f  Diff: %9.5f", 
      //        i, deltaI, xR, (xR-deltaI));
      // printf(".");
      comps.Add(fi);
      derivs.Add(dF);
      zeros->SetPoint(i-1, xR, dF->Eval(xR));
      diff->SetPoint(i-1, i, xR - deltaI);

      return true;
    }
    Double_t FindRoot(TGraph* g, Double_t xMin, Double_t xMax) { 
      WrappedGraph* wG = new WrappedGraph(g);
      ROOT::Math::RootFinder rfb(ROOT::Math::RootFinder::kBRENT);
      rfb.SetFunction(*wG, xMin, xMax); 
      rfb.Solve();
      Double_t xR = rfb.Root();
      return xR;
    }
    void DrawInPad(TVirtualPad* body) { 
      body->SetTopMargin(0.01);
      body->SetRightMargin(0.01);
      body->Divide(2,1);
      
      TVirtualPad* p = body->cd(1);
      p->SetTopMargin(0.01);
      p->SetRightMargin(0.01);
      p->Divide(1,2,0,0);
      
      TVirtualPad* q = p->cd(1);
      q->SetRightMargin(0.01);
      q->SetGridx();
      q->SetLogy();
      src->Draw();
      Double_t fmax = src->GetHistogram()->GetMaximum();
      Double_t fmin = src->GetHistogram()->GetMinimum();
      fmin          = TMath::Max(fmin, fmax/1e8);
      src->GetHistogram()->SetMinimum(fmin);
      TLatex*  ltx  = new TLatex(q->GetLeftMargin()+.02, 
         1-q->GetTopMargin()-.02,
         Form("#Delta_{p}=%f, #xi=%f, #sigma=%f", 
              Delta(), Xi(), Sigma()));
      // Printf("%s: Max=%f scale=%f", src->GetName(), fmax, 1/fmax);
      ltx->SetTextFont(42);
      ltx->SetNDC();
      ltx->SetTextAlign(13);
      ltx->Draw();

      q = p->cd(2);
      q->SetRightMargin(0.01);
      q->SetGridx();
      
      Double_t x = 1-q->GetRightMargin()-.02;
      Double_t y = 1-q->GetTopMargin()-.1;
      for (Int_t i = 1; i <= comps.GetEntries(); i++) { 
  p->cd(1);
  TF1* fi = static_cast<TF1*>(comps.At(i-1));
  fi->Draw("same");

  Double_t deltaR = zeros->GetX()[i-1];
  Double_t deltaI = i * (Delta() + Xi() * TMath::Log(i));

  Color_t col = fi->GetLineColor();
  VerticalLine(1,col,deltaI,fmin,fmax)->Draw();
  VerticalLine(3,col,deltaR,fmin,fmax)->Draw();
  
  ltx = new TLatex(x,y,Form("#Delta_{p,%2d}=%8.4f %8.4f",
          i,deltaI,deltaR));
  ltx->SetTextAlign(33);
  ltx->SetTextFont(42);
  ltx->SetNDC();
  ltx->Draw();
  y -= ltx->GetTextSize() + .01;

  p->cd(2);

  TGraph* df = static_cast<TGraph*>(derivs.At(i-1));
  df->Draw(i == 1 ? "alp" : "lp");
  df->GetHistogram()->GetXaxis()->SetRangeUser(XMin(), XMax()); 
  VerticalLine(1,col,deltaI,-.6,1.1)->Draw();
  VerticalLine(3,col,deltaR,-.6,1.1)->Draw();
      }
      p->cd(2);
      zeros->Draw("p");
      
      gStyle->SetOptFit(111111);
      gStyle->SetStatY(.6);
      p = body->cd(2);
      p->SetRightMargin(0.02);
      p->SetTopMargin(0.02);
      diff->Draw("alp");
      diff->GetHistogram()->GetListOfFunctions()->Add(fit);
      diff->GetHistogram()->SetXTitle("N_{particle}");
      diff->GetHistogram()->SetYTitle("#Delta_{r}-#Delta_{p}");
      p->Clear();
      diff->Draw("alp");

      body->Modified();
      body->Update();
      body->cd();
    }
    void WriteOut(TDirectory* dir) { 
      dir->cd();
      src->Write();
      zeros->Write();
      diff->Write();
      comps.Write(comps.GetName(),  TObject::kSingleKey);
      derivs.Write(derivs.GetName(),  TObject::kSingleKey);
      fit->Write();
    }
  };
  //__________________________________________________________________
  TGraphErrors*   cs;
  TGraphErrors*   ps;
  TGraph2DErrors* c2s;
  TGraph2DErrors* p2s;
  TMultiGraph*    mcs;
  TMultiGraph*    mps;
  
  virtual TestF::Test* MakeTest() { return new Test(); }
  virtual const char* Prefix() const { return "shift"; }
  virtual TNtuple* MakeNTuple () 
  { 
    return new TNtuple("shift","Shift scan",
           "Delta:xi:sigma:chi2nu:p:ep:c:ec");
  }
  virtual void PreLoop(UShort_t mode, Int_t n) 
  { 
    if (mode == 0) {
      c2s = new TGraph2DErrors(n*n);c2s->SetName("cs");
      p2s = new TGraph2DErrors(n*n);p2s->SetName("ps");
      mcs = new TMultiGraph("mcs", "C as function of #sigma/#xi");
      mps = new TMultiGraph("mps", "P as function of #sigma/#xi");
      c2s->SetDirectory(0);
      p2s->SetDirectory(0);
      c2s->SetMarkerStyle(20);
      p2s->SetMarkerStyle(20);
      c2s->SetTitle("c(#xi,#sigma)");
      p2s->SetTitle("p(#xi,#sigma)");
    }
    else {
      const char* var = (mode == 1 ? "#sigma" : "#xi");
      cs = new TGraphErrors(n); cs->SetName("cs");
      ps = new TGraphErrors(n); ps->SetName("ps");
      cs->SetMarkerStyle(20);
      ps->SetMarkerStyle(20);
      cs->SetTitle(Form("c(%s)", var));
      ps->SetTitle(Form("p(%s)", var));
    }
    Printf("%-8s | %-8s | %-8s | %-8s | %-8s %-9s | %-8s %-9s",
     "Delta_p", "xi", "sigma", "chi^2/nu", "c", "error", "p", "error");
  }
  void FillMG(TMultiGraph* mg,
        Int_t        i, 
        Int_t        j, 
        Int_t        n, 
        Double_t     u, 
        Double_t     y, 
        Double_t     e) 
  { 
    TList* l = mg->GetListOfGraphs();
    TGraphErrors* ge = (l ? static_cast<TGraphErrors*>(l->At(i)) : 0);
    if (!ge) { 
      ge = new TGraphErrors(n);
      ge->SetName(Form("%s%02d", mg->GetName(), i));
      ge->SetTitle(Form(Form("%d", i)));
      ge->SetMarkerColor(AliLandauGaus::GetIColor(i));
      ge->SetLineColor(AliLandauGaus::GetIColor(i));
      ge->SetMarkerStyle(20);
      mg->Add(ge);
    }
    ge->SetPoint(j, u, y);
    ge->SetPointError(j, 0, e);
  }
  virtual void Step(UShort_t     mode, 
        Int_t        i, 
        Int_t        j, 
        Int_t        n, 
        TestF::Test* bt,
        TNtuple*     nt) 
  { 
    Test* t = static_cast<Test*>(bt);
    Double_t xi    = t->Xi();
    Double_t sigma = t->Sigma();
    Double_t c     = t->fit->GetParameter(0);
    Double_t p     = t->fit->GetParameter(1);
    Double_t ec    = t->fit->GetParError(0);
    Double_t ep    = t->fit->GetParError(1);
    Int_t    idx   = i * n + j;
    if (mode == 0) {
      c2s->SetPoint(idx, xi, sigma, c);
      p2s->SetPoint(idx, xi, sigma, p);
      c2s->SetPointError(idx, 0, 0, ec);
      p2s->SetPointError(idx, 0, 0, ep);
      FillMG(mcs, i, j, n, sigma/xi, c, ec);
      FillMG(mps, i, j, n, sigma/xi, p, ep);
    }
    else { 
      cs->SetPoint(idx, (mode==1 ? sigma : xi), c);
      ps->SetPoint(idx, (mode==1 ? sigma : xi), p);
      cs->SetPointError(idx, 0, ec);
      ps->SetPointError(idx, 0, ep);
    }
    nt->Fill(1,xi,sigma,t->fit->GetChisquare()/t->fit->GetNDF(),c,ec,p,ep);
  }
  void PostLoop(UShort_t mode, TCanvas* can, TNtuple* /*nt*/, TFile* out) 
  { 
    can->Clear();
    if (mode == 0) { 
      can->Divide(2,2);
      can->cd(1); c2s->Draw("TRI2Z");
      can->cd(2); p2s->Draw("TRI2Z");
      can->cd(3); mcs->Draw("APL");
      can->cd(4); mps->Draw("APL");
      c2s->GetHistogram()->SetXTitle("#xi");
      p2s->GetHistogram()->SetXTitle("#xi");
      c2s->GetHistogram()->SetYTitle("#sigma");
      p2s->GetHistogram()->SetYTitle("#sigma");
      c2s->GetHistogram()->SetZTitle("c");
      p2s->GetHistogram()->SetZTitle("p");
      mcs->GetHistogram()->SetXTitle("#sigma/#xi");
      mps->GetHistogram()->SetXTitle("#sigma/#xi");
      mcs->GetHistogram()->SetYTitle("c");
      mps->GetHistogram()->SetYTitle("p");

      out->cd();
      c2s->Write();
      p2s->Write();
      can->Write();
      mcs->Write();
      mps->Write();
    } 
    else { 
      can->Divide(2,1);
      can->cd(1); cs->Draw("APL");
      can->cd(2); ps->Draw("APL");

      out->cd();
      can->Write();
      cs->Write();
      ps->Write();
    }
    PrintCanvas(can, 0, 0);
  }
};
#endif

//====================================================================
#ifdef TEST_FITTER
#include "AliLandauGausFitter.h"
#include <TRandom.h>

struct TestFit : public TestF
{
  struct Test : public TestF::Test 
  {
    TH1* dist;
    TF1* res;
    TH1* pars;
    Bool_t fDoShift;
    Test() : TestF::Test(), dist(0), res(0), pars(0), fDoShift(true) {}
    static void FillRandom(TH1* h, TF1* f, Int_t n=1000000) { 
      TAxis* xAxis = h->GetXaxis();

      Int_t first  = xAxis->GetFirst();
      Int_t last   = xAxis->GetLast();
      Int_t nbinsx = last-first+1;

      TArrayD integral(nbinsx+1);
      integral[0] = 0;
      for (Int_t i = 1; i <= nbinsx; i++) {
  Double_t fint = f->Integral(xAxis->GetBinLowEdge(i+first-1),
            xAxis->GetBinUpEdge(i+first-1));
  integral[i] = integral[i-1] + fint;
      }

      //   - Normalize integral to 1
      if (integral[nbinsx] == 0) {
  Error("FillRandom", "Integral = zero"); 
  return;
      }
      for (Int_t i = 1 ; i <= nbinsx; i++)  
  integral[i] /= integral[nbinsx];

      //   --------------Start main loop ntimes
      for (Int_t i = 0; i < n; i++) {
  Double_t r1 = gRandom->Rndm(i);
  Int_t    j  = TMath::BinarySearch(nbinsx,&integral[0],r1);
  Double_t x  = (xAxis->GetBinLowEdge(j+first) +
           xAxis->GetBinWidth(j+first) * 
           (r1-integral[j])/(integral[j+1] - integral[j]));
  h->Fill(x);
      }
    }
    void DoRun() 
    {
      AliLandauGaus::EnableSigmaShift(1); // Enable sigma-shift
      Double_t xMin = XMin();
      Double_t xMax = XMax();
      Int_t    n    = N();
      src->SetLineColor(kBlue);
      src->SetLineStyle(2);
      dist = new TH1D("dist", "Landau-Gaus distribution",500,xMin,xMax);
      dist->SetXTitle("#Delta");
      dist->SetYTitle("dN/d#Delta");
      dist->SetFillStyle(3001);
      dist->SetFillColor(kMagenta+1);
      dist->Sumw2();
      dist->SetDirectory(0);
      FillRandom(dist, src);
      
      Int_t    peakBin = dist->GetMaximumBin();
      Double_t max     = dist->GetBinContent(peakBin);
      dist->Scale(1/max);

      AliLandauGaus::EnableSigmaShift(fDoShift); // Disable sigma-shift
      AliLandauGausFitter f(xMin, xMax, 10);
      // f.SetDebug(true);
      
      TF1* r = f.FitNParticle(dist, n);
      if (!r) return;
      r->SetName("fit");

      // Get the status for the covariance calculation: 
      // 0: Not calculatuted 
      // 1: Approximate 
      // 2: Forced positive definite 
      // 3: Accurate 
      Int_t covStatus = 
  static_cast<TFitResult*>(f.GetFitResults().At(n-1))->CovMatrixStatus();
      if (covStatus < 2) ::Warning("", "Fit may not be good");
      
      r->SetRange(xMin, xMax);
      res = static_cast<TF1*>(r->Clone());
      res->SetLineColor(kPink+10);
      res->SetLineStyle(1);
      res->SetLineWidth(3);
      dist->GetListOfFunctions()->Add(res);
      
      Double_t min     = 0.1*dist->GetMinimum();
      TF1*     old     = src; // Tuck away src
      src              = res; // Use fit now
      for (Int_t i = 1; i <= n; i++) { 
  Double_t rDelta  = Delta();
  Double_t rxi     = Xi();
  Double_t rsigma  = Sigma();
  TF1*     comp    = Fi(i);
  comp->SetLineWidth(2);
  comp->SetLineStyle(1);

  // Modifies arguments!
  AliLandauGaus::IPars(i, rDelta, rxi, rsigma);
  TLine* l = VerticalLine(1, comp->GetLineColor(), rDelta, min,
        1.1*old->Eval(rDelta));
  dist->GetListOfFunctions()->Add(comp);
  dist->GetListOfFunctions()->Add(l);
      }
      src = old; // restore 
      
    }
    void DrawInPad(TVirtualPad* p)
    {
      // Double_t scale = src->GetMaximum();
      // src->SetParameter(0, 1/scale);
      p->Clear();
      p->SetTopMargin(0.0);
      p->SetRightMargin(0.0);
      p->SetBottomMargin(0.0);
      p->SetLeftMargin(0.0);
      p->SetFillColor(kGray);
      p->Divide(2,1,0.001,0.001);
      TVirtualPad* q = p->cd(1);
      q->SetLogy();
      q->SetTopMargin(0.01);
      q->SetRightMargin(0.02);
      q->SetFillColor(kWhite);
      gStyle->SetOptStat(0);
      gStyle->SetOptFit(11111);
      gStyle->SetOptTitle(false);
      gStyle->SetStatY(1-q->GetTopMargin());

      dist->DrawCopy("hist");
      dist->DrawCopy("e same");
      AliLandauGaus::EnableSigmaShift(1); // Enable sigma-shift
      src->DrawClone("same");
      AliLandauGaus::EnableSigmaShift(fDoShift); // Disable sigma-shift
      if (!res) return;


      q = p->cd(2);
      q->SetTopMargin(0.01);
      q->SetRightMargin(0.02);
      q->SetFillColor(kWhite);
      Int_t nPar = src->GetNpar();
      pars = new TH1D("pars", "#chi^{2}/#nu & #Deltap_{i}/#deltap_{i}", 
          nPar+1, 0, nPar+1);
      pars->SetDirectory(0);
      pars->SetFillColor(kCyan+2);
      pars->SetFillStyle(3001);
      pars->SetYTitle("#Deltap_{i}/#deltap_{i}");

      TPaveStats* stats = new TPaveStats(0,0,0,0);
      stats->SetX1NDC(gStyle->GetStatX()-gStyle->GetStatW());      
      stats->SetY1NDC(gStyle->GetStatY()-gStyle->GetStatH());
      stats->SetX2NDC(gStyle->GetStatX());
      stats->SetY2NDC(gStyle->GetStatY());
      stats->SetTextFont(42);
      
      stats->SetBorderSize(1);
      stats->SetFillColor(kWhite);
      stats->AddText("");
      pars->GetListOfFunctions()->Add(stats);
      TString fmt(Form("%%s = %%%s", stats->GetFitFormat()));
      for (Int_t i = 0; i < nPar; i++) { 
  Double_t low, high;
  res->GetParLimits(i, low, high);
  Double_t pSrc = src->GetParameter(i);
  Double_t pRes = res->GetParameter(i);
  Double_t eRes = res->GetParError(i);
  Bool_t   fix  = (low == high && ((pRes == 0 && low == 1) || low==pRes));
  Double_t diff = fix ? 0 : TMath::Abs(pSrc-pRes)/eRes;
  pars->GetXaxis()->SetBinLabel(i+1,res->GetParName(i));
  pars->SetBinContent(i+1, diff);
  stats->AddText(Form(fmt.Data(), res->GetParName(i), pSrc));
      }
      pars->GetXaxis()->SetBinLabel(nPar+1,"#chi^{2}/#nu");
      pars->SetBinContent(nPar+1,res->GetChisquare() / res->GetNDF());
      pars->DrawCopy();
      q->Modified();
      q->Update();

      p->Modified();
      p->Update();
      p->cd();
      
      PrintFs(false);
    }
    void PrintFs(Bool_t full=false) {
      if (full) {
  Int_t nPar = src->GetNpar();
  Printf("%-2s %-10s | %-8s | %8s %-9s | %9.5s", 
         "#", "Name", "Source", "Fit", "Error", "delta");
  for (Int_t i = 0; i < nPar; i++) { 
    Double_t low, high;
    res->GetParLimits(i, low, high);
    Double_t pSrc = src->GetParameter(i);
    Double_t pRes = res->GetParameter(i);
    Double_t eRes = res->GetParError(i);
    Bool_t   fix  = (low == high && ((pRes == 0 && low==1) || low==pRes));
    Double_t diff = fix ? 0 : TMath::Abs(pSrc-pRes)/eRes;
    Printf("%2d %10s | %8.4f | %8.4f %9.5f | %9.5f %s", 
     i, src->GetParName(i), pSrc, pRes, eRes, diff, 
     (fix ? "fixed" : (diff < 2 ? "ok" : "bad")));
  }
  Printf("chi^2/nu = %8.4f/%-3d = %9.5f", 
         res->GetChisquare(), res->GetNDF(), 
         res->GetChisquare() / res->GetNDF());
      }
      else 
  Printf("%8.4f", res->GetChisquare() / res->GetNDF());
    }
    virtual void WriteOut(TDirectory* d) 
    {
      d->cd();
      src->Write();
      dist->Write();
      if (res) res->Write();
      if (pars) pars->Write();
    }
    
  };
  //__________________________________________________________________
  TestF::Test* MakeTest() { return new Test(); } 
  const char* Prefix() const { return "fit"; }
  void PreLoop(UShort_t,Int_t) {
    Printf("%-8s | %-8s | %-8s | %-8s ",  "Delta_p", "xi", "sigma", "chi^2/nu");
  }
  void Step(UShort_t,Int_t,Int_t,Int_t,TestF::Test*,TNtuple*) {}
  void PostLoop(UShort_t,TCanvas* c,TNtuple*,TFile*) { PrintCanvas(c,0,0); }
};
#endif
// 
// EOF
//