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

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 FitterUtils::PlotShape2D(RooDataSet& originDataSet, RooDataSet& genDataSet, RooAbsPdf& shape, string plotsfile, string canvName, RooRealVar& B_plus_M, RooRealVar& misPT)
{
   //**************Prepare TFile to save the plots

   TFile f2(plotsfile.c_str(), "UPDATE");

   //**************Plot Signal Zero Gamma

   TH2F* th2fKey = (TH2F*)shape.createHistogram("th2Shape", B_plus_M, Binning(20), YVar(misPT, Binning(20)));
   cout<<genDataSet.sumEntries()<<endl;
   TH2F* th2fGen = (TH2F*)genDataSet.createHistogram("th2fGen", B_plus_M, Binning(20), YVar(misPT, Binning(20)));

   RooPlot* plotM = B_plus_M.frame();
   originDataSet.plotOn(plotM);
   shape.plotOn(plotM);

   RooPlot* plotMisPT = misPT.frame();
   originDataSet.plotOn(plotMisPT);
   shape.plotOn(plotMisPT);

   TCanvas canv(canvName.c_str(), canvName.c_str(), 800, 800);
   canv.Divide(2,2);
   canv.cd(1); th2fGen->Draw("lego");
   canv.cd(2); th2fKey->Draw("surf");
   canv.cd(3); plotM->Draw();
   canv.cd(4); plotMisPT->Draw();

   canv.Write();

   f2.Close();
}

void reweight_eta_1d(TH1F* weight_eta,TH1F* weight_etao,TH2F*weight_etan,TH2F* weight_eta2o,TH2F* weight_etanr,TH2F*weight_eta2,RooDataSet **dset, RooDataSet *dsetdestination, int numvar){

  if (!(*dset)) return;

  TH1F *hnum = new TH1F("hnum","hnum",n_etabins_forreweighting,etabins_forreweighting);
 TH1F *hden = new TH1F("hden","hden",n_etabins_forreweighting,etabins_forreweighting);
//  TH1F *hnum = new TH1F("hnum","hnum",25,0.,2.5);
//  TH1F *hden = new TH1F("hden","hden",25,0.,2.5);
  hnum->Sumw2();
  hden->Sumw2();

  const char* etaname=Form("rooeta%d",numvar);

  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue(etaname)),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue(etaname)),dsetdestination->store()->weight(i));
  }
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());

  hnum->Divide(hden);
  TH1F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_etarew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float eta = args.getRealValue(etaname);
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(eta)));
    if(debug){ 
		weight_eta->Fill(neww);	
		weight_etao->Fill(oldw);
		weight_etan->Fill(h->FindBin(fabs(eta)),neww);	
		weight_eta2o->Fill(h->FindBin(fabs(eta)),oldw);
		if(oldw!=0 && neww!=0)weight_etanr->Fill(h->FindBin(fabs(eta)),oldw/neww);
		else {weight_etanr->Fill(-10,1);}
	   // weight_pt2->Fill(pt,neww/oldw);
	   if(oldw!=0 && neww!=0)weight_eta2->Fill(fabs(eta),oldw/neww);
	   else {weight_eta2->Fill(-10,1);}
	}  
	  newdset->add(args,neww);
	  
	  }


  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;

  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Eta 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;

  delete old_dset;

};

void MakeSpinPlots::DrawSpinBackground(TString tag, TString mcName,bool signal){
  bool drawSM = (smName!="" && smName!=mcName);

  TCanvas cv;
  double thisN  = ws->data(mcName+"_Combined")->reduce(TString("evtcat==evtcat::")+tag)->sumEntries();
  float norm = thisN; //607*lumi/12.*thisN/(totEB+totEE);
  cout << norm <<endl;
  if(signal) norm = ws->data(Form("Data_%s_%s_sigWeight",tag.Data(),mcName.Data()))->sumEntries();
  RooPlot *frame = ws->var("cosT")->frame(0,1,5);

  RooDataSet* bkgWeight = (RooDataSet*)ws->data(Form("Data_%s_%s_bkgWeight",tag.Data(),mcName.Data()));
  RooDataSet* tmp = (RooDataSet*)ws->data("Data_Combined")->reduce(TString("((mass>115 && mass<120) || (mass>130 && mass<135)) && evtcat==evtcat::")+tag);
  tmp->plotOn(frame,RooFit::Rescale(norm/tmp->sumEntries()));
  cout << "b" <<endl;
  ws->pdf(Form("%s_FIT_%s_cosTpdf",mcName.Data(),tag.Data()))->plotOn(frame,RooFit::LineColor(kGreen),RooFit::Normalization(norm/tmp->sumEntries()));
  if(drawSM) ws->pdf(Form("%s_FIT_%s_cosTpdf",smName.Data(),tag.Data()))->plotOn(frame,RooFit::LineColor(kRed),RooFit::Normalization(norm/tmp->sumEntries()));
  cout << "c   " <<bkgWeight <<endl;

  bkgWeight->plotOn(frame,RooFit::Rescale(norm/bkgWeight->sumEntries()),RooFit::MarkerColor(kBlue) );  
  if(signal){
    cout << "d" <<endl;
      
    ws->data(Form("Data_%s_%s_sigWeight",tag.Data(),mcName.Data()))->plotOn(frame,RooFit::MarkerStyle(4));
  }
  cout << "d" <<endl;
  
  frame->SetMaximum(frame->GetMaximum()*(signal?0.8:0.4)*norm/tmp->sumEntries());
  frame->SetMinimum(-1*frame->GetMaximum());
  TLegend l(0.6,0.2,0.95,0.45);
  l.SetFillColor(0);
  l.SetBorderSize(0);
  l.SetHeader(tag);
  l.AddEntry(frame->getObject(0),"Data m#in [115,120]#cup[130,135]","p");
  l.AddEntry(frame->getObject(1),mcName,"l");
  if(drawSM) l.AddEntry(frame->getObject(2),"SM Higgs","l");
  l.AddEntry(frame->getObject(2+drawSM),"background weighted Data","p");
  if(signal) l.AddEntry(frame->getObject(3+drawSM),"signal weighted Data","p");

  cout << "e" <<endl;

  frame->Draw();
  l.Draw("SAME");
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_%s%s_%s_%s.png",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data(),tag.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/C/CosThetaDist_%s%s_%s_%s.C",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data(),tag.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_%s%s_%s_%s.pdf",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data(),tag.Data()) );
}

