gROOT->Reset();
float Mn = 939.566;  // neutron mass
float m1 = 938.272;  // 1  -  proton 
float m2 = .511;     // 2  -  electron
float m3 = 0.;       // 3  -  neutrino
float Q,dm;
float K1min,K1max,K2min,K2max,K3min,K3max;
float flat ( float, float);
float K3p( float ), K3m( float );

void uniform(){
  double x1, x2;
  float dx, dy;
  int NMAX=10000,k;

  dm = Mn - m1;
  Q  = Mn - m1 - m2 - m3;
  K1min = 0.;
  K1max = Q*(Mn - m1 + m2 + m3)/(2*Mn);
  K2min = 0.;
  K2max = Q*(Mn - m2 + m1 + m3)/(2*Mn);
  K3min = 0.;
  K3max = Q*(Mn - m3 + m2 + m1)/(2*Mn);

  TCanvas *c2 = new TCanvas("c2", "Uniform Distribution",200,10,900,900);
  c2->Divide(1,2);
  c2->cd(1);

  dx = 0.2*(K1max - K1min);
  dy = 0.2*(K3max - K3min);

  TF2 *f1 = new TF2("f1", fun2, K1min-dx, K1max+dx, K3min-dy, K3max+dy, 0 );
  TH2D *h1 = new TH2D("h1", "flat function", 100, K1min-dx, K1max+dx, 100, K3min-dy, K3max+dy);
  h1->Eval(f1);

  h1->SetStats(kFALSE);
//   h1->SetContour(64);
  gPad->SetTheta(50.);
  gPad->SetPhi(160.);
  gStyle->SetPalette(1,0);
  h1->Draw("surf1z");
  TH2D *h2 = new TH2D("h2", "flat distribution", 100, K1min-dx, K1max+dx, 100, K3min-dy, K3max+dy);
  for (k=0; k< NMAX; k++){
    h1->GetRandom2( x1, x2 );
    h2->Fill(x1,x2);
  }
  c2->cd(2); 
  h2->SetStats(kFALSE);
  h2->Draw("box");
  return;
}

double fun2 ( double *x, double *par){

  return (double ) flat( (float ) x[0], (float) x[1] );
}

float flat(float x1, float x2){

  if ( x1 < K1min   || x1 > K1max   ) return 0;
  if ( x2 < K3m(x1) || x2 > K3p(x1) ) return 0;
  return 1.;
}

float K3p(float K1){
  float l1,l2,l3,r1,r2,res,s23;

  s23 = dm*dm - 2*Mn*K1;

  l1  = K1 - 2*dm;
  l1  = 4*Mn*Mn*K1*( 2*Mn + l1 );

  l2  = ( dm*dm - (m2 + m3)*(m2 + m3) );
  l2  = l2 - 2*Mn*K1;

  l3  = ( dm*dm - (m2 - m3)*(m2 - m3) );
  l3  = l3 - 2*Mn*K1;

  r1  = Mn*( 2*K1 + (1 + m1/Mn)*dm);
  r1  = dm*dm - r1;
 
  r2  = dm*dm + m2*m2 - m3*m3;
  r2 = r2 - 2*Mn*K1;

  if(l1 <= 0) {
//     printf("l1 <= 0 %e\n",l1);
    l1 = 0;
  }
  if(l2 <= 0) {
//     printf("l2 <= 0 %e\n",l2);
    l2 = 0;
  }
  if(l3 <= 0) {
//     printf("l3 <= 0 %e\n",l3);
    l3 = 0;
  }

  res = r1*r2/(2*s23) + sqrt(l1*l2*l3)/(2*s23);

  res = dm*(Mn + m1) + m3*m3 - m2*m2 - 2*Mn*m3 + res;
  res = res/(2*Mn) - m3;

  return res;
}

float K3m(float K1){
  float l1,l2,l3,r1,r2,res,s23;

  s23 = dm*dm - 2*Mn*K1;

  l1  = K1 - 2*dm;
  l1  = 4*Mn*Mn*K1*( 2*Mn + l1 );

  l2  = ( dm*dm - (m2 + m3)*(m2 + m3) );
  l2  = l2 - 2*Mn*K1;

  l3  = ( dm*dm - (m2 - m3)*(m2 - m3) );
  l3  = l3 - 2*Mn*K1;

  r1  = Mn*( 2*K1 + (1 + m1/Mn)*dm);
  r1  = dm*dm - r1;
 
  r2  = dm*dm + m2*m2 - m3*m3;
  r2 = r2 - 2*Mn*K1;

  if(l1 <= 0) {
//     printf("l1 <= 0 %e\n",l1);
    l1 = 0;
  }
  if(l2 <= 0) {
//     printf("l2 <= 0 %e\n",l2);
    l2 = 0;
  }
  if(l3 <= 0) {
//     printf("l3 <= 0 %e\n",l3);
    l3 = 0;
  }

  res = r1*r2/(2*s23) - sqrt(l1*l2*l3)/(2*s23);

  res = dm*(Mn + m1) + m3*m3 - m2*m2 - 2*Mn*m3 + res;
  res = res/(2*Mn) - m3;

  return res;
}