void genr8_tree_eta(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];
   Double_t pi=3.14159;
   Int_t j,jj;
   
   #define npts 8;
   Int_t maxevents=10000;
   // Int_t maxevents=10;

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

   Float_t k1, k2, k3, q1, q2, q3;
   Float_t k10, k20, k30, q10, q20, q30;   // set values for nominal Pb-glass

   // define structure

   struct Events_t {

      // generated values
      Float_t Kg1_px;
      Float_t Kg1_py;
      Float_t Kg1_pz;
      Float_t Kg1_E;
      Float_t Kg2_px;
      Float_t Kg2_py;
      Float_t Kg2_pz;
      Float_t Kg2_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 mg1g22;
      Float_t t;
      Float_t xtgt;
      Float_t ytgt;
      Float_t ztgt;

      // reconstructed values
      Float_t insertk1;
      Float_t insertk2;
      Float_t insertk3;
      Float_t insertq1;
      Float_t insertq2;
      Float_t insertq3;
      Float_t RKg1_px;
      Float_t RKg1_py;
      Float_t RKg1_pz;
      Float_t RKg1_E;
      Float_t RKg2_px;
      Float_t RKg2_py;
      Float_t RKg2_pz;
      Float_t RKg2_E;
      Float_t Rprot_px;
      Float_t Rprot_py;
      Float_t Rprot_pz;
      Float_t Rprot_E;
      Float_t REb;
      Float_t Rpb;
      Float_t Rmb;
      Float_t Rmt;
      Float_t Rmg1g22;
      Float_t Rt;
      Float_t Rxtgt;
      Float_t Rytgt;
      Float_t Rztgt;
      };
   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=11;
   sprintf(filename,"eta_%dgev.root",Ebeam);
   TFile *tfile = new TFile(filename,"RECREATE");
   TTree *geneta = new TTree("geneta","eta generation using genr8");
   geneta->Branch("Events",&Events.Kg1_px,
           "Kg1_px/F:Kg1_py:Kg1_pz:Kg1_E:Kg2_px:Kg2_py:Kg2_pz:Kg2_E:prot_px:prot_py:prot_pz:prot_E:Eb:pb:mb:mt:mg1g22:t:xtgt:ytgt:ztgt:insertk1:insertk2:insertk3:insertq1:insertq2:insertq3:RKg1_px/F:RKg1_py:RKg1_pz:RKg1_E:RKg2_px:RKg2_py:RKg2_pz:RKg2_E:Rprot_px:Rprot_py:Rprot_pz:Rprot_E:REb:Rpb:Rmb:Rmt:Rmg1g22:Rt:Rxtgt:Rytgt:Rztgt");
   // sprintf(filename,"eta_mixed_%dgev.list",Ebeam);  // used for mixed events
   sprintf(filename,"eta_%dgev.list",Ebeam);
   printf ("filename=%s\n",filename);
   FILE *in1;
   in1 = fopen(filename,"r");

   jj=0;
   while(fgets(string, 256, in1)!=NULL && jj < maxevents) {
      jj++;
      // 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.Kg1_px = px;
           Events.Kg1_py = py;
           Events.Kg1_pz = pz;
           Events.Kg1_E = E;
         }
         if (partid == 2) {
           Events.Kg2_px = px;
           Events.Kg2_py = py;
           Events.Kg2_pz = pz;
           Events.Kg2_E = E;
         }
         if (partid == 3) {
           Events.prot_px = px;
           Events.prot_py = py;
           Events.prot_pz = pz;
           Events.prot_E = E;
         }
      }

      // for mixed events use generated beam energy
      // Events.Eb = Events.Kg1_pz+Events.Kg2_pz+Events.prot_pz;
      Events.Eb = Ebeam;

      Events.mb = 0;
      Events.pb = Events.Eb;
      Events.mt = 0.938;
      Float_t eta_px = Events.Kg1_px + Events.Kg2_px;
      Float_t eta_py = Events.Kg1_py + Events.Kg2_py;
      Float_t eta_pz = Events.Kg1_pz + Events.Kg2_pz;
      Float_t eta_E = Events.Kg1_E + Events.Kg2_E;
      Events.mg1g22 = (Events.Kg1_E+Events.Kg2_E)*(Events.Kg1_E+Events.Kg2_E) -
                    (Events.Kg1_px+Events.Kg2_px)*(Events.Kg1_px+Events.Kg2_px) -
                    (Events.Kg1_py+Events.Kg2_py)*(Events.Kg1_py+Events.Kg2_py) -
                    (Events.Kg1_pz+Events.Kg2_pz)*(Events.Kg1_pz+Events.Kg2_pz);
      Events.t = (eta_E-Events.Eb)*(eta_E-Events.Eb) -
                    (eta_px)*(eta_px) -
                    (eta_py)*(eta_py) -
                    (eta_pz-Events.Eb)*(eta_pz-Events.Eb);

      TRandom1 *r = new TRandom1();   // randomize interaction vertex

      Double_t collimator_radius = 0.35/2;
      Double_t xtgt0 = 0;
      Double_t ytgt0 = 0;
      Double_t ztgt0 = 65.;
      Double_t zfcal = 625.3;
      Double_t target_length=30;
      // Double_t target_length=3;
      Double_t random_phi=r->Uniform()*2*pi;
      Double_t radius=r->Uniform();
      Events.xtgt=xtgt0 + collimator_radius*sqrt(radius)*sin(random_phi);
      Events.ytgt=ytgt0 + collimator_radius*sqrt(radius)*cos(random_phi);
      Events.ztgt=ztgt0 + target_length*(r->Uniform()-0.5); 

      // printf ("runno=%d, eventno=%d, nparticles=%d, Eb=%f, mg1g22=%f t=%f\n",runno,eventno,nparticles,Events.Eb,Events.mg1g22,Events.t);

      // fill reconstructed quantities

      // assume sigma_E = k1*sqrt(E) + c2 + c3*E;

      // photon reconstruction
      // Float_t the0 = 0.7; // 3 blocks
      Float_t the0 = 1.15; // 5 blocks (includes fiducial block)
      Float_t the1 = 3.2;
      Float_t the2a = 10;
      Float_t the2b = 11;
      Float_t the3 = 110;

      Double_t    g1_the = (180/pi)*atan2(sqrt(Events.Kg1_px*Events.Kg1_px+Events.Kg1_py*Events.Kg1_py),Events.Kg1_pz);

      Double_t e_factor=1.0;
      Double_t x_factor=1.0;
      Double_t efcal_factor=1.0;
      Double_t xfcal_factor=1.0;

       // Pb-glass nominal
       Events.insertk1 = e_factor*0.05; 
       Events.insertk2 = e_factor*0.; 
       Events.insertk3 = e_factor*0.02;
       Events.insertq1 = x_factor*0.64;
       Events.insertq2 = x_factor*0.;
       Events.insertq3 = x_factor*0.; 

       // Pb-glass Ashot parameters
       /*Events.insertk1 = e_factor*0.0481; 
       Events.insertk2 = e_factor*0.; 
       Events.insertk3 = e_factor*0.0299;
       Events.insertq1 = x_factor*0.522;
       Events.insertq2 = x_factor*0.;
       Events.insertq3 = x_factor*0.039; */

       // Rad hard glass nominal
       /*Events.insertk1 = e_factor*0.06; 
       Events.insertk2 = e_factor*0.; 
       Events.insertk3 = e_factor*0.024;
       Events.insertq1 = x_factor*0.64;
       Events.insertq2 = x_factor*0.;
       Events.insertq3 = x_factor*0.; */

       // PbWO4 Ashot
       /*Events.insertk1 = e_factor*0.0117;
       Events.insertk2 = e_factor*0.018;
       Events.insertk3 = e_factor*0.0115;
       Events.insertq1 = x_factor*0.274;
       Events.insertq2 = x_factor*0.;
       Events.insertq3 = x_factor*0.;*/


       // Shashlyk 110x110
       /*Events.insertk1 = e_factor*0.032;
       Events.insertk2 = e_factor*0.003;
       Events.insertk3 = e_factor*0.016;
       Events.insertq1 = x_factor*1.54;
       Events.insertq2 = x_factor*0.;
       Events.insertq3 = x_factor*0.31;*/


       // Shashlyk 55x55
       /*Events.insertk1 = e_factor*0.032;
       Events.insertk2 = e_factor*0.003;
       Events.insertk3 = e_factor*0.016;
       Events.insertq1 = x_factor*1.54/2;
       Events.insertq2 = x_factor*0.;
       Events.insertq3 = x_factor*0.31/2;*/


       // Pb-glass nominal
       k10 = 0.05;
       k20 = 0;
       k30 = 0.02;
       q10 = 0.64;
       q20 = 0;
       q30 = 0;

        // Pb-glass Ashot parameters
        /*k10 = 0.0481;
        k20 = 0.;
        k30 = 0.0299; 
        q10 = 0.522;
        q20 = 0.;
        q30 = 0.039; */

     
      
      if (g1_the < the1) {
          k1 = Events.insertk1;
          k2 = Events.insertk2;
          k3 = Events.insertk3;
          q1 = Events.insertq1;
          q2 = Events.insertq2;
          q3 = Events.insertq3;
       }
       elseif (g1_the < the2b) {
         k1 = k10*efcal_factor;
         k2 = k20*efcal_factor;
         k3 = k30*efcal_factor;
         q1 = q10*xfcal_factor;
         q2 = q20*xfcal_factor;
         q3 = q30*xfcal_factor;

       }
       elseif (g1_the < the3) {
          // k1 = 0.;
          k1 = 0.054;
          k2 = 0.;
          k3 = 0.023;
          // q1 = 0.;
          q1 = 0.5;
          q2 = 0.;
          q3 = 0.;
       }

      Double_t mean_E = Events.Kg1_E;
      Double_t sigma_E = sqrt(k1*sqrt(mean_E)*k1*sqrt(mean_E) + k2*k2 + k3*mean_E*k3*mean_E);
      if (g1_the > the0 && g1_the < the2a || g1_the > the2b && g1_the < the3)) {
          Events.RKg1_E=   r->Gaus(mean_E,sigma_E);
       }
      else {
         Events.RKg1_E = 0;
      }

      // printf ("\n1)  E=%f px=%f py=%f pz=%f \n",Events.RKg1_E,Events.RKg1_px,Events.RKg1_py,Events.RKg1_pz);
      Events.RKg1_px = Events.Kg1_px*Events.RKg1_E/mean_E;
      Events.RKg1_py = Events.Kg1_py*Events.RKg1_E/mean_E;
      Events.RKg1_pz = Events.Kg1_pz*Events.RKg1_E/mean_E;
      // printf ("2)  E=%f px=%f py=%f pz=%f \n",Events.RKg1_E,Events.RKg1_px,Events.RKg1_py,Events.RKg1_pz);

      Float_t Dist = zfcal-Events.ztgt;   
      Float_t Dist0 = zfcal-ztgt0;
      if (Events.RKg1_pz > 0) {
          Double_t mean_x = Dist*Events.RKg1_px/Events.RKg1_pz+Events.xtgt;
          Double_t sigma_x = sqrt((q1/sqrt(mean_E))*(q1/sqrt(mean_E)) + (q2/mean_E)*(q2/mean_E) + q3*q3);
          Double_t mean_y = Dist*Events.RKg1_py/Events.RKg1_pz+Events.ytgt;
          Double_t sigma_y = sqrt((q1/sqrt(mean_E))*(q1/sqrt(mean_E)) + (q2/mean_E)*(q2/mean_E) + q3*q3);
          // printf ("3) x=%f y=%f z=%f \n",mean_x,mean_y,Dist);

          // use center of target in the reconstruction
          Double_t Rx = r->Gaus(mean_x,sigma_x);
          Double_t Ry = r->Gaus(mean_y,sigma_y);
          Double_t Rr = sqrt(Rx*Rx + Ry*Ry);
          Events.RKg1_px = (Events.RKg1_pz/Dist0) * Rx;
          Events.RKg1_py = (Events.RKg1_pz/Dist0) * Ry;
          Events.RKg1_pz = sqrt(Events.RKg1_E*Events.RKg1_E - Events.RKg1_px*Events.RKg1_px - Events.RKg1_py*Events.RKg1_py);
          Double_t Rcost = cos(atan2(sqrt(Events.RKg1_px*Events.RKg1_px + Events.RKg1_py*Events.RKg1_py),Events.RKg1_pz));
          //printf ("\n1)  E=%f Rx=%f Ry=%f Rr=%f Rcost=%f  px=%f py=%f pz=%f\n",Events.RKg1_E,Rx-mean_x,Ry-mean_y,Rr,Rcost,Events.RKg1_px,Events.RKg1_py,Events.RKg1_pz);
          /*Events.RKg1_pz = Events.RKg1_E*Rcost;
          Events.RKg1_px = (Events.RKg1_pz/Dist0) * Rx;
          Events.RKg1_py = (Events.RKg1_pz/Dist0) * Ry;
          printf ("2)  E=%f Rx=%f Ry=%f Rr=%f Rcost=%f px=%f py=%f pz=%f\n",Events.RKg1_E,Rx,Ry,Rr,Rcost,Events.RKg1_px,Events.RKg1_py,Events.RKg1_pz);*/
      }
      else {
         Events.RKg1_px = 0;
         Events.RKg1_py = 0;
         Events.RKg1_pz = 0;
         Events.RKg1_E = 0;
      }

     // printf ("E1=%g, sigma=%g, x=%g, sigmax=%g, y=%g, sigmay=%g\n",mean_E,sigma_E,mean_x,sigma_x,mean_y,sigma_y);

      Double_t g2_the = (180/pi)*atan2(sqrt(Events.Kg2_px*Events.Kg2_px+Events.Kg2_py*Events.Kg2_py),Events.Kg2_pz);
      if (g2_the < the1) {
          k1 = Events.insertk1;
          k2 = Events.insertk2;
          k3 = Events.insertk3;
          q1 = Events.insertq1;
          q2 = Events.insertq2;
          q3 = Events.insertq3;
       }
       elseif (g2_the < the2b) {
         k1 = k10*efcal_factor;
         k2 = k20*efcal_factor;
         k3 = k30*efcal_factor;
         q1 = q10*xfcal_factor;
         q2 = q20*xfcal_factor;
         q3 = q30*xfcal_factor;
       }
       elseif (g1_the < the3) {
          // k1 = 0.;
          k1 = 0.054;
          k2 = 0.;
          k3 = 0.023;
          // q1 = 0;
          q1 = 0.5;
          q2 = 0.;
          q3 = 0.;
       }

      Double_t mean_E = Events.Kg2_E;
      Double_t sigma_E = sqrt(k1*sqrt(mean_E)*k1*sqrt(mean_E) + k2*k2 + k3*mean_E*k3*mean_E);
      if (g2_the > the0 && g2_the < the2a || g2_the > the2b && g2_the < the3)) {
          Events.RKg2_E=   r->Gaus(mean_E,sigma_E);
       }
      else {
         Events.RKg2_E = 0;
      }
      Events.RKg2_px = Events.Kg2_px*Events.RKg2_E/mean_E;
      Events.RKg2_py = Events.Kg2_py*Events.RKg2_E/mean_E;
      Events.RKg2_pz = Events.Kg2_pz *Events.RKg2_E/mean_E;  

      if (Events.RKg2_pz > 0) {
          Double_t mean_x = Dist*Events.RKg2_px/Events.RKg2_pz+Events.xtgt;
          Double_t sigma_x = sqrt((q1/sqrt(mean_E))*(q1/sqrt(mean_E)) + (q2/mean_E)*(q2/mean_E) + q3*q3);
          Double_t mean_y = Dist*Events.RKg2_py/Events.RKg2_pz+Events.ytgt;
          Double_t sigma_y = sqrt((q1/sqrt(mean_E))*(q1/sqrt(mean_E)) + (q2/mean_E)*(q2/mean_E) + q3*q3);
          // printf ("3) x=%f y=%f z=%f \n",mean_x,mean_y,Dist);
          Events.RKg2_px = (Events.RKg2_pz/Dist0) * r->Gaus(mean_x,sigma_x);
          Events.RKg2_py = (Events.RKg2_pz/Dist0) * r->Gaus(mean_y,sigma_y);
          Events.RKg2_pz = sqrt(Events.RKg2_E*Events.RKg2_E - Events.RKg2_px*Events.RKg2_px - Events.RKg2_py*Events.RKg2_py);
          // printf ("3)  E=%f px=%f py=%f pz=%f \n",Events.RKg2_E,Events.RKg2_px,Events.RKg2_py,Events.RKg2_pz);
      }
      else {
         Events.RKg2_px = 0;
         Events.RKg2_py = 0;
         Events.RKg2_pz = 0;
         Events.RKg2_E = 0;
      }


     // printf ("E2=%g, sigma=%g, x=%g, sigmax=%g, y=%g, sigmay=%g\n\n",mean_E,sigma_E,mean_x,sigma_x,mean_y,sigma_y);

      // proton 
      Double_t k4 = 0.02;
      Double_t k5 = 0.02;
      
      Double_t mean_p = sqrt(Events.prot_px*Events.prot_px+Events.prot_py*Events.prot_py+Events.prot_pz*Events.prot_pz);
      Double_t sigma_p = k4*mean_p + k5*mean_p*mean_p;
      Double_t Rprot_p =   r->Gaus(mean_p,sigma_p);
      Events.Rprot_px = Events.prot_px*Rprot_p/mean_p;
      Events.Rprot_py = Events.prot_py*Rprot_p/mean_p;
      Events.Rprot_pz = Events.prot_pz*Rprot_p/mean_p;
      Events.Rprot_E = sqrt(Rprot_p*Rprot_p + Events.mt*Events.mt);

      // beam energy
      Double_t kk1 = 0;
      Double_t kk2 = 0.;
      Double_t kk3 = 0.001;

      Double_t mean_E = Events.Eb;
      Double_t sigma_E = kk1*sqrt(mean_E) + kk2 + kk3*mean_E;
      Events.REb =   r->Gaus(mean_E,sigma_E);
      Events.Rpb = Events.REb;
      Events.Rmb = Events.mb ;
      Events.Rmt = Events.mt;

      Float_t Reta_px = Events.RKg1_px + Events.RKg2_px;
      Float_t Reta_py = Events.RKg1_py + Events.RKg2_py;
      Float_t Reta_pz = Events.RKg1_pz + Events.RKg2_pz;
      Float_t Reta_E = Events.RKg1_E + Events.RKg2_E;
      Events.Rmg1g22 = (Events.RKg1_E+Events.RKg2_E)*(Events.RKg1_E+Events.RKg2_E) -
                    (Events.RKg1_px+Events.RKg2_px)*(Events.RKg1_px+Events.RKg2_px) -
                    (Events.RKg1_py+Events.RKg2_py)*(Events.RKg1_py+Events.RKg2_py) -
                    (Events.RKg1_pz+Events.RKg2_pz)*(Events.RKg1_pz+Events.RKg2_pz);
      Events.Rt = (Reta_E-Events.REb)*(Reta_E-Events.REb) -
                    (Reta_px)*(Reta_px) -
                    (Reta_py)*(Reta_py) -
                    (Reta_pz-Events.REb)*(Reta_pz-Events.REb);

      // target 
      Events.Rxtgt=xtgt0;
      Events.Rytgt=ytgt0;
      Events.Rztgt=ztgt0;    // nominal target center

      geneta->Fill();
    }
    geneta->Print();
    

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

}