C++ TF1::GetNpar方法代码示例

TF1 *GausBF::GetBf(Int_t is1, Int_t is9) const
{
  TF1 *func = new TF1("fbf", FB(), -200, 200);
  for (Int_t i = 0; i < func->GetNpar() && i < Np; i++)
    func->SetParameter(i, GetPar(is1, is9, i));

  return func;
}

TF1* GetFit(const char *address, int lyr, double xinit) {
  const char *bgr = Form("[6]*TMath::Exp([7]*(x-%f))",xinit);
  const char *la1 = "(1-[3]-[4]-[5])*TMath::Landau(x,[1],[2],1)";
  const char *la2 = "[3]*TMath::Landau(x,2*[1]+1.4*[2],2.0*[2],1)";
  const char *la3 = "[4]*TMath::Landau(x,3*[1]+3.3*[2],3.0*[2],1)";
  const char *la4 = "[5]*TMath::Landau(x,4*[1]+5.6*[2],4.0*[2],1)";
  TF1 *ret = new TF1("fit_H1",
                     Form("[0]*(%s+%s+%s+%s)+%s",la1,la2,la3,la4,bgr) );
  ret->SetParName(0,"A");
  ret->SetParName(1,"lambda");
  ret->SetParName(2,"sigma");
  ret->SetParName(3,"f2");
  ret->SetParameter(0,1e4);  ret->SetParLimits(0,1e2,1e7);
  ret->SetParameter(1,21);   ret->SetParLimits(1,13,25.0);
  ret->SetParameter(2,3.0);  ret->SetParLimits(2,1.5,5.0);
  if(lyr>1) {
    ret->SetParameter(1,15);   ret->SetParLimits(1,2,30);
    ret->SetParameter(2,5.0);  ret->SetParLimits(2,0.8,10.0);
  }
  ret->SetParameter(3,0.20); ret->SetParLimits(3,0,0.40);
  ret->SetParameter(4,0.0);
  ret->SetParameter(5,0.0);
  ret->SetParameter(6,500); ret->SetParLimits(6,1e2,1e7);
  ret->SetParameter(7,-1);  ret->SetParLimits(7,-10,-0.1);

  ret->SetLineColor(kRed-3);

  ifstream infit;
  infit.open( address );
  double tmp;
  bool found=false;
  for(int n=0; n!=ret->GetNpar(); ++n) {
    infit >> tmp;
    if(!infit.good()) break;
    ret->SetParameter(n,tmp);
    infit >> tmp;
    found = true;
  }
  infit.close();
  if(found) {
    cout << " Previous fit results found" << endl;
    ret->SetParLimits(4,0,0.40);
    ret->SetParLimits(5,0,0.40);
  } else {
    ret->SetParLimits(4,+1,-1);
    ret->SetParLimits(5,+1,-1);
  }
  return ret;
}

void fit() {
  FILE *ofile;
  ofile = fopen("xsect-integrated-me.txt","w");
  TFile *_file0 = TFile::Open("h3maker-hn.root","update");
  _file0->Delete("*_f;*");
  TH2 *h2xsect = new TH2("hq2wXsect","Q^2:W",32,1.6,3.2,7,1.5,5.1);
  Double_t qbinedges[] = { 1.5, 1.6, 1.8, 2.1, 2.4, 2.76, 3.3, 5.1 };
  h2xsect->GetYaxis()->Set(7,qbinedges);
  TH3 *h3 = (TH3*)_file0->Get("hq2wmmp");
  int qbins = h3->GetZaxis()->GetNbins();
  int wbins = h3->GetYaxis()->GetNbins();
  fprintf(ofile, "W\tQ2\txsect\terror\tpol4p0\tpol4p1\tpol4p2\tpol4p3\tpol4p4\tgN\tgM\tgS\n");
  for (int iq = 0; iq < qbins; iq++) {
    TString hsn = TString::Format("hs%d",iq);
    THStack *hs = (THStack*)_file0->Get(hsn.Data());
    TIter next(hs->GetHists());
    //while (TObject *obj = next()) {
    //TH1 *h = (TH1*)obj;
    while (TH1 *h = (TH1*)next()) {
      float *wq = getwq(h);
      float wval = wq[0];
      float qval = wq[1];
      fitmmp(h);
      TH1 *htmp = (TH1*)h->Clone("hbgsubtracted");
      TF1 *fbg = (TF1*)h->GetListOfFunctions()->FindObject("fbg");
      htmp->Add(fbg,-1);
      double N = htmp->Integral(34,43);
      double qwidth = h3->GetZaxis()->GetBinWidth(iq+1);
      int wbin = h3->GetYaxis()->FindBin(wval);
      double wwidth = h3->GetYaxis()->GetBinWidth(wbin);
      double xsect = N/(0.891*wwidth*qwidth*19.844);
      double err2 = 0;
      for (int immp = 34; immp < 44; immp++) err2 += htmp->GetBinError(immp)*htmp->GetBinError(immp);
      //fprintf(ofile, "%.3f\t%.3f\t%.0f\t%.0f",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      fprintf(ofile, "%.3f\t%.3f\t%.3e\t%.3e",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      TF1 *ftmp = (TF1*)h->GetListOfFunctions()->At(0);
      int npar = ftmp->GetNpar();
      for (int ipar = 0; ipar < npar; ipar++) fprintf(ofile, "\t%.3e", ftmp->GetParameter(ipar));
      fprintf(ofile, "\n");
    }
    hsn.Append("_f");
    _file0->WriteObject(hs,hsn.Data());
    delete hs;
  }
  fclose(ofile);
  delete _file0;
}

