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

// 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;
         }         
      }
   }

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

//.........这里部分代码省略.........
    channelNames.insert("SF_CZpeak_1j"+string(!do2011?"_2012":""));
    channelNames.insert("OF_CfrecZpeak_1j"+string(!do2011?"_2012":""));
    channelNames.insert("OF_CZpeak_1j"+string(!do2011?"_2012":""));
    channelNames.insert("OF_mainControl_1j"+string(!do2011?"_2012":""));
    channelNames.insert("OF_topbox_1j"+string(!do2011?"_2012":""));

    channelNames.insert("ee_signalLike1_2j"+string(!do2011?"_2012":""));
    channelNames.insert("SF_topbox_2j"+string(!do2011?"_2012":""));
  }
  else {
    cout << "Channel " << channel << " not defined. Please check!" << endl;
    exit(1);
  }

  // bool fix = 1;
  stringstream inFileName;

  inFileName << "workspaces/" << inFolderName << "/" << mass << ".root";
  TFile f(inFileName.str().c_str());
  
  RooWorkspace* w = (RooWorkspace*)f.Get("combWS");
  if (!w) w = (RooWorkspace*)f.Get("combined");
  
  RooDataSet* data = (RooDataSet*)w->data("combData");
  if (!data) data = (RooDataSet*)w->data("obsData");
  
  ModelConfig* mc = (ModelConfig*)w->obj("ModelConfig");
  
  RooRealVar* weightVar = w->var("weightVar");
  
  RooRealVar* mu = (RooRealVar*)mc->GetParametersOfInterest()->first();
  if (!mu) mu = w->var("SigXsecOverSM");

  const RooArgSet* mc_obs = mc->GetObservables();
  const RooArgSet* mc_nuis = mc->GetNuisanceParameters();
  const RooArgSet* mc_globs = mc->GetGlobalObservables();
  const RooArgSet* mc_poi = mc->GetParametersOfInterest();

  RooArgSet nuis = *mc_nuis;
  RooArgSet antiNuis = *mc_nuis;

  RooArgSet globs = *mc_globs;
  RooArgSet antiGlobs = *mc_globs;

  RooArgSet allParams;

  RooSimultaneous* simPdf = (RooSimultaneous*)mc->GetPdf();
  RooCategory* cat = (RooCategory*)&simPdf->indexCat();

  RooArgSet nuis_tmp = nuis;
  RooArgSet fullConstraints = *simPdf->getAllConstraints(*mc_obs,nuis_tmp,false);

  vector<string> foundChannels;
  vector<string> skippedChannels;  

  cout << "Getting constraints" << endl;
  map<string, RooDataSet*> data_map;
  map<string, RooAbsPdf*> pdf_map;
  RooCategory* decCat = new RooCategory("dec_channel","dec_channel");
  // int i = 0;
  TIterator* catItr = cat->typeIterator();
  RooCatType* type;
  RooArgSet allConstraints;
  while ((type = (RooCatType*)catItr->Next())) {
    RooAbsPdf* pdf =  simPdf->getPdf(type->GetName());

void runQ(const char* inFileName,
	    const char* wsName = "combined",
	    const char* modelConfigName = "ModelConfig",
	    const char* dataName = "obsData",
	    const char* asimov0DataName = "asimovData_0",
	    const char* conditional0Snapshot = "conditionalGlobs_0",
	    const char* asimov1DataName = "asimovData_1",
	    const char* conditional1Snapshot = "conditionalGlobs_1",
	    const char* nominalSnapshot = "nominalGlobs",
	    string smass = "130",
	    string folder = "test")
{
  double mass;
  stringstream massStr;
  massStr << smass;
  massStr >> mass;

  bool errFast = 0;
  bool goFast = 1;
  bool remakeData = 1;
  bool doRightSided = 1;
  bool doInj = 0;
  bool doObs = 1;
  bool doMedian = 1;

  TStopwatch timer;
  timer.Start();

  TFile f(inFileName);
  RooWorkspace* ws = (RooWorkspace*)f.Get(wsName);
  if (!ws)
  {
    cout << "ERROR::Workspace: " << wsName << " doesn't exist!" << endl;
    return;
  }
  ModelConfig* mc = (ModelConfig*)ws->obj(modelConfigName);
  if (!mc)
  {
    cout << "ERROR::ModelConfig: " << modelConfigName << " doesn't exist!" << endl;
    return;
  }
  RooDataSet* data = (RooDataSet*)ws->data(dataName);
  if (!data)
  {
    cout << "ERROR::Dataset: " << dataName << " doesn't exist!" << endl;
    return;
  }





  mc->GetNuisanceParameters()->Print("v");

  RooNLLVar::SetIgnoreZeroEntries(1);
  ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
  ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
  ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(1);
  cout << "Setting max function calls" << endl;
  //ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(20000);
  RooMinimizer::SetMaxFunctionCalls(10000);

  ws->loadSnapshot("conditionalNuis_0");
  RooArgSet nuis(*mc->GetNuisanceParameters());

  RooRealVar* mu = (RooRealVar*)mc->GetParametersOfInterest()->first();



  if (string(mc->GetPdf()->ClassName()) == "RooSimultaneous" && remakeData)
  {
    RooSimultaneous* simPdf = (RooSimultaneous*)mc->GetPdf();
    double min_mu;
    data = makeData(data, simPdf, mc->GetObservables(), mu, mass, min_mu);
  }







  RooDataSet* asimovData0 = (RooDataSet*)ws->data(asimov0DataName);
  if (!asimovData0)
  {
    cout << "Asimov data doesn't exist! Please, allow me to build one for you..." << endl;
    makeAsimovData(mc, true, ws, mc->GetPdf(), data, 1);
    ws->Print();
    asimovData0 = (RooDataSet*)ws->data("asimovData_0");
  }

  RooDataSet* asimovData1 = (RooDataSet*)ws->data(asimov1DataName);
  if (!asimovData1)
  {
    cout << "Asimov data doesn't exist! Please, allow me to build one for you..." << endl;
    makeAsimovData(mc, true, ws, mc->GetPdf(), data, 0);
    ws->Print();
    asimovData1 = (RooDataSet*)ws->data("asimovData_1");
  }
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   for(Int_t i=0; i<parameterScan->numEntries(); ++i){
      tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
      //cout <<"get threshold"<<endl;
      double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
      double poiVal = tmpPoint->getRealValue(firstPOI->GetName()) ;
      histOfThresholds->Fill(poiVal,arMax);
   }
   TCanvas* c1 = new TCanvas();
   c1->Divide(2);
   c1->cd(1);
   histOfThresholds->SetMinimum(0);
   histOfThresholds->Draw();
   c1->cd(2);

   // -------------------------------------------------------
   // Now we generate the expected bands and power-constraint

   // First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
   RooAbsReal* nll = mc->GetPdf()->createNLL(*data);
   RooAbsReal* profile = nll->createProfile(*mc->GetParametersOfInterest());
   firstPOI->setVal(0.);
   profile->getVal(); // this will do fit and set nuisance parameters to profiled values
   RooArgSet* poiAndNuisance = new RooArgSet();
   if(mc->GetNuisanceParameters())
      poiAndNuisance->add(*mc->GetNuisanceParameters());
   poiAndNuisance->add(*mc->GetParametersOfInterest());
   w->saveSnapshot("paramsToGenerateData",*poiAndNuisance);
   RooArgSet* paramsToGenerateData = (RooArgSet*) poiAndNuisance->snapshot();
   cout << "\nWill use these parameter points to generate pseudo data for bkg only" << endl;
   paramsToGenerateData->Print("v");


   RooArgSet unconditionalObs;
   unconditionalObs.add(*mc->GetObservables());
   unconditionalObs.add(*mc->GetGlobalObservables()); // comment this out for the original conditional ensemble

   double CLb=0;
   double CLbinclusive=0;

   // Now we generate background only and find distribution of upper limits
   TH1F* histOfUL = new TH1F("histOfUL","",100,0,firstPOI->getMax());
   histOfUL->GetXaxis()->SetTitle("Upper Limit (background only)");
   histOfUL->GetYaxis()->SetTitle("Entries");
   for(int imc=0; imc<nToyMC; ++imc){

      // set parameters back to values for generating pseudo data
      //    cout << "\n get current nuis, set vals, print again" << endl;
      w->loadSnapshot("paramsToGenerateData");
      //    poiAndNuisance->Print("v");

      RooDataSet* toyData = 0;
      // now generate a toy dataset
      if(!mc->GetPdf()->canBeExtended()){
         if(data->numEntries()==1)
            toyData = mc->GetPdf()->generate(*mc->GetObservables(),1);
         else
            cout <<"Not sure what to do about this model" <<endl;
      } else{
         //      cout << "generating extended dataset"<<endl;
         toyData = mc->GetPdf()->generate(*mc->GetObservables(),Extended());
      }

      // generate global observables
      // need to be careful for simpdf
      //    RooDataSet* globalData = mc->GetPdf()->generate(*mc->GetGlobalObservables(),1);

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

   // profll.SetReuseNLL(mOptimize);
   // slrts.SetReuseNLL(mOptimize);
   // ropl.SetReuseNLL(mOptimize);

   AsymptoticCalculator::SetPrintLevel(printLevel);

   HypoTestCalculatorGeneric *  hypoCalc = 0;
   // note here Null is B and Alt is S+B
   if (calcType == 0) hypoCalc = new  FrequentistCalculator(*data, *sbModel, *bModel);
   else if (calcType == 1) hypoCalc= new  HybridCalculator(*data, *sbModel, *bModel);
   else if (calcType == 2) hypoCalc= new  AsymptoticCalculator(*data, *sbModel, *bModel);

   if (calcType == 0) 
       ((FrequentistCalculator*)hypoCalc)->SetToys(ntoys, ntoys/nToysRatio);
   if (calcType == 1) 
       ((HybridCalculator*)hypoCalc)->SetToys(ntoys, ntoys/nToysRatio);
   if (calcType == 2 ) { 
      if (testStatType == 3) ((AsymptoticCalculator*) hypoCalc)->SetOneSidedDiscovery(true);  
      if (testStatType != 2 && testStatType != 3)  
         Warning("StandardHypoTestDemo","Only the PL test statistic can be used with AsymptoticCalculator - use by default a two-sided PL");
      

   }


   // check for nuisance prior pdf in case of nuisance parameters 
   if (calcType == 1 && (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters() )) {
         RooAbsPdf * nuisPdf = 0; 
         if (nuisPriorName) nuisPdf = w->pdf(nuisPriorName);
         // use prior defined first in bModel (then in SbModel)
         if (!nuisPdf)  { 
            Info("StandardHypoTestDemo","No nuisance pdf given for the HybridCalculator - try to deduce  pdf from the   model");
            if (bModel->GetPdf() && bModel->GetObservables() ) 
               nuisPdf = RooStats::MakeNuisancePdf(*bModel,"nuisancePdf_bmodel");
            else 
               nuisPdf = RooStats::MakeNuisancePdf(*sbModel,"nuisancePdf_sbmodel");
         }   
         if (!nuisPdf ) {
            if (bModel->GetPriorPdf())  { 
               nuisPdf = bModel->GetPriorPdf();
               Info("StandardHypoTestDemo","No nuisance pdf given - try to use %s that is defined as a prior pdf in the B model",nuisPdf->GetName());            
            }
            else { 
               Error("StandardHypoTestDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified or can be derived");
               return;
            }
         }
         assert(nuisPdf);
         Info("StandardHypoTestDemo","Using as nuisance Pdf ... " );
         nuisPdf->Print();
      
         const RooArgSet * nuisParams = (bModel->GetNuisanceParameters() ) ? bModel->GetNuisanceParameters() : sbModel->GetNuisanceParameters();
         RooArgSet * np = nuisPdf->getObservables(*nuisParams);
         if (np->getSize() == 0) { 
            Warning("StandardHypoTestDemo","Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range");
         }
         delete np;
      
         ((HybridCalculator*)hypoCalc)->ForcePriorNuisanceAlt(*nuisPdf);
         ((HybridCalculator*)hypoCalc)->ForcePriorNuisanceNull(*nuisPdf);
   }

   // hypoCalc->ForcePriorNuisanceAlt(*sbModel->GetPriorPdf());
   // hypoCalc->ForcePriorNuisanceNull(*bModel->GetPriorPdf());

result fit_toy(RooWorkspace* wspace, int n, const RooArgSet* globals) {
    RooRandom::randomGenerator()->SetSeed(0);
    // TFile f(filename);
    // RooWorkspace *wspace = (RooWorkspace*)f.Get("combined");
    ModelConfig* model = (ModelConfig*)wspace->obj("ModelConfig");

    RooAbsPdf* pdf;
    pdf = model->GetPdf();

    RooAbsPdf* top_constraint = (RooAbsPdf*)wspace->obj("top_ratio_constraint");
    RooAbsPdf* vv_constraint = (RooAbsPdf*)wspace->obj("vv_ratio_constraint");
    RooAbsPdf* top_vv_constraint_sf = (RooAbsPdf*)wspace->obj("top_vv_ratio_sf_constraint");
    RooAbsPdf* top_vv_constraint_of = (RooAbsPdf*)wspace->obj("top_vv_ratio_of_constraint");


    // generate constraint global observables
    RooRealVar *nom_top_ratio = (RooRealVar*)wspace->obj("nom_top_ratio");
    nom_top_ratio->setRange(0, 100);
    RooRealVar *nom_vv_ratio = (RooRealVar*)wspace->obj("nom_vv_ratio");
    nom_vv_ratio->setRange(0,100);
    RooRealVar *nom_top_vv_ratio_sf = (RooRealVar*)wspace->obj("nom_top_vv_ratio_sf");
    nom_top_vv_ratio_sf->setRange(0,100);
    RooRealVar *nom_top_vv_ratio_of = (RooRealVar*)wspace->obj("nom_top_vv_ratio_of");
    nom_top_vv_ratio_of->setRange(0,100);

    RooDataSet *nom_top_generated = top_constraint->generateSimGlobal(RooArgSet(*nom_top_ratio), 1);
    nom_top_ratio->setVal(((RooRealVar*)nom_top_generated->get(0)->find("nom_top_ratio"))->getVal());

    RooDataSet *nom_vv_generated = vv_constraint->generateSimGlobal(RooArgSet(*nom_vv_ratio), 1);
    nom_vv_ratio->setVal(((RooRealVar*)nom_vv_generated->get(0)->find("nom_vv_ratio"))->getVal());

    RooDataSet *nom_top_vv_sf_generated = top_vv_constraint_sf->generateSimGlobal(RooArgSet(*nom_top_vv_ratio_sf), 1);
    nom_top_vv_ratio_sf->setVal(((RooRealVar*)nom_top_vv_sf_generated->get(0)->find("nom_top_vv_ratio_sf"))->getVal());

    RooDataSet *nom_top_vv_of_generated = top_vv_constraint_of->generateSimGlobal(RooArgSet(*nom_top_vv_ratio_of), 1);
    nom_top_vv_ratio_of->setVal(((RooRealVar*)nom_top_vv_of_generated->get(0)->find("nom_top_vv_ratio_of"))->getVal());

    NumEventsTestStat* dummy = new NumEventsTestStat(*pdf);

    ToyMCSampler* mc = new ToyMCSampler(*dummy, 1);
    mc->SetPdf(*pdf);
    mc->SetObservables(*model->GetObservables());
    mc->SetGlobalObservables(*globals);
    mc->SetNuisanceParameters(*model->GetNuisanceParameters());
    mc->SetParametersForTestStat(*model->GetParametersOfInterest());
    mc->SetNEventsPerToy(n);

    RooArgSet constr;
    constr.add(*(model->GetNuisanceParameters()));
    RemoveConstantParameters(&constr);

    RooDataSet* toy_data = (RooDataSet*)mc->GenerateToyData(*const_cast<RooArgSet*>(model->GetSnapshot()));

    RooFitResult *res = pdf->fitTo(*toy_data, Constrain(constr), PrintLevel(0), Save(),
                                               Range("fitRange"), InitialHesse(),
                                               ExternalConstraints(RooArgSet(*top_constraint, *vv_constraint, *top_vv_constraint_sf, *top_vv_constraint_of)));
    result yield = get_results(wspace, res);
    yield.of.generated_sum.val = toy_data->sumEntries("(channelCat==channelCat::of) & (obs_x_of>120)");
    yield.sf.generated_sum.val = toy_data->sumEntries("(channelCat==channelCat::sf) & (obs_x_sf>120)");

    delete mc;
    delete dummy;
    // f.Close();

    return yield;
}

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); 
//.........这里部分代码省略.........

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


   double nSigmaToVary=5.;
   double muVal=0;
   bool doFit=false;

   // -------------------------------------------------------
   // 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;
   }

   // -------------------------------------------------------
   // now use the profile inspector

   RooRealVar* obs = (RooRealVar*)mc->GetObservables()->first();
   TList* list = new TList();


   RooRealVar * firstPOI = dynamic_cast<RooRealVar*>(mc->GetParametersOfInterest()->first());

   firstPOI->setVal(muVal);
   //  firstPOI->setConstant();
   if(doFit){
      mc->GetPdf()->fitTo(*data);
   }

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


   mc->GetNuisanceParameters()->Print("v");
   int  nPlotsMax = 1000;
   cout <<" check expectedData by category"<<endl;
   RooDataSet* simData=NULL;
   RooSimultaneous* simPdf = NULL;
   if(strcmp(mc->GetPdf()->ClassName(),"RooSimultaneous")==0){
      cout <<"Is a simultaneous PDF"<<endl;
      simPdf = (RooSimultaneous *)(mc->GetPdf());
   } else {
      cout <<"Is not a simultaneous PDF"<<endl;
   }


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

//.........这里部分代码省略.........
   // 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;
  }

  mc->Print();
  /////////////////////////////////////////////////////////////
  // Now find the upper limit based on the asymptotic results
  ////////////////////////////////////////////////////////////
  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  ProfileLikelihoodCalculator plc(*data,*mc);
  LikelihoodInterval* interval = plc.GetInterval();
  double plcUpperLimit = interval->UpperLimit(*firstPOI);
  delete interval;
  cout << "\n\n--------------------------------------"<<endl;
  cout <<"Will generate sampling distribution at " << firstPOI->GetName() << " = " << plcUpperLimit <<endl;
  int nPOI = mc->GetParametersOfInterest()->getSize();
  if(nPOI>1){
    cout <<"not sure what to do with other parameters of interest, but here are their values"<<endl;
    mc->GetParametersOfInterest()->Print("v");
  }

  /////////////////////////////////////////////
  // create thte test stat sampler
  ProfileLikelihoodTestStat ts(*mc->GetPdf());

  // to avoid effects from boundary and simplify asymptotic comparison, set min=-max
  if(allowNegativeMu)
    firstPOI->setMin(-1*firstPOI->getMax());

  // temporary RooArgSet
  RooArgSet poi;
  poi.add(*mc->GetParametersOfInterest());

  // create and configure the ToyMCSampler
  ToyMCSampler sampler(ts,nToyMC);
  sampler.SetPdf(*mc->GetPdf());
  sampler.SetObservables(*mc->GetObservables());
  sampler.SetGlobalObservables(*mc->GetGlobalObservables());
  if(!mc->GetPdf()->canBeExtended() && (data->numEntries()==1)){
    cout << "tell it to use 1 event"<<endl;
    sampler.SetNEventsPerToy(1);
  }
  firstPOI->setVal(plcUpperLimit); // set POI value for generation
  sampler.SetParametersForTestStat(*mc->GetParametersOfInterest()); // set POI value for evaluation

  if (useProof) {
     ProofConfig pc(*w, nworkers, "",false);
     sampler.SetProofConfig(&pc); // enable proof
  }

  firstPOI->setVal(plcUpperLimit);
  RooArgSet allParameters;
  allParameters.add(*mc->GetParametersOfInterest());
  allParameters.add(*mc->GetNuisanceParameters());
  allParameters.Print("v");

  SamplingDistribution* sampDist = sampler.GetSamplingDistribution(allParameters);
  SamplingDistPlot plot;
  plot.AddSamplingDistribution(sampDist);
  plot.GetTH1F(sampDist)->GetYaxis()->SetTitle(Form("f(-log #lambda(#mu=%.2f) | #mu=%.2f)",plcUpperLimit,plcUpperLimit));
  plot.SetAxisTitle(Form("-log #lambda(#mu=%.2f)",plcUpperLimit));

  TCanvas* c1 = new TCanvas("c1");
  c1->SetLogy();
  plot.Draw();
  double min = plot.GetTH1F(sampDist)->GetXaxis()->GetXmin();
  double max = plot.GetTH1F(sampDist)->GetXaxis()->GetXmax();

  TF1* f = new TF1("f",Form("2*ROOT::Math::chisquared_pdf(2*x,%d,0)",nPOI),min,max);
  f->Draw("same");
  c1->SaveAs("standard_test_stat_distribution.pdf");

}