C++ TH1F::GetBinError方法代码示例

Double_t calc_asymmetry_syst(TString syst)
{
	TFile *f1 = new TFile("histos/unfolded_syst_"+syst+".root");
	TFile *f3 = new TFile("histos/"+sample+"/efficiency.root");

	TH1F *hunf = (TH1F*)f1->Get("unfolded");
	TH1F *hgen_presel = (TH1F*)f3->Get("hgen_presel");
	//TH1F *hgen_presel_rebin = (TH1F*)f3->Get("hgen_presel_rebin");
	
	TH1F *hasy= (TH1F*)f1->Get("asymmetry");
	//TH1F *hStatErr = (TH1F*)f1->Get("staterr");
	
	// unfolded in bins of generated
	TH1F *hunf_rebin_width = new TH1F(var_y+"_unf",var_y+"_unf",bin_x,list_x);
	for(Int_t i = 1; i <= bin_x; i++) {
		hunf_rebin_width->SetBinContent(i,hunf->GetBinContent(i));
		hunf_rebin_width->SetBinError(i,hunf->GetBinError(i));
	}

	Double_t asy_gen = asymmetry(hgen_presel);
	Double_t asy_unf = hasy->GetMean();

	//cout << endl << asy_gen << " " << asy_unf << endl;

	Double_t diff = asy_gen - asy_unf;
	//Double_t uncertainty = diff/asy_gen;
	Double_t uncertainty = diff;
	cout << uncertainty << endl;
	//cout << hStatErr->GetMean() << endl;

	return uncertainty;
}

TH1F* shape_histos(TH1F* hin, const TString datacard, const TString name)
/*
  use the proper histograms and errors including shpae uncertainties as provided by combine
*/
{
  TH1F* hout = (TH1F*)hin->Clone(); hout->Clear();
  TFile* mlfit = new TFile("fitresults/mlfit.root", "READ");
  TH1F* shape = (TH1F*)mlfit->Get(TString("shapes_fit_s/").Append(datacard).Append("/").Append(name)); // currently problems with data and hioggsprocesses -> different name

  if(shape==0) std::cout << " No histogram found for " << name << std::endl;

  //  for(int i=0; i<hout->GetNbinsX(); ++i)
  for(int i=1; i<hout->GetNbinsX()+1; i++)
  {

    Float_t Norig = hin->GetBinContent(i)*hin->GetBinWidth(i);
    Float_t Nshape = shape->GetBinContent(i);
    Float_t scale_err = (Nshape==0) ? 1 : Norig/Nshape;

    //    hout->SetBinContent(i,shape->GetBinContent(i));
    hout->SetBinContent(i,hin->GetBinContent(i)*hin->GetBinWidth(i));
    hout->SetBinError(i,shape->GetBinError(i)*scale_err);
  }
  mlfit->Close();
  return hout;
}

TH1F * divideHistosForRatio(TH1F * sig, TH1F * bkg) {

  TH1F * cloneS = (TH1F *) sig->Clone( sig->GetTitle() );
  TH1F * cloneB = (TH1F *) bkg->Clone( bkg->GetTitle() );

  /*
  for( Int_t bin = 1; bin <= cloneB->GetNbinsX(); bin++ )
    cloneB->SetBinError( bin, 0. );    
  
  cloneS->Divide( cloneS, cloneB, 1, 1, "" );
  */

  for( Int_t bin = 1; bin <= cloneS->GetNbinsX(); bin++ ){

    Double_t sv = cloneS->GetBinContent(bin);
    Double_t se = cloneS->GetBinError(bin);

    Double_t bv = cloneB->GetBinContent(bin);

    Double_t rat = -1.;
    Double_t err =  0.;

    if( sv != 0. && bv != 0. ) {
      rat = sv / bv;
      // err = sqrt( sv / (bv*bv) );
      err = sqrt( se / (bv*bv) );
    }

    cloneS->SetBinContent( bin, rat );
    cloneS->SetBinError( bin, err );
  }

  return cloneS;
}