TF1* GetFit(const char *run, const char *inname, const char *outname) {
  ifstream infit;
  infit.open( Form("%s/fit/%s.dat",run,inname) );
  double tmp;
  bool found;
  infit >> alp >> tmp >> lda >> elda >> sgm >> esgm >> f2 >> tmp >> f3 >> tmp >> f4;
  infit.close();
  cout << alp << endl;
  cout << lda << endl;
  cout << sgm << endl;
  cout << f2 << endl;
  cout << f3 << endl;
  cout << f4 << endl;

  double f1 = 1-f2-f3-f4;
  if(1-f2<0.01) return NULL;  
  const char *la1 = Form("%f*TMath::Landau(x,1.00*%f*[0],1.0*%f*[1],1)",f1,lda,sgm);
  const char *la2 = Form("%f*TMath::Landau(x,2.14*%f*[0],2.0*%f*[1],1)",f2,lda,sgm);
  const char *la3 = Form("%f*TMath::Landau(x,3.33*%f*[0],3.0*%f*[1],1)",f3,lda,sgm);
  const char *la4 = Form("%f*TMath::Landau(x,4.55*%f*[0],4.0*%f*[1],1)",f4,lda,sgm);
  TF1 *ret = new TF1("fit_H1", Form("%f*(%s+%s+%s+%s)",alp,la1,la2,la3,la4) );
  ret->SetParName(0,"sLDA");
  ret->SetParName(1,"SSGM");
  ret->SetParameter(0,1.0);  ret->SetParLimits(0,1e-3,2.0);
  ret->SetParameter(1,1.0);  ret->SetParLimits(1,1e-3,2.0);
  ret->SetLineColor(kRed-3);
  infit.open( Form("%s/fit/%s.dat",run,outname) );
  tmp;
  found=false;
  for(int n=0; n!=ret->GetNpar(); ++n) {
    infit >> tmp;
    if(!infit.good()) break;
    ret->SetParameter(n,tmp);
    infit >> tmp;
    found = true;
  }
  infit.close();
  if(found) {
    cout << " Previous fit results found" << endl;
  }
  cout << lda << endl;
  return ret;
}

void GACalibManager::ReadFile(string fname)
{
  try{
    ifstream fin(fname.c_str(), ios::in);
    if(!fin){
      cout << "[GACalibManager-E]: cannot open parameter file: " << fname << endl;
      throw 1; 
    }
    Int_t f_type;
    istringstream iss;
    string line;
    getline(fin, line);
    iss.str(line);
    iss >> f_type;
    Double_t xmin, xmax;
    getline(fin, line);
    iss.str(line);
    iss >> xmin;
    iss >> xmax;
    GACalibFunctions *fGACalibFunction = new GACalibFunctions();
    fGACalibrator = new GACalibrator*[n_hist];
    TF1* Func = fGACalibFunction->CalibFunc((CalFuncType)f_type);
    n_par = Func->GetNpar();
    Double_t **par;
    par = new Double_t*[n_hist];
	getline(fin, line);
    for(int i=0; i<n_hist; i++){
		par[i] = new Double_t[n_par];
		for(int j=0; j<n_par; j++)
			fin >> par[i][j];
		fGACalibrator[i] = new GACalibrator((CalFuncType)f_type, par[i], i);
		fGACalibrator[i]->SetRange(xmin,xmax);
      }
  }
  catch(bad_alloc){
    cout << "[GACalibManager-E]: Exception bad_alloc." << endl;
    throw 1;
  }
}

Double_t sr_fitError(Double_t *xx, Double_t *p) {

  assert(_sr_fitError_func);
  assert(_sr_fitError_emat);
  double x = *xx;
  double k = p[0];
  TF1 *f = _sr_fitError_func;
  int n = f->GetNpar();
  TMatrixD &emat = (*_sr_fitError_emat);
  assert(emat.GetNrows()==n);
  assert(emat.GetNcols()==n);
  
  vector<double> df(n);
  for (int i = 0; i != n; ++i) {

    double p = f->GetParameter(i);
    double dp = 0.1*f->GetParError(i);
    f->SetParameter(i, p+dp);
    double fup = f->Eval(x);
    f->SetParameter(i, p-dp);
    double fdw = f->Eval(x);
    f->SetParameter(i, p);
    df[i] = (dp ? (fup - fdw) / (2.*dp) : 0);
  }

  double sumerr2(0);
  for (int i = 0; i != n; ++i) {
    for (int j = 0; j != n; ++j) {
      sumerr2 += emat[i][j]*df[i]*df[j];
    }
  }

  double err = sqrt(sumerr2);

  return (f->Eval(x) + k*err);
}

//.........这里部分代码省略.........
        
        TF1* fvn_combinemassvnfit = new TF1(Form("fvn_combinemassvnfit_%d",i), crystalball_function_v2, fit_range_low, fit_range_high, 15);
        
        fmass_combinemassvnfit->SetLineColor(2);
        fmass_combinemassvnfit->SetLineWidth(1);
        
        fvn_combinemassvnfit->SetLineColor(2);
        fvn_combinemassvnfit->SetLineWidth(1);

        ROOT::Math::WrappedMultiTF1 wfmass_combinemassvnfit(*fmass_combinemassvnfit,1);
        ROOT::Math::WrappedMultiTF1 wfvn_combinemassvnfit(*fvn_combinemassvnfit,1);
        
        ROOT::Fit::DataOptions opt;
        ROOT::Fit::DataRange range_massfit;

        range_massfit.SetRange(fit_range_low,fit_range_high);
        ROOT::Fit::BinData datamass(opt,range_massfit);
        ROOT::Fit::FillData(datamass, h_data);
        
        ROOT::Fit::DataRange range_vnfit;
        range_vnfit.SetRange(fit_range_low,fit_range_high);
        ROOT::Fit::BinData datavn(opt,range_vnfit);
        ROOT::Fit::FillData(datavn, vn_data);
        
        ROOT::Fit::Chi2Function chi2_B(datamass, wfmass_combinemassvnfit);
        ROOT::Fit::Chi2Function chi2_SB(datavn, wfvn_combinemassvnfit);
        
        GlobalChi2_poly3bkg_floatwidth globalChi2(chi2_B, chi2_SB);

        ROOT::Fit::Fitter fitter;
        
        const int Npar = 15;
        double par0[Npar];
        for( int ipar = 0; ipar < f->GetNpar(); ipar++ ) par0[ipar] = f->GetParameter(ipar);
        par0[11] = 0.01;
        par0[12] = 0.10;
        par0[13] = 0.05;
        par0[14] = 0.01;

        
        fitter.Config().SetParamsSettings(Npar,par0);
        // fix parameter
        fitter.Config().ParSettings(0).Fix();
        fitter.Config().ParSettings(1).Fix();
        fitter.Config().ParSettings(2).Fix();
        fitter.Config().ParSettings(3).Fix();
        fitter.Config().ParSettings(4).Fix();
        fitter.Config().ParSettings(5).Fix();
        fitter.Config().ParSettings(6).Fix();
        fitter.Config().ParSettings(7).Fix();
        fitter.Config().ParSettings(8).Fix();
        fitter.Config().ParSettings(9).Fix();
        fitter.Config().ParSettings(10).Fix();
        
        fitter.Config().MinimizerOptions().SetPrintLevel(0);
        fitter.Config().SetMinimizer("Minuit2","Migrad");

        fitter.FitFCN(Npar,globalChi2,0,datamass.Size()+datavn.Size(),true);
        ROOT::Fit::FitResult result = fitter.Result();
        result.Print(std::cout);
        
        fmass_combinemassvnfit->SetFitResult( result, iparmassfit_poly3bkg_floatwidth);
        fmass_combinemassvnfit->SetRange(range_massfit().first, range_massfit().second);
        fmass_combinemassvnfit->SetLineColor(kRed);
        h_data->GetListOfFunctions()->Add(fmass_combinemassvnfit);
        //c->cd();

