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

//______________________________________________________________________________
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);

}

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< 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;
}

void unbinFitosc_pts()
{
  TMinuit *gMinuit = new TMinuit(4);
  gMinuit -> SetFCN(mll_fit_pts);
  
  Double_t arglist[10];
  Double_t vstart[10];
  Double_t step[10];

  Double_t p0=0,p1=1,p2=2,p3=3;
  Int_t ierflg=0;
 
  Double_t par0;
  Double_t err0;

  // Set to 0.5 for likelihood, 1 for chisquare

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

  // Set initial values of parameters

  vstart[0] = oscpar0_init;

  step[0] = oscerr0_step;

  gMinuit->mnparm(0, "mistag", vstart[0], step[0], 0., 0., ierflg);


  // Do minimization
  arglist[0] = 3000.;
  arglist[1] = 0.1;
  gMinuit->mnexcm("CALL FCN", &p1, 1, ierflg);
  gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
  gMinuit->mnexcm("CALL FCN", &p3, 1, ierflg);

  // Get parameters 
  gMinuit->GetParameter(0,par0,err0);


  // Export the fitted parameters so we can graph them.

  fitpar=par0;

  fiterr=err0;


  RooComplex ct = RooMath::ComplexErrFunc(2.,1.);
  Double_t ctre=ct.re();
  Double_t ctim=ct.im();
  printf("ct  %10.5f  %10.5f  \n",ctre,ctim);

  delete gMinuit;

}

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;
}

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

 cout << " the content of data[2] is = " << Pcdata[2] << "+-"<<Pedata[2] << endl ;
 cout << " the content of pythia[2] is = " << Pcpythia[2] << "+-"<<Pepythia[2] << endl ;
 cout << " the content of pythia[2] is = " << Pcpompyt[2] << "+-"<<Pepompyt[2] << endl ;


 // Here  we go: the minuit show
 
  TMinuit *gMinuit = new TMinuit(NPAR);  //initialize TMinuit with a maximum of NPAR params
  gMinuit->SetFCN(fcn);
  
  Double_t arglist[NPAR]; // ???
  Int_t ierflg = 0;