//Make a histogram from a TF1
TH1F* makeHist(TF1 *iFit,TH1F *iH,std::string iName) { 
  TH1F *lH = (TH1F*) iH->Clone(iName.c_str());
  for(int i0 = 0; i0 < lH->GetNbinsX()+1; i0++) lH->SetBinContent(i0,iFit->Eval(lH->GetXaxis()->GetBinCenter(i0)));
  for(int i0 = 0; i0 < lH->GetNbinsX()+1; i0++) lH->SetBinContent(i0,lH->GetBinContent(i0)*lH->GetXaxis()->GetBinWidth(i0));
  for(int i0 = 0; i0 < lH->GetNbinsX()+1; i0++) lH->SetBinError (i0,lH->GetBinError   (i0)*lH->GetXaxis()->GetBinWidth(i0));
  return lH;
}

void DYjets_WeightedYield(TTree* tr, TCut cut, float& val, float& err)
{
  TH1F* h = new TH1F("h","h",1,-3.0,3.0);
//  std::cout<<"tr->GetEntries()="<<tr->GetEntries(cut)<<std::endl;
  tr->Draw("phoSCEta>>h",cut,"goff");
  val = h->GetBinContent(1);
  err = h->GetBinError(1);
  delete h;
}// end of DYjets_WeightedYield

TH1F* fixrange(TH1F* old, int numB) {

  float x1, x2;
  string name = old->GetName();

  if (name.find("Ht")!=string::npos) {
    x1 = 30.;
    x2 = 500.;
    if (numB==2) x2 = 400.;
  } else if (name.find("jet_pt")!=string::npos) {
    x1 = 30.;
    x2 = 300.;
    if (numB==2) {
      if (name.find("first")!=string::npos) x2 = 200.;
      if (name.find("second")!=string::npos) x2 = 120.;
    }
  } else if (name.find("pt_Z")!=string::npos) {
    x1 = 0.;
    x2 = 300.;
    if (numB==2) x2 = 230.;
  } else {
    x1 = old->GetXaxis()->GetBinCenter(1);
    x2 = old->GetXaxis()->GetBinCenter(old->GetNbinsX());
  }

  int nx = old->GetXaxis()->FindBin(x2)-old->GetXaxis()->FindBin(x1)+1;

  x1 = old->GetXaxis()->GetBinLowEdge(old->GetXaxis()->FindBin(x1));
  x2 = old->GetXaxis()->GetBinUpEdge(old->GetXaxis()->FindBin(x2));

  TH1F* tmp = new TH1F("tmp",old->GetTitle(),nx,x1,x2);
  tmp->Sumw2();

  tmp->GetXaxis()->SetTitle(old->GetXaxis()->GetTitle());
  tmp->GetYaxis()->SetTitle(old->GetYaxis()->GetTitle());

  for (int i=0;i<=old->GetNbinsX()+1;i++) {
    int ii = tmp->GetXaxis()->FindBin(old->GetXaxis()->GetBinCenter(i));
    float c1 = tmp->GetBinContent(ii);
    float e1 = tmp->GetBinError(ii);
    float c2 = old->GetBinContent(i);
    float e2 = old->GetBinError(i);

    tmp->SetBinContent(ii,c1+c2);
    tmp->SetBinError(ii,TMath::Sqrt(e1*e1+e2*e2));
  }

  tmp->SetEntries(old->GetEntries());

  old->Delete();
  tmp->SetName(name.c_str());

  return tmp;
}

float stat(TString mass){
  TString samplename = "TTbar_Powheg";
  if (mass == "1725") samplename += ".root";
  else samplename += "_mtop" + mass + ".root"; 
  TFile *f = TFile::Open(path + samplename);
  TH1F* h;
  f->GetObject("H_Yields_ElMu", h);
  float y = h->GetBinError(5); 
  delete f;
  return y*Lumi;
}

/**
 * Saving Histogram Data
 */
void ExportHistogram(TFile* f,const char* histKey,const char* outputfile){
  // Root Setup
  TH1F* h = (TH1F*) f->Get(histKey);

  // File Setup
  ofstream out;
  out.open(outputfile);
  out<<"Bin,Value,Error"<<endl;
  for (int i = 1; i < h->GetNbinsX()-1; i++){
    out<<h->GetBinCenter(i)<<","<<h->GetBinContent(i)<<","<<h->GetBinError(i)<<endl;
  }
  out.close();
}

