void Bfit()
{
   ifstream in;
   in.open("fieldAtRandomPoints");      // CMSSW Tosca field data
   ofstream out;
   out.open("outField");   
   
    TFile *file = new TFile("ntuple.root", "recreate");
    TNtuple *ntuple = new TNtuple("ntuple","Data from file fieldAtRandomPoints"
    				 ,"r,phi,z");


// Initialize TMinuit with a maximum of 10 params
   TMinuit *gMinuit = new TMinuit(10);
   gMinuit->SetFCN(fcn);

   Double_t arglist[10];
   Int_t ierflg = 0;

   arglist[0] = 1;
   gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);


// Set starting values, step sizes and limits for parameters
   Double_t pini[7] = {4.5922,15.241,0.00088268,2.2468, 0.0031074,2.156,0.868};
   Double_t step[7] = {0.01,   0.01, 0.0001,    0.01,   0.001,    0.01, 0.01};
   Double_t llim[7] = {0.0,    0.00, 0.0000,    1.20,   0.000,    1.00, 0.25};
   Double_t ulim[7] = {7.0,   20.00, 0.0200,    2.70,   0.050,    2.40, 1.00};

// Read input data from the file fieldAtRandomPoints
   for (Int_t i=0; i<np; i++) {
     in>>r[i]>>phi[i]>>z[i]>>Br[i]>>Bphi[i]>>Bz[i];
     pnt[0]=r[i];
     pnt[1]=phi[i];
     pnt[2]=z[i];
     fBrphiz(pnt, fldE, pini);
    
     out<<i<<setprecision(5)<<setw(10)<<0.1*r[i]<<setprecision(2)<<setw(5)<<phi[i]
           <<setprecision(7)<<setw(12)<<0.1*z[i]<<"  "
        <<Br[i]<<"  "<<Bz[i]<<"  "<<fldE[0]<<"  "<<fldE[2]<<endl;
   
     ntuple->Fill(r[i], phi[i], z[i]);
   }
   
   file->Write();

   gMinuit->mnparm(0,"radius",pini[0],step[0],llim[0],ulim[0],ierflg);
   gMinuit->mnparm(1,"length",pini[1],step[1],llim[1],ulim[1],ierflg);
   gMinuit->mnparm(2,"normr ",pini[2],step[2],llim[2],ulim[2],ierflg);
   gMinuit->mnparm(3,"posir ",pini[3],step[3],llim[3],ulim[3],ierflg);
   gMinuit->mnparm(4,"normz ",pini[4],step[4],llim[4],ulim[4],ierflg);
   gMinuit->mnparm(5,"posiz ",pini[5],step[5],llim[5],ulim[5],ierflg);
   gMinuit->mnparm(6,"width ",pini[6],step[6],llim[6],ulim[6],ierflg);

// Now ready for minimization step
//   gMinuit->mnseek();
   arglist[0] = 1000;
   arglist[1] = 0.1;
   gMinuit->mnexcm("SEEK", arglist ,2,ierflg);
   gMinuit->mnexcm("SIMPLEX", arglist ,2,ierflg);
   gMinuit->mnexcm("MINOS", arglist ,2,ierflg);

// Print results
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
   gMinuit->mnprin(3,amin);

}