//void RunToyScan5(TString fileName, double startVal, double stopVal, TString outFile) {
void frequentist(TString fileName) {
  cout << "Starting frequentist " << time(NULL) << endl;
  double startVal = 0;
  double stopVal = 200;
  TString outFile = "";

  int nToys = 1 ;
  int nscanpoints = 2 ;

  /*
  gROOT->LoadMacro("RooBetaPdf.cxx+") ;
  gROOT->LoadMacro("RooRatio.cxx+") ;
  gROOT->LoadMacro("RooPosDefCorrGauss.cxx+") ;
  */

  // get relevant objects out of the "ws" file

  TFile *file = TFile::Open(fileName);
  if(!file){
    cout <<"file not found" << endl;
    return;
  } 

  RooWorkspace* w = (RooWorkspace*) file->Get("workspace");
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  ModelConfig* mc = (ModelConfig*) w->obj("S+B_model");
  RooAbsData* data = w->data("data");

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

  RooRealVar* myPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  myPOI->setRange(0, 1000.);

  ModelConfig* bModel = (ModelConfig*) w->obj("B_model");
  ModelConfig* sbModel = (ModelConfig*) w->obj("S+B_model");

  ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
  profll.SetPrintLevel(2);
  profll.SetOneSided(1);
  TestStatistic * testStat = &profll;

  HypoTestCalculatorGeneric *  hc = 0;
  hc = new FrequentistCalculator(*data, *bModel, *sbModel);
  
  ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
  toymcs->SetMaxToys(10000);
  toymcs->SetNEventsPerToy(1);
  toymcs->SetTestStatistic(testStat);


  ((FrequentistCalculator *)hc)->SetToys(nToys,nToys);
  
  HypoTestInverter calc(*hc);
  calc.SetConfidenceLevel(0.95);
  calc.UseCLs(true);
  //calc.SetVerbose(true);
  calc.SetVerbose(2);

  cout << "About to set fixed scan " << time(NULL) << endl;
  calc.SetFixedScan(nscanpoints,startVal,stopVal);
  cout << "About to do inverter " << time(NULL) << endl;
  HypoTestInverterResult * res_toysCLs_calculator = calc.GetInterval();

  cout << "CLs = " << res_toysCLs_calculator->UpperLimit() 
	    << "   CLs_exp = " << res_toysCLs_calculator->GetExpectedUpperLimit(0) 
	    << "   CLs_exp(-1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(-1) 
	    << "   CLs_exp(+1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(1) << endl ;

  /*
  // dump results string to output file
  ofstream outStream ;
  outStream.open(outFile,ios::app) ;
  
  outStream << "CLs = " << res_toysCLs_calculator->UpperLimit() 
	    << "   CLs_exp = " << res_toysCLs_calculator->GetExpectedUpperLimit(0) 
	    << "   CLs_exp(-1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(-1) 
	    << "   CLs_exp(+1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(1) << endl ;
  
  outStream.close() ;
  */


  cout << "End of frequentist " << time(NULL) << endl;
  return ;

}

