C++ RooAbsReal::createProfile方法代码示例

//.........这里部分代码省略.........
   //  fixed_again.add( *ws1->var("nsig2_pp") );
   // fixed_again.add( *ws1->var("sigma1") );
   //  fixed_again.add( *ws1->var("nbkg_pp") );
   // fixed_again.add( *ws1->var("npow") );
   // fixed_again.add( *ws1->var("muPlusPt") );
   // fixed_again.add( *ws1->var("muMinusPt") );
   // fixed_again.add( *ws1->var("mscale_hi") );
   // fixed_again.add( *ws1->var("mscale_pp") );
   //  
   // ws1->Print();
   cout << "99999" << endl;

   // create signal+background Model Config
   RooStats::ModelConfig sbHypo("SbHypo");
   sbHypo.SetWorkspace( *ws1 );
   sbHypo.SetPdf( *ws1->pdf("joint") );
   sbHypo.SetObservables( obs );
   sbHypo.SetGlobalObservables( globalObs );
   sbHypo.SetParametersOfInterest( poi );
   sbHypo.SetNuisanceParameters( nuis );
   sbHypo.SetPriorPdf( *ws1->pdf("step") ); // this is optional

   // ws1->Print();
   /////////////////////////////////////////////////////////////////////
   RooAbsReal * pNll = sbHypo.GetPdf()->createNLL( *data,NumCPU(10) );
   cout << "111111" << endl;
   RooMinuit(*pNll).migrad(); // minimize likelihood wrt all parameters before making plots
   cout << "444444" << endl;
   RooPlot *framepoi = ((RooRealVar *)poi.first())->frame(Bins(10),Range(0.,0.2),Title("LL and profileLL in raa3"));
   cout << "222222" << endl;
   pNll->plotOn(framepoi,ShiftToZero());
   cout << "333333" << endl;
   
   RooAbsReal * pProfile = pNll->createProfile( globalObs ); // do not profile global observables
   pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
   pProfile->plotOn(framepoi,LineColor(kRed));
   framepoi->SetMinimum(0);
   framepoi->SetMaximum(3);
   TCanvas *cpoi = new TCanvas();
   cpoi->cd(); framepoi->Draw();
   cpoi->SaveAs("cpoi.pdf");

   ((RooRealVar *)poi.first())->setMin(0.);
   RooArgSet * pPoiAndNuisance = new RooArgSet("poiAndNuisance");
   // pPoiAndNuisance->add(*sbHypo.GetNuisanceParameters());
   // pPoiAndNuisance->add(*sbHypo.GetParametersOfInterest());
   pPoiAndNuisance->add( nuis );
   pPoiAndNuisance->add( poi );
   sbHypo.SetSnapshot(*pPoiAndNuisance);

   RooPlot* xframeSB = pObs->frame(Title("SBhypo"));
   data->plotOn(xframeSB,Cut("dataCat==dataCat::hi"));
   RooAbsPdf *pdfSB = sbHypo.GetPdf();
   RooCategory *dataCat = ws1->cat("dataCat");
   pdfSB->plotOn(xframeSB,Slice(*dataCat,"hi"),ProjWData(*dataCat,*data));
   TCanvas *c1 = new TCanvas();
   c1->cd(); xframeSB->Draw();
   c1->SaveAs("c1.pdf");

   delete pProfile;
   delete pNll;
   delete pPoiAndNuisance;
   ws1->import( sbHypo );
   /////////////////////////////////////////////////////////////////////
   RooStats::ModelConfig bHypo = sbHypo;
   bHypo.SetName("BHypo");

