C++ ModelConfig::SetSnapshot方法代码示例

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

    //Set Nuisnaces

    RooArgSet nuis(*newworkspace->var("prime_lumi"), *newworkspace->var("prime_eff"), *newworkspace->var("prime_rho"), *newworkspace->var("bprime"));

    // priors (for Bayesian calculation)
    newworkspace->factory("Uniform::prior_signal(sigma)"); // for parameter of interest
    newworkspace->factory("Uniform::prior_bg_b(bprime)"); // for data driven nuisance parameter
    newworkspace->factory("PROD::prior(prior_signal,prior_bg_b)"); // total prior


    //Observed data is pulled from histogram.
    //TFile *data_file = new TFile(data_file_name);
    TFile *data_file = new TFile(data_file_name);
    TH2D *data_hist = (TH2D *)data_file->Get(data_hist_name_in_file);
    RooDataHist *pData = new RooDataHist("data", "data", obs, data_hist);
    newworkspace->import(*pData);

    // Now, we will draw our data from a RooDataHist.
    /*TFile *data_file = new TFile(data_file_name);
    TTree *data_tree = (TTree *) data_file->Get(data_hist_name_in_file);
    RooDataSet *pData = new RooDataSet("data", "data", data_tree, obs);
    newworkspace->import(*pData);*/


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

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


    // background-only model
    // use the same PDF as s+b, with sig=0
    // POI value under the background hypothesis
    // (We will set the value to 0 later)

    Double_t poiValueForBModel = 0.0;
    ModelConfig* pBModel = new ModelConfig(*(RooStats::ModelConfig *)newworkspace->obj("SbModel"));
    pBModel->SetName("BModel");
    pBModel->SetWorkspace(*newworkspace);
    newworkspace->import(*pBModel);

    // find global maximum with the signal+background model
    // with conditional MLEs for nuisance parameters
    // and save the parameter point snapshot in the Workspace
    //  - safer to keep a default name because some RooStats calculators
    //    will anticipate it
    RooAbsReal * pNll = pSbModel->GetPdf()->createNLL(*pData);
    RooAbsReal * pProfile = pNll->createProfile(RooArgSet());
    pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
    RooArgSet * pPoiAndNuisance = new RooArgSet();
    if(pSbModel->GetNuisanceParameters())
        pPoiAndNuisance->add(*pSbModel->GetNuisanceParameters());
    pPoiAndNuisance->add(*pSbModel->GetParametersOfInterest());
    cout << "\nWill save these parameter points that correspond to the fit to data" << endl;
    pPoiAndNuisance->Print("v");
    pSbModel->SetSnapshot(*pPoiAndNuisance);
    delete pProfile;
    delete pNll;
    delete pPoiAndNuisance;


    // Find a parameter point for generating pseudo-data
    // with the background-only data.
    // Save the parameter point snapshot in the Workspace
    pNll = pBModel->GetPdf()->createNLL(*pData);
    pProfile = pNll->createProfile(poi);
    ((RooRealVar *)poi.first())->setVal(poiValueForBModel);
    pProfile->getVal(); // this will do fit and set nuisance parameters to profiled values
    pPoiAndNuisance = new RooArgSet();
    if(pBModel->GetNuisanceParameters())
        pPoiAndNuisance->add(*pBModel->GetNuisanceParameters());
    pPoiAndNuisance->add(*pBModel->GetParametersOfInterest());
    cout << "\nShould use these parameter points to generate pseudo data for bkg only" << endl;
    pPoiAndNuisance->Print("v");
    pBModel->SetSnapshot(*pPoiAndNuisance);
    delete pProfile;
    delete pNll;
    delete pPoiAndNuisance;

    // save workspace to file
    newworkspace->writeToFile("ws_twobin.root");

    // clean up
    delete newworkspace;
    delete pData;
    delete pSbModel;
    delete pBModel;


} // ----- end of tutorial ----------------------------------------

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

  // 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);
  pWs->import(*pBModel);


  // find global maximum with the signal+background model
  // with conditional MLEs for nuisance parameters
  // and save the parameter point snapshot in the Workspace
  //  - safer to keep a default name because some RooStats calculators
  //    will anticipate it
  RooAbsReal * pNll = pSbModel->GetPdf()->createNLL(*pData);
  RooAbsReal * pProfile = pNll->createProfile(RooArgSet());
  pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
  RooArgSet * pPoiAndNuisance = new RooArgSet();
  if(pSbModel->GetNuisanceParameters())
    pPoiAndNuisance->add(*pSbModel->GetNuisanceParameters());
  pPoiAndNuisance->add(*pSbModel->GetParametersOfInterest());
  cout << "\nWill save these parameter points that correspond to the fit to data" << endl;
  pPoiAndNuisance->Print("v");
  pSbModel->SetSnapshot(*pPoiAndNuisance);
  delete pProfile;
  delete pNll;
  delete pPoiAndNuisance;

  // Find a parameter point for generating pseudo-data
  // with the background-only data.
  // Save the parameter point snapshot in the Workspace
  pNll = pBModel->GetPdf()->createNLL(*pData);
  pProfile = pNll->createProfile(poi);
  ((RooRealVar *)poi.first())->setVal(poiValueForBModel);
  pProfile->getVal(); // this will do fit and set nuisance parameters to profiled values
  pPoiAndNuisance = new RooArgSet();
  if(pBModel->GetNuisanceParameters())
    pPoiAndNuisance->add(*pBModel->GetNuisanceParameters());
  pPoiAndNuisance->add(*pBModel->GetParametersOfInterest());
  cout << "\nShould use these parameter points to generate pseudo data for bkg only" << endl;
  pPoiAndNuisance->Print("v");
  pBModel->SetSnapshot(*pPoiAndNuisance);
  delete pProfile;
  delete pNll;
  delete pPoiAndNuisance;

  // inspect workspace
  pWs->Print();

  // save workspace to file
  pWs->writeToFile("ws_twobin.root");

  // clean up
  delete pWs;
  delete pData;
  delete pSbModel;
  delete pBModel;

} // ----- end of tutorial ----------------------------------------

