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

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

float getR(TMinuit& m) {
  //get results
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  m.mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  return amin;
}

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

//.........这里部分代码省略.........
      if( tmp_errpara[j] == 0. ) tmp_errpara[j] = par[j]*.1;      
    }
    //do minos if fit sucessed.

  }
  if (ierflg != 0 )  {
    printf(" *********** Fit failed! ************\n");
    gMinuit->GetParameter(0, para[0],errpara[0]);
    gMinuit->GetParameter(1, para[1],errpara[1]);
    para[0]=0.; errpara[0]=0.;

    c1->cd();
    c1->Draw();  
    //gPad->SetLogy();
    hdata->SetNdivisions(505,"XY");
    hdata->SetXTitle("comb. ISO (GeV)");
    hdata->SetYTitle("Entries");
    hdata->SetTitle("");
    hdata->SetMarkerStyle(8);
    hdata->SetMinimum(0.);
    if ( hdata->GetMaximum()<10.) hdata->SetMaximum(15.);
    else hdata->SetMaximum(hdata->GetMaximum()*1.25);
    if ( strcmp(EBEE,"EE")==0 &&ptbin == 15 ) hdata->SetMaximum(hdata->GetMaximum()*1.25);
   
    hdata->Draw("phe");  

    return fitted;    
  }

  
  // Print results
//   Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  gMinuit->mnprin(1,amin);  
  gMinuit->mnmatu(1);
  printf(" ========= happy ending !? =========================== \n");
  
  printf("FCN =  %3.3f \n", amin);

  //use new PDF form
  double tmppar[12];
  for(int ii=0; ii<9; ii++){
    tmppar[ii] = para[ii+2];
    fmcsigfit->SetParameter(ii,tmppar[ii]);
    fbkgfit->SetParameter(ii,tmppar[ii]);
  }

  c101->cd(1);
  
  //fmcsigfit->SetParameters(tmppar);
  //fmcsigfit->SetParameter(2,0.1);
  //fmcsigfit->SetLineStyle(2);

  fmcsigfit->Draw("same");
  c101->cd(2);

  fbkgfit->SetParameter(4,fbkgfit->GetParameter(4)*fmcbkg->Integral(-1., 20.)/fbkgfit->Integral(-1., 20.));
  fbkgfit->Draw("same");

  char fname[100];
  sprintf(fname,"plots/template_Ifit%s_%d.pdf",EBEE,ptbin);
  c101->SaveAs(fname);


  f11->SetParameters(tmppar);

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

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

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

//.........这里部分代码省略.........
        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
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  gMinuit->mnprin(1,amin);  
  gMinuit->mnmatu(1);
  printf(" ========= happy ending !? =========================== \n");
  
  printf("FCN =  %3.3f \n", amin);

  double yerr[100];
  for(int i=0;i<100;i++){
    yerr[i] = 0.;
  }

  hsig->Scale(para[0]);
  hbkg->Scale(para[1]);
  TH1D *hfit = (TH1D*)hbkg->Clone();
  hfit->Add(hsig);

  hsig->SetLineColor(4);
  hsig->SetLineWidth(2);
//   hsig->SetFillColor(5);
//   hsig->SetFillStyle(3001);

//   hbkg->SetLineWidth(2);
  // plot
  c1->Draw();  
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0); 
  //gPad->SetLogy();
  hdata->SetLineColor(1);

  hdata->SetXTitle(xTitle.c_str());
  hdata->SetYTitle(yTitle.c_str());
  hdata->SetTitle("");

// -- 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 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 ;
    }
  
  cout << " The final parameter values are: "<< fpar[0] << " , " << fpar[1] << endl;

  Float_t SD_frac = (fpar[1]*Hpompyt->Integral())/Hdata->Integral();
  //  Float_t NSD_events = fpar[0]*Hpythia->Integral();
  // Float_t Data_events = Hdata->Integral();

  cout << " The fraction of SD events is = " << SD_frac << endl;


  // }
  // void fill_minEE_class::show(){


  //Show the final plot

  MyC->cd(1);

  Hpythia->Scale(fpar[0]); 
  Hpompyt->Scale(fpar[1]); 

  Hdata->SetMarkerColor(2);
  Hdata->GetXaxis()->SetTitle("Min(EE-,EE+) Energy [GeV] ");
  Hdata->GetYaxis()->SetTitle("Entries");
  Hpompyt->SetLineColor(4);
  Hpythia->SetLineColor(6);
  Hpompyt->SetLineStyle(2);
  Hpythia->SetLineStyle(2);
   
  Float_t DMax = 1.2*Hdata->GetMaximum();
  Hdata->SetMaximum(DMax);
  Hdata->Draw("p");

  Hpythia->Draw("SAMES");
  Hpompyt->Draw("SAMES");
  
  HSum->Add(Hpythia);
  HSum->Add(Hpompyt);
  HSum->SetLineColor(1);
  HSum->SetLineStyle(1);
  HSum->Draw("SAMES");

}

