void dalitz(void)
{
// 
// plot the number of fits to a cdc segment that result in prob > prob_cut (currently set to 0.01).
//

#include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
gStyle->SetOptFit(kFALSE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   Int_t j,jj;
   
   #define npts 8;
   Int_t ngen=10000;

   Double_t xmin=2;
   Double_t xmax=4;
   Double_t ymin=0.5;
   Double_t ymax=3;

   //
   TCanvas *c1 = new TCanvas("c1","c1 Dalitz Plot",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   // plot pulse function

   Double_t mt = 0.938272;   // target units are GeV
   Double_t mb = 0;          // beam
   Double_t m1 = 0.938272;   // product 1
   Double_t m2 = 0.497614;   // product 2
   Double_t m3 = 0.497614;   // product 3
   Double_t Eb = 3;         // beam energy
   Double_t s = mt*mt + mb*mb + 2*Eb*mt; // total c.m. energy squared
   sprintf (string,"Target mass=%f\nBeam mass=%f\nProduct 1 mass=%f\nProduct 2 mass=%f\nProduct 3 mass=%f\nBeam energy=%f\ns=%f\n",mt,mb,m1,m2,m3,Eb,s);
   printf ("%s",string);

   xmin = (m1+m2)*(m1+m2);
   xmax = (sqrt(s)-m3)*(sqrt(s)-m3);

   TF1 *Dcontour_max = new TF1("Dcontour_max",Dcontour_max,xmin,xmax,4);
   Dcontour_max->SetParameters(s,m1,m2,m3);
   TF1 *Dcontour_min = new TF1("Dcontour_min",Dcontour_min,xmin,xmax,4);
   Dcontour_min->SetParameters(s,m1,m2,m3);
   TGraph *contour = new TGraph (Dcontour_max);

   contour->SetTitle("");
   contour->GetXaxis()->SetRangeUser(xmin,xmax);
   contour->GetYaxis()->SetRangeUser(ymin,ymax);
   contour->GetXaxis()->SetTitleSize(0.05);
   contour->GetYaxis()->SetTitleSize(0.05);
   contour->GetYaxis()->SetTitleOffset(1.5);
   contour->GetXaxis()->SetTitle("m_{Kp}^{2}");
   contour->GetYaxis()->SetTitle("m_{KK}^{2}");
   contour->GetXaxis()->SetNdivisions(5);
   contour->SetMarkerColor(1);
   contour->SetMarkerStyle(21);
   contour->Draw("AL");

   // Dcontour_max->SetLineColor(1);
   // Dcontour_max->Draw("same");
   Dcontour_min->SetLineColor(1);
   Dcontour_min->Draw("same");
   /*sprintf(string,"Dcontour_max: );
   printf("string=%s\n",string);*/

   Double_t mphi=1.020;
   Double_t mphi2=mphi*mphi;
   TLine *lphi = new TLine(xmin,mphi2,xmax,mphi2);
   lphi->SetLineColor(2);
   lphi->Draw("same");
   Double_t mtheta=1.54;
   Double_t mtheta2=mtheta*mtheta;
   TLine *ltheta = new TLine(mtheta2,ymin,mtheta2,ymax);
   ltheta->SetLineColor(2);
   ltheta->Draw("same");

   sprintf(string,"m1=%f\n",m1);
   t1 = new TLatex(0.4,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   // t1->Draw();

   TCanvas *c2 = new TCanvas("c2","c2 Dalitz Plot",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   c2->SetGridx();
   c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   TF1 *Pm12s = new TF1("Pm12s",Pm12s,sqrt(xmin),sqrt(xmax),4);
   sprintf (string,"s=%f\nProduct 1 mass=%f\nProduct 2 mass=%f\nProduct 3 mass=%f\n",s,m1,m2,m3);
   printf ("%s",string);
   Pm12s->SetParameters(s,m1,m2,m3);
   Pm12s->Draw();

   //

   c1->SaveAs("dalitz.eps");
   c1->SaveAs("dalitz.gif");

}
Double_t Dcontour_max (Double_t *x, Double_t *par) {

   Double_t Ms=par[0];
   Double_t m1=par[1];
   Double_t m2=par[2];
   Double_t m3=par[3];
   Double_t M=sqrt(Ms);
   Double_t m12s=x[0];

   Double_t M=sqrt(Ms);
   Double_t m1s=m1*m1;
   Double_t m2s=m2*m2;
   Double_t m3s=m3*m3;
   Double_t m12=sqrt(m12s);

   char string[256];

   if (m12s<(m1+m2)*(m1+m2) | m12s>(M-m3)*(M-m3)) {
      // sprintf (string,"m12=%f, M=%f, m1=%f, m2=%f, m3=%f\n",m12,M,m1,m2,m3);
      // printf ("string=%s",string);
      return -1;
      }

   Double_t E2 = (m12s-m1s+m2s)/(2*m12);
   Double_t E3 = (Ms-m12s-m3s)/(2*m12);
   Double_t m23s = (E2+E3)*(E2+E3)-(sqrt(E2*E2-m2s) - sqrt(E3*E3-m3s))*(sqrt(E2*E2-m2s) - sqrt(E3*E3-m3s));
   // Double_t m23s = (E2+E3)*(E2+E3)-(sqrt(E2*E2-m2s) + sqrt(E3*E3-m3s))*(sqrt(E2*E2-m2s) + sqrt(E3*E3-m3s));
   
   /*sprintf (string,"Dcontour=%f b=%f c=%f DeE=%f E=%f\n",pde,b,c,DeE,E);
   printf ("string=%s",string);*/

   return m23s;
}
Double_t Dcontour_min (Double_t *x, Double_t *par) {

   Double_t Ms=par[0];
   Double_t m1=par[1];
   Double_t m2=par[2];
   Double_t m3=par[3];
   Double_t M=sqrt(Ms);
   Double_t m12s=x[0];

   Double_t M=sqrt(Ms);
   Double_t m1s=m1*m1;
   Double_t m2s=m2*m2;
   Double_t m3s=m3*m3;
   Double_t m12=sqrt(m12s);

   char string[256];

   if (m12s<(m1+m2)*(m1+m2) | m12s>(M-m3)*(M-m3)) {
      // sprintf (string,"m12=%f, M=%f, m1=%f, m2=%f, m3=%f\n",m12,M,m1,m2,m3);
      // printf ("string=%s",string);
      return -1;
      }

   Double_t E2 = (m12s-m1s+m2s)/(2*m12);
   Double_t E3 = (Ms-m12s-m3s)/(2*m12);
   // Double_t m23s = (E2+E3)*(E2+E3)-(sqrt(E2*E2-m2s) - sqrt(E3*E3-m3s))*(sqrt(E2*E2-m2s) - sqrt(E3*E3-m3s));
   Double_t m23s = (E2+E3)*(E2+E3)-(sqrt(E2*E2-m2s) + sqrt(E3*E3-m3s))*(sqrt(E2*E2-m2s) + sqrt(E3*E3-m3s));
   
   /*sprintf (string,"Dcontour=%f b=%f c=%f DeE=%f E=%f\n",pde,b,c,DeE,E);
   printf ("string=%s",string);*/

   return m23s;
}
Double_t Pm12s (Double_t *x, Double_t *par) {

   // generate phase space distribution in m12s=m12*m12
   // See Hagedorn Eq. 7-43

   char string[256];

   Double_t mstars=par[0];
   Double_t m1=par[1];
   Double_t m2=par[2];
   Double_t m3=par[3];
   Double_t m12=x[0];
   Double_t mstar=sqrt(mstars);

   // sprintf(string,"m1=%f, m2=%f, m3=%f, m12=%f, mstar=%f\n",m1,m2,m3,m12,mstar);
   // printf ("Pm12s: %s",string);

   Double_t xin[1];
   Double_t parin[2];

   xin[0] = mstar;
   parin[0] = m12;
   parin[1] = m3;
   Double_t pm12 = sqrt(p2mstar(xin,parin));

   xin[0] = m12;
   parin[0] = m1;
   parin[1] = m2;
   Double_t pmj = sqrt(p2mstar(xin,parin));

   // sprintf(string,"pm12=%f, pmj=%f\n",pm12,pmj);
   // printf ("Pm12s: %s",string);

   Double_t distr = pm12 * pmj;

   return distr;

}
Double_t p2mstar (Double_t *x, Double_t *par) {

   // See Hagedorn Eq. 7-43

   Double_t m1=par[0];
   Double_t m2=par[1];
   Double_t mstar=x[0];

   Double_t mstars=mstar*mstar;
   Double_t m1s=m1*m1;
   Double_t m2s=m2*m2;

   char string[256];

   if (mstar < m1+m2) {
      // sprintf (string,"mstar=%f, m1=%f, m2=%f\n",mstar,m1,m2);
      // printf ("string=%s",string);
      return -1;
      }

   Double_t p2 = (mstars-(m1+m2)*(m1+m2))*(mstars - (m1-m2)*(m1-m2))/(4*mstars);
   
   /*sprintf (string,"Dcontour=%f b=%f c=%f DeE=%f E=%f\n",pde,b,c,DeE,E);
   printf ("string=%s",string);*/

   return p2;
}