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

//.........这里部分代码省略.........
      if(isSingle11==0)continue;
      MCBu_DTFNoFix_eta_Prime_MF11=MCEta_Mass11[0];
      MC11FlatTree->Fill();
    }
    
    MC11FlatTree->Write();
    MC11FlatOut->Close();
  }
  
  //_____________________________________________LOAD ROODATASETS___________________________________

  TFile* MCFlatInput12= new TFile("MCMinimalFile12.root");
  TTree* MCFlatInputTree12=(TTree*)MCFlatInput12->Get("DecayTree");

  TFile* MCFlatInput11= new TFile("MCMinimalFile11.root");
  TTree* MCFlatInputTree11=(TTree*)MCFlatInput11->Get("DecayTree");

  RooRealVar MCBMass("Bu_DTF_MF","Bu_DTF_MF",5000.0,5600.0);
  RooRealVar MCEtaMass("eta_prime_MM","eta_prime_MM",700.0,1200.0);
  RooRealVar BDT_response("BDT_response","BDT_response",-1.0,1.0);
  RooRealVar gamma_CL("gamma_CL","gamma_CL",0.1,1.0);
  RooArgSet Args(MCBMass,MCEtaMass,BDT_response,gamma_CL);

  RooDataSet* MCData12 = new RooDataSet("MCData12","MCData12",Args,Import(*MCFlatInputTree12));
  
  std::cout <<" Data File 12 Loaded"<<std::endl;
  
  RooDataSet* MCData11 = new RooDataSet("MCData11","MCData11",Args,Import(*MCFlatInputTree11));

  std::cout<<" Data File 11 loaded"<<std::endl;

  RooDataSet* MCDataAll= new RooDataSet("MCDataAll","MCDataAll",Args);

  MCDataAll->append(*MCData12);
  MCDataAll->append(*MCData11);
  
  RooPlot* massFrame = MCBMass.frame(Title("Data Import Check"),Bins(50));
  MCDataAll->plotOn(massFrame);
  
  RooPlot *BDTFrame = BDT_response.frame(Title("BDT Cut Check"),Bins(50));
  MCDataAll->plotOn(BDTFrame);
  TCanvas C;
  C.Divide(2,1);
  C.cd(1);
  massFrame->Draw();
  C.cd(2);
  BDTFrame->Draw();
  C.SaveAs("ImportChecks.eps");

  //________________________________MAKE MCROODATACATEGORIES__________________________________

  RooDataSet* MCBData=(RooDataSet*)MCDataAll->reduce(RooArgSet(MCBMass));
  MCBData->Print("v");
  
  RooDataSet* MCEtaData=(RooDataSet*)MCDataAll->reduce(RooArgSet(MCEtaMass));
  MCEtaData->Print("v");

  RooCategory MCMassType("MCMassType","MCMassType") ;
  MCMassType.defineType("B") ;
  MCMassType.defineType("Eta") ;
  
  // Construct combined dataset in (x,sample)
  RooDataSet MCcombData("MCcombData","MC combined data",Args,Index(MCMassType),Import("B",*MCBData),Import("Eta",*MCEtaData));

  
  //=============================================== MC FIT MODEL===================================

