void radial_acceptance(void)
{
// Plot acceptance of "close-packed" (2nd layer) vs radius.

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

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.2);
gStyle->SetPadBottomMargin(0.15);
// gStyle->SetOptFit(1111);
//

   char string[256];
   Int_t j,jj;
   Double_t kpi=3.14159;

   //
   // define canvas 

   TCanvas *c1 = new TCanvas("c1","Plot geometrical acceptance for CDC straw configurations",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   c1->SetGridx();
   c1->SetGridy();
   c1->Divide(1,2);
   c1->cd(1);
   c1_1->SetBorderMode(0);
   c1_1->SetFillColor(0);
   c1_1->SetGridx();
   c1_1->SetGridy();

   Double_t rmin=10;
   Double_t rmax=55;
   Double_t ymin=0;
   Double_t ymax=20;
   Double_t rstraw=0.8;

   TF1 *acc = new TF1("acc","100*sqrt(3)*[0]/x",rmin,rmax);
   acc->SetParameter(0,rstraw);
   acc->SetLineColor(2);
   // gPad->SetLogy(1);

   acc->SetTitle("");
   // acc->GetXaxis()->SetRangeUser(rmin,rmax);
   acc->GetYaxis()->SetRangeUser(ymin,ymax);
   acc->GetXaxis()->SetTitleSize(0.05);
   acc->GetXaxis()->SetLabelSize(0.04);
   acc->GetYaxis()->SetTitleSize(.05);
   acc->GetYaxis()->SetLabelSize(0.04);
   // acc->GetYaxis()->SetTitleOffset(1.5);
   acc->GetXaxis()->SetTitle("Radius (cm)");
   acc->GetYaxis()->SetTitle("Geometrical Acceptance Loss (%)");
   acc->GetXaxis()->SetNdivisions(10);
   acc->GetYaxis()->SetNdivisions(10);
   acc->Draw();

   sprintf (string,"Close-packed geometry\n");
   t1 = new TLatex(0.3,0.80,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();
   sprintf (string,"Paired layer\n");
   t1 = new TLatex(0.3,0.75,string);
   t1->SetNDC();
   t1->SetTextSize(0.035);
   t1->Draw();
   Double_t mean = acc->Integral(rmin,rmax)/(rmax-rmin);
   TLine *average = new TLine(rmin,mean,rmax,mean);
   average->SetLineWidth(2);
   average->Draw();

   // compute radius of paired layer
   
   Double_t Circ=2*kpi*rmax;
   Int_t maxtubes = floor(Circ/(2*rstraw));
   printf ("kpi=%f, Circ=%f, maxtubes=%d\n",kpi,Circ,maxtubes);
   Double_t Circ=2*kpi*rmin;
   Int_t mintubes = ceil(Circ/(2*rstraw));
   printf ("rstraw=%f, Circ=%f, mintubes=%d, maxtubes=%d\n",rstraw,Circ,mintubes,maxtubes);
   Int_t ntubes = mintubes;
   Double_t r0 = rmin;

   Int_t klayer;
   #define maxlayer 15;
   Double_t Rlayer[maxlayer];
   Double_t layer_radius[maxlayer];
   for (klayer=1;klayer<maxlayer;klayer++) {
      Double_t r0=2*rstraw*ntubes/(2*kpi);
      Double_t Dphi = 2*kpi/ntubes;
      r1 = r0*cos(Dphi/2) + sqrt(3)*rstraw;
      rtube = rstraw*(r1/r0);
      printf ("layer=%d,ntubes=%d, r1=%f, Dphi=%f, rtube=%f\n",2*klayer-1,ntubes,r1,Dphi,rtube);
      Rlayer[klayer]=r1;
      layer_radius[klayer]=rtube;
      Double_t Circ=2*kpi*r1;
      ntubes = ceil(Circ/(2*rstraw)) + 7;
   }

   TGraph *rad = new TGraph(klayer,Rlayer,layer_radius);
   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   c1_2->SetGridx();
   c1_2->SetGridy();
   rad->SetTitle("");
   rad->GetXaxis()->SetRangeUser(rmin,rmax);
   rad->GetYaxis()->SetRangeUser(0.8*rstraw,1.2*rstraw);
   rad->GetXaxis()->SetTitleSize(0.05);
   rad->GetYaxis()->SetTitleSize(0.05);
   // rad->GetYaxis()->SetTitleOffset(1.5);
   rad->GetYaxis()->SetTitle("Paired tube radius (cm)");
   rad->GetXaxis()->SetTitle("Layer radius (cm)");
   rad->SetMarkerColor(2);
   rad->SetMarkerStyle(21);
   rad->Draw("AP");

   c1->cd(2);
   c1_2->SetBorderMode(0);
   c1_2->SetFillColor(0);
   c1_2->SetGridx();
   c1_2->SetGridy();
   
   
   sprintf(string,"radial_acceptance.eps");
   c1->SaveAs(string);
   sprintf(string,"radial_acceptance.gif");
   c1->SaveAs(string);
}