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

void GetFWHM(TF1* floo, double& mean, double& fwhm) {
	TF1* fgt = (TF1*) gROOT->FindObject("fgt");
	if (fgt != 0) delete fgt;
	fgt = new TF1("fgt", gaustail, 0, 6000, 5);
	for (int i=0;i<5;i++)
		fgt->SetParameter(i,floo->GetParameter(i));
	double mean = fgt->GetMaximumX();
	double max = fgt->GetMaximum();
	double left = fgt->GetX(max/2.,0,mean);
	double right = fgt->GetX(max/2.,mean,6000);
	fwhm = right - left;
	cout << "Found " << mean << " " << fwhm << endl;
}

double acceptance2(double accrecoil, double n, 
		   double fbeta=0.27, double frecoil=0.71)
{						
  TF1* f = (TF1*)(gROOT->GetFunction("acc"));
  if(!f)
    f = new TF1("acc","acceptance([0],[1],x)",0,1);
  f->SetParameters(n,frecoil);
  double fcut = f->GetX(accrecoil);
  return acceptance(n,fbeta,fcut);
}

Double_t beta2_for_5sigma(const Int_t fNPts, const Double_t fRangeMin, const Double_t fRangeMax, const Double_t fN_sig_100, const Double_t fN_bkg_100, const Double_t fSigma_N_bkg, const string& fTitle) 
{
 
 Double_t x[fNPts], y[fNPts];

 Double_t step = (fRangeMax - fRangeMin)/(fNPts-1);
   
 for (Int_t i = 0; i < fNPts; i++) {
 
  x[i] = fRangeMin + step*i;
  y[i] = ScP(x[i], fN_sig_100, fN_bkg_100, fSigma_N_bkg);
 }
 
 TCanvas *c_temp = new TCanvas("c_temp","",1120,800);
 c_temp->cd();

 string title = fTitle + ";#beta^{2};S_{cP}";

 TH2F *bg_temp = new TH2F("bg_temp",title.c_str(), 100, fRangeMin, fRangeMax, 100, 0.8*y[0], 1.2*y[fNPts-1]);
 bg_temp->SetStats(kFALSE);
 bg_temp->SetTitleOffset(1.,"X");
 bg_temp->SetTitleOffset(1.,"Y");
 bg_temp->Draw();
 
 TGraph *scP = new TGraph(fNPts, x, y);
 scP->SetMarkerSize(1.);
 scP->SetMarkerStyle(24);
 scP->SetMarkerColor(kRed);
 scP->Draw("P"); 
 scP->Fit("pol2");

 TF1 *fit = (TF1*)scP->GetFunction("pol2");
 Double_t beta2 = fit->GetX(5);
 Double_t N_sig = beta2*fN_sig_100;
 Double_t N_s_b = N_sig + fN_bkg_100;
 
 cout<<">> beta2 for 5 sigma discovery for "<<fTitle<<" = "<<beta2<<"\n";
 cout<<">> ** N_sig = "<<N_sig<<"\n";
 cout<<">> ** N_bkg = "<<fN_bkg_100<<"\n";
 cout<<">> ** N_s_b = "<<N_s_b<<"\n";
 
 if(fSigma_N_bkg != 0) c_temp->SaveAs((fTitle + "_significance_beta_sys.png").c_str());
 else c_temp->SaveAs((fTitle + "_significance_beta.png").c_str());

 delete scP;
 delete bg_temp;
 delete c_temp;
 
 return beta2;
}

