C++ RooFitResult::covarianceMatrix方法代码示例

MCMCInterval * Tprime::GetMcmcInterval(ModelConfig mc,
                                       double conf_level,
                                       int n_iter,
                                       int n_burn,
                                       double left_side_tail_fraction,
                                       int n_bins) {
    //
    // Bayesian MCMC calculation using arbitrary ModelConfig
    // Want an efficient proposal function, so derive it from covariance
    // matrix of fit
    //

    RooAbsData * _data = data;
    //RooAbsData * _data = pWs->data("obsData");
    //RooStats::ModelConfig * _mc = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");
    RooStats::ModelConfig * _mc = GetModelConfig();
    _mc->Print();

    //RooFitResult * fit = pWs->pdf("model_tprime")->fitTo(*_data,Save());
    RooFitResult * fit = _mc->GetPdf()->fitTo(*_data,Save());
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
    ph.SetCacheSize(100);
    ProposalFunction * pf = ph.GetProposalFunction();

    //delete pf;
    //pf = new SequentialProposal();

    MCMCCalculator mcmc( *_data, mc );
    mcmc.SetConfidenceLevel(conf_level);
    mcmc.SetNumIters(n_iter);          // Metropolis-Hastings algorithm iterations
    mcmc.SetProposalFunction(*pf);
    mcmc.SetNumBurnInSteps(n_burn); // first N steps to be ignored as burn-in
    mcmc.SetLeftSideTailFraction(left_side_tail_fraction);
    mcmc.SetNumBins(n_bins);

    //mcInt = mcmc.GetInterval();
    try {
        mcInt = mcmc.GetInterval();
    } catch ( std::length_error &ex) {
        mcInt = 0;
    }

    //std::cout << "!!!!!!!!!!!!!! interval" << std::endl;
    if (mcInt == 0) std::cout << "No interval found!" << std::endl;

    delete fit;
    delete pf;

    return mcInt;
}