void reweight_pt_1d(TH1F* weight_pt, TH1F* weight_pto,TH2F*weight_ptn,TH2F*weight_pt2o,TH2F* weight_ptnr, TH2F* weight_pt2,RooDataSet **dset, RooDataSet *dsetdestination, int numvar){

  if (!(*dset)) return;
///numerator and denominator
  TH1F *hnum = new TH1F("hnum","hnum",n_ptbins_forreweighting,ptbins_forreweighting);
  TH1F *hden = new TH1F("hden","hden",n_ptbins_forreweighting,ptbins_forreweighting);
  hnum->Sumw2();
  hden->Sumw2();

  const char* ptname=Form("roopt%d",numvar);
// RooAbsData->get*() Load a given row of data
 //getRealValue Get value of a RooAbsReal stored in set with given name. If none is found, value of defVal is returned.
  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue(ptname)),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue(ptname)),dsetdestination->store()->weight(i));
  }
  
 //normalize to one 
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());
  hnum->Divide(hden);
  TH1F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_ptrew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float pt = args.getRealValue(ptname);
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(pt)));
    if(debug){
		weight_pt->Fill(neww);	
		weight_pto->Fill(oldw);
		weight_pt2o->Fill(h->FindBin(fabs(pt)),oldw);
		weight_ptn->Fill(h->FindBin(fabs(pt)),neww);	
		if(oldw!=0 && neww!=0)weight_ptnr->Fill(h->FindBin(fabs(pt)),oldw/neww);
		else {weight_ptnr->Fill(-10,1);}
		if(oldw!=0 && neww!=0)weight_pt2->Fill(pt,oldw/neww);
		else {weight_pt2->Fill(-10,1);}
     }
	 newdset->add(args,neww);
  }

  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;

  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Pt 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;
  delete old_dset;
};

//2d reweighting of rho and its sigma
void reweight_rhosigma(TH1F* weight_rho, TH1F* weight_rhoo,TH2F*weight_rhon,TH2F*weight_rho2o,TH2F* weight_rhonr, TH2F* weight_rho2,TH2F*weight_sigman,TH2F*weight_sigma2o,TH2F* weight_sigmanr, TH2F* weight_sigma2,RooDataSet **dset, RooDataSet *dsetdestination, bool deleteold){
  if (!(*dset)) return;
//  TH2F *hnum = new TH2F("hnum","hnum",n_rhobins_forreweighting,rhobins_forreweighting,n_sigmabins_forreweighting,sigmabins_forreweighting);
//  TH2F *hden = new TH2F("hden","hden",n_rhobins_forreweighting,rhobins_forreweighting,n_sigmabins_forreweighting,sigmabins_forreweighting);
  TH2F *hnum = new TH2F("hnum","hnum",100,0,100,20,0,20);
  TH2F *hden = new TH2F("hden","hden",100,0,100,20,0,20);
  hnum->Sumw2();
  hden->Sumw2();
  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue("roorho")),fabs((*dset)->get(i)->getRealValue("roosigma")),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue("roorho")),fabs(dsetdestination->get(i)->getRealValue("roosigma")),dsetdestination->store()->weight(i));
  }

  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());
//data/MC
  hnum->Divide(hden);
  TH2F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_rhosigmarew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float rho = args.getRealValue("roorho");
    float sigma = args.getRealValue("roosigma");
    float neww = oldw*h->GetBinContent(h->FindBin(rho,sigma));
	if(debug){
		weight_rho->Fill(neww);
		weight_rhoo->Fill(oldw);
		weight_rho2o->Fill(h->GetXaxis()->FindBin(rho),oldw);	
		weight_rhon->Fill(h->GetXaxis()->FindBin(rho),neww);	
		if(oldw!=0)weight_rhonr->Fill(h->GetXaxis()->FindBin(rho),oldw/neww);
		else {weight_rhonr->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
		if(oldw!=0)weight_rho2->Fill(rho,oldw/neww);
		weight_sigma2o->Fill(h->GetYaxis()->FindBin(sigma),oldw);	
		weight_sigman->Fill(h->GetYaxis()->FindBin(sigma),neww);	
		if(oldw!=0)weight_sigmanr->Fill(h->GetYaxis()->FindBin(sigma),oldw/neww);
		else {weight_sigmanr->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
		if(oldw!=0)weight_sigma2->Fill(sigma,oldw/neww);
	}
		newdset->add(args,neww);
	  }
	  newdset->SetName((*dset)->GetName());
	  newdset->SetTitle((*dset)->GetTitle());
	  delete hnum; delete hden;
	  RooDataSet *old_dset = *dset;
	  *dset=newdset;
	  std::cout << "RhoSigma2D rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;
	  if (deleteold) delete old_dset;
	};

void MakeSpinPlots::DrawSpinSubTotBackground(TString mcName,bool signal){
  bool drawSM = (smName!="" && smName!=mcName);

  TCanvas cv;
  double thisN  = ws->data(mcName+"_Combined")->sumEntries();
  float norm = thisN;


  if(signal) norm = ws->var(Form("Data_%s_FULLFIT_Nsig",mcName.Data()))->getVal();
  RooPlot *frame = ws->var("cosT")->frame(0,1,10);

  RooDataSet* tmp = (RooDataSet*)ws->data(Form("Data_Combined"))->reduce("(mass>115 && mass<120) || (mass>130 && mass<135)");
  tmp->plotOn(frame,RooFit::Rescale(norm/tmp->sumEntries()));

  ws->pdf(Form("%s_FIT_cosTpdf",mcName.Data()))->plotOn(frame,RooFit::LineColor(kGreen),RooFit::Normalization(norm/tmp->sumEntries()));
  if(drawSM)  ws->pdf(Form("%s_FIT_cosTpdf",smName.Data()))->plotOn(frame,RooFit::LineColor(kRed),RooFit::Normalization(norm/tmp->sumEntries()));
  if(signal){
    RooDataHist *h = (RooDataHist*)ws->data( Form("Data_%s_Combined_bkgSub_cosT",mcName.Data()) );
    h->plotOn(frame,RooFit::MarkerStyle(4));
    std::cout << "Nsig: " << h->sumEntries() << std::endl;
  }

  
  frame->SetMaximum(frame->GetMaximum()*(signal?2.:1.2)*norm/tmp->sumEntries());
  frame->SetMinimum(-1*frame->GetMaximum());
  TLegend l(0.6,0.2,0.95,0.45);
  l.SetFillColor(0);
  l.SetBorderSize(0);
  l.SetHeader("Combined");
  l.AddEntry(frame->getObject(0),"Data m#in [115,120]#cup[130,135]","p");
  l.AddEntry(frame->getObject(1),mcName,"l");
  if(drawSM) l.AddEntry(frame->getObject(2),"SM Higgs","l");
  if(signal) l.AddEntry(frame->getObject(2+drawSM),"bkg-subtracted Data","p");
  
  frame->Draw();
  l.Draw("SAME");
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_SimpleSub_%s%s_%s.png",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/C/CosThetaDist_SimpleSub_%s%s_%s.C",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_SimpleSub_%s%s_%s.pdf",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data()) );
}

