void kinematics_2body(void)
{
// 
// Input the ascii file with particle masses in beam energy to define kinematics
//    reads incident energy and masses from kinematics_2body.input
//    Original: Oct 14, 2009
//
// Output root tree containing variables of interest for the scattering process in kinematics_2body.root
//
// At the moment no histograms are produced by ths program, but could be easily added.
//

#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 title[80];
   char option[80];
   char name[10];
   Int_t j,jj;
   
   #define npts 8;
   #define pi 3.14159265;
   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 Eb;
      Float_t Pb;
      Float_t mb;
      Float_t mt;
      Float_t Escat;
      Float_t Pscat;
      Float_t beta_scat;
      Float_t mscat;
      Float_t Erecoil;
      Float_t Precoil;
      Float_t beta_recoil;
      Float_t mrecoil;
      Float_t theta_scat;   // scat in degrees
      Float_t theta_recoil; // recoil in degrees
      Float_t t;
      Float_t W;
      };
   Events_t Events;

   // open output files

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

   sprintf(filename,"kinematics_2body.input");
   // printf ("filename=%s\n",filename);
   FILE *in1;
   in1 = fopen(filename,"r");

   Float_t Eproj;
   Float_t mproj;
   Float_t mtarg;
   Float_t mrec;
   Float_t msc;
   Float_t thmin, thmax;
   Float_t Esc;
   Float_t Prec2;
   Float_t Prec;
   Float_t cosrec;
   Int_t nevents;

   if (fgets(title, 80, in1)!=NULL) strtok(title,"\n");  //remove carriage return from end of title
   if (fgets(option, 80, in1)!=NULL) strtok(option,"\n"); //remove carriage return from end of option
   fscanf (in1,"%f%f%f%f%f%f%f%d",&mproj,&mtarg,&msc,&mrec,&Eproj,&thmin,&thmax,&nevents);


      printf ("\nInput File=%s\ntitle=%s\noption=%s\nmprojectile=%f\nmtarget=%f\nmscattered=%f\nmrecoil=%f\nEproj=%f\nthmin,thmax=%f,%f\nnevents=%d\n\n",
               filename,title,option,mproj,mtarg,msc,mrec,Eproj,thmin,thmax,nevents);

   // create a tree 

   sprintf(filename,"kinematics_2body.root");
   TFile *tfile = new TFile(filename,"RECREATE");
   TTree *twobody = new TTree("twobody",title);
   twobody->Branch("Events",&Events.Eb,
           "Eb/F:Pb:mb:mt:Escat:Pscat:beta_scat:mscat:Erecoil:Precoil:beta_recoil:mrecoil:theta_scat:theta_recoil:t:W");

   // calculate threshold energy

   Float_t Ethresh = ((msc+mrec)*(msc+mrec) - mproj*mproj - mtarg*mtarg)/(2*mtarg);

   printf ("The energy threshold for this reaction is %f\n\n",Ethresh);
   if (Eproj < Ethresh) {
     printf ("*** kinematics_2body: Eproj=%f < Ethreshold=%f***\n\n",Eproj,Ethresh);
     exit (1);
     }
      
   Events.Eb = Eproj;
   Events.mb = mproj;
   Float_t Pproj = sqrt(Eproj*Eproj - mproj*mproj);
   Events.Pb = Pproj;
   Events.mt = mtarg;
   Events.mscat = msc;
   Events.mrecoil = mrec;
   Float_t s = (Events.Eb+Events.mt)*(Events.Eb+Events.mt) - Events.Pb*Events.Pb;
   // Float_t sp = (mproj*mproj + mtarg*mtarg +2*Eproj*mtarg);
   // printf ("s=%f, s2=%f\n",s,sp);
   Events.W = sqrt(s);
   Float_t A = (Eproj*mtarg)/(Eproj+mtarg) + 0.5*(mtarg*mtarg + mproj*mproj + msc*msc - mrec*mrec)/(Eproj+mtarg);

   // compute some limiting angles

   Float_t E3cm = (s+msc*msc-mrec*mrec)/(2*sqrt(s));
   Float_t p3cm = sqrt(E3cm*E3cm-msc*msc);
   Float_t E1cm = (s+mproj*mproj-mtarg*mtarg)/(2*sqrt(s));
   Float_t p1cm = sqrt(E1cm*E1cm-mproj*mproj);
   Float_t beta_3cm = p3cm/E3cm;
   Float_t beta = Events.Pb/(Events.Eb+mtarg);
   Float_t gamma = (Events.Eb+mtarg)/Events.W;
   Float_t tantheta_max;
   Float_t theta_max;
   if (beta_3cm < beta) {
       tantheta_max = (1/gamma)*beta_3cm/sqrt(beta*beta - beta_3cm*beta_3cm);       
       theta_max = atan(tantheta_max)*180/pi;
       }
   else {
       tantheta_max = 0;
       theta_max = 180;
   }

   // printf ("p3cm=%f, beta_3cm=%f, beta=%f, gamma=%f, tantheta_max=%f, theta_max=%f\n",p3cm,beta_3cm,beta ,gamma,tantheta_max,theta_max);

   Float_t arg1 = (mproj*mproj-msc*msc-mtarg*mtarg+mrec*mrec)/(2*sqrt(s));
   Float_t t0 = arg1*arg1 - (p1cm-p3cm)*(p1cm-p3cm);
   Float_t t1 = arg1*arg1 - (p1cm+p3cm)*(p1cm+p3cm);
   if (thmin<=0 && thmax <=0) {
      thmin = 0;
      thmax = theta_max;
   }

   printf("Angle kinematic limits thmin=%f, thmax=%f\n",thmin,thmax);
   printf("t kinematic limits are t0=%f, t1=%f\n\n",t0,t1);

   for (j=0;j<nevents+1;j++) {
 
      Events.theta_scat = thmin + j*(thmax-thmin)/nevents;
      Float_t costhe = cos(Events.theta_scat*pi/180);
      Float_t alpha = Pproj*costhe/(Eproj+mtarg);
      Float_t arg = 1 - (msc/A)*(msc/A)*(1-alpha*alpha);
      if (arg >= 0) {
         Esc = A*(1+alpha*sqrt(arg))/(1-alpha*alpha);
         }
      else {
         Esc = msc;
      }
      Psc = sqrt(Esc*Esc - msc*msc);
      Prec2 = Psc*Psc - 2*Pproj*Psc*costhe+Pproj*Pproj;
      if (Prec2 >= 0) {
         Prec = sqrt(Prec2);
         cosrec = (Pproj-Psc*costhe)/Prec;
         }
      else {
         Prec = 0;
         cosrec = 1;
      }
      Float_t therec = acos(cosrec);

      Events.Escat = Esc;
      Events.Pscat = Psc;
      Events.beta_scat = Psc/Esc;
      Events.Erecoil = sqrt(Prec*Prec+mrec*mrec);
      Events.Precoil = Prec;
      Events.beta_recoil = Prec/Events.Erecoil;
      Events.theta_recoil = therec*180/pi;
      Events.t = -2*mrec*(Events.Erecoil-mrec);
   
      // printf ("event=%d, theta_scat=%f, t=%f\n",j,Events.theta_scat,Events.t);
     
      twobody->Fill();
   }
      twobody->Print();

   //
   /*TCanvas *c1 = new TCanvas("c1","Kinematics_2body 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("kinematics_2body.eps");
   c1->SaveAs("kinematics_2body.gif");*/

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

}