//_________________________________________________
void TestJeffreysGaussSigma(){
  // this one is VERY sensitive
  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
  //   and you get really bizzare shapes
  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
  //   and the PDF falls off too fast at high sigma
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
  w.factory("n[100,.1,2000]");
  w.factory("ExtendPdf::p(g,n)");
  //  w.var("sigma")->setConstant();
  w.var("mu")->setConstant();
  w.var("n")->setConstant();
  w.var("x")->setBins(301);

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;


  //  w.defineSet("poi","mu,sigma");
  //w.defineSet("poi","mu,sigma,n");
  w.defineSet("poi","sigma");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();
  //  return;

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","sqrt(2.)/sigma",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("sigma")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

void JeffreysPriorDemo(){
  RooWorkspace w("w");
  w.factory("Uniform::u(x[0,1])");
  w.factory("mu[100,1,200]");
  w.factory("ExtendPdf::p(u,mu)");

  //  w.factory("Poisson::pois(n[0,inf],mu)");

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
  //  RooDataHist* asimov2 = w.pdf("pois")->generateBinned(*w.var("n"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  w.defineSet("poi","mu");
  w.defineSet("obs","x");
  //  w.defineSet("obs2","n");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",10);

  //  JeffreysPrior pi2("jeffreys2","jeffreys",*w.pdf("pois"),*w.set("poi"),*w.set("obs2"));

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1./sqrt(mu)",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  //  pi.plotOn(plot, Normalization(1,RooAbsReal::Raw),Precision(.1));
  pi.plotOn(plot);
  //  pi2.plotOn(plot,LineColor(kGreen),LineStyle(kDotted));
  test->plotOn(plot,LineColor(kRed));
  plot->Draw();

}

//_________________________________________________
void TestJeffreysGaussMean(){
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
  w.factory("n[10,.1,200]");
  w.factory("ExtendPdf::p(g,n)");
  w.var("sigma")->setConstant();
  w.var("n")->setConstant();

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  //  w.defineSet("poi","mu,sigma");
  w.defineSet("poi","mu");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

MCMCInterval * TwoBody::GetMcmcInterval_OldWay(ModelConfig mc,
					double conf_level,
					int n_iter,
					int n_burn,
					double left_side_tail_fraction,
					int n_bins){
  // use MCMCCalculator  (takes about 1 min)
  // Want an efficient proposal function, so derive it from covariance
  // matrix of fit
  
  RooFitResult* fit = ws->pdf("model")->fitTo(*data,Save());
  ProposalHelper ph;
  ph.SetVariables((RooArgSet&)fit->floatParsFinal());
  ph.SetCovMatrix(fit->covarianceMatrix());
  ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
  ph.SetCacheSize(100);
  ProposalFunction* pf = ph.GetProposalFunction();
  
  MCMCCalculator mcmc( *data, mc );
  mcmc.SetConfidenceLevel(conf_level);
  mcmc.SetNumIters(n_iter);          // Metropolis-Hastings algorithm iterations
  mcmc.SetProposalFunction(*pf);
  mcmc.SetNumBurnInSteps(n_burn); // first N steps to be ignored as burn-in
  mcmc.SetLeftSideTailFraction(left_side_tail_fraction);
  mcmc.SetNumBins(n_bins);
  
//mcInt = mcmc.GetInterval();
  try {
    mcInt = mcmc.GetInterval();
  } catch ( std::length_error &ex) {
    mcInt = 0;
  }

  //std::cout << "!!!!!!!!!!!!!! interval" << std::endl;
  if (mcInt == 0) std::cout << "No interval found!" << std::endl;
  
  return mcInt;
}

//.........这里部分代码省略.........
    RooRealVar sg("sg", "sg", sgVal_);
    RooRealVar a("a", "a", aVal_);
    RooRealVar n("n", "n", nVal_);

    RooCBShape CB("CB","CB",*mZ1,bwMean,sg,a,n);
    RooRealVar f("f","f", fVal_);

    RooRealVar mean("mean","mean",meanVal_);
    RooRealVar sigma("sigma","sigma",sigmaVal_);
    RooRealVar f1("f1","f1",f1Val_);

    RooGenericPdf RelBW("RelBW","1/( pow(mZ1*mZ1-bwMean*bwMean,2)+pow(mZ1,4)*pow(bwGamma/bwMean,2) )", RooArgSet(*mZ1,bwMean,bwGamma) );

    RooAddPdf RelBWxCB("RelBWxCB","RelBWxCB", RelBW, CB, f);
    RooGaussian gauss("gauss","gauss",*mZ1,mean,sigma);
    RooAddPdf RelBWxCBxgauss("RelBWxCBxgauss","RelBWxCBxgauss", RelBWxCB, gauss, f1);

    RooProdPdf *PDFRelBWxCBxgauss;
    PDFRelBWxCBxgauss = new RooProdPdf("PDFRelBWxCBxgauss","PDFRelBWxCBxgauss", 
                                     RooArgList(gauss1, gauss2, RelBWxCBxgauss) );
    if(p4sZ1ph_.size()==1)    
      PDFRelBWxCBxgauss = new RooProdPdf("PDFRelBWxCBxgauss","PDFRelBWxCBxgauss", 
                                     RooArgList(gauss1, gauss2, gaussph1, RelBWxCBxgauss) );
    if(p4sZ1ph_.size()==2)
      PDFRelBWxCBxgauss = new RooProdPdf("PDFRelBWxCBxgauss","PDFRelBWxCBxgauss", 
                                     RooArgList(gauss1, gauss2, gaussph1, gaussph2, RelBWxCBxgauss) );

    // observable set
    RooArgSet *rastmp;
      rastmp = new RooArgSet(*pT1RECO,*pT2RECO);
    if(p4sZ1ph_.size()==1)
      rastmp = new RooArgSet(*pT1RECO,*pT2RECO,*pTph1RECO);
    if(p4sZ1ph_.size()>=2)
      rastmp = new RooArgSet(*pT1RECO,*pT2RECO,*pTph1RECO,*pTph2RECO);

    RooDataSet* pTs = new RooDataSet("pTs","pTs", *rastmp);
    pTs->add(*rastmp); 

    //RooAbsReal* nll;
    //nll = PDFRelBWxCBxgauss->createNLL(*pTs);
    //RooMinuit(*nll).migrad();

    RooFitResult* r = PDFRelBWxCBxgauss->fitTo(*pTs,RooFit::Save(),RooFit::PrintLevel(-1));
    const TMatrixDSym& covMatrix = r->covarianceMatrix();
   
    const RooArgList& finalPars = r->floatParsFinal();
    for (int i=0 ; i<finalPars.getSize(); i++){
        TString name = TString(((RooRealVar*)finalPars.at(i))->GetName());

        if(debug_) cout<<"name list of RooRealVar for covariance matrix "<<name<<endl;

    }

    int size = covMatrix.GetNcols();
    //TMatrixDSym covMatrixTest_(size);
    covMatrixZ1_.ResizeTo(size,size);
    covMatrixZ1_ = covMatrix;   

    if(debug_) cout<<"save the covariance matrix"<<endl;
    
    l1 = pT1->getVal()/RECOpT1; l2 = pT2->getVal()/RECOpT2;
    double pTerrZ1REFIT1 = pT1->getError(); double pTerrZ1REFIT2 = pT2->getError();

    pTerrsZ1REFIT_.push_back(pTerrZ1REFIT1);
    pTerrsZ1REFIT_.push_back(pTerrZ1REFIT2);

    if(p4sZ1ph_.size()>=1){

      if(debug_) cout<<"set refit result for Z1 fsr photon 1"<<endl;

      lph1 = pTph1->getVal()/RECOpTph1;
      double pTerrZ1phREFIT1 = pTph1->getError();
      if(debug_) cout<<"scale "<<lph1<<" pterr "<<pTerrZ1phREFIT1<<endl;  
   
      pTerrsZ1phREFIT_.push_back(pTerrZ1phREFIT1);

    } 
    if(p4sZ1ph_.size()==2){

      lph2 = pTph2->getVal()/RECOpTph2;
      double pTerrZ1phREFIT2 = pTph2->getError();
      pTerrsZ1phREFIT_.push_back(pTerrZ1phREFIT2);

    }

    //delete nll;
    delete r;
    delete mZ1;
    delete pT1; delete pT2; delete pTph1; delete pTph2;
    delete pT1RECO; delete pT2RECO; delete pTph1RECO; delete pTph2RECO;
    delete ph1v3Dph2; delete p1v3Dph1; delete p2v3Dph1; delete p1v3Dph2; delete p2v3Dph2;
    delete PDFRelBWxCBxgauss;
    delete pTs;
    delete rastmp;

    if(debug_) cout<<"end Z1 refit"<<endl;

    return 0;

}

//.........这里部分代码省略.........
//     MyLimit result(plInt->IsInInterval(exp_sig),
    MyLimit result(exp_sig_val>lowLim&&exp_sig_val<uppLim,lowLim,uppLim);
    // std::cout << "isIn " << result << std::endl;
    delete plInt;
//     delete modelConfig;
    return result;
  }

  // use FeldmaCousins (takes ~20 min)  
  if ( method == FeldmanCousinsMethod ) {
    FeldmanCousins fc(*data, modelConfig);
    fc.SetConfidenceLevel(0.95);
    //number counting: dataset always has 1 entry with N events observed
    fc.FluctuateNumDataEntries(false); 
    fc.UseAdaptiveSampling(true);
    fc.SetNBins(100);
    PointSetInterval* fcInt = NULL;
    fcInt = (PointSetInterval*) fc.GetInterval(); // fix cast
    double lowLim = fcInt->LowerLimit(*wspace->var("s"));
    double uppLim = fcInt->UpperLimit(*wspace->var("s"));
//     double exp_sig_val = wspace->var("s")->getVal();
    cout << "Feldman Cousins interval on s = [" << lowLim << " " << uppLim << endl;
    // std::cout << "isIn " << result << std::endl;
    MyLimit result(exp_sig_val>lowLim&&exp_sig_val<uppLim,
		   fcInt->LowerLimit(*wspace->var("s")),fcInt->UpperLimit(*wspace->var("s")));
    delete fcInt;
    return result;
  }


  // use BayesianCalculator (only 1-d parameter of interest, slow for this problem)  
  if ( method == BayesianMethod ) {
    BayesianCalculator bc(*data, modelConfig);
    bc.SetConfidenceLevel(0.95);
    bc.SetLeftSideTailFraction(0.5);
    SimpleInterval* bInt = NULL;
    if( wspace->set("poi")->getSize() == 1)   {
      bInt = bc.GetInterval();
      if ( draw ) {
	TCanvas* c = new TCanvas("Bayesian");
	// the plot takes a long time and print lots of error
	// using a scan it is better
	bc.SetScanOfPosterior(50);
	RooPlot* bplot = bc.GetPosteriorPlot();
	bplot->Draw();
      }
      cout << "Bayesian interval on s = [" << 
	bInt->LowerLimit( ) << ", " <<
	bInt->UpperLimit( ) << "]" << endl;
      // std::cout << "isIn " << result << std::endl;
      MyLimit result(bInt->IsInInterval(exp_sig),
		     bInt->LowerLimit(),bInt->UpperLimit());
      delete bInt;
      return result;
    } else {
    cout << "Bayesian Calc. only supports on parameter of interest" << endl;
    return MyLimit();
    }
  }


  // use MCMCCalculator  (takes about 1 min)
  // Want an efficient proposal function, so derive it from covariance
  // matrix of fit
  if ( method == MCMCMethod ) {
    RooFitResult* fit = wspace->pdf("model")->fitTo(*data,Save());
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(kTRUE); // auto-create mean vars and add mappings
    ph.SetCacheSize(100);
    ProposalFunction* pf = ph.GetProposalFunction();
    
    MCMCCalculator mc(*data, modelConfig);
    mc.SetConfidenceLevel(0.95);
    mc.SetProposalFunction(*pf);
    mc.SetNumBurnInSteps(100); // first N steps to be ignored as burn-in
    mc.SetNumIters(100000);
    mc.SetLeftSideTailFraction(0.5); // make a central interval
    MCMCInterval* mcInt = NULL;
    mcInt = mc.GetInterval();
    MCMCIntervalPlot mcPlot(*mcInt); 
    mcPlot.Draw();
    cout << "MCMC interval on s = [" << 
      mcInt->LowerLimit(*wspace->var("s") ) << ", " <<
      mcInt->UpperLimit(*wspace->var("s") ) << "]" << endl;
    // std::cout << "isIn " << result << std::endl;
    MyLimit result(mcInt->IsInInterval(exp_sig),
		   mcInt->LowerLimit(*wspace->var("s")),mcInt->UpperLimit(*wspace->var("s")));
    delete mcInt;
    return result;
  }
  

  t.Print();

