void genr8_read(void)
{
// 
// Input the ascii file produced by genr8 and generate a tree with the four vectors.
//

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

gROOT->Reset();

gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);;
// gStyle->SetOptStat(kTRUE);
gStyle->SetOptFit(kFALSE);
// gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
// gStyle->SetFillColor(0);
//

   char string[256];
   char filename[132];
   char name[10];
   char cut[132];
   Int_t j,jj;
   
   #define npts 55;


   // define structure

   struct Events_t {
      Float_t Kp_px;
      Float_t Kp_py;
      Float_t Kp_pz;
      Float_t Kp_E;
      Float_t Km_px;
      Float_t Km_py;
      Float_t Km_pz;
      Float_t Km_E;
      Float_t prot_px;
      Float_t prot_py;
      Float_t prot_pz;
      Float_t prot_E;
      Float_t Eb;
      Float_t pb;
      Float_t mb;
      Float_t mt;
      Float_t mKK2;
      Float_t mKp2;
      Float_t t_K;
      Float_t t_phi;
      };
   Events_t Events;

   // open output files

   Int_t runno, eventno, nparticles;
   Int_t partid, q;
   Float_t mass,px,py,pz,E;

  // define histograms


   Double_t xmin=2;
   Double_t xmax=3;
   Double_t ymin=0.8;
   Double_t ymax=1.8;
   Int_t nbinsx = 100;
   Int_t nbinsy = 100;

   Int_t Ebeam=2500;
   Float_t Emin=1.8;
   Float_t Emax=2.1;
   Float_t tcut=0.45;
   Float_t norm=70;

   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   // printf ("cut=%s\n",cut);
   // printf ("Add title=%s\n",strcat(cut,"title "));
   
   TH1F *H1Eb =   new TH1F("H1Eb","Beam Energy (GeV)",nbinsx,1.5,2.5);
   TH1F *H1Eb_tcut =   new TH1F("H1Eb_tcut",strcat(cut,"Beam Energy (GeV)"),nbinsx,1.5,2.5);
   TH1F *H1tK =   new TH1F("H1tK","-t_{K}(GeV^{2})",nbinsx,0,2);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH1F *H1tK_tcut =   new TH1F("H1tK_tcut",strcat(cut,"-t_{K}(GeV^{2})"),nbinsx,0,2);
   TH1F *H1tphi =   new TH1F("H1tphi","-t_{#phi} (GeV^{2}) ",nbinsx,0,2);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH1F *H1tphi_tcut =   new TH1F("H1tphi_tcut",strcat(cut,"-t_{#phi} (GeV^{2})"),nbinsx,0,2);
   TH1F *H1mKK=   new TH1F("H1mKK","mKK (GeV)",nbinsx,0.9,1.1);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH1F *H1mKK_tcut =   new TH1F("H1mKK_tcut",strcat(cut,"mKK (GeV)"),nbinsx,0.9,1.1);
   TH1F *H1mKp=   new TH1F("H1mKp","mKp (GeV)",nbinsx,1.4,1.9);
   TH1F *H1mKp_copy=   new TH1F("H1mKp","mKp (GeV)",nbinsx,1.4,1.9);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH1F *H1mKp_tcut=   new TH1F("H1mKp_tcut",strcat(cut,"mKp (GeV)"),nbinsx,1.4,1.9);
   TH1F *H1mKp_tcut_copy=   new TH1F("H1mKp_tcut",strcat(cut,"mKp (GeV)"),nbinsx,1.4,1.9);
   TH1F *H1mKK2=   new TH1F("H1mKK2","mKK^{2} (GeV^{2})",nbinsx,0.9,1.1);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH1F *H1mKK2_tcut =   new TH1F("H1mKK2_tcut",strcat(cut,"mKK^{2} (GeV^{2})"),nbinsx,0.9,1.1);
   TH1F *H1mKp2 =   new TH1F("H1mKp2","mKp^{2}(GeV^{2})",nbinsx,xmin,xmax);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH1F *H1mKp2_tcut =   new TH1F("H1mKp2_tcut",strcat(cut,"mKp^{2}(GeV^{2})"),nbinsx,xmin,xmax);
   TH2F *H2mKK2mKp2 =   new TH2F("H2mKK2mKp2","mKK vs mKp",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
   sprintf (cut,"[-t_{K}<%.2f GeV^{2}]  ",tcut);
   TH2F *H2mKK2mKp2_tcut =   new TH2F("H2mKK2mKp2_tcut",strcat(cut,"mKK vs mKp"),nbinsx,xmin,xmax,nbinsy,ymin,ymax);

   // read root tree 

   sprintf(filename,"phi%d.root",Ebeam);
   printf ("filename = %s",filename);
   TFile *tfile = new TFile(filename);
   TTree *genphi = (TTree *) gROOT->FindObject("genphi");

   genphi->SetBranchAddress("Events",&Events.Kp_px);
   
   Int_t nentries = genphi->GetEntries();
   for (Int_t j=0; j<nentries; j++){
      genphi->GetEntry(j);
      // printf ("nentries=%d, j=%d, Eb=%f, mKK2=%f, mKp2=%f t_K=%f t_phi=%f\n",nentries,j,Events.Eb,Events.mKK2,Events.mKp2,Events.t_K,Events.t_phi);

      // fill histograms between Emin and Emax

      if (Events.Eb >= Emin && Events.Eb <= Emax) {
         H1Eb->Fill(Events.Eb);
         H1tK->Fill(-Events.t_K);
         H1tphi->Fill(-Events.t_phi);
         H1mKK->Fill(sqrt(Events.mKK2));
         H1mKp->Fill(sqrt(Events.mKp2));
         H1mKp_copy->Fill(sqrt(Events.mKp2));
         H1mKK2->Fill(Events.mKK2);
         H1mKp2->Fill(Events.mKp2);
         H2mKK2mKp2->Fill(Events.mKp2,Events.mKK2);
  
      // fill histograms with tcut
         if (-Events.t_K < tcut) {
            H1Eb_tcut->Fill(Events.Eb);
            H1tK_tcut->Fill(-Events.t_K);
            H1tphi_tcut->Fill(-Events.t_phi);
            H1mKK_tcut->Fill(sqrt(Events.mKK2));
            H1mKp_tcut->Fill(sqrt(Events.mKp2));
            H1mKp_tcut_copy->Fill(sqrt(Events.mKp2));
            H1mKK2_tcut->Fill(Events.mKK2);
            H1mKp2_tcut->Fill(Events.mKp2);
            H2mKK2mKp2_tcut->Fill(Events.mKp2,Events.mKK2);
         }
       }
      }
   genphi->Print();
    

   // variables for dalitz plot

   Double_t mt = 0.938;   // target units are GeV
   Double_t mb = 0;          // beam
   Double_t m1 = 0.938;   // product 1
   Double_t m2 = 0.497;   // product 2
   Double_t m3 = 0.497;   // product 3
   Double_t Eb = Ebeam;         // beam energy
   Eb = Eb/1000;
   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);

   //
   TCanvas *c1 = new TCanvas("c1","c1 Genr8_Read Canvas",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(2,2);

   c1->cd(1);  
   c1_1->SetGridx();
   c1_1->SetGridy();
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);

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


   // TGraph *thresh = new TGraph (npts,combins,threshold);
   // TLegend *leg = new TLegend(0.15,0.70,0.35,0.9);
   // leg->AddEntry(thresh,"Combinations","p");

   sprintf(string,"Eb=%.2f GeV\n",Eb);
   if (Ebeam == 2500) sprintf(string,"Eb=%.2f-%.2f GeV\n",Emin,Emax);
   t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   // t1->Draw();
   H1Eb->SetTitle(string);
   // H1Eb->GetXaxis()->SetRangeUser(xmin,xmax);
   // H1Eb->GetYaxis()->SetRangeUser(ymin,ymax);
   H1Eb->GetXaxis()->SetTitleSize(0.05);
   H1Eb->GetYaxis()->SetTitleSize(0.05);
   H1Eb->GetYaxis()->SetTitleOffset(1.5);
   H1Eb->GetXaxis()->SetTitle("Beam Energy (GeV)");
   // H1Eb->GetYaxis()->SetTitle("m_{KK}^{2}");
   H1Eb->GetXaxis()->SetNdivisions(5);
   H1Eb->SetMarkerColor(1);
   H1Eb->SetMarkerStyle(21);
   H1Eb->Draw();

   c1->cd(2); 
   c1_2->SetGridx();
   c1_2->SetGridy();
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   H1tK->GetXaxis()->SetNdivisions(5);
   H1tK->Draw();
   c1->cd(3);
   c1_3->SetGridx();
   c1_3->SetGridy();
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);
   H1tphi->GetXaxis()->SetNdivisions(5);
   H1tphi->Draw();
   c1->cd(4);
   c1_4->SetGridx();
   c1_4->SetGridy();
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);
   H1mKK->GetXaxis()->SetNdivisions(5);
   H1mKK->Draw();


   //
   TCanvas *c3 = new TCanvas("c3","c3 Genr8_Read Canvas",200,10,700,700);
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   c3->Divide(2,2);

   c3->cd(1);  
   c3_1->SetGridx();
   c3_1->SetGridy();
   c3_1->SetBorderMode(0);
   c3_1->SetFillColor(0);

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


   // TGraph *thresh = new TGraph (npts,combins,threshold);
   // TLegend *leg = new TLegend(0.15,0.70,0.35,0.9);
   // leg->AddEntry(thresh,"Combinations","p");

   sprintf(string,"[-t_{K} < %.2f GeV^{2}]  Eb=%.2f GeV \n",tcut,Eb);
   if (Ebeam == 2500) sprintf(string,"[-t_{K} < %.2f GeV^{2}]  Eb=%.2f-%.2f GeV \n",tcut,Emin,Emax);
   t1 = new TLatex(0.5,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   // t1->Draw();
   H1Eb_tcut->SetTitle(string);
   // H1Eb_tcut->GetXaxis()->SetRangeUser(xmin,xmax);
   // H1Eb_tcut->GetYaxis()->SetRangeUser(ymin,ymax);
   H1Eb_tcut->GetXaxis()->SetTitleSize(0.05);
   H1Eb_tcut->GetYaxis()->SetTitleSize(0.05);
   H1Eb_tcut->GetYaxis()->SetTitleOffset(1.5);
   H1Eb_tcut->GetXaxis()->SetTitle("Beam Energy (GeV)");
   // H1Eb_tcut->GetYaxis()->SetTitle("m_{KK}^{2}");
   H1Eb_tcut->GetXaxis()->SetNdivisions(5);
   H1Eb_tcut->SetMarkerColor(1);
   H1Eb_tcut->SetMarkerStyle(21);
   H1Eb_tcut->Draw();

   c3->cd(2); 
   c3_2->SetGridx();
   c3_2->SetGridy();
   c3_2->SetBorderMode(0);
   c3_2->SetFillColor(0);
   H1tK_tcut->GetXaxis()->SetNdivisions(5);
   H1tK_tcut->Draw();
   c3->cd(3);
   c3_3->SetGridx();
   c3_3->SetGridy();
   c3_3->SetBorderMode(0);
   c3_3->SetFillColor(0);
   H1tphi_tcut->GetXaxis()->SetNdivisions(5);
   H1tphi_tcut->Draw();
   c3->cd(4);
   c3_4->SetGridx();
   c3_4->SetGridy();
   c3_4->SetBorderMode(0);
   c3_4->SetFillColor(0);
   H1mKK_tcut->GetXaxis()->SetNdivisions(5);
   H1mKK_tcut->Draw();


   TCanvas *c2 = new TCanvas("c2","c2 Genr8_Read Canvas",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->Divide(2,2);

   c2->cd(1);  
   c2_1->SetGridx();
   c2_1->SetGridy();
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);

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

   sprintf(string,"Eb=%.2f GeV\n",Eb);
   if (Ebeam == 2500) sprintf(string,"Eb=%.2f-%.2f GeV\n",Emin,Emax);
   t1 = new TLatex(0.45,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   // xmin = (m1+m2)*(m1+m2);
   // xmax = (sqrt(s)-m3)*(sqrt(s)-m3);
   // draw kinematic boundaries

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

   Dcontour_max->SetMarkerSize(0.5);
   Dcontour_max->SetMarkerStyle(20);
   if (Ebeam != 2500) Dcontour_max->Draw("psame");

   Dcontour_min->SetMarkerSize(0.5);
   Dcontour_min->SetMarkerStyle(20);
   Dcontour_min->SetLineColor(1);
   if (Ebeam != 2500) Dcontour_min->Draw("psame");
   /*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 *ltheta2 = new TLine(mtheta2,ymin,mtheta2,ymax);
   ltheta2->SetLineColor(2);
   ltheta2->Draw("same");

   c2->cd(2);  
   c2_2->SetGridx();
   c2_2->SetGridy();
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   H1mKK2->GetXaxis()->SetNdivisions(5);
   H1mKK2->Draw();
   c2->cd(3);  
   c2_3->SetGridx();
   c2_3->SetGridy();
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   H1mKp2->Draw();
   c2->cd(4);  
   c2_4->SetGridx();
   c2_4->SetGridy();
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);
   H1mKp->GetXaxis()->SetNdivisions(5);
   H1mKp->Draw();
   Float_t y1=0;
   Float_t y2=H1mKp->GetMaximum();
   TLine *ltheta = new TLine(mtheta,y1,mtheta,y2);
   ltheta->SetLineColor(2);
   ltheta->Draw("same");


   TCanvas *c4 = new TCanvas("c4","c4 Genr8_Read Canvas",200,10,700,700);
   c4->SetBorderMode(0);
   c4->SetFillColor(0);
   c4->Divide(2,2);

   c4->cd(1);  
   c4_1->SetGridx();
   c4_1->SetGridy();
   c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);

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

   sprintf(string,"Eb=%.2f GeV\n",Eb);
   if (Ebeam == 2500) sprintf(string,"Eb=%.2f-%.2f GeV\n",Emin,Emax);
   t1 = new TLatex(0.45,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   // xmin = (m1+m2)*(m1+m2);
   // xmax = (sqrt(s)-m3)*(sqrt(s)-m3);
   // draw kinematic boundaries

   // 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);

   Dcontour_max->SetMarkerSize(0.5);
   Dcontour_max->SetMarkerStyle(20);
   if (Ebeam != 2500) Dcontour_max->DrawCopy("psame");

   Dcontour_min->SetMarkerSize(0.5);
   Dcontour_min->SetMarkerStyle(20);
   Dcontour_min->SetLineColor(1);
   if (Ebeam != 2500) Dcontour_min->DrawCopy("psame");
   /*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");

   c4->cd(2);  
   c4_2->SetGridx();
   c4_2->SetGridy();
   c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   H1mKK2_tcut->GetXaxis()->SetNdivisions(5);
   H1mKK2_tcut->Draw();
   c4->cd(3);  
   c4_3->SetGridx();
   c4_3->SetGridy();
   c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);
   H1mKp2_tcut->Draw();
   c4->cd(4);  
   c4_4->SetGridx();
   c4_4->SetGridy();
   c4_4->SetBorderMode(0);
   c4_4->SetFillColor(0);
   H1mKp_tcut->GetXaxis()->SetNdivisions(5);
   H1mKp_tcut->Draw();
   Float_t y1=0;
   Float_t y2=H1mKp_tcut->GetMaximum();
   TLine *ltheta = new TLine(mtheta,y1,mtheta,y2);
   ltheta->SetLineColor(2);
   ltheta->Draw("same");

   // plot measured mass distribution

   Float_t data_counts[npts]={-1,0,0,0,3,10,11,18,15,36,31,39,35,45,50,43,52,51,76,74,62,51,45,43,44,49,38,41,34,34,45,48,36,35,26,29,38,36,31,24,29,21,15,22,10,12,11,16,12,8,6,2,8,4,-1};
   Float_t data_xerrors[npts];
   Float_t data_yerrors[npts];
   Float_t data_mKp[npts];
   Float_t data_delta=0.0075;
   Float_t data_offset=1.4;
   for (j=0;j<npts;j++) {
       data_mKp[j]=data_delta*(j-0.5) + data_offset;
       data_xerrors[j] = 0.0001;
       data_yerrors[j] = sqrt(data_counts[j]);
       // printf ("j=%d, data_mKp=%f \n",j,data_mKp[j]);
       // printf ("j=%d, data_mKp=%f, data_counts=%f\n",j,data_mKp[j],data_counts[j]);
       }

   TGraph *datamKp = new TGraphErrors(npts,data_mKp,data_counts,data_xerrors,data_yerrors);
   TCanvas *c5 = new TCanvas("c5","c5 Genr8_Read Canvas",200,10,700,700);
   c5->SetBorderMode(0);
   c5->SetFillColor(0);

   Float_t x1=1.4;
   Float_t x2=1.8;
   Float_t ydata1=0;
   Float_t ydata2=100;
   datamKp->SetTitle("");
   datamKp->GetXaxis()->SetRangeUser(x1,x2);
   datamKp->GetYaxis()->SetRangeUser(ydata1,ydata2);
   datamKp->GetXaxis()->SetTitleSize(0.05);
   datamKp->GetYaxis()->SetTitleSize(0.05);
   datamKp->GetYaxis()->SetTitleOffset(1.5);
   datamKp->GetXaxis()->SetTitle("m_{pK} (GeV)");
   datamKp->GetYaxis()->SetTitle("Counts");
   datamKp->GetXaxis()->SetNdivisions(5);
   // datamKp->GetYaxis()->SetNdivisions(5);
   datamKp->SetMarkerColor(2);
   datamKp->SetMarkerStyle(21);
   datamKp->Draw("Ap");

   // normalize simulation

   H1mKp_tcut_copy->Scale(norm/H1mKp_tcut_copy->GetMaximum());
   H1mKp_tcut_copy->SetLineColor(2);
   H1mKp_tcut_copy->Draw("same");
   H1mKp_copy->Scale(norm/H1mKp_copy->GetMaximum());
   H1mKp_copy->SetLineColor(4);
   // H1mKp_copy->Draw("same");

   sprintf(string,"Eb=%.2f GeV\n",Eb);
   if (Ebeam == 2500) sprintf(string,"E_{b}=%.2f-%.2f GeV\n",Emin,Emax);
   t1 = new TLatex(0.45,0.8,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"tcut=%.2f GeV\n",tcut);
   if (Ebeam == 2500) sprintf(string,"-t_{K} < %.2f GeV^{2}\n",tcut);
   t1 = new TLatex(0.45,0.74,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   
   


   if (Ebeam != 2500 ) {
      sprintf(filename,"genr8_read_c1_%d.eps",Ebeam);
      c1->SaveAs(filename);
      sprintf(filename,"genr8_read_c1_%d.gif",Ebeam);
      c1->SaveAs(filename);
      sprintf(filename,"genr8_read_c2_%d.eps",Ebeam);
      c2->SaveAs(filename);
      sprintf(filename,"genr8_read_c2_%d.gif",Ebeam);
      c2->SaveAs(filename);
      sprintf(filename,"genr8_read_c3_%d.eps",Ebeam);
      c3->SaveAs(filename);
      sprintf(filename,"genr8_read_c3_%d.gif",Ebeam);
      c3->SaveAs(filename);
      sprintf(filename,"genr8_read_c4_%d.eps",Ebeam);
      c4->SaveAs(filename);
      sprintf(filename,"genr8_read_c4_%d.gif",Ebeam);
      c4->SaveAs(filename);
      }
    else {
      Int_t E1 = (Emin+0.0005)*1000;
      Int_t E2 = (Emax+0.0005)*1000;
      Int_t icut = (tcut+0.005)*100;
      printf ("tcut=%f, icut=%d\n",tcut,icut);
      sprintf(filename,"genr8_read_c1_%d-%d_tcut%d.eps",E1,E2,icut);
      c1->SaveAs(filename);
      sprintf(filename,"genr8_read_c1_%d-%d_tcut%d.gif",E1,E2,icut);
      c1->SaveAs(filename);
      sprintf(filename,"genr8_read_c2_%d-%d_tcut%d.eps",E1,E2,icut);
      c2->SaveAs(filename);
      sprintf(filename,"genr8_read_c2_%d-%d_tcut%d.gif",E1,E2,icut);
      c2->SaveAs(filename);
      sprintf(filename,"genr8_read_c3_%d-%d_tcut%d.eps",E1,E2,icut);
      c3->SaveAs(filename);
      sprintf(filename,"genr8_read_c3_%d-%d_tcut%d.gif",E1,E2,icut);
      c3->SaveAs(filename);
      sprintf(filename,"genr8_read_c4_%d-%d_tcut%d.eps",E1,E2,icut);
      c4->SaveAs(filename);
      sprintf(filename,"genr8_read_c4_%d-%d_tcut%d.gif",E1,E2,icut);
      c4->SaveAs(filename);
      }

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