void hrstrans3(){

//    LeRose SNAKE scaled by 1.4
    THRSTrans *trans = new THRSTrans( 0.1746, -0.1385, -0.1281, 0.050178, 0.037056, THRSTrans::kStd);
//
//
//    // Fit to data optics
//    THRSTrans *trans = new THRSTrans(  0.178259, -0.137781, -0.128674, 0.050178, 0.037056, THRSTrans::kStd);

    /*

//    trans->GetElement(5)->Print();
   trans->ShowOutput();
   trans->ShowAcc();
   trans->ShowFocalLengths();
   return;

   */

    // LeRose Translated Numbers
    //THRSTrans *trans = new THRSTrans( 0.1861, -0.1415, -0.1375, 0.5, 0.5, THRSTrans::kStd);
//   trans->ShowOutput();
  // return;
    // LeRose Numbers Fit to Transport matrix elements, PREX
//    THRSTrans *trans = new THRSTrans( 1.34761e-01, -1.38523e-01, -2.13742e-01, 0.02075 };
//
//

//     THRSTrans *trans = new THRSTrans(  1.344e-01, -1.431e-01, -1.855e-01, 0.005763, 0.000014, THRSTrans::kPREX);
//     trans->SetSeptumPsi(  1.22603, 0.929424);


    // HRS Database fit forcing 
//    THRSTrans *trans = new THRSTrans(     1.38796e-01, -1.39940e-01, -1.73034e-01, 0.02075, THRSTrans::kPREX);

//     trans->ShowOutput(4);

//     trans->ShowOutput();
//     trans->ShowFocalLengths();
//     trans->ShowAcc();
//   return;


    TMinuit *gMinuit = new TMinuit(5);
    gMinuit->SetFCN(fcn);
    Double_t arglist[10];
    Int_t ierflg = 0;

    arglist[0] = 1;
    gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
    Double_t vstart[5] = {
    // Standard tune from data
    0.178615, -0.137474, -0.128384, 0.043667, 0.036814};

    //PREX Tune
       

    Double_t step[5] = {0.001 , 0.001 , 0.001, 0.001, 0.001};

    gMinuit->mnparm(0, "q1", vstart[0], step[0], 0,0,ierflg);
    gMinuit->mnparm(1, "q2", vstart[1], step[1], 0,0,ierflg);
    gMinuit->mnparm(2, "q3", vstart[2], step[2], 0, 0,ierflg);
    gMinuit->mnparm(3, "K1", vstart[3], step[3], 0.0,1.5,ierflg);
    gMinuit->mnparm(4, "K2", vstart[4], step[4], 0.0,1.5,ierflg);


    arglist[0] = 5000;
    arglist[1] = 1.;
    gMinuit->mnexcm("SIMPLEX", arglist ,2,ierflg);
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    gMinuit->mnimpr();
    gMinuit->mnimpr();
    double q1, q2, q3, psi1, psi2, e1;

    gMinuit->GetParameter(0, q1, e1);
    gMinuit->GetParameter(1, q2, e1);
    gMinuit->GetParameter(2, q3, e1);
    gMinuit->GetParameter(3, psi1, e1);
    gMinuit->GetParameter(4, psi2, e1);


    THRSTrans *result = new THRSTrans( q1, q2, q3, psi1, psi2);

    result->ShowOutput();
    result->GetOptics()->Print();
    result->ShowAcc();
    result->ShowFocalLengths();

    return;
}