void LL(){

  //y0 = 0.000135096401209 sigma_y0 = 0.000103896581837 x0 = 0.000446013873443 sigma_x0 =1.81384394011e-06
  //0.014108652249 0.0168368471049 0.0219755396247 0.000120423865262 1.5575931164 1.55759310722 3.41637854038
  //0.072569437325 0.084063541977 0.0376693978906 0.000284216132439 0.51908074913 0.519080758095 1.12037749267
 // double d = 0.014108652249;
 //  double sd = 0.0168368471049;
 //  double mc = 0.0219755396247;
 //  double smc = 0.000120423865262;
 //  double r0 = d/mc;

  double d = 0.072569437325;
  double sd =  0.084063541977;
  double mc =  0.0376693978906;
  double smc =  0.00028421613243;
  double r0 = d/mc;

  RooRealVar x("x","x",mc*0.9,mc*1.1);
  RooRealVar x0("x0","x0",mc);
  RooRealVar sx("sx","sx",smc);

  RooRealVar r("r","r",r0,0.,5.);
  RooRealVar y0("y0","y0",d); 
  RooRealVar sy("sy","sy",sd); 
  
  RooProduct rx("rx","rx",RooArgList(r,x));

  RooGaussian g1("g1","g1",x,x0,sx);
  RooGaussian g2("g2","g2",rx,y0,sy);

  RooProdPdf LL("LL","LL",g1,g2);

  RooArgSet obs(x0,y0); //observables
  RooArgSet poi(r); //parameters of interest
  RooDataSet data("data", "data", obs);
  data.add(obs); //actually add the data


  RooFitResult* res = LL.fitTo(data,RooFit::Minos(poi),RooFit::Save(),RooFit::Hesse(false));
  if(res->status()==0) {
    r.Print();
    x.Print();
    cout << r.getErrorLo() << " " << r.getErrorHi() << endl;
  } else {
    cout << "Likelihood maximization failed" << endl;
  }
  
  RooAbsReal* nll = LL.createNLL(data); 
  RooPlot* frame = r.frame();
  RooAbsReal* pll = nll->createProfile(poi);
  pll->plotOn(frame);//,RooFit::LineColor(ROOT::kRed));
  frame->Draw();

  r.setVal(0.);
  cout << pll->getVal() << endl; 

  return;
    
    


}

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

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

