C++ THStack::GetStack方法代码示例

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

//.........这里部分代码省略.........
   histo[iName][iCut]->SetFillColor(vColor[iName]);
   histo[iName][iCut]->SetLineWidth(2);
   histo[iName][iCut]->SetFillStyle(3001);
   
   bool isSig = false;
   for (std::vector<std::string>::const_iterator itSig = SignalName.begin(); itSig != SignalName.end(); itSig++){
    if (reduced_name_samples.at(iName) == *itSig) isSig = true;
   }
   
   if (!isSig && reduced_name_samples.at(iName) != "DATA") {
    hs[iCut]->Add(histo[iName][iCut]);
   }
   else {
    if (!isSig) {
     histo[iName][iCut]->SetMarkerStyle(20);
     histo[iName][iCut]->SetMarkerSize(1);
     histo[iName][iCut]->SetMarkerColor(kBlack);
     histo[iName][iCut]->SetLineColor(kBlack);
     histo[iName][iCut]->SetFillColor(kBlack);
     histo[iName][iCut]->SetLineWidth(2);
     histo[iName][iCut]->SetFillStyle(3001);  
    }
    else {
     histo[iName][iCut]->SetMarkerStyle(21);
     histo[iName][iCut]->SetMarkerSize(1);
     histo[iName][iCut]->SetLineWidth(2);
     histo[iName][iCut]->SetFillStyle(3001);  
    }
   }
  }
  ///==== histo sum MC ====    
  ///==== Add systrematic error ====
  AddError(hs[iCut],LumiSyst);
  histoSumMC[iCut] = ((TH1F*)(hs[iCut]->GetStack()->Last()));

  std::cout << " MC / DATA[" << iCut << "] = "<< histoSumMC[iCut]->Integral() << " / " << histo[numDATA][iCut]->Integral() << " = " << (histo[numDATA][iCut]->Integral() ? histoSumMC[iCut]->Integral()/ histo[numDATA][iCut]->Integral() : 0) << std::endl;
  
  ///==== legend ====
  if (!LegendBuilt){
   for (unsigned int iName=0; iName<reduced_name_samples.size(); iName++){
    leg->AddEntry(histo[iName][iCut],reduced_name_samples.at(iName).c_str(),"pf");    
    LegendBuilt = true;
   }
  }
 }
 std::cout << std::endl << std::endl;
 
 ///==== calculate number of events after each step of the analysis ====
 //  [iName][iCut]
 hsTrend = new THStack("Trend","Trend");
 
 for (unsigned int iCut = 0; iCut<vCut.size(); iCut++){
  TString nameTHTrendPie = Form("%d_Trend_Pie",iCut);
  hTrendPie[iCut] = new TPie (nameTHTrendPie,nameTHTrendPie,reduced_name_samples.size());
 }
 
 for (unsigned int iName=0; iName<reduced_name_samples.size(); iName++){
  TString nameTHTrend = Form("%s_Trend",reduced_name_samples.at(iName).c_str());
  hTrend[iName] = new TH1F (nameTHTrend,nameTHTrend,vCut.size()+1,0,vCut.size()+1);
  hTrend[iName]->GetXaxis()->SetTitle("Selections");

  if (iName == numDATA) {
   hTrend[iName]->SetMarkerStyle(20);
   hTrend[iName]->SetMarkerSize(1);
   hTrend[iName]->SetMarkerColor(kBlack);
   hTrend[iName]->SetLineColor(kBlack);

//.........这里部分代码省略.........
       //fixed bin histo
       else hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xLow,xHigh);
       hist[theHistCounter] -> Sumw2();
       hist[theHistCounter] -> SetFillColor(col[theHistCounter]);
       hist[theHistCounter] -> SetFillStyle(1001);
       theLegends.push_back(*listOfSamples.at(iSample)->Legend());
    }

    //fill the histogram
    TString thisScale = Form("%f *", *(listOfSamples.at(iSample)->Scale()));
    if (isFirstOfSerie) tree[iSample] -> Draw(myVar + " >> " + TString(hist[theHistCounter] -> GetName()),thisScale + "evt_weight*kf_weight*pu_weight" + myCut);
    else tree[iSample] -> Draw(myVar + " >>+ " + TString(hist[theHistCounter] -> GetName()),thisScale + "evt_weight*kf_weight*pu_weight" + myCut);
    
    //add the histogram to the stack if the last of the series:
    //either last sample or ~ sample followed by non ~ sample
    if (isLastOfSerie) {
       if (xlowVec != 0) {
         for (int iBin = 1; iBin <= nBins; iBin++) hist[theHistCounter]->SetBinError  (iBin,hist[theHistCounter]->GetBinError(iBin)/hist[theHistCounter]->GetBinWidth(iBin));
         for (int iBin = 1; iBin <= nBins; iBin++) hist[theHistCounter]->SetBinContent(iBin,hist[theHistCounter]->GetBinContent(iBin)/hist[theHistCounter]->GetBinWidth(iBin));
       }
       hs -> Add(hist[theHistCounter]);
       theHistCounter++;
    }
    
  }//end loop on samples

  //Fix the legend
  for (int iHisto = theHistCounter-1; iHisto >= 0; iHisto--) {
    leg -> AddEntry(hist[iHisto], theLegends[iHisto], "f");   
  }
  
  //get the maximum to properly set the frame
  float theMax = hdata -> GetBinContent(hdata -> GetMaximumBin()) + hdata -> GetBinError(hdata -> GetMaximumBin());
  TH1* theMCSum = (TH1*) hs->GetStack()->Last();
  float theMaxMC = theMCSum->GetBinContent(theMCSum->GetMaximumBin()) + theMCSum->GetBinError(theMCSum->GetMaximumBin());
  if (theMaxMC > theMax) theMax = theMaxMC;
  
  //prepare the ratio band and plot
  TH1* theMCRatioBand = makeRatioBand(theMCSum);
  TH1* theRatioPlot = makeRatioPlot(hdata,theMCSum);
    
  TCanvas* can = new TCanvas();
  can -> SetLogy(isLog);
  
  TPad *pad1 = new TPad("pad1","top pad",0,0.30,1,1);
  pad1->SetBottomMargin(0.02);
  pad1->SetLeftMargin(0.13);
  pad1->Draw();
  TPad *pad2 = new TPad("pad2","bottom pad",0,0.0,1,0.30);
  pad2->SetTopMargin(0.02);
  pad2->SetLeftMargin(0.13);
  pad2->SetBottomMargin(0.4);
  pad2->SetGridy();
  pad2->Draw();
  
  pad1->cd();
  hs->Draw("hist");
  hdata->Draw("same,pe");
  if (drawSignal) hsignal->Draw("same,hist");
  if (drawLegend) leg->Draw("same");
  //hs->GetXaxis()->SetTitle(myAxisNameX);
  hs->GetYaxis()->SetTitle(myAxisNameY);
  hs->GetXaxis()->SetLabelSize(0.04);
  hs->GetYaxis()->SetLabelSize(0.04);
  hs->GetXaxis()->SetLabelOffset(0.025);
  hs->GetYaxis()->SetLabelOffset(0.035);

