C++ ModelConfig类代码示例

void MultiNetwork::init(const ModelConfig& config,
                        ParamInitCallback callback,
                        const std::vector<ParameterType>& parameterTypes,
                        bool useGpu) {
  CHECK_GT(config.sub_models_size(), 1) << "sub_models_size should GT 1";
  // check submodel[0] is root
  CHECK_EQ("root", config.sub_models(0).name())
      << "sub_models(0) should be root";
  // ignore root
  subNetworks_.resize(config.sub_models_size() - 1);
  // base class
  NeuralNetwork::init(config, callback, parameterTypes, useGpu);
  // sub networks
  for (int i = 1; i < config.sub_models_size(); ++i) {
    std::string subModelName = config.sub_models(i).name();
    if (FLAGS_parallel_nn) {
      subNetworks_[i - 1] = std::unique_ptr<ParallelNeuralNetwork>(
          new ParallelNeuralNetwork(subModelName, this));
    } else {
      subNetworks_[i - 1] = std::unique_ptr<NeuralNetwork>(
          NeuralNetwork::newNeuralNetwork(subModelName, this));
    }
    subNetworks_[i - 1]->init(config);
  }
}

void upper_limit_Bayesian_MCMC(Model* model,double confidence,int Niters){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating upper limit with the MCMC method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  RooWorkspace* wspace = new RooWorkspace("wspace");
  ModelConfig* modelConfig = new ModelConfig("bayes");
  modelConfig->SetWorkspace(*wspace);
  modelConfig->SetPdf(*model->get_complete_likelihood());
  modelConfig->SetPriorPdf(*model->get_POI_prior());
  modelConfig->SetParametersOfInterest(*model->get_POI_set());
  //  modelConfig->SetNuisanceParameters();
  modelConfig->SetNuisanceParameters(*model->get_nuisance_set());


  //configure the calculator
  //model->Print();

  cout<<" POI size "<<model->get_POI_set()->getSize()<<endl; 
  RooRealVar* firstPOI = (RooRealVar*) modelConfig->GetParametersOfInterest()->first();


  MCMCCalculator mcmccalc(*model->get_data(), *modelConfig );
  mcmccalc.SetTestSize(1-confidence);
  mcmccalc.SetLeftSideTailFraction(0);
  mcmccalc.SetNumIters(Niters);
  MCMCInterval* interval = mcmccalc.GetInterval();
  double ul = interval->UpperLimit(*firstPOI);
  cout<<"UpperLimit: "<<ul<<endl;

}

//____________________________________
void DoHypothesisTest(RooWorkspace* wks){


  // Use a RooStats ProfileLikleihoodCalculator to do the hypothesis test.
  ModelConfig model;
  model.SetWorkspace(*wks);
  model.SetPdf("model");

  //plc.SetData("data");

  ProfileLikelihoodCalculator plc;
  plc.SetData( *(wks->data("data") ));

  // here we explicitly set the value of the parameters for the null.
  // We want no signal contribution, eg. mu = 0
  RooRealVar* mu = wks->var("mu");
//   RooArgSet* nullParams = new RooArgSet("nullParams");
//   nullParams->addClone(*mu);
  RooArgSet poi(*mu);
  RooArgSet * nullParams = (RooArgSet*) poi.snapshot();
  nullParams->setRealValue("mu",0);


  //plc.SetNullParameters(*nullParams);
  plc.SetModel(model);
  // NOTE: using snapshot will import nullparams
  // in the WS and merge with existing "mu"
  // model.SetSnapshot(*nullParams);

  //use instead setNuisanceParameters
  plc.SetNullParameters( *nullParams);



  // We get a HypoTestResult out of the calculator, and we can query it.
  HypoTestResult* htr = plc.GetHypoTest();
  cout << "-------------------------------------------------" << endl;
  cout << "The p-value for the null is " << htr->NullPValue() << endl;
  cout << "Corresponding to a signifcance of " << htr->Significance() << endl;
  cout << "-------------------------------------------------\n\n" << endl;


}