//.........这里部分代码省略.........
  
   MaxLikelihoodEstimateTestStat maxll(*sbModel->GetPdf(),*poi); 
   NumEventsTestStat nevtts;

   AsymptoticCalculator::SetPrintLevel(mPrintLevel);
  
   // create the HypoTest calculator class 
   HypoTestCalculatorGeneric *  hc = 0;
   if (type == 0) hc = new FrequentistCalculator(*data, *bModel, *sbModel);
   else if (type == 1) hc = new HybridCalculator(*data, *bModel, *sbModel);
   // else if (type == 2 ) hc = new AsymptoticCalculator(*data, *bModel, *sbModel, false, mAsimovBins);
   // else if (type == 3 ) hc = new AsymptoticCalculator(*data, *bModel, *sbModel, true, mAsimovBins);  // for using Asimov data generated with nominal values 
   else if (type == 2 ) hc = new AsymptoticCalculator(*data, *bModel, *sbModel, false );
   else if (type == 3 ) hc = new AsymptoticCalculator(*data, *bModel, *sbModel, true );  // for using Asimov data generated with nominal values 
   else {
      Error("StandardHypoTestInvDemo","Invalid - calculator type = %d supported values are only :\n\t\t\t 0 (Frequentist) , 1 (Hybrid) , 2 (Asymptotic) ",type);
      return 0;
   }
  
   // set the test statistic 
   TestStatistic * testStat = 0;
   if (testStatType == 0) testStat = &slrts;
   if (testStatType == 1 || testStatType == 11) testStat = &ropl;
   if (testStatType == 2 || testStatType == 3 || testStatType == 4) testStat = &profll;
   if (testStatType == 5) testStat = &maxll;
   if (testStatType == 6) testStat = &nevtts;

   if (testStat == 0) { 
      Error("StandardHypoTestInvDemo","Invalid - test statistic type = %d supported values are only :\n\t\t\t 0 (SLR) , 1 (Tevatron) , 2 (PLR), 3 (PLR1), 4(MLE)",testStatType);
      return 0;
   }
  
  
   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   if (toymcs && (type == 0 || type == 1) ) { 
      // look if pdf is number counting or extended
      if (sbModel->GetPdf()->canBeExtended() ) { 
         if (useNumberCounting)   Warning("StandardHypoTestInvDemo","Pdf is extended: but number counting flag is set: ignore it ");
      }
      else { 
         // for not extended pdf
         if (!useNumberCounting  )  { 
            int nEvents = data->numEntries();
            Info("StandardHypoTestInvDemo","Pdf is not extended: number of events to generate taken  from observed data set is %d",nEvents);
            toymcs->SetNEventsPerToy(nEvents);
         }
         else {
            Info("StandardHypoTestInvDemo","using a number counting pdf");
            toymcs->SetNEventsPerToy(1);
         }
      }

      toymcs->SetTestStatistic(testStat);
    
      if (data->isWeighted() && !mGenerateBinned) { 
         Info("StandardHypoTestInvDemo","Data set is weighted, nentries = %d and sum of weights = %8.1f but toy generation is unbinned - it would be faster to set mGenerateBinned to true\n",data->numEntries(), data->sumEntries());
      }
      toymcs->SetGenerateBinned(mGenerateBinned);
  
      toymcs->SetUseMultiGen(mOptimize);
    
      if (mGenerateBinned &&  sbModel->GetObservables()->getSize() > 2) { 
         Warning("StandardHypoTestInvDemo","generate binned is activated but the number of ovservable is %d. Too much memory could be needed for allocating all the bins",sbModel->GetObservables()->getSize() );
      }

      // set the random seed if needed

//.........这里部分代码省略.........
            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);
   else hc = new HybridCalculator(*data, *bModel, *sbModel);

   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   //=== DEBUG
   ///// toymcs->SetWS( w ) ;
   //=== DEBUG
   toymcs->SetNEventsPerToy(1);
   toymcs->SetTestStatistic(testStat);
   if (optimize) toymcs->SetUseMultiGen(true);


   if (type == 1) { 
      HybridCalculator *hhc = (HybridCalculator*) hc;
      hhc->SetToys(ntoys,ntoys); 

      // check for nuisance prior pdf 
      if (bModel->GetPriorPdf() && sbModel->GetPriorPdf() ) {
         hhc->ForcePriorNuisanceAlt(*bModel->GetPriorPdf());
         hhc->ForcePriorNuisanceNull(*sbModel->GetPriorPdf());
      }
      else {
         if (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters() ) {
            Error("RA2bHypoTestInvDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified");
            return 0;
         }
      }
   } 
   else 
      ((FrequentistCalculator*) hc)->SetToys(ntoys,ntoys); 

   // Get the result
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);


   TStopwatch tw; tw.Start(); 

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

   // -------------------------------------------------------
   // 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);*/
   /* ((ToyMCSampler*) fc.GetTestStatSampler())->SetGenerateBinned(true); */
   ToyMCSampler*  toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
   ProfileLikelihoodTestStat* testStat = dynamic_cast<ProfileLikelihoodTestStat*>(toymcsampler->GetTestStatistic());
   testStat->SetOneSided(true);

   // 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
   // However, the test statistic has to be installed on the workers
   // so either turn off PROOF or include the modified test statistic
   // in your `$ROOTSYS/roofit/roostats/inc` directory,
   // add the additional line to the LinkDef.h file,
   // and recompile root.
   if (useProof) {
      ProofConfig pc(*w, nworkers, "", false);
      toymcsampler->SetProofConfig(&pc); // enable proof
   }

   if(mc->GetGlobalObservables()){
      cout << "will use global observables for unconditional ensemble"<<endl;
      mc->GetGlobalObservables()->Print();
      toymcsampler->SetGlobalObservables(*mc->GetGlobalObservables());

//.........这里部分代码省略.........
  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
  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);
  // ((ToyMCSampler*) fc.GetTestStatSampler())->SetGenerateBinned(true);
  ToyMCSampler*  toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler(); 
  ProfileLikelihoodTestStat* testStat = dynamic_cast<ProfileLikelihoodTestStat*>(toymcsampler->GetTestStatistic());
  testStat->SetOneSided(true);


  // test speedups:
  testStat->SetReuseNLL(true);
  //  toymcsampler->setUseMultiGen(true); // not fully validated

  // 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, 4, "",false); 
  if(mc->GetGlobalObservables()){
    cout << "will use global observables for unconditional ensemble"<<endl;
    mc->GetGlobalObservables()->Print();
    toymcsampler->SetGlobalObservables(*mc->GetGlobalObservables());
  }
  toymcsampler->SetProofConfig(&pc);	// enable proof