//   delete modelConfig;
  return MyLimit();

}

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

        mZ = new RooFormulaVar("mZ", "TMath::Sqrt(2*@0+2*@1+2*@2+2*@3+2*@4+2*@[email protected]*@[email protected]*@7)", RooArgList(p1D2,p1Dph1,p2Dph1,p1Dph2,p2Dph2,ph1Dph2, m1, m2));
        RelBW = new RooGenericPdf("RelBW","1/( pow(mZ*mZ-bwMean*bwMean,2)+pow(mZ,4)*pow(bwGamma/bwMean,2) )", RooArgSet(*mZ,bwMean,bwGamma) );
//        PDFRelBW = new RooProdPdf("PDFRelBW", "PDFRelBW", RooArgList(gauss1_lep, gauss2_lep, gauss1_gamma, gauss2_gamma, *RelBW));

        }

     //true shape
     RooRealVar sg("sg", "sg", sgVal_);
     RooRealVar a("a", "a", aVal_);
     RooRealVar n("n", "n", nVal_);

     RooCBShape CB("CB","CB",*mZ,bwMean,sg,a,n);
     RooRealVar f("f","f", fVal_);

     RooRealVar mean("mean","mean",meanVal_);
     RooRealVar sigma("sigma","sigma",sigmaVal_);
     RooRealVar f1("f1","f1",f1Val_);

     RooAddPdf *RelBWxCB;
     RelBWxCB = new RooAddPdf("RelBWxCB","RelBWxCB", *RelBW, CB, f);
     RooGaussian *gauss;
     gauss = new RooGaussian("gauss","gauss",*mZ,mean,sigma);
     RooAddPdf *RelBWxCBxgauss;
     RelBWxCBxgauss = new RooAddPdf("RelBWxCBxgauss","RelBWxCBxgauss", *RelBWxCB, *gauss, f1);

     RooProdPdf *PDFRelBWxCBxgauss;
     PDFRelBWxCBxgauss = new RooProdPdf("PDFRelBWxCBxgauss","PDFRelBWxCBxgauss", 
                                     RooArgList(gauss1_lep, gauss2_lep, *RelBWxCBxgauss) );


    //make fit
    RooArgSet *rastmp;
    rastmp = new RooArgSet(pTRECO1_lep, pTRECO2_lep);
