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

void performClosure(RooRealVar *mass, RooAbsPdf *pdf, RooDataSet *data, string closurename, double wmin=110., double wmax=130., double slow=110., double shigh=130., double step=0.002) {
  
  // plot to perform closure test
  cout << "Performing closure test..." << endl; 
  double nbins = (wmax-wmin)/step;
  TH1F *h = new TH1F("h","h",int(floor(nbins+0.5)),wmin,wmax);
  if (data){
    pdf->fillHistogram(h,RooArgList(*mass),data->sumEntries());
    h->Scale(2*h->GetNbinsX()/double(binning_));
  }
  else {
    pdf->fillHistogram(h,RooArgList(*mass));
  }
  int binLow = h->FindBin(slow);
  int binHigh = h->FindBin(shigh)-1;
  TH1F *copy = new TH1F("copy","c",binHigh-binLow,h->GetBinLowEdge(binLow),h->GetBinLowEdge(binHigh+1));
  for (int b=0; b<copy->GetNbinsX(); b++) copy->SetBinContent(b+1,h->GetBinContent(b+1+binLow));
  double areaCov = 100*h->Integral(binLow,binHigh)/h->Integral();
 
  // style
  h->SetLineColor(kBlue);
  h->SetLineWidth(3);
  h->SetLineStyle(7);
  copy->SetLineWidth(3);
  copy->SetFillColor(kGray);
  
  TCanvas *c = new TCanvas();
  if (data){
    RooPlot *plot = mass->frame(Bins(binning_),Range("higgsRange"));
    plot->addTH1(h,"hist");
    plot->addTH1(copy,"same f");
    if (data) data->plotOn(plot);
    pdf->plotOn(plot,Normalization(h->Integral(),RooAbsReal::NumEvent),NormRange("higgsRange"),Range("higgsRange"),LineWidth(1),LineColor(kRed),LineStyle(kDashed));
    plot->Draw();
    c->Print(closurename.c_str());
  }
  else {
    RooPlot *plot = mass->frame(Bins(binning_),Range("higgsRange"));
    h->Scale(plot->getFitRangeBinW()/h->GetBinWidth(1));
    copy->Scale(plot->getFitRangeBinW()/h->GetBinWidth(1));
    pdf->plotOn(plot,LineColor(kRed),LineWidth(3));
    plot->Draw();
    h->Draw("hist same");
    copy->Draw("same f");
    c->Print(closurename.c_str());
  }
  cout << "IntH: [" << h->GetBinLowEdge(binLow) << "-" << h->GetBinLowEdge(binHigh+1) << "] Area = " << areaCov << endl;
  delete c;
  delete copy;
  delete h;
}

  TH1F * getpdgs(std::string const & dr, std::string const & we) {
    
    TCanvas * canvas_temp = new TCanvas("temp");
    tree->Draw((dr+">>h").c_str(), we.c_str());
    TH1F * hist = (TH1F*)gDirectory->Get("h");
    int const min = hist->GetBinLowEdge(hist->FindFirstBinAbove(0));
    int const max = hist->GetBinLowEdge(hist->FindLastBinAbove(0)) +
      hist->GetBinWidth(hist->FindLastBinAbove(0)) + 1;
    int const binno = max - min;
    delete canvas_temp;
    delete hist;

    canvas_temp = new TCanvas("temp");
    tree->Draw(Form((dr+">>h(%d,%d,%d)").c_str(), binno, min, max), we.c_str());
    ofile->cd();
    delete canvas_temp; 
    TH1F * histb = (TH1F*)gDirectory->Get("h");

    if(histb->GetEntries() == 0) {
      delete histb;
      return nullptr;
    }    

    return histb;
    
  }

// get FWHHM
vector<double> getFWHM(RooRealVar *mass, RooAbsPdf *pdf, RooDataSet *data, double wmin=110., double wmax=130., double step=0.0004) {
 
  cout << "Computing FWHM...." << endl;
  double nbins = (wmax-wmin)/step;
  TH1F *h = new TH1F("h","h",int(floor(nbins+0.5)),wmin,wmax);
  if (data){
    pdf->fillHistogram(h,RooArgList(*mass),data->sumEntries());
  }
  else {
    pdf->fillHistogram(h,RooArgList(*mass));
  }
  
  double hm = h->GetMaximum()*0.5;
  double low = h->GetBinCenter(h->FindFirstBinAbove(hm));
  double high = h->GetBinCenter(h->FindLastBinAbove(hm));

  cout << "FWHM: [" << low << "-" << high << "] Max = " << hm << endl;
  vector<double> result;
  result.push_back(low);
  result.push_back(high);
  result.push_back(hm);
  result.push_back(h->GetBinWidth(1));
  
  delete h;
  return result;
}

