C++ TMinuit::SetFCN方法代码示例

//______________________________________________________________________________
vector< vector<double> > cSpline(int nPoints, int npar, vector <double> xData, vector <double> yData, vector <double> yErrorData, double stepSpline, double start, double step)
{
	//Populate the global variables//-> DONE IN THE MAIN
	xData_GLOB = xData;//-> DONE IN THE MAIN
	yData_GLOB = yData;//-> DONE IN THE MAIN
	yErrorData_GLOB = yErrorData;//-> DONE IN THE MAIN

	//Initialize Minuit
	vector<double> vstart, vstep;
	for(int i=0; i<npar; i++) 	//set starting values and step sizes for parameters
	{
		vstart.push_back(start);
		vstep.push_back(step);
	}

	TMinuit *myMinuit = new TMinuit(npar);  //initialize TMinuit with a maximum of npar (5) -> PASSED AS ARGUMENT
	myMinuit->SetFCN(fcn);//-> DONE IN THE MAIN
	myMinuit->SetPrintLevel(-1);//No output: -1, output:1//-> DONE IN THE MAIN

	double arglist[10];//-> DONE IN THE MAIN
	int ierflg = 0;//-> DONE IN THE MAIN
	arglist[0] = 1;//-> DONE IN THE MAIN
	myMinuit->mnexcm("SET ERR", arglist, 1, ierflg);//-> DONE IN THE MAIN

	for (int i = 0; i < npar; i++) {//-> DONE IN THE MAIN
		stringstream ss;//-> DONE IN THE MAIN
		ss<<"a"<<i;//-> DONE IN THE MAIN
		myMinuit->mnparm(i, ss.str().c_str(), vstart.at(i), vstep.at(i), 0, 0, ierflg);//-> DONE IN THE MAIN
	}//-> DONE IN THE MAIN

	//Perform the Minuit fit
	arglist[0] = 500;//-> DONE IN THE MAIN
	arglist[1] = 1.;//-> DONE IN THE MAIN
	myMinuit->mnexcm("MIGRAD", arglist, 2, ierflg); //minimization

	//Retrieve best-fit parameters
	vector< double > bestFitParams, e_bestFitParams;
	for(int i=0; i<npar; i++)
	{
		double par, epar;
		myMinuit->GetParameter(i, par, epar); //retrieve best fit parameters
		bestFitParams.push_back(par);
		e_bestFitParams.push_back(epar);
	}

	//Store the best-fit spline in a TGraph
	gsl_spline_init (spline_GLOB, &xData[0], &bestFitParams[0], xData.size()); //initialize the spline

	int nPointsSpline = int ((xData[nPoints-1]-xData[0])/stepSpline); //calculate the number of points of the spline

	vector< vector<double> > cSplineValues;
	cSplineValues.push_back( vector< double > ());//x
	cSplineValues.push_back( vector< double > ());//y
	for (int i= 0; i < nPointsSpline; i++){
		cSplineValues[0].push_back(xData[0]+i*stepSpline);
		cSplineValues[1].push_back(gsl_spline_eval (spline_GLOB, cSplineValues[0].back(), acc_GLOB));
	}

	return cSplineValues;
}

std::vector<double> run_fit(double parSRS){
    TMinuit minuit;
    double arglist[10];
    int iflag;
    arglist[0] = -1;
    minuit.mnexcm("SET PRINT",arglist,1,iflag);
    minuit.mnexcm("SET NOW",arglist,0,iflag);
    minuit.SetFCN(fcn);
    minuit.mnexcm("SET STR",arglist,1,iflag);
    minuit.mnparm(0,"parSRS",parSRS,0.1,-100,100,iflag);

    arglist[0] = 10000;
    arglist[1] = 0.01;
    minuit.mnexcm("MIGRAD",arglist,2,iflag);
    minuit.mnexcm("HESSE",arglist,0,iflag);

    double fmin, fedm,errdef;
    int npari, nparx, istat;
    minuit.mnstat(fmin,fedm,errdef,npari,nparx,istat);

    double val,err;
    std::vector<double> ret_ary;
    for(int p = 0; p < 1; p++) {
        minuit.GetParameter(p, val, err);
        ret_ary.push_back(val);
    }
    ret_ary.push_back(fmin);

    return ret_ary;
}