CalcMtEff(string outStr=""){
gStyle->SetOptStat(0);  ///to avoid the stat. on the plots
char tempname[200];
int W = 600;
int H = 600;
int H_ref = 600;
int W_ref = 800;
float T = 0.08*H_ref;
float B = 0.12*H_ref;
float L = 0.12*W_ref;
float R = 0.04*W_ref;

TCanvas* c1 = new TCanvas("name","name",10,10,W,H);
c1->SetFillColor(0);
c1->SetBorderMode(0);
c1->SetFrameFillStyle(0);
c1->SetFrameBorderMode(0);
c1->SetLeftMargin( L/W );
c1->SetRightMargin( R/W );
c1->SetTopMargin( T/H );
c1->SetBottomMargin( B/H );
c1->SetTickx(0);
c1->SetTicky(0);



Float_t legendX1 = .7; //.50;
Float_t legendX2 = .9; //.70;
Float_t legendY1 = .70; //.65;
Float_t legendY2 = .85;
TLegend* catLeg1 = new TLegend(legendX1,legendY1,legendX2,legendY2);
catLeg1->SetTextSize(0.032);
catLeg1->SetTextFont(42);
catLeg1->SetTextSize(0.04);
catLeg1->SetFillColor(0);
catLeg1->SetLineColor(0);
catLeg1->SetBorderSize(0);

TH1D * thist_tt, * thist_wj;
TH1D * thist_tt2, * thist_wj2;
TH1D * thist_t, * thist_t2;
THStack * stack;

TFile * after_tt = new TFile("Stack/Elog410_MTSelMinus_HadTauEstimation_TTbar_stacked.root","R");
TFile * before_tt = new TFile("Stack/Elog410_BeoreMT_HadTauEstimation_TTbar_stacked.root","R");
TFile * after_wj = new TFile("Stack/Elog410_MTSelMinus_HadTauEstimation_WJet_stacked.root","R");
TFile * before_wj = new TFile("Stack/Elog410_BeoreMT_HadTauEstimation_WJet_stacked.root","R");
TFile * after_t = new TFile("Stack/Elog410_MTSelMinus_HadTauEstimation_T_stacked.root","R");
TFile * before_t = new TFile("Stack/Elog410_BeoreMT_HadTauEstimation_T_stacked.root","R");

sprintf(tempname,"MtEff_%s_.root",outStr.c_str());
TFile * outFile = new TFile(tempname,"RECREATE");

////////
// ttbar
////////

stack = (THStack *) after_tt->Get("searchH")->Clone("after");
thist_tt = (TH1D *) stack->GetStack()->Last();
stack = (THStack *) before_tt->Get("searchH")->Clone("before");
thist_tt2 = (TH1D *) stack->GetStack()->Last(); 

TH1D * MtCutEff_tt = static_cast<TH1D*>(thist_tt->Clone("MtCutEff_tt"));
MtCutEff_tt->Divide(thist_tt,thist_tt2,1,1,"B");
//
//one time drawing options
//
double XUp = 46. , maxVal=2.;
  MtCutEff_tt->SetMaximum(maxVal);
    MtCutEff_tt->SetTitle("");
    MtCutEff_tt->GetXaxis()->SetLabelFont(42);
    MtCutEff_tt->GetXaxis()->SetLabelOffset(0.007);
    MtCutEff_tt->GetXaxis()->SetLabelSize(0.04);
    MtCutEff_tt->GetXaxis()->SetTitleSize(0.04);
    MtCutEff_tt->GetXaxis()->SetTitleOffset(0.9);
    MtCutEff_tt->GetXaxis()->SetTitleFont(42);
    MtCutEff_tt->GetYaxis()->SetLabelFont(42);
    MtCutEff_tt->GetYaxis()->SetLabelOffset(0.007);
    MtCutEff_tt->GetYaxis()->SetLabelSize(0.04);
    MtCutEff_tt->GetYaxis()->SetTitleSize(0.04);
    MtCutEff_tt->GetYaxis()->SetTitleOffset(1.25);
    MtCutEff_tt->GetYaxis()->SetTitleFont(42);
    MtCutEff_tt->GetXaxis()->SetTitle("search bins");
    MtCutEff_tt->GetYaxis()->SetTitle("MT efficiency");
    MtCutEff_tt->GetXaxis()->SetRangeUser(1.,XUp);
    MtCutEff_tt->SetMaximum(maxVal);
//
//end of drawing options
//
MtCutEff_tt->SetLineColor(1);
MtCutEff_tt->Draw("same");
MtCutEff_tt->Write();

////////
// wjet
////////
stack = (THStack *) after_wj->Get("searchH")->Clone("after");
thist_wj = (TH1D *) stack->GetStack()->Last();
stack = (THStack *) before_wj->Get("searchH")->Clone("before");
thist_wj2 = (TH1D *) stack->GetStack()->Last();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // set pad size
  canvas_up->SetPad(0., 1 - up_height, 0.97, 1.);
  canvas_dw->SetPad(0., dw_height_offset, 0.97, dw_height+dw_height_offset);
  canvas_up->SetFrameFillColor(0);
  canvas_up->SetFillColor(0);
  canvas_dw->SetFillColor(0);
  canvas_dw->SetFrameFillColor(0);
  canvas_dw->SetBottomMargin(0.25);
  
  // set top margin 0 for bottom figure
  canvas_dw->SetTopMargin(0);
  
  // draw top figure
  canvas_up->cd();

  TH1D * GenHist, * EstHist,* thist;
  TH1D * histTemplate;
  THStack * tempstack;

  double HT_x_max=2500.;
  double HT_x_min=400.;
  double MHT_x_max=1000.;
  double NJet_x_max=15.;
  double NBtag_x_max=4.;
  double search_x_max=19.;
  double search_x_min=1.;


  sprintf(tempname,"%s",histname.c_str());
  if(sample.find("stacked")!=string::npos){
   
    EstHist=(TH1D*)EstFile->Get(tempname)->Clone();
    tempstack=(THStack*)GenFile->Get(tempname)->Clone();   
    GenHist=(TH1D*) tempstack->GetStack()->Last();
    
    EstHistD=(TH1D*) EstFile->Get(tempname)->Clone();
    tempstack=(THStack*)GenFile->Get(tempname)->Clone();
    GenHistD=(TH1D*) tempstack->GetStack()->Last();
    
  }
  else{
    EstHist=(TH1D*) EstFile->Get(tempname)->Clone();
    GenHist=(TH1D*) GenFile->Get(tempname)->Clone();
    EstHistD=(TH1D*) EstFile->Get(tempname)->Clone();
    GenHistD=(TH1D*) GenFile->Get(tempname)->Clone();
  }

  GenHist->SetLineColor(2);
  EstHist->SetLineColor(4);
  //GenHist->GetXaxis()->SetLabelFont(42);
  //GenHist->GetXaxis()->SetLabelOffset(0.007);
  //GenHist->GetXaxis()->SetLabelSize(0.04);
  GenHist->GetXaxis()->SetTitleSize(0.05);
  GenHist->GetXaxis()->SetTitleOffset(0.9);
  GenHist->GetXaxis()->SetTitleFont(42);
  //GenHist->GetYaxis()->SetLabelFont(42);
  //GenHist->GetYaxis()->SetLabelOffset(0.007);
  //GenHist->GetYaxis()->SetLabelSize(0.04);
  GenHist->GetYaxis()->SetTitleSize(0.05);
  GenHist->GetYaxis()->SetTitleOffset(1.25);
  GenHist->GetYaxis()->SetTitleFont(42);

  //EstHist->GetXaxis()->SetLabelFont(42);
  //EstHist->GetXaxis()->SetLabelOffset(0.007);
  //EstHist->GetXaxis()->SetLabelSize(0.04);
  EstHist->GetXaxis()->SetTitleSize(0.05);

OptRes Optimization(string var = "topness", string signal = "T2tt_850_100", string region = "baseline", string channel = "allChannels", string filename = "plots/plotsProducer/1DStack.root"){

	OptRes res;
	
	//TFile* fplots = new TFile("plots/plotsProducer/1DStack.root","READ");
	TFile* fplots = new TFile(filename.c_str(),"READ");
	//string channel = "allChannels";
	//string region = "baseline";
	//string var = "topness";
	//string var = "topness_m5";
	//string signal = "T2tt_850_100";
	
	
	//retrieve canvas
	string cname = channel+"/"+region+"/"+var;
	TCanvas* c = fplots->Get(cname.c_str());
	//c->Draw();
	
	//signal
	string splot_name = "v:"+var+"|p:"+signal+"|r:"+region+"|c:"+channel+"|t:1DEntries";
	TH1F* h_sig = (TH1F*) c->GetPrimitive(splot_name.c_str());
	//sig->Draw();
	
	//////////////////////////////////////////////
	// Search the Stack that contains all bkg
	//////////////////////////////////////////////
	vector<TH1F*> h_bkgs;
	//c->GetListOfPrimitives()->Print();
	TIter next(c->GetListOfPrimitives());
	TObject* h;
	while ((h = (TH1F*) next())){
	//obj->Draw(next.GetOption());
		if(h->InheritsFrom("THStack")){
			THStack* stack = (THStack*) h;
			for(int i=0;i<stack->GetStack()->GetEntries();i++){
				TH1F* hist = (TH1F*) stack->GetStack()->At(i);
				//cout<<hist->IsA()->GetName()<<endl;
				//cout<<hist->GetName()<<" "<<hist->Integral()<<endl;
				//search tt_2l
				//if(string(hist->GetName()).find("W+jets")!=string::npos) continue;
				if(string(hist->GetName()).find("ttbar_2l")!=string::npos){
					//cout<<"mean = "<<hist->GetMean()<<" rms = "<<hist->GetRMS()<<" kurtosis = "<<hist->GetKurtosis()<<" skewness = "<<hist->GetSkewness()<<endl;
				}
				//else continue;
				//conly save the last one because plots are cumulative !!!
				// that's cheat ....
				if(i==stack->GetStack()->GetEntries()-1) h_bkgs.push_back(hist);
			}
			//cout<<((THStack*)h)->GetStack()->GetEntries()<<endl;
			//TH1* h_bkg = ((TH1*)(h->GetStack()->Last()));
		}
	}
	
	
	////////////////////////////////////////////////
	//     Loop over all bins and compute eff/sig
	////////////////////////////////////////////////
	
	res.h_eff_sig = (TH1F*) h_sig->Clone();	
	res.h_eff_bkg = (TH1F*) h_sig->Clone();	
	res.h_Zbi = (TH1F*) h_sig->Clone();	
	TH1F* h_ROC = new TH1F("h_ROC","ROC curve",10,0,1);

	//Compute the integral
	float integ_sig = res.h_eff_sig->Integral(0,res.h_eff_sig->GetNbinsX()+1);
	float integ_bkg = 0;
	for(unsigned ih=0;ih<h_bkgs.size();ih++){
		integ_bkg+=h_bkgs[ih]->Integral(0,h_bkgs[ih]->GetNbinsX()+1);
	}
	///////////////////////

	//-- efficiencoy for signal
	for(int i=1;i<res.h_eff_sig->GetNbinsX()+1;i++){
		float eff = 0;
		if(integ_sig!=0)  eff = res.h_eff_sig->Integral(i,res.h_eff_sig->GetNbinsX()+1)/integ_sig;
		res.h_eff_sig->SetBinContent(i,eff);
		res.h_eff_sig->SetBinError(i,0);
	}
	//-- efficiency for bkg
	for(int i=1;i<res.h_eff_sig->GetNbinsX()+1;i++){
		float eff = 0;
		for(unsigned ih=0;ih<h_bkgs.size();ih++){
			eff += h_bkgs[ih]->Integral(i,h_bkgs[ih]->GetNbinsX()+1);
		}
		if(integ_bkg!=0)  eff /= integ_bkg;
		res.h_eff_bkg->SetBinContent(i,eff);
		res.h_eff_bkg->SetBinError(i,0);
	}
	//-- compute significance
	//vector<pair<float,float> > roc_v;
	float* roc_sig = new float[res.h_eff_sig->GetNbinsX()];
	float* roc_bkg = new float[res.h_eff_sig->GetNbinsX()];
	int imax = 0;
	for(int i=1;i<res.h_eff_sig->GetNbinsX()+1;i++){
		float Ns = res.h_eff_sig->GetBinContent(i)*integ_sig;
		//float Ns = res.h_eff_sig->GetBinContent(i);
		//float Nb = res.h_eff_bkg->GetBinContent(i);
		float Nb = res.h_eff_bkg->GetBinContent(i)*integ_bkg;
		float signif = 0;
		if(Nb!=0) signif = Ns/sqrt(Nb);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    treedata[iDatas] = (TTree*) infile -> Get(listOfDatasets.at(iDatas)->Name()->Data());

    //fill the histogram
    if ( iDatas == 0 ) treedata[iDatas] -> Draw(myVar + " >> hdata","hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
    else treedata[iDatas] -> Draw(myVar + " >>+ hdata","hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
    
    if ( iDatas == 0 ) leg -> AddEntry(hdata, "DATA (19.8 fb^{-1})", "pl");    
    
  }//end loop on datasets
  if (xlowVec != 0) {
    for (int iBin = 1; iBin <= nBins; iBin++) hdata->SetBinError  (iBin,hdata->GetBinError(iBin)/hdata->GetBinWidth(iBin));
    for (int iBin = 1; iBin <= nBins; iBin++) hdata->SetBinContent(iBin,hdata->GetBinContent(iBin)/hdata->GetBinWidth(iBin));
  }
       
  int theHistCounter = 0;
  // loop through the samples and produce the plots
  for (UInt_t iSample=0; iSample < listOfSamples.size(); iSample++) {

    //determine if the histo is first of the series
    bool isFirstOfSerie = (*listOfSamples.at(iSample)->Legend()).CompareTo(" ");
    bool isLastOfSerie = false;
    if (iSample == listOfSamples.size() - 1) isLastOfSerie = true;
    if (iSample < listOfSamples.size() - 1 && (*listOfSamples.at(iSample+1)->Legend()).CompareTo(" ") != 0) isLastOfSerie = true;
    
    //get the tree
    tree[iSample] = (TTree*) infile -> Get(listOfSamples.at(iSample)->Name()->Data());
    //if sample first of the list create a new histogram
    if (isFirstOfSerie) {
       TString thisHistName = "h_" + *(listOfSamples.at(iSample)->Name());
       //variable bin histo
       if (xlowVec != 0) hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xlowVec);
       //fixed bin histo
       else hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xLow,xHigh);
       hist[theHistCounter] -> Sumw2();
       hist[theHistCounter] -> SetFillColor(col[theHistCounter]);
       hist[theHistCounter] -> SetFillStyle(1001);
       theLegends.push_back(*listOfSamples.at(iSample)->Legend());
       cout << *listOfSamples.at(iSample)->Legend() << " ";
    }

    //fill the histogram
    
    TString thisScale = Form("%f *", *(listOfSamples.at(iSample)->Scale()));
    if (isFirstOfSerie) tree[iSample] -> Draw(myVar + " >> " + TString(hist[theHistCounter] -> GetName()),thisScale + "hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
    else tree[iSample] -> Draw(myVar + " >>+ " + TString(hist[theHistCounter] -> GetName()),thisScale + "hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
        
    //add the histogram to the stack if the last of the series:
    //either last sample or ~ sample followed by non ~ sample
    if (isLastOfSerie) {
       if (xlowVec != 0) {
         for (int iBin = 1; iBin <= nBins; iBin++) hist[theHistCounter]->SetBinError  (iBin,hist[theHistCounter]->GetBinError(iBin)/hist[theHistCounter]->GetBinWidth(iBin));
         for (int iBin = 1; iBin <= nBins; iBin++) hist[theHistCounter]->SetBinContent(iBin,hist[theHistCounter]->GetBinContent(iBin)/hist[theHistCounter]->GetBinWidth(iBin));
       }
       hs -> Add(hist[theHistCounter]);
       cout << hist[theHistCounter] -> GetBinContent(1) << " +/- " <<  hist[theHistCounter] -> GetBinError(1) << endl;
       if ( listOfSamples.at(iSample)->Name()->Contains("wgptg130") ) {
         nSig +=  hist[theHistCounter] -> GetBinContent(1) ;
         nSigErrSq +=  pow(hist[theHistCounter] -> GetBinError(1),2) ;
       }
       else {
         nBkg +=  hist[theHistCounter] -> GetBinContent(1) ;
         nBkgErrSq +=  pow(hist[theHistCounter] -> GetBinError(1),2) ;
       }
       theHistCounter++;
    }
    
  }//end loop on samples

  //Fix the legend
  for (int iHisto = theHistCounter-1; iHisto >= 0; iHisto--) {
    leg -> AddEntry(hist[iHisto], theLegends[iHisto], "f");   
  }
  
  //get the maximum to properly set the frame
  float theMax = hdata -> GetBinContent(hdata -> GetMaximumBin()) + hdata -> GetBinError(hdata -> GetMaximumBin());
  TH1* theMCSum = (TH1*) hs->GetStack()->Last();
  cout << "total MC " << theMCSum -> GetBinContent(1) << " +/- " << theMCSum -> GetBinError(1) << endl;
  float theMaxMC = theMCSum->GetBinContent(theMCSum->GetMaximumBin()) + theMCSum->GetBinError(theMCSum->GetMaximumBin());
  if (theMaxMC > theMax) theMax = theMaxMC;
  
    
  if (isBlind) cout << "total DATA " << hdata -> GetBinContent(1) << " +/- " <<  hdata -> GetBinError(1) << endl;

  //compute the scale factor
  float nData =  hdata -> GetBinContent(1);
  float nDataErr =  hdata -> GetBinError(1);
  newScaleFactor = (nData - nBkg)/(nSig/defScaleFactor);
  float newScaleFactorErr = sqrt( (pow(nDataErr,2) + nBkgErrSq)/(nData - nBkg)/(nData - nBkg) + nSigErrSq/nSig/nSig ) * newScaleFactor;
  cout << "\n The scale factor is " << newScaleFactor << " +/- " << newScaleFactorErr << endl;
  
  //cleanup the memory allocation
  delete theMCSum;
  delete hs;
  delete leg;
  delete hdata;
  infile -> Close();
  delete infile;
  
  return;
}

//.........这里部分代码省略.........
  outputFile->cd();

  //--//
  THStack *hs = new THStack("hs","");
  // hs->Add(tauPt_Stage1_WJetsToLNu);
  hs->Add(total_diboson);
  hs->Add(total_DY);
  hs->Add(total_ST);
  hs->Add(total_TT);
  //  hs->Add(tauPt_Stage1_DYJetsToLL_M50);
  hs->Add(total_QCD);
  //  hs->Add(total_ZJets);
  hs->Add(total_WJets);
  //tauPt_Stage1_->SetTitle("");
    /*
  eff1->GetYaxis()->SetTitle("Events");
  //eff1->SetMarkerStyle(0);
  //eff1->SetMarkerColor(kBlack);
  eff1->SetLineColor(kBlack);
  eff1->SetLineWidth(2);
  */
  TCanvas* my_canvas = new TCanvas("canvas","canvas",800,600);
  my_canvas->cd();
  //  hs->Draw();
  gPad->SetLogy();
  hs->Draw("HIST");                                                                                                                                        
  hs->SetMaximum(100000);
  hs->SetMinimum(0.1);
  hs->GetXaxis()->SetRangeUser(0, 1000);
  //  hs->GetXaxis()->SetLimits(40, 3200);
  hs->GetXaxis()->SetTitle("#tau pT [GeV]");
  hs->GetYaxis()->SetTitle("Events");

  TH1F* total =  (TH1F*)hs->GetStack()->Last()->Clone();
  // total->SetFillStyle(3004);                                                                                                                                      
  // total->SetFillColor(kGray+1); 
  //  total->SetMarkerColor(0);
  // total->Draw("SAME E2");
  //hs->GetStack()->Last()->Draw("same E");
  // hs->Draw("SAME HIST");

  //  tauPt_Stage1_Run2016E->Draw("SAME E0");
  tauPt_Stage1_Run2016all->Draw("SAME E0");
  TH1F* mydata =  (TH1F*)tauPt_Stage1_Run2016all->Clone();

  //tauPt_Stage1_Run2016C->Draw("SAME E0");
  //tauPt_Stage1_Run2016CDE->Draw("SAME E0");
  //
  tauPt_Stage1_Wprime_M4000->Draw("SAME HIST");
  //  hs->SetOption("HIST L");
  
  TLatex* CMS_text = new TLatex(0.20,0.90,"CMS");
  CMS_text->SetNDC();
  CMS_text->SetTextSize(0.05);
  CMS_text->SetTextAngle(0);
  CMS_text->Draw("same");
  TLatex* CMS_text_2 = new TLatex(0.20,0.85,"Preliminary");
  CMS_text_2->SetNDC();
  CMS_text_2->SetTextFont(42);
  CMS_text_2->SetTextSize(0.05);
  CMS_text_2->SetTextAngle(0);
  CMS_text_2->Draw("same");    

  TLatex* lumiText = new TLatex(0.92,0.975,"35.9 fb^{-1} (13 TeV)");
  lumiText->SetNDC();
  lumiText->SetTextFont(42);

//.........这里部分代码省略.........
      if (MCpoissonUncertainty_flag == 1) {

	hMC[j]->SetBinError(i,sqrt(hMC[j]->GetBinContent(i)));
	
      }

    }

    hMC[j]->SetFillColor(histColor[j]);
    stackNorm += hMC[j]->Integral();

  }

  // loop again on MC histograms to scale them and then fill the thstack

  for (Int_t j = 0; j < nFiles; j++) {

    if (data0_noData1 == 0 && MCnormalizedToData_flag != 0) {

      Double_t dataNorm = hdata->Integral();

      if (binDensity_flag != 0) hMC[j]->Scale(dataNorm/stackNorm,"width");
      else hMC[j]->Scale(dataNorm/stackNorm);

    } else if (binDensity_flag != 0) hMC[j]->Scale(1.0,"width");  // option width divides by bin width and manages the correct error setting

    hMCstack->Add(hMC[j]);

  }

  if (data0_noData1 == 0) {

    cout << "Events in data: " << hdata->Integral() << endl;
    cout << "Events in MC   : " << ((TH1D*) hMCstack->GetStack()->Last())->Integral() << endl;

  }
  
  if (data0_noData1 == 0 && binDensity_flag != 0) hdata->Scale(1.0,"width");
  //if (data0_noData1 == 0 && MCnormalizedToData_flag != 0)

  // now here we go with the canvas

  TH1D * ratioplot = NULL; // will use it for the ratio plots

  TPad *subpad_1 = NULL;  // will use it to access specific subpad in canvas
  TPad *subpad_2 = NULL; 

  TCanvas *c;
  if (data0_noData1 == 0) c = new TCanvas(canvasName.c_str(), (var + " distribution").c_str(), 700, 700);
  else c = new TCanvas(canvasName.c_str(), (var + " distribution").c_str());
  TLegend *leg = new TLegend(0.7,0.6,0.99,0.94);  

  // if there are data, split canvas to draw the dta/MC ratio plot

  if (data0_noData1 == 0) {

    subpad_1 = new TPad("pad_1","",0.0,0.28,1.0,1.0);
    if (yAxisLog_flag) subpad_1->SetLogy();
    //subpad_1->SetBottomMargin(0);
    subpad_2 = new TPad("pad_2","",0.0,0.0,1.0,0.32);
    subpad_2->SetGridy();
    //subpad_2->SetTopMargin(0);
    subpad_2->SetBottomMargin(0.3);
    subpad_1->Draw();
    subpad_2->Draw();

////////////////////////////////////////////////////////////////////////////////////////
// Draw and save the canvas
// This bit shouldn't need to be updated except the region and channel labels
// As of Jan 10,2013 - the known issues are:
// ++ Error bands don't show correctly
// ++ Data histogram needs to be converted to TGraphAsymmErrors and 
//    plotted w/ poisson errors
// ++ Ratio plot needs to divide the Data TGraphAsymmErrors and
//    TGraphAsymmErrors for the total SM (by turning off the errors since they are
//    shown explicitly)
////////////////////////////////////////////////////////////////////////////////////////
void PlotUnifier(
                  TH1F* histos[9],                   // Array of pointers to histograms
                  TGraphAsymmErrors* smErrorBand,    // Error band for the top pad
                  TGraphAsymmErrors* ratioErrorBand  // Error band for the bottom pad
                )
{
  ////////////////////////////////////////////////////////////////////////////////////////
  // Set Histogram Properties, build the stack and make the ratio plot
  histos[0]->SetMarkerStyle(20)     ; histos[0]->SetMarkerSize(1.2); 
  histos[1]->SetFillColor(kOrange-2); histos[4]->SetLineWidth (2);
  histos[2]->SetFillColor(kAzure-9 ); histos[1]->SetLineWidth (2);
  histos[3]->SetFillColor(kSpring+1); histos[2]->SetLineWidth (2);
  histos[4]->SetFillColor(kAzure+4 ); histos[3]->SetLineWidth (2);
  histos[5]->SetFillColor(kGray    ); histos[5]->SetLineWidth (2);
  histos[6]->SetFillColor(kYellow-9); histos[6]->SetLineWidth (2);
  histos[7]->SetLineColor(kMagenta ); histos[7]->SetLineWidth (3); histos[7]->SetLineStyle(2);
  histos[8]->SetLineColor(kGreen+1 ); histos[8]->SetLineWidth (3); histos[8]->SetLineStyle(2);
  
  ////////////////////////////////////////////////////////////////////////////////////////
  // Here comes the stack 
  THStack* sm = new THStack("sm","sm");
  sm->Add(histos[6]);
  sm->Add(histos[5]);
  sm->Add(histos[4]);
  sm->Add(histos[3]);
  sm->Add(histos[2]);
  sm->Add(histos[1]);

  ////////////////////////////////////////////////////////////////////////////////////////
  // Define the legend and fill in the stack order
  TLegend* legendA        = new TLegend(0.49,0.65,0.74,0.85);
  TLegend* legendB        = new TLegend(0.62,0.65,0.87,0.85);
  legendA->SetTextFont(42);    
  legendB->SetTextFont(42);    
  legendA->SetTextSize(0.03);
  legendB->SetTextSize(0.03);
  legendA->SetBorderSize(0);
  legendB->SetBorderSize(0);
  legendA->SetFillColor(0);
  legendB->SetFillColor(0);
  legendA->AddEntry(histos[0],"Data"               ,"p");
  legendA->AddEntry(histos[1],"Z+jets"             ,"f");
  legendA->AddEntry(histos[2],"WW"                 ,"f");
  legendA->AddEntry(histos[3],"t#bar{t}+Wt"        ,"f");
  legendA->AddEntry(histos[4],"ZV"                 ,"f");
  legendB->AddEntry(histos[5],"Non-prompt leptons" ,"f");
  legendB->AddEntry(histos[6],"Higgs"              ,"f");
  TH1F* histoError = new TH1F("histoError"         ,"histoError",1,0,1);
  histoError->SetFillStyle(3004);
  histoError->SetFillColor(kBlack);
  legendB->AddEntry(histoError,"Bkg. Uncert."      ,"f");
  legendB->AddEntry(histos[7],"Signal 1"           ,"l");
  legendB->AddEntry(histos[8],"Signal 2"           ,"l");

  ////////////////////////////////////////////////////////////////////////////////////////
  // Ratio plot 
  TH1F* ratio      = (TH1F*) histos[0]->Clone(); ratio->Reset(); 
  TH1F* stackHisto = (TH1F*) sm->GetStack()->Last();
  ratio->Divide(histos[0],stackHisto);
  ratio->SetMarkerSize(1.2);
  ratio->SetMarkerStyle(20);
  ratio->SetLineColor(kBlack);
  ratio->SetLineWidth(2);

  ratio->GetXaxis()->SetTitle("Observable");
  ratio->GetYaxis()->SetTitle("Data/SM");
  ratio->GetXaxis()->SetLabelSize(0.13);
  ratio->GetXaxis()->SetLabelOffset(0.02);
  ratio->GetXaxis()->SetTitleSize(0.14);
  ratio->GetXaxis()->SetTitleOffset(1.2);

  ratio->GetYaxis()->SetRangeUser(0,2);
  ratio->GetYaxis()->SetLabelSize(0.13);
  ratio->GetYaxis()->SetLabelOffset(0.0125);
  ratio->GetYaxis()->SetTitleSize(0.14);
  ratio->GetYaxis()->SetTitleOffset(0.5);
  ratio->GetYaxis()->SetNdivisions(5);

  ////////////////////////////////////////////////////////////////////////////////////////
  // Error bands
  smErrorBand->SetMarkerSize(0); 
  smErrorBand->SetFillStyle(3004);
  smErrorBand->SetFillColor(kBlack);
  ratioErrorBand->SetMarkerSize(0); 
  ratioErrorBand->SetFillStyle(3004);
  ratioErrorBand->SetFillColor(kBlack);
 
  ////////////////////////////////////////////////////////////////////////////////////////
  // Draw and save the canvas
//.........这里部分代码省略.........