C++ RooDataSet::numEntries方法代码示例

void trig_pt()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  TFile *f_MC= new TFile("TnP_WptCutToTrig_MCptminus.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/Tnp_WptCut_to_Mu15_eta2p1_pt/cnt_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &ptMC=pointMC["pt"],&effMC = pointMC["efficiency"];
    XMC[i]=ptMC.getVal();
    XMCerrL[i]=-ptMC.getAsymErrorLo();
    XMCerrH[i]=ptMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(11,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->Draw("AP");
  //grMC->Draw("psame");
  myCan->SaveAs("trig_McMinus_pt.png");
  myCan->SaveAs("trig_McMinus_pt.eps");

}

void glbToId_eta()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  TFile *f_MC= new TFile("TnP_GlbToID_MCetaplus_WptTight2012_eta.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/WptTight2012_eta/fit_eff");
  //RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/WptTight2012_eta/cnt_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &etaMC=pointMC["eta"],&effMC = pointMC["efficiency"];
    XMC[i]=etaMC.getVal();
    XMCerrL[i]=-etaMC.getAsymErrorLo();
    XMCerrH[i]=etaMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(Nbin,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->Draw("AP");
  //grMC->Draw("psame");
  myCan->SaveAs("glbToId_MCplus_eta.png");
  myCan->SaveAs("glbToId_MCplus_eta.eps");

}

void sc_wptCut_P_et()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  /************************
  TFile *f_RD= new TFile("TnP_Z_Trigger_RDpt.root","read");
  RooDataSet *dataset = (RooDataSet*)f_RD->Get("tpTree/Track_To_TightCombRelIso_Mu15_eta2p1_pt/fit_eff");
  cout<<"ntry: "<<dataset->numEntries()<<endl;
  double X[11],XerrL[11],XerrH[11],Y[11],YerrLo[11],YerrHi[11];
  for(int i(0); i<dataset->numEntries();i++)
  {
    const RooArgSet &point=*dataset->get(i);
    RooRealVar &pt=point["pt"],&eff = point["efficiency"];
    X[i]=pt.getVal();
    XerrL[i]=-pt.getAsymErrorLo();
    XerrH[i]=pt.getAsymErrorHi();
    Y[i]=eff.getVal();
    YerrLo[i]=-eff.getAsymErrorLo();
    YerrHi[i]=eff.getAsymErrorHi();
  }
  gr=new TGraphAsymmErrors(11,X,Y,XerrL,XerrH,YerrLo,YerrHi);
  gr->Draw("AP");
***************************/
  TFile *f_MC= new TFile("efficiency-mc-SCToPfElectron_et_P.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("SuperClusterToPFElectron/SCtoWptCut_efficiency/cnt_eff");
  //RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/Track_with_TightCombRelIso_to_Mu15_eta2p1_pt/fit_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;

  double XMC[binSize],XMCerrL[binSize],XMCerrH[binSize],YMC[binSize],YMCerrLo[binSize],YMCerrHi[binSize];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &ptMC=pointMC["probe_sc_et"],&effMC = pointMC["efficiency"];
    XMC[i]=ptMC.getVal();
    XMCerrL[i]=-ptMC.getAsymErrorLo();
    XMCerrH[i]=ptMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(binSize,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  //myCan->SetLogx();
  grMC->Draw("AP");
  //grMC->Draw("psame");
  TLine *myLine=new TLine(25,0,25,1);
  myLine->Draw("same");
  myCan->SaveAs("sc_wptCut_P_et.png");
  myCan->SaveAs("sc_wptCut_P_et.eps");

}

RooDataHist * genHistFromModelPdf(const char * name, RooAbsPdf *model,  RooRealVar *var,    double ScaleLumi,  int range, int rebin, int seed ) {
  double genEvents =  model->expectedEvents(*var);
  TRandom3 *rndm = new TRandom3();
  rndm->SetSeed(seed);
  double nEvt = rndm->PoissonD( genEvents) ;
  int intEvt = ( (nEvt- (int)nEvt) >= 0.5) ? (int)nEvt +1 : int(nEvt);
  RooDataSet * data = model->generate(*var ,   intEvt   );
  cout<< " expected events for " << name << " = "<< genEvents << endl; 
  cout<< " data->numEntries() for name " << name << " == " << data->numEntries()<< endl;
  // cout<< " nEvt from PoissonD for" << name << " == " << nEvt<< endl;
  //cout<< " cast of nEvt  for" << name << " == " << intEvt<< endl; 
  RooAbsData *binned_data = data->binnedClone();
  TH1 * toy_hist = binned_data->createHistogram( name, *var, Binning(range/rebin )  );
  for(int i = 1; i <= toy_hist->GetNbinsX(); ++i) {
    toy_hist->SetBinError( i,  sqrt( toy_hist->GetBinContent(i)) );
    if(toy_hist->GetBinContent(i) == 0.00) {
      cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;   
    }
    if(toy_hist->GetBinContent(i) < 0.1) {
      toy_hist->SetBinContent(i, 0.0);
      toy_hist->SetBinError(i, 0.0);
      cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;   
    }  
  }
  RooDataHist * toy_rooHist = new RooDataHist(name, name , RooArgList(*var), toy_hist );
  return toy_rooHist; 
}

pair<double,double> datEvents(RooWorkspace *work, int m_hyp, int cat, bool spin=false){
  
  vector<double> result;
  RooDataSet *data = (RooDataSet*)work->data(Form("data_mass_cat%d",cat));
  double evs = data->numEntries();
  double evsPerGev;
  if (!spin) evsPerGev = data->sumEntries(Form("CMS_hgg_mass>=%4.1f && CMS_hgg_mass<%4.1f",double(m_hyp)-0.5,double(m_hyp)+0.5));
  else evsPerGev = data->sumEntries(Form("mass>=%4.1f && mass<%4.1f",double(m_hyp)-0.5,double(m_hyp)+0.5));
  return pair<double,double>(evs,evsPerGev);
}

void wspaceread_backgrounds(int channel = 1)
{

	gSystem->AddIncludePath("-I$ROOFITSYS/include");
	gROOT->ProcessLine(".L ~/tdrstyle.C");
	
	string schannel;
	if (channel == 1) schannel = "4mu";
	if (channel == 2) schannel = "4e";
	if (channel == 3) schannel = "2mu2e";
	std::cout << "schannel = " << schannel << std::endl;

	// R e a d   w o r k s p a c e   f r o m   f i l e
	// -----------------------------------------------
	
	// Open input file with workspace (generated by rf14_wspacewrite)
	char infile[192];
	sprintf(infile,"/scratch/hep/ntran/dataFiles/HZZ4L/datasets/datasets/%s/ZZAnalysisTree_ZZTo4L_lowmass.root",schannel.c_str());
	TFile *f = new TFile(infile) ;
	char outfile[192];
	sprintf( outfile, "figs/pdf_%s_bkg_highmass.eps", schannel.c_str() );
	//f->ls();
	
	
	RooDataSet* set = (RooDataSet*) f->Get("data");
	RooArgSet* obs = set->get() ;
	obs->Print();
	RooRealVar* CMS_zz4l_mass = (RooRealVar*) obs->find("CMS_zz4l_mass") ;
	
	for (int i=0 ; i<set->numEntries() ; i++) { 
		set->get(i) ; 
		//cout << CMS_zz4l_mass->getVal() << " = " << set->weight() << endl ; 
	} 
	
	gSystem->Load("PDFs/RooqqZZPdf_cxx.so");
	//gSystem->Load("PDFs/RooggZZPdf_cxx.so");

	// LO contribution

	//RooRealVar m4l("m4l","m4l",100.,1000.);
	RooRealVar a1("a1","a1",224.,100.,1000.);
	RooRealVar a2("a2","a2",-209.,-1000.,1000.);
	RooRealVar a3("a3","a3",121.,20.,1000.);
	RooRealVar a4("a4","a4",-0.022,-10.,10.);
	RooRealVar b1("b1","b1",181.,100.,1000.);
	RooRealVar b2("b2","b2",707.,0.,1000.);
	RooRealVar b3("b3","b3",60.,20.,1000.);
	RooRealVar b4("b4","b4",0.04,-10.,10.);
	RooRealVar b5("b5","b5",5.,0.,1000.);
	RooRealVar b6("b6","b6",0.,-10.,10.);
	RooRealVar frac_bkg("frac_bkg","frac_bkg",0.5,0.,1.);
	
	//a1.setConstant(kTRUE);
	//a2.setConstant(kTRUE);
	//a3.setConstant(kTRUE);
	//a4.setConstant(kTRUE);
	//b1.setConstant(kTRUE);
	//b2.setConstant(kTRUE);
	//b3.setConstant(kTRUE);
	//b4.setConstant(kTRUE);
	//b5.setConstant(kTRUE);
	//b6.setConstant(kTRUE);

	RooqqZZPdf bkg_qqzz("bkg_qqzz","bkg_qqzz",*CMS_zz4l_mass,a1,a2,a3,a4,b1,b2,b3,b4,b5,b6,frac_bkg);

	RooFitResult *r = bkg_qqzz.fitTo( *set, SumW2Error(kTRUE) );//, Save(kTRUE), SumW2Error(kTRUE)) ;

	// Plot Y
	RooPlot* frameM4l = CMS_zz4l_mass->frame(Title("M4L"),Bins(100)) ;
	set->plotOn(frameM4l) ;
	bkg_qqzz.plotOn(frameM4l) ;
	
	TCanvas *c = new TCanvas("c","c",800,600);
	c->cd();
	frameM4l->Draw();
	
	
	
	/*
	// Retrieve workspace from file
	RooWorkspace* w = (RooWorkspace*) f->Get("workspace") ;
	
	w->Print();
	
	///*
	RooRealVar* CMS_zz4l_mass = w->var("CMS_zz4l_mass") ;
	RooAbsPdf* background_nonorm = w->pdf("background_nonorm") ;
	//RooAbsData* backgroundData = w->data("backgroundData") ;
	RooAbsData* data_bkg_red = w->data("data_bkg_red") ;
	
	RooArgSet* obs = data_bkg_red->get() ; 
	RooRealVar* xdata = obs->find(CMS_zz4l_mass.GetName()) ; 
	for (int i=0 ; i<data_bkg_red->numEntries() ; i++) { 
		data_bkg_red->get(i) ; 
		cout << xdata->getVal() << " = " << data_bkg_red->weight() << endl ; 
	} 
	std::cout << "nEntries = " << data_bkg_red->numEntries() << std::endl;
	obs->Print();

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

//.........这里部分代码省略.........
           && fabs( zeeTree->fEle2Eta) < 2.5 )) continue;

    //*************************************************************************
    //pt bins and eta bins
    //*************************************************************************
    Int_t Ele1PtBin = -1;
    Int_t Ele1EtaBin = -1;
    Int_t Ele2PtBin = -1;
    Int_t Ele2EtaBin = -1;
    if (ele1pt > 10 && ele1pt < 20) Ele1PtBin = 0;
    else if (ele1pt < 30) Ele1PtBin = 1;
    else if (ele1pt < 40) Ele1PtBin = 2;
    else Ele1PtBin = 3;
    if (ele2pt > 10 && ele2pt < 20) Ele2PtBin = 0;
    else if (ele2pt < 30) Ele2PtBin = 1;
    else if (ele2pt < 40) Ele2PtBin = 2;
    else Ele2PtBin = 3;
    if (fabs(zeeTree->fEle1SCEta) < 1.0) Ele1EtaBin = 0;
    else if (fabs(zeeTree->fEle1SCEta) < 1.479) Ele1EtaBin = 1;
    else Ele1EtaBin = 2;
    if (fabs(zeeTree->fEle2SCEta) < 1.0) Ele2EtaBin = 0;
    else if (fabs(zeeTree->fEle2SCEta) < 1.479) Ele2EtaBin = 1;
    else Ele2EtaBin = 2;

    if (CategoryBin == 0) { 
      if (!(Ele1EtaBin == 0 && Ele2EtaBin == 0)) continue; 
    }
    else if (CategoryBin == 1) {
      if (!(Ele1EtaBin == 1 && Ele2EtaBin == 1)) continue; 
    }
    else if (CategoryBin == 2) {
      if (!(Ele1EtaBin == 2 && Ele2EtaBin == 2)) continue;
    }
    
    //*************************************************************************
    // restrict range of mass
    //*************************************************************************
    double zMass = (ele1FourVector+ele2FourVector).M();
    if (zMass < minMass || zMass > maxMass) continue;

    //*************************************************************************
    //set mass variable
    //*************************************************************************
    zMassArgSet.setRealValue("mass", zMass);

    data->add(zMassArgSet);
  }

  cout << "dataset size: " << data->numEntries() << endl;


  RooRealVar *cbBias = (RooRealVar*)w->var("#Deltam_{CB}");
  RooRealVar *cbSigma = (RooRealVar*)w->var("sigma_{CB}");
  RooRealVar *cbCut   = (RooRealVar*)w->var("a_{CB}");
  RooRealVar *cbPower = (RooRealVar*)w->var("n_{CB}");

//   // Now if it's a restricted fit, fix values of cbCut and cbPower to MC values.
//   if (isRestricted) {
//     cbCut.setConstant(kTRUE);
//     cbPower.setConstant(kTRUE);
//   }

  // Mass model for signal electrons p.d.f.
  RooAddPdf *model = (RooAddPdf*)w->pdf("model");


  TCanvas* c = new TCanvas("c","c", 0,0,800,600);

  //========================== Plotting  ============================
  //Create a frame
  RooPlot* plot = mass->frame(Range(minMass,maxMass),Bins(nbins));
  // Add data and model to canvas
  plot->SetTitle("");
  plot->GetYaxis()->SetTitleOffset(1.4);
  plot->GetXaxis()->SetTitle("m_{ee} (GeV/c^{2})");

  data->plotOn(plot);
  model->plotOn(plot);
//   model->paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(*cbBias, *cbSigma, *cbCut, *cbPower)), Layout(0.60,0.90,0.90));
  model->paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(*cbBias, *cbSigma, *cbCut, *cbPower)), Layout(0.12,0.38,0.60));
  plot->getAttText()->SetTextSize(.025);
  plot->Draw();

  // Print Fit Values
  TLatex *tex = new TLatex();
  tex->SetNDC();
  tex->SetTextSize(.04);
  tex->SetTextFont(2);
  tex->DrawLatex(0.195,0.775, "Run 2012A/B");
  tex->Draw();
//   tex->SetTextSize(0.022);
//   tex->DrawLatex(0.195, 0.75, "Z #rightarrow ee^{+}");
//   tex->SetTextSize(0.024);
//   tex->DrawLatex(0.645, 0.59, Form("BW Mean = %.2f GeV/c^{2}", bwMean.getVal()));
//   tex->DrawLatex(0.645, 0.54, Form("BW #sigma = %.2f GeV/c^{2}", bwWidth.getVal()));
  c->Update();
  c->SaveAs((outFilename + ".gif").c_str());


}

//____________________________________
void DoSPlot(RooWorkspace* ws){
  std::cout << "Calculate sWeights" << std::endl;

  RooAbsPdf* model = ws->pdf("model");
  RooRealVar* nsig = ws->var("nsig");
  RooRealVar* nBbkg = ws->var("nBbkg");
  RooRealVar* nbkg = ws->var("nbkg");
  RooRealVar* nbkg2 = ws->var("nbkg2");
  RooDataSet* data = (RooDataSet*) ws->data("data");

  // fit the model to the data.
  model->fitTo(*data, Extended() );

  RooMsgService::instance().setSilentMode(true);

  // Now we use the SPlot class to add SWeights to our data set
  // based on our model and our yield variables
  RooStats::SPlot* sData = new RooStats::SPlot("sData","An SPlot",
					       *data, model, RooArgList(*nsig,*nBbkg,*nbkg,*nbkg2) );


  // Check that our weights have the desired properties

  std::cout << "Check SWeights:" << std::endl;

  std::cout << std::endl <<  "Yield of sig is " 
	    << nsig->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("nsig") << std::endl;

  std::cout << std::endl <<  "Yield of Bbkg is " 
	    << nBbkg->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("nBbkg") << std::endl;

  std::cout << std::endl <<  "Yield of bkg is " 
	    << nbkg->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("nbkg") << std::endl;

  std::cout << std::endl <<  "Yield of bkg2 is " 
	    << nbkg2->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("nbkg2") << std::endl;

  cout << endl;   cout << endl;   cout << endl;
  float sum20=0;
  float sum50=0;
  float sum100=0;
  float sum200=0;
  float sum300=0;
  float sum600=0;
  float sum900=0;
  float sum1200=0;
  float total=0;

  // saving weights into a file
  ofstream myfile;
  myfile.open ("weights.txt");
  // plot the weight event by event with the Sum of events values as cross-check
  for(Int_t i=0; i < data->numEntries(); i++) {
      //myfile << sData->GetSWeight(i,"nsig") << " " << sData->GetSWeight(i,"nBbkg") << " " << sData->GetSWeight(i,"nbkg") << " " << sData->GetSWeight(i,"nbkg2") << endl;  
      //myfile << sData->GetSWeight(i,"nsig") <<endl;
    myfile << (unsigned int) data->get(i)->getRealValue("run")
             << " " << (unsigned int) data->get(i)->getRealValue("event")
	   << " " << (float) data->get(i)->getRealValue("FourMu_Mass")
             << " " << sData->GetSWeight(i,"nsig")
             << endl;
     // std::cout << "nsig Weight   " << sData->GetSWeight(i,"nsig") 
     //		<< "   nBbkg Weight   " << sData->GetSWeight(i,"nBbkg")
     //		<< "   nbkg Weight   " << sData->GetSWeight(i,"nbkg")
     //		<< "   nbkg2 Weight  " << sData->GetSWeight(i,"nbkg2")
//		<< "   Total Weight   " << sData->GetSumOfEventSWeight(i) 
//		<< std::endl;
      total+=sData->GetSWeight(i,"nsig");         
      if(i<20) sum20+=sData->GetSWeight(i,"nsig");
      if(i<50) sum50+=sData->GetSWeight(i,"nsig");
      if(i<100) sum100+=sData->GetSWeight(i,"nsig");
      if(i<200) sum200+=sData->GetSWeight(i,"nsig");
      if(i<300) sum300+=sData->GetSWeight(i,"nsig");
      if(i<600) sum600+=sData->GetSWeight(i,"nsig");
      if(i<900) sum900+=sData->GetSWeight(i,"nsig");
      if(i<1200) sum1200+=sData->GetSWeight(i,"nsig");

    }
  myfile.close();

  std::cout << std::endl;

  std::cout<<"Sum of the sWeights is: "<<total<<std::endl;
  std::cout<<"Sum of the first 20 sWeights is: "<<sum20<<std::endl;
  std::cout<<"Sum of the first 50 sWeights is: "<<sum50<<std::endl;
  std::cout<<"Sum of the first 100 sWeights is: "<<sum100<<std::endl;
  std::cout<<"Sum of the first 200 sWeights is: "<<sum200<<std::endl;
  std::cout<<"Sum of the first 300 sWeights is: "<<sum300<<std::endl;
  std::cout<<"Sum of the first 600 sWeights is: "<<sum600<<std::endl;
  std::cout<<"Sum of the first 900 sWeights is: "<<sum900<<std::endl;
  std::cout<<"Sum of the first 1200 sWeights is: "<<sum1200<<std::endl;
  std::cout<<"Total # of events: "<<data->numEntries()<<std::endl;

  // import this new dataset with sWeights
  std::cout << "import new dataset with sWeights" << std::endl;
  ws->import(*data, Rename("dataWithSWeights"));
//.........这里部分代码省略.........

// The actual job
void backgroundFits_qqzz_1Dw(int channel, int sqrts, int VBFtag)
{
  if(sqrts==7)return;
  TString schannel;
  if      (channel == 1) schannel = "4mu";
  else if (channel == 2) schannel = "4e";
  else if (channel == 3) schannel = "2e2mu";
  else cout << "Not a valid channel: " << schannel << endl;

  TString ssqrts = (long) sqrts + TString("TeV");

  cout << "schannel = " << schannel << "  sqrts = " << sqrts << " VBFtag = " << VBFtag << endl;

  TString outfile;
  if(VBFtag<2) outfile = "CardFragments/qqzzBackgroundFit_" + ssqrts + "_" + schannel + "_" + Form("%d",int(VBFtag)) + ".txt";
  if(VBFtag==2) outfile = "CardFragments/qqzzBackgroundFit_" + ssqrts + "_" + schannel + ".txt";
  ofstream of(outfile,ios_base::out);
  of << "### background functions ###" << endl;


  gSystem->AddIncludePath("-I$ROOFITSYS/include");
  gROOT->ProcessLine(".L ../CreateDatacards/include/tdrstyle.cc");
  setTDRStyle(false);
  gStyle->SetPadLeftMargin(0.16);

  TString filepath;
  if (sqrts==7) {
    filepath = filePath7TeV;
  } else if (sqrts==8) {
    filepath = filePath8TeV;
  }

  TChain* tree = new TChain("SelectedTree");
  tree->Add( filepath+ "/" + (schannel=="2e2mu"?"2mu2e":schannel) + "/HZZ4lTree_ZZTo*.root");


  RooRealVar* MC_weight = new RooRealVar("MC_weight","MC_weight",0.,2.) ; 
  RooRealVar* ZZMass = new RooRealVar("ZZMass","ZZMass",100.,1000.);
  RooRealVar* NJets30 = new RooRealVar("NJets30","NJets30",0.,100.);
  RooArgSet ntupleVarSet(*ZZMass,*NJets30,*MC_weight);
  RooDataSet *set = new RooDataSet("set","set",ntupleVarSet,WeightVar("MC_weight"));

  Float_t myMC,myMass;
  Short_t myNJets;
  int nentries = tree->GetEntries();

  tree->SetBranchAddress("ZZMass",&myMass);
  tree->SetBranchAddress("MC_weight",&myMC);
  tree->SetBranchAddress("NJets30",&myNJets);

  for(int i =0;i<nentries;i++) {
    tree->GetEntry(i);
    if(VBFtag==1 && myNJets<2)continue;
    if(VBFtag==0 && myNJets>1)continue;

    ntupleVarSet.setRealValue("ZZMass",myMass);
    ntupleVarSet.setRealValue("MC_weight",myMC);
    ntupleVarSet.setRealValue("NJets30",(double)myNJets);

    set->add(ntupleVarSet, myMC);
  }

  double totalweight = 0.;
  double totalweight_z = 0.;
  for (int i=0 ; i<set->numEntries() ; i++) { 
    //set->get(i) ; 
    RooArgSet* row = set->get(i) ;
    //row->Print("v");
    totalweight += set->weight();
    if (row->getRealValue("ZZMass") < 200) totalweight_z += set->weight();
  } 
  cout << "nEntries: " << set->numEntries() << ", totalweight: " << totalweight << ", totalweight_z: " << totalweight_z << endl;

  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
	
  //// ---------------------------------------
  //Background
  RooRealVar CMS_qqzzbkg_a0("CMS_qqzzbkg_a0","CMS_qqzzbkg_a0",115.3,0.,200.);
  RooRealVar CMS_qqzzbkg_a1("CMS_qqzzbkg_a1","CMS_qqzzbkg_a1",21.96,0.,200.);
  RooRealVar CMS_qqzzbkg_a2("CMS_qqzzbkg_a2","CMS_qqzzbkg_a2",122.8,0.,200.);
  RooRealVar CMS_qqzzbkg_a3("CMS_qqzzbkg_a3","CMS_qqzzbkg_a3",0.03479,0.,1.);
  RooRealVar CMS_qqzzbkg_a4("CMS_qqzzbkg_a4","CMS_qqzzbkg_a4",185.5,0.,200.);
  RooRealVar CMS_qqzzbkg_a5("CMS_qqzzbkg_a5","CMS_qqzzbkg_a5",12.67,0.,200.);
  RooRealVar CMS_qqzzbkg_a6("CMS_qqzzbkg_a6","CMS_qqzzbkg_a6",34.81,0.,100.);
  RooRealVar CMS_qqzzbkg_a7("CMS_qqzzbkg_a7","CMS_qqzzbkg_a7",0.1393,0.,1.);
  RooRealVar CMS_qqzzbkg_a8("CMS_qqzzbkg_a8","CMS_qqzzbkg_a8",66.,0.,200.);
  RooRealVar CMS_qqzzbkg_a9("CMS_qqzzbkg_a9","CMS_qqzzbkg_a9",0.07191,0.,1.);
  RooRealVar CMS_qqzzbkg_a10("CMS_qqzzbkg_a10","CMS_qqzzbkg_a10",94.11,0.,200.);
  RooRealVar CMS_qqzzbkg_a11("CMS_qqzzbkg_a11","CMS_qqzzbkg_a11",-5.111,-100.,100.);
  RooRealVar CMS_qqzzbkg_a12("CMS_qqzzbkg_a12","CMS_qqzzbkg_a12",4834,0.,10000.);
  RooRealVar CMS_qqzzbkg_a13("CMS_qqzzbkg_a13","CMS_qqzzbkg_a13",0.2543,0.,1.);
	
  if (channel == 1){
    ///* 4mu
    CMS_qqzzbkg_a0.setVal(103.854);
    CMS_qqzzbkg_a1.setVal(10.0718);
    CMS_qqzzbkg_a2.setVal(117.551);
    CMS_qqzzbkg_a3.setVal(0.0450287);
    CMS_qqzzbkg_a4.setVal(185.262);
//.........这里部分代码省略.........

// The actual job
void backgroundFits_ggzz_1Dw(int channel, int sqrts, int VBFtag)
{
  TString schannel;
  if      (channel == 1) schannel = "4mu";
  else if (channel == 2) schannel = "4e";
  else if (channel == 3) schannel = "2e2mu";
  else cout << "Not a valid channel: " << schannel << endl;

  TString ssqrts = (long) sqrts + TString("TeV");

  cout << "schannel = " << schannel << "  sqrts = " << sqrts << " VBFtag = "<< VBFtag << endl;

  TString outfile;
  outfile = "CardFragments/ggzzBackgroundFit_" + ssqrts + "_" + schannel + "_" + Form("%d",int(VBFtag)) + ".txt";
  ofstream of(outfile,ios_base::out);

  gSystem->AddIncludePath("-I$ROOFITSYS/include");
  gROOT->ProcessLine(".L ../CreateDatacards/include/tdrstyle.cc");
  setTDRStyle(false);
  gStyle->SetPadLeftMargin(0.16);
	
  TString filepath;filepath.Form("AAAOK/ZZ%s/ZZ4lAnalysis.root",schannel.Data());
  TFile *f = TFile::Open(filepath);
  TTree *tree = f->Get("ZZTree/candTree");

  RooRealVar* MC_weight = new RooRealVar("MC_weight","MC_weight",0.,2.) ; 
  RooRealVar* ZZMass = new RooRealVar("ZZMass","ZZMass",100,100.,1000.);
  RooRealVar* NJets30 = new RooRealVar("NJets30","NJets30",0.,5.);
  RooArgSet ntupleVarSet(*ZZMass,*NJets30,*MC_weight);
  RooDataSet *set = new RooDataSet("set","set",ntupleVarSet,WeightVar("MC_weight"));
  //RooArgSet ntupleVarSet(*ZZMass,*NJets30);  
  //RooDataSet *set = new RooDataSet("set","set",ntupleVarSet);

  Float_t myMC,myMass;
  Int_t myNJets;
  int nentries = tree->GetEntries();

  Float_t myPt,myJetPt,myJetEta,myJetPhi,myJetMass,myFisher;
  Int_t myExtralep,myBJets;
  tree->SetBranchAddress("ZZMass",&myMass);
  tree->SetBranchAddress("genHEPMCweight",&myMC);
  tree->SetBranchAddress("nCleanedJetsPt30",&myNJets);
  tree->SetBranchAddress("ZZPt",&myPt);
  tree->SetBranchAddress("nExtraLep",&myExtralep);
  tree->SetBranchAddress("nCleanedJetsPt30BTagged",&myBJets);
  tree->SetBranchAddress("DiJetDEta",&myFisher);

  for(int i =0;i<nentries;i++) {
    tree->GetEntry(i);
    if(myMass<100.)continue;
    int cat = category(myExtralep,myPt, myMass,myNJets, myBJets,/* jetpt, jeteta, jetphi, jetmass,*/myFisher);
    if(VBFtag != cat )continue;

    ntupleVarSet.setRealValue("ZZMass",myMass);
    ntupleVarSet.setRealValue("MC_weight",myMC);
    ntupleVarSet.setRealValue("NJets30",(double)cat);

    set->add(ntupleVarSet, myMC);
  }

  //RooRealVar* ZZLD = new RooRealVar("ZZLD","ZZLD",0.,1.);
  //char cut[10];
  //sprintf(cut,"ZZLD>0.5");
  //RooDataSet* set = new RooDataSet("set","set",tree,RooArgSet(*ZZMass,*MC_weight,*ZZLD),cut,"MC_weight");

  double totalweight = 0.;
  for (int i=0 ; i<set->numEntries() ; i++) { 
    set->get(i) ; 
    totalweight += set->weight();
    //cout << CMS_zz4l_mass->getVal() << " = " << set->weight() << endl ; 
  } 
  cout << "nEntries: " << set->numEntries() << ", totalweight: " << totalweight << endl;
		
  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
	
  //// ---------------------------------------
  //Background
  RooRealVar CMS_qqzzbkg_a0("CMS_qqzzbkg_a0","CMS_qqzzbkg_a0",115.3,0.,200.);
  RooRealVar CMS_qqzzbkg_a1("CMS_qqzzbkg_a1","CMS_qqzzbkg_a1",21.96,0.,200.);
  RooRealVar CMS_qqzzbkg_a2("CMS_qqzzbkg_a2","CMS_qqzzbkg_a2",122.8,0.,200.);
  RooRealVar CMS_qqzzbkg_a3("CMS_qqzzbkg_a3","CMS_qqzzbkg_a3",0.03479,0.,1.);
  RooRealVar CMS_qqzzbkg_a4("CMS_qqzzbkg_a4","CMS_qqzzbkg_a4",185.5,0.,200.);
  RooRealVar CMS_qqzzbkg_a5("CMS_qqzzbkg_a5","CMS_qqzzbkg_a5",12.67,0.,200.);
  RooRealVar CMS_qqzzbkg_a6("CMS_qqzzbkg_a6","CMS_qqzzbkg_a6",34.81,0.,100.);
  RooRealVar CMS_qqzzbkg_a7("CMS_qqzzbkg_a7","CMS_qqzzbkg_a7",0.1393,0.,1.);
  RooRealVar CMS_qqzzbkg_a8("CMS_qqzzbkg_a8","CMS_qqzzbkg_a8",66.,0.,200.);
  RooRealVar CMS_qqzzbkg_a9("CMS_qqzzbkg_a9","CMS_qqzzbkg_a9",0.07191,0.,1.);
	
  RooggZZPdf_v2* bkg_ggzz = new RooggZZPdf_v2("bkg_ggzz","bkg_ggzz",*ZZMass,
					      CMS_qqzzbkg_a0,CMS_qqzzbkg_a1,CMS_qqzzbkg_a2,CMS_qqzzbkg_a3,CMS_qqzzbkg_a4,
					      CMS_qqzzbkg_a5,CMS_qqzzbkg_a6,CMS_qqzzbkg_a7,CMS_qqzzbkg_a8,CMS_qqzzbkg_a9);
	
  //// ---------------------------------------
	
  RooFitResult *r1 = bkg_ggzz->fitTo( *set, Save(kTRUE), SumW2Error(kTRUE) );//, Save(kTRUE), SumW2Error(kTRUE)) ;

  cout << endl;
  cout << "------- Parameters for " << schannel << " sqrts=" << sqrts << endl;
  cout << "  a0_bkgd = " << CMS_qqzzbkg_a0.getVal() << endl;
//.........这里部分代码省略.........

int main() {

  TFile *tf = TFile::Open("root/MassFitResult.root");
  RooWorkspace *w = (RooWorkspace*)tf->Get("w");

  RooDataSet *data = (RooDataSet*)w->data("Data2012HadronTOS");
  //w->loadSnapshot("bs2kstkst_mc_pdf_fit");

  //RooRealVar *bs2kstkst_l      = new RooRealVar("bs2kstkst_l"      , "", -5., -20., 0.);
  //RooConstVar *bs2kstkst_zeta  = new RooConstVar("bs2kstkst_zeta" , "", 0.);
  //RooConstVar *bs2kstkst_fb    = new RooConstVar("bs2kstkst_fb"   , "", 0.);
  //RooRealVar *bs2kstkst_sigma  = new RooRealVar("bs2kstkst_sigma"  , "", 15, 10, 100);
  //RooRealVar *bs2kstkst_mu     = new RooRealVar("bs2kstkst_mu"     , "", 5350, 5400 );
  //RooRealVar *bs2kstkst_a      = new RooRealVar("bs2kstkst_a"      , "", 2.5,0,10);
  //RooRealVar *bs2kstkst_n      = new RooRealVar("bs2kstkst_n"      , "", 2.5,0,10);
  //RooRealVar *bs2kstkst_a2     = new RooRealVar("bs2kstkst_a2"     , "", 2.5,0,10);
  //RooRealVar *bs2kstkst_n2     = new RooRealVar("bs2kstkst_n2"     , "", 2.5,0,10);

  //RooIpatia2 *sig = new RooIpatia2("sig","",*w->var("B_s0_DTF_B_s0_M"), *bs2kstkst_l, *bs2kstkst_zeta, *bs2kstkst_fb, *bs2kstkst_sigma, *bs2kstkst_mu, *bs2kstkst_a, *bs2kstkst_n, *bs2kstkst_a2, *bs2kstkst_n2);
  //RooAbsPdf *sig = (RooAbsPdf*)w->pdf("bs2kstkst_mc_pdf");
  RooIpatia2 *sig = new RooIpatia2("bs2kstkst_mc_pdf","bs2kstkst_mc_pdf",*w->var("B_s0_DTF_B_s0_M"),*w->var("bs2kstkst_l"),*w->var("bs2kstkst_zeta"),*w->var("bs2kstkst_fb"),*w->var("bs2kstkst_sigma"),*w->var("bs2kstkst_mu"),*w->var("bs2kstkst_a"),*w->var("bs2kstkst_n"),*w->var("bs2kstkst_a2"),*w->var("bs2kstkst_n2"));

  RooAbsPdf *bkg = (RooAbsPdf*)w->pdf("bkg_pdf_HadronTOS2012");

  RooRealVar *sY = (RooRealVar*)w->var("bs2kstkst_y_HadronTOS2012");
  RooRealVar *bY = (RooRealVar*)w->var("bkg_y_HadronTOS2012");

  cout << sig << bkg << sY << bY << endl;

  RooAddPdf *pdf = new RooAddPdf("test","test", RooArgList(*sig,*bkg), RooArgList(*sY,*bY) );

  pdf->fitTo(*data, Extended() );

  // my sw
  double syVal = sY->getVal();
  double byVal = bY->getVal();

  // loop events
  int numevents = data->numEntries();

  sY->setVal(0.);
  bY->setVal(0.);

  RooArgSet *pdfvars = pdf->getVariables();

  vector<double> fsvals;
  vector<double> fbvals;

  for ( int ievt=0; ievt<numevents; ievt++ ) {

    RooStats::SetParameters(data->get(ievt), pdfvars);

    sY->setVal(1.);
    double f_s = pdf->getVal( RooArgSet(*w->var("B_s0_DTF_B_s0_M")) );
    fsvals.push_back(f_s);
    sY->setVal(0.);

    bY->setVal(1.);
    double f_b = pdf->getVal( RooArgSet(*w->var("B_s0_DTF_B_s0_M")) );
    fbvals.push_back(f_b);
    bY->setVal(0.);

    //cout << f_s << " " << f_b << endl;

  }

  TMatrixD covInv(2,2);
  covInv[0][0] = 0;
  covInv[0][1] = 0;
  covInv[1][0] = 0;
  covInv[1][1] = 0;

  for ( int ievt=0; ievt<numevents; ievt++ ) {
    data->get(ievt);
    double dsum=0;
    dsum += fsvals[ievt] * syVal;
    dsum += fbvals[ievt] * byVal;

    covInv[0][0] += fsvals[ievt]*fsvals[ievt] / (dsum*dsum);
    covInv[0][1] += fsvals[ievt]*fbvals[ievt] / (dsum*dsum);
    covInv[1][0] += fbvals[ievt]*fsvals[ievt] / (dsum*dsum);
    covInv[1][1] += fbvals[ievt]*fbvals[ievt] / (dsum*dsum);

  }

  covInv.Print();

  cout << covInv.Determinant() << endl;

  TMatrixD covMatrix(TMatrixD::kInverted,covInv);

  covMatrix.Print();

  RooStats::SPlot *sD = new RooStats::SPlot("sD","sD",*data,pdf,RooArgSet(*sY,*bY),RooArgSet(*w->var("eventNumber")));
}

void IntervalExamples()
{

   // Time this macro
   TStopwatch t;
   t.Start();


   // set RooFit random seed for reproducible results
   RooRandom::randomGenerator()->SetSeed(3001);

   // make a simple model via the workspace factory
   RooWorkspace* wspace = new RooWorkspace();
   wspace->factory("Gaussian::normal(x[-10,10],mu[-1,1],sigma[1])");
   wspace->defineSet("poi","mu");
   wspace->defineSet("obs","x");

   // specify components of model for statistical tools
   ModelConfig* modelConfig = new ModelConfig("Example G(x|mu,1)");
   modelConfig->SetWorkspace(*wspace);
   modelConfig->SetPdf( *wspace->pdf("normal") );
   modelConfig->SetParametersOfInterest( *wspace->set("poi") );
   modelConfig->SetObservables( *wspace->set("obs") );

   // create a toy dataset
   RooDataSet* data = wspace->pdf("normal")->generate(*wspace->set("obs"),100);
   data->Print();

   // for convenience later on
   RooRealVar* x = wspace->var("x");
   RooRealVar* mu = wspace->var("mu");

   // set confidence level
   double confidenceLevel = 0.95;

   // example use profile likelihood calculator
   ProfileLikelihoodCalculator plc(*data, *modelConfig);
   plc.SetConfidenceLevel( confidenceLevel);
   LikelihoodInterval* plInt = plc.GetInterval();

   // example use of Feldman-Cousins
   FeldmanCousins fc(*data, *modelConfig);
   fc.SetConfidenceLevel( confidenceLevel);
   fc.SetNBins(100); // number of points to test per parameter
   fc.UseAdaptiveSampling(true); // make it go faster

   // Here, we consider only ensembles with 100 events
   // The PDF could be extended and this could be removed
   fc.FluctuateNumDataEntries(false);

   // Proof
   //  ProofConfig pc(*wspace, 4, "workers=4", kFALSE);    // proof-lite
   //ProofConfig pc(w, 8, "localhost");    // proof cluster at "localhost"
   //  ToyMCSampler* toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
   //  toymcsampler->SetProofConfig(&pc);     // enable proof

   PointSetInterval* interval = (PointSetInterval*) fc.GetInterval();


   // example use of BayesianCalculator
   // now we also need to specify a prior in the ModelConfig
   wspace->factory("Uniform::prior(mu)");
   modelConfig->SetPriorPdf(*wspace->pdf("prior"));

   // example usage of BayesianCalculator
   BayesianCalculator bc(*data, *modelConfig);
   bc.SetConfidenceLevel( confidenceLevel);
   SimpleInterval* bcInt = bc.GetInterval();

   // example use of MCMCInterval
   MCMCCalculator mc(*data, *modelConfig);
   mc.SetConfidenceLevel( confidenceLevel);
   // special options
   mc.SetNumBins(200);        // bins used internally for representing posterior
   mc.SetNumBurnInSteps(500); // first N steps to be ignored as burn-in
   mc.SetNumIters(100000);    // how long to run chain
   mc.SetLeftSideTailFraction(0.5); // for central interval
   MCMCInterval* mcInt = mc.GetInterval();

   // for this example we know the expected intervals
   double expectedLL = data->mean(*x)
      + ROOT::Math::normal_quantile(  (1-confidenceLevel)/2,1)
      / sqrt(data->numEntries());
   double expectedUL = data->mean(*x)
      + ROOT::Math::normal_quantile_c((1-confidenceLevel)/2,1)
      / sqrt(data->numEntries()) ;

   // Use the intervals
   std::cout << "expected interval is [" <<
      expectedLL << ", " <<
      expectedUL << "]" << endl;

   cout << "plc interval is [" <<
      plInt->LowerLimit(*mu) << ", " <<
      plInt->UpperLimit(*mu) << "]" << endl;

   std::cout << "fc interval is ["<<
      interval->LowerLimit(*mu) << " , "  <<
      interval->UpperLimit(*mu) << "]" << endl;

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

void embeddedToysVBF_1DKD(int nEvts=50, int nToys=100,
		  sample mySample = kScalar_fa3_0){
  

  RooRealVar* kd = new RooRealVar("psMELA","psMELA",0,1);
  kd->setBins(50);
  RooPlot* kdframe1 = kd->frame();
  
  // 0- template
  TFile f1("ggH2j_KDdistribution_ps.root", "READ"); 
  TH2F *h_KD_ps = (TH2F*)f1.Get("h_KD");
  h_KD_ps->SetName("h_KD_ps");
  RooDataHist rdh_KD_ps("rdh_KD_ps","rdh_KD_ps",RooArgList(*kd),h_KD_ps);
  RooHistPdf pdf_KD_ps("pdf_KD_ps","pdf_KD_ps",RooArgList(*kd),rdh_KD_ps); 

  // 0+ template
  TFile f2("ggH2j_KDdistribution_sm.root", "READ"); 
  TH2F *h_KD_sm = (TH2F*)f2.Get("h_KD");
  h_KD_sm->SetName("h_KD_sm");
  RooDataHist rdh_KD_sm("rdh_KD_sm","rdh_KD_sm",RooArgList(*kd),h_KD_sm);
  RooHistPdf pdf_KD_sm("pdf_KD_sm","pdf_KD_sm",RooArgList(*kd),rdh_KD_sm); 

  RooRealVar rrv_fa3("fa3","fa3",0.5,0.,1.);  //free parameter of the model
  RooAddPdf model("model","ps+sm",pdf_KD_ps,pdf_KD_sm,rrv_fa3);  
  rrv_fa3.setConstant(kFALSE);

 TCanvas* c = new TCanvas("modelPlot","modelPlot",400,400);
  rdh_KD_ps.plotOn(kdframe1,LineColor(kBlack));
  pdf_KD_ps.plotOn(kdframe1,LineColor(kBlack));
  rdh_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  pdf_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  model.plotOn(kdframe1,LineColor(kRed));
  kdframe1->Draw();
  c->SaveAs("model1DPlot.png");
  c->SaveAs("model1DPlot.eps");

  double fa3Val=-99;
  if (mySample == kScalar_fa3_0)
    fa3Val=0.;
  else if (mySample == kScalar_fa3_1)
    fa3Val=1;
  else{
    cout<<"fa3Val not correct!"<<endl;
      return 0;
  }

  TCanvas* c = new TCanvas("modelPlot","modelPlot",400,400);
  rdh_KD_ps.plotOn(kdframe1,LineColor(kBlack));
  pdf_KD_ps.plotOn(kdframe1,LineColor(kBlack));
  rdh_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  pdf_KD_sm.plotOn(kdframe1,LineColor(kBlue));
  model.plotOn(kdframe1,LineColor(kRed));
  kdframe1->Draw();

  
  TChain* myChain = new TChain("newTree");
  myChain->Add(inputFileNames[mySample]);
  
  if(!myChain || myChain->GetEntries()<=0) {
    cout<<"error in the tree"<<endl;
    return 0;
  }
  
   RooDataSet* data = new RooDataSet("data","data",myChain,RooArgSet(*kd),"");

    cout << "Number of events in data: " << data->numEntries() << endl;
  
  // initialize tree to save toys to 
  TTree* results = new TTree("results","toy results");
  
  double fa3,fa3Error, fa3Pull;
  double fa2,fa2Error, fa2Pull;
  double phia2, phia2Error, phia2Pull;
  double phia3, phia3Error, phia3Pull;

  results->Branch("fa3",&fa3,"fa3/D");
  results->Branch("fa3Error",&fa3Error,"fa3Error/D");
  results->Branch("fa3Pull",&fa3Pull,"fa3Pull/D");

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

  RooDataSet* toyData;
  int embedTracker=0;
  RooArgSet *tempEvent;

  RooFitResult *toyfitresults;
  RooRealVar *r_fa3;

  for(int i = 0 ; i<nToys ; i++){
    cout <<i<<"<-----------------------------"<<endl;
    //if(toyData) delete toyData;
    toyData = new RooDataSet("toyData","toyData",RooArgSet(*kd));

    if(nEvts+embedTracker > data->sumEntries()){
      cout << "Playground::generate() - ERROR!!! Playground::data does not have enough events to fill toy!!!!  bye :) " << endl;
      toyData = NULL;
      abort();
      return 0;
    }

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

void track_pt(const int charge)
{
  if (charge==1)
  char  ch[20] = "plus"; 
  else if (charge==-1)
  char ch[20] = "minus";
  ofstream txtfile;
  char txtfname[100];
  char histfname[100];
  sprintf(txtfname,"pt_%s.txt",ch);
  sprintf(histfname,"pt_%s.png",ch);
  txtfile.open(txtfname);
  txtfile << fixed << setprecision(4);
  TCanvas *myCan=new TCanvas("myCan","myCan",800,600);
  gStyle->SetLineWidth(2.);
  gStyle->SetLabelSize(0.04,"xy");
  gStyle->SetTitleSize(0.05,"xy");
  gStyle->SetTitleOffset(1.1,"x");
  gStyle->SetTitleOffset(1.2,"y");
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadRightMargin(0.1);
  gStyle->SetPadBottomMargin(0.16);
  gStyle->SetPadLeftMargin(0.12);

  myCan->SetGrid();
  TLegend* Lgd = new TLegend(.8, .25,.9,.35);
  if (charge==1){
    TFile *f_MC= new TFile("TnP_Tracking_dr030e030_MCptplus.root","read");
    TFile *f_RD= new TFile("TnP_Tracking_dr030e030_RDptplus.root","read");
  }else if (charge==-1){
    TFile *f_MC= new TFile("TnP_Tracking_dr030e030_MCptminus.root","read");
    TFile *f_RD= new TFile("TnP_Tracking_dr030e030_RDptminus.root","read");
  }
  line = new TLine(15,1,80,1);
  line->SetLineStyle(2);
  line->SetLineWidth(3);
  
  TPaveText *title = new TPaveText(.1,1,.95,.95,"NDC");
  title->SetFillStyle(0);
  title->SetBorderSize(0);
  title->SetTextSize(0.04);
  title->AddText("CMS Preliminary, 18.9 pb^{-1} at #sqrt{s}=8 TeV");

  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTreeSta/eff_pt_dr030e030/fit_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin],ErrMC[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &ptMC=pointMC["pt"],&effMC = pointMC["efficiency"];
    XMC[i]=ptMC.getVal();
    XMCerrL[i]=-ptMC.getAsymErrorLo();
    XMCerrH[i]=ptMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
    ErrMC[i]=TMath::Max(YMCerrLo[i],YMCerrHi[i]);
  }
  grMC=new TGraphAsymmErrors(Nbin,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->SetMarkerStyle(21);
  grMC->SetMinimum(0.7);
  grMC->SetMaximum(1.11);
  grMC->GetXaxis()->SetNdivisions(505);
  grMC->GetXaxis()->SetTitle("Muon p_{T} [GeV/c]");
  grMC->GetYaxis()->SetTitle("Tracking Efficiency");
  grMC->Draw("AP");

  RooDataSet *datasetRD = (RooDataSet*)f_RD->Get("tpTreeSta/eff_pt_dr030e030/fit_eff");
  double XRD[Nbin],XRDerrL[Nbin],XRDerrH[Nbin],YRD[Nbin],YRDerrLo[Nbin],YRDerrHi[Nbin],ErrRD[Nbin];
  for(int i(0); i<datasetRD->numEntries();i++)
  {
    const RooArgSet &pointRD=*datasetRD->get(i);
    RooRealVar &ptRD=pointRD["pt"],&effRD = pointRD["efficiency"];
    XRD[i]=ptRD.getVal();
    XRDerrL[i]=-ptRD.getAsymErrorLo();
    XRDerrH[i]=ptRD.getAsymErrorHi();
    YRD[i]=effRD.getVal();
    YRDerrLo[i]=-effRD.getAsymErrorLo();
    YRDerrHi[i]=effRD.getAsymErrorHi();
    ErrRD[i]=TMath::Max(YRDerrLo[i],YRDerrHi[i]);
  }

  double SF[Nbin],SFerr[Nbin],SFerrL[Nbin],SFerrH[Nbin];
  for(int i(0); i<datasetRD->numEntries();i++)
  {
    SF[i]=YRD[i]/YMC[i];
    SFerrL[i]=YRD[i]*sqrt(YMCerrLo[i]*YMCerrLo[i]/(YMC[i]*YMC[i])+YRDerrLo[i]*YRDerrLo[i]/(YRD[i]*YRD[i]))/YMC[i];
    SFerrH[i]=YRD[i]*sqrt(YMCerrHi[i]*YMCerrHi[i]/(YMC[i]*YMC[i])+YRDerrHi[i]*YRDerrHi[i]/(YRD[i]*YRD[i]))/YMC[i];
    SFerr[i]=TMath::Max(SFerrL[i],SFerrH[i]);
    txtfile << i << "\t" << YMC[i] << "\t" << ErrMC[i] << "\t" << YRD[i] << "\t" << ErrRD[i] << "\t" << SF[i] << "\t" << SFerr[i] << endl;
  }

  grRD=new TGraphAsymmErrors(Nbin,XRD,YRD,XRDerrL,XRDerrH,YRDerrLo,YRDerrHi);
  grRD->SetLineColor(kBlack);
  grRD->SetMarkerColor(kBlack);
  
  grSF=new TGraphAsymmErrors(Nbin,XRD,SF,0,0,SFerrL,SFerrH);
  grSF->SetLineColor(8);
//.........这里部分代码省略.........

void glbToId_eta()
{
  ofstream txtfile;
  char txtfname[100];
  char histfname[100];
  sprintf(txtfname,"eta_plus.txt");
  sprintf(histfname,"eta_plus.png");
  //sprintf(txtfname,"eta_minus.txt");
  //sprintf(histfname,"eta_minus.png");
  txtfile.open(txtfname);
  txtfile << fixed << setprecision(4);
  TCanvas *myCan=new TCanvas("myCan","myCan",800,600);
  gStyle->SetLineWidth(2.);
  gStyle->SetLabelSize(0.04,"xy");
  gStyle->SetTitleSize(0.05,"xy");
  gStyle->SetTitleOffset(1.1,"x");
  gStyle->SetTitleOffset(1.2,"y");
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadRightMargin(0.1);
  gStyle->SetPadBottomMargin(0.16);
  gStyle->SetPadLeftMargin(0.12);

  myCan->SetGrid();
  TLegend* Lgd = new TLegend(.8, .25,.9,.35);
  
  TFile *f_MCndof2= new TFile("TnP_GlbToID_MCetaplus_Wpteta_eta.root","read");
  TFile *f_MCndof4= new TFile("TnP_GlbToID_MCetaplus_ndof4_Wpteta_eta.root","read");

  //TFile *f_MCndof2= new TFile("TnP_GlbToID_MCetaminus_Wpteta_eta.root","read");
  //TFile *f_MCndof4= new TFile("TnP_GlbToID_MCetaminus_ndof4_Wpteta_eta.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MCndof2->Get("tpTree/Wpteta_eta/fit_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin],ErrMC[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &etaMC=pointMC["eta"],&effMC = pointMC["efficiency"];
    XMC[i]=etaMC.getVal();
    XMCerrL[i]=-etaMC.getAsymErrorLo();
    XMCerrH[i]=etaMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
    ErrMC[i]=TMath::Max(YMCerrLo[i],YMCerrHi[i]);
  }

  grMC=new TGraphAsymmErrors(Nbin,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->SetMinimum(0.5);
  grMC->SetMaximum(1.1);
  grMC->GetXaxis()->SetNdivisions(505);
  grMC->GetXaxis()->SetTitle("Muon #eta");
  grMC->GetYaxis()->SetTitle("ID+ISO Efficiency");
  grMC->Draw("AP");

  RooDataSet *datasetRD = (RooDataSet*)f_MCndof4->Get("tpTree/Wpteta_eta/fit_eff");
  double XRD[Nbin],XRDerrL[Nbin],XRDerrH[Nbin],YRD[Nbin],YRDerrLo[Nbin],YRDerrHi[Nbin],ErrRD[Nbin];
  for(int i(0); i<datasetRD->numEntries();i++)
  {
    const RooArgSet &pointRD=*datasetRD->get(i);
    RooRealVar &etaRD=pointRD["eta"],&effRD = pointRD["efficiency"];
    XRD[i]=etaRD.getVal();
    XRDerrL[i]=-etaRD.getAsymErrorLo();
    XRDerrH[i]=etaRD.getAsymErrorHi();
    YRD[i]=effRD.getVal();
    YRDerrLo[i]=-effRD.getAsymErrorLo();
    YRDerrHi[i]=effRD.getAsymErrorHi();
    ErrRD[i]=TMath::Max(YRDerrLo[i],YRDerrHi[i]);
  }
  txtfile << "Bins \t MC ndof>2\t\t\t MC ndof>4 \t\t\t Scale Factor ndof4/ndof2 "  << endl;
  for(int i(0); i<datasetRD->numEntries();i++)
  {
    txtfile << i << "\t" << YMC[i] << "+/-" << ErrMC[i] << "\t\t" << YRD[i] << "+/-" << ErrRD[i] << "\t\t" << YRD[i]/YMC[i] << "+/-" << YRD[i]*sqrt(ErrMC[i]*ErrMC[i]/(YMC[i]*YMC[i])+ErrRD[i]*ErrRD[i]/(YRD[i]*YRD[i]))/YMC[i] << endl;
  }

  grRD=new TGraphAsymmErrors(Nbin,XRD,YRD,XRDerrL,XRDerrH,YRDerrLo,YRDerrHi);
  grRD->SetLineColor(kBlack);
  grRD->SetMarkerColor(kBlack);
  Lgd->AddEntry(grMC,"MC ndof>2","pl");
  Lgd->AddEntry(grRD,"MC ndof>4","pl");
  Lgd->SetFillStyle(0);
  Lgd->Draw();
  grRD->Draw("PSAME");
  
  myCan->SaveAs(histfname);
  txtfile.close();
}