//.........这里部分代码省略.........
    RooGaussian HiggsMassConstraint("HiggsMassConstraint","",mu,HiggsMass,HiggsMassError);


    RooAddPdf fitModel(tag+"_sb_model","",RooArgList( *ws->pdf("b_"+tag), sig ),RooArgList(Nbkg,Nsig));

    RooFitResult* sbres;
    RooAbsReal* nll;
    if(constrainMu) {
      fitModel.fitTo(data,RooFit::Strategy(0),RooFit::Extended(kTRUE),RooFit::ExternalConstraints(RooArgSet(HiggsMassConstraint)));
      sbres = fitModel.fitTo(data,RooFit::Strategy(2),RooFit::Save(kTRUE),RooFit::Extended(kTRUE),RooFit::ExternalConstraints(RooArgSet(HiggsMassConstraint)));
      nll = fitModel.createNLL(data,RooFit::NumCPU(4),RooFit::Extended(kTRUE),RooFit::ExternalConstraints(RooArgSet(HiggsMassConstraint)));
    } else {
      fitModel.fitTo(data,RooFit::Strategy(0),RooFit::Extended(kTRUE));
      sbres = fitModel.fitTo(data,RooFit::Strategy(2),RooFit::Save(kTRUE),RooFit::Extended(kTRUE));
      nll = fitModel.createNLL(data,RooFit::NumCPU(4),RooFit::Extended(kTRUE));
    }
    sbres->SetName(tag+"_sb_fitres");
    ws->import(*sbres);
    ws->import(fitModel);

    RooPlot *fmgg = mgg.frame();
    data.plotOn(fmgg);
    fitModel.plotOn(fmgg);
    ws->pdf("b_"+tag+"_ext")->plotOn(fmgg,RooFit::LineColor(kRed),RooFit::Range("Full"),RooFit::NormRange("Full"));
    fmgg->SetName(tag+"_frame");
    ws->import(*fmgg);
    delete fmgg;


    RooMinuit(*nll).migrad();

    RooPlot *fNs = Nsig.frame(0,25);
    fNs->SetName(tag+"_Nsig_pll");
    RooAbsReal *pll = nll->createProfile(Nsig);
    //nll->plotOn(fNs,RooFit::ShiftToZero(),RooFit::LineColor(kRed));
    pll->plotOn(fNs);
    ws->import(*fNs);

    delete fNs;

    RooPlot *fmu = mu.frame(125,132);
    fmu->SetName(tag+"_mu_pll");
    RooAbsReal *pll_mu = nll->createProfile(mu);
    pll_mu->plotOn(fmu);
    ws->import(*fmu);

    delete fmu;

  }

  RooArgSet weights("weights");
  RooArgSet pdfs_bonly("pdfs_bonly");
  RooArgSet pdfs_b("pdfs_b");

  RooRealVar minAIC("minAIC","",1E10);
  //compute AIC stuff
  for(auto t = tags.begin(); t!=tags.end(); t++) {
    RooAbsPdf *p_bonly = ws->pdf("bonly_"+*t);
    RooAbsPdf *p_b = ws->pdf("b_"+*t);
    RooFitResult *sb = (RooFitResult*)ws->obj(*t+"_bonly_fitres");
    RooRealVar k(*t+"_b_k","",p_bonly->getParameters(RooArgSet(mgg))->getSize());
    RooRealVar nll(*t+"_b_minNll","",sb->minNll());
    RooRealVar Npts(*t+"_b_N","",blind_data->sumEntries());
    RooFormulaVar AIC(*t+"_b_AIC","2*@0+2*@1+2*@1*(@1+1)/(@[email protected])",RooArgSet(nll,k,Npts));
    ws->import(AIC);
    if(AIC.getVal() < minAIC.getVal()) {

//.........这里部分代码省略.........
   RooDataSet* parameterScan = (RooDataSet*) fc.GetPointsToScan();
   RooArgSet* tmpPoint;

   // make a histogram of parameter vs. threshold
   TH1F* histOfThresholds = new TH1F("histOfThresholds","",
                                       parameterScan->numEntries(),
                                       firstPOI->getMin(),
                                       firstPOI->getMax());
   histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
   histOfThresholds->GetYaxis()->SetTitle("Threshold");

   // loop through the points that were tested and ask confidence belt
   // what the upper/lower thresholds were.
   // For FeldmanCousins, the lower cut off is always 0
   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

//.........这里部分代码省略.........
    TIterator* iter = paramSet2->createIterator();
    TObject* var = iter->Next();
    while (var != 0) {
        counter++;
        ParamName = var->GetName();
        vParam = w->var(ParamName);
        ParamValue = vParam->getVal();
        ParamError = vParam->getError();
        ParamLimitLow = vParam->getMin();
        ParamLimitHigh = vParam->getMax();
        cout << ParamName << " has value " << ParamValue << " with error: " << ParamError << " and limits: " << ParamLimitLow << " to " << ParamLimitHigh << endl << endl;

        if (ParamError == 0) {  //Skipping fixed parameters
            cout << "Parameter was fixed, skipping its fitting" << endl;
            cout << endl << "DONE WITH " << counter << " PARAMETER OUT OF " << Nparams << endl << endl;
            var = iter->Next();
            continue;
        }

        // determining fit range: Nsigma sigma on each side unless it would be outside of parameter limits
        if ((ParamValue - Nsigma * ParamError) > ParamLimitLow) {
            FitRangeLow = (ParamValue - Nsigma * ParamError);
        }
        else {
            FitRangeLow = ParamLimitLow;
        }

        if ((ParamValue + Nsigma * ParamError) < ParamLimitHigh) {
            FitRangeHigh = (ParamValue + Nsigma * ParamError);
        }
        else {
            FitRangeHigh = ParamLimitHigh;
        }


        // P l o t    p l a i n   l i k e l i h o o d   a n d   C o n s t r u c t   p r o f i l e   l i k e l i h o o d
        // ---------------------------------------------------
        RooPlot* frame1;
        RooAbsReal* pll=NULL;

        if (Nbins != 0) {
            frame1 = vParam->frame(Bins(Nbins), Range(FitRangeLow, FitRangeHigh), Title(TString::Format("LL and profileLL in %s", ParamName.Data())));
            nll->plotOn(frame1, ShiftToZero());

            pll = nll->createProfile(*vParam);
            // Plot the profile likelihood
            pll->plotOn(frame1, LineColor(kRed), RooFit::Precision(-1));
        }
        else { //Skip profile likelihood
            frame1 = vParam->frame(Bins(10), Range(FitRangeLow, FitRangeHigh), Title(TString::Format("LL and profileLL in %s", ParamName.Data())));
            nll->plotOn(frame1, ShiftToZero());
        }

        // D r a w   a n d   s a v e   p l o t s
        // -----------------------------------------------------------------------

        // Adjust frame maximum for visual clarity
        frame1->SetMinimum(0);
        frame1->SetMaximum(20);

        TCanvas* c = new TCanvas("CLikelihoodResult", "CLikelihoodResult", 800, 600);
        c->cd(1);
        gPad->SetLeftMargin(0.15);
        frame1->GetYaxis()->SetTitleOffset(1.4);
        frame1->Draw();
        TLegend* leg = new TLegend(0.70, 0.70, 0.95, 0.88, "");
        leg->SetFillColor(kWhite);
        leg->SetBorderSize(0);
        leg->SetTextSize(0.035);
        TLegendEntry *le1 = leg->AddEntry(nll, "Plain likelihood", "l");
        le1->SetLineColor(kBlue);
        le1->SetLineWidth(3);
        TLegendEntry *le2 = leg->AddEntry(pll, "Profile likelihood", "l");
        le2->SetLineColor(kRed);
        le2->SetLineWidth(3);
        leg->Draw("same");

        //Save plot
        TString StrippedName = TString(Filename(Filename.Last('/')+1,Filename.Length()));
        StrippedName = StrippedName.ReplaceAll(".root","");
        cout << StrippedName << endl;
        gSystem->mkdir(Form("%s/root/%s", Outputdir.Data(), StrippedName.Data()), kTRUE);
        c->SaveAs(Form("%s/root/%s/Likelihood_scan_%s.root", Outputdir.Data(), StrippedName.Data(), ParamName.Data()));
        gSystem->mkdir(Form("%s/pdf/%s", Outputdir.Data(), StrippedName.Data()), kTRUE);
        c->SaveAs(Form("%s/pdf/%s/Likelihood_scan_%s.pdf", Outputdir.Data(), StrippedName.Data(), ParamName.Data()));
        gSystem->mkdir(Form("%s/png/%s", Outputdir.Data(), StrippedName.Data()), kTRUE);
        c->SaveAs(Form("%s/png/%s/Likelihood_scan_%s.png", Outputdir.Data(), StrippedName.Data(), ParamName.Data()));


        delete c;
        delete frame1;
        if (pll) delete pll;

        cout << endl << "DONE WITH " << counter << " PARAMETER OUT OF " << Nparams << endl << endl;
        //if (counter == 2){ break; } //Exit - for testing
        var = iter->Next();
    }  // End of the loop

    return true;
}

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;


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

//.........这里部分代码省略.........
  data->add( *pObs );

  // import dataset into workspace
  pWs->import(*data);

  // create set of global observables (need to be defined as constants)
  pWs->var("glob_lumi")->setConstant(true);
  pWs->var("glob_efficiency")->setConstant(true);
  pWs->var("glob_nbkg")->setConstant(true);
  RooArgSet globalObs("global_obs");
  globalObs.add( *pWs->var("glob_lumi") );
  globalObs.add( *pWs->var("glob_efficiency") );
  globalObs.add( *pWs->var("glob_nbkg") );

  // create set of parameters of interest (POI)
  RooArgSet poi("poi");
  poi.add( *pWs->var("xsec") );
  
  // create set of nuisance parameters
  RooArgSet nuis("nuis");
  nuis.add( *pWs->var("beta_lumi") );
  nuis.add( *pWs->var("beta_efficiency") );
  nuis.add( *pWs->var("beta_nbkg") );

  // create signal+background Model Config
  RooStats::ModelConfig sbHypo("SbHypo");
  sbHypo.SetWorkspace( *pWs );
  sbHypo.SetPdf( *pWs->pdf("model") );
  sbHypo.SetObservables( obs );
  sbHypo.SetGlobalObservables( globalObs );
  sbHypo.SetParametersOfInterest( poi );
  sbHypo.SetNuisanceParameters( nuis );
  sbHypo.SetPriorPdf( *pWs->pdf("prior") ); // this is optional

  // fix all other variables in model:
  // everything except observables, POI, and nuisance parameters
  // must be constant
  pWs->var("lumi_nom")->setConstant(true);
  pWs->var("efficiency_nom")->setConstant(true);
  pWs->var("nbkg_nom")->setConstant(true);
  pWs->var("lumi_kappa")->setConstant(true);
  pWs->var("efficiency_kappa")->setConstant(true);
  pWs->var("nbkg_kappa")->setConstant(true);
  RooArgSet fixed("fixed");
  fixed.add( *pWs->var("lumi_nom") );
  fixed.add( *pWs->var("efficiency_nom") );
  fixed.add( *pWs->var("nbkg_nom") );
  fixed.add( *pWs->var("lumi_kappa") );
  fixed.add( *pWs->var("efficiency_kappa") );
  fixed.add( *pWs->var("nbkg_kappa") );
  
  // set parameter snapshot that corresponds to the best fit to data
  RooAbsReal * pNll = sbHypo.GetPdf()->createNLL( *data );
  RooAbsReal * pProfile = pNll->createProfile( globalObs ); // do not profile global observables
  pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
  RooArgSet * pPoiAndNuisance = new RooArgSet("poiAndNuisance");
  pPoiAndNuisance->add(*sbHypo.GetNuisanceParameters());
  pPoiAndNuisance->add(*sbHypo.GetParametersOfInterest());
  sbHypo.SetSnapshot(*pPoiAndNuisance);
  delete pProfile;
  delete pNll;
  delete pPoiAndNuisance;

  // import S+B ModelConfig into workspace
  pWs->import( sbHypo );

  // create background-only Model Config from the S+B one
  RooStats::ModelConfig bHypo = sbHypo;
  bHypo.SetName("BHypo");
  bHypo.SetWorkspace(*pWs);

  // set parameter snapshot for bHypo, setting xsec=0
  // it is useful to understand how this block of code works
  // but you can also use it as a recipe to make a parameter snapshot
  pNll = bHypo.GetPdf()->createNLL( *data );
  RooArgSet poiAndGlobalObs("poiAndGlobalObs");
  poiAndGlobalObs.add( poi );
  poiAndGlobalObs.add( globalObs );
  pProfile = pNll->createProfile( poiAndGlobalObs ); // do not profile POI and global observables
  ((RooRealVar *)poi.first())->setVal( 0 );  // set xsec=0 here
  pProfile->getVal(); // this will do fit and set nuisance parameters to profiled values
  pPoiAndNuisance = new RooArgSet( "poiAndNuisance" );
  pPoiAndNuisance->add( nuis );
  pPoiAndNuisance->add( poi );
  bHypo.SetSnapshot(*pPoiAndNuisance);
  delete pProfile;
  delete pNll;
  delete pPoiAndNuisance;

  // import model config into workspace
  pWs->import( bHypo );

  // print out the workspace contents
  pWs->Print();

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

  return;
}

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

      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.