// internal routine to run the inverter
HypoTestInverterResult *
RooStats::HypoTestInvTool::RunInverter(RooWorkspace * w,
                                       const char * modelSBName, const char * modelBName, 
                                       const char * dataName, int type,  int testStatType, 
                                       bool useCLs, int npoints, double poimin, double poimax, 
                                       int ntoys,
                                       bool useNumberCounting,
                                       const char * nuisPriorName ){

   std::cout << "Running HypoTestInverter on the workspace " << w->GetName() << std::endl;
  
   w->Print();
  
  
   RooAbsData * data = w->data(dataName); 
   if (!data) { 
      Error("StandardHypoTestDemo","Not existing data %s",dataName);
      return 0;
   }
   else 
      std::cout << "Using data set " << dataName << std::endl;
  
   if (mUseVectorStore) { 
      RooAbsData::setDefaultStorageType(RooAbsData::Vector);
      data->convertToVectorStore() ;
   }
  
  
   // get models from WS
   // get the modelConfig out of the file
   ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);
   ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);
  
   if (!sbModel) {
      Error("StandardHypoTestDemo","Not existing ModelConfig %s",modelSBName);
      return 0;
   }
   // check the model 
   if (!sbModel->GetPdf()) { 
      Error("StandardHypoTestDemo","Model %s has no pdf ",modelSBName);
      return 0;
   }
   if (!sbModel->GetParametersOfInterest()) {
      Error("StandardHypoTestDemo","Model %s has no poi ",modelSBName);
      return 0;
   }
   if (!sbModel->GetObservables()) {
      Error("StandardHypoTestInvDemo","Model %s has no observables ",modelSBName);
      return 0;
   }
   if (!sbModel->GetSnapshot() ) { 
      Info("StandardHypoTestInvDemo","Model %s has no snapshot  - make one using model poi",modelSBName);
      sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );
   }
  
   // case of no systematics
   // remove nuisance parameters from model
   if (noSystematics) { 
      const RooArgSet * nuisPar = sbModel->GetNuisanceParameters();
      if (nuisPar && nuisPar->getSize() > 0) { 
         std::cout << "StandardHypoTestInvDemo" << "  -  Switch off all systematics by setting them constant to their initial values" << std::endl;
         RooStats::SetAllConstant(*nuisPar);
      }
      if (bModel) { 
         const RooArgSet * bnuisPar = bModel->GetNuisanceParameters();
         if (bnuisPar) 
            RooStats::SetAllConstant(*bnuisPar);
      }
   }
  
   if (!bModel || bModel == sbModel) {
      Info("StandardHypoTestInvDemo","The background model %s does not exist",modelBName);
      Info("StandardHypoTestInvDemo","Copy it from ModelConfig %s and set POI to zero",modelSBName);
      bModel = (ModelConfig*) sbModel->Clone();
      bModel->SetName(TString(modelSBName)+TString("_with_poi_0"));      
      RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
      if (!var) return 0;
      double oldval = var->getVal();
      var->setVal(0);
      bModel->SetSnapshot( RooArgSet(*var)  );
      var->setVal(oldval);
   }
   else { 
      if (!bModel->GetSnapshot() ) { 
         Info("StandardHypoTestInvDemo","Model %s has no snapshot  - make one using model poi and 0 values ",modelBName);
         RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
         if (var) { 
            double oldval = var->getVal();
            var->setVal(0);
            bModel->SetSnapshot( RooArgSet(*var)  );
            var->setVal(oldval);
         }
         else { 
            Error("StandardHypoTestInvDemo","Model %s has no valid poi",modelBName);
            return 0;
         }         
      }
   }

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