void MCStudy()
{
    
    string fileName = "h4l_Template_2012_lowmass_unconstrained_new_v00.root";
    
    const char * wsName      = "combined";    //"w";
    const char * modelSBName = "ModelConfig"; //"modelConfig";
    const char * dataName    = dataname.c_str();
    
    TFile* file = TFile::Open(fileName);

    // get the workspace out of file
    RooWorkspace* w = (RooWorkspace*) file->Get(wsName);
    if (!wsName) {
        cout << "workspace not found" << endl;
        return -1.0;
    }
    
    // get the modelConfig out of the file
    ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);
    
    // get the data out of file
    //RooAbsData* data = w->data(dataName);
    
    // get toy MC data
    TFile *fhist = new TFile("../ToyMC/toyMC_ggF125.root");
    TH1F *hsum = (TH1F*) fhist->Get("sumtoy1");
    
    RooRealVar m4l("m4l","m4l",120,130);
    RooDataHist* data = new RooDataHist("data","data",x,hsum);
    
    // get pdf
    RooAbsPdf* model = sbModel->GetPdf();
    
    RooPlot* xframe = m4l.frame();
    data->plotOn(xframe);
    model->fitTo(*data);
    model->plotOn(xframe,LineColor(kRed));
    
    TCanvas *c = new TCanvas("c","c");
    xframe->Draw();
    
}

void upper_limit_Bayesian(Model* model,double confidence){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating upper limit with the Bayesian method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  RooWorkspace* wspace = new RooWorkspace("wspace");
  ModelConfig* modelConfig = new ModelConfig("bayes");
  modelConfig->SetWorkspace(*wspace);
  modelConfig->SetPdf(*model->get_complete_likelihood());
  modelConfig->SetPriorPdf(*model->get_POI_prior());
  modelConfig->SetParametersOfInterest(*model->get_POI_set());
  //modelConfig->SetNuisanceParameters(*model->get_nuisance_set());



  //configure the calculator
  //model->Print();

  cout<<" POI size "<<model->get_POI_set()->getSize()<<endl; 
  BayesianCalculator bcalc(*model->get_data(), *modelConfig);
  //BayesianCalculator bcalc(*model->get_data(),*model->get_complete_likelihood(),*model->get_POI_set(),*model->get_POI_prior(),model->get_nuisance_set());
  //BayesianCalculator bcalc(*model->get_data(),*model->get_complete_likelihood(),*model->get_POI_set(),*model->get_POI_prior(),0);
 

  bcalc.SetLeftSideTailFraction(0); //for upper limit

  //get the interval
  bcalc.SetConfidenceLevel(confidence);
  cout<<"Calculating"<<endl;
  SimpleInterval* interval = bcalc.GetInterval();
  double ul=interval->UpperLimit();
  std::cout <<confidence <<"% CL upper limit: "<< ul<<endl;

  TCanvas *c1=new TCanvas;
  bcalc.SetScanOfPosterior(100);
  RooPlot * plot = bcalc.GetPosteriorPlot();
  plot->Draw(); 
  c1->SaveAs("bayesian_PosteriorPlot.png");

}

double upper_limit_FC(Model* model,double confidence){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating upper limit with the FC method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  RooWorkspace* wspace = new RooWorkspace("wspace");

  ModelConfig* modelConfig = new ModelConfig("FC");
  modelConfig->SetWorkspace(*wspace);
  modelConfig->SetPdf(*model->get_complete_likelihood());
  modelConfig->SetPriorPdf(*model->get_POI_prior());
  modelConfig->SetParametersOfInterest(*model->get_POI_set());
  //modelConfig->SetParametersOfInterest(*wspace->set("poi"));
  //modelConfig->SetNuisanceParameters(*wspace->set("nuis"));
  // modelConfig->SetNuisanceParameters();
  modelConfig->SetNuisanceParameters(*model->get_nuisance_set());

  RooDataSet* data = model->get_data();
  RooArgSet* poi= model->get_POI_set();
  //configure the calculator
  //model->Print();



  cout<<" POI size "<<model->get_POI_set()->getSize()<<endl; 

  // use FeldmaCousins (takes ~20 min)  
  FeldmanCousins fc(*data, *modelConfig);
  fc.SetConfidenceLevel(0.95); //0.9 central limit=0.95 upper limit
  //fc.SetTestSize(.1); // set size of test
  //number counting: dataset always has 1 entry with N events observed
  fc.FluctuateNumDataEntries(false); 
  fc.UseAdaptiveSampling(true);

  fc.SetNBins(200);
  PointSetInterval* fcInt = NULL;
  //ConfInterval* interval = 0;


  RooRealVar* firstPOI = (RooRealVar*) modelConfig->GetParametersOfInterest()->first();

  //  if(doFeldmanCousins){ // takes 7 minutes
  fcInt = (PointSetInterval*) fc.GetInterval(); // fix cast
  //xs}
  //interval = (PointSetInterval*) fc.GetInterval();

  cout<<" ["<<fcInt->LowerLimit( *firstPOI ) << ", " <<
  fcInt->UpperLimit( *firstPOI ) << "]" << endl;
  cout<<" ["<<fcInt->LowerLimit( *firstPOI ) << ", " <<
  fcInt->UpperLimit( *firstPOI ) << "]" << endl;
    
  double ul=fcInt->UpperLimit( *firstPOI );
  //double ul=interval->UpperLimit( *firstPOI );
  return ul;
}