vector<double> sigEvents(RooWorkspace *work, int m_hyp, int cat, string altSigFileName, int spin=false){

  RooWorkspace *tempWork;
  if (altSigFileName!=""){
    TFile *temp = TFile::Open(altSigFileName.c_str());
    tempWork = (RooWorkspace*)temp->Get("cms_hgg_workspace");
  }
  else {
    tempWork = work;
  }
  vector<double> result;
  RooDataSet *ggh = (RooDataSet*)tempWork->data(Form("sig_ggh_mass_m%d_cat%d",m_hyp,cat));
  RooDataSet *vbf = (RooDataSet*)tempWork->data(Form("sig_vbf_mass_m%d_cat%d",m_hyp,cat));
  RooDataSet *wzh = (RooDataSet*)tempWork->data(Form("sig_wzh_mass_m%d_cat%d",m_hyp,cat));
  RooDataSet *tth = (RooDataSet*)tempWork->data(Form("sig_tth_mass_m%d_cat%d",m_hyp,cat));
  
  double total;
  if (!spin) {
    total = ggh->sumEntries()+vbf->sumEntries()+wzh->sumEntries()+tth->sumEntries();
    result.push_back(total);
    result.push_back(100*ggh->sumEntries()/total);
    result.push_back(100*vbf->sumEntries()/total);
    result.push_back(100*wzh->sumEntries()/total);
    result.push_back(100*tth->sumEntries()/total);
  }
  else {
    ggh = (RooDataSet*)tempWork->data(Form("sig_mass_m%d_cat%d",m_hyp,cat));
    total = ggh->sumEntries();
    result.push_back(total);
    result.push_back(100.);
    result.push_back(0.);
    result.push_back(0.);
    result.push_back(0.);
  }

  return result;
}

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

  cout << "Number of events in data sig: " << data->numEntries() << endl;
  cout << "Number of events in data bkg: " << data_bkg->numEntries() << endl;
  
  // Initialize tree to save toys to 
  TTree* results = new TTree("results","toy results");
  
  double fa3,fa3Error, fa3Pull;
  double sigFrac,sigFracError, sigFracPull;
  double significance;

  results->Branch("fa3",&fa3,"fa3/D");
  results->Branch("fa3Error",&fa3Error,"fa3Error/D");
  results->Branch("fa3Pull",&fa3Pull,"fa3Pull/D");
  results->Branch("sigFrac",&sigFrac,"sigFrac/D");
  results->Branch("sigFracError",&sigFracError,"sigFracError/D");
  results->Branch("sigFracPull",&sigFracPull,"sigFracPull/D");
  results->Branch("significance",&significance,"significance/D");

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

  RooDataSet* toyData;
  RooDataSet* toyData_bkgOnly;
  int embedTracker=nEvts*counter;
  int embedTracker_bkg=TMath::Ceil(nEvts/3.75*counter);
  RooArgSet *tempEvent;

  RooFitResult *toyfitresults;
  RooFitResult *toyfitresults_sigBkg;
  RooFitResult *toyfitresults_bkgOnly;
  RooRealVar *r_fa3;
  RooRealVar *r_sigFrac;

  for(int i = 0 ; i<nToys ; i++){
    cout <<i<<"<-----------------------------"<<endl;
    //if(toyData) delete toyData;
    toyData = new RooDataSet("toyData","toyData",RooArgSet(*kd));
    toyData_bkgOnly = new RooDataSet("toyData_bkgOnly","toyData_bkgOnly",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;
    }
    if(nEvts+embedTracker_bkg > data_bkg->sumEntries()){
      cout << "Playground::generate() - ERROR!!! Playground::data does not have enough events to fill toy!!!!  bye :) " << endl;
      toyData = NULL;
      abort();
      return 0;
    }

    for(int iEvent=0; iEvent<nEvts; iEvent++){
      if(iEvent==1)
	cout << "generating event: " << iEvent << " embedTracker: " << embedTracker << endl;
      tempEvent = (RooArgSet*) data->get(embedTracker);
      toyData->add(*tempEvent);
      embedTracker++;
    }
    if(bkg){
      for(int iEvent=0; iEvent<nEvts/3.75; iEvent++){
	if(iEvent==1)
	  cout << "generating bkg event: " << iEvent << " embedTracker bkg: " << embedTracker_bkg << endl;
	tempEvent = (RooArgSet*) data_bkg->get(embedTracker_bkg);
	toyData->add(*tempEvent);
	toyData_bkgOnly->add(*tempEvent);
	embedTracker_bkg++;
      }
    }

    if(bkg)
      toyfitresults =model.fitTo(*toyData,Save());
    else
      toyfitresults =modelSignal.fitTo(*toyData,Save());

    //cout<<toyfitresults<<endl;
    r_fa3 = (RooRealVar *) toyfitresults->floatParsFinal().find("fa3");

    fa3 = r_fa3->getVal();
    fa3Error = r_fa3->getError();
    fa3Pull = (r_fa3->getVal() - fa3Val) / r_fa3->getError();
    if(sigFloating){
      r_sigFrac = (RooRealVar *) toyfitresults->floatParsFinal().find("BoverTOT");
      sigFrac = 1-r_sigFrac->getVal();
      sigFracError = r_sigFrac->getError();
      sigFracPull = (1-r_sigFrac->getVal() - sigFracVal) / r_sigFrac->getError();
    }
    // fill TTree
    results->Fill();
  }

  char nEvtsString[100];
  sprintf(nEvtsString,"_%iEvts_%iiter",nEvts, counter);

  // write tree to output file (ouputFileName set at top)
  TFile *outputFile = new TFile("embeddedToys1DKD_fa3Corr_WithBackgDetEffects_"+sampleName[mySample]+nEvtsString+".root","RECREATE");
  results->Write();
  outputFile->Close();

}