// Called within makePlot when making the pu plot
// Basically adds in the +/- 5% systematic error bars
TGraphAsymmErrors makePUPlot( THStack stack, TString anaType, const float lumi ) {

  // Got the stack, now sum the =/- 5% histograms and make two new stacks
  // Make dummy legend etc.
  THStack p5Stack("Background MC p5","");
  THStack m5Stack("Background MC m5","");
  TLegend *legend= new TLegend(0.1,0.1,0.1,0.1);

  addBackgroundHistos( anaType, "nRecoPV_p5", p5Stack, legend, lumi );
  addBackgroundHistos( anaType, "nRecoPV_m5", m5Stack, legend, lumi );

  // Get last histogram in stacks i.e. sum of all backgrounds
  TH1F * hist = new TH1F( *(TH1F*)(stack.GetStack()->Last()) );
  TH1F * p5Hist = new TH1F( *(TH1F*)p5Stack.GetStack()->Last() );
  TH1F * m5Hist = new TH1F( *(TH1F*)m5Stack.GetStack()->Last() );

  // Create structures to make TGraphAssymErrors
  Int_t n=p5Hist->GetNbinsX();
  Double_t x[n]; // x values
  Double_t y[n]; // y values
  Double_t exl[n]; // low x error
  Double_t eyl[n]; // low y error
  Double_t exh[n]; // high x error
  Double_t eyh[n]; // high y error

  // Loop over histos and get values
  for ( int bin=1; bin < n+1; bin++ ) {
    x[bin-1]=hist->GetBinCenter( bin );
    y[bin-1]=hist->GetBinContent( bin );
    exl[bin-1]=hist->GetBinWidth( bin )/2;
    eyl[bin-1]=p5Hist->GetBinContent( bin )-y[bin-1];
    exh[bin-1]=hist->GetBinWidth( bin )/2;
    eyh[bin-1]=y[bin-1]-m5Hist->GetBinContent( bin );
  }


  // Make TGraphAssymErrors
  TGraphAsymmErrors gr(n,x,y,exl,exh,eyl,eyh);

  return gr;
}

void FitterCBForSignificance(){



  //  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_gg_TuneD6T_Emine2013.root");

  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_qg_TuneD6T_Emine2013.root");



  //  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_Qstar_2012_D6T_ak5_fat30_save.root");
  //  TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_RSGraviton_2012_D6T_ak5_GGtoGG_fat30_save.root");
  //TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_RSGraviton_2012_D6T_ak5_QQtoQQ_fat30_save.root");

  //TFile* _shapes = new TFile("/afs/cern.ch/user/m/mgouzevi/scratch0/CMGTools/CMSSW_4_2_8/src/StatTools/BayesianDijetFit/Results/Resonance_Shapes_Qstar_ak5_fat30.root");
  for (int i = 1; i < 2; i++){
	   TH1F* shape = (TH1F*) _shapes->Get(Form("h_qg_%d;1", 1000+i*100));
	   TH1F* shape_largeBinning = new TH1F("shape_large_binning", "", NBINS, BOUNDARIES);

	   for (int j = 1; j < shape->GetNbinsX()+1; j++){
	     double binCenter = shape->GetBinCenter(j);
	     double weight = shape->GetBinContent(j);
	     double binWidth = shape->GetBinWidth(j);
	     shape_largeBinning->Fill(binCenter, weight);
	     //	     cout << "j = " << j << " binWidth = " << binWidth << " weight = " << weight << endl;
	     shape->SetBinContent(j, weight/binWidth);
	   }

	   for (int j = 1; j < shape_largeBinning->GetNbinsX()+1; j++){
	     double weight = shape_largeBinning->GetBinContent(j);
	     double binWidth = shape_largeBinning->GetBinWidth(j);
	     shape_largeBinning->SetBinContent(j, weight/binWidth*1000000);
	   }

	   shape->Draw();

	   shape_largeBinning->SetLineColor(kRed);
	   shape_largeBinning->Draw("SAME");

	   //   
	   FitHistWithCBShape(shape_largeBinning, 1000.+i*100);		
	}

	  cout << "mass\tsMean\tsSigma\tsAlphaHigh\tsAlphaLow\tsNHigh\tsNLow\tsFrac"<<endl;

	for (int i = 0; i < 40; i++){
	
	  cout << Masses[i] << "\t" << Means[i] << "\t" << Sigmas[i] << "\t"
	       << alphaHights[i] << "\t\t" << alphaLows[i] << "\t"
	       << nHighs[i] << "\t\t" << nLows[i] << "\t" << fracs[i] << endl;

	}
}