double StandardFrequentistDiscovery(
   const char* infile = "",
   const char* workspaceName = "channel1",
   const char* modelConfigNameSB = "ModelConfig",
   const char* dataName = "obsData",
   int toys = 1000,
   double poiValueForBackground = 0.0,
   double poiValueForSignal = 1.0
) {

   // The workspace contains the model for s+b. The b model is "autogenerated"
   // by copying s+b and setting the one parameter of interest to zero.
   // To keep the script simple, multiple parameters of interest or different
   // functional forms of the b model are not supported.

   // for now, assume there is only one parameter of interest, and these are
   // its values:

   /////////////////////////////////////////////////////////////
   // 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_channel1_GammaExample_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 -1;
#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 -1;
   } 


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

   TStopwatch *mn_t = new TStopwatch;
   mn_t->Start();

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

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

   // get the data 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 -1.0;
   }


   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   firstPOI->setVal(poiValueForSignal);
   mc->SetSnapshot(*mc->GetParametersOfInterest());
   // create null model
   ModelConfig *mcNull = mc->Clone("ModelConfigNull");
   firstPOI->setVal(poiValueForBackground);
   mcNull->SetSnapshot(*(RooArgSet*)mcNull->GetParametersOfInterest()->snapshot());



   // ----------------------------------------------------
   // Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
   // to use simultaneously with ToyMCSampler
   ProfileLikelihoodTestStat* plts =  new ProfileLikelihoodTestStat(*mc->GetPdf());
   plts->SetOneSidedDiscovery(true);
   plts->SetVarName( "q_{0}/2" );
   
//.........这里部分代码省略.........

