void read_poisson_field(void)
{
// 
// Input the ascii file produced by genr8 and generate a tree with the four vectors.
//

#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];
   Int_t j,jj,nline;

 
   #define mpts 5;

   Float_t x1=92.456;
   Float_t z1=254;
   Float_t Bdata1[mpts]={43.40,56.40,65.70,60.00,47.60};
   Float_t ydata1[mpts]={186.69,101.6,5.08,-87.63,-183.515};
   Float_t x2=187.96;
   Float_t z2=254;
   Float_t Bdata2[mpts]={36.10,43.50,48.30,45.40,32.80};
   Float_t ydata2[mpts]={184.15,100.965,19.685,-86.995,-170.815};
   Float_t x3=92.456;
   Float_t z3=203.2;
   Float_t Bdata3[mpts]={57.20,79.20,96.10,87.80,67.60};
   Float_t ydata3[mpts]={186.69,101.6,5.08,-87.63,-183.515};
   Float_t x4=187.96;
   Float_t z4=203.2;
   Float_t Bdata4[mpts]={46.50,58.30,65.50,61.80,43.80};
   Float_t ydata4[mpts]={184.15,100.965,19.685,-86.995,-170.815};
   Float_t x5=92.456;
   Float_t z5=152.4;
   Float_t Bdata5[mpts]={72.60,113.30,145.60,126.80,92.60};
   Float_t ydata5[mpts]={186.69,101.6,5.08,-87.63,-183.515};
   Float_t x6=187.96;
   Float_t z6=152.4;
   Float_t Bdata6[mpts]={62.70,75.90,88.30,82.60,62.00};
   Float_t ydata6[mpts]={184.15,100.965,19.685,-86.995,-170.815};

   #define npts 75552;

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

   // variables

   Float_t R[npts];
   Float_t z[npts];
   Float_t Br[npts];
   Float_t Bz[npts];
   Float_t Btot[npts];
   
   //  array to place Btot values at "fixed" x and R as a function of y.

   Float_t Bslice1[npts];
   Float_t yslice1[npts];
   Float_t Bslice2[npts];
   Float_t yslice2[npts];
   Float_t Bslice3[npts];
   Float_t yslice3[npts];
   Float_t Bslice4[npts];
   Float_t yslice4[npts];
   Float_t Bslice5[npts];
   Float_t yslice5[npts];
   Float_t Bslice6[npts];
   Float_t yslice6[npts];

   // open input file

   FILE *in1;
   sprintf(filename,"GLUEX_SOL_04_CHIMN_TEST_COIL_2_MOD_2_1500A_SF7_FRINGE.TXT");
   printf ("filename=%s\n",filename);
   in1 = fopen(filename,"r");
   nline = 0;
   
  // read and ignore header 

   for (j=0;jj<33;jj++) {
      nline++;
      fgets(string, 256, in1);
      printf ("line=%d string=%s",nline,string);
      }

  printf ("\n\n End of Header Information\n\n");

   Float_t edge_of_iron= 94;
   nline=0;
   while(fgets(string, 256, in1)!=NULL) {
      // printf ("line=%d string=%s",nline,string);
      sscanf (string,"%g %g %g %g %g ",&R[nline],&z[nline],&Br[nline],&Bz[nline],&Btot[nline]);
     z[nline] = z[nline]  - edge_of_iron;
      if (nline < 10 || nline > npts -10) {
          printf ("nline=%d, R=%g, z=%g, Br=%g, Bz=%g, Btot=%g\n",nline,R[nline],z[nline],Br[nline],Bz[nline],Btot[nline]);
          }
      nline++;
      }

    // loop over all data and select Btot for fixed values of x and z, as a function of y


   Float_t delta=0.95;
   Int_t mpts1=0;
   for (jj=0;jj<nline;jj++) {
          Float_t zz=z[jj];
          Float_t RR=R[jj];
          Float_t arg = RR*RR - x1*x1;
          if (fabs(z1-zz) < delta && arg > 0) {
              yslice1[mpts1] = sqrt(arg);
              Bslice1[mpts1] = Btot[jj];
              mpts1++;
              yslice1[mpts1] = -sqrt(arg);
              Bslice1[mpts1] = Btot[jj];
              mpts1++;
              }
   }


   Float_t delta=0.95;
   Int_t mpts2=0;
   for (jj=0;jj<nline;jj++) {
          Float_t zz=z[jj];
          Float_t RR=R[jj];
          Float_t arg = RR*RR - x2*x2;
          if (fabs(z2-zz) < delta && arg > 0) {
              yslice2[mpts2] = sqrt(arg);
              Bslice2[mpts2] = Btot[jj];
              mpts2++;
              yslice2[mpts2] = -sqrt(arg);
              Bslice2[mpts2] = Btot[jj];
              mpts2++;
              }
   }


   Float_t delta=0.95;
   Int_t mpts3=0;
   for (jj=0;jj<nline;jj++) {
          Float_t zz=z[jj];
          Float_t RR=R[jj];
          Float_t arg = RR*RR - x3*x3;
          if (fabs(z3-zz) < delta && arg > 0) {
              yslice3[mpts3] = sqrt(arg);
              Bslice3[mpts3] = Btot[jj];
              mpts3++;
              yslice3[mpts3] = -sqrt(arg);
              Bslice3[mpts3] = Btot[jj];
              mpts3++;
              }
   }


   Float_t delta=0.95;
   Int_t mpts4=0;
   for (jj=0;jj<nline;jj++) {
          Float_t zz=z[jj];
          Float_t RR=R[jj];
          Float_t arg = RR*RR - x4*x4;
          if (fabs(z4-zz) < delta && arg > 0) {
              yslice4[mpts4] = sqrt(arg);
              Bslice4[mpts4] = Btot[jj];
              mpts4++;
              yslice4[mpts4] = -sqrt(arg);
              Bslice4[mpts4] = Btot[jj];
              mpts4++;
              }
   }


   Float_t delta=0.95;
   Int_t mpts5=0;
   for (jj=0;jj<nline;jj++) {
          Float_t zz=z[jj];
          Float_t RR=R[jj];
          Float_t arg = RR*RR - x5*x5;
          if (fabs(z5-zz) < delta && arg > 0) {
              yslice5[mpts5] = sqrt(arg);
              Bslice5[mpts5] = Btot[jj];
              mpts5++;
              yslice5[mpts5] = -sqrt(arg);
              Bslice5[mpts5] = Btot[jj];
              mpts5++;
              }
   }


   Float_t delta=0.95;
   Int_t mpts6=0;
   for (jj=0;jj<nline;jj++) {
          Float_t zz=z[jj];
          Float_t RR=R[jj];
          Float_t arg = RR*RR - x6*x6;
          if (fabs(z6-zz) < delta && arg > 0) {
              yslice6[mpts6] = sqrt(arg);
              Bslice6[mpts6] = Btot[jj];
              mpts6++;
              yslice6[mpts6] = -sqrt(arg);
              Bslice6[mpts6] = Btot[jj];
              mpts6++;
              }
   }
              
   //
   TCanvas *c1 = new TCanvas("c1","Read_Poisson_Field Canvas",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

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

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

    xmin = -200;
    xmax = 200;
    ymin = 20;
    ymax = 100;
   
   TGraph *b1 = new TGraph (mpts1,yslice1,Bslice1);
   TLegend *leg1 = new TLegend(0.15,0.80,0.75,0.9);
   leg1->AddEntry(b1,"b1: x=36.4 in, z=100 in","p");
   b1->SetTitle("");
   b1->GetXaxis()->SetRangeUser(xmin,xmax);
   b1->GetYaxis()->SetRangeUser(ymin,ymax);
   b1->GetXaxis()->SetTitleSize(0.05);
   b1->GetYaxis()->SetTitleSize(0.05);
   b1->GetYaxis()->SetTitleOffset(1.5);
   b1->GetXaxis()->SetTitle("Vertical Position (cm)");
   b1->GetYaxis()->SetTitle("Btot (G)");
   b1->GetXaxis()->SetNdivisions(5);
   b1->SetMarkerColor(2);
   b1->SetMarkerStyle(20);
   b1->SetMarkerSize(0.5);
   b1->Draw("Ap");
   leg1->Draw();

   TGraph *d1 = new TGraph(mpts,ydata1,Bdata1);
   d1->SetMarkerColor(4);
   d1->SetMarkerStyle(20);
   d1->SetMarkerSize(0.7);
   d1->Draw("Samep");

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

    xmin = -200;
    xmax = 200;
    ymin = 20;
    ymax = 100;
   
   TGraph *b2 = new TGraph (mpts2,yslice2,Bslice2);
   TLegend *leg2 = new TLegend(0.15,0.80,0.75,0.9);
   leg2->AddEntry(b2,"b2: x=74 in, z=100 in","p");
   b2->SetTitle("");
   b2->GetXaxis()->SetRangeUser(xmin,xmax);
   b2->GetYaxis()->SetRangeUser(ymin,ymax);
   b2->GetXaxis()->SetTitleSize(0.05);
   b2->GetYaxis()->SetTitleSize(0.05);
   b2->GetYaxis()->SetTitleOffset(1.5);
   b2->GetXaxis()->SetTitle("Vertical Position (cm)");
   b2->GetYaxis()->SetTitle("Btot (G)");
   b2->GetXaxis()->SetNdivisions(5);
   b2->SetMarkerColor(2);
   b2->SetMarkerStyle(20);
   b2->SetMarkerSize(0.5);
   b2->Draw("Ap");
   leg2->Draw();

   TGraph *d2 = new TGraph(mpts,ydata2,Bdata2);
   d2->SetMarkerColor(4);
   d2->SetMarkerStyle(20);
   d2->SetMarkerSize(0.7);
   d2->Draw("Samep");

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

    xmin = -200;
    xmax = 200;
    ymin = 40;
    ymax = 120;
   
   TGraph *b3 = new TGraph (mpts3,yslice3,Bslice3);
   TLegend *leg3 = new TLegend(0.15,0.80,0.75,0.9);
   leg3->AddEntry(b3,"b3: x=36.4 in, z=80 in","p");
   b3->SetTitle("");
   b3->GetXaxis()->SetRangeUser(xmin,xmax);
   b3->GetYaxis()->SetRangeUser(ymin,ymax);
   b3->GetXaxis()->SetTitleSize(0.05);
   b3->GetYaxis()->SetTitleSize(0.05);
   b3->GetYaxis()->SetTitleOffset(1.5);
   b3->GetXaxis()->SetTitle("Vertical Position (cm)");
   b3->GetYaxis()->SetTitle("Btot (G)");
   b3->GetXaxis()->SetNdivisions(5);
   b3->SetMarkerColor(2);
   b3->SetMarkerStyle(20);
   b3->SetMarkerSize(0.5);
   b3->Draw("Ap");
   leg3->Draw();

   TGraph *d3 = new TGraph(mpts,ydata3,Bdata3);
   d3->SetMarkerColor(4);
   d3->SetMarkerStyle(20);
   d3->SetMarkerSize(0.7);
   d3->Draw("Samep");

   c1->cd(4);
   c1_4->SetGridx();
   c1_4->SetGridy();
   c1_4->SetBorderMode(0);
   c1_4->SetFillColor(0);

    xmin = -200;
    xmax = 200;
    ymin = 40;
    ymax = 120;
   
   TGraph *b4 = new TGraph (mpts4,yslice4,Bslice4);
   TLegend *leg4 = new TLegend(0.15,0.80,0.75,0.9);
   leg4->AddEntry(b4,"b4: x=74 in, z=80 in","p");
   b4->SetTitle("");
   b4->GetXaxis()->SetRangeUser(xmin,xmax);
   b4->GetYaxis()->SetRangeUser(ymin,ymax);
   b4->GetXaxis()->SetTitleSize(0.05);
   b4->GetYaxis()->SetTitleSize(0.05);
   b4->GetYaxis()->SetTitleOffset(1.5);
   b4->GetXaxis()->SetTitle("Vertical Position (cm)");
   b4->GetYaxis()->SetTitle("Btot (G)");
   b4->GetXaxis()->SetNdivisions(5);
   b4->SetMarkerColor(2);
   b4->SetMarkerStyle(20);
   b4->SetMarkerSize(0.5);
   b4->Draw("Ap");
   leg4->Draw();

   TGraph *d4 = new TGraph(mpts,ydata4,Bdata4);
   d4->SetMarkerColor(4);
   d4->SetMarkerStyle(20);
   d4->SetMarkerSize(0.7);
   d4->Draw("Samep");

   c1->cd(5);
   c1_5->SetGridx();
   c1_5->SetGridy();
   c1_5->SetBorderMode(0);
   c1_5->SetFillColor(0);

    xmin = -200;
    xmax = 200;
    ymin = 40;
    ymax = 200;
   
   TGraph *b5 = new TGraph (mpts5,yslice5,Bslice5);
   TLegend *leg5 = new TLegend(0.15,0.80,0.75,0.9);
   leg5->AddEntry(b5,"b5: x=36.4 in, z=60 in","p");
   b5->SetTitle("");
   b5->GetXaxis()->SetRangeUser(xmin,xmax);
   b5->GetYaxis()->SetRangeUser(ymin,ymax);
   b5->GetXaxis()->SetTitleSize(0.05);
   b5->GetYaxis()->SetTitleSize(0.05);
   b5->GetYaxis()->SetTitleOffset(1.5);
   b5->GetXaxis()->SetTitle("Vertical Position (cm)");
   b5->GetYaxis()->SetTitle("Btot (G)");
   b5->GetXaxis()->SetNdivisions(5);
   b5->SetMarkerColor(2);
   b5->SetMarkerStyle(20);
   b5->SetMarkerSize(0.5);
   b5->Draw("Ap");
   leg5->Draw();

   TGraph *d5 = new TGraph(mpts,ydata5,Bdata5);
   d5->SetMarkerColor(4);
   d5->SetMarkerStyle(20);
   d5->SetMarkerSize(0.7);
   d5->Draw("Samep");

   c1->cd(6);
   c1_6->SetGridx();
   c1_6->SetGridy();
   c1_6->SetBorderMode(0);
   c1_6->SetFillColor(0);

    xmin = -200;
    xmax = 200;
    ymin = 40;
    ymax = 120;
   
   TGraph *b6 = new TGraph (mpts6,yslice6,Bslice6);
   TLegend *leg6 = new TLegend(0.15,0.80,0.75,0.9);
   leg6->AddEntry(b6,"b6: x=74 in, z=60 in","p");
   b6->SetTitle("");
   b6->GetXaxis()->SetRangeUser(xmin,xmax);
   b6->GetYaxis()->SetRangeUser(ymin,ymax);
   b6->GetXaxis()->SetTitleSize(0.05);
   b6->GetYaxis()->SetTitleSize(0.05);
   b6->GetYaxis()->SetTitleOffset(1.5);
   b6->GetXaxis()->SetTitle("Vertical Position (cm)");
   b6->GetYaxis()->SetTitle("Btot (G)");
   b6->GetXaxis()->SetNdivisions(5);
   b6->SetMarkerColor(2);
   b6->SetMarkerStyle(20);
   b6->SetMarkerSize(0.5);
   b6->Draw("Ap");
   leg6->Draw();

   TGraph *d6 = new TGraph(mpts,ydata6,Bdata6);
   d6->SetMarkerColor(4);
   d6->SetMarkerStyle(20);
   d6->SetMarkerSize(0.7);
   d6->Draw("Samep");


   /*sprintf(string,"m1=%f\n",m1);
   t1 = new TLatex(0.4,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();*/

   //

   c1->SaveAs("read_poisson_field.pdf");
   c1->SaveAs("read_poisson_field.png");

   fclose (in1);

}