TestELossDist.C

Go to the documentation of this file.

#include <iomanip>
#ifndef __CINT__
#include <iostream>
#include <TH1.h>
#include <TH2.h>
#include <TRandom.h>
#include <TCanvas.h>
#include <TArrayD.h>
#include <TStyle.h>
// #include <TFitResult.h>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TMath.h>
#include <THStack.h>
#include <TList.h>
#include <TLatex.h>
#include <TProfile.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#else
class TH1;
class TH2;
class TArrayD;
class TRandom;
class TCanvas;
class TF1;
#endif

//___________________________________________________________________
static Double_t landauGaus1(Double_t* xp, Double_t* pp);
static Double_t landauGausN(Double_t* xp, Double_t* pp);
static Double_t landauGausI(Double_t* xp, Double_t* pp);

//====================================================================
struct Function
{
  enum { 
    kC, 
    kDelta, 
    kXi, 
    kSigma, 
    kN, 
    kA 
  };
  static const Double_t fgkMPShift;
  static const Double_t fgkInvRoot2Pi;
  static const Double_t fgkConvNSigma;
  static const Double_t fgkConvNSteps;
  static Double_t  Landau(Double_t x, Double_t delta, Double_t xi)
  {
    return TMath::Landau(x, delta - xi * fgkMPShift, xi);
  }
  static Double_t  LandauGaus(Double_t x,  Double_t delta, 
            Double_t xi, Double_t sigma) 
  {
    Double_t deltap = delta - xi * fgkMPShift;
    Double_t xlow   = x - fgkConvNSigma * sigma;
    Double_t xhigh  = x + fgkConvNSigma * sigma;
    Double_t step   = (xhigh - xlow) / fgkConvNSteps;
    Double_t sum    = 0;

    for (Int_t i = 0; i < fgkConvNSteps / 2; i++) { 
      Double_t x1 = xlow  + (i - .5) * step;
      Double_t x2 = xhigh - (i - .5) * step;

      Double_t s1 = 
  TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma);
      Double_t s2 = 
  TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma);
      sum += s1 + s2;
    }
    Double_t ret = step * sum * fgkInvRoot2Pi / sigma;
    return ret;
  }
  static Double_t ILandauGaus(Int_t i, Double_t x, Double_t delta, 
            Double_t xi, Double_t sigma)
  {
    if (i <= 1) return LandauGaus(x, delta, xi, sigma);
    Double_t di      = i;
    Double_t deltai  = i * (delta + xi * TMath::Log(di));
    Double_t xii     = i * xi;
    Double_t sigmai  = TMath::Sqrt(di) * sigma;
    
    if (sigmai < 1e-10) return Landau(x, deltai, xii);
    Double_t ret =  LandauGaus(x, deltai, xii, sigmai);;
    // Info("ILandauGaus", "Fi(%f;%d,%f,%f,%f)->%f",
    //      x, i, deltai, xii, sigmai, ret);
    return ret;
  }
  static Double_t IdLandauGausdPar(Int_t    i,     Double_t x, 
           UShort_t ipar,  Double_t dPar, 
           Double_t delta, Double_t xi, 
           Double_t sigma)
  {
    if (dPar == 0) return 0;
    Double_t dp      = dPar;
    Double_t d2      = dPar / 2;
    Double_t deltaI  =  i * (delta + xi * TMath::Log(i));
    Double_t xiI     =  i * xi;
    Double_t si      =  TMath::Sqrt(Double_t(i));
    Double_t sigmaI  =  si*sigma;
    Double_t y1      = 0;
    Double_t y2      = 0;
    Double_t y3      = 0;
    Double_t y4      = 0;
    switch (ipar) {
    case 0: 
      y1 = ILandauGaus(i, x, deltaI+i*dp, xiI, sigmaI);
      y2 = ILandauGaus(i, x, deltaI+i*d2, xiI, sigmaI);
      y3 = ILandauGaus(i, x, deltaI-i*d2, xiI, sigmaI);
      y4 = ILandauGaus(i, x, deltaI-i*dp, xiI, sigmaI);
      break;
    case 1: 
      y1 = ILandauGaus(i, x, deltaI, xiI+i*dp, sigmaI);
      y2 = ILandauGaus(i, x, deltaI, xiI+i*d2, sigmaI);
      y3 = ILandauGaus(i, x, deltaI, xiI-i*d2, sigmaI);
      y4 = ILandauGaus(i, x, deltaI, xiI-i*dp, sigmaI);
      break;
    case 2: 
      y1 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*dp);
      y2 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*d2);
      y3 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*d2);
      y4 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*dp);
      break;
    default:
      return 0;
    } 
    
    Double_t d0  = y1 - y4;
    Double_t d1  = 2 * (y2 - y3);
    
    Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;
    
    return g;
  }
  static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
            Double_t sigma, Int_t n, Double_t* a)
  {
    Double_t res = LandauGaus(x, delta, xi, sigma);
    for (Int_t i = 2; i <= n; i++) 
      res += a[i-2] * ILandauGaus(i, x, delta, xi, sigma);
    return res;
  }
  Double_t NEstimate(Double_t x, Double_t delta, Double_t xi, 
         Double_t sigma, Int_t n, Double_t* a)
  {
    Double_t num = LandauGaus(x, delta, xi, sigma);
    Double_t den = num;
    for (Int_t i = 2; i <= n; i++) {
      Double_t f =  ILandauGaus(i, x, delta, xi, sigma);
      num        += a[i-2] * f * i;
      den        += a[i-2] * f;
    }
    if (den < 1e-4) return 0;
    return num / den;
  }
  Double_t DEstimatePar(Double_t x, Int_t ipar, Double_t delta, Double_t xi, 
      Double_t sigma, UShort_t n, Double_t* a, 
      Double_t dPar=0.001) 
  {
    Double_t dp    = dPar;
    Double_t d2    = dPar / 2;
    Double_t y1    = 0;
    Double_t y2    = 0;
    Double_t y3    = 0;
    Double_t y4    = 0;
    switch (ipar) { 
    case Function::kC: 
    case Function::kN: 
      return 0;
    case Function::kDelta:
      y1 = NEstimate(x, delta+dp, xi, sigma, n, a);
      y2 = NEstimate(x, delta+d2, xi, sigma, n, a);
      y3 = NEstimate(x, delta-d2, xi, sigma, n, a);
      y4 = NEstimate(x, delta-dp, xi, sigma, n, a);
      break;
    case Function::kXi:
      y1 = NEstimate(x, delta, xi+dp, sigma, n, a);
      y2 = NEstimate(x, delta, xi+d2, sigma, n, a);
      y3 = NEstimate(x, delta, xi-d2, sigma, n, a);
      y4 = NEstimate(x, delta, xi-dp, sigma, n, a);
      break;
    case Function::kSigma:
      y1 = NEstimate(x, delta, xi, sigma+dp, n, a);
      y2 = NEstimate(x, delta, xi, sigma+d2, n, a);
      y3 = NEstimate(x, delta, xi, sigma-d2, n, a);
      y4 = NEstimate(x, delta, xi, sigma-dp, n, a);
      break;
    default: {
      Int_t j = ipar-kA;
      if (j+1 > n) return 0;
      Double_t aa = a[j];
      a[j] = aa+dp; y1 = NEstimate(x, delta, xi, sigma, n, a);
      a[j] = aa+d2; y2 = NEstimate(x, delta, xi, sigma, n, a);
      a[j] = aa-d2; y3 = NEstimate(x, delta, xi, sigma, n, a);
      a[j] = aa-dp; y4 = NEstimate(x, delta, xi, sigma, n, a);
      a[j] = aa;
    }
      break;
    }
    Double_t d0  = y1 - y4;
    Double_t d1  = 2 * (y2 - y3);
    
    Double_t g   = 1/(2*dp) * (4*d1 - d0) / 3;    
    return g;
  }
  static TF1* MakeFunc(Int_t n, Double_t c, Double_t delta, 
           Double_t xi, Double_t sigma, Double_t* a, 
           Double_t xmin, Double_t xmax)
  {
    Int_t nPar = kN + n;
    TF1* f = new TF1(Form("landGaus%d",n), &landauGausN, xmin, xmax, nPar);
    f->SetNpx(500);
    f->SetParName(kC,      "C");
    f->SetParName(kDelta,  "#Delta_{p}");
    f->SetParName(kXi,     "#xi");
    f->SetParName(kSigma,  "#sigma");
    f->SetParName(kN,      "N");
    f->SetParameter(kC,     c);
    f->SetParameter(kDelta, delta);
    f->SetParameter(kXi,    xi);
    f->SetParameter(kSigma, sigma);
    f->FixParameter(kN,     n);
    f->SetParLimits(kDelta, xmin, 4);
    f->SetParLimits(kXi,    0.00, 4);
    f->SetParLimits(kSigma, 0.01, 0.11);

    for (Int_t i = 2; i <= n; i++) { 
      Int_t j = i-2;
      f->SetParName(kA+j, Form("a_{%d}", i));
      f->SetParameter(kA+j, a[j]);
      f->SetParLimits(kA+j, 0, 1);
    }
    return f;
  }
  static TF1* MakeIFunc(Int_t i, Double_t c, Double_t delta, 
      Double_t xi, Double_t sigma,
      Double_t xmin, Double_t xmax)
  {
    Int_t nPar = 5;
    TF1* f = new TF1(Form("landGausI%d",i), &landauGausI, xmin, xmax, nPar);
    f->SetNpx(100);
    f->SetParName(kC,      "C");
    f->SetParName(kDelta,  "#Delta_{p}");
    f->SetParName(kXi,     "#xi");
    f->SetParName(kSigma,  "#sigma");
    f->SetParName(kN,      "N");
    f->SetParameter(kC,     c);
    f->SetParameter(kDelta, delta);
    f->SetParameter(kXi,    xi);
    f->SetParameter(kSigma, sigma);
    f->FixParameter(kN,     i);
    return f;
  }

  // --- Object code -------------------------------------------------
  Function(Int_t n, Double_t c, Double_t delta, Double_t xi, Double_t sigma, 
     Double_t* a, Double_t xmin, Double_t xmax)
    : fF(MakeFunc(n, c, delta, xi, sigma, a, xmin, xmax))
  {}
  Function(TF1* f) : fF(f) {}
  TF1* GetF1() { return fF; }
  Double_t Evaluate(Double_t x) const
  {
    return fF->Eval(x);
  }
  Double_t Evaluate(Double_t x, Double_t& e) const
  {
    Double_t delta     = GetDelta();
    Double_t xi        = GetXi();
    Double_t sigma     = GetSigma();
    Double_t dFdDelta2 = 0;
    Double_t dFdXi2    = 0;
    Double_t dFdSigma2 = 0;
    e              = 0;
    for (Int_t i = 1; i <= GetN(); i++) { 
      Double_t a      = GetA(i);
      Double_t dDelta = a*IdLandauGausdPar(i, x, 1, 0.001, delta, xi, sigma); 
      Double_t dXi    = a*IdLandauGausdPar(i, x, 2, 0.001, delta, xi, sigma);
      Double_t dSigma = a*IdLandauGausdPar(i, x, 3, 0.001, delta, xi, sigma);
      Double_t dFda   = a*ILandauGaus(i, x, delta, xi, sigma);
      dFdDelta2 += dDelta * dDelta;
      dFdXi2    += dXi    * dXi;
      dFdSigma2 += dSigma * dSigma;
      e += TMath::Power(dFda * GetEA(i),2);
    }      
    Double_t edelta = GetEDelta();
    Double_t exi    = GetEXi();
    Double_t esigma = GetESigma();
    e += (dFdDelta2 * edelta * edelta + 
    dFdXi2    * exi    * exi    + 
    dFdSigma2 * esigma * esigma);
    e = TMath::Sqrt(e);
    return fF->Eval(x);
  }

  Double_t EstimateNParticles(Double_t x, Short_t maxN=-1)
  {
    UShort_t n     = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
    Double_t delta = GetDelta();
    Double_t xi    = GetXi();
    Double_t sigma = GetSigma();
    return NEstimate(x, delta, xi, sigma, n, GetAs());
  }
  Double_t EstimateNParticles(Double_t x, Double_t& e, Short_t maxN=-1)
  {
    
    UShort_t  n      = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
    Double_t  delta  = GetDelta();
    Double_t  xi     = GetXi();
    Double_t  sigma  = GetSigma();
    Double_t* a      = GetAs();
    Double_t  dDelta = (DEstimatePar(x,Function::kDelta,delta,xi,sigma,n,a,\
             0.01 * GetEDelta()) * GetEDelta());
    Double_t  dXi    = (DEstimatePar(x,Function::kXi,   delta,xi,sigma,n,a,
             0.01 * GetEXi())    * GetEXi());
    Double_t  dSigma = (DEstimatePar(x,Function::kSigma,delta,xi,sigma,n,a,
             0.01 * GetESigma()) * GetESigma());
    e                = dDelta * dDelta + dXi * dXi + dSigma * dSigma;
    for (Int_t i = 2; i <= n; i++) { 
      Double_t dAi = (DEstimatePar(x, Function::kA+i-2,delta,xi,sigma,n,a,
           0.01 * GetEA(i)) * GetEA(i));
      e += dAi * dAi;
    }
    e = TMath::Sqrt(e);
    return NEstimate(x, GetDelta(), GetXi(), GetSigma(), n, GetAs());
  }
  Double_t RandomEstimateNParticles(Double_t x, Short_t maxN=-1)
  {
    UShort_t n     = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
    TArrayD  p(n);
    Double_t delta = GetDelta();
    Double_t xi    = GetXi();
    Double_t sigma = GetSigma();
    
    if (Evaluate(x) < 1e-4) return 0;

    // std::cout << "RandomEstimateNParticles(" << x << "):";
    for (Int_t i = 1; i <= n; i++) { 
      Double_t a = GetA(i);
      Double_t f = ILandauGaus(i, x, delta, xi, sigma);
      p[i-1]     = a * f;
      // std::cout << p[i-1] << ",";
    }
    Double_t r = gRandom->Uniform(p.GetSum());
    Double_t l = 0;
    // std::cout << "s=" << p.GetSum() << ",r=" << r << ",";
    for (Int_t i = 1; i <= n; i++) { 
      if (r > l && r <= l + p[i-1]) {
  // std::cout << "l=" << l << ",l+p=" << l+p[i-1] << "-> " << i 
  //  << std::endl;
  return i;
      }
      l += p[i-1];
    }
    return 0;
  }
  void Draw(Option_t* option="") { fF->Draw(option); }
  Double_t GetC() const { return fF->GetParameter(kC); }
  Double_t GetDelta() const { return fF->GetParameter(kDelta); }
  Double_t GetXi() const { return fF->GetParameter(kXi); }
  Double_t GetSigma() const { return fF->GetParameter(kSigma); }
  UShort_t GetN() const { return fF->GetParameter(kN); }
  Double_t GetA(Int_t i) const { return i == 1 ? 1 : fF->GetParameter(kA+i-2);}
  Double_t* GetAs() const { return &(fF->GetParameters()[kA]);}
  Double_t GetEC() const { return fF->GetParError(kC); }
  Double_t GetEDelta() const { return fF->GetParError(kDelta); }
  Double_t GetEXi() const { return fF->GetParError(kXi); }
  Double_t GetESigma() const { return fF->GetParError(kSigma); }
  UShort_t GetEN() const { return 0; }
  Double_t GetEA(Int_t i) const { return i == 1 ? 0 : fF->GetParError(kA+i-2);}
  Double_t* GetEAs() const { return &(fF->GetParErrors()[kA]);}
  TF1* fF;
};
#ifndef __CINT__
const Double_t Function::fgkMPShift    = -0.22278298;
const Double_t Function::fgkInvRoot2Pi = 1. / TMath::Sqrt(2*TMath::Pi());
const Double_t Function::fgkConvNSigma = 5;
const Double_t Function::fgkConvNSteps = 100;
#endif
  