/*
    if(input.nFsr == 1) {
      rastmp = new RooArgSet(pTRECO1_lep, pTRECO2_lep, pTRECO1_gamma);
      }

    if(input.nFsr == 2) {
      rastmp = new RooArgSet(pTRECO1_lep, pTRECO2_lep, pTRECO1_gamma, pTRECO2_gamma);
      }
*/
    RooDataSet* pTs = new RooDataSet("pTs","pTs", *rastmp);
    pTs->add(*rastmp);

    RooFitResult* r;
    if (mass4lRECO_ > 140) {

       r = PDFRelBW->fitTo(*pTs,RooFit::Save(),RooFit::PrintLevel(-1));

       } else {

              r = PDFRelBWxCBxgauss->fitTo(*pTs,RooFit::Save(),RooFit::PrintLevel(-1));

              }
    //save fit result
    const TMatrixDSym& covMatrix = r->covarianceMatrix();
    const RooArgList& finalPars = r->floatParsFinal();

    for (int i=0 ; i<finalPars.getSize(); i++){
 
        TString name = TString(((RooRealVar*)finalPars.at(i))->GetName());
        if(debug_) cout<<"name list of RooRealVar for covariance matrix "<<name<<endl;

    }

    int size = covMatrix.GetNcols();
    output.covMatrixZ.ResizeTo(size,size);
    output.covMatrixZ = covMatrix;
    
    output.pT1_lep = pTMean1_lep.getVal();
    output.pT2_lep = pTMean2_lep.getVal();
    output.pTErr1_lep = pTMean1_lep.getError();
    output.pTErr2_lep = pTMean2_lep.getError();
/*
    if (input.nFsr >= 1) {

       output.pT1_gamma = pTMean1_gamma.getVal();
       output.pTErr1_gamma = pTMean1_gamma.getError();
    
       }

    if (input.nFsr == 2) {

       output.pT2_gamma = pTMean2_gamma.getVal();
       output.pTErr2_gamma = pTMean2_gamma.getError();

       }
*/
    delete rastmp;
    delete pTs;
    delete PDFRelBW;
    delete mZ;
    delete RelBW;
    delete RelBWxCB;
    delete gauss;
    delete RelBWxCBxgauss;
    delete PDFRelBWxCBxgauss;
}

