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;

void contorno(){
  float K1,dK1,dx,dy;
  int np,k,NPT=200;
  float x[200],y[200];

  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);
  printf("\nDalitz plot for beta decay\n");
  printf("Nuclei masses; Mn %6.3f  m1 %6.3f \n",Mn,m1);
  printf("Q-Value %6.3f   dM  %6.3f\n\n",Q,dm);
  printf("Kinetic Energy range for particle 1 (proton):         0. < T1 < %7.6f\n",K1max);
  printf("Kinetic Energy range for particle 2 (electron):       0. < T2 < %6.4f\n",K2max);
  printf("Kinetic Energy range for particle 3 (neutrino):       0. < T3 < %6.4f\n",K3max);
  printf("\n\n");

  TCanvas *c1 = new TCanvas("c1", "Dalitz plot for neutron beta Decay",100,10,700,500);
  
  c1->SetGridx();                  // show x grid
  c1->SetGridy();                  // show y grid
  
  dx = 0.2*(K1max - K1min);
  dy = 0.2*(K3max - K3min);
  np = NPT/2;
  dK1 = (K1max - K1min)/(np-1);
  for(k=0; k<np; k++){
    x[k] = K1min + dK1*k;
    y[k] = K3p(x[k]);
    x[NPT-k-1] = x[k];
    y[NPT-k-1] = K3m(x[k]);
  } 

  TGraph *gr1 = new TGraph(NPT,x,y);
  gr1->SetLineColor(kBlack);
  gr1->SetLineWidth(2);
  gr1->SetFillColor(kYellow);
  gr1->SetFillStyle(1001);
  gr1->SetTitle("Dalitz Plot: Enu vs Ep");
  gr1->Draw("ALF");
  gr1->GetXaxis()->SetTitle("proton Kinetic Energy");
  gr1->GetXaxis()->SetLimits( K1min-dx, K1max+dx );
  gr1->Draw("L");

  return;
}

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;
}