// internal routine to run the inverter
HypoTestInverterResult *  RunInverter(RooWorkspace * w, const char * modelSBName, const char * modelBName, 
                                      const char * dataName, int type,  int testStatType, 
                                      int npoints, double poimin, double poimax, 
                                      int ntoys, bool useCls ) 
{

   std::cout << "Running HypoTestInverter on the workspace " << w->GetName() << std::endl;

   w->Print();


   RooAbsData * data = w->data(dataName); 
   if (!data) { 
      Error("RA2bHypoTestDemo","Not existing data %s",dataName);
      return 0;
   }
   else 
      std::cout << "Using data set " << dataName << std::endl;

   
   // get models from WS
   // get the modelConfig out of the file
   ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);
   ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);

   if (!sbModel) {
      Error("RA2bHypoTestDemo","Not existing ModelConfig %s",modelSBName);
      return 0;
   }
   // check the model 
   if (!sbModel->GetPdf()) { 
      Error("RA2bHypoTestDemo","Model %s has no pdf ",modelSBName);
      return 0;
   }
   if (!sbModel->GetParametersOfInterest()) {
      Error("RA2bHypoTestDemo","Model %s has no poi ",modelSBName);
      return 0;
   }
   if (!sbModel->GetParametersOfInterest()) {
      Error("RA2bHypoTestInvDemo","Model %s has no poi ",modelSBName);
      return 0;
   }
   if (!sbModel->GetSnapshot() ) { 
      Info("RA2bHypoTestInvDemo","Model %s has no snapshot  - make one using model poi",modelSBName);
      sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );
   }


   if (!bModel || bModel == sbModel) {
      Info("RA2bHypoTestInvDemo","The background model %s does not exist",modelBName);
      Info("RA2bHypoTestInvDemo","Copy it from ModelConfig %s and set POI to zero",modelSBName);
      bModel = (ModelConfig*) sbModel->Clone();
      bModel->SetName(TString(modelSBName)+TString("_with_poi_0"));      
      RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
      if (!var) return 0;
      double oldval = var->getVal();
      var->setVal(0);
      bModel->SetSnapshot( RooArgSet(*var)  );
      var->setVal(oldval);
   }
   else { 
      if (!bModel->GetSnapshot() ) { 
         Info("RA2bHypoTestInvDemo","Model %s has no snapshot  - make one using model poi and 0 values ",modelBName);
         RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
         if (var) { 
            double oldval = var->getVal();
            var->setVal(0);
            bModel->SetSnapshot( RooArgSet(*var)  );
            var->setVal(oldval);
         }
         else { 
            Error("RA2bHypoTestInvDemo","Model %s has no valid poi",modelBName);
            return 0;
         }         
      }
   }


   SimpleLikelihoodRatioTestStat slrts(*sbModel->GetPdf(),*bModel->GetPdf());
   if (sbModel->GetSnapshot()) slrts.SetNullParameters(*sbModel->GetSnapshot());
   if (bModel->GetSnapshot()) slrts.SetAltParameters(*bModel->GetSnapshot());

   // ratio of profile likelihood - need to pass snapshot for the alt
   RatioOfProfiledLikelihoodsTestStat 
      ropl(*sbModel->GetPdf(), *bModel->GetPdf(), bModel->GetSnapshot());
   ropl.SetSubtractMLE(false);
   
   //MyProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   if (testStatType == 3) profll.SetOneSided(1);
   if (optimize) profll.SetReuseNLL(true);

   TestStatistic * testStat = &slrts;
   if (testStatType == 1) testStat = &ropl;
   if (testStatType == 2 || testStatType == 3) testStat = &profll;
  
   
   HypoTestCalculatorGeneric *  hc = 0;
   if (type == 0) hc = new FrequentistCalculator(*data, *bModel, *sbModel);