TH1F * getSystErrForRatio(TH1F * bkg) {

  TH1F * cloneB = (TH1F *) bkg->Clone( bkg->GetTitle() );

  for ( Int_t bin = 1; bin <= cloneB->GetNbinsX(); bin++ ) {
    Double_t relError = ( cloneB->GetBinError(bin) / cloneB->GetBinContent(bin) );

    cloneB->SetBinContent( bin, 1. );
    cloneB->SetBinError( bin, relError );
  }

  return cloneB;
}

//*************************************************************************************************
//Normalize Hist
//*************************************************************************************************
TH1F* NormalizeHist(TH1F *originalHist) {
  TH1F* hist = (TH1F*)originalHist->Clone((string(originalHist->GetName())+"_normalized").c_str());
  Double_t norm = 0;
  hist->SetTitle("");
  for (UInt_t b=0; int(b)<hist->GetXaxis()->GetNbins()+2; ++b) {
    norm += hist->GetBinContent(b);
  }
  for (UInt_t b=0; int(b)<hist->GetXaxis()->GetNbins()+2; ++b) {
    hist->SetBinContent(b,hist->GetBinContent(b) / norm);
    hist->SetBinError(b,hist->GetBinError(b) / norm);
  }

  return hist;
}

RooHistPdf* GetRooHistPdf(TFile* fin, string region, string process, string varname, RooRealVar* var, float& norm, TH1F*& histo, bool do_mcstat, bool do_norm, float input_norm){
    TH1F* h  =  GetHisto(fin,region,process,varname, norm, do_norm, input_norm);
    histo = h;
    if(do_mcstat)
    {
        //randomisation of the histo ...
        for(int i=1;i<=h->GetNbinsX();i++)
        {
            h->SetBinContent(i,randomnessGenerator->Gaus(h->GetBinContent(i),h->GetBinError(i)));
        }
    }
    string rdhname = "rdh_"+region+"_"+process;
    RooDataHist *rdh  = new RooDataHist(rdhname.c_str(),rdhname.c_str(),RooArgList(*var),Import(*h));
    string pdfname = "pdf_"+region+"_"+process;
    RooHistPdf *pdf  = new  RooHistPdf(pdfname.c_str(),pdfname.c_str(),RooArgSet(*var),*rdh);
    return pdf;
}

// rebin the spectra
TH1F *rebin_yaxian(TH1F *h, char *histName)
{
  TH1F *hRebin = new TH1F(Form("%s_rebin",h->GetName()),Form("rebin %s",h->GetTitle()),nbins_jetPtBin,boundaries_jetPtBin);
  for (int i=1;i<=h->GetNbinsX();i++)
    {
      double val=h->GetBinContent(i);
      double valErr=h->GetBinError(i);
      int binNum = hRebin->FindBin(h->GetBinCenter(i));
      double val1 = hRebin->GetBinContent(binNum);
      double valErr1 = hRebin->GetBinError(binNum);
      hRebin->SetBinContent(binNum,val+val1);
      hRebin->SetBinError(binNum,sqrt(valErr1*valErr1+valErr*valErr));
    }
  cleanup(hRebin);
  hRebin->SetName(histName);
  return hRebin;
}

//Rebin the histogram
TH1F* rebin(TH1F* iH,int iNBins,double *iAxis) { 
  std::string lTmp = "tmp"; //Added to avoid Root output errors
  TH1F *lH = new TH1F(lTmp.c_str(),lTmp.c_str(),iNBins,iAxis);
  for(int i0 = 0; i0 < iH->GetNbinsX()+1; i0++) {
    int    lNBin = lH->GetXaxis()->FindBin(iH->GetXaxis()->GetBinCenter(i0));
    double lVal  = iH->GetBinContent(i0);
    double lErr  = iH->GetBinError  (i0);
    double lOldV = lH->GetBinContent(lNBin);
    double lOldE = lH->GetBinError  (lNBin);
    lH->SetBinContent(lNBin,lVal+lOldV);
    lH->SetBinError  (lNBin,sqrt(lOldE*lOldE+lErr*lErr));
  }
  std::string lName2 = iH->GetName();
  std::string fine_binning = "_fine_binning";
  lName2.replace(lName2.find(fine_binning),fine_binning.length(),"");
  lH->SetName (lName2.c_str());
  lH->SetTitle(lName2.c_str());
  delete iH;
  return lH;
}

