void solenoid_track(void)
{
// 
//  Track particle in uniform (solenoid) field and output trajectory
//
// Output root tree containing trajectory information in solenoid_track.root
//
//

// #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 10001;
   #define pi 3.14159265;
   Double_t radians=pi/180.;


   Double_t p=0.1;
   Double_t m=0.000511;
   Double_t theta_scat=10;  // theta_scat in degrees
   Double_t phi_scat=0;        // phi_scat in degrees
   Double_t B=2;
   Double_t omega_e = B*1.758e2;  // units are rad ns-1 for electrons
   Double_t me = 0.000511;
   Double_t c= 30;    // units are cm / ns

   Double_t E = sqrt(p*p + m*m);
   Double_t beta = p/E;
   Double_t p_perp = p*sin(theta_scat*radians);
   Double_t p_par = p*cos(theta_scat*radians);
   Double_t beta_par=p_par/E;
   Double_t beta_perp = p_perp/E;
   Double_t R=p_perp*100/(0.3*B);  // R is in cm
   Double_t omega = -omega_e * me / E;

   Double_t xmin=-2*R;
   Double_t xmax=2*R;
   Double_t ymin=-2*R;
   Double_t ymax=2*R;
   Double_t zmin=0;
   Double_t zmax=300;
   Double_t tmin=0;
   Double_t tmax=zmax/(beta*c);
   Double_t rperp_min=0;
   Double_t rperp_max=50;

   Double_t the_min=0;
   Double_t the_max=5*theta_scat;


   printf("Momentum=%f, Mass=%f, E=%f,Beta=%f, Theta_scat=%f, p_perp=%f, p_par=%f, omega=%f, R=%f\n",p,m,E,beta,theta_scat,p_perp,p_par,omega,R);

   // Loop over trajectory points

   Double_t t[npts];
   Double_t x[npts];
   Double_t y[npts];
   Double_t z[npts];
   Double_t vx[npts];
   Double_t vy[npts];
   Double_t vz[npts];
   Double_t theta[npts];
   Double_t phi[npts];
   Double_t theta_track[npts];
   Double_t phi_track[npts];
   Double_t rperp[npts];
   Double_t rperp_infp[npts];


   for (j=0;j<npts;j++) {
       t[j] = tmin + (tmax-tmin)*j/(npts-1);
       z[j] = beta_par * c * t[j];
       x[j] = R*sin(omega*t[j]);
       y[j] = R*cos(omega*t[j]);

       rperp_infp[j] = z[j]*tan(1.0*radians);
       rperp[j] = sqrt(x[j]*x[j]+(y[j]-R)*(y[j]-R));
       theta_track[j] = atan2(rperp[j] ,z[j])/radians;
       phi_track[j] = atan2(x[j], y[j]-R)/radians;
       // printf ("j=%d, t=%f, z=%f, x=%f, y=%f, theta=%f, phi=%f\n",j,t[j],z[j],x[j],y[j],theta_track[j],phi_track[j]);
       }

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

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->Divide(2,2);

   TGraph *zgr = new TGraph (npts+1,t,z);
   TGraph *xgr = new TGraph (npts+1,t,x);
   TGraph *ygr = new TGraph (npts+1,t,y);
   TGraph *rperp_infpgr = new TGraph (npts+1,t,rperp_infp);
   TGraph *rperp_infp_grz = new TGraph (npts+1,z,rperp_infp);
   TGraph *rperpgr = new TGraph (npts+1,t,rperp);
   TGraph *rperp_grz= new TGraph (npts+1,z,rperp);
   TGraph *thegr = new TGraph (npts+1,t,theta_track);
   TGraph *phigr = new TGraph (npts+1,t,phi_track);
   TLegend *leg = new TLegend(0.15,0.70,0.35,0.9);
   leg->AddEntry(zgr,"z","p");
   leg->AddEntry(xgr,"x","p");
   leg->AddEntry(ygr,"y","p");

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

   zgr->SetTitle("");
   zgr->GetXaxis()->SetRangeUser(tmin,tmax);
   zgr->GetYaxis()->SetRangeUser(zmin,zmax);
   zgr->GetXaxis()->SetTitleSize(0.05);
   zgr->GetYaxis()->SetTitleSize(0.05);
   zgr->GetYaxis()->SetTitleOffset(1.5);
   zgr->GetXaxis()->SetTitle("time (ns)");
   zgr->GetYaxis()->SetTitle("z (cm)");
   zgr->GetYaxis()->SetNdivisions(505);
   zgr->GetXaxis()->SetNdivisions(505);
   zgr->SetMarkerColor(1);
   zgr->SetMarkerStyle(20);
   zgr->SetMarkerSize(0.3);
   zgr->Draw("Ap");

   sprintf(string,"p=%.1f GeV, #theta=%.1f deg\n",p,theta_scat);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"mass=%.3f GeV\n",m);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);

   xgr->SetTitle("");
   xgr->GetXaxis()->SetRangeUser(tmin,tmax);
   xgr->GetYaxis()->SetRangeUser(xmin,xmax);
   xgr->GetXaxis()->SetTitleSize(0.05);
   xgr->GetYaxis()->SetTitleSize(0.05);
   xgr->GetYaxis()->SetTitleOffset(1.5);
   xgr->GetXaxis()->SetTitle("time (ns)");
   xgr->GetYaxis()->SetTitle("x (cm)");
   xgr->GetYaxis()->SetNdivisions(505);
   xgr->GetXaxis()->SetNdivisions(505);
   xgr->SetMarkerColor(1);
   xgr->SetMarkerStyle(20);
   xgr->SetMarkerSize(0.3);
   xgr->Draw("Ap");

   sprintf(string,"p=%.1f GeV, #theta=%.1f deg\n",p,theta_scat);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"mass=%.3f GeV\n",m);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


   c1->cd(3);
   c1_3->SetBorderMode(0);
   c1_3->SetFillColor(0);

   ygr->SetTitle("");
   ygr->GetXaxis()->SetRangeUser(tmin,tmax);
   ygr->GetYaxis()->SetRangeUser(ymin,ymax);
   ygr->GetXaxis()->SetTitleSize(0.05);
   ygr->GetYaxis()->SetTitleSize(0.05);
   ygr->GetYaxis()->SetTitleOffset(1.5);
   ygr->GetXaxis()->SetTitle("time (ns)");
   ygr->GetYaxis()->SetTitle("y (cm)");
   ygr->GetYaxis()->SetNdivisions(505);
   ygr->GetXaxis()->SetNdivisions(505);
   ygr->SetMarkerColor(1);
   ygr->SetMarkerStyle(20);
   ygr->SetMarkerSize(0.3);
   ygr->Draw("Ap");

   sprintf(string,"p=%.1f GeV, #theta=%.1f deg\n",p,theta_scat);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"mass=%.3f GeV\n",m);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   // leg->Draw();


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

   c2->SetGridx();
   c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->Divide(2,2);

   TLegend *leg2 = new TLegend(0.15,0.70,0.35,0.9);
   leg2->AddEntry(thegr,"#theta","p");
   leg2->AddEntry(phigr,"#phi","p");

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

   thegr->SetTitle("");
   thegr->GetXaxis()->SetRangeUser(tmin,tmax);
   thegr->GetYaxis()->SetRangeUser(the_min,the_max);
   thegr->GetXaxis()->SetTitleSize(0.05);
   thegr->GetYaxis()->SetTitleSize(0.05);
   thegr->GetYaxis()->SetTitleOffset(1.5);
   thegr->GetXaxis()->SetTitle("time (ns)");
   thegr->GetYaxis()->SetTitle("#theta of track (degrees)");
   thegr->GetYaxis()->SetNdivisions(505);
   thegr->GetXaxis()->SetNdivisions(505);
   thegr->SetMarkerColor(1);
   thegr->SetMarkerStyle(20);
   thegr->SetMarkerSize(0.3);
   thegr->Draw("Ap");

   sprintf(string,"p=%.1f GeV, #theta=%.1f deg\n",p,theta_scat);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"mass=%.3f GeV\n",m);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c2->cd(2);
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);

   rperpgr->SetTitle("");
   rperpgr->GetXaxis()->SetRangeUser(tmin,tmax);
   rperpgr->GetYaxis()->SetRangeUser(rperp_min,rperp_max);
   rperpgr->GetXaxis()->SetTitleSize(0.05);
   rperpgr->GetYaxis()->SetTitleSize(0.05);
   rperpgr->GetYaxis()->SetTitleOffset(1.5);
   rperpgr->GetXaxis()->SetTitle("time (ns)");
   rperpgr->GetYaxis()->SetTitle("Distance from beam axis (cm)");
   rperpgr->GetYaxis()->SetNdivisions(505);
   rperpgr->GetXaxis()->SetNdivisions(505);
   rperpgr->SetMarkerColor(2);
   rperpgr->SetMarkerStyle(20);
   rperpgr->SetMarkerSize(0.3);
   rperpgr->Draw("Ap");
   rperp_infpgr->SetMarkerColor(4);
   rperp_infpgr->SetMarkerStyle(20);
   rperp_infpgr->SetMarkerSize(0.3);
   rperp_infpgr->Draw("samep");

   sprintf(string,"p=%.1f GeV, #theta=%.1f deg\n",p,theta_scat);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"mass=%.3f GeV\n",m);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c2->cd(3);
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);

   rperp_grz->SetTitle("");
   rperp_grz->GetXaxis()->SetRangeUser(zmin,zmax);
   rperp_grz->GetYaxis()->SetRangeUser(rperp_min,rperp_max);
   rperp_grz->GetXaxis()->SetTitleSize(0.05);
   rperp_grz->GetYaxis()->SetTitleSize(0.05);
   rperp_grz->GetYaxis()->SetTitleOffset(1.5);
   rperp_grz->GetXaxis()->SetTitle("Distance z along beam axis (cm)");
   rperp_grz->GetYaxis()->SetTitle("Distance r from beam axis (cm)");
   rperp_grz->GetYaxis()->SetNdivisions(505);
   rperp_grz->GetXaxis()->SetNdivisions(505);
   rperp_grz->SetMarkerColor(2);
   rperp_grz->SetMarkerStyle(20);
   rperp_grz->SetMarkerSize(0.3);
   rperp_grz->Draw("Ap");
   rperp_infp_grz->SetMarkerColor(4);
   rperp_infp_grz->SetMarkerStyle(20);
   rperp_infp_grz->SetMarkerSize(0.3);
   rperp_infp_grz->Draw("samep");

   // leg2->Draw();

   sprintf(string,"p=%.1f GeV, #theta=%.1f deg\n",p,theta_scat);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"mass=%.3f GeV\n",m);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   sprintf(filename,"solenoid_track_c1_%.1f_%.1f_%.3f.pdf",p,theta_scat,m);
   printf("%s\n",filename);
   // c1->SaveAs(filename);
   sprintf(filename,"solenoid_track_c1_%.1f_%.1f_%.3f.png",p,theta_scat,m);
   printf("%s\n",filename);
   // c1->SaveAs(filename);
   sprintf(filename,"solenoid_track_c2_%.1f_%.1f_%.3f.pdf",p,theta_scat,m);
   printf("%s\n",filename);
   c2->SaveAs(filename);
   sprintf(filename,"solenoid_track_c2_%.1f_%.1f_%.3f.png",p,theta_scat,m);
   printf("%s\n",filename);
   c2->SaveAs(filename);

}