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

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

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

void AverageMW()
{
   // include NuTeV?
   if( inclNuTeV ) nVals = 5;
   else            nVals = 4;

   // number of parameters
   const int npar = 1;
   
   // create minuit pointer
   TMinuit *minuit = new TMinuit( npar );
   
   // set FCN function
   minuit->SetFCN( &fcn );
 
   // -1=quiet, 0=normal, 1=verbose
   minuit->SetPrintLevel( 0 );

   // start values
   Double_t startVal = 80;
   Double_t fitStep  = 0.1;
   Double_t limitMin = 70;
   Double_t limitMax = 90;
   minuit->DefineParameter( 0, "MW", startVal, fitStep, limitMin, limitMax);

   // minimize with MIGRAD
   Int_t ierr = 0;
   Double_t args[2];
   args[0] = 0; // maximum function calls (0=default)
   args[1] = 0; // tolerance at minimum (0=default)
   minuit->mnexcm( "MIGrad", args, 2, ierr );
   if ( ierr != 0 )
      cout << "Warning: Maybe fit didn't converge!" << endl;

   // fit results
   Double_t fitval_MW[1], fiterr_MW[1]; 
   minuit->GetParameter( 0, fitval_MW[0], fiterr_MW[0] );

   cout << "\n\n*************************************************" << endl;
   cout << "          chi2Min = " << chi2Min  << endl;
   cout << "          n_dof   = " << nVals-1  << endl;
   cout << "          p-value = " << TMath::Prob(chi2Min, nVals-1)  << endl;
   cout << "          MW      = " << fitval_MW[0] << " +- " << fiterr_MW[0] << endl;
   cout << "*************************************************" << endl;

   // make plot
   gStyle->SetOptTitle(0);
   gStyle->SetOptStat(0);
   gStyle->SetOptFit(0);
   TCanvas* c = new TCanvas( "c", "Mean of W boson mass", 0, 0, 850, 500 );
   c->SetGridy();
   TH2F *frame = new TH2F("frame", "Mean of W boson mass", 1, 80, 80.6, 5, 0, 5. );
   frame->SetLineColor(0);
   frame->SetTickLength(0,"Y");
   frame->SetXTitle("M_{W} [GeV]");  
   frame->GetYaxis()->CenterLabels( 1 );
   frame->GetYaxis()->SetNdivisions( frame->GetNbinsY()+10, 1 );
   frame->GetYaxis()->SetBinLabel( 1, "" ); // no labels
   frame->Draw();
   
   // make labels
   MakeLabels(frame);
   
   // draw fit value
   DrawFitValue( frame, fitval_MW[0], fiterr_MW[0] );
      
   // draw measurements
   for( Int_t i=0; i<5; i++ ){
      DrawMeasurement( frame, 4-i, val_MW[i], err_MW[i] ); 
   }

   // Redraw axis
   frame->Draw("sameaxis");
   
}

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;


	}


  
}