void Fit(TString SIGNAL,TString HISTO,double scaleSGN)
{
  gROOT->ForceStyle();
  TFile *fDat = TFile::Open("Histo_flatTree_data_tmva"+SIGNAL+".root");
  TFile *fBkg = TFile::Open("Histo_flatTree_qcd_weights_tmva"+SIGNAL+".root");
  TFile *fSgn = TFile::Open("Histo_flatTree_"+SIGNAL+"_weights_tmva"+SIGNAL+".root");
  
  TH1 *hDat = (TH1*)fDat->Get(HISTO);
  TH1 *hBkgRaw = (TH1*)fBkg->Get(HISTO);
  TH1 *hSgn = (TH1*)fSgn->Get(HISTO);
  TH1 *hDat_JESlo = (TH1*)fDat->Get(HISTO+"_JESlo");
  TH1 *hBkgRaw_JESlo = (TH1*)fBkg->Get(HISTO+"_JESlo");
  TH1 *hSgn_JESlo = (TH1*)fSgn->Get(HISTO+"_JESlo");
  TH1 *hDat_JESup = (TH1*)fDat->Get(HISTO+"_JESup");
  TH1 *hBkgRaw_JESup = (TH1*)fBkg->Get(HISTO+"_JESup");
  TH1 *hSgn_JESup = (TH1*)fSgn->Get(HISTO+"_JESup");
  
  TH1F *hBkg = (TH1F*)hBkgRaw->Clone("Bkg");
  TH1F *hBkg_JESlo = (TH1F*)hBkgRaw_JESlo->Clone("Bkg_JESlo");
  TH1F *hBkg_JESup = (TH1F*)hBkgRaw_JESup->Clone("Bkg_JESup");
  
  hBkg->Smooth(2);
  hBkg_JESlo->Smooth(2);
  hBkg_JESup->Smooth(2);
  hSgn->Smooth(2);
  hSgn_JESlo->Smooth(2);
  hSgn_JESup->Smooth(2);
  
  double lumi = 4967;
  hBkg->Scale(lumi);
  hBkg_JESlo->Scale(lumi);
  hBkg_JESup->Scale(lumi);
  double k_factor = hDat->Integral()/hBkg->Integral();
  double k_factor_JESlo = hDat->Integral()/hBkg_JESlo->Integral();
  double k_factor_JESup = hDat->Integral()/hBkg_JESup->Integral();
  hBkg->Scale(k_factor);
  cout<<"Signal entries = "<<hSgn->GetEntries()<<endl;
  hSgn->Scale(lumi/scaleSGN);
  hBkg_JESlo->Scale(k_factor_JESlo);
  hSgn_JESlo->Scale(lumi/scaleSGN);
  hBkg_JESup->Scale(k_factor_JESup);
  hSgn_JESup->Scale(lumi/scaleSGN);
  hSgn_JESlo->Scale(hSgn->Integral()/hSgn_JESlo->Integral());
  hSgn_JESup->Scale(hSgn->Integral()/hSgn_JESup->Integral());
  
  TH1 *hBkg_STATlo = (TH1*)hBkg->Clone(HISTO+"_STATlo");
  TH1 *hSgn_STATlo = (TH1*)hSgn->Clone(HISTO+"_STATlo");
  TH1 *hBkg_STATup = (TH1*)hBkg->Clone(HISTO+"_STATup");
  TH1 *hSgn_STATup = (TH1*)hSgn->Clone(HISTO+"_STATup");
  
  float y1,e1;
  for(int i=0;i<hBkg->GetNbinsX();i++) {
    y1 = hBkg->GetBinContent(i+1);
    e1 = hBkg->GetBinError(i+1);
    hBkg_STATlo->SetBinContent(i+1,y1-e1);
    hBkg_STATup->SetBinContent(i+1,y1+e1);
    y1 = hSgn->GetBinContent(i+1);
    e1 = hSgn->GetBinError(i+1);
    hSgn_STATlo->SetBinContent(i+1,y1-e1);
    hSgn_STATup->SetBinContent(i+1,y1+e1);
  }
  hBkg_STATlo->Scale(hBkg->Integral()/hBkg_STATlo->Integral());
  hBkg_STATup->Scale(hBkg->Integral()/hBkg_STATup->Integral());
  hSgn_STATlo->Scale(hSgn->Integral()/hSgn_STATlo->Integral());
  hSgn_STATup->Scale(hSgn->Integral()/hSgn_STATup->Integral());
  
  double xMIN = hBkg->GetBinLowEdge(1);
  double xMAX = hBkg->GetBinLowEdge(hBkg->GetNbinsX()+1);
  double xMIN2 = hDat->GetBinLowEdge(hDat->FindFirstBinAbove(0.5));
  double xMAX2 = hDat->GetBinLowEdge(hDat->FindLastBinAbove(0.5)+1);
  RooRealVar x("x","x",xMIN2,xMAX2);
  
  RooDataHist data("data","dataset with x",x,hDat);
  RooDataHist bkg("qcd","bkg with x",x,hBkg);
  RooDataHist sgn("signal","sgn with x",x,hSgn);
  
  RooHistPdf bkgPDF("bkgPDF","bkgPDF",x,bkg,0);
  RooHistPdf sgnPDF("sgnPDF","sgnPDF",x,sgn,0);
  
  RooRealVar f("f","f",0,0.,1.);
  
  RooAddPdf model("model","model",RooArgList(sgnPDF,bkgPDF),RooArgList(f));
  
  RooFitResult* r = model.fitTo(data,Save());
  r->Print("v");
  double N = hDat->Integral();
  double B = hBkg->Integral();
  double S = hSgn->Integral();
  double m = f.getVal();
  double e = f.getError();
  cout<<"k-factor = "<<k_factor<<endl;
  cout<<N<<" "<<B<<" "<<S<<endl;
  cout<<"Total cross section =       "<<N/lumi<<" pb"<<endl;
  cout<<"Model cross section =       "<<S/lumi<<" pb"<<endl;
  cout<<"Fitted signal strength =    "<<m*N/S<<endl;
  cout<<"Fitted signal error =       "<<e*N/S<<endl;
  double p = 0.95;
  double xup = (N/S)*(m+sqrt(2.)*e*TMath::ErfInverse((1-p)*TMath::Erf(m/e)+p));
  cout<<"Bayesian Upper limit =      "<<xup<<endl;
  RooPlot* frame1 = x.frame();  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
#endif
  }

  /*
    mass plot before and after fit
  */
  TCanvas* canv = MakeCanvas("canv", "histograms", 600, 600);
  canv->cd();
  if(log){ canv->SetLogy(1); }