//====================================================================
struct Fitter
{
  // --- Object code ------------------------------------------------
  Fitter(Int_t nMinus) 
    : fNMinus(nMinus)
  {
  }
  TF1* FitFirst(TH1* dist)
  {
    fFunctions.Clear();
    
    // Get upper edge 
    Int_t    midBin = dist->GetMaximumBin();
    Double_t midE   = dist->GetBinLowEdge(midBin);
    
    // Get low edge 
    Int_t    minBin = midBin - fNMinus;
    Double_t minE  = dist->GetBinCenter(minBin);
    Double_t maxE  = dist->GetBinCenter(midBin+2*fNMinus);
    Double_t a[]   = { 0 };
    TF1* f = Function::MakeFunc(1, 0.5, midE, 0.07, 0.1, a, minE, maxE);

    // Do the fit 
    dist->Fit(f, "RNQS", "", minE, maxE);
    fFunctions.AddAtAndExpand(f, 0);
    return f;
  }
  Function* FitN(TH1* dist, Int_t n)
  {
    if (n == 1) return new Function(FitFirst(dist));
    TF1*        f1 = static_cast<TF1*>(fFunctions.At(0));
    if (!f1) { 
      f1 = FitFirst(dist);
      if (!f1) { 
  ::Warning("FitN", "First fit missing");
  return 0;
      }
    }
    
    Double_t delta1 = f1->GetParameter(Function::kDelta);
    Double_t xi1    = f1->GetParameter(Function::kXi);
    Double_t maxEi  = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
    Double_t minE   = f1->GetXmin();
    
    TArrayD a(n-1);
    for (UShort_t i = 2; i <= n; i++) a.fArray[i-2] = (n==2? 0.05 : 0.000001);
    
    TF1* fn = Function::MakeFunc(n, f1->GetParameter(Function::kC), 
         delta1, xi1, 
         f1->GetParameter(Function::kSigma), 
         a.fArray, minE, maxEi);
    
    dist->Fit(fn, "RSNQ", "", minE, maxEi);
    dist->Fit(fn, "RSNQM", "", minE, maxEi);
    fFunctions.AddAtAndExpand(fn, n-1);
    
    return new Function(fn);
  }
  Int_t     fNMinus;
  TObjArray fFunctions;
};