TF1* GeneralCorrectionFunction(
								double ChannelRangeMin, double ChannelRangeMax,double ChannelPeakPos ,
								TH2D *h2DInputForCorrection ,
								int LineIndex,
								TString InputType="T10DeriMaxEnergy",
								TString CorrNumber="1",
								double XRangeMin=0, double XRangeMax=300,
								TString FitFuncCorr="pol2",
								double FitCorrRangeMin=10, double FitCorrRangeMax= 270,
								double TresholdForCorrection=10,
								TString FitFuncSlicesString="gaus(0)+[3]+gaus(4)"
								
								
)
{
	h2DInputForCorrection->GetYaxis()->SetRangeUser(ChannelRangeMin,ChannelRangeMax);
	char buf[60];
	sprintf(buf, "hMaxCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	TH1D *hMaxPosManually=new TH1D(buf,"",h2DInputForCorrection->GetNbinsX(),h2DInputForCorrection->GetXaxis()->GetXmin(),h2DInputForCorrection->GetXaxis()->GetXmax());
	sprintf(buf, "hMaxFitCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	TH1D *hMaxPosManuallyFit=new TH1D(buf,"",h2DInputForCorrection->GetNbinsX(),h2DInputForCorrection->GetXaxis()->GetXmin(),h2DInputForCorrection->GetXaxis()->GetXmax());
	sprintf(buf, "hGausSigmaCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	TH1D *hGausSigmaManually=new TH1D(buf,"",h2DInputForCorrection->GetNbinsX(),h2DInputForCorrection->GetXaxis()->GetXmin(),h2DInputForCorrection->GetXaxis()->GetXmax());
	for(int binX = h2DInputForCorrection->GetXaxis()->FindBin(XRangeMin);binX <= h2DInputForCorrection->GetXaxis()->FindBin(XRangeMax);binX++)
	//for(int binX = h2DInputForCorrection->GetXaxis()->FindBin(200);binX <= h2DInputForCorrection->GetXaxis()->FindBin(200);binX++)
	{
		TH1D *hProfileY =h2DInputForCorrection->ProjectionY("_py",binX,binX);	
		double MaxValue=hProfileY->GetBinCenter(hProfileY->GetMaximumBin());
		
		//hMaxPosManually->SetBinContent(binX, MaxValue);
		//h2DInputForCorrection
		//cout <<hProfileY->GetEntries()<<endl;
		//TF1* FitFuncSlices = new TF1("FitFuncSlices","gaus(0)+[3]",MaxValue-20,MaxValue+20);
		//cout << TMath::Max(MaxValue-20,double(ChannelRangeMin)) << "\t" << TMath::Min(MaxValue+20,double(ChannelRangeMax)) << "\t"<<endl;
		TF1* FitFuncGausSlices = new TF1("FitFuncGausSlices","gaus(0)",TMath::Max(MaxValue-20,double(ChannelRangeMin)),TMath::Min(MaxValue+20,double(ChannelRangeMax)));
		FitFuncGausSlices->SetParameters(hProfileY->GetBinContent(hProfileY->GetMaximumBin()),MaxValue,4);
		
		hProfileY->Fit(FitFuncGausSlices,"RNIQ");
		TF1* FitFuncSlices = new TF1("FitFuncSlices",FitFuncSlicesString.Data(),TMath::Max(MaxValue-3*FitFuncGausSlices->GetParameter(2),double(ChannelRangeMin)),TMath::Min(MaxValue+3*FitFuncGausSlices->GetParameter(2),double(ChannelRangeMax)));
		FitFuncSlices->SetParameters(FitFuncGausSlices->GetParameter(0),FitFuncGausSlices->GetParameter(1),FitFuncGausSlices->GetParameter(2),10,10,FitFuncGausSlices->GetParameter(1)-5,5);
		
		FitFuncSlices->SetParLimits(0,FitFuncGausSlices->GetParameter(0)*0.8,FitFuncGausSlices->GetParameter(0)*1.5);
		FitFuncSlices->SetParLimits(1,TMath::Max(FitFuncGausSlices->GetParameter(1)-10,double(ChannelRangeMin)),TMath::Min(FitFuncGausSlices->GetParameter(1)+10,double(ChannelRangeMax)));
		FitFuncSlices->SetParLimits(2,0,FitFuncGausSlices->GetParameter(2)*2);
		FitFuncSlices->SetParLimits(3,0,500);
		FitFuncSlices->SetParLimits(4,0,FitFuncGausSlices->GetParameter(0)*0.3);
		
		FitFuncSlices->SetParLimits(5,TMath::Max(FitFuncGausSlices->GetParameter(1)-10,double(ChannelRangeMin)),TMath::Min(MaxValue-1,double(ChannelRangeMax)));
		FitFuncSlices->SetParLimits(6,0,10);
		hProfileY->Fit(FitFuncSlices,"RINQ");
		//hProfileY->DrawCopy();
		
		//cout <<MaxValue<<"  " << FitFuncSlices->GetParameter(1) << "   " << FitFuncSlices->GetParError(1) <<endl;
		//cout <<MaxValue<<"  " << FitFuncSlices->GetParameter(1) << "   " << FitFuncSlices->GetMaximumX() <<endl;
		hMaxPosManually->SetBinContent(binX, (FitFuncSlices->GetParameter(1))/ChannelPeakPos);
		hMaxPosManually->SetBinError(binX, FitFuncSlices->GetParError(1)/ChannelPeakPos);
		hGausSigmaManually->SetBinContent(binX, FitFuncSlices->GetParameter(2));
		hGausSigmaManually->SetBinError(binX, FitFuncSlices->GetParError(2));
		if(FitFuncSlices->GetParameter(2)<TresholdForCorrection && FitFuncSlices->GetParError(2)<5)
		{
			hMaxPosManuallyFit->SetBinContent(binX, (FitFuncSlices->GetParameter(1))/ChannelPeakPos);
			hMaxPosManuallyFit->SetBinError(binX, FitFuncSlices->GetParError(1)/ChannelPeakPos);
		}
		//hSpectrumTDeriMax1090Rel_EnergyChannel_MaxPosManually->SetBinError(binX, FitFuncSlices->GetParameter(2)/ChannelPeakPos);
		hProfileY->Delete();
		//cin.ignore();
	}
	//write histos to file
		
	//sprintf(buf, "hMaxCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	hMaxPosManually->Write(0,TObject::kOverwrite);
	//sprintf(buf, "hGausSigmaCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	hGausSigmaManually->Write(0,TObject::kOverwrite);
	hMaxPosManuallyFit->Write(0,TObject::kOverwrite);
	sprintf(buf, "funcCorr%s_%sNorm_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	//fit corr function and write it to file
	TF1 *funcCorrNorm = new TF1(buf,FitFuncCorr.Data(),FitCorrRangeMin,FitCorrRangeMax);
	funcCorrNorm->SetParameters(1,0,-0);
	
	funcCorrNorm->SetParLimits(0,0.8,1);
	//funcCorrNorm->SetParLimits(2,-1E5,0);
	//if(LineIndex==2)
	//	hMaxPosManuallyFit->Fit(funcCorrNorm,"RBI");
	//else
		hMaxPosManuallyFit->Fit(funcCorrNorm,"RBQI");
	sprintf(buf, "funcCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	TF1 *funcCorr = new TF1(buf,FitFuncCorr.Data(),FitCorrRangeMin,FitCorrRangeMax);
	for(int i= 0; i<funcCorr->GetNpar();i++)
	{
		funcCorr->SetParameter(i,funcCorrNorm->GetParameter(i)*ChannelPeakPos);
	}
	//sprintf(buf, "funcCorr%s_%sNorm_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	funcCorrNorm->Write(0,TObject::kOverwrite);
	//sprintf(buf, "funcCorr%s_%s_%d",CorrNumber.Data(),InputType.Data(),LineIndex);
	funcCorr->Write(0,TObject::kOverwrite);
	h2DInputForCorrection->GetYaxis()->UnZoom();
	return funcCorr;
}

Double_t GausBF::Gfit(Int_t is1, Int_t is9, Double_t rgt, TGraph *gr,
		                                          TGraph *gb)
{
  if ((is1 == 59 && is9 == 98) || (is1 == 59 && is9 == 109) ||
      (is1 == 73 && is9 == 98) || (is1 == 86 && is9 ==  97) ||
      (is1 == 86 && is9 == 98) || (is1 == 97 && is9 ==  98)) return -1;

  AMSPoint p1 = GetP1(is1);
  AMSPoint p9 = GetP9(is9);

  if (rgt == 0) {
    TGraph grb; if (!gb) gb = &grb;
    TrProp trp(p1, p9-p1, rgt);
    for (Int_t i = 0; i < 200; i++) {
      Double_t z = -200+(i+0.5)*2;
      trp.Propagate(z);
      gb->SetPoint(i, z, trp.GuFld(trp.GetP0()).x());
    }

    TF1 *func = Bfit(gb);

    Double_t *par = GausBF::Head()->GetPar(is1, is9);
    for (Int_t i = 0; i < func->GetNpar() && i < Np; i++)
      par[i] = func->GetParameter(i);

    if (par[0] == 0) cout << "Failed: " << is1 << " " << is9 << endl;

    for (Int_t i = 0; i < Np/3; i++)
      if (par[i*3] != 0) {
	if (par[i*3+2] < 0) par[i*3+2] *= -1;
	if (par[i*3+2] < 1 || TMath::Abs(par[i*3+1]) > 200) {
	  par[i*3] = par[i*3+1] = 0;
	             par[i*3+2] = 1;
	}
      }

    Double_t err = 1e-3;
    Double_t csq = 0;
    for (Int_t i = 0; i < gb->GetN(); i++) {
      Double_t d = gb->GetY()[i]-func->Eval(gb->GetX()[i]);
      csq += d*d/err/err;
    }
    return csq;
  }

  TGraph grp; if (!gr) gr = &grp;
  TrProp trp(p1, p9-p1, rgt);

  AMSPoint pnt;
  AMSDir   dir;
  if (_Zref == 0) _Zref = TkDBc()->GetZlayerAJ(1);

  for (Int_t i = 0; i < 9; i++) {
    Double_t z = TkDBc()->GetZlayerAJ(i+1);
    trp.Propagate(z);
    gr->SetPoint(i, z, trp.GetP0y());
    if (i == 0) { pnt = trp.GetP0(); dir = trp.GetDir(); }
  }

  TF1 *func = GausBF::Head()->GetPr(is1, is9, rgt); 
  gr->Fit(func, "q0");

  for (Int_t i = 0; i < Np; i++)
    GausBF::Head()->SetPar(is1, is9, i, func->GetParameter(i+5));

  Double_t csq = 0;
  for (Int_t i = 0; i < gr->GetN(); i++) {
    Double_t res = (gr->GetY()[i]-func->Eval(gr->GetX()[i]))*1e4;
//  cout << Form("%6.1f %7.3f %7.3f  %6.1f", gr->GetX()[i], gr->GetY()[i],
//		 func->Eval(gr->GetX()[i]), res) << endl;
    gr->GetY()[i] = res;
    csq += res*res;
  }
  delete func;

  return csq/gr->GetN();
}

void fit(float bgpar2smudge=1.0) {
  TF1 *fsig = new TF1("fsig",&d_sig,0.4,2,5);
  FILE *ofile;
  ofile = fopen("xsect-integrated-me.txt","w");
  TFile *_file0 = TFile::Open("h3maker-hn.root","update");
  _file0->Delete("*_f;*");
  TH2 *h2xsect = new TH2("hq2wXsect","Q^2:W",32,1.6,3.2,7,1.5,5.1);
  Double_t qbinedges[] = { 1.5, 1.6, 1.8, 2.1, 2.4, 2.76, 3.3, 5.1 };
  h2xsect->GetYaxis()->Set(7,qbinedges);
  TH3 *h3 = (TH3*)_file0->Get("hq2wmmp");
  int qbins = h3->GetZaxis()->GetNbins();
  //int wbins = h3->GetYaxis()->GetNbins();
  fprintf(ofile, "W\tQ2\txsect\terror\tpol4p0\tpol4p1\tpol4p2\tpol4p3\tpol4p4\tgN\tgM\tgS\tstepx0\tstepx1\txsectFn\n");
  for (int iq = 0; iq < qbins; iq++) {
    TString hsn = TString::Format("hs%d",iq);
    THStack *hs = (THStack*)_file0->Get(hsn.Data());
    TIter next(hs->GetHists());
    //while (TObject *obj = next()) {
    //TH1 *h = (TH1*)obj;
    while (TH1 *h = (TH1*)next()) {
      float *wq = getwq(h);
      float wval = wq[0];
      float qval = wq[1];
      delete [] wq;
      int wbin = h3->GetYaxis()->FindBin(wval);
      float wlow = h3->GetYaxis()->GetBinLowEdge(wbin);
      float step_x0 = sqrt(wlow*wlow+MASS_P*MASS_P-2*wlow*MASS_P);
      float whigh = h3->GetYaxis()->GetBinLowEdge(wbin+1);
      float step_x1 = sqrt(whigh*whigh+MASS_P*MASS_P-2*whigh*MASS_P);
      fitmmp(h,step_x0,step_x1,wval);
      TH1 *htmp = (TH1*)h->Clone("hbgsubtracted");
      TF1 *fbg = (TF1*)h->GetListOfFunctions()->FindObject("fbg");
      htmp->Add(fbg,-1);
      double N = htmp->Integral(34,43);
      double qwidth = h3->GetZaxis()->GetBinWidth(iq+1);
      //int wbin = h3->GetYaxis()->FindBin(wval);
      double wwidth = h3->GetYaxis()->GetBinWidth(wbin);

      TF1 *ftmp = (TF1*)h->GetListOfFunctions()->At(0);
      fsig->SetParameter(0,ftmp->GetParameter(5));
      fsig->SetParameter(1,ftmp->GetParameter(6));
      fsig->SetParameter(2,ftmp->GetParameter(7));
      fsig->SetParameter(3,step_x0);
      fsig->SetParameter(4,step_x1);
      fsig->SetLineWidth(2);
      fsig->SetLineColor(kBlue+1);
      h->GetListOfFunctions()->Add((TF1*)fsig->Clone("fsig"));
      //fsig->Print();
      double Nfn = 0;
      for (int b = 1; b < h->GetNbinsX(); b++) {
        double x = h->GetXaxis()->GetBinCenter(b);
        Nfn += fsig->Eval(x);
      }
      //printf("**** %.3e\t\%.3e\n",Nfn,N);
      double xsect = N/(0.891*wwidth*qwidth*19.844);
      double xsectFn = Nfn/(0.891*wwidth*qwidth*19.844);
      double err2 = 0;
      for (int immp = 34; immp < 44; immp++) err2 += htmp->GetBinError(immp)*htmp->GetBinError(immp);
      //fprintf(ofile, "%.3f\t%.3f\t%.0f\t%.0f",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      fprintf(ofile, "%.3f\t%.3f\t%.3e\t%.3e",wval,qval,xsect/(1e6), sqrt(err2)/(1e6));
      int npar = ftmp->GetNpar();
      for (int ipar = 0; ipar < npar; ipar++) fprintf(ofile, "\t%.3e", ftmp->GetParameter(ipar));
      fprintf(ofile,"\t%.3e",xsectFn/(1e6));
      fprintf(ofile, "\n");
    }
    hsn.Append("_f");
    _file0->WriteObject(hs,hsn.Data());
    delete hs;
  }
  fclose(ofile);
  delete _file0;
}

//.........这里部分代码省略.........
	histo_massfit->Fit(Form("f_%s_%d",cfgname.Data(),ipt),"L q","",fit_range_low,fit_range_high);
	histo_massfit->Fit(Form("f_%s_%d",cfgname.Data(),ipt),"L q","",fit_range_low,fit_range_high);
	histo_massfit->Fit(Form("f_%s_%d",cfgname.Data(),ipt),"L m","",fit_range_low,fit_range_high);

	//begin combine fit
    TF1* fmass_combinemassvnfit = new TF1(Form("fmass_combinemassvnfit_%s_%d",cfgname.Data(),ipt),"[0]*([5]*([4]*TMath::Gaus(x,[1],[2]*(1.0 +[6]))/(sqrt(2*3.14159)*[2]*(1.0 +[6]))+(1-[4])*TMath::Gaus(x,[1],[3]*(1.0 +[6]))/(sqrt(2*3.14159)*[3]*(1.0 +[6])))+(1-[5])*TMath::Gaus(x,[8],[7]*(1.0 +[6]))/(sqrt(2*3.14159)*[7]*(1.0 +[6]))) + [9] + [10]*x + [11]*x*x + [12]*x*x*x", fit_range_low, fit_range_high);

	TF1* fvn_combinemassvnfit = new TF1(Form("fvn_combinemassvnfit_%s_%d",cfgname.Data(),ipt), "( ( [0]*([5]*([4]*TMath::Gaus(x,[1],[2]*(1.0 +[6]))/(sqrt(2*3.14159)*[2]*(1.0 +[6]))+(1-[4])*TMath::Gaus(x,[1],[3]*(1.0 +[6]))/(sqrt(2*3.14159)*[3]*(1.0 +[6])))+(1-[5])*TMath::Gaus(x,[8],[7]*(1.0 +[6]))/(sqrt(2*3.14159)*[7]*(1.0 +[6]))) ) / ( [0]*([5]*([4]*TMath::Gaus(x,[1],[2]*(1.0 +[6]))/(sqrt(2*3.14159)*[2]*(1.0 +[6]))+(1-[4])*TMath::Gaus(x,[1],[3]*(1.0 +[6]))/(sqrt(2*3.14159)*[3]*(1.0 +[6])))+(1-[5])*TMath::Gaus(x,[8],[7]*(1.0 +[6]))/(sqrt(2*3.14159)*[7]*(1.0 +[6]))) + [9] + [10]*x + [11]*x*x + [12]*x*x*x ) ) * [13] + ( 1.0 - ( ( [0]*([5]*([4]*TMath::Gaus(x,[1],[2]*(1.0 +[6]))/(sqrt(2*3.14159)*[2]*(1.0 +[6]))+(1-[4])*TMath::Gaus(x,[1],[3]*(1.0 +[6]))/(sqrt(2*3.14159)*[3]*(1.0 +[6])))+(1-[5])*TMath::Gaus(x,[8],[7]*(1.0 +[6]))/(sqrt(2*3.14159)*[7]*(1.0 +[6]))) ) / ( [0]*([5]*([4]*TMath::Gaus(x,[1],[2]*(1.0 +[6]))/(sqrt(2*3.14159)*[2]*(1.0 +[6]))+(1-[4])*TMath::Gaus(x,[1],[3]*(1.0 +[6]))/(sqrt(2*3.14159)*[3]*(1.0 +[6])))+(1-[5])*TMath::Gaus(x,[8],[7]*(1.0 +[6]))/(sqrt(2*3.14159)*[7]*(1.0 +[6]))) + [9] + [10]*x + [11]*x*x + [12]*x*x*x ) ) ) * ( [14] + [15] * x + [16] * x * x)", fit_range_low, fit_range_high);

	ROOT::Math::WrappedMultiTF1 wfmass_combinemassvnfit(*fmass_combinemassvnfit,1);
	ROOT::Math::WrappedMultiTF1 wfvn_combinemassvnfit(*fvn_combinemassvnfit,1);

	ROOT::Fit::DataOptions opt;
	ROOT::Fit::DataRange range_massfit;
	// set the data range
	range_massfit.SetRange(fit_range_low,fit_range_high);
	ROOT::Fit::BinData datamass(opt,range_massfit);
	ROOT::Fit::FillData(datamass, histo);

	ROOT::Fit::DataRange range_vnfit;
	range_vnfit.SetRange(fit_range_low,fit_range_high);
	ROOT::Fit::BinData datavn(opt,range_vnfit);
	ROOT::Fit::FillData(datavn, h_vnvsmass);

	ROOT::Fit::Chi2Function chi2_B(datamass, wfmass_combinemassvnfit);
	ROOT::Fit::Chi2Function chi2_SB(datavn, wfvn_combinemassvnfit);

	GlobalChi2_poly3bkg_floatwidth_poly2bkg globalChi2(chi2_B, chi2_SB);

	ROOT::Fit::Fitter fitter;

	const int Npar = 17;
	double par0[Npar];
	for( int ipar = 0; ipar < f->GetNpar(); ipar++ )
		par0[ipar] = f->GetParameter(ipar);
	par0[13] = 0.01;
	par0[14] = 0.10;
	par0[15] = 0.05;
	par0[16] = 0.05;

	// create before the parameter settings in order to fix or set range on them
	fitter.Config().SetParamsSettings(Npar,par0);
	// fix parameter
	fitter.Config().ParSettings(2).Fix();
	fitter.Config().ParSettings(3).Fix();
	fitter.Config().ParSettings(4).Fix();
	fitter.Config().ParSettings(5).Fix();
	fitter.Config().ParSettings(7).Fix();
	fitter.Config().ParSettings(8).Fix();
	// set limits on the third and 4-th parameter
	fitter.Config().ParSettings(1).SetLimits(1.7, 2.0);
	//fitter.Config().ParSettings(12).SetStepSize(0.005);
	//fitter.Config().UseWeightCorrection();

	fitter.Config().MinimizerOptions().SetPrintLevel(0);
	fitter.Config().SetMinimizer("Minuit2","Migrad");

	// fit FCN function directly
	// (specify optionally data size and flag to indicate that is a chi2 fit)
	//fitter.FitFCN(Npar,globalChi2,0,datamass.Size()+datavn.Size(),false);
	fitter.FitFCN(Npar,globalChi2,0,datamass.Size()+datavn.Size(),true);
	ROOT::Fit::FitResult result = fitter.Result();
	result.Print(std::cout);

	fmass_combinemassvnfit->SetFitResult( result, iparmassfit_poly3bkg_floatwidth_poly2bkg);
	fmass_combinemassvnfit->SetRange(range_massfit().first, range_massfit().second);

//.........这里部分代码省略.........
      gPad->SetLogx();
      h3->SetMaximum(imethod==0 ? 0.05 : (itype!=2 ? 0.1 : 0.25));
      h3->SetMinimum(imethod==0 ? -0.05 : (itype!=2 ? -0.4 : -0.25));
      h3->SetYTitle(Form("k_{FSR} = dR/d#alpha (%s)",ct));
      h3->DrawClone("AXIS");
      tex->DrawLatex(0.20,0.85,texlabel[cm]);
      tex->DrawLatex(0.20,0.80,"|#eta| < 1.3");
      TLegend *leg = tdrLeg(0.60,0.75,0.90,0.90);

      for (int  isample = 0; isample != nsamples; ++isample) {

	const char *cs = samples[isample];
	TGraphErrors *gk = gkmap[ct][cm][cs]; assert(gk);
	
	leg->AddEntry(gk,texlabel[cs],"P");

	// Fit each sample separately for pT balance
	if (true) {

	  TF1 *fk = new TF1(Form("fk_%s_%s_%s",ct,cm,cs),
			    "[0]+[1]*log(0.01*x)+[2]*pow(log(0.01*x),2)",
			    30,1300);
	  fk->SetParameters(-0.25,-0.5);
	  fk->SetLineColor(gk->GetLineColor());
	  gk->Fit(fk, "QRN");

	  tex->SetTextColor(fk->GetLineColor());
	  tex->DrawLatex(0.55,0.27-0.045*isample,
			 Form("#chi^{2}/NDF = %1.1f / %d",
			      fk->GetChisquare(), fk->GetNDF()));
	  tex->SetTextColor(kBlack);

	  // Error band
	  const int n = fk->GetNpar();
	  TMatrixD emat(n,n);
	  gMinuit->mnemat(emat.GetMatrixArray(), n);
	  TF1 *fke = new TF1(Form("fk_%s_%s_%s",ct,cm,cs),
			     sr_fitError, 30, 1300, 1);
	  _sr_fitError_func = fk;
	  _sr_fitError_emat = &emat;

	  fke->SetLineStyle(kSolid);
	  fke->SetLineColor(fk->GetLineColor()-10);
	  fke->SetParameter(0,-1);
	  fke->DrawClone("SAME");
	  fke->SetParameter(0,+1);
	  fke->DrawClone("SAME");

	  fk->DrawClone("SAME");
	  gk->DrawClone("SAME Pz");

	  // Store soft radiation corrections in fsr subdirectory
	  assert(fin->cd(ct));
	  assert(gDirectory->cd(bin));
	  if (!gDirectory->FindObject("fsr")) gDirectory->mkdir("fsr");
	  assert(gDirectory->cd("fsr"));

	  TH1D *hk = (TH1D*)(isample==0 ? hpt2->Clone() : hpt1->Clone());
	  hk->SetName(Form("hkfsr_%s_%s",cm,cs));
	  TProfile *ppt = (isample==0 ? ppt2 : ppt1);
	  if (isample==3) { ppt = ppt4; } // pas-v6
	  for (int i = 1; i != hk->GetNbinsX()+1; ++i) {
	    double pt = ppt->GetBinContent(i);
	    if (pt>30 && pt<1300) {
	      hk->SetBinContent(i, fk->Eval(pt));
	      hk->SetBinError(i, fabs(fke->Eval(pt)-fk->Eval(pt)));

//***############## main fitting Fxn ################ *****//
void FitPlotAndSave( char *Ifile ){
 
/** Plot Options***/	
//gROOT->Reset();
// gROOT->Clear();
gROOT->SetStyle("Plain") ;
gROOT->SetBatch(kFALSE);
gStyle->SetOptTitle(1);
gStyle->SetOptStat(0);
gStyle->SetOptFit(1);
gStyle->SetStatX(.89);
gStyle->SetStatY(.89) ;
gStyle->SetStatBorderSize(0);
//gStyle->SetOptStat(1111111)
gStyle->SetCanvasColor(kWhite);   // background is no longer mouse-dropping white
gStyle->SetPalette(1);        // blue to red false color palette. Use 9 for b/w
gStyle->SetCanvasBorderMode(0);     // turn off canvas borders
gStyle->SetPadBorderMode(0);
gStyle->SetPaintTextFormat("5.2f");  // What precision to put numbers if plotted with "TEXT"

// For publishing:
gStyle->SetLineWidth(2);
gStyle->SetTextSize(1.1);
gStyle->SetLabelSize(0.06,"xy");
gStyle->SetTitleSize(0.08,"xy");
gStyle->SetTitleOffset(1.2,"x");
gStyle->SetTitleOffset(1.0,"y");
gStyle->SetPadTopMargin(0.1);
gStyle->SetPadRightMargin(0.1);
gStyle->SetPadBottomMargin(0.16);
gStyle->SetPadLeftMargin(0.12);
TGaxis::SetMaxDigits(1); // Set Axis to be of the form 0.11 10^N


       TFile *ifile  = new TFile(Ifile);
        
	TF1 *fitFcn  = new TF1("fitFcn", mygaus, FitLowRange, FitHighRange, 3 );
	fitFcn->SetNpx(500);
	fitFcn->SetLineWidth(4);
	fitFcn->SetLineStyle(5);
	fitFcn->SetLineColor(kBlue);
        cout <<" Calling Fitting Fxntion" << endl;
	TH1F*h_Seed_TimeEBEB = (TH1F*)ifile->Get("EBEB/seed time");
	if(h_Seed_TimeEBEB == 0){ std::cout  <<"!! Histogram Does not exist!!" << std::endl; throw 1;}
        
	h_Seed_TimeEBEB->SetTitle("Seed Time[ns]");   
        h_Seed_TimeEBEB->SetMarkerStyle(20);
        h_Seed_TimeEBEB->SetMarkerSize(0.8);
        h_Seed_TimeEBEB->SetStats(1);
        h_Seed_TimeEBEB->SetTitleSize(0.08, "x");   
        h_Seed_TimeEBEB->SetTitleOffset(1.0, "x");    
        h_Seed_TimeEBEB->SetTitleSize(0.06, "y"); 
        h_Seed_TimeEBEB->SetTitleOffset(0.95, "y");    
        h_Seed_TimeEBEB->SetYTitle("Number of Seeds/0.05ns"); 
        h_Seed_TimeEBEB->SetXTitle("t_{seed}[ns]"); 
        h_Seed_TimeEBEB->GetXaxis()->SetRangeUser(FitLowRange, FitHighRange);   

       /** Set parms as parms of Fit Fxn **/
	fitFcn->SetParameters(500, h_Seed_TimeEBEB->GetMean(), h_Seed_TimeEBEB->GetRMS() );
	fitFcn->SetParNames("CONST", "#mu(ns)", "#sigma(ns)");
	h_Seed_TimeEBEB->Fit("fitFcn", "LL"); /**Fit with improved LL**/
	std::cout << "Printing Fit Parameters for EBEB ......   " << std::endl;
        printf("Integral of function in EBEB = %g\n", fitFcn->Integral( FitLowRange, FitHighRange));

        //*** retrive fit results***//
        int npar = fitFcn->GetNpar();
        TVirtualFitter *fit = TVirtualFitter::GetFitter();
        fit->PrintResults(2,0.);
        TMatrixD *CovMatrix = new TMatrixD ( npar, npar, fit->GetCovarianceMatrix() );
	CovMatrix->Print();
        TCanvas *c1 = new TCanvas("c1","EB-EB",200,10,800,900);
     	c1->SetGridx();
     	c1->SetGridy();
	c1->GetFrame()->SetFillColor(21);
	c1->GetFrame()->SetBorderMode(-1);
	c1->GetFrame()->SetBorderSize(5);
        /* c1->Divide(2,1);  */
	c1->cd();
	h_Seed_TimeEBEB->Draw();
	fitFcn->Draw("sames");
        c1->SetLogy(0);
	//  draw the legend
    	TLegend *leg = new TLegend(0.15,0.72,0.3,0.85);
       	leg->SetTextFont(72);
        leg->SetTextSize(0.04);
        leg->AddEntry(h_Seed_TimeEBEB,"EB","lpe");
        leg->AddEntry(fitFcn,"GAUS","l");
        leg->Draw();
	c1->SaveAs("Seed_Time_DoubleElectron_Run2012A-EB-EB.png");
}       

//.........这里部分代码省略.........
  fjes->SetLineWidth(2);
  fjes->SetRange(10.,3500.);
  fjes->Draw("SAME");
  
  //TF1 *ft = new TF1("ft","1-[0]-[1]*pow(x,[2]) + ([3]+[4]*log(x))/x",30,2200);
  //ft->SetParameters(0,0.05,-0.5,1,0.1);
  //ft->FixParameter(3,0);

  // Logarithmic sigmoid
  //TF1 *ft = new TF1("ft","[0]+(1-[0])/(1. + exp(-(log(x)-log(abs([1])))"
  //	       "/(log(abs([2])+abs([1]))-log(abs([1])))))", 30,2200);
  //ft->SetParameters(0.98, 150, 50);
  TF1 *ft = new TF1("ft","[0]+(1-[0])/(1. + exp(-(log(x)-[1])/[2]))",30,2200);
  //ft->SetParameters(0.98,log(145),log(190)-log(145));
  //ft->SetParameters(0.982,4.967,0.271);
  //ft->SetParameters(0.976,5.040,0.370); // ENDCAP
  //ft->SetParameters(0.985,5.0,0.3);
  ft->SetParameters(0.985,5.025,0.3);
  //ft->FixParameter(1,5.03); // semi-weighted average of BCD and EF
  //ft->FixParameter(2,0.395); // combined fit to BCD+EF / G+H 

  // ( 12.9*5.055+6.8*5.000)/(12.9+6.8)
  ft->FixParameter(1,5.036); // semi-weighted average of BCD/GH and EF/GH
  // ( 12.9*0.344 + 6.8*0.455)/(12.9+6.8)
  ft->FixParameter(2,0.391); // combined fit to BCD+EF / GH 

  // Log-sigmoid + powerlaw
  //TF1 *ft = new TF1("ft","[0]+(1-[0])/(1. + exp(-(log(x)-[1])/[2]))"
  //	       "*(1-[3]*pow(x,[4]))",30,2200);
  //ft->SetParameters(0.982,4.967,0.271,0.1,-0.2);
  // Double powerlaw
  //TF1 *ft = new TF1("ft","[4]-[0]*pow(x,[1])-[2]*pow(x,[3])",30,2200);
  //ft->SetParameters(0.05,-0.15,0.01,-0.3,1);
  

  mg->Fit(ft,"RN");
  ft->SetLineColor(kBlue);
  ft->SetLineWidth(2);
  ft->SetRange(10.,3500.);
  ft->Draw("SAME");

  // Map multijet with response ratio
  if (gmpf) { // we have multijet available
    TGraphErrors *gmpf2 = (TGraphErrors*)gmpf->Clone("gmpf2");
    gmpf2->SetMarkerColor(kBlack);//kGray+1);
    gmpf2->SetLineColor(kBlack);//kGray+1);
    for (int i = 0; i != gmpf->GetN(); ++i) {
      if (mjvsjes) {
	gmpf2->SetPoint(i, 0.4*gmpf->GetX()[i],
			fjes->Eval(gmpf->GetX()[i])/gmpf->GetY()[i]);
	gmpf2->SetPointError(i, 0.4*gmpf->GetEX()[i],
			     gmpf->GetEY()[i]);
      }
      else {
	gmpf2->SetPoint(i, 0.4*gmpf->GetX()[i],
			ft->Eval(gmpf->GetX()[i])/gmpf->GetY()[i]);
	gmpf2->SetPointError(i, 0.4*gmpf->GetEX()[i],
			     gmpf->GetEY()[i]);
      }
    }
    gmpf2->Draw("SAMEPz");
  } // multijet

  tex->SetTextColor(kBlue);
  tex->DrawLatex(0.50,0.85,Form("#chi^{2} / NDF = %1.1f / %d",
				ft->GetChisquare(),
				ft->GetNDF()));
  tex->SetTextColor(kBlack);
  tex->SetTextSize(0.040);
  tex->DrawLatex(0.50,0.80,Form("(#chi^{2} / NDF = %1.1f / %d)",
				fjes->GetChisquare(),
				fjes->GetNDF()));


  tex->SetTextColor(kBlue-9);
  tex->SetTextSize(0.030);
  tex->DrawLatex(0.20,0.25,ft->GetExpFormula());
  tex->DrawLatex(0.20,0.20,
		 Form("p_{0}=%1.3f#pm%1.3f"
		      ", p_{1}=%1.3f#pm%1.3f"
		      ", p_{2}=%1.3f#pm%1.3f",
		      ft->GetParameter(0),ft->GetParError(0),
		      ft->GetParameter(1),ft->GetParError(1),
		      ft->GetParameter(2),ft->GetParError(2)));
  if (ft->GetNpar()>3)
    tex->DrawLatex(0.20,0.17,
		   Form("p_{3}=%1.3f#pm%1.3f"
			", p_{4}=%1.3f#pm%1.3f",
			ft->GetParameter(3),ft->GetParError(3),
			ft->GetParameter(4),ft->GetParError(4)));

  c1->SaveAs(Form("pdf/%s.pdf",s.c_str()));

  for (int i = 0; i != ft->GetNpar(); ++i) {
    cout << Form("%s%1.4g",i==0 ? "{" : ", ",ft->GetParameter(i));
  }
  cout << "}" << endl;
    

}