//.........这里部分代码省略.........
	
	TTree* ic5t = (TTree*)inf->Get("icPu5JetAnalyzer/t");
	TTree* t = (TTree*)inf->Get("hltanalysis/HltTree");
	
	ic5t->SetBranchAddress("jtpt",jtpt_ic5);
	ic5t->SetBranchAddress("jteta",jteta_ic5);
	ic5t->SetBranchAddress("jtphi",jtphi_ic5);
	ic5t->SetBranchAddress("nref",&njets_ic5);
	
	TTree* ak3t = (TTree*)inf->Get("akPu3PFJetAnalyzer/t");
	
	ak3t->SetBranchAddress("jtpt",jtpt_ak3);
	ak3t->SetBranchAddress("jteta",jteta_ak3);
	ak3t->SetBranchAddress("jtphi",jtphi_ak3);
	ak3t->SetBranchAddress("nref",&njets_ak3);
	
	t->SetBranchAddress("Run",&run);
	t->SetBranchAddress("Event",&evt);
	t->SetBranchAddress("hiBin",&bin);
	
	TString outname = "JetSkimTree.root";
			
	int Nevents = t->GetEntries();
	for(int iev = 0; iev < Nevents; ++iev){
		t->GetEntry(iev);
		ic5t->GetEntry(iev);
		ak3t->GetEntry(iev);
		
		for(int i = 0; i < njets_ak3; ++i){
			//if(jtpt_ak3[i] < cutPtJet) continue;
			if(fabs(jteta_ak3[i]) > cutEtaJet) continue;
			hak3pujet->Fill(jtpt_ak3[i]);
							}//1st for loop
		
		int id=0;
		for(int i = 0; i < njets_ak3; ++i){
			if(fabs(jteta_ic5[i]) < cutEtaJet)
				id =1;
			
		}// 2nd for loop
		if(id == 1){
				for(int i = 0; i < njets_ak3; ++i){
					if(jtpt_ic5[i] < cutPtJet) continue;
					//if(jtpt_ak3[i] < cutPtJet) continue;
					if(fabs(jteta_ak3[i]) > cutEtaJet) continue;
					hak3pujet_cut->Fill(jtpt_ak3[i]);
				}
			
	}// if loop
	}//evt loop
	
	
	TCanvas *c1 = new TCanvas("c1", "",46,477,700,509);
	c1->Range(-125,-11.1492,1125,3.24528);
	c1->SetBorderSize(0);
	c1->SetBorderMode(0);
	c1->SetLeftMargin(0.16);
	c1->SetTopMargin(0.06);
	c1->SetBottomMargin(0.17);
	c1->SetTicks(1,1);
	
	
	
	gROOT->LoadMacro("erf.C");
		
	TGraphAsymmErrors *dividedthingy= new TGraphAsymmErrors;
	dividedthingy->BayesDivide(hak3pujet_cut,hak3pujet);
	
	
	float xstart=20.0;
	float xend=300.0;
	
	
	TF1 *fitfcn = new TF1("fitfcn",erf,xstart,xend,2); //<============
	fitfcn->SetParameter(0,40); //<============= 
	TFitResultPtr results =dividedthingy->Fit("fitfcn","VSR");
	
	cout << "Chi2: " << results->Chi2() << " / " << results->Ndf() << " NDF" << endl;
	
	float turnon=fitfcn->GetX(0.99,10,300);
	cout << "99% at: " << turnon << endl;
	char turnontext[50]="";
	sprintf(turnontext," ak5pu at %4.1f",turnon);
	
	
	
	dividedthingy->GetFunction("fitfcn")->SetLineColor(3);
	dividedthingy->SetLineColor(1);
	dividedthingy->Draw("Ap");
	
	
	
	TText *t2 = new TText(95,0.08,turnontext);
	t2->SetTextSize(0.055);
	t2->SetTextColor(1);
	t2->Draw();
		
		
	
}