void runQuickOptAllAntStepsLDBHPOL() {
  
  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");
  
  double startVal=0;               
  double stepSize=0.01;       
  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();
  fGeomTool->useKurtAnita3Numbers(1);
  gSystem->CompileMacro("quickOptAllAntStepsLDBHPOL.C","k");

  fillArrays(eventNumberIndex, thetaWaveIndex, phiWaveIndex, antIndex1, antIndex2, maxCorrTimeIndex, adjacent);
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(192);
  myMin->SetObjectFit(quickOptAllAntStepsLDBHPOL);
  myMin->SetFCN(iHateRoot);
 //   myMin->SetMaxIterations(2);
  //myMin->SetErrorDef(1000);

//   int ierflg;
//   double eps[1] = {2.};
//   myMin->mnexcm("SET EPS", eps, 1, ierflg);
  //setArray();

  double minValCableDelays[MAX_ANTENNAS] ={0};
  double maxValCableDelays[MAX_ANTENNAS] ={0};

  Double_t deltaR[MAX_ANTENNAS] = {0};
  Double_t deltaRErr[MAX_ANTENNAS] = {0};
  Double_t deltaZ[MAX_ANTENNAS] = {0};
  Double_t deltaZErr[MAX_ANTENNAS] = {0};
  Double_t deltaPhi[MAX_ANTENNAS] = {0};
  Double_t deltaPhiErr[MAX_ANTENNAS] = {0};
  Double_t deltaCableDelays[MAX_ANTENNAS] = {0};
  Double_t deltaCableDelaysErr[MAX_ANTENNAS] = {0};

  //  Double_t globalDeltaz = 0.242575;
  Double_t globalDeltaz = 0.03;
  
  //Min global   deltaZ[] = -0.242575 +/- 0.00198264

  for(int y = 0; y <MAX_ANTENNAS; y++){
 
    deltaR[y] = 0;//-0.05;
    deltaCableDelays[y] = 0. - (y==47)*0.35;
    minValCableDelays[y] = deltaCableDelays[y] - 0.15;
    maxValCableDelays[y] = deltaCableDelays[y] + 0.15;
    deltaZ[y] = 0;//(y<16)*(-globalDeltaz) + (y>15)*globalDeltaz;

    char name[30];
    sprintf(name,"r%d",y);
    myMin->DefineParameter(y, name, deltaR[y], stepSize, -0.3, 0.3);
    sprintf(name,"z%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS, name, deltaZ[y], stepSize, -0.3, 0.3);
    sprintf(name,"phi%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS*2, name, deltaPhi[y], stepSize, -0.3, 0.3);
    sprintf(name,"cable%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS*3, name, deltaCableDelays[y], stepSize, minValCableDelays[y], maxValCableDelays[y]);
  }



  for(int y = 0; y <MAX_ANTENNAS; y++){
    myMin->FixParameter(y); // fixed R 
    myMin->FixParameter(y+MAX_ANTENNAS); // fixed Z
    myMin->FixParameter(y+MAX_ANTENNAS*2); // fixed phi
    //    myMin->FixParameter(y+MAX_ANTENNAS*3); // fixed t

  }

  myMin->FixParameter(MAX_ANTENNAS*3); // fixed t0

  //*********MINUIT METHOD*******************
//   myMin->SetPrintLevel(-1);
  myMin->Migrad();
  //  int error_flag;
  //  myMin->mnexcm("MINOS",0,0,error_flag);  



  for(int u = 0; u <MAX_ANTENNAS; u++){
    myMin->GetParameter(u+MAX_ANTENNAS*3,deltaCableDelays[u],deltaCableDelaysErr[u]);
    cout << " deltaCableDelays[" << u << "] = " << deltaCableDelays[u]  << " +/- " << deltaCableDelaysErr[u] << endl;

    //    myMin->GetParameter(u+MAX_ANTENNAS,deltaZ[u],deltaZErr[u]);
//.........这里部分代码省略.........

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

//      gMinuit->GetParameter(1, fbk1, fbk1err );
      //cout << "Fsig = " << fsig << " +- " << fsigerr << endl;
      //cout << "Fbk1 = " << fbk1 << " +- " << fbk1err << endl;
      //cout << "Chi2/degree of freedom = " << chi2 <<endl;
	fbk1 = 1 - fsig ; 
//      TCanvas* c1 = new TCanvas("c1","",500,500);

  TH2F *htmp2 = new TH2F("htmp2","",100, 0., 0.025, 100, 0., data->GetMaximum()*1.25);

  htmp2->SetNdivisions(505,"XY");
  htmp2->SetXTitle("SigmaIetaIeta");
  htmp2->SetYTitle("A.U.");
  htmp2->SetLineColor(1); 

//     htmp2->Draw();


  TH1D* signal_display = (TH1D*)signal_pos->Clone();
  signal_display->SetName("signal_display");
  Double_t scale_sig = signal_display->Integral();
  signal_display->Scale(fsig/scale_sig);
  signal_display->SetFillStyle(3001);

      
  TH1D* background_display1 = (TH1D*)background_pos1->Clone();
  background_display1->SetName("background_display1");
  Double_t scale_background1 = background_display1->Integral();
  background_display1->Scale(fbk1/scale_background1);
  background_display1->SetFillStyle(3001);

/*

//.........这里部分代码省略.........
  
  arglist[0] = fitstra; 
  minuit->mnexcm("SET STR",arglist,1,ierflg);

  bool dofit = true; 
  
  if( dofit){
    minuit->Migrad();
    if (!minuit->fCstatu.Contains("CONVERGED")) {
      mcfitStatus = 1; //first try not converged.
      minuit->Migrad();
    }else{
      mcfitStatus = 0;
      
    }
  }
    
  double ftest[1000]; 

  double atest[1000];
  
  //minuit->GetParameter(0,fitpar[0],fitparErr[0]);

  if (!minuit->fCstatu.Contains("CONVERGED")) {
    printf("No convergence at fitting, routine Fit \n");
    printf("Minuit return string: %s \n",minuit->fCstatu.Data());
    mcfitStatus = 2;  //2nd try not converged.
  }else{
    mcfitStatus = -1;
  }
  fAmin = minuit->fAmin;

  for(int j=0; j<npar; j++){
    minuit->GetParameter(j,fitpar[j],fitparErr[j]);
    mcfitpar[j] = fitpar[j];
    mcfitparErr[j] = fitparErr[j];
  }
  
  
  double par[10];
  
  for(int j=0;j<npar;j++){
    par[j] = fitpar[j];
  }
  
  
  bool printScan = false; 
  //bool printScan = true; 
  par[0] = 0.0076; 
  
  if(printScan && testpaircat>=0){
    double  fmin = function2(par,true);
    
    
    float x1 = fitpar[testpaircat] - 0.1;
    float x2 = fitpar[testpaircat] + 0.1;
    float teststep = 0.01;
    if(smearMethod=="uncorrSmear"){
      x1 =  fitpar[testpaircat] - 0.005 >0 ? fitpar[testpaircat] - 0.005: 0;
      x2 =  fitpar[testpaircat] +0.005;
      teststep = 0.001;
    }
    
    double dmin1 = 1;
    double dmin2 = 1;
    double fitparErrLow[10];

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

//.........这里部分代码省略.........
 
  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");
//       printf("         ---------------------------------------------------------\n");
      
//       arglist[0] = 200; // number of iteration
//       arglist[1] = 1;
//       gMinuit->mnexcm("MINOS", arglist ,2,ierflg);
//       printf("         --------------------------------------------------------- \n");
//       printf("         Done Minos.  ierr = %d \n", ierflg);
//       Double_t amin;
//       gMinuit->mnprin(1,amin);
    }
  else {
    printf(" *********** Fit failed! ************\n");
    gMinuit->GetParameter(0, para[0],errpara[0]);
    gMinuit->GetParameter(1, para[1],errpara[1]);
    para[0]=0.; errpara[0]=0.;
  }

  
  // Print results

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

  int testpaircat = testcat; 
  
  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){
    minuit->FixParameter(1);
    minuit->FixParameter(2);
    minuit->FixParameter(3);

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

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

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

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

    minuit->FixParameter(4);
    minuit->FixParameter(5);
    minuit->FixParameter(6);
    
  }  
  
  arglist[0] = 0.0001;
  minuit->mnexcm("SET EPS",arglist,1,ierflg);
  
  arglist[0] = 1; 
  minuit->mnexcm("SET STR",arglist,1,ierflg);
  
  bool dofit = true; 
  
  if( dofit){
    minuit->Migrad();
    if (!minuit->fCstatu.Contains("CONVERGED")) {
      mcfitStatus = 1; //first try not converged.
      minuit->Migrad();
    }else{
      mcfitStatus = 0;
      
    }
    //minuit->mnmnos();
  }
  
  if (!minuit->fCstatu.Contains("CONVERGED")) {
    printf("No convergence at fitting, routine Fit \n");
    printf("Minuit return string: %s \n",minuit->fCstatu.Data());
    mcfitStatus = 2;  //2nd try not converged.
  }else{
    mcfitStatus = -1;
  }
  fAmin = minuit->fAmin;

  for(int j=0; j<npar; j++){
    minuit->GetParameter(j,fitpar[j],fitparErr[j]);
    mcfitpar[j] = fitpar[j];
    mcfitparErr[j] = fitparErr[j];
  }
  
  
  
  
  mytree->Fill();
  mytree->Write();
  fnew->Write();
  fnew->Close();

  
  
}