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

//.........这里部分代码省略.........
    double result = (ks_score>chi2_score) ? ks_score : chi2_score;
    
      href->SetNormFactor(Nevents_new);     	
      hnew->SetNormFactor(Nevents_new);
    //hnew->SetNormFactor(1);

    // ensure that the peaks of both histograms will be shown by making a dummy histogram
    float Nentries_ref = href->GetEntries();
    float Nentries_new = hnew->GetEntries();
    float XaxisMin_ref = 0, XaxisMax_ref = 0, YaxisMin_ref = 0, YaxisMax_ref = 0;
    float XaxisMin_new = 0, XaxisMax_new = 0, YaxisMin_new = 0, YaxisMax_new = 0;
    if (Nentries_ref>0) YaxisMax_ref = (href->GetMaximum()+TMath::Sqrt(href->GetMaximum()))*(Nentries_new/Nentries_ref);
    if (Nentries_new>0) YaxisMax_new = (hnew->GetMaximum()+TMath::Sqrt(hnew->GetMaximum()));

    XaxisMin_ref = href->GetXaxis()->GetXmin()>NewMin  ? href->GetXaxis()->GetXmin() : NewMin;
    XaxisMax_ref = href->GetXaxis()->GetXmax()<=NewMax ? href->GetXaxis()->GetXmax() : NewMax;
    YaxisMax_ref = (YaxisMax_ref>=YaxisMax_new) ? YaxisMax_ref : YaxisMax_new;

    if (TMath::Abs(XaxisMin_ref - XaxisMax_ref)<1E-6)
      {
	XaxisMin_ref = 0;
	XaxisMax_ref = 1;
      }
    
    TH1F *hdumb = new TH1F("hdumb","", rebinning, XaxisMin_ref, XaxisMax_ref);
    hdumb->SetMinimum(1E-1); //--For Rick
    hdumb->SetMaximum(1.05*YaxisMax_ref);
    //    if (href->GetMaximum() < hnew->GetMaximum())
    //  href->SetAxisRange(0, 1.1 * hnew->GetMaximum(), "Y");
        
    // set drawing options on the reference histogram
    href->SetStats(0);
    href->SetLineWidth(2);
    href->SetLineColor(14);
    href->SetMarkerColor(14);
    href->SetFillColor(17);
    //href->SetFillStyle(3004);
    href->GetXaxis()->SetTitle(name.c_str());
    href->GetYaxis()->SetTitle("Entries");
    href->GetYaxis()->SetTitleOffset(1.5);

    // set drawing options on the new histogram
    hnew->SetStats(0);
    hnew->SetLineWidth(2);
    hnew->SetFillStyle(3001);
    // set drawing options on the dummy histogram
    hdumb->SetStats(0);
    hdumb->GetXaxis()->SetTitle(name.c_str());
    hdumb->GetXaxis()->SetLabelSize(0.5 * hdumb->GetXaxis()->GetTitleSize());
    hdumb->GetYaxis()->SetTitle("Entries");
    hdumb->GetYaxis()->SetTitleOffset(1.5);
    hdumb->GetYaxis()->SetLabelSize(0.5 * hdumb->GetXaxis()->GetTitleSize());
    
    stringstream ss_title;
    ss_title.precision(5);
    if (ks_score>chi2_score)
      ss_title << "KS Score = " << ks_score;
    else
      ss_title << "Chi^2 Score = " << chi2_score;
    TText canvas_title(0.1,0.97,ss_title.str().c_str());


    // determine if test is a "pass" or a "fail"
    if (result <= threshold) {

      canvas_title.SetTextColor(kRed);

//.........这里部分代码省略.........
         << delim << setw(left) << Form("%.3f %s %.4f",effdata2,pm,effdata2err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effdata3,pm,effdata3err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effdata4,pm,effdata4err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effdata5,pm,effdata5err) << setw(4)
         << delimend << endl;

    cout << delimstart << setw(10) << "mc" << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effmc1,pm,effmc1err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effmc2,pm,effmc2err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effmc3,pm,effmc3err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effmc4,pm,effmc4err) << setw(4)
         << delim << setw(left) << Form("%.3f %s %.4f",effmc5,pm,effmc5err) << setw(4)
         << delimend << endl;

    cout << delimstart << setw(10) << "data/mc" << setw(4)
         << delim << setw(left) << Form("%.2f %s %.2f",ratio1,pm,ratio1err) << setw(4)
         << delim << setw(left) << Form("%.2f %s %.2f",ratio2,pm,ratio2err) << setw(4)
         << delim << setw(left) << Form("%.2f %s %.2f",ratio3,pm,ratio3err) << setw(4)
         << delim << setw(left) << Form("%.2f %s %.2f",ratio4,pm,ratio4err) << setw(4)
         << delim << setw(left) << Form("%.2f %s %.2f",ratio5,pm,ratio5err) << setw(4)
         << delimend << endl;

    //--------------------------------------
    // draw stuff
    //--------------------------------------

    hdata->Scale(1.0/hdata->Integral());
    hmc->Scale(1.0/hmc->Integral());

    if( log ) hmc->GetYaxis()->SetRangeUser(0.0001,5);
    else      hmc->GetYaxis()->SetRangeUser(0.0,1);

    hmc->GetXaxis()->SetTitle(xtitle);
    hmc->SetLineColor(2);
    hmc->SetMarkerColor(2);
    hmc->DrawNormalized("hist");
    hmc->DrawNormalized("sameE1");
    hdata->SetLineColor(4);
    hdata->SetMarkerColor(4);
    hdata->Draw("sameE1");

    grdata->SetLineColor(6);
    grmc->SetLineColor(7);
    //grdata->Draw("sameP");
    //grmc->Draw("sameP");

    TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
    leg->AddEntry(hdata , "data" , "lp");
    leg->AddEntry(hmc   , "MC"   , "lp");
    leg->SetBorderSize(0);
    leg->SetFillColor(0);
    leg->Draw();

    TLatex *t = new TLatex();
    t->SetNDC();

    if( TString(plottitle).Contains("el") ) t->DrawLatex(0.6,0.6,"electrons");
    if( TString(plottitle).Contains("mu") ) t->DrawLatex(0.6,0.6,"muons");

    if( TString(plottitle).Contains("0j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 0");
    if( TString(plottitle).Contains("1j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 1");
    if( TString(plottitle).Contains("2j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 2");
    if( TString(plottitle).Contains("3j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 3");
    if( TString(plottitle).Contains("4j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 4");

    //--------------------------------------

void frac_fit(){
  double Ndata []  = {315239.0, 40524.1, 8625.2};
  double NdataE[]  = {   579.8,   202.8,   95.4};
  double NW []     = {258150.0,    0.38,   14.8};
  double NWE[]     = {   502.8,    0.38,    3.1};
  double NZ []     = { 19473.9, 37494.8,  181.4};
  double NZE[]     = {    61.7,    75.4,    5.2};
  double Ntt []    = { 33156.3,  1637.8, 8214.0};
  double NttE[]    = {    32.9,     8.2,   18.3};
  double Nother [] = { 20934.8,   606.9,  470.8};
  double NotherE[] = {   494.4,     5.3,   10.7};

  TH1F *W     = new TH1F("W"    ,"W"    , 3, 0.5, 3.5);
  TH1F *Z     = new TH1F("Z"    ,"Z"    , 3, 0.5, 3.5);
  TH1F *tt    = new TH1F("tt"   ,"tt"   , 3, 0.5, 3.5);
  TH1F *other = new TH1F("other","other", 3, 0.5, 3.5);
  TH1F *data  = new TH1F("data" ,"data" , 3, 0.5, 3.5);
	
  TH1F *weightW  = new TH1F("W_weight" , "W_weight" , 3, 0.5, 3.5);	
  TH1F *weightZ  = new TH1F("Z_weight" , "Z_weight" , 3, 0.5, 3.5);	
  TH1F *weighttt = new TH1F("tt_weight", "tt_weight", 3, 0.5, 3.5);	

  for(int i=0; i<3; i++){
    double Nw_prime     = TMath::Power(NW[i]/NWE[i]        , 2); 
    double NZ_prime     = TMath::Power(NZ[i]/NZE[i]        , 2); 
    double Ntt_prime    = TMath::Power(Ntt[i]/NttE[i]      , 2); 
    //double Nother_prime = TMath::Power(Nother[i]/NotherE[i], 2); 
    W ->SetBinContent(i+1, Nw_prime );
    Z ->SetBinContent(i+1, NZ_prime );
    tt->SetBinContent(i+1, Ntt_prime);

    weightW -> SetBinContent(i+1, NW[i] /W ->GetBinContent(i+1));   
    weightZ -> SetBinContent(i+1, NZ[i] /Z ->GetBinContent(i+1));   
    weighttt-> SetBinContent(i+1, Ntt[i]/tt->GetBinContent(i+1));   
    //    other   -> SetBinContent(i+1, Nother[i] );
    data    -> SetBinContent(i+1, Ndata[i]);
  }

  W ->SetLineColor(kRed  );
  Z ->SetLineColor(kBlue );
  tt->SetLineColor(kGreen);

	for(int i=1; i<=3; i++){
	//	cout << "data: " << data->GetBinContent(i) << " W " << W->GetBinContent(i) << " Z " << Z->GetBinContent(i) << " tt " << tt->GetBinContent(i) << endl; 	
	//	cout << "data weight=1, W weight "  << weightW->GetBinContent(i) << " Z weight: " << weightZ->GetBinContent(i) << " weight tt " << weighttt->GetBinContent(i) << endl; 
	}
	
  TObjArray *mc= new TObjArray(4);
  mc->Add(W);
  mc->Add(Z);
  mc->Add(tt);
  //  mc->Add(other);
  TFractionFitter *fit = new TFractionFitter(data,mc);
  //fit->GetFitter()->FixParameter(3);
	fit->SetWeight(0,weightW);
	fit->SetWeight(1,weightZ);
	fit->SetWeight(2,weighttt);
	//fit->Constrain(0,0,2);
	//fit->Constrain(1,0,2);
	//fit->Constrain(2,0,2);

  fit->Fit();
  double R[3]; 

  for(int i=0; i<3; i++){
    double x=0;
    double xE=0;
    fit->GetResult(i,x,xE);
    double N=0; 
    if (i==0) N = NW[i] ; 
    if (i==1) N = NZ[i] ;
    if (i==2) N = Ntt[i];
	  cout << "x " << x << endl; 
	  cout << " initial fraction: " << Ndata[i]/N << endl; 
    R[i]=x*(Ndata[i]/N);
    cout << "R: " << R[i] << " +/- " << (xE/x)*R[i] << endl; 

  }
}

void bbtt_upg_em(std::string var,int nbins, double xmin, double xmax,std::string xtitle, std::string ytitle, double sigscale=1)
{

  TFile *outDC = new TFile("hh_em_inputs.root","RECREATE");

  SetStyle(); gStyle->SetLineStyleString(11,"20 10");
  TH1::SetDefaultSumw2(1);
 
  //std::string dir = "/data/blue/Bacon/029a/Upgrade/merged_talk_jun30/";
  std::string dir = "/afs/cern.ch/work/j/jlawhorn/public/ntuples/";
  
  std::stringstream scale; scale << sigscale;
  
  //Cut definitions
  double luminosity = 3000;
  std::stringstream lumi; lumi << luminosity;
  std::string objcut = "(tauCat1==3 && tauCat2==2 && ptTau1>20 && ptTau2>20 && tauIso2<0.4 && tauIso1<0.4 && (bTag1==2||bTag1==3||bTag1==6||bTag1==7) && (bTag2==1||bTag2==3||bTag2==6||bTag2==7) && ptB1>20 && ptB2>20 && sqrt( (etaTau1-etaTau2)**2 + (phiTau1-phiTau2)**2 )>0.4)";
  std::string jetcut = objcut+"*(mTT>20 && mTT<90)*(mBB1>70 && mBB1<140)*(mt2pileup>100)*(bdtVal>-0.05)";
  //signal region
  std::string mccut = jetcut+"*eventWeight*"+lumi.str();
  std::string sigcut = jetcut+"*eventWeight*"+lumi.str();
  std::string zjetcut = jetcut+"*eventWeight*(eventType==4)*"+lumi.str();
  std::string wjetcut = jetcut+"*eventWeight*(eventType==3)*"+lumi.str();
  std::string ewkcut = jetcut+"*eventWeight*(eventType!=1)*"+lumi.str();
  //--------------------------------------------------------------------------
  
  //Get the trees
  TTree *hhtree = load(dir+"HHToTTBB_14TeV.root"); 
  TTree *tttree = load(dir+"tt.root"); 
  //TTree *vbfhtree = load(dir+"VBFToTT_14TeV_phase2.root");
  //TTree *gfhtree = load(dir+"H.root");
  //TTree *vjettree = load(dir+"Vjets.root");
  //TTree *ewktree = load(dir+"diboson.root");
  
  //-------------------------------------------------------------------------
  
  //Get histograms
  TCanvas *canv0 = MakeCanvas("canv", "histograms", 600, 600);
  canv0->cd();
  std::string vardraw;
  /*  TH1F *Ztt = new TH1F("DY","",nbins,xmin,xmax);
  vardraw = var+">>"+"DY";
  vjettree->Draw(vardraw.c_str(),zjetcut.c_str());
  InitHist(Ztt  , xtitle.c_str(), ytitle.c_str(), TColor::GetColor(248,206,104), 1001);*/
  TH1F *ttbar = new TH1F("TTbar","",nbins,xmin,xmax);
  vardraw = var+">>"+"TTbar";
  tttree->Draw(vardraw.c_str(),mccut.c_str());
  InitHist(ttbar, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(155,152,204), 1001);
  /*  TH1F *wjets = new TH1F("Wjets","",nbins,xmin,xmax);
  vardraw = var+">>"+"Wjets";
  vjettree->Draw(vardraw.c_str(),wjetcut.c_str());
  InitHist(wjets, xtitle.c_str(), ytitle.c_str(), TColor::GetColor(222,90,106), 1001);
  TH1F *ewk = new TH1F("Ewk","",nbins,xmin,xmax);
  vardraw = var+">>"+"Ewk";
  ewktree->Draw(vardraw.c_str(),ewkcut.c_str());
  InitHist(ewk, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  TH1F *vbfh = new TH1F("VBFH","",nbins,xmin,xmax);
  vardraw = var+">>"+"VBFH";
  vbfhtree->Draw(vardraw.c_str(),mccut.c_str());
  InitHist(vbfh, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(250,202,255), 1001);
  TH1F *ggh = new TH1F("GGH","",nbins,xmin,xmax);
  vardraw = var+">>"+"GGH";
  gfhtree->Draw(vardraw.c_str(),mccut.c_str());
  InitHist(ggh, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(250,202,255), 1001);*/
  TH1F *smhh = new TH1F("SMhh","",nbins,xmin,xmax);
  vardraw = var+">>"+"SMhh";
  hhtree->Draw(vardraw.c_str(),sigcut.c_str());
  InitSignal(smhh);
  smhh->SetLineColor(kBlack);
  delete canv0;
  //---------------------------------------------------------------------------
  //Print out the yields
  Double_t error=999;
  //ofstream outfile;
  //outfile.open("yields.txt");
  //outfile << "Yields for the signal region." << std::endl;
  cout << "SM hh   "  << smhh->IntegralAndError(0,smhh->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  /*  outfile << "SM h   "  << smhh->IntegralAndError(0,smhh->GetNbinsX(),error) << "+/-" << error << endl;
      outfile << "Ztt    "  << Ztt->IntegralAndError(0,Ztt->GetNbinsX(),error) << "+/-" << error << endl;*/
  cout << "ttbar    "  << ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error) << "+/-";
  cout << error << endl; error=999;
  /*  outfile << "ewk    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
      outfile << "wjets    "  << wjets->IntegralAndError(0,wjets->GetNbinsX(),error) << "+/-" << error << endl;*/
  //--------------------------------------------------------------------------
  //outfile.close();
  outDC->cd();
  TDirectory* lTD = outDC->mkdir("emu");
  outDC->cd(lTD->GetPath());
  ttbar->SetName("data_obs");
  ttbar->SetTitle("data_obs");
  ttbar->Write();
  /*Ztt->SetName("ZTT");
  Ztt->SetTitle("ZTT");
  Ztt->Write();*/
  ttbar->SetName("TT");
  ttbar->SetTitle("TT");
  ttbar->Write();
  /*wjets->SetName("W");
  wjets->SetTitle("W");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      hDeltaPhiBack12->Fill((Phi1-Phi2)*TMath::RadToDeg());
      hThetaRelBack12->Fill(thetarel12);
    }

    if(TMath::Abs(DeltaT13-corr13) < 500){
      hDeltaTheta13->Fill((Theta1-Theta3)*TMath::RadToDeg());
      hDeltaPhi13->Fill((Phi1-Phi3)*TMath::RadToDeg());
      hThetaRel13->Fill(thetarel13);
    }
    else if(TMath::Abs(DeltaT13-corr13) > 1000 && TMath::Abs(DeltaT13-corr13) < 6000){
      hDeltaThetaBack13->Fill((Theta1-Theta3)*TMath::RadToDeg());
      hDeltaPhiBack13->Fill((Phi1-Phi3)*TMath::RadToDeg());
      hThetaRelBack13->Fill(thetarel13);
    }
  }
  
  h->SetStats(0);

  hDeltaThetaBack12->Sumw2();
  hDeltaPhiBack12->Sumw2();
  hThetaRelBack12->Sumw2();
  hDeltaThetaBack12->Scale(0.1);
  hDeltaPhiBack12->Scale(0.1);
  hThetaRelBack12->Scale(0.1);
  hDeltaThetaBack13->Sumw2();
  hDeltaPhiBack13->Sumw2();
  hThetaRelBack13->Sumw2();
  hDeltaThetaBack13->Scale(0.1);
  hDeltaPhiBack13->Scale(0.1);
  hThetaRelBack13->Scale(0.1);

  Float_t val,eval;
  TCanvas *c1=new TCanvas();
  TF1 *ff = new TF1("ff","[0]*[4]/[2]/sqrt(2*TMath::Pi())*TMath::Exp(-(x-[1])*(x-[1])*0.5/[2]/[2]) + [3]*[4]/6/[2]");
  ff->SetParName(0,"signal");
  ff->SetParName(1,"mean");
  ff->SetParName(2,"sigma");
  ff->SetParName(3,"background");
  ff->SetParName(4,"bin width");
  ff->SetParameter(0,42369);
  ff->SetParameter(1,0);
  ff->SetParLimits(2,10,1000);
  ff->SetParameter(2,150); // fix witdh if needed
  ff->SetParameter(3,319);
  ff->FixParameter(4,20000./nbint); // bin width
  
  ff->SetNpx(1000);
  
  h->Fit(ff);
  
  val = ff->GetParameter(2);
  eval = ff->GetParError(2);
  
  printf("significance = %f\n",ff->GetParameter(0)/sqrt(ff->GetParameter(0) + ff->GetParameter(3)));

  h->Draw();
  
  TF1 *func1 = (TF1 *)  h->GetListOfFunctions()->At(0);
  
  func1->SetLineColor(2);
  h->SetLineColor(4);
  
  TPaveText *text = new TPaveText(1500,(h->GetMinimum()+(h->GetMaximum()-h->GetMinimum())*0.6),9500,h->GetMaximum());
  text->SetFillColor(0);
  sprintf(title,"width = %5.1f #pm %5.1f",func1->GetParameter(2),func1->GetParError(2));
  text->AddText(title);
  sprintf(title,"signal (S) = %5.1f #pm %5.1f",func1->GetParameter(0),func1->GetParError(0));
  text->AddText(title);
  sprintf(title,"background (B) (3#sigma) = %5.1f #pm %5.1f",func1->GetParameter(3),func1->GetParError(3));
  text->AddText(title);
  sprintf(title,"significance (S/#sqrt{S+B}) = %5.1f",func1->GetParameter(0)/sqrt(func1->GetParameter(0)+func1->GetParameter(3)));
  text->AddText(title);
  
  text->SetFillStyle(0);
  text->SetBorderSize(0);
  
  text->Draw("SAME");
  
  printf("n_day = %f\nn_dayGR = %f\n",nsec*1./86400,nsecGR*1./86400);

  text->AddText(Form("rate = %f #pm %f per day",func1->GetParameter(0)*86400/nsecGR,func1->GetParError(0)*86400/nsecGR));

  TFile *fo = new TFile("output-SAVO-010203.root","RECREATE");
  h->Write();
  h2->Write();
  hDeltaTheta12->Write();
  hDeltaPhi12->Write();
  hThetaRel12->Write();
  hDeltaThetaBack12->Write();
  hDeltaPhiBack12->Write();
  hThetaRelBack12->Write();
  hDeltaTheta13->Write();
  hDeltaPhi13->Write();
  hThetaRel13->Write();
  hDeltaThetaBack13->Write();
  hDeltaPhiBack13->Write();
  hThetaRelBack13->Write();
  fo->Close();
  
}

//____________________________________
void rs_bernsteinCorrection(){

  // set range of observable
  Double_t lowRange = -1, highRange =5;

  // make a RooRealVar for the observable
  RooRealVar x("x", "x", lowRange, highRange);

  // true model
  RooGaussian narrow("narrow","",x,RooConst(0.), RooConst(.8));
  RooGaussian wide("wide","",x,RooConst(0.), RooConst(2.));
  RooAddPdf reality("reality","",RooArgList(narrow, wide), RooConst(0.8));

  RooDataSet* data = reality.generate(x,1000);

  // nominal model
  RooRealVar sigma("sigma","",1.,0,10);
  RooGaussian nominal("nominal","",x,RooConst(0.), sigma);

  RooWorkspace* wks = new RooWorkspace("myWorksspace");

  wks->import(*data, Rename("data"));
  wks->import(nominal);

  // The tolerance sets the probability to add an unecessary term.
  // lower tolerance will add fewer terms, while higher tolerance
  // will add more terms and provide a more flexible function.
  Double_t tolerance = 0.05; 
  BernsteinCorrection bernsteinCorrection(tolerance);
  Int_t degree = bernsteinCorrection.ImportCorrectedPdf(wks,"nominal","x","data");

  cout << " Correction based on Bernstein Poly of degree " << degree << endl;

  RooPlot* frame = x.frame();
  data->plotOn(frame);
  // plot the best fit nominal model in blue
  nominal.fitTo(*data,PrintLevel(-1));
  nominal.plotOn(frame);

  // plot the best fit corrected model in red
  RooAbsPdf* corrected = wks->pdf("corrected");  
  corrected->fitTo(*data,PrintLevel(-1));
  corrected->plotOn(frame,LineColor(kRed));

  // plot the correction term (* norm constant) in dashed green
  // should make norm constant just be 1, not depend on binning of data
  RooAbsPdf* poly = wks->pdf("poly");  
  poly->plotOn(frame,LineColor(kGreen), LineStyle(kDashed));
  
  // this is a switch to check the sampling distribution
  // of -2 log LR for two comparisons:
  // the first is for n-1 vs. n degree polynomial corrections
  // the second is for n vs. n+1 degree polynomial corrections
  // Here we choose n to be the one chosen by the tolerance
  // critereon above, eg. n = "degree" in the code.
  // Setting this to true is takes about 10 min.
  bool checkSamplingDist = false;

  TCanvas* c1 = new TCanvas();
  if(checkSamplingDist) {
    c1->Divide(1,2);
    c1->cd(1);
  }
  frame->Draw();

  if(checkSamplingDist) {
    // check sampling dist
    TH1F* samplingDist = new TH1F("samplingDist","",20,0,10);
    TH1F* samplingDistExtra = new TH1F("samplingDistExtra","",20,0,10);
    int numToyMC = 1000;
    bernsteinCorrection.CreateQSamplingDist(wks,"nominal","x","data",samplingDist, samplingDistExtra, degree,numToyMC);
    
    c1->cd(2);
    samplingDistExtra->SetLineColor(kRed);
    samplingDistExtra->Draw();
    samplingDist->Draw("same");
  }
}

//.........这里部分代码省略.........
    model->SetBinContent(ibin+1, model->GetBinContent(ibin+1)>0 ? model->GetBinContent(ibin+1)*model->GetBinWidth(ibin+1) : 0.01);
    model->SetBinError  (ibin+1, CONVERVATIVE_CHI2 ? 0. : model->GetBinError  (ibin+1)*model->GetBinWidth(ibin+1));
    test1->SetBinContent(ibin+1, test1->GetBinContent(ibin+1)*test1->GetBinWidth(ibin+1));
    test1->SetBinError  (ibin+1, test1->GetBinError  (ibin+1)*test1->GetBinWidth(ibin+1));
  }
  double chi2prob = test1->Chi2Test      (model,"PUW");        std::cout << "chi2prob:" << chi2prob << std::endl;
  double chi2ndof = test1->Chi2Test      (model,"CHI2/NDFUW"); std::cout << "chi2ndf :" << chi2ndof << std::endl;
  double ksprob   = test1->KolmogorovTest(model);              std::cout << "ksprob  :" << ksprob   << std::endl;
  double ksprobpe = test1->KolmogorovTest(model,"DX");         std::cout << "ksprobpe:" << ksprobpe << std::endl;  

  std::vector<double> edges;
  TH1F* zero = (TH1F*)ref->Clone("zero"); zero->Clear();
  TH1F* rat1 = (TH1F*)data->Clone("rat1"); 
  for(int ibin=0; ibin<rat1->GetNbinsX(); ++ibin){
    rat1->SetBinContent(ibin+1, Ztt->GetBinContent(ibin+1)>0 ? data->GetBinContent(ibin+1)/Ztt->GetBinContent(ibin+1) : 0);
    rat1->SetBinError  (ibin+1, Ztt->GetBinContent(ibin+1)>0 ? data->GetBinError  (ibin+1)/Ztt->GetBinContent(ibin+1) : 0);
    zero->SetBinContent(ibin+1, 0.);
    zero->SetBinError  (ibin+1, Ztt->GetBinContent(ibin+1)>0 ? Ztt ->GetBinError  (ibin+1)/Ztt->GetBinContent(ibin+1) : 0);
  }
  for(int ibin=0; ibin<rat1->GetNbinsX(); ++ibin){
    if(rat1->GetBinContent(ibin+1)>0){
      edges.push_back(TMath::Abs(rat1->GetBinContent(ibin+1)-1.)+TMath::Abs(rat1->GetBinError(ibin+1)));
      // catch cases of 0 bins, which would lead to 0-alpha*0-1
      rat1->SetBinContent(ibin+1, rat1->GetBinContent(ibin+1)-1.);
    }
  }
  float range = 0.1;
  std::sort(edges.begin(), edges.end());
  if (edges[edges.size()-2]>0.1) { range = 0.2; }
  if (edges[edges.size()-2]>0.2) { range = 0.5; }
  if (edges[edges.size()-2]>0.5) { range = 1.0; }
  if (edges[edges.size()-2]>1.0) { range = 1.5; }
  if (edges[edges.size()-2]>1.5) { range = 2.0; }
  rat1->SetLineColor(kBlack);
  rat1->SetFillColor(kGray );
  rat1->SetMaximum(+range);
  rat1->SetMinimum(-range);
  rat1->GetYaxis()->CenterTitle();
  rat1->GetYaxis()->SetTitle("#bf{Data/MC-1}");
  rat1->GetXaxis()->SetTitle("#bf{m_{#tau#tau} [GeV]}");
  rat1->Draw();
  zero->SetFillStyle(  3013);
  zero->SetFillColor(kBlack);
  zero->SetLineColor(kBlack);
  zero->SetMarkerSize(0.1);
  zero->Draw("e2histsame");
  canv0->RedrawAxis();

  TPaveText* stat1 = new TPaveText(0.20, 0.76+0.061, 0.32, 0.76+0.161, "NDC");
  stat1->SetBorderSize(   0 );
  stat1->SetFillStyle(    0 );
  stat1->SetTextAlign(   12 );
  stat1->SetTextSize ( 0.05 );
  stat1->SetTextColor(    1 );
  stat1->SetTextFont (   62 );
  stat1->AddText(TString::Format("#chi^{2}/ndf=%.3f,  P(#chi^{2})=%.3f", chi2ndof, chi2prob));
  //stat1->AddText(TString::Format("#chi^{2}/ndf=%.3f,  P(#chi^{2})=%.3f, P(KS)=%.3f", chi2ndof, chi2prob, ksprob));
  stat1->Draw();

  /*
    Ratio After fit over Prefit
  */
  TCanvas *canv1 = MakeCanvas("canv1", "histograms", 600, 400);
  canv1->SetGridx();
  canv1->SetGridy();
  canv1->cd();

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


  tree_bkg1->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg1(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg1 =(TH1F*)gPad->GetPrimitive("hbkg1");
  tree_bkg2->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg2(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg2 =(TH1F*)gPad->GetPrimitive("hbkg2");
  tree_bkg3->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg3(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && genmatch && DPHICut));
  TH1F *hbkg3 =(TH1F*)gPad->GetPrimitive("hbkg3");
     
tree_bkg4->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg4(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut)); 
  TH1F *hbkg4 =(TH1F*)gPad->GetPrimitive("hbkg4");
  tree_bkg5->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg5(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && genmatch && DPHICut));
  TH1F *hbkg5 =(TH1F*)gPad->GetPrimitive("hbkg5");
  tree_bkg6->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg6(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg6 =(TH1F*)gPad->GetPrimitive("hbkg6");  
  tree_bkg7->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg7(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg7 =(TH1F*)gPad->GetPrimitive("hbkg7");  
  tree_bkg8->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg8(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg8 =(TH1F*)gPad->GetPrimitive("hbkg8");  
  tree_bkg9->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg9(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut)); 
  TH1F *hbkg9 =(TH1F*)gPad->GetPrimitive("hbkg9");  
  tree_bkg10->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg10(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg10 =(TH1F*)gPad->GetPrimitive("hbkg10");  
 tree_bkg11->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg11(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg11 =(TH1F*)gPad->GetPrimitive("hbkg11");  
 tree_bkg12->Draw("(t1pfmetPhiCorr(t1pfmet,t1pfmetPhi,t1pfmetSumEt,0))>>hbkg12(10,-4,4)","weight*(mgg<115||mgg>135)"*(mggmax && mggmin && pt1Cut && pt2Cut && METCut && eveto  && DPHICut));
  TH1F *hbkg12 =(TH1F*)gPad->GetPrimitive("hbkg12");  



 

  hdata->SetMarkerColor(kBlack);
  hdata->SetLineColor(kBlack);
  hdata->SetMarkerStyle(20);
  
  h1->SetLineColor(kRed+3);
  h2->SetLineColor(kRed+1);
  h3->SetLineColor(kRed);
  h4->SetLineColor(kPink+2);
  h5->SetLineColor(kPink+4);
 h6->SetLineColor(kPink+7);
 h7->SetLineColor(kMagenta+2);
 h1->SetLineWidth(2);
 h2->SetLineWidth(2);
 h3->SetLineWidth(2);
 h4->SetLineWidth(2);
 h5->SetLineWidth(2);
 h6->SetLineWidth(2);
 h7->SetLineWidth(2);


  
  
  THStack *hs=new THStack("hs","");
  hbkg7->SetFillColor(kGreen+2);
  hbkg6->SetFillColor(kGreen);
  hbkg8->SetFillColor(kYellow);
  hbkg4->SetFillColor(kOrange);
  hbkg9->SetFillColor(kOrange+7);
  hbkg10->SetFillColor(kOrange+4);
  hbkg11->SetFillColor(kCyan);
  hbkg12->SetFillColor(kCyan+1);
  hbkg5->SetFillColor(kBlue+2);
  hbkg2->SetFillColor(kBlue);  
  hbkg3->SetFillColor(kMagenta-2);

void plotDavid(TH1F* bkgT, TH1F* sigT, TH1F* dataT, std::string name){
  
  gROOT->SetBatch();
  system("mkdir -p plots/ada/david");
  system("mkdir -p plots/grad/david");
  system("mkdir -p plots/ada/diff");
  system("mkdir -p plots/grad/diff");

  std::string bdt;
  TString str = dataT->GetName();
  if (str.Contains("ada")) bdt="ada";
  else if (str.Contains("grad")) bdt="grad";
  else std::cout << "Error find BDT type" << std::endl;
  assert (str.Contains("ada") || str.Contains("grad"));

  if (bkgT->GetNbinsX() != sigT->GetNbinsX() || sigT->GetNbinsX() != dataT->GetNbinsX()) std::cout << "Plot problem: calling plot for histograms with different number of bins" << std::endl;
  assert(bkgT->GetNbinsX() == sigT->GetNbinsX() || sigT->GetNbinsX() == dataT->GetNbinsX());
  
  TH1F *bkg = linearBin(bkgT);
  TH1F *sig2 = linearBin(sigT);
  TH1F *sig5 = (TH1F*)sig2->Clone();
  TH1F *sig10 = (TH1F*)sig2->Clone();
  TH1F *data = linearBin(dataT);
  TH1F *diff = (TH1F*)data->Clone();
  diff->Add(bkg,-1);

  gStyle->SetOptStat(0);
  gROOT->SetStyle("Plain");
  gROOT->ForceStyle();

  TCanvas *canv = new TCanvas();
  bkg->SetLineColor(kBlue);
  bkg->SetFillColor(kBlue-9);
  sig2->SetLineColor(kRed);
  sig2->SetLineStyle(3);
  sig2->Scale(1.);
  sig5->SetLineColor(kRed);
  sig5->SetLineStyle(7);
  sig5->Scale(3.);
  sig10->SetLineColor(kRed);
  sig10->Scale(5.);
  data->SetMarkerStyle(20);
  bkg->SetTitle(("Background, data and signal distributions for "+bdt+" "+name).c_str());
  sig2->SetTitle(("Background, data and signal distributions for "+bdt+" "+name).c_str());
  data->SetTitle(("Background, data and signal distributions for "+bdt+" "+name).c_str());
  bkg->GetXaxis()->SetTitle("");
  sig2->GetXaxis()->SetTitle("");
  data->GetXaxis()->SetTitle("");
  bkg->GetYaxis()->SetTitle("Events");
  sig2->GetYaxis()->SetTitle("Events");
  data->GetYaxis()->SetTitle("Events");

  TLegend *leg = new TLegend(0.45,0.6,0.85,0.85);
  leg->SetLineColor(0);
  leg->SetFillColor(0);
  leg->AddEntry(bkg,"Background","f");
  leg->AddEntry(data,"Data","lep");
  leg->AddEntry(sig10,"Signal (1, 3, 5 #times SM)","l");
  TPaveText *txt = new TPaveText(0.2,0.1,0.4,0.35,"NDC");
  txt->SetFillColor(0);
  txt->SetLineColor(0);
  txt->AddText("#int L = 4.76 fb^{-1}");

  bkg->GetYaxis()->SetRangeUser(1.0,2.*(data->GetMaximum()));

  bkg->Draw("e2");
  sig2->Draw("same hist");
  sig5->Draw("same hist");
  sig10->Draw("same hist");
  data->Draw("same e");
  leg->Draw("same");
  txt->Draw("same");

  canv->SetLogy();
  canv->Print(("plots/"+bdt+"/david/"+name+".png").c_str(),"png");
  canv->SetLogy(false);
  canv->Clear();
  TLegend *leg2 = new TLegend(0.45,0.6,0.85,0.85);
  leg2->SetFillColor(0);
  leg2->SetLineColor(0);
  sig10->GetYaxis()->SetRangeUser((-1*sig10->GetMaximum())+10,sig10->GetMaximum()+20);
  sig10->GetXaxis()->SetTitle("");
  sig10->GetYaxis()->SetTitle("Events");
  sig10->SetTitle(Form("Data, background difference compared to signal %s %s",bdt.c_str(),name.c_str()));
  diff->SetMarkerStyle(20);
  diff->GetXaxis()->SetTitle("");
  diff->GetYaxis()->SetTitle("Events");
  leg2->AddEntry(diff,"Data - background model","lep");
  leg2->AddEntry(sig10,"Signal (1, 3, 5 #times SM)","l");
  TF1 *line = new TF1("line","0.0",0.0,sig10->GetNbinsX()+1);
  line->SetLineColor(kBlue);
  sig10->Draw("hist");
  sig5->Draw("same hist");
  sig2->Draw("same hist");
  line->Draw("same");
  diff->Draw("p same");
  leg2->Draw("same");
  txt->Draw("same");
  canv->Print(("plots/"+bdt+"/diff/"+name+".png").c_str(),"png");

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

void drawFigure4a(TString var="Njets") 
{
  gInterpreter->ExecuteMacro("GoodStyle.C");
  gROOT->LoadMacro("tdrstyle.C");                                                                                                                           
  setTDRStyle();    

  TString variable = "p_{T,max}^{#font[12]{l}}";
  std::ostringstream filename, madname;
  filename<<"rootfiles/all_unfolding_"<<var<<".root";
  TFile* file = new TFile(filename.str().c_str(), "read");
  madname<<"hGenXs"<<var<<"_1";
  std::cout<<madname.str().c_str()<<std::endl;
  TH1F* xsValue          = (TH1F*)(file->Get("hComb_diff")->Clone("xsValue"));
  TH1F* xsValue_Madgraph = (TH1F*)(file->Get(madname.str().c_str())->Clone("xsValue_Madgraph"));
  TH1F* xsValue_MCnlo    = (TH1F*)(file->Get("mcfm_tot")->Clone("xsValue_MCnlo"));

  // Set the data errors- I don't need this because I already have complete error in my plot
  //----------------------------------------------------------------------------
  
  // Data cosmetics
  //----------------------------------------------------------------------------
  xsValue->SetLineWidth(1);
  xsValue->SetMarkerSize(_msize);
  xsValue->SetMarkerStyle(kFullCircle);
  xsValue->SetMarkerColor(kBlack);
  xsValue->SetLineColor(kBlack);
  xsValue->SetFillStyle(1001);
  xsValue->SetFillColor(kWhite);
  // Madgraph cosmetics
  //----------------------------------------------------------------------------
  xsValue_Madgraph->SetFillColor(kOrange);
  //xsValue_Madgraph->SetFillColor(kWhite);
  xsValue_Madgraph->SetFillStyle(1001);
  xsValue_Madgraph->SetLineColor(kOrange+7);
  xsValue_Madgraph->SetLineWidth(1);
  xsValue_Madgraph->SetMarkerColor(kOrange+7);
  xsValue_Madgraph->SetMarkerSize(_msize);
  xsValue_Madgraph->SetMarkerStyle(21);


  // MCNLO cosmetics
  //----------------------------------------------------------------------------
  xsValue_MCnlo->SetFillColor(kAzure-9);
   //xsValue_MCnlo->SetFillColor(kWhite);
  xsValue_MCnlo->SetFillStyle(1001);
  xsValue_MCnlo->SetLineColor(kAzure);
  xsValue_MCnlo->SetLineWidth(1);
  xsValue_MCnlo->SetMarkerColor(kAzure);
  xsValue_MCnlo->SetMarkerSize(_msize);
  //  xsValue_MCnlo->SetMarkerStyle(21);
  xsValue_MCnlo->SetMarkerStyle(24);

  //  TCanvas * c1=new TCanvas("c1", "c1");
  
//xsValue_MCnlo->Draw("pey0");
   

  // Set the canvas and pads
  //----------------------------------------------------------------------------
  TCanvas* canvas = new TCanvas("wwxs", "wwxs", 600, 600);
  //  TCanvas* canvas = new TCanvas("wwxs", "wwxs"); 

  //defalut
  //TCanvas* canvas = new TCanvas("wwxs", "wwxs", 600, 850);
  //  TPad* pad1 = new TPad("pad1", "pad1", 0, 0.55, 1, 1.000);
  //TPad* pad2 = new TPad("pad2", "pad2", 0, 0.39, 1, 0.552);
  //TPad* pad3 = new TPad("pad3", "pad3", 0, 0.23, 1, 0.392);

  TPad* pad1 = new TPad("pad1", "pad1", 0, 0.49, 1, 1.000);
  TPad* pad2 = new TPad("pad2", "pad2", 0, 0.33, 1, 0.492);
  TPad* pad3 = new TPad("pad3", "pad3", 0, 0, 1, 0.332);
  
  pad1->SetTopMargin(0.09);
  pad2->SetTopMargin(0);
  pad3->SetTopMargin(0);

  pad1->SetBottomMargin(0);
  pad2->SetBottomMargin(0);
  //  pad3->SetBottomMargin(0.15);
  pad3->SetBottomMargin(0.45);

  pad1->SetLeftMargin(0.16);
  pad2->SetLeftMargin(0.16);
  pad3->SetLeftMargin(0.16);

  pad1->SetRightMargin(0.06);
  pad2->SetRightMargin(0.06);
  pad3->SetRightMargin(0.06);



  // pad1
  //----------------------------------------------------------------------------
  pad1->Draw();
  pad1->cd();
  pad1->SetLogy();


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

//.........这里部分代码省略.........
	pad01 = new TPad("pad01","pad01",0.0000,0.8333,0.2000,1.0000,0,0,0);
	pad02 = new TPad("pad02","pad02",0.2000,0.8333,0.4000,1.0000,0,0,0);
	pad03 = new TPad("pad03","pad03",0.4000,0.8333,0.6000,1.0000,0,0,0);
	pad04 = new TPad("pad04","pad04",0.6000,0.8333,0.8000,1.0000,0,0,0);
	pad05 = new TPad("pad05","pad05",0.0000,0.6666,0.2000,0.8333,0,0,0);
	pad06 = new TPad("pad06","pad06",0.2000,0.6666,0.4000,0.8333,0,0,0);
	pad07 = new TPad("pad07","pad07",0.4000,0.6666,0.6000,0.8333,0,0,0);
	pad08 = new TPad("pad08","pad08",0.6000,0.6666,0.8000,0.8333,0,0,0);
	pad09 = new TPad("pad09","pad09",0.0000,0.3333,0.3333,0.6666,0,0,0);
	pad10 = new TPad("pad10","pad10",0.3333,0.3333,0.6666,0.6666,0,0,0);
	pad11 = new TPad("pad11","pad11",0.6666,0.3333,1.0000,0.6666,0,0,0);
	pad12 = new TPad("pad12","pad12",0.0000,0.0000,0.3333,0.3333,0,0,0);
	pad13 = new TPad("pad13","pad13",0.3333,0.0000,0.6666,0.3333,0,0,0);
	pad14 = new TPad("pad14","pad14",0.6666,0.0000,1.0000,0.3333,0,0,0);
	pad15 = new TPad("pad15","pad15",0.8000,0.8333,1.0000,1.0000,0,0,0);
	pad16 = new TPad("pad16","pad16",0.8000,0.6666,1.0000,0.8333,0,0,0);
	pad01->Draw();pad02->Draw();pad03->Draw();pad04->Draw();pad05->Draw();pad06->Draw();pad07->Draw();pad08->Draw();pad09->Draw();pad10->Draw();pad11->Draw();pad12->Draw();pad13->Draw();pad14->Draw();pad15->Draw();pad16->Draw();

	// Everything below here makes graphs for each section of the canvas

	pad01->cd();
	TString titledpHe = "dp Cut from 3He Runs";
	cout << "Drawing " << titledpHe << "..." << endl;
	TH1F *HedpNoCut = new TH1F("HedpNoCut",titledpHe,400,-0.08,0.08);
	TH1F *HedpCut = new TH1F("HedpCut",titledpHe,400,-0.08,0.08);
	chainHeR->Draw("ExTgtCor_R.dp>>HedpNoCutR(400,-0.08,0.08)", "", "");
	chainHeR->Draw("ExTgtCor_R.dp>>HedpCutR(400,-0.08,0.08)", goldR, "");
	chainHeL->Draw("ExTgtCor_L.dp>>HedpNoCutL(400,-0.08,0.08)", "", "");
	chainHeL->Draw("ExTgtCor_L.dp>>HedpCutL(400,-0.08,0.08)", goldL, "");
	HedpNoCut->Add(HedpNoCutR,HedpNoCutL);
	HedpNoCut->SetTitle(titledpHe);
	HedpNoCut->Draw();
	HedpCut->Add(HedpCutR,HedpCutL);
	HedpCut->SetLineColor(kBlack);
	HedpCut->SetFillColor(kViolet);
	HedpCut->Draw("same");

	pad05->cd();
	TString titledpC = "dp Cut from 12C Runs";
	cout << "Drawing " << titledpC << "..." << endl;
	TH1F *CdpNoCut = new TH1F("CdpNoCut",titledpC,400,-0.8,0.8);
	TH1F *CdpCut = new TH1F("CdpCut",titledpC,400,-0.8,0.8);
	chainCR->Draw("ExTgtCor_R.dp>>CdpNoCutR(400,-0.08,0.08)", "", "");
	chainCR->Draw("ExTgtCor_R.dp>>CdpCutR(400,-0.08,0.08)", goldR, "");
	chainCL->Draw("ExTgtCor_L.dp>>CdpNoCutL(400,-0.08,0.08)", "", "");
	chainCL->Draw("ExTgtCor_L.dp>>CdpCutL(400,-0.08,0.08)", goldL, "");
	CdpNoCut->Add(CdpNoCutR,CdpNoCutL);
	CdpNoCut->SetTitle(titledpC);
	CdpNoCut->Draw();
	CdpCut->Add(CdpCutR,CdpCutL);
	CdpCut->SetLineColor(kBlack);
	CdpCut->SetFillColor(kViolet);
	CdpCut->Draw("same");

	pad02->cd();
	TString titleTargetHe = "Target Cut for 3He Runs";
	cout << "Drawing " << titleTargetHe << "..." << endl;
	TH1F *HeReactZNoCut = new TH1F("HeReactZNoCut",titleTargetHe,400,-0.3,0.3);
	TH1F *HeReactZCut = new TH1F("HeReactZCut",titleTargetHe,400,-0.3,0.3);
	chainHeR->Draw("ReactPt_R.z>>HeReactZNoCutR(400,-0.3,0.3)", "", "");
	chainHeR->Draw("ReactPt_R.z>>HeReactZCutR(400,-0.3,0.3)", targetR, "");
	chainHeL->Draw("ReactPt_L.z>>HeReactZNoCutL(400,-0.3,0.3)", "", "");
	chainHeL->Draw("ReactPt_L.z>>HeReactZCutL(400,-0.3,0.3)", targetL, "");
	HeReactZNoCut->Add(HeReactZNoCutR,HeReactZNoCutL);
	HeReactZNoCut->SetTitle(titleTargetHe);
	HeReactZNoCut->Draw();

//.........这里部分代码省略.........
    // pp (average Run-2/3 parameterization)
    // weight = 2.0*pi*pt * (546.2/42.2) / pow(exp(-0.2664*pt) + pt/0.7703,8.718);
    //    weight = weight*(1+0.6362/exp(7.736*pt));

    ptpion->AddBinContent(ibin,weight);
    intpion = intpion + weight*binwidth;

    // if ( pt<1.5 ) weight = 2.0 * pi * pt * otherhagedorn;
    // if ( pt>=1.5 && pt<=5.5 ) 
    //   weight = 2.0 * pi * pt * (trans*otherhagedorn + (1.-trans)*otherpower);
    // if ( pt>5.5 ) weight = 2.0*pi*pt*otherpower;

    weight = 2.0*pi*pt * (450.8/42.2) / 
      pow(exp(-0.5553*sqrt(mtother*mtother-pimass*pimass)) + 
	  sqrt(mtother*mtother-pimass*pimass)/0.6467,8.046);

    ptothermt->AddBinContent(ibin,weight);
    intothermt = intothermt + weight*binwidth;

    weight = norm*InitializePtHydro(omass,t_fo,beta,pt);
    ptother->AddBinContent(ibin,weight);
    intother = intother + weight*binwidth;
  }

  cout << "Integral pion    (dN/dy): " << intpion << endl;
  cout << "Integral other   (dN/dy): " << intother << endl;
  cout << "Integral othermt (dN/dy): " << intothermt << endl;

  norm = intpion/intothermt;
  ptothermt->Scale(norm);
 
  for (int ibin=1; ibin<=nbins; ibin++ )
  {
    weight = ptothermt->GetBinContent(ibin)/ptpion->GetBinContent(ibin);
    ratiomt->AddBinContent(ibin,weight);
  }

  norm = expected/ratiomt->GetMaximum();
  ratiomt->Scale(norm);
  ptothermt->Scale(norm);
  cout << "meson/pion ratio should be set to: " << norm << endl;

  norm = expected*intothermt/intother;
  ptother->Scale(norm);
 
  for (int ibin=1; ibin<=nbins; ibin++ )
  {
    weight = ptother->GetBinContent(ibin)/ptpion->GetBinContent(ibin);
    ratio->AddBinContent(ibin,weight);
  }

  ptpion->SetXTitle("p_{t} [GeV/c]");
  ptpion->SetYTitle("(1/2#pi p_{t})dN/dp_{t}dy [(c/GeV)^{2}]");
  ptpion->SetLineWidth(5);
  ptpion->SetLineColor(1);
  ptother->SetLineWidth(5);
  ptother->SetLineColor(2);
  ptothermt->SetLineWidth(5);
  ptothermt->SetLineColor(4);
  ptpion->Draw("");
  ptothermt->Draw("same");
  ptother->Draw("same");

  TLatex *text_pion0 = new TLatex(3.2,1.5,"#pi^{0}");
  text_pion0->SetTextColor(1);
  text_pion0->Draw();
  TLatex *text_pionm = new TLatex(3.2,1.1,"#eta m_{t} scaled");
  text_pionm->SetTextColor(4);
  text_pionm->Draw();
  TLatex *text_pionp = new TLatex(3.2,0.7,"#eta hydro");
  text_pionp->SetTextColor(2);
  text_pionp->Draw();

   c2 = new TCanvas("c2","A Simple Graph with error bars",200,10,700,500);

   c2->SetFillColor(10);
   c2->SetGrid();
   c2->SetLogy(0);
   c2->SetTicks();
 
   c2->cd();

  ratiomt->SetXTitle("p_{t} [GeV/c]");
  ratiomt->SetYTitle("#eta / #pi^{0}");
  ratio->SetLineWidth(5);
  ratio->SetLineColor(2);
  ratiomt->SetLineWidth(5);
  ratiomt->SetLineColor(4);
  ratiomt->Draw();
  ratio->Draw("same");

  TLatex *text_pionm = new TLatex(3.2,0.42,"#eta m_{t} scaled");
  text_pionm->SetTextColor(4);
  text_pionm->Draw();
  TLatex *text_pionp = new TLatex(3.2,0.35,"#eta hydro");
  text_pionp->SetTextColor(2);
  text_pionp->Draw();

  return;
}

void BkgDemo_Diplot(TCanvas* c, TH1F** h, TLegend* l1, TLegend* l2){
		//Make a split canvas, comparison of methods and tag above, bkg subtracted tag below
//	printf("\nDo BkgDemo_Diplot for canvas %s\n",c->GetName());
//	printf("integrals stage1: tag %f usb %f lsb %f\n",h[1]->Integral(),h[5]->Integral(),h[6]->Integral());

	
	TLatex * TEX_CMSPrelim;
	if(preliminary) TEX_CMSPrelim = new TLatex(0.177136,0.953368,"CMS Preliminary"); 
	else TEX_CMSPrelim = new TLatex(0.177136,0.953368,"CMS"); 
	PrettyLatex(TEX_CMSPrelim,0.03);
	TLatex * TEX_E_TeV = new TLatex(0.800251,0.953368,"#sqrt{s} = 8 TeV");
	PrettyLatex(TEX_E_TeV,0.03);
	TLatex * TEX_lumi_fb = new TLatex(0.621859,0.953368,Form("#intL dt = %.1f fb^{-1}",Integrated_Luminosity_Data));
	PrettyLatex(TEX_lumi_fb,0.03);
	string sp1 = string(c->GetName())+"_p1";
	string sp2 = string(c->GetName())+"_p2";
    TPad *p1 = new TPad((char*)sp1.c_str(),"",0.,0.3,1.,1.);
    p1->SetBottomMargin(0);
    p1->cd();
//	c->Divide(1,2);
//	c->cd(1);
	PrettyHist(h[3],kRed);
	PrettyHist(h[5],kBlue);
	PrettyHist(h[1]);
	PrettyHist(h[6],kGreen);
	//PrettyHist(h[3],kRed);
	PrettyHist(h[7],kTeal);
PrettyBlock2(h[5],kBlue,3354,2);
	//PrettyBlock(h[5],kBlue,string("//thatch"));//PrettyMarker(h[5],kBlue,4);
	PrettyMarker(h[1]);
PrettyBlock2(h[6],kGreen,3345,2);
	//PrettyBlock(h[6],kGreen,string("\\thatch"));//PrettyMarker(h[6],kGreen,4);
	//PrettyMarker(h[3],kRed);
	PrettyMarker(h[7],kTeal);
//	h[5]->Scale(mHwidth/sidebandwidth);//lsb scaled
//	h[6]->Scale(mHwidth/sidebandwidth);//usb scaled
//	printf("integrals stage2: tag %f usb %f lsb %f\n",h[1]->Integral(),h[5]->Integral(),h[6]->Integral());
//	printf("ranges maximi before setrange tag %f lsb %f, usb %f\n",h[1]->GetMaximum(),h[5]->GetMaximum(),h[6]->GetMaximum());
//	printf("integrals stage3: tag %f usb %f lsb %f\n",h[1]->Integral(),h[5]->Integral(),h[6]->Integral());
//	printf("ranges maximi before after tag %f lsb %f usb %f\n",h[1]->GetMaximum(),h[5]->GetMaximum(),h[6]->GetMaximum());
	h[3]->SetFillStyle(0);//open rectangle
	h[3]->SetLineColor(kRed);
	h[3]->SetLineWidth(4);

	h[6]->SetMinimum(0.0);
	//float linmax = h[6]->GetMaximum();
	//float linmin = h[6]->GetMinimum();
	playNiceWithLegend(h[6],0.30,0.0);
	playNiceWithLegend(h[3],0.30,0.0);
	playNiceWithLegend(h[5],0.30,0.0);
	playNiceWithLegend(h[1],0.30,0.0);
	SameRange(h[3],h[1]);
	SameRange(h[3],h[5],h[6]);
	SameRange(h[3],h[1]);
	h[3]->Draw("e2p");
	h[5]->Draw("e2psame");
	h[6]->Draw("e2psame");
	h[3]->Draw("e2psame");
	//if(showTag) h[1]->Draw("e1psame");//tag
	TPad *p2 = new TPad((char*)sp2.c_str(),"",0.,0.125,1.,0.3);
    p2->SetTopMargin(0);
    p2->cd();
//	c->cd(2);
	TAxis * x = h[7]->GetXaxis();
	TAxis * y = h[7]->GetYaxis();
	float fontsize = 0.25;
	float fontsizeY = 0.10;
	x->SetTitleSize(fontsize);
	y->SetTitleSize(fontsizeY);
	x->SetLabelSize(fontsize);
	y->SetLabelSize(fontsizeY);

	h[7]->GetYaxis()->SetRangeUser(0.,2.);
	h[7]->Draw("e1p");
	TLine *OneLine = new TLine(x->GetXmin(),1.0,x->GetXmax(),1.0);
	OneLine->SetLineColor(kBlack);
	OneLine->SetLineWidth(2);
	OneLine->SetLineStyle(7);//dashed.
	OneLine->Draw("same");
	h[7]->Draw("e1psame");
	
	p1->cd();
//  c->cd(1);
	if(showTag) l1->AddEntry(h[1],"Higgs Mass Region");
        TH1F * box = new TH1F("box","asdf",1,0,1);
        box->SetMarkerColor(kRed);
        box->SetMarkerStyle(25);
        box->SetMarkerSize(2);
        box->SetLineColor(0);
	l1->AddEntry(box,"Data Driven Background");
	l1->AddEntry(h[5],"Lower Mass Sideband");
	l1->AddEntry(h[6], "Upper Mass Sideband");
	l1->Draw("same");
        TEX_CMSPrelim->Draw("same");
        TEX_E_TeV->Draw("same");
        TEX_lumi_fb->Draw("same");
	p2->cd();
		//c->cd(2);
	l2->SetTextSize(fontsize*0.8);
	l2->AddEntry(h[7],"Lower/Upper Sideband Ratio");
//.........这里部分代码省略.........

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

		/*Debug*/ if(printlevel > 0) cout << "Build mass plots" << endl;
	
		//diphoton mass plots, one for each data plot, for each topology. 
	if(makediphoMassPlot){
	for (int iTop = 0; iTop<nEventTopologies; iTop++) {
		/*Debug*/ if(printlevel > 1) cout << "first loop, iTop = "<<iTop <<endl;
		string canvName = string("c_mgg")+s_EventTopology_v2[iTop]+"_unsliced";
		/*Debug*/ if(printlevel > 1) cout << "made tempCanv with name "<<canvName  <<endl;
		TCanvas* tempCanv = newTCanvas((char*)canvName.c_str()); //char*error
		/*Debug*/ if(printlevel > 1) cout << "Enter the merky deptsh of diphoMassPlot"<<endl;
		diphoMassPlot(h_mGG_unsliced[Data][s_EventTopology[iTop]],tempCanv,mgg_fit_curve[iTop]);
		/*Debug*/ if(printlevel > 1) cout << "fin dinphoMass plot, try to write to fplots" <<endl;
		tempCanv->Write();
		if(saveImages) SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggifpdfjpgepspng");//save as all types
		//if(saveImages) SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggifpdfeps");//save as all types
	}
	}
	

		// %%%%%%%%%%%%%%%%%%%%%%%% Kinematics and backgrounds %%%%%%%%%%%%%%%%%
//     ___  __        ___                   ___  __     __    
//    / _ )/ /_____ _/ _ \___ __ _  ___    / _ \/ /__  / /____
//   / _  /  '_/ _ `/ // / -_)  ' \/ _ \  / ___/ / _ \/ __(_-<
//  /____/_/\_\\_, /____/\__/_/_/_/\___/ /_/  /_/\___/\__/___/
//            /___/                                           

	/*Debug*/ if(printlevel > 0) cout << "Build Kin Var plots type 1" << endl;
		//Arrangement 1
        TH1F * box = new TH1F("box","asdf",1,0,1);
        box->SetMarkerColor(kRed);
        box->SetMarkerStyle(25);
        box->SetMarkerSize(2);
        box->SetLineColor(0);
		///Demonstrate background. Plot two side bands with their ratio as a sub-plot. include a switchable Tag.


	if(makeBkgDemo){
	for (int iTop = 0; iTop<nEventTopologies; iTop++){
		for (int kKinVar = 0; kKinVar<nKinemVars; kKinVar++){
			/*Debug*/ if(printlevel >2) cout << "iTop="<<iTop<<" kKinVar="<<kKinVar<<endl;
			if (s_EventTopology[iTop].compare("metCut")==0 && s_KinemVars[kKinVar].compare("MET")==0 ){
				/*Debug*/ if(printlevel >2) cout << "skipping met plot for met cut"<<endl;
				continue;
			}
			string canvName = string("BkgFull_")+s_KinemVars[kKinVar]+s_EventTopology_v2[iTop];
			/*Debug*/ if(printlevel >6) {TH1F** h = KinVarHistMap[Data][s_EventTopology[iTop]][s_KinemVars[kKinVar]];
			/*Debug*/ if(printlevel >6) printf("Scan1 Ratio bin1 for %s %s  = %f\n",s_EventTopology[iTop].c_str(),s_KinemVars[kKinVar].c_str(),h[7]->GetBinContent(1));}
			TCanvas* tempCanv = newTCanvas((char*)canvName.c_str());
			BkgDemo_Diplot(tempCanv,
					 KinVarHistMap[Data][s_EventTopology[iTop]][s_KinemVars[kKinVar]],
						   makeL1(0.2072864,0.7203704,0.4070352,0.9203704),
						   makeL2(0.51005,0.706436,0.708543,0.902669));
			//tempCanv->Write();
			//if(saveImages)SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggif");//save as all types
			//if(saveImages)SaveCanvas(tempCanv,plotsdir+tempCanv->GetName(),"ggifpdfeps");//save as all types
		}//edn for each kinematic varriable
	}//end for each topology
	}

//     ___  __       ___  ___  __     __    
//    / _ )/ /_____ <  / / _ \/ /__  / /____
//   / _  /  '_/ _ `/ / / ___/ / _ \/ __(_-<
//  /____/_/\_\\_, /_/ /_/  /_/\___/\__/___/
//            /___/                         

//.........这里部分代码省略.........
    for (int j=0; j<nBi2East; j++) {
      meanBi2 += resBi2East[j];
    }
    for (int j=0; j<nBi2West; j++) {
      meanBi2 += resBi2West[j];
    }
    meanBi2 = meanBi2 / (double) (nBi2East+nBi2West);
    
    double sigmaBi2 = 0.;
    for (int j=0; j<nBi2East; j++) {
      sigmaBi2 += (resBi2East[j] - meanBi2)*(resBi2East[j] - meanBi2);
    }
    for (int j=0; j<nBi2West; j++) {
      sigmaBi2 += (resBi2West[j] - meanBi2)*(resBi2West[j] - meanBi2);
    }
    sigmaBi2 = sigmaBi2 / (double) (nBi2East+nBi2West);
    sigmaBi2 = sqrt( sigmaBi2 );
    
    cout << "meanBi2 = " << meanBi2 << endl;
    cout << "     sigma = " << sigmaBi2 << endl;
    errEnv << "meanBi2 = " << meanBi2 << endl;
    errEnv << "sigma = " << sigmaBi2 << endl;
  }

  cout << endl << "Results of Gaussian Fits:\n";
  errEnv << endl << "Results of Gaussian Fits:\n";

  // Ce East
  c1 = new TCanvas("c1", "c1");
  c1->SetLogy(0);

  hisCeEast->SetXTitle("East Ce E_{Q} Error [keV]");
  hisCeEast->GetXaxis()->CenterTitle();
  hisCeEast->SetLineColor(1);
  hisCeEast->Draw();
  hisCeEast->Fit("gaus", "", "", -8.0, 8.0);
  cout << "meanCeEast = " << gaus->GetParameter(1) << endl;
  cout << "     sigma = " << gaus->GetParameter(2)  << endl;
  errEnv << "meanCeEast = " << gaus->GetParameter(1) << endl;
  errEnv << "sigma = " << gaus->GetParameter(2) << endl;

  //In East
  c2 = new TCanvas("c2", "c2");
  c2->SetLogy(0);

  hisInEast->SetXTitle("East In E_{Q} Error [keV]");
  hisInEast->GetXaxis()->CenterTitle();
  hisInEast->SetLineColor(1);
  hisInEast->Draw();
  hisInEast->Fit("gaus", "", "", -8.0, 8.0);
  cout << "meanInEast = " << gaus->GetParameter(1) << endl;
  cout << "     sigma = " << gaus->GetParameter(2)  << endl;
  errEnv << "meanInEast = " << gaus->GetParameter(1) << endl;
  errEnv << "sigma = " << gaus->GetParameter(2) << endl;

  // Sn East
  c3 = new TCanvas("c3", "c3");
  c3->SetLogy(0);

  hisSnEast->SetXTitle("East Sn E_{Q} Error [keV]");
  hisSnEast->GetXaxis()->CenterTitle();
  hisSnEast->SetLineColor(1);
  hisSnEast->Draw();
  hisSnEast->Fit("gaus", "", "", -25.0, 25.0);
  cout << "meanSnEast = " << gaus->GetParameter(1) << endl;
  cout << "     sigma = " << gaus->GetParameter(2)  << endl;