//.........这里部分代码省略.........
    hnlo_cteq->SetBinError(jk, 0.);
    hnlo_mstw->SetBinError(jm, 0.);
    hnlo_nnpdf->SetBinError(jn, 0.);
    hnlo_hera->SetBinError(jh, 0.);
  }
  */

  // initial fit of the NLO curve
  TF1 *fnlo = new TF1("fnlo",
		      "[0]*exp([1]/x)*pow(x,[2])"
		      "*pow(1-x*cosh([4])/3500.,[3])", 10., 1000.);
  fnlo->SetParameters(2e14,-18,-5.2,8.9,0.5*ieta);
  fnlo->FixParameter(4,0.5*ieta);

  hnlo = hmc;
  hnlo->Fit(fnlo,"QRN");
  
  // Graph of theory points with centered bins
  const double minerr = 0.02;
  TGraphErrors *gnlo = new TGraphErrors(0);
  TGraphErrors *gnlo2 = new TGraphErrors(0); // above + minerr
  TGraphErrors *gnlocut = new TGraphErrors(0); // only up to expected pT
  gnlo->SetName("gnlo");
  gnlo2->SetName("gnlo2");
  gnlocut->SetName("gnlocut");
  for (int i = 1; i != hnlo->GetNbinsX()+1; ++i) {

    double y = hnlo->GetBinContent(i);
    double dy = hnlo->GetBinError(i);
    double ptmin = hnlo->GetBinLowEdge(i);
    double ptmax = hnlo->GetBinLowEdge(i+1);

    double y0 = fnlo->Integral(ptmin, ptmax) / (ptmax - ptmin);
    double x = fnlo->GetX(y0, ptmin, ptmax);

    int n = gnlo->GetN();
    tools::SetPoint(gnlo, n, x, y, 0, dy);
    tools::SetPoint(gnlo2, n, x, y, 0, tools::oplus(dy, minerr*y));
    if (y*(ptmax-ptmin) > 0.0001) { // for 10/fb
      int m = gnlocut->GetN();
      tools::SetPoint(gnlocut, m, x, y, 0, 0);
    }
  }

  gnlo2->Fit(fnlo,"QRN");

  // Divide graph with fit to check stability
  TGraphErrors *gnlofit = new TGraphErrors(0);
  gnlofit->SetName("gnlofit");
  for (int i = 0; i != gnlo->GetN(); ++i) {
    
    double x, y, ex, ey;
    tools::GetPoint(gnlo, i, x, y, ex, ey);
    double f = fnlo->Eval(x);
    tools::SetPoint(gnlofit, i, x, y/f, ex, ey/f);
  }

  // Rebin theory to match data bins
  TH1D *hnlo0 = hnlo;
  hnlo0->SetName("hnlo0");
  hnlo = tools::Rebin(hnlo0, hpt);
  hnlo->SetName("hnlo");

  // Same for earlier 2010 and 2011 studies to get a matching theory set
  // (moved from PDF4LHC to CT10 since 2010)
  if (false) {

void fitTF1(TCanvas *canvas, TH1F h, double XMIN_, double XMAX_, double dX_, double params[], Color_t LC_=kBlack) { //double& FWHM_, double& x0_, double& x1_, double& x2_, double& y0_, double& y1_, double& INT_, double& YIELD_) {
//TCanvas* fitTF1(TH1F h, double HSCALE_, double XMIN_, double XMAX_) {
//TF1* fitTF1(TH1F h, double HSCALE_, double XMIN_, double XMAX_) {
	gROOT->ForceStyle();
	RooMsgService::instance().setSilentMode(kTRUE);
	for(int i=0;i<2;i++) RooMsgService::instance().setStreamStatus(i,kFALSE);

	//TCanvas *canvas = new TCanvas(TString::Format("canvas_%s",h.GetName()),TString::Format("canvas_%s",h.GetName()),1800,1200);

	RooDataHist *RooHistFit = 0;
	RooAddPdf *model = 0;

	RooWorkspace w = RooWorkspace("w","workspace");

	RooRealVar x("mbbReg","mbbReg",XMIN_,XMAX_);
	RooRealVar kJES("CMS_scale_j","CMS_scale_j",1,0.9,1.1);
	RooRealVar kJER("CMS_res_j","CMS_res_j",1,0.8,1.2);
	kJES.setConstant(kTRUE);
	kJER.setConstant(kTRUE);

	TString hname = h.GetName();

	RooHistFit = new RooDataHist("fit_"+hname,"fit_"+hname,x,&h);


	RooRealVar YieldVBF = RooRealVar("yield_"+hname,"yield_"+hname,h.Integral());
	RooRealVar m("mean_"+hname,"mean_"+hname,125,100,150);
	RooRealVar s("sigma_"+hname,"sigma_"+hname,12,3,30);
	RooFormulaVar mShift("mShift_"+hname,"@0*@1",RooArgList(m,kJES));
	RooFormulaVar sShift("sShift_"+hname,"@0*@1",RooArgList(m,kJER));
	RooRealVar a("alpha_"+hname,"alpha_"+hname,1,-10,10);
	RooRealVar n("exp_"+hname,"exp_"+hname,1,0,100);
	RooRealVar b0("b0_"+hname,"b0_"+hname,0.5,0.,1.);
	RooRealVar b1("b1_"+hname,"b1_"+hname,0.5,0.,1.);
	RooRealVar b2("b2_"+hname,"b2_"+hname,0.5,0.,1.);
	RooRealVar b3("b3_"+hname,"b3_"+hname,0.5,0.,1.);
	
	RooBernstein bkg("signal_bkg_"+hname,"signal_bkg_"+hname,x,RooArgSet(b0,b1,b2));
	RooRealVar fsig("fsig_"+hname,"fsig_"+hname,0.7,0.0,1.0);
	RooCBShape sig("signal_gauss_"+hname,"signal_gauss_"+hname,x,mShift,sShift,a,n);

	model = new RooAddPdf("signal_model_"+hname,"signal_model_"+hname,RooArgList(sig,bkg),fsig);

	//RooFitResult *res = model->fitTo(*RooHistFit,RooFit::Save(),RooFit::SumW2Error(kFALSE),"q");
	model->fitTo(*RooHistFit,RooFit::Save(),RooFit::SumW2Error(kFALSE),"q");

	//res->Print();
	//model->Print();
	
	canvas->cd();
	canvas->SetTopMargin(0.1);
	RooPlot *frame = x.frame();
// no scale
	RooHistFit->plotOn(frame);
	model->plotOn(frame,RooFit::LineColor(LC_),RooFit::LineWidth(2));//,RooFit::LineStyle(kDotted));
	model->plotOn(frame,RooFit::Components(bkg),RooFit::LineColor(LC_),RooFit::LineWidth(2),RooFit::LineStyle(kDashed));
// with scale
//	RooHistFit->plotOn(frame,RooFit::Normalization(HSCALE_,RooAbsReal::NumEvent));
//	model->plotOn(frame,RooFit::Normalization(HSCALE_,RooAbsReal::NumEvent),RooFit::LineWidth(1));
//	model->plotOn(frame,RooFit::Components(bkg),RooFit::LineColor(kBlue),RooFit::LineWidth(1),RooFit::LineStyle(kDashed),RooFit::Normalization(HSCALE_,RooAbsReal::NumEvent));
	 
	frame->GetXaxis()->SetLimits(50,200);
	frame->GetXaxis()->SetNdivisions(505);
	frame->GetXaxis()->SetTitle("M_{b#bar{b}} (GeV)");
	frame->GetYaxis()->SetTitle("Events");
	frame->Draw();
	h.SetFillColor(kGray);
	h.Draw("hist,same");
	frame->Draw("same");
	gPad->RedrawAxis();

	TF1 *tmp = model->asTF(x,fsig,x);
	//tmp->Print();

	double y0_ = tmp->GetMaximum();
	double x0_ = tmp->GetMaximumX();
	double x1_ = tmp->GetX(y0_/2.,XMIN_,x0_);
	double x2_ = tmp->GetX(y0_/2.,x0_,XMAX_);
	double FWHM_ = x2_-x1_;
	double INT_ = tmp->Integral(XMIN_,XMAX_);
	double YIELD_= YieldVBF.getVal();
	double y1_ = dX_*0.5*y0_*(YieldVBF.getVal()/tmp->Integral(XMIN_,XMAX_));

	params[0] = x0_;
	params[1] = x1_;
	params[2] = x2_;
	params[3] = y0_;
	params[4] = y1_;
	params[5] = FWHM_;
	params[6] = INT_;
	params[7] = YIELD_;

	//cout<<"Int = "<<tmp->Integral(XMIN_,XMAX_)<<", Yield = "<<YieldVBF.getVal()<<", y0 = "<<y0_<<", y1 = "<<y1_ <<", x0 = "<< x0_ << ", x1 = "<<x1_<<", x2 = "<<x2_<<", FWHM = "<<FWHM_<<endl;

	TLine ln = TLine(x1_,y1_,x2_,y1_);
	ln.SetLineColor(kMagenta+3);
	ln.SetLineStyle(7);
	ln.SetLineWidth(2);
	ln.Draw();
	
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  grapherrors_myStyle(gr,"gr",2,30,1,30,1001,-1111,-1111,510,510,21,36,1.,"M_{t} (GeV)","M_{J/#psi+l} (GeV)");
  gr->Draw("AP");
  cn_calib->Update();
  fit->SetLineColor(30);
  fit_up->SetLineColor(29); fit_up->SetLineStyle(2); fit_up->SetLineWidth(2);
  fit_up->Draw("same");
  fit_down->SetLineColor(29); fit_down->SetLineStyle(2); fit_down->SetLineWidth(2);
  fit_down->Draw("same");
  leg_myStyle(leg_calib);
  leg_calib->SetTextAlign(22);
  leg_calib->Draw();  
  channel_tex->Draw("same");
  cms_myStyle(lumi, false);
  TString outBlind = outdir;
  if (binned < 2)
    outBlind += "BlindedUnbinnedCalibration";
  if (binned == 2)
    outBlind += "BlindedBinnedCalibration";
  cn_calib->SaveAs(outBlind+".pdf");
  cn_calib->SaveAs(outBlind+".C");
  cn_calib->SaveAs(outBlind+".jpg");
  cn_calib->SaveAs(outBlind+".eps");

  // Unblinded analysis
  double *data_fit = new double[2];
  if (binned == 0)
    data_fit = unbinnedFit(indir+fileData, xlim, 0., channel_tex, outdir, lumi);
  if (binned > 0)
    data_fit = binnedFit(indir+fileData, xlim, 0., channel_tex, outdir, lumi);
  double mean_gaus = data_fit[0];
  double err_gaus = data_fit[1];
  TString mjpsil_res = TString::Format("M_{J/#psi+l} = (%3.1f #pm %3.1f) GeV", mean_gaus, err_gaus);

  double x_horiz[4]={gr->GetXaxis()->GetXmin(),fit->GetX(mean_gaus+err_gaus,0,250),fit->GetX(mean_gaus-err_gaus,0,250),gr->GetXaxis()->GetXmin()};
  double y_horiz[4]={mean_gaus+err_gaus,mean_gaus+err_gaus,mean_gaus-err_gaus,mean_gaus-err_gaus};
  TPolyLine *horiz = new TPolyLine(4,x_horiz,y_horiz);
  TLine *horiz_sup = new TLine(x_horiz[0], y_horiz[0], x_horiz[1], y_horiz[1]);
  TLine *horiz_med = new TLine(gr->GetXaxis()->GetXmin(), mean_gaus, fit->GetX(mean_gaus,0,250), mean_gaus);
  TLine *horiz_inf = new TLine(x_horiz[2], y_horiz[2], x_horiz[3], y_horiz[3]);
  double x_vert_min[4]={fit_up->GetX(mean_gaus-err_gaus,0,250),fit_up->GetX(mean_gaus-err_gaus,0,250),x_horiz[2],x_horiz[2]};
  double y_vert_min[4]={gr->GetYaxis()->GetXmin(),mean_gaus-err_gaus,mean_gaus-err_gaus,gr->GetYaxis()->GetXmin()};
  TPolyLine *vert_min = new TPolyLine(4,x_vert_min,y_vert_min);
  double x_vert_mean[4]={x_horiz[2],x_horiz[2],x_horiz[1],x_horiz[1]};
  double y_vert_mean[4]={gr->GetYaxis()->GetXmin(),y_horiz[2],y_horiz[1],gr->GetYaxis()->GetXmin()};
  TPolyLine *vert_mean = new TPolyLine(4,x_vert_mean,y_vert_mean);
  double x_vert_max[4]={x_horiz[1],x_horiz[1],fit_down->GetX(mean_gaus+err_gaus,0,250),fit_down->GetX(mean_gaus+err_gaus,0,250)};
  double y_vert_max[4]={gr->GetYaxis()->GetXmin(),mean_gaus+err_gaus,mean_gaus+err_gaus,gr->GetYaxis()->GetXmin()};
  TPolyLine *vert_max = new TPolyLine(4,x_vert_max,y_vert_max);
  TLine *vert1 = new TLine(x_vert_min[0],y_vert_min[0],x_vert_min[1],y_vert_min[1]);
  TLine *vert2 = new TLine(x_vert_mean[0],y_vert_mean[0],x_vert_mean[1],y_vert_mean[1]);
  TLine *vert3 = new TLine(fit->GetX(mean_gaus,0,250),y_vert_mean[0],fit->GetX(mean_gaus,0,250),mean_gaus);
  TLine *vert4 = new TLine(x_vert_mean[3],y_vert_mean[3],x_vert_mean[2],y_vert_mean[2]);
  TLine *vert5 = new TLine(x_vert_max[3],y_vert_max[3],x_vert_max[2],y_vert_max[2]);
  TString mt_res; 
  double pm_fit = (x_vert_mean[3]-x_vert_mean[0])/2.;
  double p_fit = x_vert_max[3]-x_vert_mean[3];
  double m_fit = x_vert_mean[0]-x_vert_min[0];
  if (fabs(p_fit-m_fit)<0.1) 
    mt_res= TString::Format("M_{t} = (%0.1f #pm %0.1f #pm %0.1f) GeV", fit->GetX(mean_gaus,0,250), pm_fit, fabs(m_fit));
  else 
    mt_res= TString::Format("M_{t} = (%0.1f #pm %0.1f #splitline{_{+ %0.1f}}{_{- %0.1f}}) GeV", fit->GetX(mean_gaus,0,250),pm_fit, fabs(p_fit), fabs(m_fit));
  double *mtop_res = new double[2];
  mtop_res[0] = fit->GetX(mean_gaus,0,250);
  mtop_res[1] = pm_fit + max(fabs(p_fit), fabs(m_fit));

  TLegend *leg_res = new TLegend(0.2,0.75,0.65,0.87,NULL,"brNDC");