//====================================================================
static Double_t landauGaus1(Double_t* xp, Double_t* pp)
{
  Double_t x     = xp[0];
  Double_t c     = pp[Function::kC];
  Double_t delta = pp[Function::kDelta];
  Double_t xi    = pp[Function::kXi];
  Double_t sigma = pp[Function::kSigma];
    
  return c * Function::LandauGaus(x, delta, xi, sigma);
}
static Double_t landauGausN(Double_t* xp, Double_t* pp)
{
  Double_t  x     = xp[0];
  Double_t  c     = pp[Function::kC];
  Double_t  delta = pp[Function::kDelta];
  Double_t  xi    = pp[Function::kXi];
  Double_t  sigma = pp[Function::kSigma];
  Int_t     n     = Int_t(pp[Function::kN]);
  Double_t* a     = &pp[Function::kA];
    
  return c * Function::NLandauGaus(x, delta, xi, sigma, n, a);
}
static Double_t landauGausI(Double_t* xp, Double_t* pp)
{
  Double_t  x     = xp[0];
  Double_t  c     = pp[Function::kC];
  Double_t  delta = pp[Function::kDelta];
  Double_t  xi    = pp[Function::kXi];
  Double_t  sigma = pp[Function::kSigma];
  Int_t     i     = Int_t(pp[Function::kN]);
    
  return c * Function::ILandauGaus(i, x, delta, xi, sigma);
}