void rs101_limitexample()
{
    // --------------------------------------
    // An example of setting a limit in a number counting experiment with uncertainty on background and signal

    // to time the macro
    TStopwatch t;
    t.Start();

    // --------------------------------------
    // The Model building stage
    // --------------------------------------
    RooWorkspace* wspace = new RooWorkspace();
    wspace->factory("Poisson::countingModel(obs[150,0,300], sum(s[50,0,120]*ratioSigEff[1.,0,3.],b[100]*ratioBkgEff[1.,0.,3.]))"); // counting model
    //  wspace->factory("Gaussian::sigConstraint(ratioSigEff,1,0.05)"); // 5% signal efficiency uncertainty
    //  wspace->factory("Gaussian::bkgConstraint(ratioBkgEff,1,0.1)"); // 10% background efficiency uncertainty
    wspace->factory("Gaussian::sigConstraint(gSigEff[1,0,3],ratioSigEff,0.05)"); // 5% signal efficiency uncertainty
    wspace->factory("Gaussian::bkgConstraint(gSigBkg[1,0,3],ratioBkgEff,0.2)"); // 10% background efficiency uncertainty
    wspace->factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)"); // product of terms
    wspace->Print();

    RooAbsPdf* modelWithConstraints = wspace->pdf("modelWithConstraints"); // get the model
    RooRealVar* obs = wspace->var("obs"); // get the observable
    RooRealVar* s = wspace->var("s"); // get the signal we care about
    RooRealVar* b = wspace->var("b"); // get the background and set it to a constant.  Uncertainty included in ratioBkgEff
    b->setConstant();

    RooRealVar* ratioSigEff = wspace->var("ratioSigEff"); // get uncertain parameter to constrain
    RooRealVar* ratioBkgEff = wspace->var("ratioBkgEff"); // get uncertain parameter to constrain
    RooArgSet constrainedParams(*ratioSigEff, *ratioBkgEff); // need to constrain these in the fit (should change default behavior)

    RooRealVar * gSigEff = wspace->var("gSigEff");     // global observables for signal efficiency
    RooRealVar * gSigBkg = wspace->var("gSigBkg");  // global obs for background efficiency
    gSigEff->setConstant();
    gSigBkg->setConstant();

    // Create an example dataset with 160 observed events
    obs->setVal(160.);
    RooDataSet* data = new RooDataSet("exampleData", "exampleData", RooArgSet(*obs));
    data->add(*obs);

    RooArgSet all(*s, *ratioBkgEff, *ratioSigEff);

    // not necessary
    modelWithConstraints->fitTo(*data, RooFit::Constrain(RooArgSet(*ratioSigEff, *ratioBkgEff)));

    // Now let's make some confidence intervals for s, our parameter of interest
    RooArgSet paramOfInterest(*s);

    ModelConfig modelConfig(wspace);
    modelConfig.SetPdf(*modelWithConstraints);
    modelConfig.SetParametersOfInterest(paramOfInterest);
    modelConfig.SetNuisanceParameters(constrainedParams);
    modelConfig.SetObservables(*obs);
    modelConfig.SetGlobalObservables( RooArgSet(*gSigEff,*gSigBkg));
    modelConfig.SetName("ModelConfig");
    wspace->import(modelConfig);
    wspace->import(*data);
    wspace->SetName("w");
    wspace->writeToFile("rs101_ws.root");



    // First, let's use a Calculator based on the Profile Likelihood Ratio
    //ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest);
    ProfileLikelihoodCalculator plc(*data, modelConfig);
    plc.SetTestSize(.05);
    ConfInterval* lrinterval = plc.GetInterval();  // that was easy.

    // Let's make a plot
    TCanvas* dataCanvas = new TCanvas("dataCanvas");
    dataCanvas->Divide(2,1);

    dataCanvas->cd(1);
    LikelihoodIntervalPlot plotInt((LikelihoodInterval*)lrinterval);
    plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S");
    plotInt.Draw();

    // Second, use a Calculator based on the Feldman Cousins technique
    FeldmanCousins fc(*data, modelConfig);
    fc.UseAdaptiveSampling(true);
    fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed
    fc.SetNBins(100); // number of points to test per parameter
    fc.SetTestSize(.05);
    //  fc.SaveBeltToFile(true); // optional
    ConfInterval* fcint = NULL;
    fcint = fc.GetInterval();  // that was easy.

    RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true));

    // Third, use a Calculator based on Markov Chain monte carlo
    // Before configuring the calculator, let's make a ProposalFunction
    // that will achieve a high acceptance rate
    ProposalHelper ph;
    ph.SetVariables((RooArgSet&)fit->floatParsFinal());
    ph.SetCovMatrix(fit->covarianceMatrix());
    ph.SetUpdateProposalParameters(true);
    ph.SetCacheSize(100);
    ProposalFunction* pdfProp = ph.GetProposalFunction();  // that was easy

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