void PIDPedes()
{
  TH1F* Energy;
  Char_t Buff[256];

  for(Int_t ch=0; ch<24; ch++)
  {
    sprintf(Buff, "PID_Energy%d", ch);
    Energy = (TH1F*)gROOT->FindObject(Buff);
    printf("%4.0f\n", Energy->GetMaximumBin()*Energy->GetBinWidth(Energy->GetMaximumBin()));
  }
}

void PIDTimes()
{
  TH1F* Time;
  Char_t Buff[256];

  for(Int_t ch=0; ch<24; ch++)
  {
    sprintf(Buff, "PID_Time%d", ch);
    Time = (TH1F*)gROOT->FindObject(Buff);
    printf("%5.1f\n", Time->GetMaximumBin()*Time->GetBinWidth(Time->GetMaximumBin()) - 200.0);
  }
}

void VetoTimes()
{
  TH1F* Time;
  Char_t Buff[256];

  for(Int_t ch=0; ch<384; ch++)
  {
    sprintf(Buff, "Veto_Time%d", ch);
    Time = (TH1F*)gROOT->FindObject(Buff);
    printf("%6.1f\n", (Time->GetMaximumBin()*Time->GetBinWidth(Time->GetMaximumBin()) - 500.0)/0.05000);
  }
}

void LadderTimes()
{
  TH1F* Time;
  Double_t par[4];
  Double_t Max;
  Char_t Buff[256];

  for(Int_t ch=0; ch<352; ch++)
  {
    sprintf(Buff, "Ladder_Time%d", ch);
    Time = (TH1F*)gROOT->FindObject(Buff);
    Max = Time->GetMaximumBin()*Time->GetBinWidth(Time->GetMaximumBin()) - 1000.0;
   
    TF1* gauss = new TF1("gauss", "gaus", Max-8, Max+8);
    Time->Fit(gauss, "RQ+");
    gauss->GetParameters(&par[0]);
    
    printf("%5.2f\n", par[1]/0.117710);
  }
}

void countZs(const char* filename, int run1=1, int run2=0, float lumi=0, bool add=false) {
   TFile *file = new TFile(filename);
   if (!file) return;
   
   TCanvas *c1 = new TCanvas();
   c1->cd();

   TH1F *h = (TH1F*) file->Get("Run summary/DiMuonHistograms/GlbGlbMuon_HM");
   h->GetYaxis()->SetRangeUser(0,200);
   h->Draw();

   TF1 *f = new TF1("f","gaus(0)+pol2(3)",70,110);
   f->SetParName(0,"norm");
   f->SetParName(1,"mass");
   f->SetParName(2,"width");
   f->SetParameters(1.16028e+02,9.07785e+01,2.11556e+00,-2.17600e+01,9.67443e-01,-6.98382e-03);
   f->SetParLimits(0,0,1e6);
   f->SetParLimits(1,85,95);
   f->SetParLimits(2,1.5,5);
   h->Fit(f,"","",70,110);
   h->GetYaxis()->SetRangeUser(0,200);
   c1->Update();

   TF1 *f2 = new TF1("f2","gaus(0)",70,110);
   f2->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2));
   int integral = (int) (f2->Integral(70,110) / h->GetBinWidth(2));

   TLatex *tl = new TLatex();
   tl->DrawLatexNDC(0.1,0.95,Form("pPb 8TeV [%i-%i, %.1f nb^{-1}]",run1,run2,lumi));
   tl->DrawLatexNDC(0.18,0.75,Form("%i Z bosons",integral));
   tl->DrawLatexNDC(0.18,0.7,Form("(#sigma = %.1f nb)",integral/lumi));

   if (!add) c1->SaveAs(Form("plotZ_%i_%i.pdf",run1,run2));
   else c1->SaveAs("plotZ.gif+100");

   file->Close();
}