void statTest(double mu_pe, double mu_hyp, ModelConfig *mc , RooDataSet *data ){

    int nToyMC = 5;
    // set roofit seed
    RooRandom::randomGenerator()->SetSeed();

    cout << endl;
    cout << endl;
    cout << "Will generate " << nToyMC << " pseudo-experiments for : " << endl;
    cout << " - mu[pseudo-data] = " << mu_pe  << endl;
    cout << " - mu[stat-test]   = " << mu_hyp << endl;
    cout << endl;

    // Check number of POI (for Wald approx)
    RooArgSet *ParamOfInterest = (RooArgSet*) mc->GetParametersOfInterest();
    int nPOI = ParamOfInterest->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");
    }
    RooRealVar* firstPOI    = (RooRealVar*) ParamOfInterest->first(); 
    RooAbsPdf *simPdf = (mc->GetPdf());
    //PrintAllParametersAndValues( *mc->GetGlobalObservables() );
    //PrintAllParametersAndValues( *mc->GetObservables() );
    firstPOI->setVal(0.0); // FIXME

    //simPdf->fitTo( *data, Hesse(kTRUE), Minos(kTRUE), PrintLevel(1) );
    simPdf->fitTo( *data );

    // set up the sampler
    ToyMCSampler sampler;
    sampler.SetPdf(*mc->GetPdf());
    sampler.SetObservables(*mc->GetObservables());
    sampler.SetNToys(nToyMC);
    sampler.SetGlobalObservables(*mc->GetGlobalObservables());
    sampler.SetParametersForTestStat(*mc->GetParametersOfInterest());
    RooArgSet* poiset = dynamic_cast<RooArgSet*>(ParamOfInterest->Clone());


    // only unconditional fit
    MinNLLTestStat *minNll = new MinNLLTestStat(*mc->GetPdf());
    minNll->EnableDetailedOutput(true);
    sampler.AddTestStatistic(minNll);

    // enable PROOF if desired
    //ProofConfig pc(*w, 8, "workers=8", kFALSE);
    //sampler.SetProofConfig(&pc);

    // evaluate the test statistics - this is where most of our time will be spent
    cout << "Generating " << nToyMC << " toys...this will take a few minutes" << endl;
    TStopwatch *mn_t = new TStopwatch; 
    mn_t->Start();
    RooDataSet* sd = sampler.GetSamplingDistributions(*poiset);
    cout << "Toy generation complete :" << endl;
    // stop timing
    mn_t->Stop();
    cout << " total CPU time: " << mn_t->CpuTime() << endl;
    cout << " total real time: " << mn_t->RealTime() << endl; 

    // now sd contains all information about our test statistics, including detailed output
    // we might eg. want to explore the results either directly, or first converting to a TTree
    // do the conversion
    TFile f("mytoys.root", "RECREATE");
    TTree *toyTree = RooStats::GetAsTTree("toyTree", "TTree created from test statistics", *sd);
    // save result to file, but in general do whatever you like
    f.cd();
    toyTree->Write();
    f.Close();