//.........这里部分代码省略.........
  RooDataSet* parameterScan = (RooDataSet*) fc.GetPointsToScan();
  RooArgSet* tmpPoint;

  // make a histogram of parameter vs. threshold
  TH1F* histOfThresholds = new TH1F("histOfThresholds","",
				    parameterScan->numEntries(),
				    firstPOI->getMin(),
				    firstPOI->getMax());
  histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
  histOfThresholds->GetYaxis()->SetTitle("Threshold");

  // loop through the points that were tested and ask confidence belt
  // what the upper/lower thresholds were.
  // For FeldmanCousins, the lower cut off is always 0
  for(Int_t i=0; i<parameterScan->numEntries(); ++i){
    tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
    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-constriant
  ////////////////////////////////////////////////////////////

  // 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");


  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
    w->loadSnapshot("paramsToGenerateData");

    // in 5.30 there is a nicer way to generate toy data  & randomize global obs
    RooAbsData* toyData = toymcsampler->GenerateToyData(*paramsToGenerateData);

    // get test stat at observed UL in observed data
    firstPOI->setVal(observedUL);
    double toyTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*toyData,tmpPOI);
    //    toyData->get()->Print("v");
    //    cout <<"obsTSatObsUL " <<obsTSatObsUL << "toyTS " << toyTSatObsUL << endl;