//====================================================================
struct Generator
{
  Double_t fDelta;
  Double_t fXi;
  Double_t fSigma;
  Int_t    fMaxN;
  TH1*     fSingle;
  TH1*     fSummed;
  TH1*     fSimple;
  TH2*     fNSum;
  TH1*     fN;
  TArrayD  fRP; // Raw probabilities
  TArrayD  fP; // Normalized probabilities 
  TF1*     fF;
  Int_t    fNObs;

  Generator(Double_t delta=.55, Double_t xi=0.06, Double_t sigma=0.02,
      Int_t maxN=5, const Double_t* a=0)
    : fDelta(delta), fXi(xi), fSigma(sigma), fMaxN(maxN),
      fSingle(0), fSummed(0), fSimple(0), fNSum(0), fRP(), fP(fMaxN), fF(0),
      fNObs(0)
  {
    fMaxN = TMath::Min(fMaxN, 5);
    const Double_t ps[] = { 1, .05, .02, .005, .002 };
    fRP.Set(5, ps);
    if (a) { 
      for (Int_t i = 0; i < maxN; i++) fRP[i] = a[i];
    }

    Double_t sum = 0; 
    for (Int_t i = 0; i < fMaxN; i++) sum += fRP[i];
    for (Int_t i = 0; i < fP.fN; i++) {
      fP[i] = fRP[i] / sum;
      if (i != 0) fP[i] += fP[i-1];
    }

    fSingle = new TH1D("single", "Single signals", 500, 0, 10);
    fSingle->SetXTitle("Signal");
    fSingle->SetYTitle("Events");
    fSingle->SetFillColor(kRed+1);
    fSingle->SetMarkerColor(kRed+1);
    fSingle->SetLineColor(kRed+1);
    fSingle->SetFillStyle(3001);
    fSingle->SetMarkerStyle(21);
    fSingle->SetDirectory(0);
    fSingle->Sumw2();

    fSummed = static_cast<TH1D*>(fSingle->Clone("summed"));
    fSummed->SetTitle("Summed signals");
    fSummed->SetFillColor(kBlue+1);
    fSummed->SetMarkerColor(kBlue+1);
    fSummed->SetLineColor(kBlue+1);
    fSummed->SetDirectory(0);

    fSimple = static_cast<TH1D*>(fSingle->Clone("summed"));
    fSimple->SetTitle("Summed signals");
    fSimple->SetFillColor(kGreen+1);
    fSimple->SetMarkerColor(kGreen+1);
    fSimple->SetLineColor(kGreen+1);
    fSimple->SetDirectory(0);

    fNSum = new TH2D("single", "Observation signals", 
         fMaxN, .5, fMaxN+.5, 100, 0, 10);
    fNSum->SetMarkerColor(kMagenta+1);
    fNSum->SetLineColor(kMagenta+1);
    fNSum->SetFillColor(kMagenta+1);
    fNSum->SetYTitle("Signal");
    fNSum->SetXTitle("N");
    fNSum->SetZTitle("Events");
    fNSum->SetDirectory(0);
    fNSum->Sumw2();

    fN = new TH1D("n", "Number of particles", fMaxN, .5, fMaxN+.5);
    fN->SetXTitle("N");
    fN->SetYTitle("Events");
    fN->SetMarkerColor(kMagenta+1);
    fN->SetLineColor(kMagenta+1);
    fN->SetFillColor(kMagenta+1);
    fN->SetMarkerStyle(22);
    fN->SetFillStyle(3001);
    fN->SetDirectory(0);

    std::cout << "Probabilities:" << std::setprecision(4) 
        << "\n  ";
    for (Int_t i = 0; i < fRP.fN; i++) 
      std::cout << std::setw(7) << fRP[i] << " ";
    std::cout << " = " << fRP.GetSum() << "\n  ";
    for (Int_t i = 0; i < fP.fN; i++) 
      std::cout << std::setw(7) << fP[i] << " ";
    std::cout << std::endl;

    Double_t aa[] = { 0 };
    fF = Function::MakeFunc(1, 1, fDelta, fXi, fSigma, aa,  0, 10);
  }
  void Clear()
  {
    fSingle->Reset();
    fSummed->Reset();
    fNSum->Reset();
    fN->Reset();
  }
  Int_t GenerateN()
  {
    Double_t rndm = gRandom->Uniform();
    Int_t    ret  = 0;
    if      (                 rndm < fP[0]) ret = 1;
    else if (rndm >= fP[0] && rndm < fP[1]) ret = 2;
    else if (rndm >= fP[1] && rndm < fP[2]) ret = 3;
    else if (rndm >= fP[2] && rndm < fP[3]) ret = 4;
    else if (rndm >= fP[3] && rndm < fP[4]) ret = 5;
    // Printf("%f -> %d", rndm, ret);
    fN->Fill(ret);
    return ret;
  }
  Double_t Generate1Signal()
  {
    Double_t ret = 0;
    if (fF) ret = fF->GetRandom();
    else { 
      Double_t rmpv = gRandom->Gaus(fDelta, fSigma);
      ret           = gRandom->Landau(rmpv, fXi);
    }
    fSingle->Fill(ret);
    fSimple->Fill(gRandom->Landau(fDelta - fXi * Function::fgkMPShift, fXi));
    return ret;
  }
  Double_t GenerateNSignal(Int_t n)
  {
    Double_t ret = 0;
    for (Int_t i = 0; i < n; i++) ret += Generate1Signal();
    fSummed->Fill(ret);
    fNSum->Fill(n, ret);
    return ret;
  }
  void Generate(Int_t nObs=1000, Bool_t reset=true)
  {
    if (reset) Clear();
    fNObs = nObs;

    for (Int_t i = 0; i < nObs; i++) {
      // if (((i+1) % (nObs/10)) == 0) 
      //  std::cout << "Event " << std::setw(6) << i << std::endl;
      Int_t n = GenerateN();
      GenerateNSignal(n);
    }
    Double_t m = fSingle->GetMaximum();
    fSingle->Scale(1. / m);
    m = fSummed->GetMaximum();
    fSummed->Scale(1. / m);
    m = fSimple->GetMaximum();
    fSimple->Scale(1. / m);

    std::cout << "Resulting probabilities:\n  ";
    for (Int_t i = 1; i <= fN->GetNbinsX(); i++) 
      std::cout << std::setprecision(4) << std::setw(7) 
    << fN->GetBinContent(i) << " ";
    std::cout << std::endl;
  }
  void PrintParameter(Double_t input, TF1* f, Int_t i)
  {
    Double_t val = 0, err = 0;
    if (i < f->GetNpar()) {
      val = f->GetParameter(i);
      err = f->GetParError(i);
    }
    std::cout << "  " << std::setw(16) << f->GetParName(i) << ": "
        << std::fixed 
        << std::setprecision(4) << std::setw(6) << input << " -> " 
        << std::setprecision(4) << std::setw(6) << val   << " +/- "
        << std::setprecision(5) << std::setw(7) << err   << " [delta:" 
        << std::setw(3) << int(100 * (input-val) / input) << "%,err:" 
        << std::setw(3) << int(100 * err / val) << "%]"
        << std::endl;
  }
  void PrintResult(TF1* f) 
  {
    Double_t chi2 = f->GetChisquare();
    Int_t    ndf  = f->GetNDF();
    std::cout << "Chi^2/NDF = " << chi2 << "/" << ndf << "=" 
        << chi2/ndf << std::endl;
    PrintParameter(1,      f, Function::kC);
    PrintParameter(fDelta, f, Function::kDelta);
    PrintParameter(fXi,    f, Function::kXi);
    PrintParameter(fSigma, f, Function::kSigma);
    PrintParameter(fMaxN,  f, Function::kN);
    for (Int_t i = 0; i < fMaxN-1; i++) 
      PrintParameter(fRP[i+1], f, Function::kA+i);
  }
  void DrawParameter(Double_t x, Double_t y, Double_t input, TF1* f, Int_t i)
  {
    Double_t val = 0, err = 0;
    if (i < f->GetNpar()) {
      val = f->GetParameter(i);
      err = f->GetParError(i);
    }

    TLatex* ltx = new TLatex(x, y, f->GetParName(i));
    ltx->SetTextSize(.04);
    ltx->SetTextFont(132);
    ltx->SetTextAlign(11);
    ltx->SetNDC();
    ltx->Draw();
    
    ltx->DrawLatex(x+.05,  y, Form("%5.3f", input));
    ltx->DrawLatex(x+.14, y, "#rightarrow");
    ltx->SetTextAlign(21);
    ltx->DrawLatex(x+.3,  y, Form("%5.3f #pm %6.4f", val, err));
    ltx->SetTextAlign(11);
    ltx->DrawLatex(x+.48,  y, Form("[#Delta=%3d%%, #delta=%3d%%]",
          int(100 * (input-val)/input), 
          int(100*err/val)));
  }
  void DrawResult(Double_t x, Double_t y, TF1* f)
  {
    Double_t chi2 = f->GetChisquare();
    Int_t    ndf  = f->GetNDF();
    TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/NDF=%5.2f/%3d=%6.3f", 
          chi2, ndf, chi2/ndf));
    ltx->SetTextSize(.04);
    ltx->SetNDC();
    ltx->SetTextFont(132);
    ltx->Draw();

    x += .05;
    y -= .035; DrawParameter(x, y, 1,      f, Function::kC);
    y -= .035; DrawParameter(x, y, fDelta, f, Function::kDelta);
    y -= .035; DrawParameter(x, y, fXi,    f, Function::kXi);
    y -= .035; DrawParameter(x, y, fSigma, f, Function::kSigma);
    y -= .035; DrawParameter(x, y, fMaxN,  f, Function::kN);
    for (Int_t i = 0; i < fMaxN-1; i++) {
      y -= .035; DrawParameter(x, y, fRP[i+1], f, Function::kA+i);
    }
  }
  void Draw(const char* type="png")
  {
    gStyle->SetOptTitle(0);
    gStyle->SetOptStat(0);
    gStyle->SetPalette(1);

    TCanvas* c = new TCanvas("c", "Generator distributions", 1200, 1000);
    c->SetFillColor(0);
    c->SetBorderSize(0);
    c->SetBorderMode(0);
    c->SetTopMargin(0.01);
    c->SetRightMargin(0.01);

    c->Divide(2, 2);

    // --- First pad:  Data (multi, single, landau) + Fit -------------
    TVirtualPad* p = c->cd(1);
    p->SetLogy();
    p->SetBorderSize(0);
    p->SetBorderMode(0);
    p->SetTopMargin(0.01);
    p->SetRightMargin(0.01);
    fSummed->Draw();
    fSingle->Draw("same");
    fSimple->Draw("same");

    if (fF) { 
      Double_t max = fF->GetMaximum();
      Double_t a[] = { 0 };
      TF1* fcopy = Function::MakeFunc(1, 
              fF->GetParameter(Function::kC) /max, 
              fF->GetParameter(Function::kDelta), 
              fF->GetParameter(Function::kXi), 
              fF->GetParameter(Function::kSigma), 
              a, 0, 10);
      fcopy->SetLineColor(2);
      fcopy->SetLineStyle(2);
      fcopy->Draw("same");
    }

    Fitter*   fitter = new Fitter(4);
    Function* fit    = fitter->FitN(fSummed, 5);
    TF1*      f      = fit->GetF1();
    f->Draw("same");
    f->SetLineColor(kBlack);
    TF1* fst = static_cast<TF1*>(fitter->fFunctions.At(0));
    fst->SetLineWidth(3);
    fst->SetLineStyle(2);
    fst->SetLineColor(kMagenta);
    fst->Draw("same");
    PrintResult(f);
    DrawResult(.25, .95, f);
    TGraphErrors* gr = new TGraphErrors(fSummed->GetNbinsX());
    gr->SetName("fitError");
    gr->SetTitle("Fit w/errors");
    gr->SetFillColor(kYellow+1);
    gr->SetFillStyle(3001);
    for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) { 
      Double_t x = fSummed->GetBinCenter(j);
      Double_t e = 0;
      Double_t y = fit->Evaluate(x, e);
      gr->SetPoint(j-1, x, y);
      gr->SetPointError(j-1, 0, e);
    }
    gr->Draw("same l e4");

    TLegend* l = new TLegend(.15, .11, .5, .4);
    l->SetBorderSize(0);
    l->SetFillColor(0);
    l->SetFillStyle(0);
    l->SetTextFont(132);
    l->AddEntry(fSimple, "L: Single-L", "lp");
    l->AddEntry(fSingle, "LG: Single-L#otimesG", "lp");
    l->AddEntry(fSummed,"NLG: Multi-L#otimesG", "lp");
    l->AddEntry(gr, "f_{N}: Fit to NLG");
    l->Draw();

    // --- Second pad:  Data (1, 2, 3, ...) --------------------------
    p = c->cd(2);
    p->SetLogy();
    p->SetBorderSize(0);
    p->SetBorderMode(0);
    p->SetTopMargin(0.01);
    p->SetRightMargin(0.01);
    p->SetGridx();
    THStack* stack = new THStack(fNSum, "y");
    TIter next(stack->GetHists());
    TH1* h = 0;
    Int_t i = 2;
    Double_t max = 0;
    TLegend* l2 = new TLegend(.6, .5, .95, .95);
    l2->SetBorderSize(0);
    l2->SetFillStyle(0);
    l2->SetFillColor(0);
    l2->SetTextFont(132);
    while ((h = static_cast<TH1*>(next()))) { 
      if (i == 2) max = h->GetMaximum();
      h->SetLineColor(i); 
      h->SetFillColor(i); 
      h->SetFillStyle(3001);
      h->Scale(1/max);
      l2->AddEntry(h, Form("%d particles", i - 1), "f");
      i++;
    }
    stack->Draw("hist nostack");
    stack->GetHistogram()->SetXTitle("Signal");
    stack->GetHistogram()->SetYTitle("Events");
    i = 2;
    for (Int_t j = 1; j <= 5; j++) {
      TF1* fi = Function::MakeIFunc(j, fRP[j-1], fDelta, fXi, fSigma, 0, 10);
      if (j == 1) max = fi->GetMaximum();
      fi->SetParameter(0, fRP[j-1]/max);
      fi->SetLineColor(i);
      fi->SetLineWidth(3);
      fi->SetLineStyle(2);
      fi->Draw("same");

      TF1* fj = Function::MakeIFunc(j, fit->GetC() * fit->GetA(j), 
            fit->GetDelta(), fit->GetXi(), 
            fit->GetSigma(), 0, 10);
      fj->SetLineColor(i);
      fj->SetLineWidth(2);
      fj->Draw("same");

      std::cout << "Component " << j << " scale=" << fi->GetParameter(0) 
    << " (" << fRP[j-1] << ")" << " " << fj->GetParameter(0) 
    << " (" << fit->GetC() << "*" << fit->GetA(j) << ")" 
    << std::endl;
      i++;
    }
    TLegendEntry* ee = l2->AddEntry("", "Input f_{i}", "l");
    ee->SetLineStyle(2);
    ee->SetLineWidth(3);
    l2->AddEntry("", "Fit f_{i}", "l");
    l2->Draw();
    // fNSum->Draw("lego2");

    // --- Third pad:  Mean multiplicity ------------------------------
    TH1* nhist1 = new TH1F("nhist1", "NHist", 5, 0.5, 5.5);
    nhist1->SetFillColor(kCyan+1);
    nhist1->SetMarkerColor(kCyan+1);
    nhist1->SetLineColor(kCyan+1);
    nhist1->SetFillStyle(3002);
    TH1* nhist2 = new TH1F("nhist2", "NHist", 5, 0.5, 5.5);
    nhist2->SetFillColor(kYellow+1);
    nhist2->SetMarkerColor(kYellow+1);
    nhist2->SetLineColor(kYellow+1);
    nhist2->SetFillStyle(3002);

    TH2* resp = new TH2F("resp", "Reponse", 100, 0, 10, 6, -.5, 5.5);
    resp->SetXTitle("Signal");
    resp->SetYTitle("# particles");
    for (Int_t j = 0; j < fNObs; j++) { 
      if (((j+1) % (fNObs / 10)) == 0) 
  std::cout << "Event # " << j+1 << std::endl;
      Double_t x = fSummed->GetRandom();
      Double_t y = fit->RandomEstimateNParticles(x);
      nhist1->Fill(y);
      resp->Fill(x, y);
    }
    TGraphErrors* graph = new TGraphErrors(fSummed->GetNbinsX());
    graph->SetName("evalWeighted");
    graph->SetTitle("Evaluated function");
    graph->SetLineColor(kBlack);
    graph->SetFillColor(kYellow+1);
    graph->SetFillStyle(3001);
    for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) { 
      Double_t x = fSummed->GetBinCenter(j);
      Double_t e = 0;
      Double_t y = fit->EstimateNParticles(x, e);
      nhist2->Fill(y, fSummed->GetBinContent(j));
      graph->SetPoint(j-1, x, y);
      graph->SetPointError(j-1, 0, e);
    }
    nhist1->Scale(1. / nhist1->GetMaximum());
    nhist2->Scale(1. / nhist2->GetMaximum());
    fN->Scale(1. / fN->GetMaximum());

    p = c->cd(3);
    p->SetLogy();
    p->SetBorderSize(0);
    p->SetBorderMode(0);
    p->SetTopMargin(0.01);
    p->SetRightMargin(0.01);
    fN->Draw();
    nhist1->Draw("same");
    nhist2->Draw("same");

    TLegend* l3 = new TLegend(.3, .7, .95, .95);
    l3->SetBorderSize(0);
    l3->SetFillStyle(0);
    l3->SetFillColor(0);
    l3->SetTextFont(132);
    l3->AddEntry(fN,  
     Form("Input n distribution: #bar{m}=%5.3f, s_{m}=%5.3f",
          fN->GetMean(), fN->GetRMS()), "f");
    l3->AddEntry(nhist1, 
     Form("Random draw from fit: #bar{m}=%5.3f, s_{m}=%5.3f", 
          nhist1->GetMean(), nhist1->GetRMS()), "f");
    l3->AddEntry(nhist2, 
     Form("Weighted evaluation of fit: #bar{m}=%5.3f, s_{m}=%5.3f", 
          nhist2->GetMean(), nhist2->GetRMS()), "f");
    l3->Draw();
    
    // --- Fourth pad: Reponse function ------------------------------
    p = c->cd(4);
    // p->SetLogy();
    p->SetBorderSize(0);
    p->SetBorderMode(0);
    p->SetTopMargin(0.01);
    p->SetRightMargin(0.01);
    p->SetGridx();
    TProfile* prof1 = fNSum->ProfileY();
    prof1->SetLineColor(kMagenta+1);
    prof1->Draw();
    TProfile* prof2 = resp->ProfileX();
    prof2->SetLineColor(kCyan+2);
    prof2->Draw("same");
    graph->SetLineWidth(2);
    graph->Draw("same LP E4");

    TLegend* l4 = new TLegend(.2, .11, .8, .4);
    l4->SetBorderSize(0);
    l4->SetFillStyle(0);
    l4->SetFillColor(0);
    l4->SetTextFont(132);
    l4->AddEntry(prof1, "Input distribution of N particles", "l");
    l4->AddEntry(prof2, "Random estimate of N particles", "l");
    l4->AddEntry(graph, "E: #sum_{i}^{N}i a_{i}f_{i}/"
     "#sum_{i}^{N}a_{i}f_{i} evaluated over NLG", "lf");
    l4->Draw();

    c->cd();
    
    if (type && type[0] != '\0') 
      c->Print(Form("TestELossDist.%s", type));
  }
};

void
TestELossDist(const char* type="png")
{
  gRandom->SetSeed(12345);
  const Double_t a1[] = { 1, .05, .02, .005, .002 }; // 7TeV pp
  const Double_t a2[] = { 1, .1, .05, .01, .005 };
  const Double_t a3[] = { 1, .2, .10, .05, .01 };

  Generator* g = new Generator(0.55, 0.06, 0.06, 5, a1);
  g->Generate(100000);
  g->Draw(type);

  
}