//.........这里部分代码省略.........
   if(ierflg!=0){
     cout << "+++++++++++++++++++++" << endl;
     cout << "+ ups.. check again +" << endl;
     cout << "+++++++++++++++++++++" << endl;
   }

   // MIGRAD: error optimiser
   if(ierflg==0 && 0){
     arglist[0] = 500000;
     gMinuit->mnexcm("MIGRAD", arglist ,1,ierflg);
     cout << " MIGRAD flag is " << ierflg << endl;
     if(ierflg!=0){
       cout << "+++++++++++++++++++++" << endl;
       cout << "+ ups.. check again +" << endl;
       cout << "+++++++++++++++++++++" << endl;
     }
   }

   // Minos: error optimiser
   if(ierflg==0 && 0){
     arglist[0] = 500000;
     gMinuit->mnexcm("MINOS", arglist ,1,ierflg);
     cout << " MINOS flag is " << ierflg << endl;
     if(ierflg!=0){
       cout << "+++++++++++++++++++++" << endl;
       cout << "+ ups.. check again +" << endl;
       cout << "+++++++++++++++++++++" << endl;
     }
   }

// Print results
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
   gMinuit->mnprin(3,amin);
 
   cout << "prob = " << TMath::Prob(amin,gEntriesUsed-gMinuit->GetNumFreePars())*100 << " %" << endl;
   cout << "ndf = " << gEntriesUsed-gMinuit->GetNumFreePars() << endl;

   
   cout << "output" << endl;
