C++ TLegend::SetFillColor方法代码示例

//.........这里部分代码省略.........
            if( TMath::Abs(c->akPu3PF.jteta[jet2index]) > 1.3 )
                continue;
            Double_t dphi = TMath::ACos(TMath::Cos(c->photon.phi[leadingIndex] - c->akPu3PF.jtphi[jet2index]));
            if( dphi < 2.7 )
                continue;

            Double_t alpha = jet3pt/leadingPt;
            Double_t ratio = jet2pt/leadingPt;

            int ptbin = returnPtBin(leadingPt);
            if(ptbin != -1)
                Ratio[ii][ptbin]->Fill(alpha,ratio);

            alphas->Fill(alpha);
            hPhotonPt[ii]->Fill(leadingPt);
            hPhotonEta[ii]->Fill( c->photon.eta[leadingIndex] );
            hJet2Pt[ii]->Fill(jet2pt);
            hJet2Eta[ii]->Fill( c->akPu3PF.jteta[jet2index] );
            hJet3Pt[ii]->Fill(jet3pt);
            hJet3Eta[ii]->Fill( c->akPu3PF.jteta[jet3index] );
        }
    }

    TCanvas *canvas[nptbins];
    for(int i = 0; i < nptbins; i++)
    {
        TH1D *hRatiod = Ratio[0][i]->ProjectionX();
        TH1D *hRatiomc = Ratio[1][i]->ProjectionX();
        RdRmcalpha[i]->Divide(hRatiod,hRatiomc,1,1);

        canvas[i] = new TCanvas();
        RdRmcalpha[i]->Draw();
    }

    TCanvas *danvas = new TCanvas();
    alphas->Draw();
    danvas->SaveAs("alpha_counts.gif");


    hPhotonPt[0]->Scale( 1./hPhotonPt[0]->GetEntries() );
    hPhotonPt[1]->Scale( 1./hPhotonPt[1]->GetEntries() );
    TCanvas *d1 = new TCanvas();
    hPhotonPt[0]->Draw();
    hPhotonPt[1]->SetMarkerColor(kRed);
    hPhotonPt[1]->Draw("same");
    d1->SetLogy();
    TLegend *leg = new TLegend(0.75,0.75,0.9,0.9);
    leg->SetFillColor(0);
    leg->AddEntry(hPhotonPt[0], "Data","P");
    leg->AddEntry(hPhotonPt[1], "MC", "P");
    leg->DrawClone();
    d1->SaveAs("photon_pt.gif");

    hPhotonEta[0]->Scale( 1./hPhotonEta[0]->GetEntries() );
    hPhotonEta[1]->Scale( 1./hPhotonEta[1]->GetEntries() );
    TCanvas *d2 = new TCanvas();
    hPhotonEta[0]->Draw();
    hPhotonEta[1]->SetMarkerColor(kRed);
    hPhotonEta[1]->Draw("same");
    leg->DrawClone();
    d2->SaveAs("photon_eta.gif");

    hJet2Pt[0]->Scale( 1./hJet2Pt[0]->GetEntries() );
    hJet2Pt[1]->Scale( 1./hJet2Pt[1]->GetEntries() );
    TCanvas *d3 = new TCanvas();
    hJet2Pt[0]->Draw();
    hJet2Pt[1]->SetMarkerColor(kRed);
    hJet2Pt[1]->Draw("same");
    d3->SetLogy();
    leg->DrawClone();
    d3->SaveAs("jet2_pt.gif");

    hJet2Eta[0]->Scale( 1./hJet2Eta[0]->GetEntries() );
    hJet2Eta[1]->Scale( 1./hJet2Eta[1]->GetEntries() );
    TCanvas *d4 = new TCanvas();
    hJet2Eta[0]->Draw();
    hJet2Eta[1]->SetMarkerColor(kRed);
    hJet2Eta[1]->Draw("same");
    leg->DrawClone();
    d4->SaveAs("jet2_eta.gif");

    hJet3Pt[0]->Scale( 1./hJet3Pt[0]->GetEntries() );
    hJet3Pt[1]->Scale( 1./hJet3Pt[1]->GetEntries() );
    TCanvas *d5 = new TCanvas();
    hJet3Pt[0]->Draw();
    hJet3Pt[1]->SetMarkerColor(kRed);
    hJet3Pt[1]->Draw("same");
    d5->SetLogy();
    leg->DrawClone();
    d5->SaveAs("jet3_pt.gif");

    hJet3Eta[0]->Scale( 1./hJet3Eta[0]->GetEntries() );
    hJet3Eta[1]->Scale( 1./hJet3Eta[1]->GetEntries() );
    TCanvas *d6 = new TCanvas();
    hJet3Eta[0]->Draw();
    hJet3Eta[1]->SetMarkerColor(kRed);
    hJet3Eta[1]->Draw("same");
    leg->DrawClone();
    d6->SaveAs("jet3_eta.gif");
}