vector<double> fit()
{
  vector<double> out(2*nbeta+1);
  
  TMinuit minu;
  //minu.SetPrintLevel(-1);
  minu.SetFCN(ch2);
  //set tolerance
  double tol=1e-16;
  int iflag;
  //minu.SetMaxIterations(10000000);
  //minu.mnexcm("SET ERR",&tol,1,iflag);
  
  for(int ib=0;ib<nbeta;ib++)
    {
      minu.DefineParameter(2*ib+0,"infi",1,0.001,0,0);
      minu.DefineParameter(2*ib+1,"coef",1,0.001,0,0);
    }
  minu.DefineParameter(2*nbeta,"M",0,0.001,0,0);
  
  minu.Migrad();
  
  double dum;
  for(int ib=0;ib<=2*nbeta;ib++) minu.GetParameter(ib,out[ib],dum);
  
  return out;
}

void fit(boot &A,boot &B,boot &C,bvec &X,bvec &Y)
{
  //copy X
  X_fit=new double[nens];
  for(int iens=0;iens<nens;iens++) X_fit[iens]=X[iens].med();
  Y_fit=new double[nens];
  err_Y_fit=new double[nens];
  
  TMinuit minu;
  minu.SetPrintLevel(-1);
  
  minu.DefineParameter(0,"A",0.0,0.0001,0,0);
  minu.DefineParameter(1,"B",0.0,0.0001,0,0);
  minu.DefineParameter(2,"C",0.0,0.0001,0,0);
  
  minu.SetFCN(chi2_wr);
  
  double C2;
  for(int iboot=0;iboot<nboot+1;iboot++)
    {
      if(iboot>0)
        minu.SetPrintLevel(-1);
      
      minu.DefineParameter(3,"a380",lat[0][iboot],0.0001,0,0);
      minu.DefineParameter(4,"a390",lat[1][iboot],0.0001,0,0);
      minu.DefineParameter(5,"a405",lat[2][iboot],0.0001,0,0);
      minu.DefineParameter(6,"a420",lat[3][iboot],0.0001,0,0);
      minu.FixParameter(3);
      minu.FixParameter(4);
      minu.FixParameter(5);
      minu.FixParameter(6);
      
      for(int iens=0;iens<nens;iens++)
        {
          Y_fit[iens]=Y.data[iens].data[iboot];
          err_Y_fit[iens]=Y.data[iens].err();
        }
      
      //minimize
      minu.Migrad();
      
      //get back parameters
      double dum;
      minu.GetParameter(0,A.data[iboot],dum);
      minu.GetParameter(1,B.data[iboot],dum);
      minu.GetParameter(2,C.data[iboot],dum);
      
      double lat_med[4]={lat[0].med(),lat[1].med(),lat[2].med(),lat[3].med()};
      if(iboot==0) C2=chi2(A.data[iboot],B[iboot],C[iboot],lat_med);
    }
  
  //calculate the chi2
  cout<<"A = ("<<A<<"), B=("<<B<<"), C=("<<C<<")"<<endl;
  cout<<"Chi2 = "<<C2<<" / "<<nens-3<<" = "<<C2/(nens-3)<<endl;
  
  delete[] X_fit;
  delete[] Y_fit;
  delete[] err_Y_fit;
}

//______________________________________________________________________________
void Ifit()
{
// The z values
	z[0]=1;
	z[1]=0.96;
	z[2]=0.89;
	z[3]=0.85;
	z[4]=0.78;
// The errors on z values
        Float_t error = 0.01;
	errorz[0]=error;
	errorz[1]=error;
	errorz[2]=error;
	errorz[3]=error;
	errorz[4]=error;
// the x values
	x[0]=1.5751;
	x[1]=1.5825;
	x[2]=1.6069;
	x[3]=1.6339;
	x[4]=1.6706;
// the y values
	y[0]=1.0642;
	y[1]=0.97685;
	y[2]=1.13168;
	y[3]=1.128654;
	y[4]=1.44016;

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

   Double_t arglist[10];
   Int_t ierflg = 0;

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

// Set starting values and step sizes for parameters
   static Double_t vstart[4] = {3, 1 , 0.1 , 0.01};
   static Double_t step[4] = {0.1 , 0.1 , 0.01 , 0.001};
   gMinuit->mnparm(0, "a1", vstart[0], step[0], 0,0,ierflg);
   gMinuit->mnparm(1, "a2", vstart[1], step[1], 0,0,ierflg);
   gMinuit->mnparm(2, "a3", vstart[2], step[2], 0,0,ierflg);
   gMinuit->mnparm(3, "a4", vstart[3], step[3], 0,0,ierflg);

// Now ready for minimization step
   arglist[0] = 500;
   arglist[1] = 1.;
   gMinuit->mnexcm("MIGRAD", 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);

}