//   gMinuit->mnexcm("SHO COV", arglist ,0,ierflg);
//   gMinuit->mnexcm("SHO COR", arglist ,0,ierflg);
//   gMinuit->mnexcm("SHO EIG", arglist ,0,ierflg);
   
   //get fit parameters
   double par[15], epar;
   for(int i=0; i<NoPar; i++){
     gMinuit->GetParameter(i,par[i],epar);
   }

   
   TFile *output = new TFile("FitParFieldStrength_out.root","RECREATE");

   int tcount = 0;
   float tmpBAngle, tmpZenith;
   float r0x[gMaxEvents], r0y[gMaxEvents];
   float bx[gMaxEvents], by[gMaxEvents];
   
   for(int i=0; i<gEntries; i++){
     if(gFieldStrengthAnt[i]>0 && gFieldStrengthNS[i]>0 && gFieldStrengthEW[i]>0){
       tmpBAngle = gBAngle[i] * TMath::Pi() / 180.;
       tmpZenith = gZenith[i] * TMath::Pi() / 180.;
       
       r0x[tcount] = -gDistanceShowerAxis[i];
       r0y[tcount] = (log( (gFieldStrengthAnt[i]/gEffBand) / ( par[0]*cos(tmpZenith+par[3])*(par[5]+sin(tmpBAngle+par[4]))*powf(10,par[2]*(gEg[i]-8))))  );

//.........这里部分代码省略.........
  gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
  printf (" -------------------------------------------- \n");
  printf("Finished.  ierr = %d \n", ierflg);

  infoBin.clear();
  infoBin_err.clear();

  double para[NPARBIN+1],errpara[NPARBIN+1];
  if ( ierflg == 0 ) 
    {
      for(int j=0; j<=NPARBIN-1;j++) {
        gMinuit->GetParameter(j, para[j],errpara[j]);
        para[NPARBIN] = dataCollBin.size();
        infoBin.push_back(para[j]);
        infoBin_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 );
      infoBin.push_back(sigCollBin.size());

    }
  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
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;
  gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  gMinuit->mnprin(1,amin);  
  gMinuit->mnmatu(1);
  printf(" ========= happy ending !? =========================== \n");
  
  printf("FCN =  %3.3f \n", amin);

  double yerr[20];
  for(int i=0;i<20;i++){
    yerr[i] = 0.;
  }

  hsig->Scale(para[0]);
  hbkg->Scale(para[1]);
  TH1F *hfit = (TH1F*)hsig->Clone();
  hfit->Add(hbkg);


  hsig->SetLineColor(1);
  hsig->SetFillColor(5);
  hsig->SetFillStyle(3001);

  hbkg->SetLineWidth(2);
  // plot
  c1->Draw();  
  //gPad->SetLogy();
  hdata->SetLineColor(1);
  hdata->SetNdivisions(505,"XY");
  hdata->SetXTitle("Iso_{ECAL}+Iso_{HCAL}+Iso_{TRK} (GeV)");
  hdata->SetYTitle("Entries");
  hdata->SetTitle("");
  hdata->SetMarkerStyle(8);
  hdata->SetMinimum(0.);

//.........这里部分代码省略.........
      TF1* g2 = new TF1("g2","gaus",0.,4000.);
      Int_t BinDown = hdE->GetXaxis()->FindBin(maxt-2.5);
      Int_t BinUp = hdE->GetXaxis()->FindBin(maxt+2.5);
      Double_t maxy = hdE->ProjectionY()->GetMaximumBin();
      Double_t maxdE = hdE->GetYaxis()->GetBinCenter(maxy);
      hdE->ProjectionY("hproj2",BinDown,BinUp)->Fit("g2", "I", "", maxdE-500., maxdE+500);
      merry = g2->GetParameter(2);
      cuty = maxdE+5*merry;
      
      cout<<"err: "<<merrx<<" "<<merry<<endl;
      cout<<"cut: "<<cutx<<" "<<cuty<<endl;
      //Ranges of fitting 
      nbinsx = hdE->GetXaxis()->FindBin(38.);     
      nbinx0 = hdE->GetXaxis()->FindBin(maxt+4*merrx);

      //I set range of fitting in y
      Int_t up = hdE->GetYaxis()->GetNbins();
      nbinsy = up;            
      Int_t nentr = hdE->GetEntries();
      for(Int_t i=0; i<nbinsy; i++){
			if(i<maxy) continue;
			if(hdE->ProjectionX("hproj",i+1,i+1)->GetEntries()<0.0005*nentr){
	  			nbinsy = i+1; 
	}
      }
      
      cout<<"x range: "<<maxt+4*merrx<<" - "<<38.<<endl;
      cout<<"Y range: "<<hdE->GetYaxis()->GetBinCenter(nbinsy)<<endl;
     
      //3Hen peak position
      ym = maxdE;
      xm =  hdE->GetXaxis()->GetBinCenter(hdE->ProjectionX()->GetMaximumBin());
      cout<<" xm, ym: "<<xm<<" "<<ym<<endl;

      //Fit to profile - start parameters
      TF1* g0 = new TF1("g0","pol2",15.,38.);
      TProfile* pp;
      pp = hdE->ProfileX();
      pp->Fit("g0","IQ","0",maxt-5*merrx,/*hdE->GetXaxis()->GetBinCenter(nbinsx)*/38.);   
      
      //initialize TMinuit with a maximum of 3 params
      ierflg = 0;
      arglist[0] = 1;
      minuit->mnexcm("SET ERR", arglist ,1,ierflg);
      
      // Set starting values and step sizes for parameters
      if(bin<5){
		vstart[0] = g0->GetParameter(2);
		vstart[1] = g0->GetParameter(1);
	//	vstart[2] = g0->GetParameter(0);
      }

      Double_t step[3] = {0.01 , 1. , 10.};
      minuit->mnparm(0, "a", vstart[0], step[0], 0,0,ierflg);
      minuit->mnparm(1, "b", vstart[1], step[1], 0,0,ierflg);
      // minuit->mnparm(2, "c", vstart[2], step[2], 0,0,ierflg);
      
      // Now ready for minimization step
      arglist[0] = 1500;
      arglist[1] = 1;
      minuit->mnexcm("MIGRAD", arglist ,2,ierflg);
      
      // Print results
      Double_t amin,edm,errdef;
      Int_t nvpar,nparx,icstat;
      minuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
      //gMinuit->mnprin(3,amin);
      
      // get parameters
      
      double pval[2],perr[2],plo[2],phi[2];
      
      TString para0,para1,para2;
      
      int istat;
      
      minuit->mnpout(0,para0,pval[0],perr[0],plo[0],phi[0],istat);
      minuit->mnpout(1,para1,pval[1],perr[1],plo[1],phi[1],istat);
      //minuit->mnpout(2,para2,pval[2],perr[2],plo[2],phi[2],istat);
      
      fFit[bin-1][el-1] = new TF1(Form("fFit_bin%02d_el%02d",bin,el), "pol2", 0., 40.); 
      fFit[bin-1][el-1]->FixParameter(2, pval[0]);
      fFit[bin-1][el-1]->FixParameter(1, pval[1]);
      fFit[bin-1][el-1]->FixParameter(0, ym - pval[0]*xm*xm - pval[1]*xm);
    }
  }

  TCanvas* cc[6];
  for(Int_t el=0; el<6; el++){
    cc[el] = new TCanvas(Form("c_el%02d",el),Form("c_el%02d",el),1350,950);
    cc[el]->Divide(4,6);

    for(Int_t bin = 0; bin<24; bin++){  
      cc[el]->cd(bin+1);
      hdEToF[el][bin]->Draw("colz");
      fFit[el][bin]->Draw("same");
    }  
  }

}