// Set starting values and step sizes for parameters

  Double_t vstart[NPAR] ;
  Double_t step[NPAR] ;
  Double_t par[NPAR] ,fpar[NPAR];
  char parName[NPAR];
  Int_t n;
  parName[0] = "Pythia";
  parName[1] = "Pompyt";

  par[0] = Hdata->Integral()/Hpythia->Integral();
  par[1] = 0.05*Hdata->Integral()/Hpompyt->Integral();


  for (n=0 ; n<NPAR ; n++)
    {
      vstart[n] = par[n]  ;
      step[n] = 0.5 ;
      sprintf(parName,"a%d",n);
      gMinuit->mnparm(n, parName, vstart[n], step[n], 0,0,ierflg);
    }

  cout << "par[0] set to = " << par[0] << " while par[1] = " << par[1] << endl;
  arglist[0] = 1;

  gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);  
  gMinuit->mnexcm("SET PRINT", arglist,1,ierflg);

  //Scan on parameter = 1
  arglist[0] = 1;
  gMinuit->mnexcm("SCAN",  arglist,1,ierflg);
  
  //Maximum number of calls
  arglist[0] = 500;
  gMinuit->mnexcm("MIGRAD", arglist,1,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);


  //Get the final fit parameters
  for (n=0 ; n<NPAR ; n++)
    {
      Double_t parameter, erro ;
      gMinuit->GetParameter(n,parameter,erro) ;
      fpar[n] = parameter ;

//_________________________________________________________________________________
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.;

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

//_________________________________________________________________________________
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);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    double fitpar[10];
    double fitparErr[10];
    
  double smearcat[4] = {1.1,1.1,1.1,1.1};
  double smearcatMax[4] = {3,3,3,3};

  double deltaEcat[4] = {0,0,0,0};
  
  
  if(smearMethod=="uncorrSmear"){
    smearcatMax[0] = 0.05;
    smearcatMax[1] = 0.05;
    smearcatMax[2] = 0.05;
    smearcatMax[3] = 0.05;

    smearcat[0] = 0.007;
    smearcat[1] = 0.01;
    smearcat[2] = 0.015;
    smearcat[3] = 0.02;
    
    if(fitdet==2){
      smearcat[0] = 0.02;
      smearcat[1] = 0.02;
      smearcat[2] = 0.02;
      smearcat[3] = 0.02;
    }
    
  }
  

  for(int j=0; j<4; j++){
    TString parname = TString (Form("smearcat%d",j));
    if(smearMethod=="corrSmear"){
      minuit->mnparm(j, parname, smearcat[j], STEPMN, 0.5,smearcatMax[j],ierflg);
    }else  if(smearMethod=="uncorrSmear"){
      minuit->mnparm(j, parname, smearcat[j], STEPMN, 0,smearcatMax[j],ierflg);
    }

    fitpar[j] = smearcat[j]; //initialized values
  }
  
  for(int j=4; j<npar; j++){
    TString parname = TString (Form("scalecat%d",j-4));
    minuit->mnparm(j, parname, deltaEcat[j-4], STEPMN, -0.03,0.03,ierflg);
    fitpar[j] = deltaEcat[j-4]; //initialized values
  }
  

  if(fixscale==1){
    for(int j=4; j<npar; j++){
      minuit->FixParameter(j);
    }
  }
  


  if(testpaircat==1){
    minuit->FixParameter(0);
    minuit->FixParameter(2);
    minuit->FixParameter(3);

    minuit->FixParameter(4);
    minuit->FixParameter(6);
    minuit->FixParameter(7);

  }else if( testpaircat==0){

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

//  for (int i=1;i<2000;i++){ //the Sig or Bkg histo has 2000 bins
   for (int i=1;i<= sigTemplate->GetNbinsX();i++){ 
    if(/* ( ( bkgTemplate2->GetBinContent(i) ) != 0) ||*/ ( ( bkgTemplate1->GetBinContent(i) ) != 0) || ( ( sigTemplate->GetBinContent(i) ) != 0 ) ){
    Number_of_points = Number_of_points+1;

    }
  }
//cout<<"Number_of_points:"<<Number_of_points<<endl;
int dof = Number_of_points-1 ;//dof=Number_of_points - free_par
                              //and in this case we normalize it so free_par =1
//cout<<"dof of this fit:"<<dof<<endl;



/////////


  TMinuit *gMinuit = new TMinuit(1); // initialize TMinuit with a maximum of 1 params
//  TMinuit *gMinuit = new TMinuit(2);  //initialize TMinuit with a maximum of 5 (1param??) params
  gMinuit->SetFCN(fcn); // sets function to minimize: fcn is Chi2 with errors on templates

  Double_t arglist[10];
  Int_t ierflg = 0; // status flag, it is 0 when ereything goes fine

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

  // set first parameter
  Double_t vstart = sigFrac_intial;
  //Double_t vstart = 0.48 ;//0.11;//intial value
  Double_t step = 0.001;
  gMinuit->mnparm(0, "fsig", vstart, step, 0,1,ierflg);

/*
  // set second parameter
    vstart = bkg1Frac_intial;
//    vstart = 0.4; //0.69;//intial value
    step = 0.001;
  gMinuit->mnparm(1, "fbk1", vstart, step, 0,1,ierflg);
*/

  // Now ready for minimization step
  arglist[0] = 1000;
  arglist[1] = 0.01;
  gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

  Double_t fsig=0;
  Double_t fsigerr=0;
  Double_t fbk1=0;
//  Double_t fbk1err=0;
  Double_t chi2 = 0;


  if ( ierflg == 0 ) 
    {

      // Print results
      Double_t amin,edm,errdef;
      Int_t nvpar,nparx,icstat;
      gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
      gMinuit->mnprin(3,amin);
      chi2 = (gMinuit->fAmin)/dof;
      gMinuit->GetParameter(0, fsig, fsigerr );  

//.........这里部分代码省略.........
    dataColl.push_back(hdata->GetBinContent(ibin));
    sigColl.push_back(hsig->GetBinContent(ibin));
    bkgColl.push_back(hbkg->GetBinContent(ibin));    
  }
  printf( " -----  Got %d, %d, %d events for fit ----- \n ", dataColl.size(), sigColl.size(), bkgColl.size() );  
  if ( dataColl.size() != sigColl.size() || sigColl.size()!=bkgColl.size() ) {
    printf(" error ...  inconsistent hit collection size \n");
    fin_data->Close();
    fin->Close();
    fin_gjet6000->Close();

    return fitted;
  }

  //--------------------------------------------------
  //init parameters for fit
  Double_t vstart[10] = {1., 1.};
  vstart[0] = sigfrac*ndata;
  vstart[1] = (1-sigfrac)*ndata;
 
  TMinuit *gMinuit = new TMinuit(NPAR);  
  gMinuit->Command("SET STR 1");
  gMinuit->SetFCN(fcn);
  Double_t arglist[10];
  Int_t ierflg = 0;
  
  arglist[0] = 1;
  gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
  arglist[0] = 1;
  gMinuit->mnexcm("SET PRINT", arglist ,1,ierflg);

  Double_t step[] = { 0.1, 0.1,};

  gMinuit->mnparm(0,  "Signal yield"  , vstart[0],  step[0], 0., ndata*2.  , ierflg);
  gMinuit->mnparm(1,  "background yield"  , vstart[1],  step[1], 0., ndata*2. , ierflg);
  
  printf(" --------------------------------------------------------- \n");
  printf(" Now ready for minimization step \n --------------------------------------------------------- \n");
  
  arglist[0] = 2000; // number of iteration
  arglist[1] = 1.;
  gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
  printf (" -------------------------------------------- \n");
  printf("Finished.  ierr = %2.2f \n", ierflg);

  info.clear();
  info_err.clear();

  double para[NPAR+1],errpara[NPAR+1];
  if ( ierflg == 0 ) 
    {
      for(int j=0; j<=NPAR-1;j++) {
        gMinuit->GetParameter(j, para[j],errpara[j]);
        para[NPAR] = dataColl.size();
        info.push_back(para[j]);
        info_err.push_back(errpara[j]);
        printf("Parameter (yeild) %d = %f +- %f\n",j,para[j],errpara[j]);
	
      }
      printf(" fitted yield %2.3f \n", (para[0]+para[1])/ndata );

      info.push_back(sigColl.size());

      //do minos if fit sucessed.
//       printf("         ---------------------------------------------------------\n");
//       printf("          Now call for minos step \n");

// -- main function
void fit_chi2(TH1D* dataInput, TH1D* sigTemplate, TH1D* bkgTemplate, std::string prefix, Double_t& sigFrac, Double_t& sigFrac_err)
{
  gStyle->SetOptStat(kFALSE);
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);

  Double_t scale=1.;

  data = (TH1D*)dataInput->Clone();
  data->SetName("data");
  data->SetLineColor(1);
  data->SetMarkerColor(1);
  data->SetXTitle("Fisher's isolation [GeV]");
  data->Sumw2();
  scale = 1.0/(Double_t)data->Integral(0,1000); 
  cout << "scale for data = " << scale << endl;
  data->Scale(scale);


  fit_result = (TH1D*)dataInput->Clone();
  fit_result->SetName("fit_result");
  fit_result->SetLineColor(8); 
  fit_result->SetMarkerColor(8);
  fit_result->SetLineStyle(2);
  fit_result->Sumw2();
  fit_result->Reset();


  signal_pos = (TH1D*)sigTemplate->Clone();
  signal_pos->SetName("signal_pos");
  signal_pos->SetLineColor(2);
  signal_pos->SetMarkerColor(2);
  signal_pos->SetFillColor(2);
  signal_pos->SetXTitle("Fisher's isolation [GeV]");
  signal_pos->Sumw2();
  scale = 1.0/(Double_t)signal_pos->Integral(0,1000); 
  cout << "scale for signal template = " << scale << endl;
  signal_pos->Scale(scale);


  background_pos = (TH1D*)bkgTemplate->Clone();
  background_pos->SetName("background_pos");
  background_pos->SetLineColor(4);
  background_pos->SetMarkerColor(4);
  background_pos->SetFillColor(4);
  background_pos->SetXTitle("Fisher's isolation [GeV]");
  background_pos->Sumw2();
  scale = 1.0/(Double_t)background_pos->Integral(0,1000); 
  cout << "scale for background template = " << scale << endl;
  background_pos->Scale(scale);



  TMinuit *gMinuit = new TMinuit(1);  //initialize TMinuit with a maximum of 5 (1param??) params
  gMinuit->SetFCN(fcn); // sets function to minimize: fcn is Chi2 with errors on templates

  Double_t arglist[10];
  Int_t ierflg = 0; // status flag, it is 0 when ereything goes fine

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


  Double_t vstart = 0.5;
  Double_t step = 0.001;
  gMinuit->mnparm(0, "fsig", vstart, step, 0,1,ierflg);



  // Now ready for minimization step
  arglist[0] = 1000;
  arglist[1] = 0.01;
  gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

  Double_t fsig=0;
  Double_t fsigerr=0;
  Double_t chi2 = 0;


  if ( ierflg == 0 ) 
    {

      // Print results
      Double_t amin,edm,errdef;
      Int_t nvpar,nparx,icstat;
      gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
      gMinuit->mnprin(3,amin);
      chi2 = gMinuit->fAmin;
      gMinuit->GetParameter(0, fsig, fsigerr);  
      cout << "Fsig = " << fsig << " +- " << fsigerr << endl;

      TCanvas* c1 = new TCanvas("c1","",500,500);

      data->Draw();
      TH1D* signal_display = (TH1D*)signal_pos->Clone();
      signal_display->SetName("signal_display");
//.........这里部分代码省略.........

//___________________________________________________________________________
Double_t* IfitBin(TH1F* dataInput, TH1F* sigTemplate, TH1F* bkgTemplate, 
	       int fit_data=1)
{

  cout << "Input files are " << dataInput->GetName() << "\t" << sigTemplate->GetName() << "\t" << bkgTemplate->GetName() << endl;

  TCanvas *c1 = new TCanvas("HF1", "Histos1", 0, 0, 600, 600);
  dataCollBin.clear();
  sigCollBin.clear();
  bkgCollBin.clear();

  Double_t* fitted = new Double_t[8];
  fitted[0] = fitted[1] = fitted[2] = fitted[3] = fitted[4] = fitted[5] = fitted[6] = fitted[7] = 0.0;

  TH1F *hsum = new TH1F();
  float ntemplate = 1.;
  float sigfrac = 0.1;
  TH1F *hsum_norm = new TH1F();
  TH1F *hdata;
  TH1F *hsig  = (TH1F*)sigTemplate->Clone();
  hsig->SetName("hsig");
  hsig->Rebin(6);

  TH1F *hbkg  = (TH1F*)bkgTemplate->Clone();
  hbkg->SetName("hbkg");
  hbkg->Rebin(6);

  float ndata=0;
  if ( fit_data>0 ) {
    hdata = (TH1F*)dataInput->Clone();
    hdata -> SetName("hdata");
    ndata = hdata->Integral();
  }else {
    hsum = (TH1F*)hsig->Clone();
    hsum->Add(hbkg,1);
    cout << "For histogram " << sigTemplate->GetName() << " and " << bkgTemplate->GetName() << " sum = " << 
      hsum->Integral() << endl;

    if (hsum->Integral()>1.) ntemplate = hsum->Integral();
    sigfrac = hsig->Integral()/ntemplate;

    hsum_norm = (TH1F*)hsum->Clone();  
    hsum_norm->Scale(1./hsum->Integral());

    hdata = (TH1F*)hsum_norm->Clone();
    hdata -> SetName("hdata");
    ndata=ntemplate;
    hdata->FillRandom(hsum_norm, ndata);
  }
  if(ndata==0) {
    printf(" ---  no events in the fit \n");
    return fitted;
  }
    
  printf(" --------- before the fit ------------- \n");
  printf("Nsig %2.3f, Nbg %2.3f, Ntemplate %3.3f \n", hsig->Integral(), hbkg->Integral(), ntemplate);
//   printf("Purity %2.3f, init size %4.3f,  test sample size %4d\n", hsig->Integral()/hsum->Integral(), hsum->Integral(), ndata);
  printf(" -------------------------------------- \n");

  hdata->Rebin(6);
  int nbins = hdata->GetNbinsX();

  hsig->Scale(1./hsig->Integral());
  hbkg->Scale(1./hbkg->Integral());  

  for (int ibin=1; ibin<=nbins; ibin++) {
    dataCollBin.push_back(hdata->GetBinContent(ibin));
    sigCollBin.push_back(hsig->GetBinContent(ibin));
    bkgCollBin.push_back(hbkg->GetBinContent(ibin));    
  }
  printf( " -----  Got %d, %d, %d events for fit ----- \n ", dataCollBin.size(),
	  sigCollBin.size(), bkgCollBin.size() );  
  if ( dataCollBin.size() != sigCollBin.size() || sigCollBin.size()!=bkgCollBin.size() ) {
    printf(" error ...  inconsistent hit collection size \n");
    return fitted;
  }

  //--------------------------------------------------
  //init parameters for fit
  Double_t vstart[10] = {1., 1.};
  vstart[0] = sigfrac*ndata;
  vstart[1] = (1-sigfrac)*ndata;
 
  TMinuit *gMinuit = new TMinuit(NPARBIN);  
  gMinuit->Command("SET STR 1");
  gMinuit->SetFCN(fcnBin);
  Double_t arglist[10];
  Int_t ierflg = 0;
  
  arglist[0] = 1;
  gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
  arglist[0] = 1;
  gMinuit->mnexcm("SET PRINT", arglist ,1,ierflg);

  Double_t step[] = { 0.1, 0.1,};

  gMinuit->mnparm(0,  "Signal yield"  , vstart[0],  step[0], 0., ndata*2.  , ierflg);
  gMinuit->mnparm(1,  "background yield"  , vstart[1],  step[1], 0., ndata*2. , ierflg);
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  printf(" -------------------------------------- \n");



  printf( " -----  Got %d, %d, %d events for fit ----- \n ", dataColl.size(),
	  sigColl.size(), bkgColl.size() );  

  //--------------------------------------------------
  //init parameters for fit
  Double_t vstart[11] = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.};
  vstart[0] = sigfrac*ndata;
  vstart[1] = (1-sigfrac)*ndata;
  for (int ii=0; ii<9; ii++) {    
    vstart[ii+2] = Para[ii]; //8 shape parameters
  }
  TMinuit *gMinuit = new TMinuit(NPAR);  
  gMinuit->Command("SET STR 1");
  gMinuit->SetFCN(fcn);
  Double_t arglist[11];
  Int_t ierflg = 0;
  
  arglist[0] = 1;
  gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
  arglist[0] = 1;
  gMinuit->mnexcm("SET PRINT", arglist ,1,ierflg);

  Double_t step[] = { 1.,1.,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,};

  for ( int ii=0; ii<9; ii++){
    printf(" para %d, %.5f, err %.5f \n", ii, Para[ii], Para_err[ii]);
  }

  float sigma = 3.;
  gMinuit->mnparm(0,  "Signal yield"  , vstart[0],  step[0], 0., ndata*2.  , ierflg);
  gMinuit->mnparm(1,  "background yield"  , vstart[1],  step[1], 0., ndata*2. , ierflg);

//   gMinuit->mnparm(2,  "constant"     , Para[0],  0.00,  Para[0], Para[0] , ierflg);
//   gMinuit->mnparm(3,  "exp tail"     , Para[1],  0.01,  Para[1]-sigma*Para_err[1], Para[1]+sigma*Para_err[1], ierflg);
//   gMinuit->mnparm(4,  "exg mean"     , Para[2],  0.01,  Para[2]-sigma*Para_err[2], Para[2]+sigma*Para_err[2], ierflg);
//   gMinuit->mnparm(5,  "exg width"    , Para[3],  0.01,  Para[3]-sigma*Para_err[3], Para[3]+sigma*Para_err[3], ierflg);
//   gMinuit->mnparm(6,  "constant"     , Para[4],  0.00,  Para[4]                  , Para[4]                  , ierflg);
//   gMinuit->mnparm(7,  "bg exp turnon", Para[5],  0.01,  Para[5]-sigma*Para_err[5], Para[5]+sigma*Para_err[5], ierflg);
//   gMinuit->mnparm(8,  "bg x offset  ", Para[6],  0.01,  Para[6]-sigma*Para_err[6], Para[6]+sigma*Para_err[6], ierflg);
//   gMinuit->mnparm(9,  "bg bend slope", Para[7],  0.01,  0.001                    , 0.1                      , ierflg);
// //   gMinuit->mnparm(10, "bg bend power", Para[8],  0.01,  Para[8]-sigma*Para_err[8], Para[8]+sigma*Para_err[8], ierflg);
//   gMinuit->mnparm(10, "bg bend power", Para[8],  0.01,  0.5                       , 5.                       , ierflg);

//   gMinuit->mnparm(2,  "constant"     , Para[0], TMath::Abs(Para[0]*0.0) ,  Para[0], Para[0], ierflg);
//   gMinuit->mnparm(3,  "exp tail"     , Para[1], TMath::Abs(Para[1]*0.01) ,  Para[1]-sigma*Para_err[1], Para[1]+sigma*Para_err[1], ierflg); 
// //   gMinuit->mnparm(3,  "exp tail"     , Para[1], TMath::Abs(Para[1]*0.1) ,  0.8    , 1.3    , ierflg);
//   gMinuit->mnparm(4,  "exg mean"     , Para[2], TMath::Abs(Para[2]*0.1) ,  0.5    , 1.0    , ierflg);
//   gMinuit->mnparm(5,  "exg width"    , Para[3], TMath::Abs(Para[3]*0.1) ,  0.25   , 0.5    , ierflg);
//   gMinuit->mnparm(6,  "constant"     , Para[4], TMath::Abs(Para[4]*0.0) ,  Para[4], Para[4], ierflg);
//   gMinuit->mnparm(7,  "bg exp turnon", Para[5], TMath::Abs(Para[5]*0.1) ,  -0.7   , -0.3   , ierflg);
//   gMinuit->mnparm(8,  "bg x offset  ", Para[6], TMath::Abs(Para[6]*0.0) ,  -0.15  , -0.05  , ierflg);
//   gMinuit->mnparm(9,  "bg bend slope", Para[7], TMath::Abs(Para[7]*0.1) ,  0.01   , 0.05   , ierflg);
//   gMinuit->mnparm(10, "bg bend power", Para[8], TMath::Abs(Para[8]*0.1) ,  0.5    , 1.5    , ierflg);

  gMinuit->mnparm(2,  "constant"     , Para[0],  0.00,  Para[0], Para[0] , ierflg);
  gMinuit->mnparm(3,  "exp tail"     , Para[1],  0.00,  Para[1]-sigma*Para_err[1], Para[1]+sigma*Para_err[1], ierflg);
  gMinuit->mnparm(4,  "exg mean"     , Para[2],  0.00,  Para[2]-sigma*Para_err[2], Para[2]+sigma*Para_err[2], ierflg);
  gMinuit->mnparm(5,  "exg width"    , Para[3],  0.00,  Para[3]-sigma*Para_err[3], Para[3]+sigma*Para_err[3], ierflg);
  gMinuit->mnparm(6,  "constant"     , Para[4],  0.00,  Para[4]                  , Para[4]                  , ierflg);
  gMinuit->mnparm(7,  "bg exp turnon", Para[5],  0.00,  Para[5]-sigma*Para_err[5], Para[5]+sigma*Para_err[5], ierflg);
  gMinuit->mnparm(8,  "bg x offset  ", Para[6],  0.00,  Para[6]-sigma*Para_err[6], Para[6]+sigma*Para_err[6], ierflg);
  gMinuit->mnparm(9,  "bg bend slope", Para[7],  0.00,  0.001                    , 0.1                      , ierflg);