/*
    TFile* tmpFile = new TFile("mytoys.root","READ");
    TTree* myTree = (TTree*)tmpFile->Get("toyTree");

    // get boundaries for histograms
    TIter nextLeaf( (myTree->GetListOfLeaves())->MakeIterator() );
    TObject* leafObj(0);
    map<TString, float> xMaxs;
    map<TString, float> xMins;
    for(int i(0); i<myTree->GetEntries(); i++) {
      myTree->GetEntry(i);
      nextLeaf = ( (myTree->GetListOfLeaves())->MakeIterator() );
      while( (leafObj = nextLeaf.Next()) ) {
        TString name(leafObj->GetName());
        float value(myTree->GetLeaf( leafObj->GetName() )->GetValue());
        if(value > xMaxs[name]) { xMaxs[name] = value; }
        if(value < xMins[name]) { xMins[name] = value; }
      } // loop over leaves
    } // loop over tree entries

    // plot everything in the tree
    myTree->GetEntry(0);
    nextLeaf = ( (myTree->GetListOfLeaves())->MakeIterator() );
    leafObj = 0;
    // make a histogram per leaf
    map<TString, TH1F*> hists;
    myTree->GetEntry(0);
    while( (leafObj = nextLeaf.Next()) ) {
      if(!leafObj) { continue; }
      //cout << leafObj->GetName() << endl;
      TString name(leafObj->GetName());
      // special ones : fit related things
//.........这里部分代码省略.........

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

   ToyMCSampler * sampler = (ToyMCSampler *)hypoCalc->GetTestStatSampler();

   if (sampler && (calcType == 0 || calcType == 1) ) { 

      // look if pdf is number counting or extended
      if (sbModel->GetPdf()->canBeExtended() ) { 
         if (useNC)   Warning("StandardHypoTestDemo","Pdf is extended: but number counting flag is set: ignore it ");
      }
      else {
         // for not extended pdf
         if (!useNC)  { 
            int nEvents = data->numEntries();
            Info("StandardHypoTestDemo","Pdf is not extended: number of events to generate taken  from observed data set is %d",nEvents);
            sampler->SetNEventsPerToy(nEvents);
         }
         else {
            Info("StandardHypoTestDemo","using a number counting pdf");
            sampler->SetNEventsPerToy(1);
         }
      }
      
      if (data->isWeighted() && !generateBinned) { 
         Info("StandardHypoTestDemo","Data set is weighted, nentries = %d and sum of weights = %8.1f but toy generation is unbinned - it would be faster to set generateBinned to true\n",data->numEntries(), data->sumEntries());
      }
      if (generateBinned)  sampler->SetGenerateBinned(generateBinned);


      // set the test statistic
      if (testStatType == 0) sampler->SetTestStatistic(slrts); 
      if (testStatType == 1) sampler->SetTestStatistic(ropl); 
      if (testStatType == 2 || testStatType == 3) sampler->SetTestStatistic(profll); 

   }

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 HypoTestInvDemo(const char * fileName ="GausModel_b.root",
                     const char * wsName = "w",
                     const char * modelSBName = "model_sb",
                     const char * modelBName = "model_b",
                     const char * dataName = "data_obs",                  
                     int type = 0,  // calculator type 
                     int testStatType = 0, // test stat type
                     int npoints = 10,   
                     int ntoys=1000,
                     bool useCls = true )
{ 
   /*
    type = 0 Freq calculator 
    type = 1 Hybrid 

    testStatType = 0 LEP
                 = 1 Tevatron 
                 = 2 PL


   */

   if (fileName==0) { 
      std::cout << "give input filename " << std::endl;
      return; 
   }
   TFile * file = new TFile(fileName); 

   RooWorkspace * w = dynamic_cast<RooWorkspace*>( file->Get(wsName) );
   if (!w) {      
      return; 
   }
   w->Print();


   RooAbsData * data = w->data(dataName); 
   if (!data) { 
      Error("HypoTestDemo","Not existing data %s",dataName);
   }

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


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

   RatioOfProfiledLikelihoodsTestStat 
   ropl(*bModel->GetPdf(), *sbModel->GetPdf(), sbModel->GetSnapshot());
   ropl.SetSubtractMLE(false);
   
   ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   profll.SetOneSided(0);

   TestStatistic * testStat = &slrts;
   if (testStatType == 1) testStat = &ropl;
   if (testStatType == 2) testStat = &profll;
  
   
   HypoTestCalculatorGeneric *  hc = 0;
   if (type == 0) hc = new FrequentistCalculator(*data, *sbModel, *bModel);
   else new HybridCalculator(*data, *sbModel, *bModel);

   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   //toymcs->SetNEventsPerToy(1);
   toymcs->SetTestStatistic(testStat);


    if (type == 1) { 
      HybridCalculator *hhc = (HybridCalculator*) hc;
      hhc->SetToys(ntoys,ntoys); 
      // hhc->ForcePriorNuisanceAlt(*pdfNuis);
      // hhc->ForcePriorNuisanceNull(*pdfNuis);
   } 
   else 
      ((FrequentistCalculator*) hc)->SetToys(ntoys,ntoys); 

  // Get the result
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);


   TStopwatch tw; tw.Start(); 
   const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
   RooRealVar *poi = (RooRealVar*)poiSet->first();

   // fit the data first
   sbModel->GetPdf()->fitTo(*data);
   double poihat  = poi->getVal();
   //poi->setVal(30);
   //poi->setError(10);


   HypoTestInverter calc(*hc);
   // GENA: for two-sided interval
   //calc.SetConfidenceLevel(0.95);
   // GENA: for 95% upper limit
//.........这里部分代码省略.........