double TwoBody::printMcmcUpperLimit( double peak, ModelConfig &_mc, std::string filename ){
  //
  // print out the upper limit on the first Parameter of Interest
  //

  std::string _legend = "[TwoBody::printMcmcUpperLimit]: ";

  char buf[1024];

  double _limit = numeric_limits<double>::max();

  if (mcInt){
    //mc.SetWorkspace(*ws);
    RooRealVar * firstPOI = (RooRealVar*) _mc.GetParametersOfInterest()->first();
    _limit = mcInt->UpperLimit(*firstPOI);
    std::cout << "\n95% upper limit on " << firstPOI->GetName() << " is : " <<
      _limit << endl;
    
    if (bMcmcConverged){
      sprintf(buf, "%7.1f   %7.6f", peak, _limit);
    }
    else{
      sprintf(buf, "# %7.1f   %7.6f  # MCMC did not converge", peak, _limit);
    }

  }
  else{
    sprintf(buf, "# MCMC did not converge");
  }

  if (filename.size()!=0){
    
    std::ofstream aFile;
    
    // append to file if exists
    aFile.open(filename.c_str(), std::ios_base::app);
    
    aFile << buf << std::endl;
    
    // close outfile here so it is safe even if subsequent iterations crash
    aFile.close();
    
  }
    
  return _limit;
}

Double_t TwoBody::GetPoiUpper(std::string channel, Double_t peak, ModelConfig &_mc){
  //
  // Estimate a good value for the upper boundary of the range of POI
  //

  std::string legend = "[TwoBody::GetPoiUpper]: ";

  double _range = -1.0;

  std::cout << legend << "doing a rough estimate of the POI range" << std::endl;
  // adding auto-channel option (multi) while
  // preserving backwards compatibility
  // if single channel
  _range = GetPoiUpperSimple(channel, peak);
  
  std::cout << legend << "crude estimate for poi range (x3): "
	    << 3.0*_range << std::endl;
  std::cout << legend 
	    << "this will be used if the profile likelihood ratio estimate fails"
	    << std::endl;
  std::cout << legend 
	    << "will try to estimate POI range better with profile likelihood ratio limit now"
	    << std::endl;

  Double_t result = 0.1;

  // estimate limit with profile likelihood ratio and
  // set the range to 3 times the limit
  
  // query intervals
  RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
  RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
  RooRealVar * _poi = (RooRealVar *)_mc.GetParametersOfInterest()->first();
  double upper_limit = GetPlrInterval(0.95, _mc)->UpperLimit( *_poi );
  RooMsgService::instance().setGlobalKillBelow(msglevel);
  
  // safety in case upper limit == 0
  if (upper_limit<std::numeric_limits<double>::min()){
    upper_limit = _range;
  }

  result = 3.0*upper_limit;

  return result;
}