void addNuisanceWithToys(std::string iFileName,std::string iChannel,std::string iBkg,std::string iEnergy,std::string iName,std::string iDir,bool iRebin=true,bool iVarBin=false,int iFitModel=1,int iFitModel1=1,double iFirst=150,double iLast=1500,std::string iSigMass="800",double iSigScale=0.1,int iNToys=1000) { 
  std::cout << "======> " << iDir << "/" << iBkg << " -- " << iFileName << std::endl;  
  if(iVarBin) std::cout << "option not implemented yet!";
  if(iVarBin) return;
  //double lFirst = 200;
  //double lLast  = 1500;
  double lFirst = iFirst;
  double lLast  = iLast;

  std::cout << "===================================================================================================================================================" <<std::endl;
  std::cout << "Using Initial fit model: " << iFitModel << ", fitting range: " << iFirst << "-" << iLast << " , using alternative fit model: " << iFitModel1 << std::endl; 
  std::cout << "===================================================================================================================================================" <<std::endl;

  TFile *lFile = new TFile(iFileName.c_str());
  TH1F  *lH0   = (TH1F*) lFile->Get((iDir+"/"+iBkg).c_str());
  TH1F  *lData = (TH1F*) lFile->Get((iDir+"/data_obs").c_str());
  TH1F  *lSig = 0;

  // for now, use bbH signal for testing in b-tag and ggH in no-btag
  if(iDir.find("_btag") != std::string::npos) lSig = (TH1F*)lFile->Get((iDir+"/bbH"+iSigMass+"_fine_binning").c_str());
  else lSig = (TH1F*)lFile->Get((iDir+"/ggH"+iSigMass+"_fine_binning").c_str());

  TH1F *lH0Clone = (TH1F*)lH0->Clone("lH0Clone");     // binning too fine as of now? start by rebinning
  TH1F *lDataClone = (TH1F*)lData->Clone("lDataClone");   
  TH1F *lSigClone = (TH1F*)lSig->Clone("lSigClone");  
 // lH0Clone->Rebin(2);
 // lDataClone->Rebin(2);
 // lSigClone->Rebin(2);

  lSig->Rebin(10);  

  //Define the fit function
  RooRealVar lM("m","m" ,0,5000);
  lM.setRange(lFirst,lLast);
  RooRealVar lA("a","a" ,50,  0.1,200);
  RooRealVar lB("b","b" ,0.0 , -10.5,10.5);
  RooRealVar lA1("a1","a1" ,50,  0.1,1000);
  RooRealVar lB1("b1","b1" ,0.0 , -10.5,10.5);

  RooDataHist *pH0  =  new RooDataHist("Data","Data" ,RooArgList(lM),lH0);
  double lNB0 = lH0->Integral(lH0->FindBin(lFirst),lH0->FindBin(lLast));
  double lNSig0 = lSig->Integral(lSig->FindBin(lFirst),lSig->FindBin(lLast));
 //lNB0=500;
// lNSig0=500;
 lSig->Scale(iSigScale*lNB0/lNSig0);                                         // scale signal to iSigScale*(Background yield), could try other options
 lNSig0 = lSig->Integral(lSig->FindBin(lFirst),lSig->FindBin(lLast));        // readjust norm of signal hist   
  
  //Generate the "default" fit model 

  RooGenericPdf *lFit  = 0; lFit = new RooGenericPdf("genPdf","exp(-m/(a+b*m))",RooArgList(lM,lA,lB));
  if(iFitModel == 1) lFit = new RooGenericPdf("genPdf","exp(-a*pow(m,b))",RooArgList(lM,lA,lB));
  if(iFitModel == 1) {lA.setVal(0.3); lB.setVal(0.5);}
  if(iFitModel == 2) lFit = new RooGenericPdf("genPdf","a*exp(b*m)",RooArgList(lM,lA,lB));
  if(iFitModel == 2) {lA.setVal(0.01); lA.setRange(0,10); }
  if(iFitModel == 3) lFit = new RooGenericPdf("genPdf","a/pow(m,b)",RooArgList(lM,lA,lB));
 
  // Generate the alternative model
  
  RooGenericPdf *lFit1  = 0; lFit1 = new RooGenericPdf("genPdf","exp(-m/(a1+b1*m))",RooArgList(lM,lA1,lB1));
  if(iFitModel1 == 1) lFit1 = new RooGenericPdf("genPdf","exp(-a1*pow(m,b1))",RooArgList(lM,lA1,lB1));
  if(iFitModel1 == 1) {lA1.setVal(0.3); lB1.setVal(0.5);}
  if(iFitModel1 == 2) lFit1 = new RooGenericPdf("genPdf","a1*exp(b1*m)",RooArgList(lM,lA1,lB1));
  if(iFitModel1 == 2) {lA1.setVal(0.01); lA1.setRange(0,10); }
  if(iFitModel1 == 3) lFit1 = new RooGenericPdf("genPdf","a1/pow(m,b1)",RooArgList(lM,lA1,lB1));
  
  //=============================================================================================================================================
  //Perform the tail fit and generate the shift up and down histograms
  //=============================================================================================================================================

  RooFitResult  *lRFit = 0;
  lRFit = lFit->fitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(lFirst,lLast),RooFit::Strategy(0)); 
  TMatrixDSym lCovMatrix   = lRFit->covarianceMatrix(); 
  TMatrixD  lEigVecs(2,2);    lEigVecs = TMatrixDSymEigen(lCovMatrix).GetEigenVectors();
  TVectorD  lEigVals(2);      lEigVals = TMatrixDSymEigen(lCovMatrix).GetEigenValues();
  cout << " Ve---> " << lEigVecs(0,0) << " -- " << lEigVecs(1,0) << " -- " << lEigVecs(0,1) << " -- " << lEigVecs(1,1) << endl;
  cout << " Co---> " << lCovMatrix(0,0) << " -- " << lCovMatrix(1,0) << " -- " << lCovMatrix(0,1) << " -- " << lCovMatrix(1,1) << endl;
  double lACentral = lA.getVal();
  double lBCentral = lB.getVal();
  lEigVals(0) = sqrt(lEigVals(0));
  lEigVals(1) = sqrt(lEigVals(1));
  cout << "===> " << lEigVals(0) << " -- " << lEigVals(1) << endl;
  
  TH1F* lH     = (TH1F*) lFit->createHistogram("fit" ,lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  lA.setVal(lACentral + lEigVals(0)*lEigVecs(0,0));
  lB.setVal(lBCentral + lEigVals(0)*lEigVecs(1,0));
  TH1F* lHUp   = (TH1F*) lFit->createHistogram("Up"  ,lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  lA.setVal(lACentral - lEigVals(0)*lEigVecs(0,0));
  lB.setVal(lBCentral - lEigVals(0)*lEigVecs(1,0));
  TH1F* lHDown = (TH1F*) lFit->createHistogram("Down",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  lA.setVal(lACentral + lEigVals(1)*lEigVecs(0,1));
  lB.setVal(lBCentral + lEigVals(1)*lEigVecs(1,1));
  TH1F* lHUp1   = (TH1F*) lFit->createHistogram("Up1",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  lA.setVal(lACentral - lEigVals(1)*lEigVecs(0,1));
  lB.setVal(lBCentral - lEigVals(1)*lEigVecs(1,1));
  TH1F* lHDown1 = (TH1F*) lFit->createHistogram("Down1",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

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

void addNuisance(std::string iFileName,std::string iChannel,std::string iBkg,std::string iEnergy,std::string iName,std::string iDir,bool iRebin=true,bool iVarBin=false,int iFitModel=1,double iFirst=150,double iLast=1500) { 
  std::cout << "======> " << iDir << "/" << iBkg << " -- " << iFileName << std::endl;  
  if(iVarBin) addVarBinNuisance(iFileName,iChannel,iBkg,iEnergy,iName,iDir,iRebin,iFitModel,iFirst,iLast);
  if(iVarBin) return;

  TFile *lFile = new TFile(iFileName.c_str());
  TH1F  *lH0   = (TH1F*) lFile->Get((iDir+"/"+iBkg).c_str());
  TH1F  *lData = (TH1F*) lFile->Get((iDir+"/data_obs").c_str());

  //Define the fit function
  RooRealVar lM("m","m" ,0,5000);   //lM.setBinning(lBinning);
  RooRealVar lA("a","a" ,50,  0.1,100);
  RooRealVar lB("b","b" ,0.0 , -10.5,10.5); //lB.setConstant(kTRUE);
  RooDataHist *pH0  =  new RooDataHist("Data","Data" ,RooArgList(lM),lH0);
  RooGenericPdf *lFit  = 0; lFit = new RooGenericPdf("genPdf","exp(-m/(a+b*m))",RooArgList(lM,lA,lB));
  if(iFitModel == 1) lFit = new RooGenericPdf("genPdf","exp(-a*pow(m,b))",RooArgList(lM,lA,lB));
  if(iFitModel == 1) {lA.setVal(0.3); lB.setVal(0.5);}
  if(iFitModel == 2) lFit = new RooGenericPdf("genPdf","a*exp(b*m)",RooArgList(lM,lA,lB));
  if(iFitModel == 3) lFit = new RooGenericPdf("genPdf","a/pow(m,b)",RooArgList(lM,lA,lB));
  RooFitResult  *lRFit = 0;
  double lFirst = iFirst;
  double lLast  = iLast;
  //lRFit = lFit->chi2FitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(lFirst,lLast));
  lRFit = lFit->fitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(lFirst,lLast),RooFit::Strategy(0)); 
  TMatrixDSym lCovMatrix   = lRFit->covarianceMatrix(); 
  TMatrixD  lEigVecs(2,2);    lEigVecs = TMatrixDSymEigen(lCovMatrix).GetEigenVectors();
  TVectorD  lEigVals(2);      lEigVals = TMatrixDSymEigen(lCovMatrix).GetEigenValues();
  cout << " Ve---> " << lEigVecs(0,0) << " -- " << lEigVecs(1,0) << " -- " << lEigVecs(0,1) << " -- " << lEigVecs(1,1) << endl;
  cout << " Co---> " << lCovMatrix(0,0) << " -- " << lCovMatrix(1,0) << " -- " << lCovMatrix(0,1) << " -- " << lCovMatrix(1,1) << endl;
  double lACentral = lA.getVal();
  double lBCentral = lB.getVal();
  lEigVals(0) = sqrt(lEigVals(0));
  lEigVals(1) = sqrt(lEigVals(1));
  cout << "===> " << lEigVals(0) << " -- " << lEigVals(1) << endl;
  
  TH1F* lH     = (TH1F*) lFit->createHistogram("fit" ,lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));
  lA.setVal(lACentral + lEigVals(0)*lEigVecs(0,0));
  lB.setVal(lBCentral + lEigVals(0)*lEigVecs(1,0));
  TH1F* lHUp   = (TH1F*) lFit->createHistogram("Up"  ,lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));
  lA.setVal(lACentral - lEigVals(0)*lEigVecs(0,0));
  lB.setVal(lBCentral - lEigVals(0)*lEigVecs(1,0));
  TH1F* lHDown = (TH1F*) lFit->createHistogram("Down",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  lA.setVal(lACentral + lEigVals(1)*lEigVecs(0,1));
  lB.setVal(lBCentral + lEigVals(1)*lEigVecs(1,1));
  TH1F* lHUp1   = (TH1F*) lFit->createHistogram("Up1",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));
  lA.setVal(lACentral - lEigVals(1)*lEigVecs(0,1));
  lB.setVal(lBCentral - lEigVals(1)*lEigVecs(1,1));
  TH1F* lHDown1 = (TH1F*) lFit->createHistogram("Down1",lM,RooFit::Binning(lH0->GetNbinsX(),lH0->GetXaxis()->GetXmin(),lH0->GetXaxis()->GetXmax()));

  std::string lNuisance1 =  iBkg+"_"+"CMS_"+iName+"1_" + iChannel + "_" + iEnergy;
  std::string lNuisance2 =  iBkg+"_"+"CMS_"+iName+"2_" + iChannel + "_" + iEnergy;
  lHUp    = merge(lNuisance1 + "Up"   ,lFirst,lH0,lHUp);
  lHDown  = merge(lNuisance1 + "Down" ,lFirst,lH0,lHDown);
  lHUp1   = merge(lNuisance2 + "Up"   ,lFirst,lH0,lHUp1);
  lHDown1 = merge(lNuisance2 + "Down" ,lFirst,lH0,lHDown1);
  lH      = merge(lH0->GetName()      ,lFirst,lH0,lH);

  if(iRebin) { 
    const int lNBins = lData->GetNbinsX();
    double *lAxis    = getAxis(lData);
    lH0     = rebin(lH0    ,lNBins,lAxis);
    lH      = rebin(lH     ,lNBins,lAxis);
    lHUp    = rebin(lHUp   ,lNBins,lAxis);
    lHDown  = rebin(lHDown ,lNBins,lAxis);
    lHUp1   = rebin(lHUp1  ,lNBins,lAxis);
    lHDown1 = rebin(lHDown1,lNBins,lAxis);
  }

  // we dont need this bin errors since we do not use them (fit tails replaces bin-by-bin error!), therefore i set all errors to 0, this also saves us from modifying the add_bbb_error.py script in which I otherwise would have to include a option for adding bbb only in specific ranges
  int lMergeBin = lH->GetXaxis()->FindBin(iFirst);
  for(int i0 = lMergeBin; i0 < lH->GetNbinsX()+1; i0++){
    lH->SetBinError  (i0,0);
    lHUp->SetBinError  (i0,0);
    lHDown->SetBinError  (i0,0);
    lHUp1->SetBinError  (i0,0);
    lHDown1->SetBinError  (i0,0);
  }


  TFile *lOutFile =new TFile("Output.root","RECREATE");
  cloneFile(lOutFile,lFile,iDir+"/"+iBkg);
  lOutFile->cd(iDir.c_str());
  lH     ->Write();
  lHUp   ->Write(); 
  lHDown ->Write(); 
  lHUp1  ->Write(); 
  lHDown1->Write(); 

  // Debug Plots
  lH0->SetStats(0);
  lH->SetStats(0);
  lHUp->SetStats(0);
  lHDown->SetStats(0);
  lHUp1->SetStats(0);
  lHDown1->SetStats(0);
  lH0    ->SetLineWidth(1); lH0->SetMarkerStyle(kFullCircle);
  lH     ->SetLineColor(kGreen);
  lHUp   ->SetLineColor(kRed);
  lHDown ->SetLineColor(kRed);
//.........这里部分代码省略.........