//............................................
// Constructor for accpetance from a ROOT Tfile
SlicedAcceptance::SlicedAcceptance( string type, string fileName,string histName, bool fluctuate, bool quiet ) :
	slices(), nullSlice(new AcceptanceSlice(0.,0.,0.)), tlow(), thigh(), beta(), _sortedSlices(false), maxminset(false), t_min(0.), t_max(0.), _hasChecked(false), _storedDecision(false)
{
	if(!quiet) cout << "Root file being used for acceptance" << endl;
	(void)type;
	if( type != "RootFile" ) {   }//do nothing for now

	string fullFileName = StringProcessing::FindFileName( fileName, quiet );

	if( !quiet ) cout << "Opening: " << fullFileName << endl;

	TFile* file = TFile::Open(TString(fullFileName));
 	if(!quiet) cout << "File " << fullFileName << " opened!" << endl;
	TH1F* histo = (TH1F*)file->Get(TString(histName));
	if(!quiet) cout << "Histo " << histName << " opened!" << endl;
	histo->Draw();
	if(!quiet) cout << "Histo " << histName << " drawn!" << endl;
//	histo->Sumw2();

	if(fluctuate){
	cout << "WARNING! You have fluctuated the acceptance." << endl;
	cout << "WARNING! This is for systematic studies only. " << endl;
	cout << "WARNING! Projections and pull fits will have a different fluctuated acceptance to the PDF you fit with." << endl;
	cout << "WARNING! ONLY USE FluctuateAcceptance:True FOR SYSTEMATIC STUDIES" << endl;
	TRandom3 * rng = new TRandom3(0);	
	//Randomly fluctuate bin contents within error: 
	for (int l = 1; l <= histo->GetNbinsX(); ++l){
	if(!quiet)	cout << "Bin content and error before: " << histo->GetBinContent(l) << "+/-" << histo->GetBinError(l);
		histo->SetBinContent(l,rng->Gaus(histo->GetBinContent(l),histo->GetBinError(l)));
	if(!quiet)	cout <<  " and after: " << histo->GetBinContent(l) << "+/-" << histo->GetBinError(l) << endl;
	}
	delete rng;
	}

	histo->Scale(1./(histo->GetBinContent(histo->GetMaximumBin())));
	histo->SetMinimum(0);


	double maxend = histo->GetBinLowEdge(histo->GetNbinsX()) + histo->GetBinWidth(histo->GetNbinsX());
	double height;
	double start;
	double end  = histo->GetBinLowEdge(histo->GetNbinsX()) + histo->GetBinWidth(histo->GetNbinsX());
	double dheight;
	for (int l = 1; l <= histo->GetNbinsX(); ++l){
		height = histo->GetBinContent(l);
		dheight = height;
		for (int n = l; n>0; n--){
			if(histo->GetBinContent(n)<height){
				dheight = height - histo->GetBinContent(n);
				cout << l << " " << n << " " << dheight << endl;
				break;
			}
		}
		start = histo->GetBinLowEdge(l);
		end = maxend;
		for (int m = l; m <= histo->GetNbinsX(); ++m){
			double thisbinheight = histo->GetBinContent(m);
			if(thisbinheight<height){
				end = histo->GetBinLowEdge(m);
				break;
			}
		}
		slices.push_back( new AcceptanceSlice( start, end, dheight ) );
	if(!quiet)	cout << start << "	" << end << "      " << dheight << endl;
	}
	histo->Delete();
//	delete histo;
	file->Close();
	delete file;

	if( !quiet ) cout << "Time Acc Slices: " << slices.size() << endl;
	if( slices.size() == 1 )
	{
		cout << "SlicedAcceptance: SERIOUS ERROR" << endl;
		exit(0);
	}
	//....done.....


	_sortedSlices = this->isSorted();

	if( _sortedSlices )
	{
		if( !quiet ) cout << "Sliced Acceptance is using sorted horizontal slices" << endl;
	}
	else
	{
		if( !quiet ) cout << "Sliced Acceptance is NOT using sorted horizontal slices" << endl;
	}
}

