void genr8_tree(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();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
//gROOT->LoadMacro("$ROOTSYS/test/libEvent.so"); 
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];
   char filename[132];
   char name[10];
   Int_t j,jj;
   
   #define npts 8;
   Int_t ngen=10000;
   Int_t nevents=10;

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

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

   // create a tree 

   Int_t Ebeam=2500;
   sprintf(filename,"phi%d.root",Ebeam);
   TFile *tfile = new TFile(filename,"RECREATE");
   TTree *genphi = new TTree("genphi","phi generation using genr8");
   genphi->Branch("Events",&Events.Kp_px,
           "Kp_px/F:Kp_py:Kp_pz:Kp_E:Km_px:Km_py:Km_pz:Km_E:prot_px:prot_py:prot_pz:prot_E:Eb:pb:mb:mt:mKK2:mKp2:t_K:t_phi:");
   sprintf(filename,"phi%d.list",Ebeam);
   printf ("filename=%s\n",filename);
   FILE *in1;
   in1 = fopen(filename,"r");

   while(fgets(string, 256, in1)!=NULL) {
      // printf ("string=%s",string);
      sscanf (string,"%d %d %d",&runno,&eventno,&nparticles);
      // printf ("runno=%d, eventno=%d, nparticles=%d\n",runno,eventno,nparticles);

      // decode lines with particle information

      for (j=0;j<nparticles;j++) {
         fgets(string, 256,in1);
         sscanf (string,"%d%s%f",&partid,name,&mass);
         // printf ("partid=%d, name=%s, mass=%f\n",partid,name,mass);
         fgets(string, 256,in1);
         sscanf (string,"%d%f%f%f%f",&q,&px,&py,&pz,&E);
         printf ("runno=%d, eventno=%d, nparticles=%d, partid=%d, name=%s, mass=%f, q=%d, px=%f, py=%f, pz=%f, E=%f\n",runno,eventno,nparticles,partid,name,mass,q,px,py,pz,E);
         if (partid == 1) {
           Events.Kp_px = px;
           Events.Kp_py = py;
           Events.Kp_pz = pz;
           Events.Kp_E = E;
         }
         if (partid == 2) {
           Events.Km_px = px;
           Events.Km_py = py;
           Events.Km_pz = pz;
           Events.Km_E = E;
         }
         if (partid == 3) {
           Events.prot_px = px;
           Events.prot_py = py;
           Events.prot_pz = pz;
           Events.prot_E = E;
         }
      }
      Events.Eb = Events.Kp_pz+Events.Km_pz+Events.prot_pz;
      Events.mb = 0;
      Events.pb = Events.Eb;
      Events.mt = 0.938;
      Float_t phi_px = Events.Kp_px + Events.Km_px;
      Float_t phi_py = Events.Kp_py + Events.Km_py;
      Float_t phi_pz = Events.Kp_pz + Events.Km_pz;
      Float_t phi_E = Events.Kp_E + Events.Km_E;
      Events.mKK2 = (Events.Kp_E+Events.Km_E)*(Events.Kp_E+Events.Km_E) -
                    (Events.Kp_px+Events.Km_px)*(Events.Kp_px+Events.Km_px) -
                    (Events.Kp_py+Events.Km_py)*(Events.Kp_py+Events.Km_py) -
                    (Events.Kp_pz+Events.Km_pz)*(Events.Kp_pz+Events.Km_pz);
      Events.mKp2 = (Events.Kp_E+Events.prot_E)*(Events.Kp_E+Events.prot_E) -
                    (Events.Kp_px+Events.prot_px)*(Events.Kp_px+Events.prot_px) -
                    (Events.Kp_py+Events.prot_py)*(Events.Kp_py+Events.prot_py) -
                    (Events.Kp_pz+Events.prot_pz)*(Events.Kp_pz+Events.prot_pz);
      Events.t_K = (Events.Km_E-Events.Eb)*(Events.Km_E-Events.Eb) -
                    (Events.Km_px)*(Events.Km_px) -
                    (Events.Km_py)*(Events.Km_py) -
                    (Events.Km_pz-Events.Eb)*(Events.Km_pz-Events.Eb);
      Events.t_phi = (phi_E-Events.Eb)*(phi_E-Events.Eb) -
                    (phi_px)*(phi_px) -
                    (phi_py)*(phi_py) -
                    (phi_pz-Events.Eb)*(phi_pz-Events.Eb);
      printf ("runno=%d, eventno=%d, nparticles=%d, Eb=%f, mKK2=%f, mKp2=%f t_K=%f t_phi=%f\n",runno,eventno,nparticles,Events.Eb,Events.mKK2,Events.mKp2,Events.t_K,Events.t_phi);
      genphi->Fill();
    }
    genphi->Print();
    

   //
   TCanvas *c1 = new TCanvas("c1","Genr8_Tree Canvas",200,10,700,700);
   c1->SetBorderMode(0);
   c1->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,"Hello=%s\n","hello");
   printf ("%s",string);

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


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

   //

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

   fclose (in1);
   tfile->Write();

}