//.........这里部分代码省略.........

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

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

  // make model
  RooWorkspace* wspace = new RooWorkspace("wspace");

  wspace->factory("Gaussian::sigCons(prime_SigEff[0,-5,5], nom_SigEff[0,-5,5], 1)");
  wspace->factory("expr::SigEff('1.0*pow(1.20,@0)',prime_SigEff)"); // // 1+-20%, 1.20=exp(20%)

  wspace->factory("Poisson::on(non[0,50], sum::splusb(prod::SigUnc(s[0,0,50],SigEff),mainb[8.8,0,50],dilep[0.9,0,20],tau[2.3,0,20],QCD[0.,0,10],MC[0.1,0,4]))");

  wspace->factory("Gaussian::mcCons(prime_rho[0,-5,5], nom_rho[0,-5,5], 1)");
  wspace->factory("expr::rho('1.0*pow(1.39,@0)',prime_rho)"); // // 1+-39%
  wspace->factory("Poisson::off(noff[0,200], prod::rhob(mainb,rho,mu_plus_e[0.74,0.01,10],1.08))");
  wspace->factory("Gaussian::mcCons2(mu_plus_enom[0.74,0.01,4], mu_plus_e, sigmatwo[.05])");

  wspace->factory("Gaussian::dilep_pred(dilep_nom[0.9,0,20], dilep, sigma3[2.2])");
  wspace->factory("Gaussian::tau_pred(tau_nom[2.3,0,20], tau, sigma4[0.5])");
  wspace->factory("Gaussian::QCD_pred(QCD_nom[0.0,0,10], QCD, sigma5[1.0])");
  wspace->factory("Gaussian::MC_pred(MC_nom[0.1,0.01,4], MC, sigma7[0.14])");

  wspace->factory("PROD::model(on,off,mcCons,mcCons2,sigCons,dilep_pred,tau_pred,QCD_pred,MC_pred)");

  RooArgSet obs(*wspace->var("non"), *wspace->var("noff"), *wspace->var("mu_plus_enom"), *wspace->var("dilep_nom"), *wspace->var("tau_nom"), "obs");
  obs.add(*wspace->var("QCD_nom"));  obs.add(*wspace->var("MC_nom"));
  RooArgSet globalObs(*wspace->var("nom_SigEff"), *wspace->var("nom_rho"), "global_obs");
  // fix global observables to their nominal values
  wspace->var("nom_SigEff")->setConstant();
  wspace->var("nom_rho")->setConstant();

  RooArgSet poi(*wspace->var("s"), "poi");
  RooArgSet nuis(*wspace->var("mainb"), *wspace->var("prime_rho"), *wspace->var("prime_SigEff"), *wspace->var("mu_plus_e"), *wspace->var("dilep"), *wspace->var("tau"), "nuis");
  nuis.add(*wspace->var("QCD"));  nuis.add(*wspace->var("MC"));


  wspace->factory("Uniform::prior_poi({s})");
  wspace->factory("Uniform::prior_nuis({mainb,mu_plus_e,dilep,tau,QCD,MC})");
  wspace->factory("PROD::prior(prior_poi,prior_nuis)");

  wspace->var("non")->setVal(8); //observed
  //wspace->var("non")->setVal(12); //expected observation
  wspace->var("noff")->setVal(7); //observed events in control region
  wspace->var("mu_plus_enom")->setVal(0.74);
  wspace->var("dilep_nom")->setVal(0.9);
  wspace->var("tau_nom")->setVal(2.3);
  wspace->var("QCD")->setVal(0.0);
  wspace->var("MC")->setVal(0.1);


  RooDataSet * data = new RooDataSet("data","",obs);
  data->add(obs);
  wspace->import(*data);


  /////////////////////////////////////////////////////
  // Now the statistical tests
  // model config
  ModelConfig* pSbModel = new ModelConfig("SbModel");
  pSbModel->SetWorkspace(*wspace);
  pSbModel->SetPdf(*wspace->pdf("model"));
  pSbModel->SetPriorPdf(*wspace->pdf("prior"));
  pSbModel->SetParametersOfInterest(poi);
  pSbModel->SetNuisanceParameters(nuis);
  pSbModel->SetObservables(obs);
  pSbModel->SetGlobalObservables(globalObs);
  wspace->import(*pSbModel);

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


  Double_t poiValueForBModel = 0.0;
  ModelConfig* pBModel = new ModelConfig(*(RooStats::ModelConfig *)wspace->obj("SbModel"));
  pBModel->SetName("BModel");
  pBModel->SetWorkspace(*wspace);
  wspace->import(*pBModel);


  RooAbsReal * pNll = pSbModel->GetPdf()->createNLL(*data);
  RooAbsReal * pProfile = pNll->createProfile(RooArgSet());
  pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
  RooArgSet * pPoiAndNuisance = new RooArgSet();
  //if(pSbModel->GetNuisanceParameters())
  //  pPoiAndNuisance->add(*pSbModel->GetNuisanceParameters());
  pPoiAndNuisance->add(*pSbModel->GetParametersOfInterest());
  cout << "\nWill save these parameter points that correspond to the fit to data" << endl;
  pPoiAndNuisance->Print("v");
  pSbModel->SetSnapshot(*pPoiAndNuisance);
  delete pProfile;
  delete pNll;
  delete pPoiAndNuisance;


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