void drawSignalBGData( const std::string& dirName, const std::string& channel, const std::string& var, const std::string& varName, const std::string& units, int nBins, float xMin, float xMax, TTree* tree_sig, TTree* tree_bg, TTree* tree_data, const std::string& fullSelection, const std::string& fullSelection_sidebands, const std::string& suffix ) {


  DrawBase* db = new DrawBase("sigbfdata" );

  TH1D* h1_data = new TH1D("data", "", nBins, xMin, xMax ); 
  TH1D* h1_sig = new TH1D("sig", "", nBins, xMin, xMax ); 
  TH1D* h1_bg = new TH1D("bg", "", nBins, xMin, xMax ); 

  h1_data->Sumw2();
  h1_sig->Sumw2();
  h1_bg->Sumw2();

  h1_sig->SetFillColor(46);
  h1_bg->SetFillColor(38);
  h1_sig->SetLineColor(46);
  h1_bg->SetLineColor(38);

  h1_sig->SetFillStyle(3004);
  h1_bg->SetFillStyle(3005);


  if( BLINDED )
    tree_data->Project( "data", var.c_str(), fullSelection_sidebands.c_str() );
  else
    tree_data->Project( "data", var.c_str(), fullSelection.c_str() );
  tree_sig->Project(  "sig",  var.c_str(), fullSelection.c_str() );
  tree_bg ->Project(  "bg",   var.c_str(), fullSelection.c_str() );


  float data_integral = h1_data->Integral();
  if( data_integral < 2 ) data_integral=3.;
  h1_sig->Scale( data_integral/h1_sig->Integral() );
  h1_bg->Scale( data_integral/h1_bg->Integral() );
  


  TGraphAsymmErrors* gr_data = fitTools::getGraphPoissonErrors( h1_data, false );
  gr_data->SetMarkerStyle(20);  
  gr_data->SetMarkerSize(1.3);  
   

  TString category_tstr(channel);
  bool isLeptonic = category_tstr.Contains("leptonic");

  std::string channelName_legend = (isLeptonic) ? "Leptonic Channel" : "Hadronic Channel";


  TLegend* legend = new TLegend( 0.55, 0.66, 0.9, 0.92, channelName_legend.c_str() );
  legend->SetTextFont(42);
  legend->SetFillColor(kWhite);
  legend->SetTextSize(0.038);
  legend->SetTextColor(kBlack);
  if( useCS )
    legend->AddEntry( gr_data, "Data Control Sample", "P" );
  else {
    if( BLINDED ) {
      legend->AddEntry( gr_data, "Sidebands Data", "P" );
    } else {
      legend->AddEntry( gr_data, "Sidebands Data", "P" );
    }
  }
  legend->AddEntry( h1_sig, "tHq (125)", "F" );
  legend->AddEntry( h1_bg, "ttH (125)", "F" );


  float yMax = 0.;
  float yMax_data = h1_data->GetMaximum();
  if( yMax_data>0. ) {
    if( yMax_data < 1.5 ) yMax=3.5;
    else if( yMax_data < 2.5 ) yMax=6.;
    else                       yMax=yMax_data*2.;
  } else {
    float yMax_mc = h1_sig->GetMaximum()*1.2;
    yMax = 1.8*yMax_mc;
  }


  TH2D* h2_axes = new TH2D( "axes", "", 10, xMin, xMax, 10, 0., yMax);
  std::string xTitle = varName;
  if( units!="" ) xTitle = xTitle + " [" + units + "]";
  h2_axes->SetXTitle(xTitle.c_str());
  h2_axes->SetYTitle( "Arbitrary Units" );

  TCanvas* c1 = new TCanvas("c1", "", 600, 600);
  c1->cd();


  TPaveText* labelTop = db->get_labelTop();

  h2_axes->Draw();
  legend->Draw("same");   
  labelTop->Draw("same");
  h1_bg->Draw("h same");
  h1_sig->Draw("h same");
  gr_data->Draw("p same");


  gPad->RedrawAxis();

//.........这里部分代码省略.........

void DrawNCLUSTER()
{
  //=========Macro generated from canvas: r1/Energy
  //=========  (Fri Apr  6 17:14:48 2012) by ROOT version5.28/00g
  TCanvas *r1 = new TCanvas("r1", "Ncluster",12,24,550,741);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetFrameLineWidth(1);
  r1->Range(0,0,1,1);
  r1->SetFillColor(0);
  r1->SetBorderMode(0);
  r1->SetBorderSize(0);
  r1->SetTickx(1);
  r1->SetTicky(1);
  r1->SetLeftMargin(0.16);
  r1->SetRightMargin(0.01);
  r1->SetTopMargin(0.0256917);
  r1->SetBottomMargin(0.07692308);
  r1->SetFrameBorderMode();
  
  // ------------>Primitives in pad: r1_1
  TPad *r1_1 = new TPad("r1_1", "Energy_1",0.02,0.37,0.95,0.99);
  r1_1->Draw();
  r1_1->cd();
  r1_1->Range(-19,0.01,95,95);
  r1_1->SetFillColor(0);
  r1_1->SetBorderMode(0);
  r1_1->SetBorderSize(2);
  r1_1->SetTickx(1);
  r1_1->SetTicky(1);
  r1_1->SetLeftMargin(0.16);
  r1_1->SetRightMargin(0.01);
  r1_1->SetTopMargin(0.02);
  r1_1->SetBottomMargin(0);
  r1_1->SetFrameBorderMode(0);
  r1_1->SetFrameBorderMode(0);
   
  TH1D *he = new TH1D("he"," ",55,0,55);
  he->SetMinimum(0.01);
  he->SetMaximum(70);
  he->SetStats(0);
  he->GetXaxis()->SetTitle("E_{beam} [GeV]");
  he->GetXaxis()->SetLabelFont(43);
  he->GetXaxis()->SetLabelSize(0);
  he->GetXaxis()->SetTitleFont(43);
  he->GetXaxis()->SetTitleSize(0); 
  he->GetYaxis()->SetTitle("<N_{cluster}>");
  he->GetYaxis()->SetLabelFont(43);
  he->GetYaxis()->SetTitleSize(30);
  he->GetYaxis()->SetLabelSize(20);
  he->GetYaxis()->SetTitleFont(43);
  he->GetYaxis()->SetTitleOffset(1.7);
  he->GetZaxis()->SetLabelFont(42);
  he->GetZaxis()->SetTitleSize(0.05);
  he->GetZaxis()->SetTitleFont(42);
  he->Draw("");
   
  std::vector<result> resultData=readTXT(std::string("ncluster_augsep2012_electron.txt"));
  TGraphErrors *gre = new TGraphErrors(resultData.size());
  gre->SetName("Graph");
  gre->SetTitle("Graph");
  gre->SetLineColor(1);
  gre->SetFillStyle(1);
  gre->SetFillColor(1);
  gre->SetLineWidth(2);
  gre->SetMarkerColor(1);
  gre->SetMarkerStyle(34);
  gre->SetMarkerSize(1.2);
  for(unsigned int i=0; i<resultData.size(); i++){
    gre->SetPoint(i,resultData.at(i).ebeam,resultData.at(i).ncluster);
    gre->SetPointError(i,0,resultData.at(i).nclusterError);
  }

  TH1F *Graph_Graph3 = new TH1F("Graph_Graph3","Graph",100,0,87.5);
  Graph_Graph3->SetMinimum(0);
  Graph_Graph3->SetMaximum(1193.483);
  Graph_Graph3->SetDirectory(0);
  Graph_Graph3->SetStats(0);
  Graph_Graph3->GetXaxis()->SetLabelFont(42);
  Graph_Graph3->GetXaxis()->SetTitleSize(0.05);
  Graph_Graph3->GetXaxis()->SetTitleFont(42);
  Graph_Graph3->GetYaxis()->SetLabelFont(42);
  Graph_Graph3->GetYaxis()->SetTitleSize(0.05);
  Graph_Graph3->GetYaxis()->SetTitleFont(42);
  Graph_Graph3->GetZaxis()->SetLabelFont(42);
  Graph_Graph3->GetZaxis()->SetTitleSize(0.05);
  Graph_Graph3->GetZaxis()->SetTitleFont(42);
  gre->SetHistogram(Graph_Graph3);
   
  gre->Draw("p");
   
  std::vector<result> resultFTFP=readTXT(std::string("ncluster_ftfp_bert_hp_electron.txt"));
  gre = new TGraphErrors(resultFTFP.size());
  gre->SetName("Graph");
  gre->SetTitle("Graph");
  gre->SetFillColor(kRed-3);
  gre->SetMarkerColor(kRed-3);
  gre->SetLineWidth(2);
  gre->SetMarkerStyle(20);
  for(unsigned int i=0; i<resultFTFP.size(); i++){
//.........这里部分代码省略.........

//.........这里部分代码省略.........

  //TCut njets    = "njets>=2";
  TCut tkisoold = "tkisoold/probe->pt()>0.1";
  TCut tkisonew = "tkisonew/probe->pt()>0.1";

  //-----------------------------------------
  // check nvtx data vs. MC
  //-----------------------------------------

  TH1F *hnvtx_mc   = new TH1F("hnvtx_mc"  ,"",30,0,30);
  TH1F *hnvtx_data = new TH1F("hnvtx_data","",30,0,30);

  hnvtx_mc->Sumw2();
  hnvtx_data->Sumw2();

  chdata->Draw("nvtx>>hnvtx_data",(eltnpcut||mutnpcut));
  chmc->Draw("nvtx>>hnvtx_mc",(eltnpcut||mutnpcut)*vtxweight);

  TCanvas *c1 = new TCanvas();
  c1->cd();
  
  hnvtx_mc->SetLineColor(2);
  hnvtx_mc->SetMarkerColor(2);
  hnvtx_data->SetLineColor(4);
  hnvtx_data->SetMarkerColor(4);
  
  hnvtx_data->GetXaxis()->SetTitle("N_{VTX}");
  hnvtx_data->DrawNormalized();
  hnvtx_mc->DrawNormalized("same");

  TLegend *leg = new TLegend(0.6,0.6,0.8,0.8);
  leg->AddEntry(hnvtx_data,"data","lp");
  leg->AddEntry(hnvtx_mc,"MC","lp");
  leg->SetFillColor(0);
  leg->SetBorderSize(0);
  leg->Draw();

  if( printplot ) c1->Print("plots/nvtx.pdf");

  //--------------------------
  // absolute track isolation
  //--------------------------
  
  bool residual = true;
  bool log      = true;
  
  char* var = "tkisonewnoveto";
  
  plotDistribution( chdata , chmc , TCut(eltnpcut)        , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "el_tkiso_0j" , printplot , residual , log );
  plotDistribution( chdata , chmc , TCut(mutnpcut)        , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "mu_tkiso_0j" , printplot , residual , log );

  plotDistribution( chdata , chmc , TCut(eltnpcut+njets1) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "el_tkiso_1j" , printplot , residual , log );
  plotDistribution( chdata , chmc , TCut(mutnpcut+njets1) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "mu_tkiso_1j" , printplot , residual , log );
  
  plotDistribution( chdata , chmc , TCut(eltnpcut+njets2) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "el_tkiso_2j" , printplot , residual , log );
  plotDistribution( chdata , chmc , TCut(mutnpcut+njets2) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "mu_tkiso_2j" , printplot , residual , log );

  plotDistribution( chdata , chmc , TCut(eltnpcut+njets3) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "el_tkiso_3j" , printplot , residual , log );
  plotDistribution( chdata , chmc , TCut(mutnpcut+njets3) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "mu_tkiso_3j" , printplot , residual , log );

  plotDistribution( chdata , chmc , TCut(eltnpcut+njets4) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "el_tkiso_4j" , printplot , residual , log );
  plotDistribution( chdata , chmc , TCut(mutnpcut+njets4) , vtxweight , var , 10 , 0 , 10 , "abs tkiso [GeV]"    , "mu_tkiso_4j" , printplot , residual , log );

  //--------------------------
  // relative track isolation
  //--------------------------

void compareSystVariable( const std::string& dirName, const std::string& channel, const std::string& dataset, TTree* tree_central, TFile* file_systUp, TFile* file_systDown, const std::string& fullSelection, const std::string& systName, const std::string& varName, const std::string& axisTitle, const std::string& units, int nBins, float xMin, float xMax, const std::string& suffix) {


  DrawBase* db = new DrawBase("syst" );

  TTree* tree_systUp = (TTree*)file_systUp->Get(dataset.c_str());
  TTree* tree_systDown = (file_systDown!=0) ? (TTree*)file_systDown->Get(dataset.c_str()) : 0;

  TH1D* h1_central = new TH1D("central", "", nBins, xMin, xMax);
  TH1D* h1_systUp = new TH1D("systUp", "", nBins, xMin, xMax);
  TH1D* h1_systDown = new TH1D("systDown", "", nBins, xMin, xMax);

  h1_central->Sumw2();
  h1_systUp->Sumw2();
  h1_systDown->Sumw2();

  tree_central->Project( "central", varName.c_str(), fullSelection.c_str() );
  tree_systUp->Project( "systUp", varName.c_str(), fullSelection.c_str() );
  if( tree_systDown!=0 )
    tree_systDown->Project( "systDown", varName.c_str(), fullSelection.c_str() );

  h1_central ->SetLineColor(kBlack);
  h1_systUp  ->SetLineColor(kGreen);
  h1_systDown->SetLineColor(kRed);

  h1_central ->SetLineWidth(2);
  h1_systUp  ->SetLineWidth(2);
  h1_systDown->SetLineWidth(2);

  //h1_systUp  ->SetLineStyle(2);
  //h1_systDown->SetLineStyle(2);


  TLegend* legend = new TLegend( 0.6, 0.66, 0.9, 0.92, dataset.c_str() );
  legend->SetTextFont(42);
  legend->SetFillColor(kWhite);
  legend->SetTextSize(0.038);
  legend->SetTextColor(kBlack);
  legend->AddEntry( h1_central, "Central", "L" );
  std::string systUp_text = systName + " +1#sigma";
  std::string systDown_text = systName + " -1#sigma";
  legend->AddEntry( h1_systUp, systUp_text.c_str(), "L" );
  if( tree_systDown!=0 )
    legend->AddEntry( h1_systDown, systDown_text.c_str(), "L" );


  float yMax = 1.6*h1_central->GetMaximum();


  TH2D* h2_axes = new TH2D( "axes", "", 10, xMin, xMax, 10, 0., yMax);
  std::string xTitle = axisTitle;
  if( units!="" ) xTitle = xTitle + " [" + units + "]";
  h2_axes->SetXTitle(xTitle.c_str());
  h2_axes->SetYTitle( "Events" );

  TCanvas* c1 = new TCanvas("c1", "", 600, 600);
  c1->cd();


  TPaveText* labelTop = db->get_labelTop();

  h2_axes->Draw();
  legend->Draw("same");   
  labelTop->Draw("same");
  h1_central->Draw("h same");
  h1_systUp->Draw("h same");
  if( tree_systDown!=0 )
    h1_systDown->Draw("h same");


  gPad->RedrawAxis();

  std::string canvasName = dirName + "/" + channel + "/syst" + systName + "_" + varName + "_" + dataset + "_" + channel;
  if( suffix!="" ) canvasName = canvasName + "_" + suffix;
  std::string canvasName_eps = canvasName + ".eps";
  std::string canvasName_png = canvasName + ".png";

  c1->SaveAs(canvasName_eps.c_str());
  c1->SaveAs(canvasName_png.c_str());


  delete db;

  delete h1_central;
  delete h1_systUp;
  delete h1_systDown;

  delete h2_axes;

  delete legend;
  delete c1;



}

void vs_PlotQCDcomp() {

  Bool_t saveC  = false;
  Bool_t diJets = true;
  Bool_t isMC   = false;

  TString vsSave;
  vsSave = "_T07w";

  TString jSave;
  if(diJets) jSave = "2j";
  else       jSave = "3j";

  TString sample;
  if(diJets) sample = "select_1ph_2jets";
  else       sample = "select_1ph_3jets";

  TString FOtag;
  if(isMC) FOtag = "_FO_CorrMC";
  //  else     FOtag = "_FO_Corr";
  else     FOtag = "_FO";

  TString sysu = "SYSTUP_";
  TString sysd = "SYSTDOWN_";

  TString outDir;
  outDir = "Plots_PhotonSusyAnalysis/QCDweights/";


  //  string inputFile1 = "/data/user/vsola/CMSSW_Releases/CMSSW_5_3_9/src/Plots_PhotonSusyAnalysis/Merged_QCD_PhotonJet_T07w_PAT/mergedHistos.root";
  //  string inputFile1 = "/data/user/vsola/CMSSW_Releases/CMSSW_5_3_9/src/Plots_PhotonSusyAnalysis/Merged_Data_V05w_PAT/mergedHistos.root";
  string inputFile1 = "/data/user/vsola/CMSSW_Releases/CMSSW_5_3_9/src/Plots_PhotonSusyAnalysis/PhotonHadReReco_22Jan2013_V05_PAT/mergedHistos.root";
  
  
  setMyTDRStyle();
  gROOT->SetStyle("mytdrStyle");
  
  //  gStyle->SetHistMinimumZero();
  //  gStyle->SetPaintTextFormat("4.2f");

  gStyle->SetHistLineWidth(2);
  gStyle->UseCurrentStyle();

  //  gStyle->SetPadTopMargin(1.0);
  //  gStyle->SetPadLeftMargin(3.2);
  //  gStyle->SetPadRightMargin(4.5);
  //  gStyle->SetPadBottomMargin(3.2);

  gROOT->ForceStyle(1);


  static const Int_t nHt  =  7;
  static const Int_t nPt  = 13;
  static const Int_t nMet = 16;
  static const Int_t nJet = 15;

  Double_t binPt[nPt+1]   = {0.,75.,90.,120.,160.,210.,260.,320.,400.,500.,650.,800.,1000.,1500.};
  Double_t binHt[nHt+1]   = {0.,400.,450.,550.,700.,900.,1200.,1500.};
  Double_t binMet[nMet+1] = {0.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,120.,160.,200.,270.,350.,500.};
  Double_t binJet[nJet+1] = {0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15};


  TLatex *as = new TLatex();
  as->SetNDC(true);
  as->SetTextColor(12);
  as->SetTextFont(43);
  as->SetTextSize(19);

  TLegend *legend = new TLegend(0.60, 0.70, 0.75, 0.85, "");
  legend->SetFillColor(10);
  legend->SetFillStyle(1001);
  legend->SetTextSize(0.04);
  legend->SetBorderSize(0);
  legend->SetShadowColor(0);

  std::string outLumi  = "CMS Work in Progress - QCD MC #sqrt{s} = 8 TeV - #geq 1 #gamma, #geq 2 j";

  TFile* f1 = TFile::Open( inputFile1.c_str() );

  //photon Pt
  TCanvas * cPt = new TCanvas("cPt","cPt");

  TH1F* hpt   = (TH1F*) f1->Get( sample+"/PreselCut_photonPt" );
  TH1F* hptFO = (TH1F*) f1->Get( sample+FOtag+"/PreselCut_photonPt" );

  TH1F *hptR = (TH1F*) hpt->Rebin(nPt,hpt->GetTitle(),binPt);
  if (hptR->GetSumw2N() == 0)
    hptR->Sumw2();

  TH1F *hptFOR = (TH1F*) hptFO->Rebin(nPt,hptFO->GetTitle(),binPt);
  if (hptFOR->GetSumw2N() == 0)
    hptFOR->Sumw2();

  TH1F* hptSU = (TH1F*) f1->Get( sample+FOtag+"/PreselCut_"+sysu+"photonPt" );
  TH1F* hptSD = (TH1F*) f1->Get( sample+FOtag+"/PreselCut_"+sysd+"photonPt" );

  TH1F *hptSUR = (TH1F*) hptSU->Rebin(nPt,hptSU->GetTitle(),binPt);
  if (hptSUR->GetSumw2N() == 0)
    hptSUR->Sumw2();

//.........这里部分代码省略.........

//.........这里部分代码省略.........
	
	mu_orig_import->Draw();
	
	TH1F* mu_orig_mine = new TH1F("mu_orig_mine","#mu distribution before reweighting;#mu;fraction of events",21,-0.5,20.5);
	TH1F* mu_RW_mine = new TH1F("mu_RW_mine","#mu distribution after reweighting;#mu;fraction of events",21,-0.5,20.5);
	TH1F* mu_data_mine = new TH1F("mu_data_mine","#mu distribution in data;#mu;fraction of events",21,-0.5,20.5);
	
	
	for(int iBin = 1; iBin <= mu_orig_mine->GetNbinsX(); iBin++)
	{
		mu_orig_mine->SetBinContent(iBin,mu_orig_import->GetBinContent(iBin));
		mu_RW_mine->SetBinContent(iBin,mu_RW_import->GetBinContent((iBin-1)*5+1));
		
		float dataMuLowEdge  = mu_orig_mine->GetBinLowEdge(iBin);
		float dataMuHighEdge = mu_orig_mine->GetBinLowEdge(iBin) + mu_orig_mine->GetBinWidth(iBin);
		int dataBinNumberLowEdge  = mu_data->GetXaxis()->FindFixBin(dataMuLowEdge);
		int dataBinNumberHighEdge = mu_data->GetXaxis()->FindFixBin(dataMuHighEdge);
		float dataBinContent = mu_data->Integral( dataBinNumberLowEdge, dataBinNumberHighEdge-1 );
		
		mu_data_mine->SetBinContent(iBin,dataBinContent);
		
		
	}
	mu_orig_mine->SetBinContent(1,mu_orig_import->GetBinContent(101));

	mu_orig_import->Scale(1./mu_orig_import->GetSumOfWeights());
	mu_RW_import->Scale(1./mu_RW_import->GetSumOfWeights());
	mu_data->Scale(1./mu_data->GetSumOfWeights());
	//mu_data->Scale(10./mu_data->GetSumOfWeights());
	
	
	TCanvas* c1 = new TCanvas("c1","c1", 800,600);
     
	mu_orig_import->Draw();
	
	c1->Print("mu_mc.png");
	c1->Print("mu_mc.eps");
	
	mu_data->GetXaxis()->SetTitle("#mu");
	mu_data->GetYaxis()->SetTitle("fraction of events");
	
	mu_data->Draw();
	
	c1->Print("mu_data.png");
	c1->Print("mu_data.eps");
	
	
	float myMax[] = {mu_RW_import->GetMaximum(),mu_data->GetMaximum()};
	
	float Max_h = *max_element(myMax,myMax+2);
	
	mu_RW_import->SetMaximum(1.2*Max_h);
	mu_RW_import->SetMarkerStyle(25);
	mu_RW_import->SetMarkerSize(1.5);
	mu_RW_import->GetXaxis()->SetTitle("True number of interactions per event");
	mu_RW_import->GetXaxis()->SetRangeUser(0,60);
	mu_RW_import->GetYaxis()->SetTitle("Fraction of events");
	
	mu_RW_import->Draw("phist");
	//~ 
	//~ //mu_data->SetLineColor(kRed);
	//~ //mu_data->Draw("same");
	mu_data->SetLineColor(kGreen);
	mu_data->SetMarkerStyle(22);
	mu_data->SetMarkerSize(1.5);
	
	mu_data->SetMarkerColor(kGreen);
	
	mu_data->Draw("phistsame");
	mu_data->Draw("same");
	mu_RW_import->Draw("same");
	mu_orig_import->SetLineColor(kRed);
	
	mu_orig_import->Draw("histsame");
	
	gPad->RedrawAxis();
	
	TLegend* leg = new TLegend(0.6,0.7,0.8,0.93);
	leg->SetBorderSize(0);
	leg->SetFillColor(0);
	leg->SetTextSize(legend_size);
	leg->AddEntry(mu_RW_import,"re-weighted MC","p");
	leg->AddEntry(mu_data,"data 2012ABCD","lp");
	leg->AddEntry(mu_orig_import,"Summer12 53x MC","l");
	
	leg->Draw();
	
	c1->Print("mu_reweight.png");
	c1->Print("mu_reweight.eps");
	
	c1->SetLogy();
	mu_RW_import->SetMaximum(1200*Max_h);
	
	c1->Print("mu_reweight_log.png");
	c1->Print("mu_reweight_log.eps");
	
	//~ 

	return 0;
}

//.........这里部分代码省略.........


             }//endcap
             

         
         }//trigger pass 

         }//Filling Yield(N_Obsereved) in Data




     

     }//Data 


     //trigger efficiency
     



     ///////////////////////////defining legend
    char Legname1[100];
    TLegend *leg[24];
    for(int k0=0;k0<24;k0++){
    sprintf(Legname1,"leg_1D%i",k0);
    leg[k0]=new TLegend(0.5,0.7,0.80,0.89);
    leg[k0]->SetTextFont(62);
    leg[k0]->SetLineColor(1);
    leg[k0]->SetLineStyle(1);
    leg[k0]->SetLineWidth(3);
    leg[k0]->SetFillColor(0);
    leg[k0]->SetFillStyle(1001);
    leg[k0]->SetShadowColor(0);
    leg[k0]->SetDrawOption(0);
    leg[k0]->SetBorderSize(0);
    leg[k0]->SetTextSize(0.03);
    }


     ///////////////////////////////ZinvMC Yield

     //TCanvas *ZinvMCYield=new TCanvas("ZinvMC","ZinvMC");
     //ZinvMCYield->SetLogy();
     h_NZinv18bin->Sumw2();
     h_NZinv18bin->SetFillColor(1);
     h_NZinv18bin->SetFillStyle(1000);
     h_NZinv18bin->GetXaxis()->SetTitle("bin Number");
     h_NZinv18bin->GetYaxis()->SetTitle("Z/Gamma Ratio");
     //h_NZinv18bin->Draw("hist");

   //////////////////////////////////////////////Calculating Zgamma Ratio start
      TPaveText *tpav1 = new TPaveText(0.1956522,0.6247818,0.729097,0.8970332,"brNDC");

    tpav1->SetBorderSize(0);
    tpav1->SetFillStyle(0);
    tpav1->SetTextAlign(11);
    tpav1->SetTextFont(42);
    tpav1->SetTextSize(0.04);
    tpav1->AddText("HT >500");
    tpav1->AddText("#gamma p_{T} > 100 ");
    tpav1->AddText("NJets >=4");
    tpav1->AddText("MHT>200");
    tpav1->AddText("Btags=0");