#if defined MSSM
  if(!log){ data->GetXaxis()->SetRange(0, data->FindBin(350)); } else{ data->GetXaxis()->SetRange(0, data->FindBin(1000)); };
#else
  data->GetXaxis()->SetRange(0, data->FindBin(350));
#endif

  data->SetNdivisions(505);
  data->SetMinimum(min);
  float maxZZ=ZZ->GetBinContent(ZZ->GetMaximumBin()); float maxdata=data->GetBinContent(data->GetMaximumBin());
  if (maxdata>maxZZ)
     data->SetMaximum(1.8*maxdata);
  else
     data->SetMaximum(1.8*maxZZ);
  data->Draw("e");

//  TH1F* errorBand = (TH1F*)ZZ ->Clone();
  TH1F* errorBand = (TH1F*)Zjets ->Clone();
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(1);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  errorBand  ->Scale(0.15);
  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
    if(errorBand->GetBinContent(idx)>0){
      std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
      break;
    }
  }
  //if(log){
    ZZ  ->Draw("histsame");
    Zjets->Draw("histsame");
    $DRAW_ERROR
#ifndef DROP_SIGNAL
    ZH_htt  ->Draw("histsame");
#endif
  //}
  //else{
//#ifndef DROP_SIGNAL
//    ggH  ->Draw("histsame");
//#endif
//    Ztt  ->Draw("histsame");
//    ttbar->Draw("histsame");
//    EWK  ->Draw("histsame");
//    Fakes->Draw("histsame");
//    $DRAW_ERROR
//  }
  data->Draw("esame");
  canv->RedrawAxis();

//  //CMSPrelim(dataset, "#tau_{e}#tau_{#mu}", 0.17, 0.835);  
  CMSPrelim(dataset, "", 0.16, 0.835);  
  TPaveText* chan     = new TPaveText(0.20, 0.74+0.061, 0.32, 0.74+0.161, "NDC");
  chan->SetBorderSize(   0 );
  chan->SetFillStyle(    0 );
  chan->SetTextAlign(   12 );
  chan->SetTextSize ( 0.05 );
  chan->SetTextColor(    1 );