void OneSidedFrequentistUpperLimitWithBands_intermediate(const char* infile = "",
					    const char* workspaceName = "combined",
					    const char* modelConfigName = "ModelConfig",
					    const char* dataName = "obsData"){


  double confidenceLevel=0.95;
  // degrade/improve number of pseudo-experiments used to define the confidence belt.  
  // value of 1 corresponds to default number of toys in the tail, which is 50/(1-confidenceLevel)
  double additionalToysFac = 1.;  
  int nPointsToScan = 30; // number of steps in the parameter of interest 
  int nToyMC = 100; // number of toys used to define the expected limit and band

  TStopwatch t;
  t.Start();
  /////////////////////////////////////////////////////////////
  // First part is just to access a user-defined file 
  // or create the standard example file if it doesn't exist
  ////////////////////////////////////////////////////////////
  const char* filename = "";
  if (!strcmp(infile,""))
    filename = "results/example_combined_GaussExample_model.root";
  else
    filename = infile;
  // Check if example input file exists
  TFile *file = TFile::Open(filename);

  // if input file was specified byt not found, quit
  if(!file && strcmp(infile,"")){
    cout <<"file not found" << endl;
    return;
  } 

  // if default file not found, try to create it
  if(!file ){
    // Normally this would be run on the command line
    cout <<"will run standard hist2workspace example"<<endl;
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;
  }

  // now try to access the file again
  file = TFile::Open(filename);
  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

  
  /////////////////////////////////////////////////////////////
  // Now get the data and workspace
  ////////////////////////////////////////////////////////////

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  // get the modelConfig out of the file
  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // make sure ingredients are found
  if(!data || !mc){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }

  cout << "Found data and ModelConfig:" <<endl;
  mc->Print();

  /////////////////////////////////////////////////////////////
  // Now get the POI for convenience
  // you may want to adjust the range of your POI
  ////////////////////////////////////////////////////////////
  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  //  firstPOI->setMin(0);
  //  firstPOI->setMax(10);

  /////////////////////////////////////////////
  // create and use the FeldmanCousins tool
  // to find and plot the 95% confidence interval
  // on the parameter of interest as specified
  // in the model config
  // REMEMBER, we will change the test statistic
  // so this is NOT a Feldman-Cousins interval
  FeldmanCousins fc(*data,*mc);
  fc.SetConfidenceLevel(confidenceLevel); 
  fc.AdditionalNToysFactor(additionalToysFac); // improve sampling that defines confidence belt
  //  fc.UseAdaptiveSampling(true); // speed it up a bit, but don't use for expectd limits
  fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
//.........这里部分代码省略.........

void exercise_3() {
  //Open the rootfile and get the workspace from the exercise_0
  TFile fIn("exercise_0.root");
  fIn.cd();
  RooWorkspace *w = (RooWorkspace*)fIn.Get("w");

  //You can set constant parameters that are known
  //If you leave them floating, the fit procedure will determine their uncertainty
  w->var("mean")->setConstant(kFALSE); //don't fix the mean, it's what we want to know the interval for!
  w->var("sigma")->setConstant(kTRUE);
  w->var("tau")->setConstant(kTRUE);
  w->var("Nsig")->setConstant(kTRUE);
  w->var("Nbkg")->setConstant(kTRUE);

  //Set the RooModelConfig and let it know what the content of the workspace is about
  ModelConfig model;
  model.SetWorkspace(*w);
  model.SetPdf("PDFtot");

  //Let the model know what is the parameter of interest
  RooRealVar* mean = w->var("mean");
  mean->setRange(120., 130.);   //this is mostly for plotting reasons
  RooArgSet poi(*mean);

  // set confidence level
  double confidenceLevel = 0.68;

  //Build the profile likelihood calculator
  ProfileLikelihoodCalculator plc; 
  plc.SetData(*(w->data("PDFtotData"))); 
  plc.SetModel(model);
  plc.SetParameters(poi);
  plc.SetConfidenceLevel(confidenceLevel);

  //Get the interval
  LikelihoodInterval* plInt = plc.GetInterval();

  //Now let's do the same for the Bayesian Calculator
  //Now we also need to specify a prior in the ModelConfig
  //To be quicker, we'll use the PDF factory facility of RooWorkspace
  //NB!! For simplicity, we are using a flat prior, but this doesn't mean it's the best choice!
  w->factory("Uniform::prior(mean)");
  model.SetPriorPdf(*w->pdf("prior"));

  //Construct the bayesian calculator
  BayesianCalculator bc(*(w->data("PDFtotData")), model);
  bc.SetConfidenceLevel(confidenceLevel);
  bc.SetParameters(poi);
  SimpleInterval* bcInt = bc.GetInterval();

  // Let's make a plot
  TCanvas dataCanvas("dataCanvas");
  dataCanvas.Divide(2,1);
  dataCanvas.cd(1);

  LikelihoodIntervalPlot plotInt((LikelihoodInterval*)plInt);
  plotInt.SetTitle("Profile Likelihood Ratio and Posterior for mH");
  plotInt.SetMaximum(3.);
  plotInt.Draw();

  dataCanvas.cd(2);
  RooPlot *bcPlot = bc.GetPosteriorPlot();
  bcPlot->Draw();

  dataCanvas.SaveAs("exercise_3.gif");

  //Now print the interval for mH for the two methods
  cout << "PLC interval is [" << plInt->LowerLimit(*mean) << ", " << 
    plInt->UpperLimit(*mean) << "]" << endl;

  cout << "Bayesian interval is [" << bcInt->LowerLimit() << ", " << 
    bcInt->UpperLimit() << "]" << endl;

}

// implementation
void TwoBinInstructional( void ){
  
  // let's time this example
  TStopwatch t;
  t.Start();

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

  // make model
  RooWorkspace * pWs = new RooWorkspace("ws");

  // derived from data
  pWs->factory("xsec[0.2,0,2]"); // POI
  pWs->factory("bg_b[10,0,50]");    // data driven nuisance

  // predefined nuisances
  pWs->factory("lumi[100,0,1000]");
  pWs->factory("eff_a[0.2,0,1]");
  pWs->factory("eff_b[0.05,0,1]");
  pWs->factory("tau[0,1]");
  pWs->factory("xsec_bg_a[0.05]"); // constant
  pWs->var("xsec_bg_a")->setConstant(1);

  // channel a (signal): lumi*xsec*eff_a + lumi*bg_a + tau*bg_b
  pWs->factory("prod::sig_a(lumi,xsec,eff_a)");
  pWs->factory("prod::bg_a(lumi,xsec_bg_a)");
  pWs->factory("prod::tau_bg_b(tau, bg_b)");
  pWs->factory("Poisson::pdf_a(na[14,0,100],sum::mu_a(sig_a,bg_a,tau_bg_b))");

  // channel b (control): lumi*xsec*eff_b + bg_b
  pWs->factory("prod::sig_b(lumi,xsec,eff_b)");
  pWs->factory("Poisson::pdf_b(nb[11,0,100],sum::mu_b(sig_b,bg_b))");

  // nuisance constraint terms (systematics)
  pWs->factory("Lognormal::l_lumi(lumi,nom_lumi[100,0,1000],sum::kappa_lumi(1,d_lumi[0.1]))");
  pWs->factory("Lognormal::l_eff_a(eff_a,nom_eff_a[0.20,0,1],sum::kappa_eff_a(1,d_eff_a[0.05]))");
  pWs->factory("Lognormal::l_eff_b(eff_b,nom_eff_b[0.05,0,1],sum::kappa_eff_b(1,d_eff_b[0.05]))");
  pWs->factory("Lognormal::l_tau(tau,nom_tau[0.50,0,1],sum::kappa_tau(1,d_tau[0.05]))");
  //pWs->factory("Lognormal::l_bg_a(bg_a,nom_bg_a[0.05,0,1],sum::kappa_bg_a(1,d_bg_a[0.10]))");

  // complete model PDF
  pWs->factory("PROD::model(pdf_a,pdf_b,l_lumi,l_eff_a,l_eff_b,l_tau)");

  // Now create sets of variables. Note that we could use the factory to
  // create sets but in that case many of the sets would be duplicated
  // when the ModelConfig objects are imported into the workspace. So,
  // we create the sets outside the workspace, and only the needed ones
  // will be automatically imported by ModelConfigs

  // observables
  RooArgSet obs(*pWs->var("na"), *pWs->var("nb"), "obs");

  // global observables
  RooArgSet globalObs(*pWs->var("nom_lumi"), *pWs->var("nom_eff_a"), *pWs->var("nom_eff_b"), 
		      *pWs->var("nom_tau"),
		      "global_obs");

  // parameters of interest
  RooArgSet poi(*pWs->var("xsec"), "poi");

  // nuisance parameters
  RooArgSet nuis(*pWs->var("lumi"), *pWs->var("eff_a"), *pWs->var("eff_b"), *pWs->var("tau"), "nuis");

  // priors (for Bayesian calculation)
  pWs->factory("Uniform::prior_xsec(xsec)"); // for parameter of interest
  pWs->factory("Uniform::prior_bg_b(bg_b)"); // for data driven nuisance parameter
  pWs->factory("PROD::prior(prior_xsec,prior_bg_b)"); // total prior

  // create data
  pWs->var("na")->setVal(14);
  pWs->var("nb")->setVal(11);
  RooDataSet * pData = new RooDataSet("data","",obs);
  pData->add(obs);
  pWs->import(*pData);
  //pData->Print();

  // signal+background model
  ModelConfig * pSbModel = new ModelConfig("SbModel");
  pSbModel->SetWorkspace(*pWs);
  pSbModel->SetPdf(*pWs->pdf("model"));
  pSbModel->SetPriorPdf(*pWs->pdf("prior"));
  pSbModel->SetParametersOfInterest(poi);
  pSbModel->SetNuisanceParameters(nuis);
  pSbModel->SetObservables(obs);
  pSbModel->SetGlobalObservables(globalObs);

  // set all but obs, poi and nuisance to const
  SetConstants(pWs, pSbModel);
  pWs->import(*pSbModel);


  // background-only model
  // use the same PDF as s+b, with xsec=0
  // POI value under the background hypothesis
  Double_t poiValueForBModel = 0.0;
  ModelConfig* pBModel = new ModelConfig(*(RooStats::ModelConfig *)pWs->obj("SbModel"));
  pBModel->SetName("BModel");
  pBModel->SetWorkspace(*pWs);
//.........这里部分代码省略.........

void StandardFeldmanCousinsDemo(const char* infile = "",
                                const char* workspaceName = "combined",
                                const char* modelConfigName = "ModelConfig",
                                const char* dataName = "obsData"){

   // -------------------------------------------------------
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
   const char* filename = "";
   if (!strcmp(infile,"")) {
      filename = "results/example_combined_GaussExample_model.root";
      bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
      // if file does not exists generate with histfactory
      if (!fileExist) {
#ifdef _WIN32
         cout << "HistFactory file cannot be generated on Windows - exit" << endl;
         return;
#endif
         // Normally this would be run on the command line
         cout <<"will run standard hist2workspace example"<<endl;
         gROOT->ProcessLine(".! prepareHistFactory .");
         gROOT->ProcessLine(".! hist2workspace config/example.xml");
         cout <<"\n\n---------------------"<<endl;
         cout <<"Done creating example input"<<endl;
         cout <<"---------------------\n\n"<<endl;
      }

   }
   else
      filename = infile;

   // Try to open the file
   TFile *file = TFile::Open(filename);

   // if input file was specified byt not found, quit
   if(!file ){
      cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
      return;
   }


   // -------------------------------------------------------
   // Tutorial starts here
   // -------------------------------------------------------

   // get the workspace out of the file
   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
   if(!w){
      cout <<"workspace not found" << endl;
      return;
   }

   // get the modelConfig out of the file
   ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

   // get the modelConfig out of the file
   RooAbsData* data = w->data(dataName);

   // make sure ingredients are found
   if(!data || !mc){
      w->Print();
      cout << "data or ModelConfig was not found" <<endl;
      return;
   }

   // -------------------------------------------------------
   // create and use the FeldmanCousins tool
   // to find and plot the 95% confidence interval
   // on the parameter of interest as specified
   // in the model config
   FeldmanCousins fc(*data,*mc);
   fc.SetConfidenceLevel(0.95); // 95% interval
   //fc.AdditionalNToysFactor(0.1); // to speed up the result
   fc.UseAdaptiveSampling(true); // speed it up a bit
   fc.SetNBins(10); // set how many points per parameter of interest to scan
   fc.CreateConfBelt(true); // save the information in the belt for plotting

   // Since this tool needs to throw toy MC the PDF needs to be
   // extended or the tool needs to know how many entries in a dataset
   // per pseudo experiment.
   // In the 'number counting form' where the entries in the dataset
   // are counts, and not values of discriminating variables, the
   // datasets typically only have one entry and the PDF is not
   // extended.
   if(!mc->GetPdf()->canBeExtended()){
      if(data->numEntries()==1)
         fc.FluctuateNumDataEntries(false);
      else
         cout <<"Not sure what to do about this model" <<endl;
   }

   // We can use PROOF to speed things along in parallel
   //  ProofConfig pc(*w, 1, "workers=4", kFALSE);
   //  ToyMCSampler*  toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
   //  toymcsampler->SetProofConfig(&pc); // enable proof


   // Now get the interval
   PointSetInterval* interval = fc.GetInterval();
   ConfidenceBelt* belt = fc.GetConfidenceBelt();
//.........这里部分代码省略.........

void OneSidedFrequentistUpperLimitWithBands(const char* infile = "",
                                            const char* workspaceName = "combined",
                                            const char* modelConfigName = "ModelConfig",
                                            const char* dataName = "obsData") {



   double confidenceLevel=0.95;
   int nPointsToScan = 20;
   int nToyMC = 200;

   // -------------------------------------------------------
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
   const char* filename = "";
   if (!strcmp(infile,"")) {
      filename = "results/example_combined_GaussExample_model.root";
      bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
      // if file does not exists generate with histfactory
      if (!fileExist) {
#ifdef _WIN32
         cout << "HistFactory file cannot be generated on Windows - exit" << endl;
         return;
#endif
         // Normally this would be run on the command line
         cout <<"will run standard hist2workspace example"<<endl;
         gROOT->ProcessLine(".! prepareHistFactory .");
         gROOT->ProcessLine(".! hist2workspace config/example.xml");
         cout <<"\n\n---------------------"<<endl;
         cout <<"Done creating example input"<<endl;
         cout <<"---------------------\n\n"<<endl;
      }

   }
   else
      filename = infile;

   // Try to open the file
   TFile *file = TFile::Open(filename);

   // if input file was specified byt not found, quit
   if(!file ){
      cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
      return;
   }


   // -------------------------------------------------------
   // Now get the data and workspace

   // get the workspace out of the file
   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
   if(!w){
      cout <<"workspace not found" << endl;
      return;
   }

   // get the modelConfig out of the file
   ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

   // get the modelConfig out of the file
   RooAbsData* data = w->data(dataName);

   // make sure ingredients are found
   if(!data || !mc){
      w->Print();
      cout << "data or ModelConfig was not found" <<endl;
      return;
   }

   // -------------------------------------------------------
   // Now get the POI for convenience
   // you may want to adjust the range of your POI

   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   /*  firstPOI->setMin(0);*/
   /*  firstPOI->setMax(10);*/

   // --------------------------------------------
   // Create and use the FeldmanCousins tool
   // to find and plot the 95% confidence interval
   // on the parameter of interest as specified
   // in the model config
   // REMEMBER, we will change the test statistic
   // so this is NOT a Feldman-Cousins interval
   FeldmanCousins fc(*data,*mc);
   fc.SetConfidenceLevel(confidenceLevel);
   /*  fc.AdditionalNToysFactor(0.25); // degrade/improve sampling that defines confidence belt*/
   /*  fc.UseAdaptiveSampling(true); // speed it up a bit, don't use for expected limits*/
   fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
   fc.CreateConfBelt(true); // save the information in the belt for plotting

   // -------------------------------------------------------
   // Feldman-Cousins is a unified limit by definition
   // but the tool takes care of a few things for us like which values
   // of the nuisance parameters should be used to generate toys.
   // so let's just change the test statistic and realize this is
   // no longer "Feldman-Cousins" but is a fully frequentist Neyman-Construction.
   /*  ProfileLikelihoodTestStatModified onesided(*mc->GetPdf());*/
   /*  fc.GetTestStatSampler()->SetTestStatistic(&onesided);*/
//.........这里部分代码省略.........

void IntervalExamples()
{

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


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

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

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

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

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

   // set confidence level
   double confidenceLevel = 0.95;

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

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

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

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

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


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

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

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

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

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

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

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

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

void PlotAll(TString wsname)
{
	char* binLabels[19] = {"60","70","80","90","100","110","120","130","140","150","160","170","180","190","200","250","300","400","1000"};	

	//get the stuff from the workspace:
	
	TFile* file=TFile::Open(wsname);
	RooWorkspace* ws = (RooWorkspace*)file->Get("combined");
	ModelConfig  *mc = (ModelConfig*)ws->obj("ModelConfig");
	RooAbsData   *data = ws->data("obsData");
	RooSimultaneous* simPdf=(RooSimultaneous*)(mc->GetPdf());
	RooAbsReal* nll=simPdf->createNLL(*data);

	// FPT 0 **************************************	
	// EM channel
	
	RooCategory* chanCat = (RooCategory*) (&simPdf->indexCat());
        TIterator* iterat = chanCat->typeIterator() ;
        RooCatType* ttype = (RooCatType*)iterat->Next();

	RooAbsPdf  *pdf_stateEM  = simPdf->getPdf(ttype->GetName()) ;
	RooArgSet  *obstmpEM  = pdf_stateEM->getObservables( *mc->GetObservables() ) ;
	
	// get EM data
       	RooAbsData *dataEM = data->reduce(Form("%s==%s::%s",chanCat->GetName(),chanCat->GetName(),ttype->GetName()));
		
	RooRealVar *obsEM     = ((RooRealVar*) obstmpEM->first());
	TString chanName1(ttype->GetName());

	// create data histogram
	TH1* hdataEM = dataEM->createHistogram("Data "+chanName1,*obsEM);
	// set errors to gaussian
        for (int ib=0 ; ib<hdataEM->GetNbinsX()+1 ; ib++) hdataEM->SetBinError(ib, sqrt(hdataEM->GetBinContent(ib)));

	double EMnorm = pdf_stateEM->expectedEvents(*obsEM);
	
	//****************************
	// ME channel
	ttype = (RooCatType*)iterat->Next();
	RooAbsPdf* pdf_stateME  = simPdf->getPdf(ttype->GetName()) ;
        RooArgSet* obstmpME  = pdf_stateME->getObservables( *mc->GetObservables() ) ;

	// get ME data
	RooAbsData *dataME = data->reduce(Form("%s==%s::%s",chanCat->GetName(),chanCat->GetName(),ttype->GetName()));	
	RooRealVar* obsME = ((RooRealVar*) obstmpME->first());
	TString chanName2(ttype->GetName());

        // create data histogram
        TH1* hdataME = dataME->createHistogram("Data "+chanName2,*obsME);
        // set errors to gaussian
        for (int ib=0 ; ib<hdataME->GetNbinsX()+1 ; ib++) hdataME->SetBinError(ib, sqrt(hdataME->GetBinContent(ib)));
        
	
	// get initial BG histogram
	//TH1* h_initial_BG_EM = pdf_stateEM->createHistogram("initial_BG_EM",*obsEM);
	//TH1* h_initial_BG_ME = pdf_stateME->createHistogram("initial_BG_ME",*obsME);
	
	double MEnorm = pdf_stateME->expectedEvents(*obsME);
	cout << "EM expected events = " << EMnorm << ", ME expected events = " << MEnorm << "." << endl;
	//h_initial_BG_EM->Scale(EMnorm);
	//h_initial_BG_ME->Scale(MEnorm);	

	// get initial gammas
	int nbins = hdataEM->GetNbinsX();
        double InitGamma[nbins];
        for (int i=0; i<nbins; i++)
        {
               	TString varname = "gamma_B0_l1pt0_bin_"+NumberToString(i);
               	InitGamma[i] = ws->var(varname)->getVal();
               	cout << "initial gamma"+NumberToString(i)+" = " << InitGamma[i] << endl;
        }
        double InitFpt = ws->var("fl1pt_l1pt0")->getVal();
        cout << "initial fpt_l1pt0 = " << InitFpt <<  endl;


	// DO THE GLOBAL FIT
	
	minimize(nll);	
       
	// get final BG histograms
	TH1* h_final_BG_EM = pdf_stateEM->createHistogram("final_BG_EM",*obsEM);
        TH1* h_final_BG_ME = pdf_stateME->createHistogram("final_BG_ME",*obsME); 
	h_final_BG_EM->Scale(EMnorm);
	h_final_BG_ME->Scale(MEnorm);
	
	// uncertainty bands
	TH1D* BuncertaintyEM = new TH1D("BuncertaintyEM","BuncertaintyEM",nbins,0,nbins);
	TH1D* BuncertaintyME = new TH1D("BuncertaintyME","BuncertaintyME",nbins,0,nbins);
	for (int i=1; i<=nbins; i++){
		double sigbEM = h_final_BG_EM->GetBinError(i);
		double bEM = h_final_BG_EM->GetBinContent(i);
		BuncertaintyEM->SetBinError(i,sigbEM); BuncertaintyEM->SetBinContent(i,bEM);
		double sigbME = h_final_BG_ME->GetBinError(i);
                double bME = h_final_BG_ME->GetBinContent(i);
                BuncertaintyME->SetBinError(i,sigbME); BuncertaintyME->SetBinContent(i,bME);
	}
	//BuncertaintyEM->SetFillStyle(3004); 
	BuncertaintyEM->SetFillColor(kGreen-9);
	BuncertaintyEM->SetLineColor(kBlack); BuncertaintyEM->SetLineStyle(2);
	//BuncertaintyME->SetFillStyle(3004); 
//.........这里部分代码省略.........