void drawSingleLimitPlot( const ZGConfig& cfg, const std::string& limitsFile, float eff, const std::string& cat, const std::string& width, bool onlyExpected ) {


  std::string axisName;
  if( cat=="fit_v0") 
    axisName = "95\% CL UL on #sigma #times BR(Q#rightarrowq#gamma) [fb]";
  else if (cat=="fit_em")
    axisName = "95\% CL UL on #sigma #times BR(Q#rightarrowqZ) [fb]";
  else {
    std::cout << "Unknown category " << cat << " !!! Exiting!" << std::endl;
    exit(1);
  }


  TGraph* gr_obs = new TGraph(0);
  TGraph* gr_exp = new TGraph(0);
  TGraphAsymmErrors* gr_exp_1sigma = new TGraphAsymmErrors(0);
  TGraphAsymmErrors* gr_exp_2sigma = new TGraphAsymmErrors(0);

  std::ifstream ifs(limitsFile.c_str());
  std::cout << "-> Opened file: " << limitsFile << std::endl;
  int iPointExp = 0;
  int iPointObs = 0;
  float lastMass = -1;
  float lastObs = -1;

  while( ifs.good() ) {

    float m, obs, exp, exp_m1s, exp_m2s, exp_p1s, exp_p2s;
    std::string s_m, s_obs, s_exp, s_exp_m1s, s_exp_m2s, s_exp_p1s, s_exp_p2s;
    ifs >> s_m >> m >> s_obs >> obs >> s_exp >> exp >> s_exp_m1s >> exp_m1s >> s_exp_m2s >> exp_m2s >> s_exp_p1s >>  exp_p1s >> s_exp_p2s >> exp_p2s;
    TString m_tstr(s_m);
    if( m_tstr.BeginsWith("#") ) continue;
    if( m==lastMass ) continue;


    float conversion = eff*cfg.lumi();

    obs    /=conversion;
    exp    /=conversion;
    exp_m1s/=conversion;
    exp_m2s/=conversion;
    exp_p1s/=conversion;
    exp_p2s/=conversion;
    std::cout << "m: " << m << " obs: " << obs << " exp: " << exp << " exp_m1s: " << exp_m1s << " exp_m2s: " << exp_m2s << " exp_p1s: " << exp_p1s << " exp_p2s: " << exp_p2s << std::endl;


    bool okForObs = true;
    if( obs*conversion>0.1 && okForObs ) {
      gr_obs       ->SetPoint( iPointObs, m, obs );
      iPointObs++;
    }

    bool okForExp = true;
    if( okForExp ) {
      gr_exp       ->SetPoint( iPointExp, m, exp );
      gr_exp_1sigma->SetPoint( iPointExp, m, exp );
      gr_exp_2sigma->SetPoint( iPointExp, m, exp );
      gr_exp_1sigma->SetPointError( iPointExp, 0., 0., exp-exp_m1s, exp_p1s-exp );
      gr_exp_2sigma->SetPointError( iPointExp, 0., 0., exp-exp_m2s, exp_p2s-exp );
      iPointExp++;
    }

    lastMass = m;
    lastObs = obs;

  }

  gr_obs->SetLineWidth(2);

  gr_exp->SetLineWidth(2);
  gr_exp->SetLineStyle(2);

  gr_exp_1sigma->SetLineWidth(0);
  gr_exp_1sigma->SetFillColor(8);

  gr_exp_2sigma->SetLineWidth(0);
  gr_exp_2sigma->SetFillColor(219);

  TCanvas* c1 = new TCanvas( "c1", "", 600, 600 );
  c1->cd();

  float yMax = (cat=="fit_v0") ? 100. : 200.;

  TH2D* h2_axes = new TH2D("axes", "", 10, 400., 1000., 10, 0., yMax );
  h2_axes->SetYTitle( axisName.c_str() );
  //h2_axes->SetYTitle( "95\% CL UL on #sigma #times BR(A#rightarrowZ#gamma#rightarrowl^{+}l^{-}#gamma) [fb]");
  h2_axes->SetXTitle( "Resonance Mass [GeV]");
  h2_axes->Draw();


  gr_exp_2sigma->Draw("E3 same");
  gr_exp_1sigma->Draw("E3 same");
  gr_exp       ->Draw("L  same");
  if( !onlyExpected )
    gr_obs       ->Draw("L  same");


  gr_exp_1sigma->SetLineWidth(2);
  gr_exp_1sigma->SetLineStyle(2);
//.........这里部分代码省略.........

void phimetphimu2b3()
{
//=========Macro generated from canvas: b3/phimetphimu2b3/b3/phimetphimu2b3
//=========  (Tue Nov 24 14:24:06 2015) by ROOT version6.02/05
   TCanvas *b3/phimetphimu2b3 = new TCanvas("b3/phimetphimu2b3", "b3/phimetphimu2b3",380,402,700,700);
   b3/phimetphimu2b3->Range(0,0,1,1);
   b3/phimetphimu2b3->SetFillColor(0);
   b3/phimetphimu2b3->SetBorderMode(0);
   b3/phimetphimu2b3->SetBorderSize(2);
   b3/phimetphimu2b3->SetLogy();
   b3/phimetphimu2b3->SetFrameBorderMode(0);
   b3/phimetphimu2b3->SetFrameBorderMode(0);
   
   THStack *b3/phimetphimu2b3 = new THStack();
   b3/phimetphimu2b3->SetName("b3/phimetphimu2b3");
   b3/phimetphimu2b3->SetTitle("b3/phimetphimu2b3");
   
   TH1F *b3/phimetphimu2b3_stack_20 = new TH1F("b3/phimetphimu2b3_stack_20","b3/phimetphimu2b3",20,-3.141593,3.141593);
   b3/phimetphimu2b3_stack_20->SetMinimum(-0);
   b3/phimetphimu2b3_stack_20->SetMaximum(-nan);
   b3/phimetphimu2b3_stack_20->SetDirectory(0);
   b3/phimetphimu2b3_stack_20->SetStats(0);

   Int_t ci;      // for color index setting
   TColor *color; // for color definition with alpha
   ci = TColor::GetColor("#000099");
   b3/phimetphimu2b3_stack_20->SetLineColor(ci);
   b3/phimetphimu2b3_stack_20->GetXaxis()->SetTitle("b3/phimetphimu2b3");
   b3/phimetphimu2b3_stack_20->GetXaxis()->SetLabelFont(42);
   b3/phimetphimu2b3_stack_20->GetXaxis()->SetLabelSize(0.035);
   b3/phimetphimu2b3_stack_20->GetXaxis()->SetTitleSize(0.035);
   b3/phimetphimu2b3_stack_20->GetXaxis()->SetTitleFont(42);
   b3/phimetphimu2b3_stack_20->GetYaxis()->SetTitle("Events/pb");
   b3/phimetphimu2b3_stack_20->GetYaxis()->SetLabelFont(42);
   b3/phimetphimu2b3_stack_20->GetYaxis()->SetLabelSize(0.035);
   b3/phimetphimu2b3_stack_20->GetYaxis()->SetTitleSize(0.035);
   b3/phimetphimu2b3_stack_20->GetYaxis()->SetTitleFont(42);
   b3/phimetphimu2b3_stack_20->GetZaxis()->SetLabelFont(42);
   b3/phimetphimu2b3_stack_20->GetZaxis()->SetLabelSize(0.035);
   b3/phimetphimu2b3_stack_20->GetZaxis()->SetTitleSize(0.035);
   b3/phimetphimu2b3_stack_20->GetZaxis()->SetTitleFont(42);
   b3/phimetphimu2b3->SetHistogram(b3/phimetphimu2b3_stack_20);
   
   
   TH1D *phimetphimu2b396 = new TH1D("phimetphimu2b396"," #Delta_{#phi}[#mu2,MET]",20,-3.141593,3.141593);

   ci = TColor::GetColor("#00cc00");
   phimetphimu2b396->SetFillColor(ci);

   ci = TColor::GetColor("#00cc00");
   phimetphimu2b396->SetLineColor(ci);

   ci = TColor::GetColor("#00cc00");
   phimetphimu2b396->SetMarkerColor(ci);
   phimetphimu2b396->SetMarkerStyle(22);
   phimetphimu2b396->GetXaxis()->SetTitle("phimetphimu2b3");
   phimetphimu2b396->GetXaxis()->SetLabelFont(42);
   phimetphimu2b396->GetXaxis()->SetLabelSize(0.035);
   phimetphimu2b396->GetXaxis()->SetTitleSize(0.035);
   phimetphimu2b396->GetXaxis()->SetTitleFont(42);
   phimetphimu2b396->GetYaxis()->SetTitle("Events/pb");
   phimetphimu2b396->GetYaxis()->SetLabelFont(42);
   phimetphimu2b396->GetYaxis()->SetLabelSize(0.035);
   phimetphimu2b396->GetYaxis()->SetTitleSize(0.035);
   phimetphimu2b396->GetYaxis()->SetTitleFont(42);
   phimetphimu2b396->GetZaxis()->SetLabelFont(42);
   phimetphimu2b396->GetZaxis()->SetLabelSize(0.035);
   phimetphimu2b396->GetZaxis()->SetTitleSize(0.035);
   phimetphimu2b396->GetZaxis()->SetTitleFont(42);
   b3/phimetphimu2b3->Add(phimetphimu2b3,"");
   
   TH1D *phimetphimu2b397 = new TH1D("phimetphimu2b397"," #Delta_{#phi}[#mu2,MET]",20,-3.141593,3.141593);

   ci = TColor::GetColor("#00ffff");
   phimetphimu2b397->SetFillColor(ci);

   ci = TColor::GetColor("#00ffff");
   phimetphimu2b397->SetLineColor(ci);

   ci = TColor::GetColor("#00ffff");
   phimetphimu2b397->SetMarkerColor(ci);
   phimetphimu2b397->SetMarkerStyle(20);
   phimetphimu2b397->GetXaxis()->SetTitle("phimetphimu2b3");
   phimetphimu2b397->GetXaxis()->SetLabelFont(42);
   phimetphimu2b397->GetXaxis()->SetLabelSize(0.035);
   phimetphimu2b397->GetXaxis()->SetTitleSize(0.035);
   phimetphimu2b397->GetXaxis()->SetTitleFont(42);
   phimetphimu2b397->GetYaxis()->SetTitle("Events/pb");
   phimetphimu2b397->GetYaxis()->SetLabelFont(42);
   phimetphimu2b397->GetYaxis()->SetLabelSize(0.035);
   phimetphimu2b397->GetYaxis()->SetTitleSize(0.035);
   phimetphimu2b397->GetYaxis()->SetTitleFont(42);
   phimetphimu2b397->GetZaxis()->SetLabelFont(42);
   phimetphimu2b397->GetZaxis()->SetLabelSize(0.035);
   phimetphimu2b397->GetZaxis()->SetTitleSize(0.035);
   phimetphimu2b397->GetZaxis()->SetTitleFont(42);
   b3/phimetphimu2b3->Add(phimetphimu2b3,"");
   
   TH1D *phimetphimu2b398 = new TH1D("phimetphimu2b398"," #Delta_{#phi}[#mu2,MET]",20,-3.141593,3.141593);

//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TGraphAsymmErrors *grmedian_cls = new TGraphAsymmErrors(nMassEff, mass, medianD);
  grmedian_cls->SetName("LimitExpectedCLs");
  TGraphAsymmErrors *gr68_cls = new TGraphAsymmErrors(nMassEff, mass, medianD, 0, 0, down68err, up68err);
  gr68_cls->SetName("Limit68CLs");
  TGraphAsymmErrors *gr95_cls = new TGraphAsymmErrors(nMassEff, mass, medianD, 0, 0, down95err, up95err);
  gr95_cls->SetName("Limit95CLs");

  // TGraphAsymmErrors *grthSM=new TGraphAsymmErrors(nMassEff1,mass1,xs,0,0,0,0);//xs_downerr,xs_uperr);
  TGraph *grthSM=new TGraph(nMassEff,mass,xs);//xs_downerr,xs_uperr);
  grthSM->SetName("SMXSection");


  // TGraphAsymmErrors *grthSM10=new TGraphAsymmErrors(nMassEff1,mass1,xs10,0,0,0,0);
  TGraph *grthSM10=new TGraph(nMassEff,mass,xs10);
  grthSM10->SetName("SMXSection_2nd");

  // double fr_left = 590.0, fr_down = 1E-5, fr_right = 2000.0, fr_up = 0.5; 
   double fr_left = 590.0, fr_down = 0.00001, fr_right = 2500.0, fr_up = 5;

  TCanvas *cMCMC = new TCanvas("c_lim_Asymptotic", "canvas with limits for Asymptotic CLs", 630, 600);
  cMCMC->cd();
  cMCMC->SetGridx(1);
  cMCMC->SetGridy(1);
  // draw a frame to define the range

  TH1F *hr = cMCMC->DrawFrame(fr_left, fr_down, fr_right, fr_up, "");
  TString VV = "ZH";
  
  hr->SetXTitle("M_{Zp} [GeV]");
  hr->SetYTitle("#sigma_{95%} [pb]"); // #rightarrow 2l2q
  hr->SetMinimum(0.0001);
  hr->SetMaximum(100.0);

  gr95_cls->SetFillColor(kYellow);
  gr95_cls->SetFillStyle(1001);//solid
  gr95_cls->SetLineStyle(kDashed);
  gr95_cls->SetLineWidth(3);
  gr95_cls->GetXaxis()->SetTitle("M_{V'} [GeV]");
  gr95_cls->GetYaxis()->SetTitle("#sigma_{95%} #times BR(V' #rightarrow " + VV + ") [pb]"); // #rightarrow 2l2q
  gr95_cls->GetXaxis()->SetRangeUser(fr_left, fr_right);

  gr95_cls->Draw("3");
  //  gr95_cls->SetMinimum(0.00001);
  //gr95_cls->SetMaximum(1000.0);
  
  //grmedian_cls->SetMinimum(0.00001);
  //grmedian_cls->SetMaximum(1000.0);
  
  gr68_cls->SetFillColor(kGreen);
  gr68_cls->SetFillStyle(1001);//solid
  gr68_cls->SetLineStyle(kDashed);
  gr68_cls->SetLineWidth(3);
  gr68_cls->Draw("3same");
  grmedian_cls->GetXaxis()->SetTitle("M_{V'} [GeV]");
  grmedian_cls->GetYaxis()->SetTitle("#sigma_{95%} #times BR(V' #rightarrow " + VV + ") [pb]"); // #rightarrow 2l2q
  grmedian_cls->SetMarkerStyle(24);//25=hollow squre
  grmedian_cls->SetMarkerColor(kBlack);
  grmedian_cls->SetLineStyle(2);
  grmedian_cls->SetLineWidth(3);

  if(obs){
    grobslim_cls->SetMarkerColor(kBlack);
    grobslim_cls->SetMarkerStyle(21);//24=hollow circle
    grobslim_cls->SetMarkerSize(1.0);
    grobslim_cls->SetLineStyle(1);
    grobslim_cls->SetLineWidth(3);

void drawGraphs(TGraphErrors* data, TGraphErrors* mc, const std::string& method, const std::string& xTitle, const std::string& yTitle, const std::string& legendTitle, double lumi, const std::string& outputName, int dataMarkerStyle = 20, int dataMarkerColor = kBlack, int mcMarkerStyle = 29, int mcMarkerColor = kBlue) {

  data->SetMarkerSize(1.5);
  data->SetMarkerStyle(dataMarkerStyle);
  data->SetMarkerColor(dataMarkerColor);
  data->SetLineColor(dataMarkerColor);

  mc->SetMarkerSize(1.5);
  mc->SetMarkerStyle(mcMarkerStyle);
  mc->SetMarkerColor(mcMarkerColor);
  mc->SetLineColor(mcMarkerColor);

  // Fit
  TF1* data_fct = nullptr;
  TF1* mc_fct = nullptr;

  TH1D* errors_data = new TH1D("errors_data", "errors", 100, 0, 1);
  errors_data->SetStats(false);
  errors_data->SetFillColor(LIGHT_GRAY);
  //errors_data->SetLineColor(kRed);
  errors_data->SetFillStyle(1001);

  TH1D* errors_mc = (TH1D*) errors_data->Clone("errors_mc");
  errors_mc->SetFillColor(LIGHT_BLUE);

  if (method == "Balancing") {
    data_fct = new TF1("data_fct", "[0] - x*x*[1]", 0, 1);
    data_fct->SetLineColor(dataMarkerColor);
    data_fct->SetLineWidth(1);
    data_fct->SetLineStyle(2);

    data->Fit(data_fct, "RQN");
    (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_data, 0.68);

    mc_fct = new TF1("mc_fct", "[0] - x*x*[1]", 0, 1);
    mc_fct->SetLineColor(mcMarkerColor);
    mc_fct->SetLineWidth(1);
    mc_fct->SetLineStyle(2);

    mc->Fit(mc_fct, "RQN");
    (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_mc, 0.68);
  } else {
    data_fct = new TF1("data_fct", "[0] + x*[1]", 0, 1);
    data_fct->SetLineColor(dataMarkerColor);
    data_fct->SetLineWidth(1);
    data_fct->SetLineStyle(2);

    data->Fit(data_fct, "RQN");
    (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_data, 0.68);

    mc_fct = new TF1("mc_fct", "[0] + x*[1]", 0, 1);
    mc_fct->SetLineColor(mcMarkerColor);
    mc_fct->SetLineWidth(1);
    mc_fct->SetLineStyle(2);

    mc->Fit(mc_fct, "RQN");
    (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_mc, 0.68);

  }

  data_fct->SetRange(0, 1);
  mc_fct->SetRange(0, 1);

  TMultiGraph* mg = new TMultiGraph();
  mg->Add(data);
  mg->Add(mc);

  TString title = TString::Format(";%s;%s", xTitle.c_str(), yTitle.c_str());
  mg->SetTitle(title);

  TCanvas* canvas = new TCanvas("canvas", "", 800, 800);

  mg->Draw("ap");

  errors_data->Draw("e3 same");
  errors_mc->Draw("e3 same");

  data_fct->Draw("same");
  mc_fct->Draw("same");

  mg->Draw("ap same");


  TLegend* legend = new TLegend(0.18, 0.18, 0.55, 0.35);
  legend->SetTextFont(42);
  legend->SetFillColor(kWhite);
  legend->SetFillStyle(0);
  legend->SetTextSize(0.035);
  legend->SetBorderSize(1);

  TString legendTitleWithPtCut = TString::Format("%s, p_{T}^{#gamma} #geq 170 GeV", legendTitle.c_str());

  legend->SetHeader(legendTitleWithPtCut);
  legend->AddEntry(data, TString::Format("%s (data)", method.c_str()), "p");
  legend->AddEntry(mc, TString::Format("%s (MC)", method.c_str()), "p");
  legend->Draw();

  TLatex tl;
  tl.SetNDC();
  tl.SetTextSize(0.035);
//.........这里部分代码省略.........

void drawCombinedGraphs(TGraphErrors* balancingData, TGraphErrors* balancingMC, TGraphErrors* mpfData, TGraphErrors* mpfMC, const std::string& xTitle, const std::string& yTitle, const std::string& legendTitle, double lumi, const std::string& outputName) {

  balancingData->SetMarkerSize(1.5);
  balancingData->SetMarkerStyle(20);
  balancingData->SetMarkerColor(kBlack);
  balancingData->SetLineColor(kBlack);

  mpfData->SetMarkerSize(1.5);
  mpfData->SetMarkerStyle(20);
  mpfData->SetMarkerColor(kRed);
  mpfData->SetLineColor(kRed);

  balancingMC->SetMarkerSize(1.5);
  balancingMC->SetMarkerStyle(29);
  balancingMC->SetMarkerColor(kBlue);
  balancingMC->SetLineColor(kBlue);

  mpfMC->SetMarkerSize(1.5);
  mpfMC->SetMarkerStyle(29);
  mpfMC->SetMarkerColor(46);
  mpfMC->SetLineColor(46);
  
  TH1D* errors_bal_data = new TH1D("errors_bal_data", "errors", 100, 0, 1);
  errors_bal_data->SetStats(false);
  errors_bal_data->SetFillColor(LIGHT_GRAY);
  errors_bal_data->SetFillStyle(1001);

  TH1D* errors_bal_mc = (TH1D*) errors_bal_data->Clone("errors_bal_mc");
  errors_bal_mc->SetFillColor(LIGHT_BLUE);

  TH1D* errors_mpf_data = new TH1D("errors_mpf_data", "errors", 100, 0, 1);
  errors_mpf_data->SetStats(false);
  errors_mpf_data->SetFillColor(LIGHT_RED);
  errors_mpf_data->SetFillStyle(1001);

  TH1D* errors_mpf_mc = (TH1D*) errors_bal_data->Clone("errors_mpf_mc");
  errors_mpf_mc->SetFillColor(LIGHT_MARRON);

  TF1* balancingData_fct = new TF1("balancingData_fct", "[0] - x*x*[1]", 0, 1);
  balancingData_fct->SetLineColor(kBlack);
  balancingData_fct->SetLineWidth(1);
  balancingData_fct->SetLineStyle(2);

  balancingData->Fit(balancingData_fct, "QRN");
  (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_bal_data, 0.68);

  TF1* balancingMC_fct = new TF1("mc_fct", "[0] - x*x*[1]", 0, 1);
  balancingMC_fct->SetLineColor(kBlue);
  balancingMC_fct->SetLineWidth(1);
  balancingMC_fct->SetLineStyle(2);

  balancingMC->Fit(balancingMC_fct, "QRN");
  (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_bal_mc, 0.68);

  TF1* mpfData_fct = new TF1("mpfData_fct", "[0] + x*[1]", 0, 1);
  mpfData_fct->SetLineColor(kRed);
  mpfData_fct->SetLineWidth(1);
  mpfData_fct->SetLineStyle(2);

  mpfData->Fit(mpfData_fct, "QRN");
  (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_mpf_data, 0.68);

  TF1* mpfMC_fct = new TF1("mc_fct", "[0] + x*[1]", 0, 1);
  mpfMC_fct->SetLineColor(46);
  mpfMC_fct->SetLineWidth(1);
  mpfMC_fct->SetLineStyle(2);

  mpfMC->Fit(mpfMC_fct, "QRN");
  (TVirtualFitter::GetFitter())->GetConfidenceIntervals(errors_mpf_mc, 0.68);

  TString balancing_ratio_legend = TString::Format("#color[4]{#scale[1]{r_{bal} = %.03f #pm %.03f}}", balancingData_fct->GetParameter(0) / balancingMC_fct->GetParameter(0), sqrt(pow(balancingData_fct->GetParError(0), 2) + pow(balancingMC_fct->GetParError(0), 2)));

  TString mpf_ratio_legend = TString::Format("#color[2]{#scale[1]{r_{MPF} = %.03f #pm %.03f}}", mpfData_fct->GetParameter(0) / mpfMC_fct->GetParameter(0), sqrt(pow(mpfData_fct->GetParError(0), 2) + pow(mpfMC_fct->GetParError(0), 2)));

  TMultiGraph* mg = new TMultiGraph();
  mg->Add(balancingData);
  mg->Add(balancingMC);
  mg->Add(mpfData);
  mg->Add(mpfMC);

  TString title = TString::Format(";%s;%s", xTitle.c_str(), yTitle.c_str());
  mg->SetTitle(title);

  TCanvas* canvas = new TCanvas("canvas", "", 800, 800);

  mg->Draw("ap");

  balancingData_fct->SetRange(0, 1);
  balancingMC_fct->SetRange(0, 1);
  mpfData_fct->SetRange(0, 1);
  mpfMC_fct->SetRange(0, 1);

  errors_bal_data->Draw("e3 same");
  errors_bal_mc->Draw("e3 same");

  errors_mpf_data->Draw("e3 same");
  errors_mpf_mc->Draw("e3 same");

  mg->Draw("ap same");

//.........这里部分代码省略.........

//.........这里部分代码省略.........
  matrixLeftInverse.Invert(det2);
  //   dumpElements(matrixLeftInverse);

  TMatrixD lambda(nbins,NELE);
  lambda= matrixLeftInverse*matrixRight;
  //   dumpElements(lambda);
  
  TMatrixD transposeLambda(NELE,nbins);
  transposeLambda.Transpose(lambda);
  //   dumpElements(transposeLambda);

  //   dumpElements(lambda);

  TVectorD combined_value(nbins);
  combined_value = lambda*measurement;
  //   dumpElements(combined_value);

  TMatrixD combined_error(nbins,nbins);
  combined_error = lambda*(errorM*transposeLambda);
  

  // after all matrix operation, now set the histogram
  for(int i=0; i<nbins;i++){

    h_combine->SetBinContent(i+1, combined_value(i));
  
    double error = combined_error(i,i);
  
    if(error>1e-10)error = sqrt(error);
    else error=1e-10;

    h_combine->SetBinError(i+1,error);

  }
  
  for(int i=1; i <= h_combine->GetNbinsX(); i++){

    if(h_combine->GetBinContent(i)>1e-10)
      cout << "Bin " << i << ": " << h_combine->GetBinContent(i) 
	   << " +- " << h_combine->GetBinError(i) << "\t"
	   << "rele = " << 
	(h_combine->GetBinError(i)/h_combine->GetBinContent(i)) 
	   << endl;

  }


  h_e      ->SetXTitle(xtitle.data());
  h_mu     ->SetXTitle(xtitle.data());
  h_combine->SetXTitle(xtitle.data());

  std::string command = "recreate";
  if(update)command="update";
  TFile* outFile = new TFile("bigmatrix_corr.root", command.data());
  h_e      ->Write();
  h_mu     ->Write();
  h_combine->Write();
  heff_e   ->Write();
  heff_mu  ->Write();
  outFile->Close();


  TCanvas* c1 = new TCanvas("c1","",500,500);
  if(logScale)
    c1->SetLogy(1);
  float x1NDC = 0.67;
  float y1NDC = 0.764831;
  float x2NDC = 0.830;
  float y2NDC = 0.908898;

  h_e->Draw("e");
  h_mu->Draw("esame");
  h_combine->Draw("e1same");

  TLegend* leg = new TLegend(x1NDC,y1NDC,x2NDC,y2NDC);  
  
  leg->SetFillColor(0);
  leg->SetFillStyle(0);
  leg->SetTextSize(0.05);
  leg->SetBorderSize(0);
  leg->AddEntry(h_combine, "combined");
  leg->AddEntry(h_e, "e");
  leg->AddEntry(h_mu, "#mu");
  leg->Draw("same");

  std::string remword  ="h_";
  std::string remword2 ="h";
  size_t pos  = eikoName.find(remword);
  if(pos!= std::string::npos)
    eikoName.replace(pos,remword2.length(),"");

  string dirName = "fig";
  gSystem->mkdir(dirName.data());

  string fileName = dirName + "/" + "final_emu" + eikoName;
  c1->Print(Form("%s.eps",fileName.data()));
  c1->Print(Form("%s.gif",fileName.data()));
  c1->Print(Form("%s.pdf",fileName.data()));

}

int main() {
  ModTDRStyle();

  std::vector<Scan> scans;
  //scans.push_back({"higgsCombinefullScan.MultiDimFit.mH125.root", "with syst.", 1, nullptr});
  //scans.push_back({"higgsCombinefastScan.MultiDimFit.mH125.root", "no syst.", 32, nullptr});
  //scans.push_back({"higgsCombinenoBBBScan.MultiDimFit.mH125.root", "no bbb syst.", 38, nullptr});
  
  // scans.push_back({"thesis/higgsCombineFullScan.MultiDimFit.mH125.root", "Stat+Syst+Theory", 1, nullptr});
  scans.push_back({"higgsCombineFullScan.MultiDimFit.mH125.root", "Mass Scan", kAzure-4, nullptr});
  // scans.push_back({"higgsCombineStatOnly.MultiDimFit.mH125.root", "Stat Only", kBlue+1, nullptr});
  //scans.push_back({"thesis/higgsCombineStatAndTh.MultiDimFit.mH125.root", "Stat+Theory", 39, nullptr});
  TCanvas c1("canvas","canvas");

  std::vector<TLine *> lines;


  TLegend *leg = new TLegend(0.37,0.65,0.73,0.9,"","brNDC");

  unsigned counter = 0;
  for (auto & sc : scans) {
    TFile f1(sc.file.c_str());
    TTree *t1 = dynamic_cast<TTree*>(f1.Get("limit"));
    double best1 = 0.0;
    sc.gr = new TGraph(ExtractGraph(t1, best1));
    auto x1 = GetCrossings(*(sc.gr), 1.0);
    TString res;
    if (x1.size() == 2) {
      double err = (x1[1]-x1[0])/2.0;
      std::cout << "Best fit is: " << best1 << " +/- " << err << std::endl;
      lines.push_back(new TLine(x1[0],0,x1[0],1.0));
      lines.back()->SetLineColor(sc.color);
      lines.back()->SetLineWidth(2);
      lines.push_back(new TLine(x1[1],0,x1[1],1.0));
      lines.back()->SetLineColor(sc.color);
      lines.back()->SetLineWidth(2);
      res = TString::Format("%.1f#pm%.1f",best1,err);
    }
    sc.gr->SetLineColor(sc.color);
    sc.gr->SetLineWidth(3);
    sc.gr->Draw(counter ? "LSAME" : "AL");
    TString leg_text = "#splitline{"+sc.label+"}{"+res+"}";
    leg->AddEntry(sc.gr, leg_text, "L");
    counter++;
  }
  // c1.cd();
  // // g1.Print();
  // g1.SetLineColor(1);
  // g1.SetLineWidth(2);
  // // g1.SetMarkerColor(7);
  // g1.Draw("AC");
  scans[0].gr->SetMaximum(4);
  scans[0].gr->GetXaxis()->SetRangeUser(110, 130);
  // scans[0].gr->GetXaxis()->SetTitle("Signal Strength, #mu");
  scans[0].gr->GetXaxis()->SetTitle("m_{H} [GeV]");
  scans[0].gr->GetYaxis()->SetTitle("-2 #Delta ln L");
  scans[0].gr->GetXaxis()->SetTitleFont(62);
  scans[0].gr->GetYaxis()->SetTitleFont(62);
  leg->SetBorderSize(1);
  leg->SetTextFont(42);
  leg->SetTextSize(0.03);
  leg->SetLineColor(0);
  leg->SetLineStyle(1);
  leg->SetLineWidth(1);
  leg->SetFillColor(0);
  leg->SetFillStyle(1001);
  leg->Draw();
  lines.push_back(new TLine(110,1,130,1));
  lines.back()->SetLineColor(2);
  for (auto l : lines) l->Draw();
  TLatex *title_latex = new TLatex();
  title_latex->SetNDC();
  title_latex->SetTextSize(0.035);
  title_latex->SetTextFont(62);
  title_latex->SetTextAlign(31);
  double height = 0.94;
  title_latex->DrawLatex(0.95,height,"19.7 fb^{-1} (8 TeV) + 4.9 fb^{-1} (7 TeV)");
  title_latex->SetTextAlign(11);
  title_latex->DrawLatex(0.17,height,"H#rightarrow#tau#tau");
  title_latex->SetTextSize(0.08);
  title_latex->DrawLatex(0.21, 0.25, "#mu#tau_{h}");
  c1.Update();
  c1.SaveAs("scan.pdf");
  return 0;
}