C++ setTDRStyle函数代码示例

void efficiencyL0_01(std::string fullPath, bool doOverview = true)
{
    const std::string outpdf("efficiencyL0_01plot");
    const int fVerbose(0);
    setTDRStyle();
    gStyle->SetOptStat(112211);
    gStyle->SetPalette(1);
    // Canvas
    if (doOverview)
	c = new TCanvas("c1","c1",1400,1200);
    else
	c = new TCanvas("c1","c1",600,600);
    const unsigned int nPadX = doOverview ? 3 : 1;
    const unsigned int nPadY = doOverview ? 2 : 1;
    c->Divide(nPadX,nPadY);
    const unsigned int nPads=nPadX*nPadY;
    // Open file
    TFile *f = TFile::Open(fullPath.c_str());
    if (f==0)
    {
	cout << "File " << fullPath << " not found -- exiting" << endl;
	return;
    }
    if(fVerbose>0)
	cout << "Succesfully opened file " << fullPath << endl;
    // Get TTrees
    TTree* tree = (TTree*) f->Get("events");
    if (tree==0)
    {
	cout << "Tree events not found -- exiting" << endl;
	return;
    }
    if(fVerbose>0) cout << "Got TTree evemts with " << tree->GetEntries() << " entries" << endl;

    // General cuts
    std::string cutacc = "genL0vtxR>1&&genL0vtxR<35&&TMath::Abs(genL0vtxZ)<100";

    // Do plots
    int canvasCdCounter(0), canvasPageCounter(0);
    const int nbinseta(5);
    const int nbinspt(5);
    const double etamax(2.6);
    const double ptmax(40.);
    //const double ptbins[] = {0,1,3,5,10,15,40};
    const double ptbins[] = {0,0.5,1,1.5,2,2.5,3,3.5,4,5,6,8,10,40};
    const int ptbins_size = (sizeof(ptbins)/sizeof(double));
    std::vector<double> ptbinvec(ptbins,ptbins+ptbins_size);
    //const double etabins[] = {0.0,0.5,1.0,1.4,1.8,2.2};
    const double etabins[] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0};
    const int etabins_size = (sizeof(etabins)/sizeof(double));
    std::vector<double> etabinvec(etabins,etabins+etabins_size);
    //const int nbinseta(2);
    //const double etamax(2.4);
    //const int nbinspt(1);
    const bool etabetrag(true);

    canvaspager(c,outpdf,nPads,canvasCdCounter,canvasPageCounter);
    doPlot2d((TPad*)c->cd(canvasCdCounter),tree,"h1", "genL0pt:TMath::Abs(genL0eta)", cutacc.c_str(), etabinvec, ptbinvec,"Generated #Lambda","|#eta(#Lambda)|","p_{T}(#Lambda)","","GeV/c");
    setTH2params((TPad*)c->cd(canvasCdCounter),(TH2F*)gDirectory->GetList()->FindObject("h1"),true);

    canvaspager(c,outpdf,nPads,canvasCdCounter,canvasPageCounter);
    doPlot2d((TPad*)c->cd(canvasCdCounter),tree,"h2", "genL0pt:TMath::Abs(genL0eta)", (cutacc+(cutacc.size()>0?"&&":"")+"L0matched==1").c_str(), etabinvec, ptbinvec,"Generated #Lambda, found in reco","|#eta(#Lambda)|","p_{T}(#Lambda)","","GeV/c");
    setTH2params((TPad*)c->cd(canvasCdCounter),(TH2F*)gDirectory->GetList()->FindObject("h2"),true);

#ifdef DOPIDTABLE
    // now we make a pidTable out of this
    PidTable *pid = new PidTable(1);
    pid->readFromHist(gDirectory, "h2", "h1");
    pid->dumpToFile("pid_lambda0.dat");
#endif

    // make a clone for the ratio
    {
	TH2F *hratio = (TH2F*)gDirectory->GetList()->FindObject("h2")->Clone("h2Dratio");
    }
    canvaspager(c,outpdf,nPads,canvasCdCounter,canvasPageCounter);
    doPlotRatio2d((TPad*)c->cd(canvasCdCounter),"#Lambda reco efficiency","h1","h2Dratio");

    // Now we do a 1D histo for each x-bin
    //plot1Dfrom2DforeachXbin((TH2F*)gDirectory->GetList()->FindObject("h1"));
    canvaspager(c,outpdf,nPads,canvasCdCounter,canvasPageCounter);
    plot1Dfrom2DforeachXbin((TPad*)c->cd(canvasCdCounter),"h2Dratio");

    // 1D ratio plots
    const double ptbinsratio[] = {0,0.5,1,1.5,2,2.5,3,3.5,4,5,6,8,12,16,40};
    const int ptbinsratio_size = (sizeof(ptbinsratio)/sizeof(double));
    std::vector<double> ptbinratiovec(ptbinsratio,ptbinsratio+ptbinsratio_size);
    const double etabinsratio[] = {0,0.2};
    const int etabinsratio_size = (sizeof(etabinsratio)/sizeof(double));
    std::vector<double> etabinratiovec(etabinsratio,etabinsratio+etabinsratio_size);

    std::string cutaccEta = (cutacc.size()>0?cutacc+"&&":"")+ "TMath::Abs(genL0eta)<2.5";
    std::string cutaccPt = (cutacc.size()>0?cutacc+"&&":"")+ "genL0pt>2&&genL0pt<40";

    canvaspager(c,outpdf,nPads,canvasCdCounter,canvasPageCounter);
    doRatioPlot((TPad*)c->cd(canvasCdCounter),tree,tree,"hpt","genL0pt","genL0pt",cutaccEta.c_str(), (cutaccEta+(cutaccEta.size()>0?"&&":"")+"L0matched==1").c_str(),ptbinratiovec,"#Lambda reco efficiency for |#eta|<2.5","p_{T}","GeV/c");
    canvaspager(c,outpdf,nPads,canvasCdCounter,canvasPageCounter);
    doRatioPlot((TPad*)c->cd(canvasCdCounter),tree,tree,"heta","TMath::Abs(genL0eta)","TMath::Abs(genL0eta)",cutaccPt.c_str(), (cutaccPt+(cutaccPt.size()>0?"&&":"")+"L0matched==1").c_str(),25,0,2.5,"#Lambda reco efficiency for 2<p_{T}<40","|#eta|","");

    // finalize current page
//.........这里部分代码省略.........