//.........这里部分代码省略.........
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }
  // make b model
  ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);


   // case of no systematics
   // remove nuisance parameters from model
   if (noSystematics) { 
      const RooArgSet * nuisPar = sbModel->GetNuisanceParameters();
      if (nuisPar && nuisPar->getSize() > 0) { 
         std::cout << "StandardHypoTestInvDemo" << "  -  Switch off all systematics by setting them constant to their initial values" << std::endl;
         RooStats::SetAllConstant(*nuisPar);
      }
      if (bModel) { 
         const RooArgSet * bnuisPar = bModel->GetNuisanceParameters();
         if (bnuisPar) 
            RooStats::SetAllConstant(*bnuisPar);
      }
   }


  if (!bModel ) {
      Info("StandardHypoTestInvDemo","The background model %s does not exist",modelBName);
      Info("StandardHypoTestInvDemo","Copy it from ModelConfig %s and set POI to zero",modelSBName);
      bModel = (ModelConfig*) sbModel->Clone();
      bModel->SetName(TString(modelSBName)+TString("B_only"));      
      RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
      if (!var) return;
      double oldval = var->getVal();
      var->setVal(0);
      //bModel->SetSnapshot( RooArgSet(*var, *w->var("lumi"))  );
      bModel->SetSnapshot( RooArgSet(*var)  );
      var->setVal(oldval);
  }
  
   if (!sbModel->GetSnapshot() || poiValue > 0) { 
      Info("StandardHypoTestDemo","Model %s has no snapshot  - make one using model poi",modelSBName);
      RooRealVar * var = dynamic_cast<RooRealVar*>(sbModel->GetParametersOfInterest()->first());
      if (!var) return;
      double oldval = var->getVal();
      if (poiValue > 0)  var->setVal(poiValue);
      //sbModel->SetSnapshot( RooArgSet(*var, *w->var("lumi") ) );
      sbModel->SetSnapshot( RooArgSet(*var) );
      if (poiValue > 0) var->setVal(oldval);
      //sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );
   }

   



   // part 1, hypothesis testing 
   SimpleLikelihoodRatioTestStat * slrts = new SimpleLikelihoodRatioTestStat(*bModel->GetPdf(), *sbModel->GetPdf());
   // null parameters must includes snapshot of poi plus the nuisance values 
   RooArgSet nullParams(*bModel->GetSnapshot());
   if (bModel->GetNuisanceParameters()) nullParams.add(*bModel->GetNuisanceParameters());
   
   slrts->SetNullParameters(nullParams);
   RooArgSet altParams(*sbModel->GetSnapshot());
   if (sbModel->GetNuisanceParameters()) altParams.add(*sbModel->GetNuisanceParameters());
   slrts->SetAltParameters(altParams);


   ProfileLikelihoodTestStat * profll = new ProfileLikelihoodTestStat(*bModel->GetPdf());

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

      pdflist.add( *allNuisancePdfs ) ;

      pdflist.Print() ;
      printf("\n") ;

      RooProdPdf* likelihood = new RooProdPdf( "likelihood", "hbb likelihood", pdflist ) ;
      likelihood->Print() ;


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


  //  printf("\n\n Running a test fit.\n\n") ;


  //  dsObserved -> Print() ;
  //  dsObserved -> printMultiline(cout, 1, kTRUE, "") ;


  //  printf("\n\n =============================================\n\n") ;
  //  likelihood -> fitTo( *dsObserved, PrintLevel(3), Hesse(0), Minos(0) ) ;
  //  printf("\n\n =============================================\n\n") ;







     //-- Set up RooStats models.

      printf("\n\n Setting up S+B model.\n\n") ;

      RooArgSet poi( *rv_sig_strength, "poi" ) ;
      RooUniform signal_prior( "signal_prior", "signal_prior", *rv_sig_strength ) ;

      ModelConfig sbModel ("SbModel");
      sbModel.SetWorkspace( workspace ) ;
      sbModel.SetPdf( *likelihood ) ;
      sbModel.SetParametersOfInterest( poi );
      sbModel.SetPriorPdf(signal_prior);
      sbModel.SetObservables( *observedParametersList );
      sbModel.SetNuisanceParameters( *allNuisances );
      sbModel.SetGlobalObservables( *globalObservables );

      workspace.Print() ;

      printf("\n\n Doing fit for S+B model.\n" ) ; fflush(stdout) ;

      RooAbsReal* pNll = sbModel.GetPdf()->createNLL(*dsObserved);
      RooAbsReal* pProfile = pNll->createProfile(RooArgSet());
      pProfile->getVal();
      RooArgSet* pPoiAndNuisance = new RooArgSet();
      pPoiAndNuisance->add(*sbModel.GetParametersOfInterest());
      if(sbModel.GetNuisanceParameters()) pPoiAndNuisance->add(*sbModel.GetNuisanceParameters());
      printf("\n\n Will save these parameter points that correspond to the fit to data.\n\n") ; fflush(stdout) ;
      pPoiAndNuisance->Print("v");
      sbModel.SetSnapshot(*pPoiAndNuisance);
      workspace.import (sbModel);

      delete pProfile ;
      delete pNll ;
      delete pPoiAndNuisance ;

      printf("\n\n Setting up BG-only model.\n\n") ;

      ModelConfig bModel (*(RooStats::ModelConfig *)workspace.obj("SbModel"));
      bModel.SetName("BModel");
      bModel.SetWorkspace(workspace);

      printf("\n\n Doing fit for BG-only model.\n" ) ; fflush(stdout) ;
      pNll = bModel.GetPdf()->createNLL(*dsObserved);
      pProfile = pNll->createProfile(*bModel.GetParametersOfInterest());
      ((RooRealVar *)(bModel.GetParametersOfInterest()->first()))->setVal(0.);
      pProfile->getVal();
      pPoiAndNuisance = new RooArgSet();
      pPoiAndNuisance->add(*bModel.GetParametersOfInterest());
      if(bModel.GetNuisanceParameters()) pPoiAndNuisance->add(*bModel.GetNuisanceParameters());
      printf("\n\n Should use these parameter points to generate pseudo data for bkg only.\n\n") ; fflush(stdout) ;
      pPoiAndNuisance->Print("v");
      bModel.SetSnapshot(*pPoiAndNuisance);
      workspace.import (bModel);

      delete pProfile ;
      delete pNll ;
      delete pPoiAndNuisance ;

      workspace.Print() ;

      printf("\n\n Saving workspace in : %s\n\n", outfile ) ;

      gSystem->Exec(" mkdir -p outputfiles " ) ;

      workspace.writeToFile( outfile ) ;




   } // build_hbb_workspace1.