// --------------------------------- //
void make_lin_fit(double & slope, double & d_slope, double & offset, double & d_offset){
    TMinuit min;
    min.SetPrintLevel(-1);
    //min.SetPrintLevel(0);
    int err = min.DefineParameter(0, "slope", slope, d_slope, 0.05, 1.0);
    assert(err==0);
    err = min.DefineParameter(1, "offset", offset, d_offset, 0.001, 0.2);
    assert(err==0);
    min.SetFCN(chi2_linear);
    min.mnmigr();
    min.GetParameter(0, slope, d_slope);
    min.GetParameter(1, offset, d_offset);
}

std::vector<double> run_fit(double parOOR, double parOEFG, double parOFTMR,
                            double parOTO, double parSRS) {
    TMinuit minuit;
    double arglist[10];
    int iflag;
    arglist[0] = -1;
    minuit.mnexcm("SET PRINT",arglist,1,iflag);
    minuit.mnexcm("SET NOW",arglist,0,iflag);
    minuit.SetFCN(fcn);
    arglist[0] = 2;
    minuit.mnexcm("SET STR",arglist,1,iflag);
    minuit.mnparm(0,"parOOR",parOOR,0.1,0,0,iflag);//or
    minuit.mnparm(1,"parEFG",parOEFG,0.1,0,0,iflag);//efgp
    minuit.mnparm(2,"parFTMR",parOFTMR,0.1,0,0,iflag);//ftr
    minuit.mnparm(3,"parTO",parOTO,0.1,0,0,iflag);//top
    minuit.mnparm(4,"parSRS",parSRS,0.1,0,0,iflag);//srs

    if (parSRS == 0)
        minuit.FixParameter(4);

    arglist[0] = 10000;
    arglist[1] = 0.0001;
    minuit.mnexcm("MIGRAD",arglist,2,iflag);
    minuit.mnexcm("HESSE",arglist,0,iflag);

    double fmin, fedm,errdef;
    int npari, nparx, istat;
    minuit.mnstat(fmin,fedm,errdef,npari,nparx,istat);

    double val,err;
    std::vector<double> ret_ary;
    for(int p = 0; p < 5; p++) {
        minuit.GetParameter(p, val, err);
        ret_ary.push_back(val);
    }
    ret_ary.push_back(fmin);

    return ret_ary;
}

Double_t Fitter::MinuitFit() {
  TMinuit *gMinuit = new TMinuit(1);
  gMinuit->SetFCN(ChiSquare);
  
  //settings 
  Double_t arglist[6];
  Int_t ierflg = 0;
  
  // Set starting values and step sizes for parameters
  gMinuit->mnparm(0, "Sig", _W1, _step, _W1min, _W1max,ierflg);
  //gMinuit->mnparm(0, "Sig", _W1, _step,0,0,ierflg);
  
  // Now ready for minimization step
  gMinuit->mnexcm("MIGRAD", arglist ,0,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);
  
  //gMinuit->mnexcm("MINOS", arglist ,0,ierflg);
  for(int i=0; i<1; i++)
    {
      gMinuit -> GetParameter(i, _W1, _W1err);
      cout<<"Parameter fit  "<<i<<" value "<<_W1<<" err "<<_W1err<<endl;
    }
  
  delete gMinuit;
  
  //calculate bkg error:
  double entries_mc = _hsig->Integral(_hdata->GetXaxis()->FindBin(_xmin), _hdata->GetXaxis()->FindBin(_xmax));
  double entries_bkg = _hbkg->Integral(_hdata->GetXaxis()->FindBin(_xmin), _hdata->GetXaxis()->FindBin(_xmax));
  _W2err = _W1err*entries_mc/entries_bkg;

  return amin;
}

void parab_fit(double ch2,jack &a,jack &b,jack &c,double *X,jvec Y,const char *path,bool fl_parab=1)
{
  ofstream out(path);
  out<<"@type xydy\n";
  
  double MX=0;
  //copy X
  nc_fit=Y.nel;
  X_fit=new double[nc_fit];
  for(int iel=0;iel<nc_fit;iel++)
    if(use[iel])
    {
      X_fit[iel]=X[iel];
      out<<X[iel]<<" "<<Y[iel]<<endl;
      if(X[iel]>MX) MX=X[iel];
    }
  Y_fit=new double[nc_fit];
  err_Y_fit=new double[nc_fit];
  
  TMinuit minu;
  minu.SetPrintLevel(0);
  
  int npars=3;
  minu.DefineParameter(0,"A",0.0,0.0001,0,0);
  minu.DefineParameter(1,"B",0.0,0.0001,0,0);
  minu.DefineParameter(2,"C",0.0,0.0001,0,0);
  if(!fl_parab)
    {
      minu.FixParameter(0);
      npars--;
    }
  minu.SetFCN(chi2_wr);
  
  for(int ijack=0;ijack<=njack;ijack++)
    {
      minu.DefineParameter(2,"C",fixed_A[ijack],0.0001,0,0);
      minu.FixParameter(2);
      if(ijack>0)
        minu.SetPrintLevel(-1);
      
      for(int iel=0;iel<nc_fit;iel++)
        {
          Y_fit[iel]=Y.data[iel].data[ijack];
          err_Y_fit[iel]=Y[iel].err();
        }
      
      //minimize
      minu.Migrad();
      
      //get back parameters
      double dum;
      minu.GetParameter(0,a.data[ijack],dum);
      minu.GetParameter(1,b.data[ijack],dum);
      minu.GetParameter(2,c.data[ijack],dum);
    }
  
  out<<"&\[email protected] xy\n";
  for(double t=0;t<MX;t+=MX/100)
    out<<t<<" "<<parab(a[njack],b[njack],c[njack],t)<<endl;
  
  out.close();
}

