#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TRandom.h"

Double_t p(Double_t x);
Double_t pp(Double_t x);

const    Int_t npx = 100;
const    Int_t npz = 4000;

const    Double_t dz  = 1. / Double_t(npz);
const    Double_t dx  = 1. / Double_t(npx);

const    Double_t eps = dz/2.;

void dglap(Double_t qmin = .2, Double_t qmax = 10.)
{

    Double_t x   = 0.;
    Double_t z   = 0.;

    TF1* g = new TF1("*g", 0., 1., 0);
    g->SetNpx(npx);
    g->SetSavedPoint(npx + 1,  0.);
    g->SetSavedPoint(npx + 2,  1.);

    TF1* kern = new TF1("*k", 0., 1., 0);
    kern->SetNpx(npz);
    kern->SetSavedPoint(npz + 1,  0.);
    kern->SetSavedPoint(npz + 2,  1.);

    for (Int_t i = 0; i < npx + 1; i++) {
      x = i * dx;
      g->SetSavedPoint(i,  TMath::Exp(- (x-1.) * (x-1) / 2.e-5 ));
    }

    Double_t tmax = 2. * TMath::Log(qmax/0.2);
    Double_t tmin = 2. * TMath::Log(qmin/0.2);
    Double_t t = tmax;

    Double_t dt = 0.05;
    Double_t kint = 0.;

    for (Int_t j = 0; j < npz; j++) {
      z = j * dz + eps;
      Double_t w1 = p(z);
      kint += w1;
    } // z loop

    Double_t integral1 = kint * dz;
    TF1* f = 0;
    while (1) {
	t-= dt;
	printf("Evolving to t = %9.3f \n", t);

	if (t < tmin) break;
	if (f) delete f;
	f = new TF1(*g);
	for (Int_t i = 0; i < npx; i++) {
	  x = dx * i;
	// branching out at t
	// kernel
	//
	  Double_t f0 = f->Eval(x);
	  Double_t f1 = f0 * (1. - dt * integral1);
	  //
	  // branching in at t
	  // kernel
	  kint = 0.;
	  for (Int_t j = 0; j < npz; j++) {
	    z = j * dz;
	    if (z < x) {
	      kern->SetSavedPoint(j, 0.);
	      continue;
	    }
	    Double_t w1  = 0.;
	    Double_t arg = 0.;
	    // some code missing here
	    //
	    //
	    //
	    //
	    //
	    kern->SetSavedPoint(j, w1);
	  } // z loop
	  f1 +=  dt * kern->Integral(x, 1.);
	  g->SetSavedPoint(i,  f1);
	} // x loop
    }// evolve
    
  TCanvas*   c2 = new TCanvas("c2");
  c2->Divide(1,1);
  c2->SetLogy();
  g->SetMinimum(1.e-5);
  g->SetMaximum(20.);    
  g->SetLineColor(4);
  g->Draw("");
  g->GetHistogram()->SetXTitle("z");
}


Double_t p(Double_t x)
{
  Double_t y;
  if (x > 0.001) {
    y = (0.2 * (x * x + (1. - x) * (1. - x)));
    //    y = (0.2 * (1. - x * x) / (1.01 - x));
  } else {
    y = 0.;
  }
  return y;
}

Double_t pp(Double_t x)
{
  Double_t y;
  if (x > 0.001) {
    y = (0.2 * (x * x + (1. - x) * (1. - x))/x);
    //  y = (0.2 * (1. - x * x) / (1.01 - x) / x);
  } else {
    y = 0.;
  }
  return y;
}