void plot_pll(TString fname="monoh_withsm_SRCR_bg11.7_bgslop-0.0_nsig0.0.root")
{
  SetAtlasStyle();



  TFile* file =  TFile::Open(fname);
  RooWorkspace* wspace = (RooWorkspace*) file->Get("wspace");

  cout << "\n\ncheck that eff and reco terms included in BSM component to make fiducial cross-section" <<endl;
  wspace->function("nsig")->Print();
  RooRealVar* reco = wspace->var("reco");
  if(  wspace->function("nsig")->dependsOn(*reco) ) {
    cout << "all good." <<endl;
  } else {
    cout << "need to rerun fit_withsm using DO_FIDUCIAL_LIMIT true" <<endl;
    return;
  }

  /*
  // DANGER
  // TEST WITH EXAGGERATED UNCERTAINTY
  wspace->var("unc_theory")->setMax(1);
  wspace->var("unc_theory")->setVal(1);
  wspace->var("unc_theory")->Print();
  */

  // this was for making plot about decoupling/recoupling approach
  TCanvas* tc = new TCanvas("tc","",400,400);
  RooPlot *frame = wspace->var("xsec_bsm")->frame();
  RooAbsPdf* pdfc = wspace->pdf("jointModeld");
  RooAbsData* data = wspace->data("data");
  RooAbsReal *nllJoint = pdfc->createNLL(*data, RooFit::Constrained()); // slice with fixed xsec_bsm
  RooAbsReal *profileJoint = nllJoint->createProfile(*wspace->var("xsec_bsm"));

  wspace->allVars().Print("v");
  pdfc->fitTo(*data);
  wspace->allVars().Print("v");
  wspace->var("xsec_bsm")->Print();
  double nllmin = 2*nllJoint->getVal();
  wspace->var("xsec_bsm")->setVal(0);
  double nll0 = 2*nllJoint->getVal();
  cout << Form("nllmin = %f, nll0 = %f, Z=%f", nllmin, nll0, sqrt(nll0-nllmin)) << endl;
  nllJoint->plotOn(frame, RooFit::LineColor(kGreen), RooFit::LineStyle(kDotted), RooFit::ShiftToZero(), RooFit::Name("nll_statonly")); // no error
  profileJoint->plotOn(frame,RooFit::Name("pll") );
  wspace->var("xsec_sm")->Print();
  wspace->var("theory")->Print();
  wspace->var("theory")->setConstant();
  profileJoint->plotOn(frame, RooFit::LineColor(kRed), RooFit::LineStyle(kDashed), RooFit::Name("pll_smfixed") );

  frame->GetXaxis()->SetTitle("#sigma_{BSM, fid} [fb]");
  frame->GetYaxis()->SetTitle("-log #lambda  ( #sigma_{BSM, fid} )");
  double temp = frame->GetYaxis()->GetTitleOffset();
  frame->GetYaxis()->SetTitleOffset( 1.1* temp );

  frame->SetMinimum(1e-7);
  frame->SetMaximum(4);


  // Legend
  double x1,y1,x2,y2;
  GetX1Y1X2Y2(tc,x1,y1,x2,y2);
  TLegend *legend_sr=FastLegend(x2-0.75,y2-0.3,x2-0.25,y2-0.5,0.045);
  legend_sr->AddEntry(frame->findObject("pll"),"with #sigma_{SM} uncertainty","L");
  legend_sr->AddEntry(frame->findObject("pll_smfixed"),"with #sigma_{SM} constant","L");
  legend_sr->AddEntry(frame->findObject("nll_statonly"),"no systematics","L");
  frame->Draw();
  legend_sr->Draw("SAME");



  // descriptive text
  vector<TString> pavetext11;
  pavetext11.push_back("#bf{#it{ATLAS Internal}}");
  pavetext11.push_back("#sqrt{#it{s}} = 8 TeV #scale[0.6]{#int}Ldt = 20.3 fb^{-1}");
  pavetext11.push_back("#it{H}+#it{E}_{T}^{miss} , #it{H #rightarrow #gamma#gamma}, #it{m}_{#it{H}} = 125.4 GeV");

  TPaveText* text11=CreatePaveText(x2-0.75,y2-0.25,x2-0.25,y2-0.05,pavetext11,0.045);
  text11->Draw();

  tc->SaveAs("pll.pdf");



  /*
  wspace->var("xsec_bsm")->setConstant(true);
  wspace->var("eff"     )->setConstant(true);
  wspace->var("mh"      )->setConstant(true);
  wspace->var("sigma_h" )->setConstant(true);
  wspace->var("lumi"    )->setConstant(true);
  wspace->var("xsec_sm" )->setVal(v_xsec_sm);
  wspace->var("eff"     )->setVal(1.0);
  wspace->var("lumi"    )->setVal(v_lumi);
  TH1* nllHist = profileJoint->createHistogram("xsec_bsm",100);
  TFile* out = new TFile("nllHist.root","REPLACE");
  nllHist->Write()
  out->Write();
  out->Close();
  */

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

void LHFit()
{
    
    double lorange = 100.;
    double hirange = 140.;
    
    // Import data
    TFile *file = new TFile("/atlas/data18a/yupan/HZZ4l2012/MiniTree/MiniTree_2013Moriond/v03/MC12a/JHUPythia8_AU2CTEQ6L1_ggH125_Spin0p_ZZ4lep.root");
    TTree *tree = (TTree*)file->Get("tree_incl_4mu");
    
    // Define variables
    float m4l = 0;
    float m34 = 0;
    float m4lerr = 0;
    float wgt = 0;
    
    // Get the number of entries in the tree
    int nevents = tree->GetEntries();
    
    tree->SetBranchStatus("*",0);
    tree->SetBranchStatus("m4l_unconstrained",1);
    tree->SetBranchStatus("mZ2_unconstrained",1);
    tree->SetBranchStatus("m4lerr_unconstrained",1);
    tree->SetBranchStatus("weight_corr",1);
    
    tree->SetBranchAddress("m4l_unconstrained",&m4l);
    tree->SetBranchAddress("mZ2_unconstrained",&m34);
    tree->SetBranchAddress("m4lerr_unconstrained",&m4lerr);
    tree->SetBranchAddress("weight_corr",&wgt);
    
    
    ///////////////
    // Build pdfs
    //////////////
    RooRealVar* mass = new RooRealVar("m4l","mass",lorange,hirange,"GeV");
    RooRealVar* mZ2 = new RooRealVar("m34","mZ2",10.,60.,"GeV");
    RooRealVar merr("m4lerr","mass err",0.1,5.0,"GeV");
    RooRealVar weight("weight","weight",0,10);

    
    float totalwgt(0);
    for (Int_t i=0; i<nevents; i++) {
        tree->GetEntry(i);
        totalwgt+=wgt;
    }
    std::cout<<"total weight = "<<totalwgt<<std::endl;
    
    
    //////////////////
    // Create ND dataset
    /////////////////
    RooDataSet signal("signal","signal",RooArgSet(*mass,*mZ2,merr,weight),"weight");
    
    std::cout<<"Reading in events from signal minitree"<<std::endl;
    for (int i=0; i<nevents; i++) {
        tree->GetEntry(i);
        mass->setVal(m4l);
        mZ2->setVal(m34);
        merr.setVal(m4lerr);
        weight.setVal(wgt/totalwgt);
        signal.add(RooArgSet(*mass,*mZ2,merr),wgt/totalwgt);
    }
    
    // Create 1D kernel estimation for signal mass
    RooKeysPdf kestmass("kestmass","kestmass",*mass,signal);
    // Create histogram of 1D kernel estimation pdf for mass
    TH1* hm = kestmass.createHistogram("hm",*mass);
    // Create mass pdf from the 1D kernel histogram
    RooDataHist* dmass = new RooDataHist("dmass","dmass",RooArgSet(*mass),hm);
    RooHistPdf* masspdf = new RooHistPdf("masspdf","masspdf",*mass,*dmass,2);
    
    
    
    //////////////////////////////
    // Construct plain likelihood
    /////////////////////////////
    
    // construct unbinned likelihood
    RooAbsReal *nll = masspdf->createNLL(signal,NumCPU(2));
    
    // Minimize likelihood w.r.t all parameters before making plots
    RooMinuit(*nll).migrad();
    
    // Plot likelihood scan mass
    RooPlot* frame1 = mass->frame(Bins(10),Title("LL and profileLL in mass"));
    
    // Construct profile likelihood in mass
    RooAbsReal *pll_mass = nll->createProfile(*mass);
    
    // Plot the profile likelihood in mass
    pll_mass->plotOn(frame1,LineColor(kRed));
    

    
    
    TCanvas *c = new TCanvas("c","c",800,600);
    gPad->SetLeftMargin(0.15);
    frame1->GetYaxis()->SetTitleOffset(1.4);
    frame1->Draw();
    c->SaveAs("testProLL.png");
//.........这里部分代码省略.........

void rf605_profilell()
{
  // C r e a t e   m o d e l   a n d   d a t a s e t 
  // -----------------------------------------------

  // Observable
  RooRealVar x("x","x",-20,20) ;

  // Model (intentional strong correlations)
  RooRealVar mean("mean","mean of g1 and g2",0,-10,10) ;
  RooRealVar sigma_g1("sigma_g1","width of g1",3) ;
  RooGaussian g1("g1","g1",x,mean,sigma_g1) ;

  RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,6.0) ;
  RooGaussian g2("g2","g2",x,mean,sigma_g2) ;

  RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
  RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;

  // Generate 1000 events
  RooDataSet* data = model.generate(x,1000) ;
  


  // C o n s t r u c t   p l a i n   l i k e l i h o o d
  // ---------------------------------------------------

  // Construct unbinned likelihood
  RooAbsReal* nll = model.createNLL(*data,NumCPU(2)) ;

  // Minimize likelihood w.r.t all parameters before making plots
  RooMinuit(*nll).migrad() ;

  // Plot likelihood scan frac 
  RooPlot* frame1 = frac.frame(Bins(10),Range(0.01,0.95),Title("LL and profileLL in frac")) ;
  nll->plotOn(frame1,ShiftToZero()) ;

  // Plot likelihood scan in sigma_g2
  RooPlot* frame2 = sigma_g2.frame(Bins(10),Range(3.3,5.0),Title("LL and profileLL in sigma_g2")) ;
  nll->plotOn(frame2,ShiftToZero()) ;



  // C o n s t r u c t   p r o f i l e   l i k e l i h o o d   i n   f r a c
  // -----------------------------------------------------------------------

  // The profile likelihood estimator on nll for frac will minimize nll w.r.t
  // all floating parameters except frac for each evaluation

  RooAbsReal* pll_frac = nll->createProfile(frac) ;

  // Plot the profile likelihood in frac
  pll_frac->plotOn(frame1,LineColor(kRed)) ;

  // Adjust frame maximum for visual clarity
  frame1->SetMinimum(0) ;
  frame1->SetMaximum(3) ;



  // C o n s t r u c t   p r o f i l e   l i k e l i h o o d   i n   s i g m a _ g 2 
  // -------------------------------------------------------------------------------

  // The profile likelihood estimator on nll for sigma_g2 will minimize nll
  // w.r.t all floating parameters except sigma_g2 for each evaluation
  RooAbsReal* pll_sigmag2 = nll->createProfile(sigma_g2) ;

  // Plot the profile likelihood in sigma_g2
  pll_sigmag2->plotOn(frame2,LineColor(kRed)) ;

  // Adjust frame maximum for visual clarity
  frame2->SetMinimum(0) ;
  frame2->SetMaximum(3) ;



  // Make canvas and draw RooPlots
  TCanvas *c = new TCanvas("rf605_profilell","rf605_profilell",800, 400);
  c->Divide(2);
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ;

  delete pll_frac ;
  delete pll_sigmag2 ;
  delete nll ;
}

Int_t Tprime::SetParameterPoints( std::string sbModelName,
                                  std::string bModelName ) {
    //
    // Fit the data with S+B model.
    // Make a snapshot of the S+B parameter point.
    // Profile with POI=0.
    // Make a snapshot of the B parameter point
    // (B model is the S+B model with POI=0
    //

    Double_t poi_value_for_b_model = 0.0;

    // get S+B model config from workspace
    RooStats::ModelConfig * pSbModel = (RooStats::ModelConfig *)pWs->obj(sbModelName.c_str());
    pSbModel->SetWorkspace(*pWs);

    // get parameter of interest set
    const RooArgSet * poi = pSbModel->GetParametersOfInterest();

    // get B model config from workspace
    RooStats::ModelConfig * pBModel = (RooStats::ModelConfig *)pWs->obj(bModelName.c_str());
    pBModel->SetWorkspace(*pWs);

    // make sure that data has been loaded
    if (!data) return -1;

    // find parameter point for global maximum with the S+B model,
    // with conditional MLEs for nuisance parameters
    // and save the parameter point snapshot in the Workspace
    RooAbsReal * nll = pSbModel->GetPdf()->createNLL(*data);
    RooAbsReal * profile = nll->createProfile(RooArgSet());
    profile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
    RooArgSet * poiAndNuisance = new RooArgSet();
    if(pSbModel->GetNuisanceParameters())
        poiAndNuisance->add(*pSbModel->GetNuisanceParameters());
    poiAndNuisance->add(*pSbModel->GetParametersOfInterest());
    pWs->defineSet("SPlusBModelParameters", *poiAndNuisance);
    pWs->saveSnapshot("SPlusBFitParameters",*poiAndNuisance);
    pSbModel->SetSnapshot(*poi);
    RooArgSet * sbModelFitParams = (RooArgSet *)poiAndNuisance->snapshot();
    cout << "\nWill save these parameter points that correspond to the fit to data" << endl;
    sbModelFitParams->Print("v");
    delete profile;
    delete nll;
    delete poiAndNuisance;
    delete sbModelFitParams;

    //

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

    return 0;
}