int main(int argc, char *argv[]){

	OptionParser(argc,argv);
	

	RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
	RooMsgService::instance().setSilentMode(true);


  
	system(Form("mkdir -p %s",outdir_.c_str()));

	vector<string> procs;
	split(infilenames_,infilenamesStr_,boost::is_any_of(","));
  
   TPython::Exec("import os,imp,re");
    const char * env = gSystem->Getenv("CMSSW_BASE") ; 
      std::string globeRt = env;
      TPython::Exec(Form("buildSMHiggsSignalXSBR = imp.load_source('*', '%s/src/flashggFinalFit/Signal/python/buildSMHiggsSignalXSBR.py')",globeRt.c_str()));
      TPython::Eval(Form("buildSMHiggsSignalXSBR.Init%dTeV()", 13));
   for (unsigned int i =0 ; i<infilenames_.size() ; i++){
	    int mH  =(int) TPython::Eval(Form("int(re.search('_M(.+?)_','%s').group(1))",infilenames_[i].c_str())); 
	   double WH_XS  =  (double)TPython::Eval(Form("buildSMHiggsSignalXSBR.getXS(%d,'%s')",mH,"WH"));
	   double ZH_XS  =  (double)TPython::Eval(Form("buildSMHiggsSignalXSBR.getXS(%d,'%s')",mH,"ZH"));
     float tot_XS = WH_XS + ZH_XS;
     float wFrac=  WH_XS /tot_XS ;
     float zFrac=  ZH_XS /tot_XS ;
      std::cout << "mass "<< mH << " wh fraction "<< WH_XS /tot_XS << ", zh fraction "<< ZH_XS /tot_XS <<std::endl; 
     TFile *infile =  TFile::Open(infilenames_[i].c_str());
	   string outname  =(string) TPython::Eval(Form("'%s'.split(\"/\")[-1].replace(\"VH\",\"WH_VH\")",infilenames_[i].c_str())); 
     TFile *outfile = TFile::Open(outname.c_str(),"RECREATE") ;
    TDirectory* saveDir = outfile->mkdir("tagsDumper");
    saveDir->cd();

    RooWorkspace *inWS = (RooWorkspace*) infile->Get("tagsDumper/cms_hgg_13TeV");
    RooRealVar *intLumi = (RooRealVar*)inWS->var("IntLumi");
    RooWorkspace *outWS = new RooWorkspace("cms_hgg_13TeV");
    outWS->import(*intLumi);
    std::list<RooAbsData*> data =  (inWS->allData()) ;
    std::cout <<" [INFO] Reading WS dataset contents: "<< std::endl;
        for (std::list<RooAbsData*>::const_iterator iterator = data.begin(), end = data.end(); iterator != end; ++iterator )  {
              RooDataSet *dataset = dynamic_cast<RooDataSet *>( *iterator );
              if (dataset) {
	            string zhname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"zh\")",dataset->GetName())); 
	            string whname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"wh\")",dataset->GetName())); 
              RooDataSet *datasetZH = (RooDataSet*) dataset->emptyClone(zhname.c_str(),zhname.c_str());
              RooDataSet *datasetWH = (RooDataSet*) dataset->emptyClone(whname.c_str(),whname.c_str());
                TRandom3 r;
                r.Rndm();
                double x[dataset->numEntries()];
                r.RndmArray(dataset->numEntries(),x);
                for (int j =0; j < dataset->numEntries() ; j++){
                    
                    if( x[j] < wFrac){
                    dataset->get(j);
                    datasetWH->add(*(dataset->get(j)),dataset->weight());
                    } else{
                    dataset->get(j);
                    datasetZH->add(*(dataset->get(j)),dataset->weight());
                    }
                }
              float w =datasetWH->sumEntries();
              float z =datasetZH->sumEntries();
           if(verbose_){
              std::cout << "Original dataset " << *dataset <<std::endl;
              std::cout << "WH       dataset " << *datasetWH <<std::endl;
              std::cout << "ZH       dataset " << *datasetZH <<std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "WH fraction (obs) : WH " << w/(w+z) <<",   ZH "<< z/(w+z) << std::endl;
              std::cout << "WH fraction (exp) : WH " << wFrac <<",   ZH "<< zFrac << std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "********************************************" <<std::endl;
              }
               outWS->import(*datasetWH);
               outWS->import(*datasetZH);
                }
             RooDataHist *datahist = dynamic_cast<RooDataHist *>( *iterator );

              if (datahist) {
	            string zhname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"zh\")",datahist->GetName())); 
	            string whname  =(string) TPython::Eval(Form("'%s'.replace(\"wzh\",\"wh\")",datahist->GetName())); 
              RooDataHist *datahistZH = (RooDataHist*) datahist->emptyClone(zhname.c_str(),zhname.c_str());
              RooDataHist *datahistWH = (RooDataHist*) datahist->emptyClone(whname.c_str(),whname.c_str());
                TRandom3 r;
                r.Rndm();
                double x[datahist->numEntries()];
                r.RndmArray(datahist->numEntries(),x);
                for (int j =0; j < datahist->numEntries() ; j++){
                    
                    datahistWH->add(*(datahist->get(j)),datahist->weight()*wFrac);
                    datahistZH->add(*(datahist->get(j)),datahist->weight()*zFrac);
                }
              float w =datahistWH->sumEntries();
              float z =datahistZH->sumEntries();
           if(verbose_){
              std::cout << "Original datahist " << *datahist <<std::endl;
              std::cout << "WH       datahist " << *datahistWH <<std::endl;
//.........这里部分代码省略.........

void FitterUtils::prepare_PDFs(string trigStr, string BDTVar, double BDTcut,
      string signalfile, string partrecofile, string combfile, string JpsiLeakfile,
      double minBMass, double maxBMass,
      string signaltree, string partrecotree, string combtree, string JpsiLeaktree)
{


   //***********Get the datasets
   TFile* fSignal = new TFile(signalfile.c_str());
   TTree* tSignal = (TTree*)fSignal->Get(signaltree.c_str());
   TFile* fPartReco = new TFile(partrecofile.c_str());
   TTree* tPartReco = (TTree*)fPartReco->Get(partrecotree.c_str());
   TFile* fComb = new TFile(combfile.c_str());
   TTree* tComb = (TTree*)fComb->Get(combtree.c_str()); 
   TFile* fJpsiLeak = new TFile(JpsiLeakfile.c_str());
   TTree* tJpsiLeak = (TTree*)fJpsiLeak->Get(JpsiLeaktree.c_str()); 


   //**********Define variables
   RooRealVar trigVar(trigStr.c_str(), trigStr.c_str(), -10, 10);
   RooRealVar BDTRooRealVar(BDTVar.c_str(), BDTVar.c_str(), -1,1);
   RooRealVar B_plus_M("B_plus_M", "M_{visible}", minBMass, maxBMass, "MeV");
   RooRealVar misPT("misPT", "p_{#perp}", 0, 5000, "MeV");
   RooRealVar B_plus_DTFM_M_zero("B_plus_DTFM_M_zero", "M_{constr}", 0, 20000, "MeV"); 
   RooRealVar e_plus_BremMultiplicity("e_plus_BremMultiplicity","e_plus_BremMultiplicity", -1,2);
   RooRealVar e_minus_BremMultiplicity("e_minus_BremMultiplicity","e_minus_BremMultiplicity", -1,2);

   RooRealVar weightPartReco("weightPartReco", "weightPartReco", 0, 10);
   RooRealVar weightLeakage("weightLeakage", "weightLeakage", 0, 10);
   RooRealVar dataMCWeightee("DataMCWeightee", "DataMCWeightee",0, 30);


   //***********Set only variables needed

   tSignal->SetBranchStatus("*", 0); tSignal->SetBranchStatus("B_plus_M", 1); tSignal->SetBranchStatus("misPT", 1); tSignal->SetBranchStatus("B_plus_DTFM_M_zero", 1); tSignal->SetBranchStatus(BDTVar.c_str(),1);
   tSignal->SetBranchStatus("e_plus_BremMultiplicity", 1); tSignal->SetBranchStatus("e_minus_BremMultiplicity", 1); tSignal->SetBranchStatus(trigStr.c_str()); tSignal->SetBranchStatus("DataMCWeightee",1);

   tPartReco->SetBranchStatus("*", 0); tPartReco->SetBranchStatus("B_plus_M", 1); tPartReco->SetBranchStatus("misPT", 1); tPartReco->SetBranchStatus("B_plus_DTFM_M_zero", 1);tPartReco->SetBranchStatus(BDTVar.c_str(),1);
   tPartReco->SetBranchStatus("e_plus_BremMultiplicity", 1); tPartReco->SetBranchStatus("e_minus_BremMultiplicity", 1); tPartReco->SetBranchStatus(trigStr.c_str()); tPartReco->SetBranchStatus("weightPartReco",1);

   tComb->SetBranchStatus("*", 0); tComb->SetBranchStatus("B_plus_M", 1); tComb->SetBranchStatus("misPT", 1); tComb->SetBranchStatus("B_plus_DTFM_M_zero", 1);tComb->SetBranchStatus(BDTVar.c_str(),1);
   tComb->SetBranchStatus("e_plus_BremMultiplicity", 1); tComb->SetBranchStatus("e_minus_BremMultiplicity", 1); tComb->SetBranchStatus(trigStr.c_str());

   tJpsiLeak->SetBranchStatus("*", 0); tJpsiLeak->SetBranchStatus("B_plus_M", 1); tJpsiLeak->SetBranchStatus("misPT", 1); 
   tJpsiLeak->SetBranchStatus("B_plus_DTFM_M_zero", 1);tJpsiLeak->SetBranchStatus(BDTVar.c_str(),1);
   tJpsiLeak->SetBranchStatus("e_plus_BremMultiplicity", 1); tJpsiLeak->SetBranchStatus("e_minus_BremMultiplicity", 1); tJpsiLeak->SetBranchStatus(trigStr.c_str()); 
   tJpsiLeak->SetBranchStatus("weightLeakage",1);


   //***********Set Binning


   RooBinning defaultMBins(floor((maxBMass-minBMass)/(40.)), B_plus_M.getMin(), B_plus_M.getMax() ); 
   RooBinning defaultMisPTBins(floor(40), misPT.getMin(), misPT.getMax()); 
   RooBinning broaderMBins(floor((maxBMass-minBMass)/(80.)), B_plus_M.getMin(), B_plus_M.getMax()); 
   RooBinning broaderMisPTBins(floor(40), misPT.getMin(), misPT.getMax()); 

   B_plus_M.setBinning( defaultMBins);
   misPT.setBinning( defaultMisPTBins );
   B_plus_M.setBinning( broaderMBins, "broaderBins");
   misPT.setBinning( broaderMisPTBins, "broaderBins" );

   B_plus_DTFM_M_zero.setBins(100);

   RooArgSet argset(BDTRooRealVar, B_plus_DTFM_M_zero, misPT,  B_plus_M, trigVar, e_plus_BremMultiplicity, e_minus_BremMultiplicity);
   RooArgSet argsetPartReco(BDTRooRealVar, B_plus_DTFM_M_zero, misPT,  B_plus_M, trigVar, e_plus_BremMultiplicity, e_minus_BremMultiplicity, weightPartReco);
   RooArgSet argsetLeakage(BDTRooRealVar, B_plus_DTFM_M_zero, misPT,  B_plus_M, trigVar, e_plus_BremMultiplicity, e_minus_BremMultiplicity, weightLeakage);
   RooArgSet argsetSignal(BDTRooRealVar, B_plus_DTFM_M_zero, misPT,  B_plus_M, trigVar, e_plus_BremMultiplicity, e_minus_BremMultiplicity, dataMCWeightee);

   cout<<"getting the datasets:"<<endl;

   RooDataSet* dataSetSignalZeroGamma;
   RooDataSet* dataSetSignalOneGamma;
   RooDataSet* dataSetSignalTwoGamma;
   RooDataSet* dataSetPartReco;
   RooDataSet* dataSetJpsiLeak;
   RooDataSet* dataSetComb;

   TFile* fw(NULL);

   string BDTCutString( ("("+BDTVar+">"+d2s(BDTcut)+")").c_str()  );

   dataSetSignalZeroGamma = new RooDataSet("dataSetSignalZeroGamma", "dataSetSignalZeroGamma", argsetSignal, Import(*tSignal), Cut(( " ("+trigStr+"  > 0.9) && "+BDTCutString+" && ((e_plus_BremMultiplicity+e_minus_BremMultiplicity) > -0.5) && ((e_plus_BremMultiplicity+e_minus_BremMultiplicity) < 0.5) && B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str()), WeightVar("DataMCWeightee")  );

   dataSetSignalOneGamma = new RooDataSet("dataSetSignalOneGamma", "dataSetSignalOneGamma", argsetSignal, Import(*tSignal), Cut(( " ("+trigStr+"  > 0.9) && "+BDTCutString+" && ((e_plus_BremMultiplicity+e_minus_BremMultiplicity) > 0.5) && ((e_plus_BremMultiplicity+e_minus_BremMultiplicity) < 1.5) && B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str()), WeightVar("DataMCWeightee")  );

   dataSetSignalTwoGamma = new RooDataSet("dataSetSignalTwoGamma", "dataSetSignalTwoGamma", argsetSignal, Import(*tSignal), Cut(( " ("+trigStr+"  > 0.9) && "+BDTCutString+" && ((e_plus_BremMultiplicity+e_minus_BremMultiplicity) > 1.5) && ((e_plus_BremMultiplicity+e_minus_BremMultiplicity) < 2.5) && B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str()), WeightVar("DataMCWeightee")  );


   dataSetPartReco = new RooDataSet("dataSetPartReco", "dataSetPartReco",  argsetPartReco, Import(*tPartReco),Cut(("("+trigStr+"  > 0.9) && "+BDTCutString+ " && B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str()), WeightVar("weightPartReco"));

   dataSetJpsiLeak = new RooDataSet("dataSetJpsiLeak", "dataSetJpsiLeak",  argsetLeakage, Import(*tJpsiLeak),Cut(("B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str()), WeightVar("weightLeakage"));

    // dataSetComb = new RooDataSet("dataSetComb", "dataSetComb", tComb, argset, ("("+trigStr+"  > 0.9) && (UBDT3R > "+d2s(BDTcut-0.03)+")  && B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str());
   dataSetComb = new RooDataSet("dataSetComb", "dataSetComb", tComb, argset, ("("+trigStr+"  > 0.9) && "+BDTCutString+"  && B_plus_M > "+d2s(minBMass)+" && B_plus_M < "+d2s(maxBMass)).c_str());

   cout<<"Number of zero: "<< dataSetSignalZeroGamma->sumEntries()<<endl;
   cout<<"Number of one: "<< dataSetSignalOneGamma->sumEntries()<<endl;
   cout<<"Number of two: "<< dataSetSignalTwoGamma->sumEntries()<<endl;
   cout<<"Number of PartReco: "<< dataSetPartReco->sumEntries()<<endl;
//.........这里部分代码省略.........

void angularDistributions_spin0(int plotIndex=0, int binning=80){

  gROOT->ProcessLine(".L  ../PDFs/RooXZsZs_5D.cxx+");
  gROOT->ProcessLine(".L  ../PDFs/RooSpinOne_7D.cxx+");  
  gROOT->ProcessLine(".L  ../PDFs/RooSpinTwo_7D.cxx+");  
  gROOT->ProcessLine(".L  ../src/AngularPdfFactory.cc+");
  gROOT->ProcessLine(".L  ../src/ScalarPdfFactory.cc+");
  gROOT->ProcessLine(".L  ../src/VectorPdfFactory.cc+");
  gROOT->ProcessLine(".L  ../src/TensorPdfFactory.cc+");

  gROOT->ProcessLine(".L  ~/tdrstyle.C");
  setTDRStyle();

  // observables
  RooRealVar* z1mass = new RooRealVar("z1mass","m_{1} [GeV]",40,110);
  RooRealVar* z2mass = new RooRealVar("z2mass","m_{2} [GeV]",1e-09,65);
  RooRealVar* hs = new RooRealVar("costhetastar","cos#theta*",-1,1);
  RooRealVar* h1 = new RooRealVar("costheta1","cos#theta_{1}",-1,1);
  RooRealVar* h2 = new RooRealVar("costheta2","cos#theta_{2}",-1,1);
  RooRealVar* Phi = new RooRealVar("phi","#Phi",-TMath::Pi(),TMath::Pi());
  RooRealVar* Phi1 = new RooRealVar("phistar1","#Phi_{1}",-TMath::Pi(),TMath::Pi());

  vector<RooRealVar*> measureables;
  measureables.push_back(z1mass);
  measureables.push_back(z2mass);
  measureables.push_back(hs);
  measureables.push_back(h1);
  measureables.push_back(h2);
  measureables.push_back(Phi);
  measureables.push_back(Phi1);

  RooRealVar* mzz = new RooRealVar("mzz","mzz",125,100,1000);

  ScalarPdfFactory* SMHiggs = new ScalarPdfFactory(z1mass,z2mass,h1,h2,Phi,mzz);
  SMHiggs->makeSMHiggs();
  SMHiggs->makeParamsConst(true);

  ScalarPdfFactory* PSHiggs = new ScalarPdfFactory(z1mass,z2mass,h1,h2,Phi,mzz);
  PSHiggs->makePSHiggs();
  PSHiggs->makeParamsConst(true);

  ScalarPdfFactory* LGHiggs = new ScalarPdfFactory(z1mass,z2mass,h1,h2,Phi,mzz);
  LGHiggs->makeLGHiggs();
  LGHiggs->makeParamsConst(true);

  ScalarPdfFactory* A2Higgs = new ScalarPdfFactory(z1mass,z2mass,h1,h2,Phi,mzz);
  A2Higgs->makeCustom(0.0,1.0,0.0,0.0,0.0,0.0);
  A2Higgs->makeParamsConst(true);

  TChain* treeSM = new TChain("angles");
  treeSM->Add("/scratch0/hep/whitbeck/OLDHOME/4lHelicity/generatorJHU_V02-01-00/SMHiggs_store/SMHiggs_125GeV.root");
  RooDataSet* dataSM = new RooDataSet("dataSM","dataSM",treeSM,RooArgSet(*z1mass,*z2mass,*hs,*h1,*h2,*Phi,*Phi1));

  TChain* treePS = new TChain("angles");
  treePS->Add("/scratch0/hep/whitbeck/OLDHOME/4lHelicity/generatorJHU_V02-01-00/psScalar_store/psScalar_125GeV.root");
  RooDataSet* data = new RooDataSet("dataPS","dataPS",treePS,RooArgSet(*z1mass,*z2mass,*hs,*h1,*h2,*Phi,*Phi1));

  TChain* treeLG = new TChain("angles");
  treeLG->Add("/scratch0/hep/whitbeck/OLDHOME/4lHelicity/generatorJHU_V02-01-00/0hPlus_store/0hPlus_125GeV.root");
  RooDataSet* dataLG = new RooDataSet("dataLG","dataLG",treeLG,RooArgSet(*z1mass,*z2mass,*hs,*h1,*h2,*Phi,*Phi1));

  double normFix;
  if( ! strcmp( measureables[plotIndex]->GetName() , "phistar1" ) ) normFix=.001/(2.*3.1415);
  else if( ! strcmp( measureables[plotIndex]->GetName() , "costhetastar" ) ) normFix=.001/2.0;
  else normFix=.001;
  
  RooPlot* plot = measureables[plotIndex]->frame(binning);
  plot->GetXaxis()->CenterTitle();
  plot->GetYaxis()->CenterTitle();
  plot->GetYaxis()->SetTitle(" ");
  plot->GetXaxis()->SetNdivisions(-505);

  cout << 1./dataSM->sumEntries() << endl;
  double rescale=(double)(1./dataSM->sumEntries());
  dataSM->plotOn(plot,MarkerColor(kRed),MarkerStyle(4),MarkerSize(1.5),LineWidth(0),XErrorSize(0),Rescale(.001),DataError(RooAbsData::None));
  SMHiggs.PDF->plotOn(plot,LineColor(kRed),LineWidth(2),Normalization( normFix ));
  cout << 1./dataPS->sumEntries() << endl;
  rescale=(double)(1./dataPS->sumEntries());
  dataPS->plotOn(plot,MarkerColor(kBlue),MarkerStyle(27),MarkerSize(1.9),XErrorSize(0),Rescale(.001),DataError(RooAbsData::None));
  PSHiggs.PDF->plotOn(plot,LineColor(kBlue),LineWidth(2),Normalization( normFix ));
  cout << 1./dataLG->sumEntries() << endl;
  rescale=(double)(1./dataLG->sumEntries());
  dataLG->plotOn(plot,MarkerColor(kGreen+3),MarkerStyle(25),MarkerSize(1.5),XErrorSize(0),Rescale(.001),DataError(RooAbsData::None));
  LGHiggs.PDF->plotOn(plot,LineColor(kGreen+3),LineWidth(2),Normalization( normFix ));

  TGaxis::SetMaxDigits(3);

  TCanvas* can =new TCanvas("can","can",600,600);

  gStyle->SetPadLeftMargin(0.05);

  char temp[150];
  sprintf(temp,"%s>>SM_histo(%i,%i,%i)",measureables[plotIndex]->GetName(),binning,(int)measureables[plotIndex]->getMin(),(int)measureables[plotIndex]->getMin());
  treeSM->Draw(temp);
  TH1F* SM_histo = (TH1F*) gDirectory->Get("SM_histo");
  sprintf(temp,"%s>>PS_histo(%i,%i,%i)",measureables[plotIndex]->GetName(),binning,(int)measureables[plotIndex]->getMin(),(int)measureables[plotIndex]->getMin());
  treePS->Draw(temp);
  TH1F* PS_histo = (TH1F*) gDirectory->Get("PS_histo");
  sprintf(temp,"%s>>LG_histo(%i,%i,%i)",measureables[plotIndex]->GetName(),binning,(int)measureables[plotIndex]->getMin(),(int)measureables[plotIndex]->getMin());
  treeLG->Draw(temp);
//.........这里部分代码省略.........

void forData(string channel, string catcut, bool removeMinor=true){

  // Suppress all the INFO message

  RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);

  // Input files and sum all backgrounds

  TChain* treeData  = new TChain("tree");
  TChain* treeZjets = new TChain("tree");

  if( channel == "ele" ){

    treeData->Add(Form("%s/data/SingleElectron-Run2015D-05Oct2015-v1_toyMCnew.root",  channel.data()));
    treeData->Add(Form("%s/data/SingleElectron-Run2015D-PromptReco-V4_toyMCnew.root", channel.data()));

  }

  else if( channel == "mu" ){

    treeData->Add(Form("%s/data/SingleMuon-Run2015D-05Oct2015-v1_toyMCnew.root",  channel.data()));
    treeData->Add(Form("%s/data/SingleMuon-Run2015D-PromptReco-V4_toyMCnew.root", channel.data()));

  }

  else return;

  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-100to200_13TeV_toyMCnew.root", channel.data()));
  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-200to400_13TeV_toyMCnew.root", channel.data()));
  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-400to600_13TeV_toyMCnew.root", channel.data()));
  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-600toInf_13TeV_toyMCnew.root", channel.data()));

  // To remove minor background contribution in data set (weight is -1)

  if( removeMinor ){

    treeData->Add(Form("%s/VV/WW_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));
    treeData->Add(Form("%s/VV/WZ_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));
    treeData->Add(Form("%s/VV/ZZ_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));
    treeData->Add(Form("%s/TT/TT_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));

  }

  // Define all the variables from the trees

  RooRealVar cat ("cat", "", 0, 2);
  RooRealVar mJet("prmass", "M_{jet}",  30.,  300., "GeV");
  RooRealVar mZH ("mllbb",   "M_{ZH}", 900., 3000., "GeV");
  RooRealVar evWeight("evweight", "", -1.e3, 1.e3);

  // Set the range in jet mass

  mJet.setRange("allRange", 30., 300.);
  mJet.setRange("lowSB",    30.,  65.);
  mJet.setRange("highSB",  135., 300.);
  mJet.setRange("signal",  105., 135.);

  RooBinning binsmJet(54, 30, 300);

  RooArgSet variables(cat, mJet, mZH, evWeight);

  TCut catCut = Form("cat==%s", catcut.c_str());
  TCut sbCut  = "prmass>30 && !(prmass>65 && prmass<135) && prmass<300";
  TCut sigCut = "prmass>105 && prmass<135";

  // Create a dataset from a tree -> to process an unbinned likelihood fitting

  RooDataSet dataSetData   ("dataSetData",    "dataSetData",    variables, Cut(catCut),           WeightVar(evWeight), Import(*treeData));
  RooDataSet dataSetDataSB ("dataSetDataSB",  "dataSetDataSB",  variables, Cut(catCut && sbCut),  WeightVar(evWeight), Import(*treeData));
  RooDataSet dataSetZjets  ("dataSetZjets",   "dataSetZjets",   variables, Cut(catCut),           WeightVar(evWeight), Import(*treeZjets));
  RooDataSet dataSetZjetsSB("dataSetZjetsSB", "dataSetZjetsSB", variables, Cut(catCut && sbCut),  WeightVar(evWeight), Import(*treeZjets));  
  RooDataSet dataSetZjetsSG("dataSetZjetsSG", "dataSetZjetsSG", variables, Cut(catCut && sigCut), WeightVar(evWeight), Import(*treeZjets));
  
  // Total events number

  float totalMcEv   = dataSetZjetsSB.sumEntries() + dataSetZjetsSG.sumEntries();
  float totalDataEv = dataSetData.sumEntries();

  RooRealVar nMcEvents("nMcEvents", "nMcEvents", 0., 99999.);
  RooRealVar nDataEvents("nDataEvents", "nDataEvents", 0., 99999.);

  nMcEvents.setVal(totalMcEv);
  nMcEvents.setConstant(true);

  nDataEvents.setVal(totalDataEv);
  nDataEvents.setConstant(true);

  // Signal region jet mass

  RooRealVar constant("constant", "constant", -0.02,  -1.,   0.);
  RooRealVar offset  ("offset",   "offset",     30., -50., 200.);
  RooRealVar width   ("width",    "width",     100.,   0., 200.);

  if( catcut == "1" ) offset.setConstant(true);
  
  RooErfExpPdf model_mJet("model_mJet", "model_mJet", mJet, constant, offset, width);
  RooExtendPdf ext_model_mJet("ext_model_mJet", "ext_model_mJet", model_mJet, nMcEvents);

  RooFitResult* mJet_result = ext_model_mJet.fitTo(dataSetZjets, SumW2Error(true), Extended(true), Range("allRange"), Strategy(2), Minimizer("Minuit2"), Save(1));

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

int main(int argc, char *argv[]){

	OptionParser(argc,argv);
	

	RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
	RooMsgService::instance().setSilentMode(true);


  
	system(Form("mkdir -p %s",outdir_.c_str()));

	vector<string> procs;
	split(infilenames_,infilenamesStr_,boost::is_any_of(","));
  
   TPython::Exec("import os,imp,re");
   for (unsigned int i =0 ; i<infilenames_.size() ; i++){
     TFile *infile =  TFile::Open(infilenames_[i].c_str());
	   string outname  =(string) TPython::Eval(Form("'%s'.split(\"/\")[-1].replace('root','_reduced.root')",infilenames_[i].c_str())); 
     TFile *outfile = TFile::Open(outname.c_str(),"RECREATE") ;
    TDirectory* saveDir = outfile->mkdir("tagsDumper");
    saveDir->cd();

    RooWorkspace *inWS = (RooWorkspace*) infile->Get("tagsDumper/cms_hgg_13TeV");
    RooRealVar *intLumi = (RooRealVar*)inWS->var("IntLumi");
    RooWorkspace *outWS = new RooWorkspace("cms_hgg_13TeV");
    outWS->import(*intLumi);
    std::list<RooAbsData*> data =  (inWS->allData()) ;
    std::cout <<" [INFO] Reading WS dataset contents: "<< std::endl;
        for (std::list<RooAbsData*>::const_iterator iterator = data.begin(), end = data.end(); iterator != end; ++iterator )  {
              RooDataSet *dataset = dynamic_cast<RooDataSet *>( *iterator );
              if (dataset) {
              RooDataSet *datasetReduced = (RooDataSet*) dataset->emptyClone(dataset->GetName(),dataset->GetName());
                TRandom3 r;
                r.Rndm();
                double x[dataset->numEntries()];
                r.RndmArray(dataset->numEntries(),x);
                int desiredEntries = floor(0.5+ dataset->numEntries()*fraction_);
                int modFraction = floor(0.5+ 1/fraction_);
                int finalEventCount=0;
                for (int j =0; j < dataset->numEntries() ; j++){
                    if( j%modFraction==0){
                      finalEventCount++;
                    }
                 }
                float average_weight= dataset->sumEntries()/finalEventCount;
                for (int j =0; j < dataset->numEntries() ; j++){
                    if( j%modFraction==0){
                    dataset->get(j);
                    datasetReduced->add(*(dataset->get(j)),average_weight);
                    }
                }
              float entriesIN =dataset->sumEntries();
              float entriesOUT =datasetReduced->sumEntries();
           if(verbose_){
              std::cout << "Original dataset " << *dataset <<std::endl;
              std::cout << "Reduced       dataset " << *datasetReduced <<std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "fraction (obs) : " << entriesOUT/entriesIN << std::endl;
              std::cout << "fraction (exp) : " << fraction_ << std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "********************************************" <<std::endl;
              }
               outWS->import(*datasetReduced);
                }
                
             RooDataHist *datahist = dynamic_cast<RooDataHist *>( *iterator );

              if (datahist) {
              RooDataHist *datahistOUT = (RooDataHist*) datahist->emptyClone(datahist->GetName(),datahist->GetName());
                TRandom3 r;
                r.Rndm();
                for (int j =0; j < datahist->numEntries() ; j++){
                    
                    datahistOUT->add(*(datahist->get(j)),datahist->weight());
                }
              float w =datahistOUT->sumEntries();
              float z =datahist->sumEntries();
           if(verbose_){
              std::cout << "Original datahist " << *datahist <<std::endl;
              std::cout << "Reduced  datahist " << *datahistOUT<<std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "WH fraction (obs) : " << w/(z) <<std::endl;
              std::cout << "WH fraction (exp) : " << fraction_ << std::endl;
              std::cout << "********************************************" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "" <<std::endl;
              std::cout << "********************************************" <<std::endl;
              }
               outWS->import(*datahistOUT);
                }
                  }
   saveDir->cd();
   outWS->Write();
   outfile->Close();
   infile->Close();
//.........这里部分代码省略.........

void Fit3D::plotFitAccuracy(
    const RooDataSet& mc_data,
    const RooFitResult& fit)
{
  fit.Print("v");

  double n_true_bs = mc_data.sumEntries("component == component::bs");
  double n_true_bd = mc_data.sumEntries("component == component::bd");
  double n_true_cw = mc_data.sumEntries("component == component::cw");
  double n_true_ww = mc_data.sumEntries("component == component::ww");
  double n_true_cn = mc_data.sumEntries("component == component::cn");

  RooRealVar* bs_fit = (RooRealVar*) fit.floatParsFinal().find("n_bs_pp");
  RooRealVar* bd_fit = (RooRealVar*) fit.floatParsFinal().find("n_bd_pp");
  RooRealVar* cw_fit = (RooRealVar*) fit.floatParsFinal().find("n_cw_pp");
  RooRealVar* ww_fit = (RooRealVar*) fit.floatParsFinal().find("n_ww_pp");
  RooRealVar* cn_fit = (RooRealVar*) fit.floatParsFinal().find("n_cn_pp");
  
  TString title("Fit Accuracy (N^{++}_{fit}-N^{++}_{true})");
  TString filename(output_path_ + "fit_accuracy_pp");
  if (!bs_fit) {
    bs_fit = (RooRealVar*) fit.floatParsFinal().find("n_bs_nn");
    bd_fit = (RooRealVar*) fit.floatParsFinal().find("n_bd_nn");
    cw_fit = (RooRealVar*) fit.floatParsFinal().find("n_cw_nn");
    ww_fit = (RooRealVar*) fit.floatParsFinal().find("n_ww_nn");
    cn_fit = (RooRealVar*) fit.floatParsFinal().find("n_cn_nn");
    title = TString("Fit Accuracy (N^{--}_{fit}-N^{--}_{true})");
    filename = TString(output_path_ + "fit_accuracy_nn");
  }
  if (!bs_fit) {
    // Error. Quit while ahead.
    cout << "Error in plotFitAccuracy(): "
         << "Cannot find fit variables. Check names are valid."
         << endl;
    return;
  }

  std::cout << n_true_bd << std::endl;
  std::cout << n_true_bs << std::endl;
  std::cout << n_true_cn << std::endl;
  std::cout << n_true_cw << std::endl;
  std::cout << n_true_ww << std::endl;
  
  TCanvas c1("c1", title, 200, 10, 700, 500);
  c1.SetGrid();
  double x[5] = {1, 2, 3, 4, 5};
  double y[5] = {
      bs_fit->getVal() - n_true_bs,
      bd_fit->getVal() - n_true_bd,
      cw_fit->getVal() - n_true_cw,
      ww_fit->getVal() - n_true_ww,
      cn_fit->getVal() - n_true_cn};
  double exl[5] = {0, 0, 0, 0, 0};
  double exh[5] = {0, 0, 0, 0, 0};
  double eyl[5] = {
      -bs_fit->getErrorLo(),
      -bd_fit->getErrorLo(),
      -cw_fit->getErrorLo(),
      -ww_fit->getErrorLo(),
      -cn_fit->getErrorLo()};
  double eyh[5] = {
      bs_fit->getErrorHi(),
      bd_fit->getErrorHi(),
      cw_fit->getErrorHi(),
      ww_fit->getErrorHi(),
      cn_fit->getErrorHi()};
  TGraphAsymmErrors* gr = new TGraphAsymmErrors(5, x, y, exl, exh, eyl, eyh);
  
  TLatex* cc_bs_label = new TLatex(gr->GetX()[0], gr->GetY()[0], " CC (B_{s})");
  TLatex* cc_bd_label = new TLatex(gr->GetX()[1], gr->GetY()[1], " CC (B_{d})");
  TLatex* cw_label = new TLatex(gr->GetX()[2], gr->GetY()[2], " CW");
  TLatex* ww_label = new TLatex(gr->GetX()[3], gr->GetY()[3], " WW");
  TLatex* cn_label = new TLatex(gr->GetX()[4], gr->GetY()[4], " CN");
  
  gr->GetListOfFunctions()->Add(cc_bs_label);
  gr->GetListOfFunctions()->Add(cc_bd_label);
  gr->GetListOfFunctions()->Add(cw_label);
  gr->GetListOfFunctions()->Add(ww_label);
  gr->GetListOfFunctions()->Add(cn_label);
  
  gr->SetTitle(title);
  gr->SetMarkerStyle(kOpenCircle);
  gr->SetMarkerColor(4);
  gr->Draw("AP");

  c1.Print(filename + ".eps");
  
  TFile accuracy_plot_file(filename + ".root", "RECREATE");
  gr->Write();
  accuracy_plot_file.Close();
  
  return;
}