void runAngleOpt() {

  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");

  double startVal=0;               
  double stepSize=0.1;       
  double minVal=-2;           
  double maxVal=2;            
  Double_t p0 = 0;

  //Load libraries. Need to have ANITA_UTIL_INSTALL_DIR/lib and ROOTSYS/lib in the LD_LIBRARY_PATH               
  gSystem->Load("libfftw3.so");
  gSystem->Load("libMathMore.so");
  gSystem->Load("libPhysics.so");  
  gSystem->Load("libGeom.so");  
  gSystem->Load("libMinuit.so");  
  gSystem->Load("libRootFftwWrapper.so");         
  gSystem->Load("libAnitaEvent.so");      
  gSystem->Load("libAnitaCorrelator.so");

  AnitaGeomTool *fGeomTool = AnitaGeomTool::Instance();
  gSystem->CompileMacro("anglePlotterOpt.C","k");
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(1);
  myMin->SetObjectFit(anglePlotterOpt);
  myMin->SetFCN(iHateRoot);
  //setArray();

  for(int u = 0; u < 1; u++){

    int middle = 16;

      for(int y = 16; y <32; y++){

	int leftOpt, rightOpt;
	fGeomTool->getThetaPartners(middle,leftOpt,rightOpt); 


	myMin->DefineParameter(0, "antNum", middle, stepSize, minVal, maxVal);
	myMin->FixParameter(0);  
    
	myMin->DefineParameter(1, "deltaT", startVal, stepSize, minVal, maxVal);

  
	Double_t deltaT,deltaTErr;
  
	//*********MINUIT METHOD*******************
	myMin->SetPrintLevel(-1);
	myMin->Migrad();   
	myMin->GetParameter(1,deltaT,deltaTErr);

  

	setValue(rightOpt,deltaT);

	//	printArray();

	//	cout << middle << "  " << rightOpt << "  " << deltaT << endl;
  
	cout << "deltaTArrayMod[" << rightOpt << "] = " << deltaT << ";" << endl;

	middle = rightOpt;

      }

  }

  
  //   myMin->DeleteArrays();

  //   myMin->DeleteArrays();

  
}