void plotResolution( char* var, float xmin, float xmax, const char * region, const char *xbinning, const char * skim, int regionCode, 
		     std::ofstream * resultsStg1,
		     std::ofstream * resultsStg2 ) {
  
  //Log file
  
  TString log_subdir("./logs/");
  TString resolution_log = TString("resolution_PVbins_") + TString(var) + TString("_") + TString(skim) + TString(".log");

  std::ofstream * logfile = new std::ofstream( (log_subdir + resolution_log).Data(), ios_base::app );

  TDatime *d1 = new TDatime();
  (*logfile) << d1->AsString() << std::endl;

  //Output path
  TString path("./slhc-plots/resolution/Taus/noPUC/isoStudies");
  
  TString varName = TString("reco") + TString(var);  // RECO(PFTau)
  TString l1varName = TString("l1g") + TString(var); // L1 Calo Upgrade (can be Stage 1 or Stage 2)
  
  gROOT->SetStyle("Plain");
  gROOT->SetBatch(false);
  gStyle->SetOptStat(0);

  // --- Use the CMS TDR style
  gROOT->ProcessLine(".L tdrStyle.C");
  setTDRStyle();
  tdrStyle->SetErrorX(0.5);
  tdrStyle->SetPadLeftMargin(0.18);
  tdrStyle->SetPadRightMargin(0.08);
  tdrStyle->SetLegendBorderSize(0);
  tdrStyle->SetTitleYOffset(1.3);
  tdrStyle->SetOptStat(0);
  tdrStyle->SetOptFit(0);
  tdrStyle->SetTitleFontSize(0.05);
  tdrStyle->SetStatStyle(0);

  TFile * f1 = new TFile("../L1Tree_Latest_MuTau_2012CD.root" );

  f1->cd();

  TTree * Stage1Taus = (TTree*)gDirectory->Get("rlxTauSt1Efficiency/Ntuple");

  TTree * Stage2Taus = (TTree*)gDirectory->Get("rlxTauSt2Efficiency/Ntuple");
  
  std::cout << " Ntuple ready: " << Stage1Taus << " " << Stage2Taus << std::endl;
  
  TCanvas * canvas = new TCanvas("asdf", "adsf", 800, 600);
  canvas->Draw();
  
  TList * Stage1Histos = new TList();
  TList * Stage2Histos = new TList();
  
  //this is to separate into different decay mode as in reco::PFTau::hadronicDecayMode
  // Enumerator:
  // -1 kNull
  //  0 kOneProng0PiZero
  //  1 kOneProng1PiZero
  //  2 kOneProng2PiZero
  //  ...
  // 10 kThreeProng0PiZero

  int decayMode[10] = {0};
  
  for( int k = 0 ; k < 1; ++k ) {
    
    double m1   = 0.0;
    double area = 0.0;
    double yMax = 0.0;
    
    int xMin = xmin;
    int xMax = xmax;
    
    char cuts[200];
    
    //this is to use the decayMode:
    sprintf(cuts, " ( l1Pt >= 25.0 && l1gMatch >= 1.0 )&& %s && nPVs > %d && nPVs < %d && decayMode >= %d", region, xMin, xMax, decayMode[k]);
    
    TString histoDraw = TString("(") + TString( l1varName ) + TString(" - ") + TString(varName) + TString(")/") + TString(varName) + TString(">> ");
    
    char st1HistoName[100];
    sprintf(st1HistoName, "htempSt1_%d", k);
    
    char st2HistoName[100];
    sprintf(st2HistoName, "htempSt2_%d", k);

    TH1D* h1T = new TH1D(st1HistoName, "", 50, -2, 2);
    h1T->Sumw2();

    TH1D* h2T = new TH1D(st2HistoName, "", 50, -2, 2);
    h2T->Sumw2();

    Stage1Taus->Draw( histoDraw + TString(st1HistoName), cuts, "");
    TH1F* h1 = (TH1F*)gDirectory->Get( st1HistoName )->Clone();
    area = h1->Integral();
    h1->Scale( 1.0/area );
    
    m1 = h1->GetMaximum();
    
    if ( m1 > yMax ) 
//.........这里部分代码省略.........

void plot(char* h, bool norm2data, TString prefix, TString RR) {

    TFile f1(prefix+"DATA"+RR+"CP.root", "read");
    TH1D* hdata = (TH1D*)f1.Get(h);
    hdata->SetDirectory(0);
    
    TFile f2(prefix+"DYM_CP.root", "read");
    TH1D* hsignal = (TH1D*)f2.Get(TString(h));
    hsignal->SetDirectory(0);
    
    TFile f3(prefix+"QCD_CP_dd.root", "read");
    TH1D* hqcd = (TH1D*)f3.Get("hqcd");
    hqcd->SetDirectory(0);

    TFile f3_mc(prefix+"QCD_CP.root", "read");
    TH1D* hqcd_mc = (TH1D*)f3_mc.Get(h);
    hqcd_mc->SetDirectory(0);

    //FIXME    
    cout << "See the implementation in 2D case!" << endl;
    TFile f4(prefix+"TT_CP.root", "read");
    TH1D* httbar = (TH1D*)f4.Get(h); //"httbar");
    httbar->SetDirectory(0);

    TFile f4_mc(prefix+"TT_CP.root", "read");
    TH1D* httbar_mc = (TH1D*)f4_mc.Get(h);
    httbar_mc->SetDirectory(0);

    TFile f5(prefix+"TAU_CP.root", "read");
    TH1D* hZtautau = (TH1D*)f5.Get(h);
    hZtautau->SetDirectory(0);
    
    TFile f6(prefix+"EWK_CP.root", "read");
    TH1D* hWleptonnu = (TH1D*)f6.Get(h);
    hWleptonnu->SetDirectory(0);
    
    TFile f7(prefix+"DIBOSON_CP.root", "read");
    TH1D* hDibosons = (TH1D*)f7.Get(h);
    hDibosons->SetDirectory(0);
    
    gStyle->SetOptStat(0);
    gStyle->SetPalette(1);
    gROOT->ProcessLine(".L ../tools/setTDRStyle.C");
    setTDRStyle();
    gROOT->SetStyle("tdrStyle");
    gROOT->ForceStyle(true);
    
    // signal
    hsignal->SetLineColor(kOrange);
    hsignal->SetFillColor(kOrange);
    
    // qcd
    hqcd->SetLineColor(9);
    hqcd->SetFillColor(9);

    // tt
    httbar->SetLineColor(kRed+2);
    httbar->SetFillColor(kRed+2);

    // tau
    hZtautau->SetLineColor(kGreen);
    hZtautau->SetFillColor(kGreen);

    // EWK
    hWleptonnu->SetLineColor(38);
    hWleptonnu->SetFillColor(38);

    //hupsilon->SetLineColor(51);
    //hupsilon->SetFillColor(51);

    // diboson
    hDibosons->SetLineColor(40);
    hDibosons->SetFillColor(40);

    //
    //The histograms come properly weighted with priors
    //This would normalize to number of events in data
    //

    if (norm2data) {
	    const double ttbarNorm = httbar->Integral();
	    const double WleptonnuNorm = hWleptonnu->Integral();
	    const double ZtautauNorm = hZtautau->Integral();
	    const double qcdNorm = hqcd_mc->Integral();
 	    const double dibosonNorm = hDibosons->Integral();
	    const double signalNorm = hsignal->Integral();
 	    const double fullMCintegral = ttbarNorm+WleptonnuNorm+ZtautauNorm+qcdNorm+dibosonNorm+signalNorm;

	    httbar->Scale(hdata->Integral()/fullMCintegral);
	    hWleptonnu->Scale(hdata->Integral()/fullMCintegral);
	    hZtautau->Scale(hdata->Integral()/fullMCintegral);
	    hqcd_mc->Scale(hdata->Integral()/fullMCintegral);
	    hsignal->Scale(hdata->Integral()/fullMCintegral);
	    hDibosons->Scale(hdata->Integral()/fullMCintegral);
            normalizeToPeak(httbar, hWleptonnu, hqcd_mc, hZtautau, hsignal, hDibosons, hdata);

    } else {
           //Just do not normalize, for utility purposes, not for the analysis 
    }

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

void plotLimit(string outputDir="./", TString inputs="", TString inputs_blinded="", TString inputXSec="", bool strengthLimit=true, bool blind=false, double energy=7, double luminosity=5.035, TString legendName="ee and #mu#mu channels")
{
   setTDRStyle();  
   gStyle->SetPadTopMargin   (0.05);
   gStyle->SetPadBottomMargin(0.12);
   gStyle->SetPadRightMargin (0.16);
   gStyle->SetPadLeftMargin  (0.14);
   gStyle->SetTitleSize(0.04, "XYZ");
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.45);
   gStyle->SetPalette(1);
   gStyle->SetNdivisions(505);
  
  //get the limits from the tree
  TFile* file = TFile::Open(inputs);
  printf("Looping on %s\n",inputs.Data());
  if(!file) return;
  if(file->IsZombie()) return;
  TFile* file_blinded = TFile::Open(inputs_blinded);
  printf("Looping on %s\n",inputs_blinded.Data());
  if(!file_blinded) return;
  if(file_blinded->IsZombie()) return;
  TTree* tree_blinded = (TTree*)file_blinded->Get("limit");
  tree_blinded->GetBranch("mh"              )->SetAddress(&Tmh      );
  tree_blinded->GetBranch("limit"           )->SetAddress(&Tlimit   );
  tree_blinded->GetBranch("limitErr"        )->SetAddress(&TlimitErr);
  tree_blinded->GetBranch("quantileExpected")->SetAddress(&TquantExp);
  TGraph* ExpLimitm2 = getLimitGraph(tree_blinded,0.025);
  TGraph* ExpLimitm1 = getLimitGraph(tree_blinded,0.160);
  TGraph* ExpLimit   = getLimitGraph(tree_blinded,0.500);
  TGraph* ExpLimitp1 = getLimitGraph(tree_blinded,0.840);
  TGraph* ExpLimitp2 = getLimitGraph(tree_blinded,0.975);
  file_blinded->Close(); 
  TTree* tree = (TTree*)file->Get("limit");
  tree->GetBranch("mh"              )->SetAddress(&Tmh      );
  tree->GetBranch("limit"           )->SetAddress(&Tlimit   );
  tree->GetBranch("limitErr"        )->SetAddress(&TlimitErr);
  tree->GetBranch("quantileExpected")->SetAddress(&TquantExp);
  TGraph* ObsLimit   = getLimitGraph(tree,-1   ); 
  file->Close(); 

  FILE* pFileSStrenght = fopen((outputDir+"SignalStrenght").c_str(),"w");
  std::cout << "Printing Signal Strenght" << std::endl;
  for(int i=0;i<ExpLimit->GetN();i++){
     double M = ExpLimit->GetX()[i];
     std::cout << "Mass: " << M << "; ExpLimit: " << ExpLimit->Eval(M) << std::endl; 
     printf("$%8.6E$ & $%8.6E$ & $[%8.6E,%8.6E]$ & $[%8.6E,%8.6E]$ \\\\\\hline\n",M, ExpLimit->Eval(M), ExpLimitm1->Eval(M), ExpLimitp1->Eval(M), ExpLimitm2->Eval(M),  ExpLimitp2->Eval(M));
     fprintf(pFileSStrenght, "$%8.6E$ & $%8.6E$ & $[%8.6E,%8.6E]$ & $[%8.6E,%8.6E]$ & $%8.6E$ \\\\\\hline\n",M, ExpLimit->Eval(M), ExpLimitm1->Eval(M), ExpLimitp1->Eval(M), ExpLimitm2->Eval(M),  ExpLimitp2->Eval(M), ObsLimit->Eval(M));
    if(int(ExpLimit->GetX()[i])%50!=0)continue; //printf("%f ",ObsLimit->Eval(M));
  }printf("\n");
  fclose(pFileSStrenght); 
 

  //get the pValue
  inputs = inputs.ReplaceAll("/LimitTree", "/PValueTree");
  file = TFile::Open(inputs);
  
  printf("Looping on %s\n",inputs.Data());
  if(!file) return;
  if(file->IsZombie()) return;
  
  tree = (TTree*)file->Get("limit");
  
  tree->GetBranch("limit"           )->SetAddress(&Tlimit   );
  
  TGraph* pValue     = getLimitGraph(tree,-1);
  
  file->Close();

  
  //make TH Cross-sections
   string suffix = outputDir;
   TGraph* THXSec   = Hxswg::utils::getXSec(outputDir); 
   scaleGraph(THXSec, 1000);  //convert cross-section to fb
   double cprime=1.0; double  brnew=0.0;
   double XSecScaleFactor = 1.0;
   if(suffix.find("_cp")!=string::npos){
     sscanf(suffix.c_str()+suffix.find("_cp"), "_cp%lf_brn%lf", &cprime, &brnew);
     XSecScaleFactor = pow(cprime,2) * (1-brnew);
   }
  //XSecScaleFactor = 0.001; //pb to fb
  scaleGraph(THXSec, XSecScaleFactor);


  string prod = "pp_SM";
  if(outputDir.find("ggH")!=std::string::npos)prod="gg";
  if(outputDir.find("qqH")!=std::string::npos)prod="qq";
  if(outputDir.find("ppH")!=std::string::npos)prod="pp";

  
  strengthLimit = false;
  if(prod=="pp_SM")strengthLimit=true;
 
  //TGraph *XSecMELA = Hxswg::utils::getXSecMELA(cprime);

  //Hxswg::utils::multiplyGraph(   ObsLimit, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitm2, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitm1, XSecMELA);
  //Hxswg::utils::multiplyGraph(   ExpLimit, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitp1, XSecMELA);
//.........这里部分代码省略.........

void DrawResponsePlot( Float_t genMean[], Float_t genSigma[], 
		       Float_t recoMean[], Float_t recoSigma[]) {


  gROOT->ProcessLine(".L mystyle.C");
  setTDRStyle();
  tdrStyle->SetErrorX(0.5);
  tdrStyle->SetPadLeftMargin(0.2);
  tdrStyle->SetPadRightMargin(0.10);
  tdrStyle->SetLegendBorderSize(0);
  tdrStyle->SetTitleYOffset(1.3);

  // plot full spectrum
  TGraphErrors *ptbalanceGen  = new TGraphErrors(nZPtBins, pt, genMean, 
						 errpt, genSigma);
  TGraphErrors *ptbalanceReco = new TGraphErrors(nZPtBins, pt, recoMean, 
						 errpt, recoSigma);
  // plot Zmumu values
  Float_t ptmm[9] = { 40.0, 60.0, 100.0, 140.0, 200.0, 250.0, 
		      330.0, 400.0, 520.0 };
  Float_t balancemm[9] = { 0.496, 0.568, 0.66, 0.71, 0.75, 0.765, 
			   0.775, 0.79, 0.81 }; 
  TGraph *ptbalancemm = new TGraph( 9, ptmm, balancemm);


  ptbalanceGen->GetXaxis()->SetTitle("p_{T}^{Z} (GeV/c)");
  ptbalanceGen->GetYaxis()->SetTitle("p_{T}^{jet} / p_{T}^{Z}");
  ptbalanceGen->SetMarkerStyle(22);
  ptbalanceGen->SetMarkerSize(1.2);
  ptbalanceGen->SetTitle("");
  ptbalanceGen->SetMinimum(0.3);
  ptbalanceReco->SetMarkerColor(2);
  ptbalanceReco->SetLineColor(2);
  ptbalanceReco->SetMarkerStyle(22);
  ptbalanceReco->SetMarkerSize(1.2);
  ptbalanceReco->SetMinimum(0.3);
  ptbalancemm->SetMarkerStyle(24);
  ptbalancemm->SetMarkerSize(1.2);
  ptbalancemm->SetMinimum(0.3);
  ptbalancemm->SetMarkerColor(4);
  ptbalancemm->SetLineColor(4);

  TCanvas c1("c1","",500,500);
  ptbalanceGen->Draw("APL");
  ptbalanceReco->Draw("PL");
  ptbalancemm->Draw("PL");
  leg_hist = new TLegend(0.6,0.7,0.89,0.89);
  leg_hist->AddEntry( ptbalanceGen, "Generator level", "l");
  leg_hist->AddEntry( ptbalanceReco,"Calorimeter level","l");
  leg_hist->AddEntry( ptbalancemm,"Z#rightarrow#mu#mu","l");
  leg_hist->SetFillColor(0);
  leg_hist->Draw();
  c1.SaveAs("PtBalanceVsPt.eps");
  c1.SaveAs("PtBalanceVsPt.gif");
  c1.SaveAs("PtBalanceVsPt.root");
  c1.Close();
  delete leg_hist;
  delete ptbalanceGen;
  delete ptbalanceReco;
  delete ptbalancemm;
}

void PlotEffGen(TString fileName = "DYJetsEff",
		// TString numName = "TPPassed",
		// TString denName = "TPAll",
		//		TString numName = "RecoPassedTag",
		//		TString denName = "RecoAllTag",
		TString numName = "RecoPassed",
		TString denName = "RecoAll",
		bool posTag = true,
		bool isMuId = true,
		bool isPrompt = true) {


  setTDRStyle();

  TString name("ElectronPt_");
  if (isMuId)
    name = "MuonPt_";

  TString LepQ("Neg");
  if (posTag)
    LepQ = "Pos";

  TString Type("NonPrompt");
  if (isPrompt)
    Type = "Prompt";

  int nEtaBins = 3;
  TString EtaBinName[3] = {"0To0p8",
			   "0p8To1p5",
			   "1p5To2p3"};
  TFile * file = new TFile(fileName+".root");

  int nPtBins = 7;
  float ptBins[8] = {10,15,20,30,40,50,70,100};

  TH1F * histos[3][2];
  for (int iEta=0; iEta<3; ++iEta) {
    TH1F * histNum = (TH1F*)file->Get(name+LepQ+Type+numName+EtaBinName[iEta]);
    histos[iEta][0] = TH1toTH1(histNum,nPtBins,ptBins,true,"_new");
    TH1F * histDen = (TH1F*)file->Get(name+LepQ+Type+denName+EtaBinName[iEta]);
    histos[iEta][1] = TH1toTH1(histDen,nPtBins,ptBins,true,"_new");
    
  }

  int color[3] = {1,2,4};
  int symbol[3] = {20,21,22};

  TGraphAsymmErrors * eff[3];
  for (int iEta=0; iEta<3; ++iEta) {
    eff[iEta] = new TGraphAsymmErrors();
    eff[iEta]->SetLineColor(color[iEta]);
    eff[iEta]->SetMarkerColor(color[iEta]);
    eff[iEta]->SetLineWidth(2);
    eff[iEta]->SetMarkerSize(1.6);
    eff[iEta]->SetMarkerStyle(symbol[iEta]);
    eff[iEta]->Divide(histos[iEta][0],histos[iEta][1],"Pois");
  }  

 
  TH2F * frame = new TH2F("frame","",2,0,100,2,0,1);
  TCanvas * canv = new TCanvas("canv","",700,700);
  frame->Draw();
  eff[0]->Draw("PS");
  eff[1]->Draw("PS");
  eff[2]->Draw("PS");
  canv->Update();

}

bool                                          // returns true if success
processStyleSection(FILE *fp,string& theline, bool& new_section)
{
  vector<string> v_tokens;

  string *sid      = NULL;
  TStyle *thestyle = NULL;

  if (gl_verbose)
    cout << "Processing style section" << endl;

  new_section=false;

  while (getLine(fp,theline,"style")) {
    if (!theline.size()) continue;
    if (theline[0] == '#') continue; // comments are welcome

    if (theline[0] == '[') {
      new_section=true;
      return true;
    }

    string key, value;
    if (!getKeyValue(theline,key,value)) continue;

    if (key == "id") {
      if (sid != NULL) {
	cerr << "no more than one id per style section allowed " << value << endl;
	break;
      }

      sid = new string(value);
      thestyle = new TStyle(*gStyle); // Assume current attributes, let user override specifics
      thestyle->SetNameTitle(sid->c_str(),sid->c_str());
      glmap_id2style.insert(pair<string,TStyle*>(*sid,thestyle));
      continue;
    }
    else if (!sid) {
      // assume the style is the global style
      thestyle = gStyle;
      sid = new string("using global"); // so as to pass this check the next time
      cerr << "No style ID defined, assuming global style" << endl;
      continue;
    }

    if (key == "style") {
      if (value == "TDR")
	setTDRStyle();
      else if (value == "Plain")
	gROOT->SetStyle("Plain");
      else 
	cerr << "unknown style name " << value << endl;
    }
    else if (key == "optstat")  {
      if (gl_verbose) cout << "OptStat = " << thestyle->GetOptStat() << endl;
      thestyle->SetOptStat(value.c_str());
      if (gl_verbose) cout << "OptStat = " << thestyle->GetOptStat() << endl;
    }
    else if (key == "opttitle")  thestyle->SetOptTitle(str2int(value));

    else if (key == "padrightmargin")  thestyle->SetPadRightMargin (str2flt(value));
    else if (key == "padleftmargin")   thestyle->SetPadLeftMargin  (str2flt(value));
    else if (key == "padtopmargin")    thestyle->SetPadTopMargin   (str2flt(value));
    else if (key == "padbottommargin") thestyle->SetPadBottomMargin(str2flt(value));
    else if (key == "padgridx")        thestyle->SetPadGridX       (str2int(value));
    else if (key == "padgridy")        thestyle->SetPadGridY       (str2int(value));

    // Set the position/size of the title box

    else if (key == "titlexndc")    thestyle->SetTitleX(str2flt(value));
    else if (key == "titleyndc")    thestyle->SetTitleY(str2flt(value));
    else if (key == "titlewndc")    thestyle->SetTitleW(str2flt(value));
    else if (key == "titlehndc")    thestyle->SetTitleH(str2flt(value));
    else if (key == "titlefont")    thestyle->SetTitleFont(str2int(value));
    else if (key == "titlebordersize") thestyle->SetTitleBorderSize(str2int(value));

    else if (key == "markercolor")  thestyle->SetMarkerColor(str2int(value));
    else if (key == "markerstyle")  thestyle->SetMarkerStyle(str2int(value));
    else if (key == "markersize")   thestyle->SetMarkerSize(str2int(value));
    else if (key == "linecolor")    thestyle->SetLineColor(str2int(value));
    else if (key == "linestyle")    thestyle->SetLineStyle(str2int(value));
    else if (key == "linewidth")    thestyle->SetLineWidth(str2int(value));
    else if (key == "fillcolor")    thestyle->SetFillColor(str2int(value));
    else if (key == "fillstyle")    thestyle->SetFillStyle(str2int(value));

    // axes
    else if (key == "xtitlesize")   thestyle->SetTitleSize(str2flt(value),"X");
    else if (key == "ytitlesize")   thestyle->SetTitleSize(str2flt(value),"Y");
    else if (key == "ztitlesize")   thestyle->SetTitleSize(str2flt(value),"Z");
    else if (key == "xtitleoffset") thestyle->SetTitleOffset(str2flt(value),"X");
    else if (key == "ytitleoffset") thestyle->SetTitleOffset(str2flt(value),"Y");
    else if (key == "ztitleoffset") thestyle->SetTitleOffset(str2flt(value),"Z");
    else if (key == "xlabeloffset") thestyle->SetLabelOffset(str2flt(value),"X");
    else if (key == "ylabeloffset") thestyle->SetLabelOffset(str2flt(value),"Y");
    else if (key == "zlabeloffset") thestyle->SetLabelOffset(str2flt(value),"Z");
    else if (key == "xtitlefont")   thestyle->SetTitleFont(str2int(value),"X");
    else if (key == "ytitlefont")   thestyle->SetTitleFont(str2int(value),"Y");
    else if (key == "ztitlefont")   thestyle->SetTitleFont(str2int(value),"Z");
    else if (key == "xlabelsize")   thestyle->SetLabelSize(str2flt(value),"X");
    else if (key == "ylabelsize")   thestyle->SetLabelSize(str2flt(value),"Y");
//.........这里部分代码省略.........

void plotEfficiency(char* var, double min, double max, int l1PtMin=0) {
  TFile*f1 =  TFile::Open("efficiency_tree_Zmumu_merged.root");
  f1->cd("muonEfficiency");
  TTree* tree1 = (TTree*) gDirectory->Get("Ntuple");

  int nBins = 100;
  if(var=="Pt") nBins = (int) ((max - min) / 5);
  if(var=="Eta") nBins = (int) ((max - min) / 0.2);



  ///// ----- Define your cuts ----------
  char l1PtMinStr[100];
  sprintf(l1PtMinStr, "l1Pt>%d", l1PtMin);
  TCut l1PtMinCut(l1PtMinStr);
  TCut cutNPV0("(nPVs>-1)");
  TCut cutNPV1("(nPVs>=0 && nPVs<15)");
  TCut cutNPV2("(nPVs>=15 && nPVs<25)");
  TCut cutNPV3("(nPVs>=25 && nPVs<200)");
  TCut l1isoCut("(abs(l1gDR)>0.2 || (l1gRegionEt/l1Pt<0.5))");


  ///// ----- Create all the histograms ----------
  TH1D* allReco = new TH1D("allReco", "", nBins, min, max);
  allReco->SetLineWidth(2);
  allReco->Sumw2();
  TH1D* allRecoNPV1 = allReco->Clone("allRecoNPV1");
  TH1D* allRecoNPV2 = allReco->Clone("allRecoNPV2");
  TH1D* allRecoNPV3 = allReco->Clone("allRecoNPV3");



  TH1D* l1iso = allReco->Clone("l1iso");
  TH1D* l1isoNPV1 = allReco->Clone("l1isoNPV1");
  TH1D* l1isoNPV2 = allReco->Clone("l1isoNPV2");
  TH1D* l1isoNPV3 = allReco->Clone("l1isoNPV3");


  l1isoNPV1->SetLineColor(4);
  l1isoNPV1->SetMarkerColor(4);
  l1isoNPV2->SetLineColor(6);
  l1isoNPV2->SetMarkerColor(6);
  l1isoNPV3->SetLineColor(2);
  l1isoNPV3->SetMarkerColor(2);

  l1iso->GetYaxis()->SetRangeUser(0.5, 1.1);
  l1iso->GetYaxis()->SetTitle("L1 Efficiency");
  char* xtitle = var;
  if(var=="Pt") xtitle = "Offline muon p_{T} (GeV)";
  if(var=="Eta") xtitle = "Offline muon #eta";
  l1iso->GetXaxis()->SetTitle(xtitle);





  ///// ----- Draw from tree ----------
  TString varName = TString("reco") + TString(var);
  tree1->Draw(varName+TString(">>allReco"),cutNPV0,"goff");
  tree1->Draw(varName+TString(">>allRecoNPV1"),cutNPV1,"goff");
  tree1->Draw(varName+TString(">>allRecoNPV2"),cutNPV2,"goff");
  tree1->Draw(varName+TString(">>allRecoNPV3"),cutNPV3,"goff");

  tree1->Draw(varName+TString(">>l1iso"),l1PtMinCut && cutNPV0 && l1isoCut,"goff");
  tree1->Draw(varName+TString(">>l1isoNPV1"),l1PtMinCut && cutNPV1 && l1isoCut,"goff");
  tree1->Draw(varName+TString(">>l1isoNPV2"),l1PtMinCut && cutNPV2 && l1isoCut,"goff");
  tree1->Draw(varName+TString(">>l1isoNPV3"),l1PtMinCut && cutNPV3 && l1isoCut,"goff");




  //// ---- Now divide by allReco histogram to obtain efficiency ----
  l1iso->Divide(allReco);
  l1isoNPV1->Divide(allRecoNPV1);
  l1isoNPV2->Divide(allRecoNPV2);
  l1isoNPV3->Divide(allRecoNPV3);


  // --- plot the efficiency histograms ------ 
  gROOT->ProcessLine(".L ~/tdrstyle.C");
  setTDRStyle();
  tdrStyle->SetErrorX(0.5);
  tdrStyle->SetPadLeftMargin(0.18);
  tdrStyle->SetPadRightMargin(0.08);
  tdrStyle->SetLegendBorderSize(0);
  tdrStyle->SetTitleYOffset(1.3);


  TCanvas* canEffIso = new TCanvas("canEffIso","",500,500);
  l1iso->SetFillColor(5);
  l1iso->Draw("e3");
  l1iso->Draw("esame");
  l1isoNPV1->Draw("psame");
  l1isoNPV2->Draw("psame");
  l1isoNPV3->Draw("psame");
  l1iso->Draw("esame");

  TLegend* legend = new TLegend(0.45,0.15,0.9,0.45);
  legend->SetFillColor(0);
  legend->AddEntry(l1iso, "All NPV", "PLE");  
//.........这里部分代码省略.........

void plotMC(){//main  
TFile *file[6];

file[0]= TFile::Open("../data/plot_T5Wg_mGl-800to1000.root");  
file[1]= TFile::Open("../data/plot_T5Wg_mGl-1050to1100.root");
file[2]= TFile::Open("../data/plot_T5Wg_mGl-1250to1300.root");
file[3]= TFile::Open("../data/plot_T5Wg_mGl-1350to1400.root");
file[4]= TFile::Open("../data/plot_T5Wg_mGl-1450to1500.root");
file[5]= TFile::Open("../data/plot_T5Wg_mGl-1550.root");

TH2F *p_SUSYMass[6];
TH1F *p_mcphotonET[6];
TH1F *p_mceleET[6];
TH1F *p_phoEBR9_true[6]; 
TH1F *p_phoEER9_true[6];
TH1F *p_phoEBHoverE_true[6];
TH1F *p_phoEEHoverE_true[6];
TH1F *p_phoEBsigma_true[6];
TH1F *p_phoEEsigma_true[6];
TH1F *p_phoEBChIso_true[6];
TH1F *p_phoEEChIso_true[6];
TH1F *p_phoEBNeuIso_true[6];
TH1F *p_phoEENeuIso_true[6];
TH1F *p_phoEBPhoIso_true[6];
TH1F *p_phoEEPhoIso_true[6];

TH1F *p_eleEBR9_true[6];
TH1F *p_eleEER9_true[6];

TH1F *p_eleEBR9_idselect[6];
TH1F *p_eleEER9_idselect[6];

TH1F *p_Mt[6];
TH1F *p_PhoELeDeltaR[6];
TH1F *p_PhoEleMass[6];
TH1F *p_selectPhoEleMass[6];

for(unsigned iF(0); iF < 6; iF++){
  p_SUSYMass[iF] = (TH2F*)file[iF]->Get("Mass");

  p_mcphotonET[iF] = (TH1F*)file[iF]->Get("photonET");
  p_mceleET[iF] = (TH1F*)file[iF]->Get("eleET");

  p_phoEBR9_true[iF] = (TH1F*)file[iF]->Get("p_phoEBR9");
  p_phoEBR9_true[iF]->GetXaxis()->SetRangeUser(0.5,1);
  p_phoEER9_true[iF] = (TH1F*)file[iF]->Get("p_phoEER9");
  p_phoEER9_true[iF]->GetXaxis()->SetRangeUser(0.5,1);
  p_phoEBHoverE_true[iF] = (TH1F*)file[iF]->Get("p_phoEBHoverE");
  p_phoEEHoverE_true[iF] = (TH1F*)file[iF]->Get("p_phoEEHoverE");
  p_phoEBsigma_true[iF] = (TH1F*)file[iF]->Get("p_phoEBsigma");
  p_phoEEsigma_true[iF] = (TH1F*)file[iF]->Get("p_phoEEsigma");
  p_phoEBChIso_true[iF] = (TH1F*)file[iF]->Get("p_phoEBChIso");
  p_phoEEChIso_true[iF] = (TH1F*)file[iF]->Get("p_phoEEChIso");
  p_phoEBNeuIso_true[iF] = (TH1F*)file[iF]->Get("p_phoEBNeuIso");
  p_phoEENeuIso_true[iF] = (TH1F*)file[iF]->Get("p_phoEENeuIso");
  p_phoEBPhoIso_true[iF] = (TH1F*)file[iF]->Get("p_phoEBPhoIso");
  p_phoEEPhoIso_true[iF] = (TH1F*)file[iF]->Get("p_phoEEPhoIso");

  p_eleEBR9_true[iF] = (TH1F*)file[iF]->Get("p_eleEBR9");
  p_eleEBR9_true[iF]->GetXaxis()->SetRangeUser(0.5,1);
  p_eleEER9_true[iF] = (TH1F*)file[iF]->Get("p_eleEER9");
  p_eleEER9_true[iF]->GetXaxis()->SetRangeUser(0.5,1);

  p_eleEBR9_idselect[iF] = (TH1F*)file[iF]->Get("p_eleEBR9_idselect");
  p_eleEBR9_idselect[iF]->GetXaxis()->SetRangeUser(0.5,1);
  p_eleEER9_idselect[iF] = (TH1F*)file[iF]->Get("p_eleEER9_idselect");
  p_eleEER9_idselect[iF]->GetXaxis()->SetRangeUser(0.5,1);

  p_Mt[iF] = (TH1F*)file[iF]->Get("Mt");
  p_PhoELeDeltaR[iF] = (TH1F*)file[iF]->Get("p_PhoELeDeltaR");
  p_PhoEleMass[iF] = (TH1F*)file[iF]->Get("p_PhoEleMass");
  p_selectPhoEleMass[iF] = (TH1F*)file[iF]->Get("p_selectPhoEleMass");
}

p_SUSYMass[0]->Add(p_SUSYMass[1]); 
p_SUSYMass[0]->Add(p_SUSYMass[2]);
p_SUSYMass[0]->Add(p_SUSYMass[3]);
p_SUSYMass[0]->Add(p_SUSYMass[4]);
p_SUSYMass[0]->Add(p_SUSYMass[5]);

setTDRStyle();    
TCanvas *hist_MC[6][19];
std::ostringstream canvas; 

for(unsigned iF(0); iF<6; iF++){
  for(unsigned iC(0); iC<19; iC++){
    canvas.str("");
    canvas << "energy" << iF << "_canvas" << iC;
    hist_MC[iF][iC] = new TCanvas(canvas.str().c_str(),canvas.str().c_str(),600,600);
  }
}

TCanvas *canMass = new TCanvas("SUSY Mass","SUSY Mass",600,600);
TCanvas *canEt = new TCanvas("Photon Et","Photon Et",600,600);
canEt->cd();
int norm(1);
int LineColor[] = {1,2,3,4,6,41};
string FileList[] = {"M_{#tilde{g}}=800 to 1000","M_{#tilde{g}}=1050 to 1100", "M_{#tilde{g}}=1250 to 1300","M_{#tilde{g}}=1350 to 1400","M_{#tilde{g}}=1450 to 1500","M_{#tilde{g}}=1550"};
TLegend *legEt =  new TLegend(0.4,0.55,0.86,0.85);
std::ostringstream histname;
//.........这里部分代码省略.........

// .x PlotEffPt.C("DYJets","RecoPassedTag","RecoAllTag","RecoPassedTag","RecoAllTag",1,"NonPrompt","Prompt",2,0.4,1.0,"Id+Iso efficiency")
void PlotEffPt(TString fileName = "DYJets_TP",
	       TString numName0 = "RecoPassedTag",
	       TString denName0 = "RecoAllTag",
	       TString numName1 = "RecoPassedTag",
	       TString denName1 = "RecoAllTag",
	       bool isMuId = true,
	       TString namePrompt0 = "NonPrompt",
	       TString namePrompt1 = "Prompt",
	       int EtaBin = 2,
	       float yMin = 0.4,
	       float yMax = 1.0,
	       TString yTitle="Id+Iso efficiency") {
  // inputs
  // fileName - the name of the RooT file produced by AnalysisMacroEff (without .root extension)
  // numName0 - name of the numerator histogram for non-prompt leptons
  // denName0 - name of the denominator histogram for non-prompt leptons
  // numName1 - name of the numerator histogram for prompt leptons
  // denName1 - name of the denominator histogram for prompt leptons
  // isMuId - false : electron , true : muon
  // namePrompt0 - name of the histogram for prompt leptons
  // namePrompt1 - name of the histogram for non-prompt leptons
  // EtaBin - eta bin  0 : 0-0.8,  1 : 0.8-1.5, 2 : 1.5-2.4
  // yMin, yMax - ranges of the y axis (efficiency)
  // yTitle : title of y axis

  setTDRStyle();

  int lepSign = 0;

  TString name("ElectronPt_");
  if (isMuId)
    name = "MuonPt_";

  TString LepQ("Neg");
  if (lepSign>0)
    LepQ = "Pos";

  TString Type[2];
  Type[0] = namePrompt0;
  Type[1] = namePrompt1;


  int nEtaBins = 3;
  TString EtaBinName[3] = {"0To0p8",
			   "0p8To1p5",
			   "1p5To2p3"};

  TString EtaLeg[3] = {"#eta = [0,0.8]","#eta = [0.8,1.5]","#eta = [1.5,2.5]"};
  if (isMuId)
    EtaLeg[2] = "#eta = [1.5,2.4]";

  TFile * file = new TFile(fileName+".root");

  int nPtBins = 7;
  float ptBins[8] = {10,15,20,30,40,50,70,100};
  TString numName[2];
  TString denName[2];
  numName[0] = numName0;
  numName[1] = numName1;
  denName[0] = denName0;
  denName[1] = denName1;  

  TH1F * histos[3][2][2];
  for (int iEta=0; iEta<3; ++iEta) {
    for (int iType=0; iType<2; ++iType) {
      TH1F * histNum = (TH1F*)file->Get(name+LepQ+Type[iType]+numName[iType]+EtaBinName[iEta]);
      if (lepSign==0) {
	TH1F * histNumX = (TH1F*)file->Get(name+"Pos"+Type[iType]+numName[iType]+EtaBinName[iEta]);
	histNum->Add(histNum,histNumX);
      }
      histos[iEta][iType][0] = TH1toTH1(histNum,nPtBins,ptBins,false,"_new");
      TH1F * histDen = (TH1F*)file->Get(name+LepQ+Type[iType]+denName[iType]+EtaBinName[iEta]);
      if (lepSign==0) {
	TH1F * histDenX = (TH1F*)file->Get(name+"Pos"+Type[iType]+denName[iType]+EtaBinName[iEta]);
	histDen->Add(histDen,histDenX);
      }
      histos[iEta][iType][1] = TH1toTH1(histDen,nPtBins,ptBins,false,"_new");
    }
  }

  int color[3] = {1,2,4};
  int symbol[3] = {20,21,22};

  TGraphAsymmErrors * eff[3][2];
  for (int iEta=0; iEta<3; ++iEta) {
    for (int iType=0; iType<2; ++iType) {
      eff[iEta][iType] = new TGraphAsymmErrors();
      eff[iEta][iType]->SetLineColor(color[iEta]);
      eff[iEta][iType]->SetMarkerColor(color[iEta]);
      eff[iEta][iType]->SetLineWidth(2);
      eff[iEta][iType]->SetMarkerSize(2.5);
      eff[iEta][iType]->SetMarkerStyle(symbol[iEta]);
      if (iType==1)
	eff[iEta][iType]->SetMarkerStyle(symbol[iEta]+4);
      eff[iEta][iType]->Divide(histos[iEta][iType][0],histos[iEta][iType][1]);
    }
  }  
 
  // **************
//.........这里部分代码省略.........

void CheckDataFakes_CR_Bjet(){
  
  TH1::SetDefaultSumw2(true);
  setTDRStyle();
  
  TString param = "BJetCR_fakeparam_";

  TH1F* data_presel_hist = new TH1F("mcclosure","mcclosure", 10,0.,10.);
  data_presel_hist->GetYaxis()->SetTitle("Event");
  float ymax = 250.;
  data_presel_hist->GetYaxis()->SetRangeUser(0., ymax);
  data_presel_hist->GetXaxis()->SetTitle("");
  data_presel_hist->GetXaxis()->SetBinLabel(1,"PtvsEta");
  data_presel_hist->GetXaxis()->SetBinLabel(2,"HtvsEta");
  data_presel_hist->GetXaxis()->SetBinLabel(3,"PtvsEta_nbjet");
  data_presel_hist->GetXaxis()->SetBinLabel(4,"HtvsEta_X");
  data_presel_hist->GetXaxis()->SetBinLabel(5,"PtvsEta_X");
  data_presel_hist->GetXaxis()->SetBinLabel(6,"PtvsEta_nbjet_X");
  data_presel_hist->GetXaxis()->SetBinLabel(7,"PtvsEta_HT_nbjet_X");
  data_presel_hist->GetXaxis()->SetBinLabel(8,"PtvsEta_HT_X");
  data_presel_hist->GetXaxis()->SetBinLabel(9,"PtvsEta_HT");
  
  
  TString path_data= "/home/jalmond/Analysis/LQanalyzer/data/output/SSElectron/HNDiElectron_data_5_3_14.root";
  TFile * fdata = new TFile(path_data);
  
  TH1F* h_data_presel= (TH1F*)fdata->Get(param + "20");
  double error_up_data = GetError(h_data_presel , true);
  double error_down_data = GetError(h_data_presel , false);
  
  
  
  TCanvas* c1 = new TCanvas(("Plot"), "Plot", 1600, 1200);
  c1->cd();

  data_presel_hist->GetYaxis()->SetTitle("Event");
  data_presel_hist->Draw();

  
  vector<TString> jetpt;
  jetpt.push_back("20");
  jetpt.push_back("30");
  jetpt.push_back("40");
  jetpt.push_back("60");

  
  TString path_mc= "/home/jalmond/Analysis/LQanalyzer/data/output/SSElectron/HNDiElectron_mc_5_3_14.root";
  TFile * fmc = new TFile(path_mc);


  TString path_cf= "/home/jalmond/Analysis/LQanalyzer/data/output/SSElectron/HNDiElectron_SKchargeflip_dilep_5_3_14.root";
  TFile * fcf = new TFile(path_cf);


  TString path_np= "/home/jalmond/Analysis/LQanalyzer/data/output/SSElectron/NP/REG2/HNDiElectron_SKnonprompt_dilep_5_3_14.root";
  TFile * fnp = new TFile(path_np);

  TLegend* legendH = new TLegend(0.6,0.6, 0.9,0.9);
  legendH->SetFillColor(kWhite);
  legendH->SetTextFont(42);
  
  //h_data_presel->Draw("histsame");

  double jx[9], jxerr[9];
  double jy[9], jyerrup[9], jyerrdown[9];
  float offset = 0.1; 
  for(unsigned int bin = 0; bin < h_data_presel->GetNbinsX(); bin++){
    
    jx[bin] = bin+ offset;  jxerr[bin] = 0.;
    jy[bin] = h_data_presel->GetBinContent(bin+1);
    jyerrup[bin] = error_up_data;
    jyerrdown[bin] = error_down_data;
  }

  TGraphAsymmErrors * gd = new TGraphAsymmErrors(9, jx, jy, jxerr, jxerr,  jyerrup,jyerrdown);
  gd->SetLineWidth(2.0);
  gd->SetMarkerSize(2.);

  gd->Draw( "9sameP" );
  

  TLine *line = new TLine( data_presel_hist->GetBinLowEdge(data_presel_hist->GetXaxis()->GetFirst()),h_data_presel->GetBinContent(1),data_presel_hist->GetBinLowEdge(data_presel_hist->GetXaxis()->GetLast()+1),h_data_presel->GetBinContent(1));

  line->SetLineColor(kBlack);
  line->SetLineWidth(2);
  line->Draw();
  
  float ymaxlb = ymax * 0.6;
  float ymaxhb = ymax ;
  for(unsigned int b = 1; b < data_presel_hist->GetNbinsX(); b++){
    //Bin dividers
    if( b < 4) ymax = ymaxlb;
    else ymax = ymaxhb;
    TLine *lineb = new TLine( float(b), 0., float(b), ymax);
    lineb->SetLineColor(kRed);
    lineb->SetLineWidth(2);
    lineb->Draw();
  }

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

void controlPlots(string strigger="Jet30"){

  gStyle->SetOptStat(0);
  gROOT->ProcessLine(".L ~/tdrstyle.C");
  setTDRStyle();
  //style modifications:
  gStyle->SetErrorX(0.5);
  gStyle->SetPadTopMargin(0.08);//.005
  gStyle->SetPadBottomMargin(0.10);//0.13
  gStyle->SetMarkerStyle(1);

  TCanvas *c = new TCanvas(("ptCutComp"+strigger).c_str(),("pt cut comparison (8 vs. 10 GeV), "+strigger).c_str(),100,300,500,500);
  TCanvas *cr = new TCanvas(("ptCutCompRat"+strigger).c_str(),("pt cut comparison (8 vs. 10 GeV) Ratio, "+strigger).c_str(),500,300,500,500);
  c->cd(); 

  //string strigger="jet30"; 
  double rebin= 5.0; 
  string file1 = "PT8vsPT10/ptcut8/histo-Data-HLT_" + strigger + "U-BBCorr_HadrJetMatchPt_OpenHLT-total.root"; 
  string file2 = "PT8vsPT10/ptcut10/histo-Data-HLT_" + strigger + "U-BBCorr_HadrJetMatchPt_OpenHLT-total.root";
  string plot="Vert_dR_2b_CORR";
  string label1 = "p_{T} > 8 GeV";   
  string label2 = "p_{T} > 10 GeV";   
  //string plot="Vert_dR_2b";  

  TFile *f1 = TFile::Open(file1.c_str(),"READ");
  TFile *f2 = TFile::Open(file2.c_str(),"READ");
  
  TH1F *hp1; f1->GetObject((plot).c_str(),hp1);
  TH1F *hp2; f2->GetObject((plot).c_str(),hp2);

  hp1->Rebin(rebin);
  hp2->Rebin(rebin);
  
  hp1->SetLineColor(1);
  hp1->SetMarkerColor(1);
  hp1->SetMarkerStyle(4);
  hp2->SetLineColor(2);
  hp2->SetMarkerColor(2);
  hp2->SetMarkerStyle(20);
  hp2->SetMarkerSize(1.13);
  

  bool hp1first = false;
  if(hp1->GetMaximum()>hp2->GetMaximum()) hp1first=true;
  double max; 
  if(hp1first){
    hp1->GetXaxis()->SetTitle("#Delta R"); 
    hp1->GetYaxis()->SetTitle("d#sigma /d#Delta R"); 
    hp1->GetYaxis()->SetTitleOffset(1.9); 
    hp1->GetXaxis()->SetRangeUser(0,3.9); 
    hp1->DrawCopy("PE");
    hp2->DrawCopy("PEsames");
    hp1->DrawCopy("PEsames");
    max = hp1->GetMaximum(); 
  }
  else{
    hp2->DrawCopy("PE");
    hp2->GetXaxis()->SetTitle("#Delta R"); 
    hp2->GetYaxis()->SetTitle("d#sigma /d#Delta R"); 
    hp2->GetYaxis()->SetTitleOffset(1.9); 
    hp2->GetXaxis()->SetRangeUser(0,3.9); 
    hp1->DrawCopy("PEsames");
    max = hp2->GetMaximum(); 
  }

  TLegend *leg = new TLegend(0.6,0.65,0.89,0.84);
  leg->SetHeader(("HLT_"+strigger+"U").c_str());
  leg->SetBorderSize(0);
  leg->SetFillColor(0);
  leg->AddEntry(hp1,label1.c_str(),"PE");
  leg->AddEntry(hp2,label2.c_str(),"PE");
  leg->Draw();

  TLatex *   tex; 
  if(strigger=="Jet15") tex = new TLatex(0,max*1.26,"CMS    #sqrt{s} = 7 TeV, L = 3.1 pb^{-1}");
  else if(strigger=="Jet30") tex = new TLatex(0,max*1.17,"CMS    #sqrt{s} = 7 TeV, L = 3.1 pb^{-1}");
  else if(strigger=="Jet50") tex = new TLatex(0,max*1.13,"CMS    #sqrt{s} = 7 TeV, L = 3.1 pb^{-1}");
  tex->SetTextSize(0.040); //0.044
  tex->SetLineWidth(2);
  tex->Draw();

  gPad->RedrawAxis(); 

  cr->cd(); 
  hp2->Rebin(2); 
  hp1->Rebin(2); 
  TH1F *hrat = hp2->Clone(); 
  hrat->Divide(hp2,hp1,1,1);
  hrat->GetXaxis()->SetRangeUser(0,3.9); 
  hrat->GetXaxis()->SetTitle("#Delta R"); 
  hrat->GetYaxis()->SetTitle("ratio"); 
  hrat->GetYaxis()->SetTitleOffset(1.7); 
  hrat->SetLineColor(1);
  hrat->SetMarkerColor(1);  
  hrat->Draw(); 
  TLatex *tex2; 
  if(strigger=="Jet15") tex2 = new TLatex(0,1.53,"CMS    #sqrt{s} = 7 TeV, L = 3.1 pb^{-1}");
  else if(strigger=="Jet30") tex2 = new TLatex(0,1.187,"CMS    #sqrt{s} = 7 TeV, L = 3.1 pb^{-1}");
  else if(strigger=="Jet50") tex2 = new TLatex(0,1.143,"CMS    #sqrt{s} = 7 TeV, L = 3.1 pb^{-1}");
  tex2->SetTextSize(0.040); //0.044
//.........这里部分代码省略.........

void macroPlot( TString name, TString nameGen, const TString & nameFile1, const TString & nameFile2, const TString & title,
                const TString & resonanceType, const int rebinX, const int rebinY, const int fitType,
                const TString & outputFileName, const bool genMass) {

  gROOT->SetBatch(true);

  //Save the graphs in a file
  TFile *outputFile = new TFile(outputFileName,"RECREATE");

  setTDRStyle();

  TGraphErrors *grM_1 = fit2DProj(name, nameFile1, 100, rebinX, rebinY, fitType, outputFile, resonanceType, 0, "_1");
  TGraphErrors *grM_2 = fit2DProj(name, nameFile2, 100, rebinX, rebinY, fitType, outputFile, resonanceType, 0, "_2");
  TGraphErrors *grM_Gen = fit2DProj(nameGen, nameFile2, 100, rebinX, rebinY, fitType, outputFile, resonanceType, 0, "");

  TCanvas *c = new TCanvas(name+"_Z",name+"_Z");
  c->SetGridx();
  c->SetGridy();
    
  grM_1->SetMarkerColor(1);
  grM_2->SetMarkerColor(2);
  if(genMass)
    grM_Gen->SetMarkerColor(4);
 
  TString xAxisTitle;

  double x[2],y[2];

  if( name.Contains("Eta") ) {
    x[0]=-3; x[1]=3;       //<------useful for reso VS eta
    xAxisTitle = "muon #eta";
  }
  else if( name.Contains("PhiPlus") || name.Contains("PhiMinus") ) {
    x[0] = -3.2; x[1] = 3.2;
    xAxisTitle = "muon #phi(rad)";
  }
  else {
    x[0] = 0.; x[1] = 200;
    xAxisTitle = "muon pt(GeV)";
  }
  if( resonanceType == "JPsi" || resonanceType == "Psi2S" ) { y[0]=0.; y[1]=6.; }
  else if( resonanceType.Contains("Upsilon") ) { y[0]=8.; y[1]=12.; }
  else if( resonanceType == "Z" ) { y[0]=80; y[1]=100; }

  // This is used to have a canvas containing both histogram points
  TGraph *gr = new TGraph(2,x,y);
  gr->SetMarkerStyle(8);
  gr->SetMarkerColor(108);
  gr->GetYaxis()->SetTitle("Res Mass(GeV)   ");
  gr->GetYaxis()->SetTitleOffset(1);
  gr->GetXaxis()->SetTitle(xAxisTitle);
  gr->SetTitle(title);

  // Text for the fits
  TPaveText * paveText1 = new TPaveText(0.20,0.15,0.49,0.28,"NDC");
  paveText1->SetFillColor(0);
  paveText1->SetTextColor(1);
  paveText1->SetTextSize(0.02);
  paveText1->SetBorderSize(1);
  TPaveText * paveText2 = new TPaveText(0.59,0.15,0.88,0.28,"NDC");
  paveText2->SetFillColor(0);
  paveText2->SetTextColor(2);
  paveText2->SetTextSize(0.02);
  paveText2->SetBorderSize(1);

  TPaveText * paveText3 = new TPaveText(0.59,0.32,0.88,0.45,"NDC");
  paveText3->SetFillColor(0);
  paveText3->SetTextColor(4);
  paveText3->SetTextSize(0.02);
  paveText3->SetBorderSize(1);
  
  /*****************PARABOLIC FIT (ETA)********************/
  if( name.Contains("Eta") ) {
    std::cout << "Fitting eta" << std::endl;
    TF1 *fit1 = new TF1("fit1",onlyParabolic,-3.2,3.2,2);
    fit1->SetLineWidth(2);
    fit1->SetLineColor(1);
    if( grM_1->GetN() > 0 ) grM_1->Fit("fit1","", "", -3, 3);
    setTPaveText(fit1, paveText1);

    TF1 *fit2 = new TF1("fit2","pol0",-3.2,3.2);
    // TF1 *fit2 = new TF1("fit2",onlyParabolic,-3.2,3.2,2);
    fit2->SetLineWidth(2);
    fit2->SetLineColor(2);
    if( grM_2->GetN() > 0 ) grM_2->Fit("fit2","", "", -3, 3);
    // grM_2->Fit("fit2","R");
    setTPaveText(fit2, paveText2);
    if(genMass){
      TF1 *fit3 = new TF1("fit3",onlyParabolic,-3.2,3.2,2);
      fit3->SetLineWidth(2);
      fit3->SetLineColor(4);
      if( grM_Gen->GetN() > 0 ) grM_Gen->Fit("fit3","", "", -3, 3);
      // grM_2->Fit("fit2","R");
      setTPaveText(fit3, paveText3);
    }
  }
  /*****************SINUSOIDAL FIT (PHI)********************/
//   if( name.Contains("Phi") ) {
//     std::cout << "Fitting phi" << std::endl;
//     TF1 *fit1 = new TF1("fit1",sinusoidal,-3.2,3.2,3);
//.........这里部分代码省略.........

void Opt_mT(){
  
  
  setTDRStyle();
  //gStyle->SetPalette(1);
  
  
  TLegend *legend = new TLegend(.7, 0.7, .9, 0.9);
  legend->SetFillColor(10);
  legend->SetBorderSize(0);
  legend->SetTextSize(0.04);

  
  TFile * fnp = new TFile("/home/jalmond/Analysis/LQanalyzer/data/output/ElectronOpt/HNLowMedHighMassOptimisation_SKnonprompt_dilep_5_3_14.root");
  TFile * fcf = new TFile("/home/jalmond/Analysis/LQanalyzer/data/output/ElectronOpt/HNLowMedHighMassOptimisation_SKchargeflip_dilep_5_3_14.root");
  TFile * fmc = new TFile("/home/jalmond/Analysis/LQanalyzer/data/output/ElectronOpt/HNLowMedHighMassOptimisation_mc_5_3_14.root");
  TFile * fdata = new TFile("/home/jalmond/Analysis/LQanalyzer/data/output/ElectronOpt/HNLowMedHighMassOptimisation_data_5_3_14.root");

  TH1F* h_npcutflow= (TH1F*)fnp->Get(("HighMassOptimise"));
  TH1F* h_cfcutflow= (TH1F*)fcf->Get(("HighMassOptimise"));
  TH1F* h_mccutflow= (TH1F*)fmc->Get(("HighMassOptimise"));
  TH1F* h_datacutflow= (TH1F*)fdata->Get(("HighMassOptimise"));
  
  std::vector<TString> masses;
  masses.push_back("40");
  masses.push_back("50");
  masses.push_back("60");
  masses.push_back("70");
  masses.push_back("80");
  masses.push_back("90");
  masses.push_back("100");
  masses.push_back("125");
  masses.push_back("150");
  masses.push_back("175");
  masses.push_back("200");
  masses.push_back("225");
  masses.push_back("250");
  masses.push_back("275");
  masses.push_back("300");
  masses.push_back("325");
  masses.push_back("350");
  masses.push_back("375");
  masses.push_back("400");
  masses.push_back("500");
  masses.push_back("600");
  masses.push_back("700");

  std::map<TString, TH1*> histmap;
  std::map<TString, int> nsigmap;
  for(unsigned int imass=0; imass < masses.size(); imass++){
    TFile * fsig = new TFile(("/home/jalmond/Analysis/LQanalyzer/data/output/ElectronOpt/HNLowMedHighMassOptimisation_SKHNee" + masses.at(imass) + "_nocut_5_3_14.root").Data());

    TH1* hsig =  (TH1F*)fsig->Get(("HighMassOptimise_sig"));
    histmap[masses.at(imass)] =  hsig;

    TH1* hnsig =   (TH1F*)fsig->Get(("eventcutflow"));
    nsigmap[masses.at(imass)] = hnsig->GetBinContent(1);
    cout << masses.at(imass) << " has nsig = " << hnsig->GetBinContent(1) << endl;
  }
  
  
  int i2=0;
  for(std::map<TString, TH1*>::iterator it = histmap.begin(); it != histmap.end(); it++, i2++){
    TCanvas* c1 = new TCanvas(("Plot"+masses.at(i2)).Data(), "Plot", 800, 600);
    TH1* h_mass_mtopt = (TH1*) h_mccutflow->Clone(("A"+masses.at(i2)).Data());
  

    float nsig=50000.;
    if(it->first.Contains("40")) nsig = 100000;
    if(it->first.Contains("60")) nsig = 100000;
    if(it->first.Contains("80")) nsig = 100000;

    for(int i=1; i <   h_npcutflow->GetNbinsX(); i++){
      float sig_cut = it->second->GetBinContent(i);
      
      float tot_bkg =  h_npcutflow->GetBinContent(i) + h_cfcutflow->GetBinContent(i)  + h_mccutflow->GetBinContent(i);
      
      float sig_eff = sig_cut/nsig;
      float bkgtmp = tot_bkg + (0.28*h_npcutflow->GetBinContent(i))* (0.28*h_npcutflow->GetBinContent(i));
      float denom= 1. + sqrt(bkgtmp);
      
      
      float punzi = 0.;
      if(bkgtmp > 0.) {
	if(denom != 0.) punzi= sig_eff/denom;
      }
      //if(bkgtmp > 25.)  punzi=0.;
      //if(sig_eff < 0.0004) punzi=0.;
      if(h_npcutflow->GetBinContent(i) < 1.)  punzi=0.;
      //cout << i << " " << punzi << endl;
      h_mass_mtopt->SetBinContent(i,punzi);
      
      //if(tot_bkg > 1) cout << "bin = " << i << ": tot_bkg = " << tot_bkg << " and data  = " << h_datacutflow->GetBinContent(i)  << " signal eff = " << 100.* sig_eff<< endl;
      //cout << i << ": h_npcutflow->GetBinContent(i) = " <<  h_npcutflow->GetBinContent(i) << " h_cfcutflow->GetBinContent(i) = " << h_cfcutflow->GetBinContent(i) << " h_mccutflow->GetBinContent(i) = " << h_mccutflow->GetBinContent(i)  << " ...   data = " << h_datacutflow->GetBinContent(i) << endl;
    } 

    h_mass_mtopt->GetYaxis()->SetTitle("Punzi"); 
    h_mass_mtopt->GetXaxis()->SetTitle("cut");
    
    h_mass_mtopt->SetLineWidth(3.);
//.........这里部分代码省略.........

void plotLimit(int signal = 0){
//signal: 0 = ZPN, 1 = ZPW, 2 = ZPXW, 3 = RSG  

setTDRStyle();

//gROOT->SetStyle("Plain");
gStyle->SetOptStat(0000000000); //this clears all the boxes and crap
gStyle->SetLegendBorderSize(1);

TGraph * limit_obs = new TGraph(3);
TGraph * limit_exp = new TGraph(3);
TGraphAsymmErrors * band_exp1 = new TGraphAsymmErrors();
TGraphAsymmErrors * band_exp2 = new TGraphAsymmErrors();

TGraph * limit_obs_2 = new TGraph(3);
TGraph * limit_exp_2 = new TGraph(3);
TGraphAsymmErrors * band_exp1_2 = new TGraphAsymmErrors();
TGraphAsymmErrors * band_exp2_2 = new TGraphAsymmErrors();
TGraph * limit_obs_3 = new TGraph(3);
TGraph * limit_exp_3 = new TGraph(3);
TGraphAsymmErrors * band_exp1_3 = new TGraphAsymmErrors();
TGraphAsymmErrors * band_exp2_3 = new TGraphAsymmErrors();
TGraph * limit_obs_4 = new TGraph(3);
TGraph * limit_exp_4 = new TGraph(3);
TGraphAsymmErrors * band_exp1_4 = new TGraphAsymmErrors();
TGraphAsymmErrors * band_exp2_4 = new TGraphAsymmErrors();
TGraph *theory = new TGraph(3);
TGraph *theory2 = new TGraph(3);
TGraph *theory3 = new TGraph(3);
TGraph *theory4 = new TGraph(3);

   theory->SetPoint(0,  0.01, 4.24671);

   theory->SetPoint(1,  0.1, 42.24246);

   theory->SetPoint(2,  0.3, 122.17487);
   theory2->SetPoint(0,  0.01, 0.17980);
   theory2->SetPoint(1,  0.1, 2.00723);
   theory2->SetPoint(2,  0.3, 6.99950);
   theory3->SetPoint(0, 0.01, 0.01659*0.1);
   theory3->SetPoint(1,  0.1, 0.23030*0.1);
   theory3->SetPoint(2,  0.3, 1.03387*0.1);
   theory4->SetPoint(0, 0.01, 0.00203*0.01);
   theory4->SetPoint(1, 0.1, 0.04254*0.01);
   theory4->SetPoint(2, 0.3, 0.25352*0.01); 

 string filename = "Limits/comb_width_1TeV.txt";
 if (signal == 1) filename= "Limits/comb_width_2TeV.txt";
 if (signal == 2) filename= "Limits/comb_width_3TeV.txt";
 if (signal == 3) filename= "Limits/comb_width_4TeV.txt";

ifstream infile(filename);

double mass, exp, obs, up1, up2, dn1, dn2;
int point = 0;

while (!infile.eof()){

  //infile >> mass >> exp >>  dn2 >> up2 >> dn1 >> up1;
  //obs = exp;	

  infile >> mass >> obs >> exp >>  dn2 >> up2 >> dn1 >> up1;

 
  double sf = 1.0;
  cout.precision(3);
  //cout << mass << " & " << obs << " & " << exp << " & " << "["<<dn1<<", "<<up1<<"] & ["<<dn2<<", "<<up2<<"] \\" << endl;
  cout << mass << " & \\textbf{" << obs*sf << "} & " << dn2*sf << " & " << dn1*sf << " & \\textbf{" << exp*sf << "} & " << up1*sf << " & " << up2*sf <<  " \\\\" << endl;


  //if (mass == 2500) sf = 0.01;
  //if (mass == 3000) sf = 0.001;
  //if (mass == 3500) sf = 0.0001;
  //if (mass == 4000) sf = 0.0001;


  limit_obs->SetPoint(point, mass, obs*sf);
  limit_exp->SetPoint(point, mass, exp*sf);
  band_exp1->SetPoint(point, mass, exp*sf);
  band_exp2->SetPoint(point, mass, exp*sf);
  
  band_exp1->SetPointEYhigh(point, up1*sf - exp*sf);
  band_exp1->SetPointEYlow(point, exp*sf - dn1*sf);
  band_exp2->SetPointEYhigh(point, up2*sf - exp*sf);
  band_exp2->SetPointEYlow(point, exp*sf - dn2*sf);
  point++;

}

ifstream infile2 ("Limits/comb_width_2TeV.txt");
point = 0;
while (!infile2.eof()){

  //infile >> mass >> exp >>  dn2 >> up2 >> dn1 >> up1;
  //obs = exp;	

  infile2 >> mass >> obs >> exp >>  dn2 >> up2 >> dn1 >> up1;

 
  double sf = 1.0;
//.........这里部分代码省略.........