//.........这里部分代码省略.........
	gStyle->SetOptFit(0);



	/** Printing out efficiencies and systematics **/

	vector <TString> sysnames;
	sysnames.push_back("Lifetime");
	sysnames.push_back("Decay Model");
	sysnames.push_back("Polarization");
	if(type=="DD") sysnames.push_back("DD vtx");
	
	/** Print efficiencies */ 

	cout << "\n\n" << xvarname << " bin " << " \t\t\t& ";
	for(unsigned s = 0; s < effnames.size(); s++) cout << effnames[s] << " \t\t\t\t& ";
	cout << "Upper \t\t\t\t& Total  \\\\" << endl;

	TGraphErrors * grtoteff = new TGraphErrors();
	TGraphErrors * grtoteff_lowSel = new TGraphErrors();
	for(int j = 1; j <= toteff->GetNbinsX(); j++)
	{
		if((xbins[j]==11 && !rel) || xbins[j]==15) continue;
		if(xbins[j]==9.1 || xbins[j-1]==10.1) continue;
		cout << fixed << setprecision(1) << "eff " << xbins[j-1] << "-" << xbins[j] << fixed << setprecision(5) << " \t & ";

		for(unsigned i = 0; i < hdefault.size(); i++)	
			cout << "$" << hdefault[i]->GetBinContent(j) << " \\pm " << hdefault[i]->GetBinError(j) << "$ \t & ";

		cout << "$" << uppereff->GetBinContent(j) << " \\pm " << uppereff->GetBinError(j) << "$ \t & ";
		cout << "$" << toteff->GetBinContent(j) << " \\pm " << toteff->GetBinError(j) << "$ \\\\ " << endl;

		grtoteff->SetPoint(j,toteff->GetBinCenter(j),toteff->GetBinContent(j));
		grtoteff->SetPointError(j,toteff->GetBinWidth(j)/2.,toteff->GetBinError(j));
		if(toteff_lowSel) grtoteff_lowSel->SetPoint(j,toteff_lowSel->GetBinCenter(j),toteff_lowSel->GetBinContent(j));
		if(toteff_lowSel) grtoteff_lowSel->SetPointError(j,toteff_lowSel->GetBinWidth(j)/2.,toteff_lowSel->GetBinError(j));
	}

	grtoteff->Write("toteff");
	if(toteff_lowSel) grtoteff_lowSel->Write("toteff_lowSel");
	if(!doSys) { delete MCfile; delete histFile; return 0; }


	/** Print sys separate in efficiency */
	
	vector < TH1F * > tmp;
	vector < vector < TH1F * > > sys_eff(4,tmp);
	vector < TH1F * > tot_sys_eff;

	cout << endl << endl << endl;
	if(!percent) cout << "\n\n" << xvarname << " bin\t\t& Value \t & Stats";
	else cout << "\n\n" << xvarname << " bin ";
	for(unsigned s = 0; s < sysnames.size(); s++) cout << "\t& " << sysnames[s];
	cout << " \\\\" << endl;

	for(int j = 1; j <= toteff->GetNbinsX(); j++)
	{
		if((xbins[j]==11 && !rel) || xbins[j]==15) continue;
		if(xbins[j]==9.1 || xbins[j-1]==10.1) continue;
		cout << "-----------------------------------------------------------------------------------------" << endl;
		cout << fixed << setprecision(1) << xbins[j-1] << "-" << xbins[j] << fixed << setprecision(3) << endl;
		cout << "-----------------------------------------------------------------------------------------" << endl;
		if(!percent) cout << fixed << setprecision(5);

		for(unsigned i = 0; i < effnames.size(); i++)
		{

//.........这里部分代码省略.........
  std::string newName = std::string(inputfile).substr(0, std::string(inputfile).find(".root"));
  canv->Print(TString::Format("%s.png", newName.c_str())); 
  canv->Print(TString::Format("%s.pdf", newName.c_str())); 
  canv->Print(TString::Format("%s.eps", newName.c_str())); 

  /*
    Ratio Data over MC
  */
  TCanvas *canv0 = MakeCanvas("canv0", "histograms", 600, 400);
  canv0->SetGridx();
  canv0->SetGridy();
  canv0->cd(); 
  TH1F* model;
  if(CONSERVATIVE_CHI2){
    if(std::string(extra) == std::string("#mu#mu")){
      model = (TH1F*)Zmm ->Clone("model");
    }
    else if(std::string(extra) == std::string("ee")){
      model = (TH1F*)Zee ->Clone("model");
    }
    else if(std::string(extra) == std::string("e#mu") && HWWBG){
      model = (TH1F*)ggH_hww ->Clone("model");
    }
    else{  
      model = (TH1F*)Ztt ->Clone("model");
    }
  }
  else{
    model = (TH1F*)errorBand->Clone("model");
  }
  TH1F* test1 = (TH1F*)data->Clone("test1"); 
  for(int ibin=0; ibin<test1->GetNbinsX(); ++ibin){
    //the small value in case of 0 entries in the model is added to prevent the chis2 test from failing
    model->SetBinContent(ibin+1, model->GetBinContent(ibin+1)>0 ? model->GetBinContent(ibin+1)*model->GetBinWidth(ibin+1) : 0.01);
    model->SetBinError  (ibin+1, CONSERVATIVE_CHI2 ? 0. : model->GetBinError  (ibin+1)*model->GetBinWidth(ibin+1));
    test1->SetBinContent(ibin+1, test1->GetBinContent(ibin+1)*test1->GetBinWidth(ibin+1));
    test1->SetBinError  (ibin+1, test1->GetBinError  (ibin+1)*test1->GetBinWidth(ibin+1));
  }
  double chi2prob = test1->Chi2Test      (model,"PUW");        std::cout << "chi2prob:" << chi2prob << std::endl;
  double chi2ndof = test1->Chi2Test      (model,"CHI2/NDFUW"); std::cout << "chi2ndf :" << chi2ndof << std::endl;
  double ksprob   = test1->KolmogorovTest(model);              std::cout << "ksprob  :" << ksprob   << std::endl;
  double ksprobpe = test1->KolmogorovTest(model,"DX");         std::cout << "ksprobpe:" << ksprobpe << std::endl;  

  std::vector<double> edges;
  TH1F* zero = (TH1F*)Ztt->Clone("zero"); zero->Clear();
  TH1F* rat1 = (TH1F*)data->Clone("rat"); 
  for(int ibin=0; ibin<rat1->GetNbinsX(); ++ibin){
    rat1->SetBinContent(ibin+1, errorBand->GetBinContent(ibin+1)>0 ? data->GetBinContent(ibin+1)/errorBand->GetBinContent(ibin+1) : 0);
    rat1->SetBinError  (ibin+1, errorBand->GetBinContent(ibin+1)>0 ? data->GetBinError  (ibin+1)/errorBand->GetBinContent(ibin+1) : 0);
    zero->SetBinContent(ibin+1, 0.);
    zero->SetBinError  (ibin+1, errorBand->GetBinContent(ibin+1)>0 ? errorBand ->GetBinError  (ibin+1)/errorBand->GetBinContent(ibin+1) : 0);
  }
  for(int ibin=0; ibin<rat1->GetNbinsX(); ++ibin){
    if(rat1->GetBinContent(ibin+1)>0){
      edges.push_back(TMath::Abs(rat1->GetBinContent(ibin+1)-1.)+TMath::Abs(rat1->GetBinError(ibin+1)));
      // catch cases of 0 bins, which would lead to 0-alpha*0-1
      rat1->SetBinContent(ibin+1, rat1->GetBinContent(ibin+1)-1.);
    }
    zero->SetBinContent(ibin+1, 0.);
  }
  float range = 0.1;
  std::sort(edges.begin(), edges.end());
  if (edges[edges.size()-2]>0.1) { range = 0.2; }
  if (edges[edges.size()-2]>0.2) { range = 0.5; }
  if (edges[edges.size()-2]>0.5) { range = 1.0; }
  if (edges[edges.size()-2]>1.0) { range = 1.5; }

//.........这里部分代码省略.........
  if(!doCent)hBFractionMC2->GetXaxis()->SetRangeUser(80,200);
  hBFractionMC2->SetMarkerSize(0);
  hBFractionMC2->SetMaximum(0.06);
  hBFractionMC2->SetMinimum(0.0);
  hBFractionMC2->Draw("hist");


  TGraphAsymmErrors *gBFractionMC2 = new TGraphAsymmErrors(hBFractionMC);
  if(!doCent){
    setMeanPt(gBFractionMC2,hBFractionMC,0,do3bin);
    gBFractionMC2->GetXaxis()->SetRangeUser(80,200);
  }



  TLatex *prel;
  if(doCent)prel= new TLatex(10,0.0615,"CMS preliminary");
  else prel= new TLatex(85,0.0615,"CMS preliminary");
  prel->Draw();

  TLatex *roots = new TLatex(147,0.0615,"#sqrt{s_{NN}} = 2.76 TeV");
  roots->Draw();

  if(!doCent){
    TLatex *csel = new TLatex(90,0.05,"Centrality 0-100%");
    csel->Draw();
  }
  TLatex *ptlabel;
  if(doCent==1) ptlabel= new TLatex(20,0.005,"80 < Jet p_{T} < 100 GeV/c");
  if(doCent==2) ptlabel= new TLatex(20,0.005,"100 < Jet p_{T} < 120 GeV/c");
  if(doCent)ptlabel->Draw();

  // to be precise we should evaluate mcStatErr for 3 bins seperately
  float mcStatErr[4] = {0.03,0.06,0.07,0.15};

  
  //TGraphAsymmErrors *gSyst = new TGraphAsymmErrors(3);
  TGraphErrors *gSyst;
  if(do3bin) gSyst= new TGraphErrors(3);
  else gSyst= new TGraphErrors(4);
  Double_t errCLratio, errMethod, totalSystErr;

  for(Int_t i=1;i<=hBFractionDataLTJP->GetNbinsX();i++) {
    gSyst->SetPoint(i-1,hBFractionDataLTJP->GetBinCenter(i),hBFractionDataLTJP->GetBinContent(i));
    cout<<" central value "<<hBFractionDataLTJP->GetBinContent(i)<<endl;
    errCLratio = abs(hBFractionDataLTJP->GetBinContent(i)-hBFractionDataLTJPFixC->GetBinContent(i));
    errMethod = max(abs(hBFractionDataLTJP->GetBinContent(i)-hBFractionData->GetBinContent(i)),abs(hBFractionDataLTJP->GetBinContent(i)-hBFractionJPdirect->GetBinContent(i)));
    double errJES = 0.14*hBFractionDataLTJP->GetBinContent(i);
    
    totalSystErr = norm(errCLratio,errMethod,errJES);
    gSyst->SetPointError(i-1,hBFractionDataLTJP->GetBinWidth(i)/2,totalSystErr);
    cout<<" sys error "<<totalSystErr<<endl;
    // add in MC template uncertainties
    float origStatErr = hBFractionDataLTJP->GetBinError(i);
    int statBin=i-1;
    if(doCent==1) statBin=0;
    if(doCent==2) statBin=1;
    float extraStatErr = mcStatErr[statBin]*hBFractionDataLTJP->GetBinContent(i);
    float totalStatErr = sqrt(origStatErr*origStatErr + extraStatErr*extraStatErr);
    hBFractionDataLTJP->SetBinError(i,totalStatErr);
    cout<<" total error "<<sqrt(totalSystErr*totalSystErr+totalStatErr*totalStatErr)<<endl;

  }

  gSyst->SetFillColor(5);
  gSyst->Draw("2");


  gBFractionMC2->Draw("Z,p,same");
  hBFractionMC2->Draw("hist,same");


  TGraphAsymmErrors *gBFractionDataLTJP2 = new TGraphAsymmErrors(hBFractionDataLTJP);
  if(!doCent)setMeanPt(gBFractionDataLTJP2,hBFractionDataLTJP,1,do3bin);
  gBFractionDataLTJP2->SetLineColor(1);
  gBFractionDataLTJP2->SetMarkerColor(1);
  gBFractionDataLTJP2->SetMarkerSize(1.5);
  gBFractionDataLTJP2->Draw("p,e1,same");


  TLegend *legFrac2 = new TLegend(0.2,0.15,0.8,0.34);
  if(doCent){
    legFrac2->SetX1(0.365);
    legFrac2->SetY1(0.657);
    legFrac2->SetX2(0.965);
    legFrac2->SetY2(0.848);
  }

  legFrac2->SetHeader("#int L dt = 150 #mub^{-1}");
  legFrac2->SetBorderSize(0);
  legFrac2->SetFillStyle(0);
  legFrac2->AddEntry(gBFractionDataLTJP2,"PbPb data","p");
  legFrac2->AddEntry(gBFractionMC2,"PYTHIA+HYDJET 1.8","lp");
  legFrac2->AddEntry(gSyst,"Exp. uncertainty","f");
  legFrac2->Draw();

 
  cBFraction2->RedrawAxis();

}

void DrawMass(int iSEL)
{
  gROOT->ForceStyle();
  TString SET[2] = {"A","B"};
  RooMsgService::instance().setSilentMode(kTRUE);
  for(int i=0;i<2;i++) {
    RooMsgService::instance().setStreamStatus(i,kFALSE);
  }
  TString SELECTION[2] = {"NOM","VBF"};
  TString MASS_RAW[2]  = {"mbb[1]","mbb[2]"};
  TString MASS_REG[2]  = {"mbbReg[1]","mbbReg[2]"};

 // TFile *inf = TFile::Open("Fit_VBFPowheg125_sel"+SELECTION[iSEL]+".root");
  TFile *inf = TFile::Open("/usb/data2/UAData/2015/flatTree_VBFPowheg125.root");
  TTree *tr  = (TTree*)inf->Get("Hbb/events");

  TH1F *hRaw = new TH1F("hRawMass","hRawMass",150,0,300);
  TH1F *hReg = new TH1F("hRegMass","hRegMass",150,0,300);

  RooRealVar x("mbb","mbb",60,170);

  TCanvas *can = new TCanvas("Mbb_sel"+SELECTION[iSEL],"Mbb_sel"+SELECTION[iSEL],900,750);
  
  TCut ct1 = TCut("triggerResult[0]==1||triggerResult[1]==1");
  TCut cs1 = TCut("jetBtag[b1[1]]>0.244 && jetBtag[b2[1]]>0.244 && jetPt[3]>40. && jetPt[2]>50. && jetPt[1]>70. && jetPt[0]>80. && dEtaqq[1]>2.5 && mqq[1]>250 && dPhibb[1]<2.0 && nLeptons==0");
  tr->Draw(MASS_RAW[iSEL]+">>hRawMass",ct1&&cs1);
  tr->Draw(MASS_REG[iSEL]+">>hRegMass",ct1&&cs1);

  hRaw->Sumw2();
  hReg->Sumw2();
  hRaw->Scale(1./(hRaw->Integral()*hRaw->GetBinWidth(1)));
  hReg->Scale(1./(hReg->Integral()*hReg->GetBinWidth(1)));

  RooDataHist *rRaw = new RooDataHist("rRaw","rRaw",x,hRaw);
  RooDataHist *rReg = new RooDataHist("rReg","rReg",x,hReg);

  RooRealVar m1("m1","m1",125,110,140); 
  RooRealVar m2("m2","m2",125,110,140);
  RooRealVar sL1("sL1","sL1",12,3,30); 
  RooRealVar sL2("sL2","sL2",12,3,30);
  RooRealVar sR1("sR1","sR1",12,3,30); 
  RooRealVar sR2("sR2","sR2",12,3,30);
  RooRealVar a1("a1","a1",1,-10,10); 
  RooRealVar a2("a2","a2",1,-10,10); 
  RooRealVar n1("n1","n1",1,0,100); 
  RooRealVar n2("n2","n2",1,0,100);        
  RooRealVar b10("b10","b10",0.5,0.,1.);
  RooRealVar b11("b11","b11",0.5,0.,1.);
  RooRealVar b12("b12","b12",0.5,0.,1.);
  RooRealVar b13("b13","b13",0.5,0.,1.);
  RooRealVar b20("b20","b20",0.5,0.,1.);
  RooRealVar b21("b21","b21",0.5,0.,1.);
  RooRealVar b22("b22","b22",0.5,0.,1.); 
  RooRealVar b23("b23","b23",0.5,0.,1.);     
        
  RooBernstein bkg1("bkg1","bkg1",x,RooArgSet(b10,b11,b12,b13));
  RooBernstein bkg2("bkg2","bkg2",x,RooArgSet(b20,b21,b22,b23));
  RooRealVar fsig1("fsig1","fsig1",0.7,0.,1.); 
  RooRealVar fsig2("fsig2","fsig2",0.7,0.,1.);     
  RooBifurGauss sig1("sig1","sig1",x,m1,sL1,sR1);
  RooBifurGauss sig2("sig2","sig2",x,m2,sL2,sR2);
  //RooCBShape sig1("sig1","sig1",x,m1,s1,a1,n1);
  //RooCBShape sig2("sig2","sig2",x,m2,s2,a2,n2);
        
  RooAddPdf *model1 = new RooAddPdf("model1","model1",RooArgList(sig1,bkg1),fsig1);
  RooAddPdf *model2 = new RooAddPdf("model2","model2",RooArgList(sig2,bkg2),fsig2);

  model1->fitTo(*rRaw,SumW2Error(kFALSE),"q");
  model2->fitTo(*rReg,SumW2Error(kFALSE),"q");

  hRaw->SetLineWidth(2);
  hReg->SetLineWidth(2);
  hRaw->SetLineColor(kBlack);
  hReg->SetLineColor(kRed+1); 
  hReg->SetFillColor(kRed-10);
  hRaw->SetMarkerStyle(21);
  hReg->SetMarkerStyle(20);
  hRaw->SetMarkerSize(1.5);
  hReg->SetMarkerSize(1.5);
  hRaw->SetMarkerColor(kBlack);
  hReg->SetMarkerColor(kRed+1);
  
  RooPlot* frame = x.frame();
  rRaw->plotOn(frame,LineColor(kBlack),LineWidth(1),MarkerColor(kBlack),MarkerStyle(21));
  model1->plotOn(frame,LineColor(kBlack),LineWidth(2));
  rReg->plotOn(frame,LineColor(kRed+1),LineWidth(1),MarkerColor(kRed+1),MarkerStyle(20));
  model2->plotOn(frame,LineColor(kRed+1),LineWidth(2));

  
  TF1 *tmp_func1 = model1->asTF(x,fsig1,x);
  TF1 *tmp_func2 = model2->asTF(x,fsig2,x);
  double y01 = tmp_func1->GetMaximum();
  double x01 = tmp_func1->GetMaximumX();
  double x11 = tmp_func1->GetX(y01/2,60,x01);
  double x21 = tmp_func1->GetX(y01/2,x01,200);
  double FWHM1 = x21-x11;
  
  double y02 = tmp_func2->GetMaximum();
  double x02 = tmp_func2->GetMaximumX();
  double x12 = tmp_func2->GetX(y02/2,60,x02);
//.........这里部分代码省略.........