rendered paste bodyvoid 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);
}