//_________________________________________________________________________________
void multipleSplinesWithHistogramsVb(int iEventLook = 163, int nEvents = 1000, int nPoints = 9, double seed = 231) 
{
	double lowerBound = 10; //bounds for random vector function
	double upperBound = 20;
	double lowerErrorBound = 1;
	double upperErrorBound = 2;

	//Load the data
	vector< vector<double> > xEvents, yEvents, yErrorEvents; //each of these is a vector of vectors (of randomized data points)
	for(int i = 0; i < nEvents; i++) 
	{
		vector <double> xData, yData, yErrorData; //temporary vectors that are only used to get random values from FillRand function
		FillRandVectors(nPoints, xData, yData, yErrorData, seed*(i+1), lowerBound, upperBound, lowerErrorBound, upperErrorBound); //populates random vectors for y values and y error vector
		xEvents.push_back(xData);
		yEvents.push_back(yData);
		yErrorEvents.push_back(yErrorData);
	}

	//Intialization of the variables
	const int npar = nPoints;
	const int orderSpline = 4;
	const int nbreak = npar+2-orderSpline;
	double stepSpline = 0.01;
	double xminBSplineWorkspace = 0;
	double xmaxBSplineWorkspace = 9;
	double startCSplineWorkspace = 15.;
	double stepCSplineWorkspace = 1.5;

	acc_GLOB = gsl_interp_accel_alloc ();
	spline_GLOB = gsl_spline_alloc (gsl_interp_cspline, nPoints);	
	bw_GLOB = gsl_bspline_alloc(orderSpline, nbreak);
	
	//B- and C-splines
	clock_t tbstart, tbstop, tcstart, tcstop;
	vector <double> timeb, timec;

	vector< vector<double> > xBSplineValues, yBSplineValues;
	vector< vector<double> > xCSplineValues, yCSplineValues; 

//Setup for the C-spline_________________________________________________________________________
	TMinuit *myMinuit = new TMinuit(npar);  //initialize TMinuit with a maximum of npar 
	myMinuit->SetFCN(fcn);
	myMinuit->SetPrintLevel(-1);//No output: -1, output:1

	double arglist[10];
	int ierflg = 0;
	arglist[0] = 1;
	myMinuit->mnexcm("SET ERR", arglist, 1, ierflg);

	//Initialize Minuit
	vector<double> vstart, vstep;
	for(int i=0; i<npar; i++) 	//set starting values and step sizes for parameters
	{
		vstart.push_back(startCSplineWorkspace);
		vstep.push_back(stepCSplineWorkspace);
	}

	for (int i = 0; i < npar; i++) {
		stringstream ss;
		ss<<"a"<<i;
		myMinuit->mnparm(i, ss.str().c_str(), vstart.at(i), vstep.at(i), 0, 0, ierflg);
	}

//Setup for the B-spline_________________________________________________________________________

	//Looping begins for the calculations of the B and C-splines for each event
	for(int i = 0; i < (int)xEvents.size(); i++)
	{
		//Populate the global variables
		xData_GLOB = xEvents.at(i);
		yData_GLOB = yEvents.at(i);
		yErrorData_GLOB = yErrorEvents.at(i);

		tbstart = clock();
		vector< vector<double> > bSplineValues = bSpline(nPoints, npar, xEvents.at(i), yEvents.at(i), yErrorEvents.at(i), stepSpline, xminBSplineWorkspace, xmaxBSplineWorkspace);
		tbstop = clock();
		timeb.push_back(((float)tbstop-(float)tbstart)/ (CLOCKS_PER_SEC/1000.) );		

		std::cout<<timeb.back()<<std::endl;

		xBSplineValues.push_back(bSplineValues.at(0));
		yBSplineValues.push_back(bSplineValues.at(1));

		tcstart = clock();
		vector< vector<double> > cSplineValues = cSpline(nPoints, npar, xEvents.at(i), stepSpline, myMinuit);
		tcstop = clock();
		timec.push_back(((float)tcstop-(float)tcstart)/ (CLOCKS_PER_SEC/1000.) );		

		xCSplineValues.push_back(cSplineValues.at(0));
		yCSplineValues.push_back(cSplineValues.at(1));
	}

	//Histograms______________________________________________________________________________________

	//Time
	int nbins = 100;
	double xlow = 0;
	double xup = 1.;

//.........这里部分代码省略.........

/*DoMunuit{{{*/
void DoMinuit(double *par)
{
	cout << "      *************************************************" << endl; 
	cout << "      *          Minimization Fitting for             *" << endl;
	cout << "      *              F(y) Function                    *" << endl;
	cout << "      *             Z. Ye 09/20/2010                  *" << endl;
	cout << "      *************************************************" << endl; 
	cout << endl;
	gSystem->Load("libMinuit");
	static const Int_t iParaNum=5;

	TMinuit *pMinuit = new TMinuit(iParaNum);  //initialize TMinuit with a maximum of iParaNum params
	pMinuit->SetFCN(myFcn);

	Double_t arglist[10];
	Int_t ierflg = 0;

	/// set print level
	/*  SET PRIntout  <level>
		Sets the print level, determining how much output will be
		produced. Allowed values and their meanings are displayed
		after a SHOw PRInt command, and are currently <level>=:
		[-1]  no output except from SHOW commands
		[0]  minimum output
		[1]  default value, normal output
		[2]  additional output giving intermediate results.
		[3]  maximum output, showing progress of minimizations. */
	arglist[0] = 1;
	pMinuit->mnexcm("SET PRIntout",arglist,1,ierflg);

	/*
	   SET NOWarnings
	   Supresses Minuit warning messages.
	   SET WARnings
	   Instructs Minuit to output warning messages when suspicious
	   conditions arise which may indicate unreliable results.
	   This is the default.
	   */
	arglist[0] = 1;
	pMinuit->mnexcm("SET NOWarnings",arglist,0,ierflg);

	//Set Name of Parameters
	TString namePar[iParaNum]={"f0","bigB","alpha","a","b"};

	// Double_t inPar[iParaNum]={Target->f0,   
	// 			    Target->B,   
	// 			    Target->alpha,
	// 			    Target->a,
	// 			    Target->b};	
	Double_t inPar[iParaNum]={par[0],par[1],par[2],par[3],par[4]};   

	//Set Stepsize of Parameters
	Double_t step[iParaNum]={ 0.00000001,0.00000001,0.00001,0.00000001,0.00000001};

	//Set Min of Parameters, value==0 for No-Bound
	// Double_t minVal[iParaNum]={ 0.0001,0.0001,100.0,0.0001,0.0001};
	Double_t minVal[iParaNum]={ 0.0,0.0,0.0,0.0,0.0};

	//Set Max of Parameters, value==0 for No-Bound
	// Double_t maxVal[iParaNum]={ 0.01,0.01,300.0,0.01,0.01};
	Double_t maxVal[iParaNum]={ 0.0,0.0,0.0,0.0,0.0};

	//Initializing Parameters
	for(int ii=0;ii<iParaNum;ii++)
	{
		pMinuit->DefineParameter(ii,namePar[ii], 
				inPar[ii], step[ii], minVal[ii], maxVal[ii]);
	}

	//Fix parameters
	//arglist[0] = 3;   pMinuit->mnexcm("FIX ", arglist ,1,ierflg);

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

	/*
	//SET LIMits  0->[parno]  1->[lolim]  2->[uplim]
	//	  arglist[0] = 0;	   //this means set all parameters with the same limit
	//	  arglist[0] = 1;	   //this means set 1st parameters with the specified limit
	//	  arglist[0] = 2;	   //this means set 2nd parameters with the sepecified limit
	arglist[0] = 0;
	arglist[1] = 0.;
	arglist[2] = 0.;
	pMinuit->mnexcm("SET LIMits",arglist,3,ierflg);
	*/

	//Set Iteration Times and Use Migrad to start minimization
	arglist[0] = 500;
	arglist[1] = 0.01;
	pMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

	//Get Results

	for(int ii = 0; ii < iParaNum; ii++)
	{
		pMinuit->GetParameter(ii,outPar[ii],err[ii]);    
		par[ii] = outPar[ii];
	}
	//Put the results into a file
//.........这里部分代码省略.........

//_________________________________________________________________________________
void integratedSplinesV4a(double seed = 231)
{
    //Load the data
    int nEvents = 1000; //number of times the data will be randomized
    int nPoints = 9;

    vector< vector<double> > xEvents, yEvents, yErrorEvents; //each of these is a vector of vectors (of randomized data points)
    for(int i = 0; i < nEvents; i++)
    {
        vector <double> xData, yData, yErrorData; //temporary vectors that are only used to get random values from FillRand function
        FillRandVectors(nPoints, xData, yData, yErrorData, seed*(i+1)); //populates random vectors for y values and y error vector
        xEvents.push_back(xData);
        yEvents.push_back(yData);
        yErrorEvents.push_back(yErrorData);
        //Populate the global variables
        xData_GLOB = xData;
        yData_GLOB = yData;
        yErrorData_GLOB = yErrorData;
    }

    //Intialization of the variables
    const int npar = nPoints;
    const int orderSpline = 4;
    const int nbreak = npar+2-orderSpline;
    double stepSpline = 0.01;
    double xminBSplineWorkspace = 0;
    double xmaxBSplineWorkspace = 9;
    double startCSplineWorkspace = 15.;
    double stepCSplineWorkspace = 1.5;

    acc_GLOB = gsl_interp_accel_alloc ();
    spline_GLOB = gsl_spline_alloc (gsl_interp_cspline, nPoints);
    gsl_bspline_workspace *bw = gsl_bspline_alloc(orderSpline, nbreak);

    //Setting up TMinuit for C-spline minimization
    TMinuit *myMinuit = new TMinuit(npar);  //initialize TMinuit with a maximum of npar (5)
    myMinuit->SetFCN(fcn);
    myMinuit->SetPrintLevel(-1);//No output: -1, output:1

    double arglist[10];
    int ierflg = 0;
    arglist[0] = 1;
    myMinuit->mnexcm("SET ERR", arglist, 1, ierflg);

    vector<double> vstart, vstep;
    for(int i=0; i < npar; i++) 	//set starting values and step sizes for parameters
    {
        vstart.push_back(startCSplineWorkspace);
        vstep.push_back(stepCSplineWorkspace);
    }

    for (int i = 0; i < npar; i++)
    {
        stringstream ss;
        ss<<"a"<<i;
        myMinuit->mnparm(i, ss.str().c_str(), vstart.at(i), vstep.at(i), 0, 0, ierflg);
    }

    //Perform the Minuit fit
    arglist[0] = 500;
    arglist[1] = 1.;

    //B and C spline loops__________________________________________________________________________
    clock_t tbstart, tbstop, tcstart, tcstop;
    vector <double> timeb, timec;

    vector< vector<double> > xBSplineValues, yBSplineValues;
    vector< vector<double> > xCSplineValues, yCSplineValues;
    for(int i = 0; i < (int)xEvents.size(); i++)
    {
        tbstart = clock();
        vector< vector<double> > bSplineValues = bSpline(nPoints, npar, xEvents.at(i), yEvents.at(i), yErrorEvents.at(i), stepSpline, xminBSplineWorkspace, xmaxBSplineWorkspace, bw);
        tbstop = clock();
        timeb.push_back(((float)tbstop-(float)tbstart)/ (CLOCKS_PER_SEC/1000.) );

        xBSplineValues.push_back(bSplineValues.at(0));
        yBSplineValues.push_back(bSplineValues.at(1));

        tcstart = clock();
        vector< vector<double> > cSplineValues = cSpline(nPoints, npar, xEvents.at(i), stepSpline, myMinuit, arglist, ierflg);
        tcstop = clock();
        timec.push_back(((float)tcstop-(float)tcstart)/ (CLOCKS_PER_SEC/1000.) );

        xCSplineValues.push_back(cSplineValues.at(0));
        yCSplineValues.push_back(cSplineValues.at(1));
    }

    //Histograms______________________________________________________________________________________

    //Time
    int nbins = 100;
    double xlow = 0;
    double xup = 5.;

    TH1D *hTimeB = new TH1D("Time","; time [ms]; number of runs", nbins, xlow, xup);
    hTimeB->SetStats(0);
    TH1D *hTimeC = new TH1D("TimeC","; time [ms]; number of runs", nbins, xlow, xup);
    hTimeC->SetLineColor(kRed);
    hTimeC->SetStats(0);
//.........这里部分代码省略.........

void runTopRingOpt() {

  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");

  double startVal=0;               
  double stepSize=0.1;       
  double minVal=-0.5;           
  double maxVal=0.5;            
  Double_t p0 = 0;

  //Load libraries. Need to have ANITA_UTIL_INSTALL_DIR/lib and ROOTSYS/lib in the LD_LIBRARY_PATH               
  gSystem->Load("libfftw3.so");
  gSystem->Load("libMathMore.so");
  gSystem->Load("libPhysics.so");  
  gSystem->Load("libGeom.so");  
  gSystem->Load("libMinuit.so");  
  gSystem->Load("libRootFftwWrapper.so");         
  gSystem->Load("libAnitaEvent.so");      
  gSystem->Load("libAnitaCorrelator.so");

  AnitaGeomTool *fGeomTool = AnitaGeomTool::Instance();
  gSystem->CompileMacro("topRingOpt.C","k");
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(150);
  myMin->SetObjectFit(topRingOpt);
  myMin->SetFCN(iHateRoot);
  //setArray();

  double startValDeltaT[16] ={0};
  double startValR[16] ={0};
  double startValPhi[16] ={0};

 //    startValDeltaT[0] =  -0.0519515;     startValR[0] =   -0.0101463;         startValPhi[0] =      -0.00473836 ;   
//     startValDeltaT[1] =     -0.0597062;   startValR[1] =      -0.02577;       startValPhi[1] =   0.00864501 ; 
//     startValDeltaT[2] =     -0.081435;     startValR[2] =     -0.000224044;    startValPhi[2] =   -0.000630649;
//     startValDeltaT[3] =     0.0118873;     startValR[3] =      0.019945;       startValPhi[3] =    0.014016;   
//     startValDeltaT[4] =      0.017917;      startValR[4] =     -0.00297559;    startValPhi[4] =    0.0224936 ;
//     startValDeltaT[5] =      0.0377119;     startValR[5] =     -0.014872;       startValPhi[5] =   0.0163349; 
//     startValDeltaT[6] =     -0.0426158;     startValR[6] =    -0.0562555;      startValPhi[6] =    0.0220065 ; 
//     startValDeltaT[7] =     -0.0221673;     startValR[7] =    -0.034104 ;      startValPhi[7] =    0.0158545 ; 
//     startValDeltaT[8] =      0.0263739;     startValR[8] =     0.00248804;      startValPhi[8] =    0.013246 ; 
//     startValDeltaT[9] =      -0.0938419;    startValR[9] =     -0.00344703;     startValPhi[9] =    -0.00718616; 
//     startValDeltaT[10] =      0.145264;      startValR[10] =     -0.0121874 ;     startValPhi[10] =     0.0156988 ;  
//     startValDeltaT[11] =      0.118105;      startValR[11] =     -0.0337033 ;      startValPhi[11] =   -0.00324182 ;
//     startValDeltaT[12] =      0.321805;      startValR[12] =     0.0134362  ;      startValPhi[12] =   -0.00190277 ;
//     startValDeltaT[13] =      0.0197693;     startValR[13] =     -0.000656063;      startValPhi[13] =   -0.0162318  ;
//     startValDeltaT[14] =      -0.115263;     startValR[14] =     0.0495637   ;      startValPhi[14] =   -0.0198119 ; 
//     startValDeltaT[15] =      -0.255707;   startValR[15] =     0.00189892  ;     startValPhi[15] =     0.0383932  ;
    


 for(int y = 0; y <16; y++){

   ofstream newfile("newSimonNumbers.txt");

   char name[30];
   sprintf(name,"r%d",y);
	myMin->DefineParameter(y, name, startValR[y], stepSize, minVal, maxVal);
   sprintf(name,"z%d",y);
	myMin->DefineParameter(y+16, name, startValDeltaT[y], stepSize, minVal, maxVal);
	   sprintf(name,"phi%d",y);
	myMin->DefineParameter(y+32, name, startValPhi[y], stepSize, minVal, maxVal);
 }


  
	Double_t deltaR[32],deltaRErr[32];
	Double_t deltaZ[32],deltaZErr[32];
	Double_t deltaPhi[32],deltaPhiErr[32];
  
	//*********MINUIT METHOD*******************
	myMin->SetPrintLevel(-1);
	myMin->Migrad();   

	for(int u = 0; u <16; u++){
	  myMin->GetParameter(u,deltaR[u],deltaRErr[u]);
	  //cout << "deltaR[" << u << "] = " << deltaR[u] ;
	
	
	  myMin->GetParameter(u+16,deltaZ[u],deltaZErr[u]);
	  //cout << " deltaZ[" << u << "] = " << deltaZ[u] ;
	
	  myMin->GetParameter(u+32,deltaPhi[u],deltaPhiErr[u]);
	  //cout << " deltaPhi[" << u << "] = " << deltaPhi[u] << ";" << endl;

	  newfile << u << "  " << deltaZ[u]<< "  " << deltaR[u]<< "  " << deltaPhi[u]<< "  " << 0 << endl;


	}


  
}

//.........这里部分代码省略.........
              cout << endl;
            }

	    gEntries++;
	    bNS = false;
	    bEW = false;
	  }
	  break;
	}
      }
    }
   }

   cout << gEntries << endl;
   if(NoChannels30+NoChannels17 > gMaxCh || gEntries > gMaxEvents) {
     cout << "constants limit reached" << endl;
     exit(1);
   }
 
   TMinuit *gMinuit = new TMinuit(10);
   gMinuit->SetPrintLevel(-1);
   Double_t arglist[10];
   Int_t ierflg = 0;
   int NoPar = 6;
   
   bool AntFit = false;
   
   if(AntFit){
   // Antenna Fit Settings
   //++++++++++++++++++++++
   cout << "Antenna Fit" << endl;
   cout << "+++++++++++" << endl;
   
   gMinuit->SetFCN(FcnParFS);
   
// no effect
  arglist[0] = 1;
  gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);

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

   gMinuit->mnparm(0, "scale",  1., 0.01, 0, 0, ierflg);
   gMinuit->mnparm(1, "r0",    400., 1, 0, 0, ierflg);
   gMinuit->mnparm(2, "Epower",  1., 0.01, 0.8, 1.2, ierflg);
//   gMinuit->mnparm(2, "Epower",  1., 0.1, 0, 0, ierflg);
   gMinuit->mnparm(3, "Zphase",  0., 0.1, -TMath::Pi()*2, TMath::Pi()*2, ierflg);
   gMinuit->mnparm(4, "Bphase",  0., 0.1, -TMath::Pi()*2, TMath::Pi()*2, ierflg);
   gMinuit->mnparm(5, "offset",  2., 0.1, 0, 0, ierflg);
//   gMinuit->FixParameter(0);
//   gMinuit->FixParameter(2);

   gMinuit->FixParameter(3);
   gMinuit->FixParameter(4);
   gMinuit->FixParameter(5);
   }
   else{
   // Channel Fit Settings
   //+++++++++++++++++++++
   cout << "Channel Fit" << endl;
   cout << "+++++++++++" << endl;

   gMinuit->SetFCN(FcnParFSChannel);
   
// no effect