void test_counting_experiment() {

//////////////////////   MODEL BUILDING    /////////////////////////////////
///////////////////////////////////////////////////////////////////////////
/*
N_s = N_tot_theory(Mass,Xsec) * Acceptance_SR * Eff_AmBe_bin_i * mu
N_b = N_Co_SR_bin_i * Norm_factor

Xesec: considered 10^-40 cm^2
Norm_factor = N_Data_CR / N_Co_CR   --> assuming no difference between Co and Data in CR and SR.
N_tot_theory(Mass,Xsec): for 225 livedays, 45kg and considering Xsec. It is a constant, no uncertainty at the moment.

---Costraint Signal
   nuissance parameter = Acceptance_SR, Eff_AmBe_bin_i
   Gauss(Acceptance_SR_obs | Acceptance_SR, err.)
   Poisson(S0_i | S_tot_SR  * Eff_AmBe_bin_i)
   
---Costraint Bkg
   nuissance parameter = N_Co_SR_bin_i, Norm_factor
   Gauss(Norm_factor_obs |  Norm_factor, err)
   Poisson(B0_i | N_Co_SR_bin_i)

---- WARNING:: convergence problems: mu_hat should always be >> 1, too small values have problem in finding minimum
	because mu is set >0. ---> Try to fix Xsec in order to have mu_hat ~ 10

*/

  RooWorkspace w("w");

  //gROOT->ProcessLine(".L retrieve_input_from_histo_NoSys.C+");
  gROOT->ProcessLine(".L retrieve_input_from_histo.C+");

  retrieve_input_from_histo(w);


// Building the model
   ModelConfig mc("ModelConfig",&w);
   mc.SetPdf(*w.pdf("model"));
   mc.SetParametersOfInterest(*w.var("mu"));

// Setting nuissance parameter
   mc.SetNuisanceParameters(*w.set("nuissance_parameter"));

// need now to set the global observable
   mc.SetGlobalObservables(*w.set("g_observables"));

   mc.SetObservables(*w.set("observables"));

// this is needed for the hypothesis tests
   mc.SetSnapshot(*w.var("mu"));


// make data set with the number of observed events
   RooDataSet data("data","", *w.set("observables"));
   data.add(*w.set("observables"));

// import data set in workspace and save it in a file
   w.import(data);

// import model in the workspace 
   w.import(mc);

   w.writeToFile("CountingModel.root", true);




w.Print();

data.Print();

/* 
cout << w.var("S_i")->getValV() << endl;//<< "   "   <<  w.var("S_i_exp")->getValV() << endl;
///////////////////////////////////////////////////////////////////////
ProfileLikelihoodCalculator pl(data,mc);
  pl.SetConfidenceLevel(0.95);
  LikelihoodInterval* interval = pl.GetInterval();

   // find the iterval on the first Parameter of Interest
  RooRealVar* firstPOI = (RooRealVar*) mc.GetParametersOfInterest()->first();

  double lowerLimit = interval->LowerLimit(*firstPOI);
  double upperLimit = interval->UpperLimit(*firstPOI);


  cout << "\n95% interval on " <<firstPOI->GetName()<<" is : ["<<
    lowerLimit << ", "<<
    upperLimit <<"] "<<endl;


  LikelihoodIntervalPlot * plot = new LikelihoodIntervalPlot(interval);
//  plot->SetRange(0,50);  // possible eventually to change ranges
  //plot->SetNPoints(50);  // do not use too many points, it could become very slow for some models
  plot->Draw("");  // use option TF1 if too slow (plot.Draw("tf1")

*/




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

void significance(RooWorkspace& w ) {

  ModelConfig* mc = (ModelConfig*)w.obj("mc");
  RooDataSet* data = (RooDataSet*)w.data("data");
  //data->Print();

  // define the S+B snapshot (this is used for computing the expected significance)
  ModelConfig* sbModel = mc->Clone();
  sbModel->SetName("S+B Model");
  RooRealVar* poi = (RooRealVar*) sbModel->GetParametersOfInterest()->first();
  poi->setVal(50);
  sbModel->SetSnapshot(*poi);

  ModelConfig * bModel = (ModelConfig*) sbModel->Clone();
  bModel->SetName("B model");
  poi->setVal(0);
  bModel->SetSnapshot(*poi);

  vector<double> masses;
  vector<double> p0values;
  vector<double> p0valuesExpected;
  vector<double> sigvalues;

  double massMin = 200;
  double massMax = 2500;
  int nbins = 100;

  // loop on the mass values 
  for ( double mass=massMin; mass<=massMax; mass += (massMax-massMin)/nbins ) {

    w.var("mass")->setVal( mass );

    // create the AsymptoticCalculator from data,alt model, null model
    AsymptoticCalculator * ac = new AsymptoticCalculator(*data, *sbModel, *bModel);
    ac->SetOneSidedDiscovery(true);  // for one-side discovery test                                      
    AsymptoticCalculator::SetPrintLevel(-1);

    // run the calculator
    HypoTestResult* asymCalcResult = ac->GetHypoTest();
    asymCalcResult->Print();
     
    double pvalue = asymCalcResult->NullPValue();
    double sigvalue = asymCalcResult->Significance();
    double expectedP0 = AsymptoticCalculator::GetExpectedPValues(asymCalcResult->NullPValue(),asymCalcResult->AlternatePValue(), 0, false);

    masses.push_back(mass);
    p0values.push_back(pvalue);
    p0valuesExpected.push_back(expectedP0);
    sigvalues.push_back(sigvalue);
    std::cout << "** Mass = " << mass << " p0-value = " << expectedP0 << " p-value = " << pvalue << " significance = " << sigvalue << std::endl;

  }

  TGraph* graph1  = new TGraph(masses.size(),&masses[0],&p0values[0]);
  TGraph* graph2  = new TGraph(masses.size(),&masses[0],&p0valuesExpected[0]);
  TGraph* graph3  = new TGraph(masses.size(),&masses[0],&sigvalues[0]);

  TCanvas* c2 = new TCanvas("c2","Significance", 900, 700);
  c2->Divide(1,2);
  c2->cd(1);
  graph1->SetMarkerStyle(10);
  //graph1->Draw("APC");
  graph1->Draw("AC");
  graph2->SetLineStyle(2);
  graph2->Draw("C");
  graph1->GetXaxis()->SetTitle("Mass [GeV]");
  graph1->GetYaxis()->SetTitle("p0 value");
  graph1->SetTitle("P-value vs Mass");
  graph1->SetMinimum(graph2->GetMinimum());
  graph1->SetLineColor(kBlue);
  graph2->SetLineColor(kRed);
  gPad->SetLogy(true);

  c2->cd(2);
  graph3->SetMarkerStyle(10);
  graph3->Draw("AC");
  graph3->SetLineStyle(1);
  graph3->SetLineColor(kRed);
  graph3->GetXaxis()->SetTitle("Mass [GeV]");
  graph3->GetYaxis()->SetTitle("Significance");
  graph3->SetTitle("Significance vs Mass");
  gPad->SetLogy(false);

  c2->SaveAs("significance.pdf");
  c2->SaveAs("significance.png");
}