int test_bfield(float r) {
  gROOT->Reset();
  float x[6], f[4];
  char  title[100];
  int   nradii;

  float radius[10];
  TH1F* h_bz[10];
  TH1F* h_br[10];

  if (r == -1) {
    nradii = 5;
    radius[0] = 0.;
    radius[1] = 24.9;
    radius[2] = 25.1;
    radius[3] = 149.;
    radius[4] = 149.1;
  }
  else if (r == -2) {
    nradii = 5;
    radius[0] = 150.;
    radius[1] = 207.;
    radius[2] = 208;
    radius[3] = 345.;
    radius[4] = 346.;
  }
  else {
    nradii = 1;
    radius[0] = r;
  }

  for (int ir=0; ir<nradii; ir++) {
    
    sprintf(title,"B(z) at r = %8f",radius[ir]);
    TH1F* h_bz[ir] = new TH1F("bz",title,100,-500,500);
    sprintf(title,"B(r) at r = %8f",radius[ir]);
    TH1F* h_br[ir] = new TH1F("br",title,100,-500,500);

    for (int i=1; i<=100; i++) {
      x[0] = radius[ir];
      x[1] = 0;
      x[2] = (i-50)*10.;
      x[3] = 0.;
      x[4] = 0.;
      x[5] = 1.;
      //    TBfield::Instance()->Bfield(f,x);
      TBfield::Instance()->Get(x,f);
      //    printf(" r, z = %f %f bx,by,bz = %f %f %f \n",x[0], x[2],f[0],f[1],f[2]);
      h_bz[ir]->SetBinContent(i,-f[2]);
      h_br[ir]->SetBinContent(i,sqrt(f[0]*f[0]+f[1]*f[1]));
    }
  }

  TCanvas* plot = new TCanvas("bf","bfield",0,0,600,800);
  if (nradii == 1) {
    plot->Divide(1,2);
    plot->cd(1);
    h_bz[0]->Draw();
    plot->cd(2);
    h_br[0]->Draw();
  }
  else {
    plot->Divide(2,5);
    for (int i=0;i<5;i++) {
      plot->cd(2*i+1);
      h_bz[i]->Draw();
      plot->cd(2*i+2);
      h_br[i]->Draw();
    }
  }
}