int main (int argc, char **argv)
{
  const char* chInFile = "ws.root";
  const char* chOutFile = "ws_gen.root";
  int numSignal = 10000;
  int numBkg = 100000;

  char option_char;
  while ( (option_char = getopt(argc,argv, "i:o:s:b:")) != EOF )
    switch (option_char)
      {
         case 'i': chInFile = optarg; break;
         case 'o': chOutFile = optarg; break;
         case 's': numSignal = atoi(optarg); break;
         case 'b': numBkg = atoi(optarg); break;
         case '?': fprintf (stderr,
                            "usage: %s [i<input file> o<output file>]\n", argv[0]);
      }

  cout << "In File = " << chInFile << endl;
  cout << "Out File = " << chOutFile << endl;
  cout << "Signal Events = " << numSignal << endl;
  cout << "Bkg Events = " << numBkg << endl;

  TFile inFile(chInFile,"READ");
  RooWorkspace* ws = (RooWorkspace*) inFile.Get("rws");
  TFile outFile(chOutFile,"RECREATE");

/*
  ws->var("tau")->setVal(1.417);
  ws->var("DG")->setVal(0.151);
  ws->var("beta")->setVal(0.25);
  ws->var("A02")->setVal(0.553);
  ws->var("A1")->setVal(0.487);
  ws->var("delta_l")->setVal(3.15);
  ws->var("fs")->setVal(0.147);
*/

//  ws->var("delta_l")->setConstant(kTRUE);
//  ws->var("delta_p")->setConstant(kTRUE);
//  ws->var("Dm")->setConstant(kTRUE);

  //*ws->var("xs") = numSignal/(numSignal+numBkg);
//  int numSignal = numEvents * ws->var("xs")->getVal();
//  int numBkg = numEvents - numSignal;

  ws->factory("Gaussian::dilutionGauss(d,0,0.276)");
  //ws->factory("SUM::dSignalPDF(xds[0.109]*dilutionGauss,TruthModel(d))");
  //ws->factory("SUM::dBkgPDF(xdb[0.109]*dilutionGauss,TruthModel(d))");
  ws->factory("SUM::dSignalPDF(xds[1]*dilutionGauss,TruthModel(d))");
  ws->factory("SUM::dBkgPDF(xdb[1]*dilutionGauss,TruthModel(d))");

/*
  ws->factory("GaussModel::xetGaussianS(et,meanGaussEtS,sigmaGaussEtS)");
  ws->factory("Decay::xerrorSignal(et,tauEtS,xetGaussianS,SingleSided]");

  ws->factory("PROD::xsignalTimeAngle(timeAngle|et,xerrorSignal");
  ws->factory("PROD::xsignal(massSignal,xsignalTimeAngle,DmConstraint)");
*/

  RooDataSet* dSignalData = ws->pdf("dSignalPDF")->generate(RooArgSet(*ws->var("d")),numSignal);
  RooDataSet *dataSignal = ws->pdf("signal")->generate(RooArgSet(*ws->var("m"),*ws->var("t"),*ws->var("et"),*ws->var("cpsi"),*ws->var("ctheta"),*ws->var("phi")), RooFit::ProtoData(*dSignalData));

  ws->factory("GaussModel::xetGaussianPR(et,meanGaussEtPR,sigmaGaussEtPR)");
  ws->factory("Decay::xerrBkgPR(et,tauEtPR,xetGaussianPR,SingleSided]");

  ws->factory("GaussModel::xetGaussianNP(et,meanGaussEtNP,sigmaGaussEtNP)");
  ws->factory("Decay::xerrBkgNP(et,tauEtNP,xetGaussianNP,SingleSided]");


  /* Time */
  ws->factory("GaussModel::xresolution(t,0,scale,et)");
  ws->factory("Decay::xnegativeDecay(t,tauNeg,xresolution,Flipped)");
  ws->factory("Decay::xpositiveDecay(t,tauPos,xresolution,SingleSided)");
  ws->factory("Decay::xpositiveLongDecay(t,tauLngPos,xresolution,SingleSided)");

  ws->factory("RSUM::xtBkgNP(xn*xnegativeDecay,xp*xpositiveDecay,xpositiveLongDecay");

/*               Promt and Non-Prompt                       */
   ws->factory("PROD::xtimeBkgNP(xtBkgNP|et,xerrBkgNP)");
   ws->factory("PROD::xtimeBkgPR(xresolution|et,xerrBkgPR)");

   ws->factory("PROD::xPrompt(massBkgPR,xtimeBkgPR,anglePR)");
   ws->factory("PROD::xNonPrompt(massBkgNP,xtimeBkgNP,angleNP)");

  ws->factory("SUM::xbackground(xprompt*xPrompt,xNonPrompt)");


  RooDataSet* dBkgData = ws->pdf("dBkgPDF")->generate(RooArgSet(*ws->var("d")),numBkg);
  RooDataSet* dataBkg  = ws->pdf("xbackground")->generate(RooArgSet(*ws->var("m"),*ws->var("t"),*ws->var("et"),*ws->var("cpsi"),*ws->var("ctheta"),*ws->var("phi")), numBkg);

  dataBkg->merge(dBkgData);
  dataSignal->SetName("dataGenSignal");
  dataBkg->SetName("dataGenBkg");
  ws->import(*dataSignal);
  ws->import(*dataBkg);

  ////ws->import(*dataBkg,RooFit::Rename("dataGenBkg"));

  dataSignal->append(*dataBkg);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  vector<double> muFabErrHigh;
  vector<double> muPaulErrHigh;
  vector<double> muChi2ErrHigh;
  vector<double> muAICErrHigh;
  muTree->Branch("jobn",&jobn);
  muTree->Branch("toyn",&toyn);
  muTree->Branch("truthModel",&truthModel);
  muTree->Branch("muFab",&muFab);
  muTree->Branch("muPaul",&muPaul);
  muTree->Branch("muChi2",&muChi2);
  muTree->Branch("muAIC",&muAIC);
  muTree->Branch("muFabErrLow",&muFabErrLow);
  muTree->Branch("muPaulErrLow",&muPaulErrLow);
  muTree->Branch("muChi2ErrLow",&muChi2ErrLow);
  muTree->Branch("muAICErrLow",&muAICErrLow);
  muTree->Branch("muFabErrHigh",&muFabErrHigh);
  muTree->Branch("muPaulErrHigh",&muPaulErrHigh);
  muTree->Branch("muChi2ErrHigh",&muChi2ErrHigh);
  muTree->Branch("muAICErrHigh",&muAICErrHigh);
  
  //TH1F *muDistFab = new TH1F("muDistFab","muDistFab",int(20*(mu_high-mu_low)),mu_low,mu_high);
  //TH1F *muDistPaul = new TH1F("muDistPaul","muDistPaul",int(20*(mu_high-mu_low)),mu_low,mu_high);
  //TH1F *muDistChi2 = new TH1F("muDistChi2","muDistChi2",int(20*(mu_high-mu_low)),mu_low,mu_high);
  //TH1F *muDistAIC = new TH1F("muDistAIC","muDistAIC",int(20*(mu_high-mu_low)),mu_low,mu_high);
  
  mass->setBins(320);
  RooDataSet *data = (RooDataSet*)bkgWS->data(Form("data_mass_cat%d",cat));
  //RooDataSet *data = (RooDataSet*)bkgWS->data(Form("data_cat%d_7TeV",cat));
  RooDataHist *dataBinned = new RooDataHist(Form("roohist_data_mass_cat%d",cat),Form("roohist_data_mass_cat%d",cat),RooArgSet(*mass),*data);
  RooDataSet *sigMC = (RooDataSet*)sigWS->data(Form("sig_ggh_mass_m%d_cat%d",expectSignalMass,cat));
  RooDataSet *sigMC_vbf = (RooDataSet*)sigWS->data(Form("sig_wzh_mass_m%d_cat%d",expectSignalMass,cat));
  RooDataSet *sigMC_wzh = (RooDataSet*)sigWS->data(Form("sig_vbf_mass_m%d_cat%d",expectSignalMass,cat));
  RooDataSet *sigMC_tth = (RooDataSet*)sigWS->data(Form("sig_tth_mass_m%d_cat%d",expectSignalMass,cat));
  sigMC->append(*sigMC_vbf);
  sigMC->append(*sigMC_wzh);
  sigMC->append(*sigMC_tth);
  //RooExtendPdf *ggh_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_ggh",cat));
  //RooExtendPdf *vbf_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_vbf",cat));
  //RooExtendPdf *wzh_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_wzh",cat));
  //RooExtendPdf *tth_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_tth",cat));
  //RooAbsPdf *sigPdf = new RooAddPdf(Form("sigpdfsmrel_cat%d_7TeV",cat),Form("sigpdfsmrel_cat%d_7TeV",cat),RooArgList(*ggh_pdf,*vbf_pdf,*wzh_pdf,*tth_pdf));
  
  if (!dataBinned || !sigMC){
    cerr << "ERROR -- one of data or signal is NULL" << endl;
    exit(1);
  }
  
  // set of truth models to throw toys from
  PdfModelBuilder toysModel;
  toysModel.setObsVar(mass);
  toysModel.setSignalModifier(mu);
  // add truth pdfs from config datfile these need to be cached
  // to throw a toy from the SB fit make sure that the cache happens at makeSBPdfs
  for (vector<pair<int,pair<string,string> > >::iterator it=toysMap.begin(); it!=toysMap.end(); it++){
    if (it->first==-1) { // this is a hyrbid toy
      throwHybridToys=true;
      vector<string> temp;
      split(temp,it->second.first,boost::is_any_of(","));
      split(switchFunc,it->second.second,boost::is_any_of(","));
      for (unsigned int i=0; i<temp.size(); i++){
        switchMass.push_back(atof(temp[i].c_str()));
      }
      continue; 
    }
    if (it->first==-2) { // this is a keys pdf toy
      double rho = lexical_cast<double>(it->second.first);

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

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

  RooDataSet *Data = (RooDataSet*)w->data("Data2011")->Clone("Data");
  Data->append( *((RooDataSet*)w->data("Data2012")) );

  RooDataSet *Bs2Kst0Kst0_MC = (RooDataSet*)w->data("Bs2Kst0Kst0_MC2011")->Clone("Bs2KstKst0_MC");
  Bs2Kst0Kst0_MC->append( *((RooDataSet*)w->data("Bs2Kst0Kst0_MC2012")) );

  RooDataSet *Bs2Kst0Kst01430_MC = (RooDataSet*)w->data("Bs2Kst0Kst01430_MC2011")->Clone("Bs2KstKst0_MC");
  Bs2Kst0Kst01430_MC->append( *((RooDataSet*)w->data("Bs2Kst0Kst01430_MC2012")) );

  RooDataSet *Bs2Kst01430Kst01430_MC = (RooDataSet*)w->data("Bs2Kst01430Kst01430_MC2011")->Clone("Bs2KstKst0_MC");
  Bs2Kst01430Kst01430_MC->append( *((RooDataSet*)w->data("Bs2Kst01430Kst01430_MC2012")) );

  RooDataSet *Bd2Kst0Kst0_MC = (RooDataSet*)w->data("Bd2Kst0Kst0_MC2011")->Clone("Bs2KstKst0_MC");
  Bd2Kst0Kst0_MC->append( *((RooDataSet*)w->data("Bd2Kst0Kst0_MC2012")) );

  RooDataSet *Bd2PhiKst0_MC = (RooDataSet*)w->data("Bd2PhiKst0_MC2011")->Clone("Bs2KstKst0_MC");
  Bd2PhiKst0_MC->append( *((RooDataSet*)w->data("Bd2PhiKst0_MC2012")) );

  RooDataSet *Bs2PhiKst0_MC = (RooDataSet*)w->data("Bs2PhiKst0_MC2011")->Clone("Bs2KstKst0_MC");
  Bs2PhiKst0_MC->append( *((RooDataSet*)w->data("Bs2PhiKst0_MC2012")) );

  RooDataSet *Bd2RhoKst0_MC = (RooDataSet*)w->data("Bd2RhoKst0_MC2011")->Clone("Bs2KstKst0_MC");
  Bd2RhoKst0_MC->append( *((RooDataSet*)w->data("Bd2RhoKst0_MC2012")) );

  RooDataSet *Lb2ppipipi_MC = (RooDataSet*)w->data("Lb2ppipipi_MC2011")->Clone("Bs2KstKst0_MC");
  Lb2ppipipi_MC->append( *((RooDataSet*)w->data("Lb2ppipipi_MC2012")) );

  RooDataSet *Lb2pKpipi_MC = (RooDataSet*)w->data("Lb2pKpipi_MC2011")->Clone("Bs2KstKst0_MC");
  Lb2pKpipi_MC->append( *((RooDataSet*)w->data("Lb2pKpipi_MC2012")) );


  w->import(*Data);
  w->import(*Bs2Kst0Kst0_MC);
  w->import(*Bs2Kst0Kst01430_MC);
  w->import(*Bs2Kst01430Kst01430_MC);
  w->import(*Bd2Kst0Kst0_MC);
  w->import(*Bd2PhiKst0_MC);
  w->import(*Bs2PhiKst0_MC);
  w->import(*Bd2RhoKst0_MC);
  w->import(*Lb2ppipipi_MC);
  w->import(*Lb2pKpipi_MC);

  RooRealVar *mass = (RooRealVar*)w->var("B_s0_DTF_B_s0_M");

  fitIpatia( w, "bs2kstkst_mc", "Bs2KstKst0_MC");

  // Make the PDF here
  RooRealVar *p1 = new RooRealVar("p1","p1",-0.002,-0.004,0.);
  RooExponential *exp = new RooExponential("exp","exp",*mass,*p1);
  //RooRealVar *m1 = new RooRealVar("m1","m1",5320,5380);
  //RooRealVar *s1 = new RooRealVar("s1","s1",1,20);
  //RooGaussian *sig = new RooGaussian("sig","sig",*mass,*m1,*s1);
  RooRealVar *m2 = new RooRealVar("m2","m2",5320,5380);
  RooRealVar *s2 = new RooRealVar("s2","s2",1,20);
  RooGaussian *sig_bd = new RooGaussian("sig_bd","sig_bd",*mass,*m2,*s2);

  //
  RooRealVar *bs2kstkst_l       = new RooRealVar( "bs2kstkst_l"    ,"", -5, -20, -1.);
  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   ,20 );
  RooRealVar *bs2kstkst_mu      = new RooRealVar( "bs2kstkst_mu"   ,"",5350  ,5380     );
  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","sig",*mass,*bs2kstkst_l,*bs2kstkst_zeta,*bs2kstkst_fb,*bs2kstkst_sigma,*bs2kstkst_mu,*bs2kstkst_a,*bs2kstkst_n,*bs2kstkst_a2,*bs2kstkst_n2);

  RooRealVar *bkg_y = new RooRealVar("bkg_y","bkg_y",10e3,10e5);
  RooRealVar *sig_y = new RooRealVar("sig_y","sig_y",0,20e3);
  RooRealVar *sig_bd_y = new RooRealVar("sig_bd_y","sig_bd_y",0,3000);

  RooArgList *pdfs = new RooArgList();
  RooArgList *yields = new RooArgList();

  pdfs->add( *exp );
  pdfs->add( *sig );
  pdfs->add( *sig_bd );

  yields->add( *bkg_y );
  yields->add( *sig_y );
  yields->add( *sig_bd_y );

  RooAddPdf *pdf = new RooAddPdf("pdf","pdf",*pdfs,*yields);

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

  RooPlot *plot = mass->frame();
  Data->plotOn(plot);
    // set fit params constant;
  pdf->plotOn(plot);

  TCanvas *c = new TCanvas();
  plot->Draw();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // define output works
  RooWorkspace *wOut = new RooWorkspace("wbkg","wbkg") ;
  
  // util;ities for the combined fit
  RooCategory     finalcat  ("finalcat",  "finalcat") ;  
  RooSimultaneous fullbkgpdf("fullbkgpdf","fullbkgpdf",finalcat);
  RooDataSet      datacomb  ("datacomb",  "datacomb",  RooArgList(*hmass,finalcat)) ;

  RooDataSet *datacombcat = new RooDataSet("data_combcat","",RooArgList(*hmass)) ;
  
  // add the 'combcat' to the list...if more than one cat
  if( catnames.size() > 1 ) {
    catnames.push_back("combcat");    
    catdesc.push_back("Combined");
  }
  
  for (UInt_t icat=0; icat<catnames.size(); ++icat) {
    TString catname = catnames.at(icat);
    finalcat.defineType(catname);
    
    // check if we're in a sub-cat or the comb-cat
    RooDataSet *data   = NULL;
    RooDataSet *inData = NULL;
    if( icat < (catnames.size() - 1) || catnames.size() == 1) { // this is NOT the last cat (which is by construction the combination)
      inData = (RooDataSet*)win->data(TString("data_mass_")+catname);
      if( !inData ) {
	std::cerr<<" ERROR: Could not find dataset > data_mass_"<<catname.Data()<<" < in input workspace."<<std::endl;
	return;
      }
      data = new RooDataSet(TString("data_")+catname,"",*hmass,Import(*inData));  // copy the dataset (why?)
      
      // append the data to the combined data...
      RooDataSet *datacat = new RooDataSet(TString("datacat")+catname,"",*hmass,Index(finalcat),Import(catname,*data)) ;
      datacomb.append(*datacat);
      datacombcat->append(*data);
      
      // normalization for this category
      RooRealVar *nbkg = new RooRealVar(TString::Format("CMS_hgg_%s_bkgshape_norm",catname.Data()),"",800.0,0.0,25e3);
      
      // we keep track of the normalizario vars only for N-1 cats, naming convetnions hystoric...
      if( catnames.size() > 2 && icat < (catnames.size() - 2) ) {
	RooRealVar* cbkg = new RooRealVar(TString::Format("cbkg%s",catname.Data()),"",0.0,0.0,1e3);
	cbkg->removeRange();
	normu_vec.push_back(cbkg);
	normList.add(*cbkg);
      }
      
      /// generate the Bernstrin polynomial (FIX-ME: add possibility ro create other models...)
      fstBernModel* theBGmodel = new fstBernModel(hmass, polorder[icat], icat, catname);            // using my dedicated class...
      
      std::cout<<" model name is "<<theBGmodel->getPdf()->GetName()<<std::endl;

      RooAbsPdf*    bkgshape   = theBGmodel->getPdf();                                              // the BG shape
      RooAbsPdf*    bkgpdf     = new RooExtendPdf(TString("bkgpdf")+catname,"",*bkgshape,*nbkg);    // the extended PDF
      
      // add the extedned PDF to the RooSimultaneous holding all models...
      fullbkgpdf.addPdf(*bkgpdf,catname);
      // store the NON-EXTENDED PDF for usgae to compute the error bands later..
      pdfShape_vec.push_back(bkgshape);
      pdf_vec     .push_back(bkgpdf);
      data_vec    .push_back(data);
      
    } else {
      data = datacombcat;   // we're looking at the last cat (by construction the combination)
      data_vec.push_back(data);
      

int main(){
  
  system("mkdir -p plots");
  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
  TFile *bkgFile = TFile::Open("comb_svn/hgg.inputbkgdata_8TeV_MVA.root");
  TFile *sigFile = TFile::Open("comb_svn/hgg.inputsig_8TeV_nosplitVH_MVA.root");
  RooWorkspace *bkgWS = (RooWorkspace*)bkgFile->Get("cms_hgg_workspace");
  RooWorkspace *sigWS = (RooWorkspace*)sigFile->Get("wsig_8TeV");

  RooRealVar *mass = (RooRealVar*)bkgWS->var("CMS_hgg_mass");
  RooRealVar *mu = new RooRealVar("mu","mu",-5.,5.);

  mass->setBins(320);
  cout << mass->getBins() << endl;
  RooDataSet *dataAll;
  
  int firstCat=1;
  int lastCat=1;
  float mu_low=-1.;
  float mu_high=3.;
  float mu_step=0.01;

  vector<pair<double,TGraph*> > minNlltrack;

  for (int cat=firstCat; cat<=lastCat; cat++){
    RooDataSet *data = (RooDataSet*)bkgWS->data(Form("data_mass_cat%d",cat));
    if (cat==firstCat) dataAll = (RooDataSet*)data->Clone("data_mass_all");
    else dataAll->append(*data);
    RooDataHist *dataBinned = new RooDataHist(Form("roohist_data_mass_cat%d",cat),Form("roohist_data_mass_cat%d",cat),RooArgSet(*mass),*data);
    RooDataSet *sigMC = (RooDataSet*)sigWS->data(Form("sig_mass_m125_cat%d",cat));

    if (!dataBinned || !sigMC){
      cerr << "ERROR -- one of data or signal is NULL" << endl;
      exit(1);
    }
    
    // Construct PDFs for this category using PdfModelBuilder
    PdfModelBuilder modelBuilder;
    modelBuilder.setObsVar(mass);
    modelBuilder.setSignalModifier(mu);
    // For Standard Analysis
    //if (cat>=0 && cat<=3) modelBuilder.addBkgPdf("Bernstein",5,Form("pol5_cat%d",cat));
    //if (cat>=4 && cat<=5) modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
    //if (cat>=6 && cat<=8) modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
    // To Profile Multiple PDFs
    if (cat==0 || cat==1 || cat==2 || cat==3){
      modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",5,Form("pol5_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",6,Form("pol6_cat%d",cat));
      /*
      modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",5,Form("pow5_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",5,Form("exp5_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",5,Form("lau5_cat%d",cat));
      */
    }
    if (cat==4 || cat==5 || cat==6 || cat==7 || cat==8) {
      modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
      /*
      modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
      */
    }
    map<string,RooAbsPdf*> bkgPdfs = modelBuilder.getBkgPdfs();
    modelBuilder.setSignalPdfFromMC(sigMC);
    modelBuilder.makeSBPdfs();
    map<string,RooAbsPdf*> sbPdfs = modelBuilder.getSBPdfs();

    modelBuilder.fitToData(dataBinned,true,true);
    modelBuilder.fitToData(dataBinned,false,true);

    modelBuilder.throwToy(Form("cat%d_toy0",cat),dataBinned->sumEntries(),true,true);

    // Profile this category using ProfileMultiplePdfs
    ProfileMultiplePdfs profiler;
    for (map<string,RooAbsPdf*>::iterator pdf=sbPdfs.begin(); pdf!=sbPdfs.end(); pdf++) {
      string bkgOnlyName = pdf->first.substr(pdf->first.find("sb_")+3,string::npos);
      if (bkgPdfs.find(bkgOnlyName)==bkgPdfs.end()){
        cerr << "ERROR -- couldn't find bkg only pdf " << bkgOnlyName << " for SB pdf " << pdf->first << endl;
        pdf->second->fitTo(*dataBinned);
        exit(1);
      }
      int nParams = bkgPdfs[bkgOnlyName]->getVariables()->getSize()-1;
      profiler.addPdf(pdf->second,2*nParams);
      //profiler.addPdf(pdf->second);
      cout << pdf->second->GetName() << " nParams=" << pdf->second->getVariables()->getSize() << " nBkgParams=" << nParams << endl;
    }
    profiler.printPdfs();
    //cout << "Continue?" << endl;
    //string bus; cin >> bus;
//.........这里部分代码省略.........

string runQuickRejig(string oldfilename, int ncats){
  
  string newfilename="TestWS.root";
  TFile *oldFile = TFile::Open(oldfilename.c_str());
  TFile *newFile = new TFile(newfilename.c_str(),"RECREATE");

  RooWorkspace *oldWS = (RooWorkspace*)oldFile->Get("cms_hgg_workspace");
  RooWorkspace *newWS = new RooWorkspace("wsig_8TeV");
  
  RooRealVar *mass = (RooRealVar*)oldWS->var("CMS_hgg_mass");

  RooDataSet *tot;

  RooRealVar *normT1 = new RooRealVar("nT1","nT1",500,0,550);
  RooRealVar *meanT1 = new RooRealVar("mT1","mT1",125,122,128);
  RooRealVar *sigmaT1 = new RooRealVar("sT1","sT1",3.1,1.,10.);
  RooGaussian *gausT1 = new RooGaussian("gT1","gT1",*mass,*meanT1,*sigmaT1);
  RooRealVar *normT2 = new RooRealVar("nT2","nT2",2,0,500);
  RooRealVar *meanT2 = new RooRealVar("mT2","mT2",125,122,128);
  RooRealVar *sigmaT2 = new RooRealVar("sT2","sT2",1.7,0.,10.);
  RooGaussian *gausT2 = new RooGaussian("gT2","gT2",*mass,*meanT2,*sigmaT2);
  RooRealVar *normT3 = new RooRealVar("nT3","nT3",2,0,500);
  RooRealVar *meanT3 = new RooRealVar("mT3","mT3",125,122,128);
  RooRealVar *sigmaT3 = new RooRealVar("sT3","sT3",1.7,0.,10.);
  RooGaussian *gausT3 = new RooGaussian("gT3","gT3",*mass,*meanT3,*sigmaT3);
  RooAddPdf *gausT = new RooAddPdf("gT","gT",RooArgList(*gausT1,*gausT2,*gausT3),RooArgList(*normT1,*normT2,*normT3));

  for (int cat=0; cat<ncats; cat++){
    RooDataSet *thisData = (RooDataSet*)oldWS->data(Form("sig_ggh_mass_m125_cat%d",cat));
    newWS->import(*thisData);
    RooRealVar *norm1 = new RooRealVar(Form("n1%d",cat),Form("n1%d",cat),500,0,550);
    RooRealVar *mean1 = new RooRealVar(Form("m1%d",cat),Form("m1%d",cat),125,122,128);
    RooRealVar *sigma1 = new RooRealVar(Form("s1%d",cat),Form("s1%d",cat),3.1,1.,10.);
    RooGaussian *gaus1 = new RooGaussian(Form("g1%d",cat),Form("g1%d",cat),*mass,*mean1,*sigma1);
    RooRealVar *norm2 = new RooRealVar(Form("n2%d",cat),Form("n2%d",cat),2,0,500);
    RooRealVar *mean2 = new RooRealVar(Form("m2%d",cat),Form("m2%d",cat),125,122,128);
    RooRealVar *sigma2 = new RooRealVar(Form("s2%d",cat),Form("s2%d",cat),1.7,0.,10.);
    RooGaussian *gaus2 = new RooGaussian(Form("g2%d",cat),Form("g2%d",cat),*mass,*mean2,*sigma2);
    RooRealVar *norm3 = new RooRealVar(Form("n3%d",cat),Form("n3%d",cat),2,0,500);
    RooRealVar *mean3 = new RooRealVar(Form("m3%d",cat),Form("m3%d",cat),125,122,128);
    RooRealVar *sigma3 = new RooRealVar(Form("s3%d",cat),Form("s3%d",cat),1.7,0.,10.);
    RooGaussian *gaus3 = new RooGaussian(Form("g3%d",cat),Form("g3%d",cat),*mass,*mean3,*sigma3);
    RooAddPdf *gaus = new RooAddPdf(Form("g%d",cat),"g",RooArgList(*gaus1,*gaus2,*gaus3),RooArgList(*norm1,*norm2,*norm3));
    gaus->fitTo(*thisData,SumW2Error(kTRUE));
    newWS->import(*gaus);
    if (cat==0) {
      tot = thisData;
      tot->SetName("sig_ggh_m125");
    }
    else tot->append(*thisData);
  }
  newWS->import(*tot);
  gausT->fitTo(*tot,SumW2Error(kTRUE));
  newWS->import(*gausT);
  newWS->Write();
  oldFile->Close();
  newFile->Close();
  delete newFile;
  delete newWS;

  return newfilename;

}