C++ RooFitResult::minNll方法代码示例

float getDLL(RooWorkspace* w, TString tag) {
  //RooFitResult::sexp_b_fitres

  RooFitResult *b = (RooFitResult*)w->obj(tag+"_b_fitres");
  RooFitResult *sb = (RooFitResult*)w->obj(tag+"_sb_fitres");

  return -2*(sb->minNll()-b->minNll());
}

void fit(Int_t i, Double_t va=-0.43, Double_t vb=0.){
    sprintf(namestr,"%d_%.2f_%.2f",i,va,vb);

    Double_t g1, g2;
    g1 = 0.1*r.Rndm() - 0.05;
    g2 =     r.Rndm() - 0.5;

    //a->setVal(va);
    //b->setVal(vb);

    //G001->setVal(g1);
    //G002->setVal(g2);

    //RooRealVar * G000 = new RooRealVar("G000", "G000",  0.5);
    //RooRealVar * G001 = new RooRealVar("G001", "G001",  g1, -5., 5.);
    //RooRealVar * G002 = new RooRealVar("G002", "G002",  g2, -5., 5.);
    //RooRealVar * G003 = new RooRealVar("G003", "G003", 0.0);//, -1., 1.);
    //RooRealVar * G004 = new RooRealVar("G004", "G004", 0.0);//, -1., 1.);
    //RooRealVar * a    = new RooRealVar("a",    "a",    va);
    //RooRealVar * b    = new RooRealVar("b",    "b",    vb);
    //RooRealVar * n    = new RooRealVar("n",    "n",    1000., -100., 5000.);
    RooRealVar G000("G000", "G000",  0.5);
    RooRealVar G001("G001", "G001",  g1, -5., 5.);
    RooRealVar G002("G002", "G002",  g2, -5., 5.);
    RooRealVar G003("G003", "G003", 0.0);//, -1., 1.);
    RooRealVar G004("G004", "G004", 0.0);//, -1., 1.);
    RooRealVar    a("a",    "a",    va);
    RooRealVar    b("b",    "b",    vb);
    RooRealVar    n("n",    "n",    1000., -100., 5000.);

    RooB2Kll pdf("pdf", "pdf", *cosTheta, G000, G001, G002, G003, G004, a, b, n);

    RooFitResult * fitresult = pdf.fitTo(*data,Save(kTRUE), Minos(kFALSE), NumCPU(4), SumW2Error(kTRUE));
    
    RooPlot * frame = cosTheta->frame();
    data->plotOn(frame);
    pdf.plotOn(frame);

    if(fitresult->minNll() == fitresult->minNll() && fitresult->minNll()>0 ) {
    fout << i << "\t" << a.getVal() << "\t" << b.getVal() << "\t" << fitresult->minNll() << "\t" << fitresult->status() << "\t" 
         << G000.getVal() << "\t" << G001.getVal() << "\t" << G002.getVal() << "\t" << G003.getVal() << "\t" << G004.getVal() << endl;
    }

    gROOT->ProcessLine(".x ~/lhcb/lhcbStyle.C");
    TCanvas c("fit","fit", 800, 800);
    frame->Draw();
    c.SaveAs("fits/fit"+TString(namestr)+".png");
    c.SaveAs("fits/fit"+TString(namestr)+".pdf");
}

void runFit(RooAbsPdf *pdf, RooDataSet *data, double *NLL, int *stat_t, int MaxTries, int mhLow, int mhHigh){

	int ntries=0;
	int stat=1;
	double minnll=10e8;
	while (stat!=0){
	  if (ntries>=MaxTries) break;
	  RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),Range(mhLow,mhHigh));
	  //RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),Range(85,110));
	  //RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),SumW2Error(kTRUE)
          stat = fitTest->status();
	  minnll = fitTest->minNll();
	  ntries++; 
	}
	*stat_t = stat;
	*NLL = minnll;
}

void FitBias(TString CAT,TString CUT,float SIG,float BKG,int NTOYS)
{
  gROOT->ForceStyle();
  
  RooMsgService::instance().setSilentMode(kTRUE);
  RooMsgService::instance().setStreamStatus(0,kFALSE);
  RooMsgService::instance().setStreamStatus(1,kFALSE);
  
  // -----------------------------------------
  TFile *fTemplates = TFile::Open("templates_"+CUT+"_"+CAT+"_workspace.root");
  RooWorkspace *wTemplates = (RooWorkspace*)fTemplates->Get("w");
  RooRealVar *x            = (RooRealVar*)wTemplates->var("mTop");
  RooAbsPdf *pdf_signal    = (RooAbsPdf*)wTemplates->pdf("ttbar_pdf_Nominal");
  RooAbsPdf *pdf_bkg       = (RooAbsPdf*)wTemplates->pdf("qcdCor_pdf"); 
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  x->setBins(250);   
  RooPlot *frame;

  TFile *outf;

  if (NTOYS > 1) { 
    outf = TFile::Open("FitBiasToys_"+CUT+"_"+CAT+".root","RECREATE");
  }

  float nSigInj,nBkgInj,nSigFit,nBkgFit,eSigFit,eBkgFit,nll;

  TTree *tr = new TTree("toys","toys");
  
  tr->Branch("nSigInj",&nSigInj,"nSigInj/F");
  tr->Branch("nSigFit",&nSigFit,"nSigFit/F");
  tr->Branch("nBkgInj",&nBkgInj,"nBkgInj/F");
  tr->Branch("nBkgFit",&nBkgFit,"nBkgFit/F");
  tr->Branch("eSigFit",&eSigFit,"eSigFit/F");
  tr->Branch("eBkgFit",&eBkgFit,"eBkgFit/F");
  tr->Branch("nll"    ,&nll    ,"nll/F");

  for(int itoy=0;itoy<NTOYS;itoy++) {
    // generate pseudodataset
    nSigInj = rnd->Poisson(SIG);
    nBkgInj = rnd->Poisson(BKG);
    RooRealVar *nSig = new RooRealVar("nSig","nSig",nSigInj);
    RooRealVar *nBkg = new RooRealVar("nBkg","nBkg",nBkgInj);
    RooAddPdf *model = new RooAddPdf("model","model",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nSig,*nBkg)); 
    RooDataSet *data = model->generate(*x,nSigInj+nBkgInj);
    
    RooDataHist *roohist = new RooDataHist("roohist","roohist",RooArgList(*x),*data);
    // build fit model
    RooRealVar *nFitSig = new RooRealVar("nFitSig","nFitSig",SIG,0,10*SIG);
    RooRealVar *nFitBkg = new RooRealVar("nFitBkg","nFitBkg",BKG,0,10*BKG);
    RooAddPdf *modelFit = new RooAddPdf("modelFit","modelFit",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nFitSig,*nFitBkg)); 
    // fit the pseudo dataset
    RooFitResult *res = modelFit->fitTo(*roohist,RooFit::Save(),RooFit::Extended(kTRUE));
    //res->Print();
    nSigFit = nFitSig->getVal();
    nBkgFit = nFitBkg->getVal();
    eSigFit = nFitSig->getError();
    eBkgFit = nFitBkg->getError();
    nll     = res->minNll();
    tr->Fill();
    if (itoy % 100 == 0) {
      cout<<"Toy #"<<itoy<<": injected = "<<nSigInj<<", fitted = "<<nSigFit<<", error = "<<eSigFit<<endl;
    }
    if (NTOYS == 1) {
      frame = x->frame();
      roohist->plotOn(frame); 
      model->plotOn(frame);
    }
  }
  if (NTOYS == 1) {
    TCanvas *can = new TCanvas("Toy","Toy",900,600);
    frame->Draw();
  }  
  else {
    outf->cd();
    tr->Write();
    outf->Close();
    fTemplates->Close();
  }  
}

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

    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()) {
      minAIC.setVal(AIC.getVal());
    }
    //aicExpSum+=TMath::Exp(-0.5*AIC.getVal()); //we will need this precomputed  for the next step
    pdfs_bonly.add(*p_bonly);
    pdfs_b.add(*p_b);
  }
  ws->import(minAIC);
  //compute the AIC weight
  float aicExpSum=0;
  for(auto t = tags.begin(); t!=tags.end(); t++) {
    RooFormulaVar *AIC = (RooFormulaVar*)ws->obj(*t+"_b_AIC");
    aicExpSum+=TMath::Exp(-0.5*(AIC->getVal()-minAIC.getVal())); //we will need this precomputed  for the next step    
  }
  std::cout << "aicExpSum: " << aicExpSum << std::endl;

  for(auto t = tags.begin(); t!=tags.end(); t++) {
    RooFormulaVar *AIC = (RooFormulaVar*)ws->obj(*t+"_b_AIC");
    RooRealVar *AICw = new RooRealVar(*t+"_b_AICWeight","",TMath::Exp(-0.5*(AIC->getVal()-minAIC.getVal()))/aicExpSum);
    if( TMath::IsNaN(AICw->getVal()) ) {AICw->setVal(0);}
    ws->import(*AICw);
    std::cout << *t << ":  " << AIC->getVal()-minAIC.getVal() << "    " << AICw->getVal() << std::endl;
    weights.add(*AICw);
  }
  RooAddPdf bonly_AIC("bonly_AIC","",pdfs_bonly,weights);
  RooAddPdf b_AIC("b_AIC","",pdfs_b,weights);

  //b_AIC.fitTo(data,RooFit::Strategy(0),RooFit::Extended(kTRUE),RooFit::Range("sideband_low,sideband_high"));
  //RooFitResult* bres = b_AIC.fitTo(data,RooFit::Strategy(2),RooFit::Save(kTRUE),RooFit::Extended(kTRUE),RooFit::Range("sideband_low,sideband_high"));

//.........这里部分代码省略.........
	}
	RooDataSet *startPars = new RooDataSet("startParsForce", "startParsForce", *w->set(parsName));
	startPars->add(*w->set(parsName));

	// set up parameters and ranges
	RooArgList *varyPars = new RooArgList();
	TIterator* it = w->set(parsName)->createIterator();
	while ( RooRealVar* p = (RooRealVar*)it->Next() )
	{
		if ( p->isConstant() ) continue;
		if ( forceVariables=="" && ( false
					|| TString(p->GetName()).BeginsWith("d") ///< use these variables
					// || TString(p->GetName()).BeginsWith("r")
					|| TString(p->GetName()).BeginsWith("k")
					|| TString(p->GetName()) == "g"
					) && ! (
						TString(p->GetName()) == "rD_k3pi"  ///< don't use these
						|| TString(p->GetName()) == "rD_kpi"
						// || TString(p->GetName()) == "dD_kpi"
						|| TString(p->GetName()) == "d_dk"
						|| TString(p->GetName()) == "d_dsk"
						))
		{
			varyPars->add(*p);
		}
		else if ( forceVariables.Contains(TString(p->GetName())+",") )
		{
			varyPars->add(*p);
		}
	}
	delete it;
	int nPars = varyPars->getSize();
	if ( debug ) cout << "Utils::fitToMinForce() : nPars = " << nPars << " => " << pow(2.,nPars) << " fits" << endl;
	if ( debug ) cout << "Utils::fitToMinForce() : varying ";
	if ( debug ) varyPars->Print();

	//////////

	r = fitToMinBringBackAngles(w->pdf(pdfName), false, printlevel);

	//////////

	int nErrors = 0;

	// We define a binary mask where each bit corresponds
	// to parameter at max or at min.
	for ( int i=0; i<pow(2.,nPars); i++ )
	{
		if ( debug ) cout << "Utils::fitToMinForce() : fit " << i << "        \r" << flush;
		setParameters(w, parsName, startPars->get(0));

		for ( int ip=0; ip<nPars; ip++ )
		{
			RooRealVar *p = (RooRealVar*)varyPars->at(ip);
			float oldMin = p->getMin();
			float oldMax = p->getMax();
			setLimit(w, p->GetName(), "force");
			if ( i/(int)pow(2.,ip) % 2==0 ) { p->setVal(p->getMin()); }
			if ( i/(int)pow(2.,ip) % 2==1 ) { p->setVal(p->getMax()); }
			p->setRange(oldMin, oldMax);
		}

		// check if start parameters are sensible, skip if they're not
		double startParChi2 = getChi2(w->pdf(pdfName));
		if ( startParChi2>2000 ){
			nErrors += 1;
			continue;
		}

		// refit
		RooFitResult *r2 = fitToMinBringBackAngles(w->pdf(pdfName), false, printlevel);

		// In case the initial fit failed, accept the second one.
		// If both failed, still select the second one and hope the
		// next fit succeeds.
		if ( !(r->edm()<1 && r->covQual()==3) ){
			delete r;
			r = r2;
		}
		else if ( r2->edm()<1 && r2->covQual()==3 && r2->minNll()<r->minNll() ){
			// better minimum found!
			delete r;
			r = r2;
		}
		else{
			delete r2;
		}
	}

	if ( debug ) cout << endl;
	if ( debug ) cout << "Utils::fitToMinForce() : nErrors = " << nErrors << endl;

	RooMsgService::instance().setGlobalKillBelow(INFO);

	// (re)set to best parameters
	setParameters(w, parsName, r);

	delete startPars;
	return r;
}

int main(int argc, char *argv[]){
 
  OptionParser(argc,argv);

  TStopwatch sw;
  sw.Start();

  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
  RooMsgService::instance().setSilentMode(true);

  system("mkdir -p plots/fTest");

  vector<string> procs;
  split(procs,procString_,boost::is_any_of(","));
  
  TFile *inFile = TFile::Open(filename_.c_str());
  RooWorkspace *inWS = (RooWorkspace*)inFile->Get("cms_hgg_workspace");
  
  RooRealVar *mass = (RooRealVar*)inWS->var("CMS_hgg_mass");
  //mass->setBins(320);
  //mass->setRange(mass_-10,mass_+10);
  //mass->setBins(20);
  RooRealVar *MH = new RooRealVar("MH","MH",mass_);
  MH->setVal(mass_);
  MH->setConstant(true);

  map<string,pair<int,int> > choices;
  map<string,vector<RooPlot*> > plotsRV;
  map<string,vector<RooPlot*> > plotsWV;

  for (unsigned int p=0; p<procs.size(); p++){
    vector<RooPlot*> tempRV;
    vector<RooPlot*> tempWV;
    for (int cat=0; cat<ncats_; cat++){
      RooPlot *plotRV = mass->frame(Range(mass_-10,mass_+10));
      plotRV->SetTitle(Form("%s_cat%d_RV",procs[p].c_str(),cat));
      tempRV.push_back(plotRV);
      RooPlot *plotWV = mass->frame(Range(mass_-10,mass_+10));
      plotWV->SetTitle(Form("%s_cat%d_WV",procs[p].c_str(),cat));
      tempWV.push_back(plotWV);
    }
    plotsRV.insert(pair<string,vector<RooPlot*> >(procs[p],tempRV));
    plotsWV.insert(pair<string,vector<RooPlot*> >(procs[p],tempWV));
  }

  vector<int> colors;
  colors.push_back(kBlue);
  colors.push_back(kRed);
  colors.push_back(kGreen+2);
  colors.push_back(kMagenta+1);
  
  for (int cat=0; cat<ncats_; cat++){
    for (unsigned int p=0; p<procs.size(); p++){
      string proc = procs[p];
      RooDataSet *dataRV = (RooDataSet*)inWS->data(Form("sig_%s_mass_m%d_rv_cat%d",proc.c_str(),mass_,cat));
      RooDataSet *dataWV = (RooDataSet*)inWS->data(Form("sig_%s_mass_m%d_wv_cat%d",proc.c_str(),mass_,cat));
      //mass->setBins(160);
      //RooDataHist *dataRV = dataRVtemp->binnedClone();
      //RooDataHist *dataWV = dataWVtemp->binnedClone();
      //RooDataSet *dataRVw = (RooDataSet*)dataRVtemp->reduce(Form("CMS_hgg_mass>=%3d && CMS_hgg_mass<=%3d",mass_-10,mass_+10)); 
      //RooDataSet *dataWVw = (RooDataSet*)dataWVtemp->reduce(Form("CMS_hgg_mass>=%3d && CMS_hgg_mass<=%3d",mass_-10,mass_+10));
      //RooDataHist *dataRV = new RooDataHist(Form("roohist_%s",dataRVtemp->GetName()),Form("roohist_%s",dataRVtemp->GetName()),RooArgSet(*mass),*dataRVtemp);
      //RooDataHist *dataWV = new RooDataHist(Form("roohist_%s",dataWVtemp->GetName()),Form("roohist_%s",dataWVtemp->GetName()),RooArgSet(*mass),*dataWVtemp);
      //RooDataSet *dataRV = stripWeights(dataRVweight,mass);
      //RooDataSet *dataWV = stripWeights(dataWVweight,mass);
      //RooDataSet *data = (RooDataSet*)inWS->data(Form("sig_%s_mass_m%d_cat%d",proc.c_str(),mass_,cat));
     
      int rvChoice=0;
      int wvChoice=0;
      
      // right vertex
      int order=1;
      int prev_order=0;
      int cache_order=0;
      double thisNll=0.;
      double prevNll=1.e6;
      double chi2=0.;
      double prob=0.;

      dataRV->plotOn(plotsRV[proc][cat]);
      //while (prob<0.8) {
      while (order<5) {
        RooAddPdf *pdf = buildSumOfGaussians(Form("cat%d_g%d",cat,order),mass,MH,order);
        RooFitResult *fitRes = pdf->fitTo(*dataRV,Save(true),SumW2Error(true));//,Range(mass_-10,mass_+10));
        double myNll=0.;
        thisNll = fitRes->minNll();
        //double myNll = getMyNLL(mass,pdf,dataRV);
        //thisNll = getMyNLL(mass,pdf,dataRV);
        //RooAbsReal *nll = pdf->createNLL(*dataRV);
        //RooMinuit m(*nll);
        //m.migrad();
        //thisNll = nll->getVal();
        //plot(Form("plots/fTest/%s_cat%d_g%d_rv",proc.c_str(),cat,order),mass_,mass,dataRV,pdf);
        pdf->plotOn(plotsRV[proc][cat],LineColor(colors[order-1]));
        chi2 = 2.*(prevNll-thisNll);
        //if (chi2<0. && order>1) chi2=0.;
        int diffInDof = (2*order+(order-1))-(2*prev_order+(prev_order-1));
        prob = TMath::Prob(chi2,diffInDof);
        cout << "\t RV: " << cat << " " << order << " " << diffInDof << " " << prevNll << " " << thisNll << " " << myNll << " " << chi2 << " " << prob << endl;
        prevNll=thisNll;
//.........这里部分代码省略.........

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


       if ( poiMinVal < 0. && poiMaxVal < 0. ) {

          printf("\n\n Automatic determination of scan range.\n\n") ;

          if ( startPoiVal <= 0. ) {
             printf("\n\n *** POI starting value zero or negative %g.  Quit.\n\n\n", startPoiVal ) ;
             return ;
          }

          poiMinVal = startPoiVal - 3.5 * sqrt(startPoiVal) ;
          poiMaxVal = startPoiVal + 6.0 * sqrt(startPoiVal) ;

          if ( poiMinVal < 0. ) { poiMinVal = 0. ; }

          printf("    Start val = %g.   Scan range:   %g  to  %g\n\n", startPoiVal, poiMinVal, poiMaxVal ) ;


       }



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


       double poiVals_scanDown[1000] ;
       double nllVals_scanDown[1000] ;

       //-- Do scan down from best value.

       printf("\n\n +++++ Starting scan down from best value.\n\n") ;

       double minNllVal(1.e9) ;

       for ( int poivi=0; poivi < npoiPoints/2 ; poivi++ ) {

          ////double poiValue = poiMinVal + poivi*(poiMaxVal-poiMinVal)/(1.*(npoiPoints-1)) ;
          double poiValue = best_poi_val - poivi*(best_poi_val-poiMinVal)/(1.*(npoiPoints/2-1)) ;

          rrv_poiValue -> setVal( poiValue ) ;
          rrv_poiValue -> setConstant( kTRUE ) ;


       //+++++++++++++++++++++++++++++++++++

          rminuit->migrad() ;
          rminuit->hesse() ;
          RooFitResult* rfr = rminuit->save() ;

       //+++++++++++++++++++++++++++++++++++


          if ( verbLevel > 0 ) { rfr->Print("v") ; }


          float fit_minuit_var_val = rfr->minNll() ;

          printf(" %02d : poi constraint = %.2f : allvars : MinuitVar, createNLL, PV, POI :    %.5f   %.5f   %.5f   %.5f\n",
                poivi, rrv_poiValue->getVal(), fit_minuit_var_val, nll->getVal(), plot_var->getVal(), new_poi_rar->getVal() ) ;
          cout << flush ;



          poiVals_scanDown[poivi] = new_poi_rar->getVal() ;
          nllVals_scanDown[poivi] = plot_var->getVal() ;

///
/// Perform the 1d Prob scan.
/// Saves chi2 values and the prob-Scan p-values in a root tree
/// For the datasets stuff, we do not yet have a MethodDatasetsProbScan class, so we do it all in
/// MethodDatasetsProbScan
/// \param nRun Part of the root tree file name to facilitate parallel production.
///
int MethodDatasetsProbScan::scan1d(bool fast, bool reverse)
{
	if (fast) return 0; // tmp

	if ( arg->debug ) cout << "MethodDatasetsProbScan::scan1d() : starting ... " << endl;

    // Set limit to all parameters.
    this->loadParameterLimits(); /// Default is "free", if not changed by cmd-line parameter


    // Define scan parameter and scan range.
    RooRealVar *parameterToScan = w->var(scanVar1);
    float parameterToScan_min = hCL->GetXaxis()->GetXmin();
    float parameterToScan_max = hCL->GetXaxis()->GetXmax();

		// do a free fit
		RooFitResult *result = this->loadAndFit(this->pdf); // fit on data
		assert(result);
    RooSlimFitResult *slimresult = new RooSlimFitResult(result,true);
		slimresult->setConfirmed(true);
		solutions.push_back(slimresult);
		double freeDataFitValue = w->var(scanVar1)->getVal();

    // Define outputfile
    system("mkdir -p root");
    TString probResName = Form("root/scan1dDatasetsProb_" + this->pdf->getName() + "_%ip" + "_" + scanVar1 + ".root", arg->npoints1d);
    TFile* outputFile = new TFile(probResName, "RECREATE");

    // Set up toy root tree
    this->probScanTree = new ToyTree(this->pdf, arg);
    this->probScanTree->init();
    this->probScanTree->nrun = -999; //\todo: why does this branch even exist in the output tree of the prob scan?

    // Save parameter values that were active at function
    // call. We'll reset them at the end to be transparent
    // to the outside.
    RooDataSet* parsFunctionCall = new RooDataSet("parsFunctionCall", "parsFunctionCall", *w->set(pdf->getParName()));
    parsFunctionCall->add(*w->set(pdf->getParName()));

    // start scan
    cout << "MethodDatasetsProbScan::scan1d_prob() : starting ... with " << nPoints1d << " scanpoints..." << endl;
    ProgressBar progressBar(arg, nPoints1d);
    for ( int i = 0; i < nPoints1d; i++ )
    {
        progressBar.progress();
        // scanpoint is calculated using min, max, which are the hCL x-Axis limits set in this->initScan()
        // this uses the "scan" range, as expected
        // don't add half the bin size. try to solve this within plotting method

        float scanpoint = parameterToScan_min + (parameterToScan_max - parameterToScan_min) * (double)i / ((double)nPoints1d - 1);
				if (arg->debug) cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() " << scanpoint << " " << parameterToScan_min << " " << parameterToScan_max << endl;

        this->probScanTree->scanpoint = scanpoint;

        if (arg->debug) cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() - scanpoint in step " << i << " : " << scanpoint << endl;

        // don't scan in unphysical region
        // by default this means checking against "free" range
        if ( scanpoint < parameterToScan->getMin() || scanpoint > parameterToScan->getMax() + 2e-13 ) {
            cout << "it seems we are scanning in an unphysical region: " << scanpoint << " < " << parameterToScan->getMin() << " or " << scanpoint << " > " << parameterToScan->getMax() + 2e-13 << endl;
            exit(EXIT_FAILURE);
        }

        // FIT TO REAL DATA WITH FIXED HYPOTHESIS(=SCANPOINT).
        // THIS GIVES THE NUMERATOR FOR THE PROFILE LIKELIHOOD AT THE GIVEN HYPOTHESIS
        // THE RESULTING NUISANCE PARAMETERS TOGETHER WITH THE GIVEN HYPOTHESIS ARE ALSO
        // USED WHEN SIMULATING THE TOY DATA FOR THE FELDMAN-COUSINS METHOD FOR THIS HYPOTHESIS(=SCANPOINT)
        // Here the scanvar has to be fixed -> this is done once per scanpoint
        // and provides the scanner with the DeltaChi2 for the data as reference
        // additionally the nuisances are set to the resulting fit values

        parameterToScan->setVal(scanpoint);
        parameterToScan->setConstant(true);

        RooFitResult *result = this->loadAndFit(this->pdf); // fit on data
        assert(result);

        if (arg->debug) {
            cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() - minNll data scan at scan point " << scanpoint << " : " << 2 * result->minNll() << ": "<< 2 * pdf->getMinNll() << endl;
        }
        this->probScanTree->statusScanData = result->status();

        // set chi2 of fixed fit: scan fit on data
        // CAVEAT: chi2min from fitresult gives incompatible results to chi2min from pdf
        // this->probScanTree->chi2min           = 2 * result->minNll();
        this->probScanTree->chi2min           = 2 * pdf->getMinNll();
        this->probScanTree->covQualScanData   = result->covQual();
        this->probScanTree->scanbest  = freeDataFitValue;

        // After doing the fit with the parameter of interest constrained to the scanpoint,
        // we are now saving the fit values of the nuisance parameters. These values will be
        // used to generate toys according to the PLUGIN method.
        this->probScanTree->storeParsScan(); // \todo : figure out which one of these is semantically the right one
//.........这里部分代码省略.........

   void ws_constrained_profile3D( const char* wsfile = "rootfiles/ws-data-unblind.root",
                                   const char* new_poi_name = "n_M234_H4_3b",
                                   int npoiPoints = 20,
                                   double poiMinVal = 0.,
                                   double poiMaxVal = 20.,
                                   double constraintWidth = 1.5,
                                   double ymax = 10.,
                                   int verbLevel=0 ) {


     gStyle->SetOptStat(0) ;

     //--- make output directory.

     char command[10000] ;
     sprintf( command, "basename %s", wsfile ) ;
     TString wsfilenopath = gSystem->GetFromPipe( command ) ;
     wsfilenopath.ReplaceAll(".root","") ;
     char outputdirstr[1000] ;
     sprintf( outputdirstr, "outputfiles/scans-%s", wsfilenopath.Data() ) ;
     TString outputdir( outputdirstr ) ;


     printf("\n\n Creating output directory: %s\n\n", outputdir.Data() ) ;
     sprintf(command, "mkdir -p %s", outputdir.Data() ) ;
     gSystem->Exec( command ) ;


     //--- Tell RooFit to shut up about anything less important than an ERROR.
      RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;



       if ( verbLevel > 0 ) { printf("\n\n Verbose level : %d\n\n", verbLevel) ; }


       TFile* wstf = new TFile( wsfile ) ;

       RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );

       if ( verbLevel > 0 ) { ws->Print() ; }






       RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;

       if ( verbLevel > 0 ) {
          printf("\n\n\n  ===== RooDataSet ====================\n\n") ;
          rds->Print() ;
          rds->printMultiline(cout, 1, kTRUE, "") ;
       }





       ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
       RooAbsPdf* likelihood = modelConfig->GetPdf() ;

       RooRealVar* rrv_mu_susy_all0lep = ws->var("mu_susy_all0lep") ;
       if ( rrv_mu_susy_all0lep == 0x0 ) {
          printf("\n\n\n *** can't find mu_susy_all0lep in workspace.  Quitting.\n\n\n") ;
          return ;
       }





       //-- do BG only.
       rrv_mu_susy_all0lep->setVal(0.) ;
       rrv_mu_susy_all0lep->setConstant( kTRUE ) ;










       //-- do a prefit.

       printf("\n\n\n ====== Pre fit with unmodified nll var.\n\n") ;

       RooFitResult* dataFitResultSusyFixed = likelihood->fitTo(*rds, Save(true),Hesse(false),Minos(false),Strategy(1),PrintLevel(verbLevel));
       int dataSusyFixedFitCovQual = dataFitResultSusyFixed->covQual() ;
       if ( dataSusyFixedFitCovQual < 2 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFixedFitCovQual ) ; return ; }
       double dataFitSusyFixedNll = dataFitResultSusyFixed->minNll() ;

       if ( verbLevel > 0 ) {
          dataFitResultSusyFixed->Print("v") ;
       }

       printf("\n\n Nll value, from fit result : %.3f\n\n", dataFitSusyFixedNll ) ;

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

//.........这里部分代码省略.........
			if ( scanUp )
			{
				scanvalue = min + (max-min)*(double)i/(double)nPoints1d + hCL->GetBinWidth(1)/2.;
				if ( scanvalue < scanStart ) continue;
				if ( scanvalue > scanStop ) break;
			}
			else
			{
				scanvalue = max - (max-min)*(double)(i+1)/(double)nPoints1d + hCL->GetBinWidth(1)/2.;
				if ( scanvalue > scanStart ) continue;
				if ( scanvalue < scanStop ) break;
			}

			// disable drag mode
			// (the improve method doesn't work with drag mode as parameter run
			// at their limits)
			if ( scanDisableDragMode ) setParameters(w, parsName, startPars->get(0));

			// set the parameter of interest to the scan point
			par->setVal(scanvalue);

			// don't scan in unphysical region
			if ( scanvalue < par->getMin() || scanvalue > par->getMax() ) continue;

			// status bar
			if ( (((int)nStep % (int)(nTotalSteps/printFreq)) == 0))
				cout << "MethodProbScan::scan1d() : scanning " << (float)nStep/(float)nTotalSteps*100. << "%   \r" << flush;

			// fit!
			RooFitResult *fr = 0;
			if ( arg->probforce )         fr = fitToMinForce(w, combiner->getPdfName());
			else if ( arg->probimprove )  fr = fitToMinImprove(w, combiner->getPdfName());
			else                          fr = fitToMinBringBackAngles(w->pdf(pdfName), false, -1);
			double chi2minScan = fr->minNll();
			if ( std::isinf(chi2minScan) ) chi2minScan=1e4; // else the toys in PDF_testConstraint don't work
			RooSlimFitResult *r = new RooSlimFitResult(fr); // try to save memory by using the slim fit result
			delete fr;
			allResults.push_back(r);
			bestMinFoundInScan = TMath::Min((double)chi2minScan, (double)bestMinFoundInScan);

			TString warningChi2Neg;
			if ( chi2minScan < 0 ){
				float newChi2minScan = chi2minGlobal + 25.; // 5sigma more than best point
				warningChi2Neg = "MethodProbScan::scan1d() : WARNING : " + title;
				warningChi2Neg += TString(Form(" chi2 negative for scan point %i: %f",i,chi2minScan));
				warningChi2Neg += " setting to: " + TString(Form("%f",newChi2minScan));
				//cout << warningChi2Neg << "\r" << flush;
				cout << warningChi2Neg << endl;
				chi2minScan = newChi2minScan;
			}

			// If we find a minimum smaller than the old "global" minimum, this means that all
			// previous 1-CL values are too high.
			if ( chi2minScan<chi2minGlobal ){
				if ( arg->verbose ) cout << "MethodProbScan::scan1d() : WARNING : '" << title << "' new global minimum found! "
																		<< " chi2minScan=" << chi2minScan << endl;
				chi2minGlobal = chi2minScan;
				// recompute previous 1-CL values
				for ( int k=1; k<=hCL->GetNbinsX(); k++ ){
					hCL->SetBinContent(k, TMath::Prob(hChi2min->GetBinContent(k)-chi2minGlobal, 1));
				}
			}

			double deltaChi2 = chi2minScan - chi2minGlobal;
			double oneMinusCL = TMath::Prob(deltaChi2, 1);

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

				// set start parameters from inner turn of the spiral
				int xStartPars, yStartPars;
				computeInnerTurnCoords(iStart, jStart, i, j, xStartPars, yStartPars, 1);
				RooSlimFitResult *rStartPars = mycurveResults2d[xStartPars-1][yStartPars-1];
				if ( rStartPars ) setParameters(w, parsName, rStartPars);

				// memory management:
				tMemory.Start(false);
				// delete old, inner fit results, that we don't need for start parameters anymore
				// for this we take the second-inner-most turn.
				int iOld, jOld;
				bool innerTurnExists = computeInnerTurnCoords(iStart, jStart, i, j, iOld, jOld, 2);
				if ( innerTurnExists ){
					deleteIfNotInCurveResults2d(mycurveResults2d[iOld-1][jOld-1]);
					mycurveResults2d[iOld-1][jOld-1] = 0;
				}
				tMemory.Stop();

				// alternative choice for start parameters: always from what we found at function call
				// setParameters(w, parsName, startPars->get(0));

				// set scan point
				float scanvalue1 = hCL2d->GetXaxis()->GetBinCenter(i);
				float scanvalue2 = hCL2d->GetYaxis()->GetBinCenter(j);
				par1->setVal(scanvalue1);
				par2->setVal(scanvalue2);

				// fit!
				tFit.Start(false);
				RooFitResult *fr;
				if ( !arg->probforce ) fr = fitToMinBringBackAngles(w->pdf(pdfName), false, -1);
				else                   fr = fitToMinForce(w, combiner->getPdfName());
				double chi2minScan = fr->minNll();
				tFit.Stop();
				tSlimResult.Start(false);
				RooSlimFitResult *r = new RooSlimFitResult(fr); // try to save memory by using the slim fit result
				tSlimResult.Stop();
				delete fr;
				allResults.push_back(r);
				bestMinFoundInScan = TMath::Min((double)chi2minScan, (double)bestMinFoundInScan);
				mycurveResults2d[i-1][j-1] = r;

				// If we find a new global minumum, this means that all
				// previous 1-CL values are too high. We'll save the new possible solution, adjust the global
				// minimum, return a status code, and stop.
				if ( chi2minScan > -500 && chi2minScan<chi2minGlobal ){
					// warn only if there was a significant improvement
					if ( arg->debug || chi2minScan<chi2minGlobal-1e-2 ){
						if ( arg->verbose ) cout << "MethodProbScan::scan2d() : WARNING : '" << title << "' new global minimum found! chi2minGlobal="
															<< chi2minGlobal << " chi2minScan=" << chi2minScan << endl;
					}
					chi2minGlobal = chi2minScan;
					// recompute previous 1-CL values
					for ( int k=1; k<=hCL2d->GetNbinsX(); k++ )
						for ( int l=1; l<=hCL2d->GetNbinsY(); l++ ){
							hCL2d->SetBinContent(k, l, TMath::Prob(hChi2min2d->GetBinContent(k,l)-chi2minGlobal, ndof));
						}
				}

				double deltaChi2 = chi2minScan - chi2minGlobal;
				double oneMinusCL = TMath::Prob(deltaChi2, ndof);

				// Save the 1-CL value. But only if better than before!
				if ( hCL2d->GetBinContent(i, j) < oneMinusCL ){
					hCL2d->SetBinContent(i, j, oneMinusCL);

//.........这里部分代码省略.........
	        if( 0 )
	        {
	            // it's strange ,if you can't draw following six figs in the big 3*2 pad,you can fist use these code,run once,don't exit the ROOT,then you can draw six figs!!!
			    data = new RooDataHist("data", "data", x, time2lastshowermuon[0]);
				sum  = new RooAddPdf("sum","sum",RooArgList(expLi9, expHe8, expIbd),RooArgList(coeLi9, coeHe8, NIbd));
		        RooPlot* mesframe = x.frame() ;
		        data->plotOn(mesframe) ;
		        c->cd(1);
		        mesframe->Draw();
	            return true;
	        }
        }
        for( int ihist=0 ; ihist<1 ; ihist++ )
        {
			double minNl[6]={0.};
			for( int j=0 ; j<6 ; j++ )
			{
	            std::cout<<"now is  "<< j+1<<endl;
	            NumIbd[j]=hh[j]->Integral(1,hh[j]->FindBin(40));
	            hh[j]->Rebin(100);
				data = new RooDataHist("data", "data", x, hh[j]);
				sum  = new RooAddPdf("sum","sum",RooArgList(expLi9, expHe8, expIbd),RooArgList(coeLi9, coeHe8, NIbd));
				for( int i=0 ; i<41 ; i++ )
				{
					Li9Frac.setVal(i*0.025);
					fitres = sum->fitTo((*data),Save(),PrintLevel(-1),Extended(kTRUE));
	                //fitres->Print();
					n98[i]=N98.getVal(0);
					in98[i]=N98.getError();
	                nIbd[i]=NIbd.getVal(0);
	                inIbd[i]=NIbd.getError();
	                rateMuValue[i]=rateMu.getVal(0);
	                rateMuErr[i]=rateMu.getError();
					minNLL[i] = fitres->minNll();
					if( minNLL[i]<minNl[j] )
					{
						minNl[j]=minNLL[i];
						minindex[j]=i;
					}
				}
				n98fit[j]=n98[minindex[j]];
				in98fit[j]=in98[minindex[j]];
	            rateMufit[j]=rateMuValue[minindex[j]];
	            if( draw6slice )
	            {
		            //for drawing fit figure
					Li9Frac.setVal(minindex[j]*0.025);
					fitres = sum->fitTo((*data),Save(),PrintLevel(-1),Extended(kTRUE));
		            //fitres->Print();
		            mesframe[j] = x.frame() ;
		            data->plotOn(mesframe[j]) ;
		            sum->plotOn(mesframe[j]);
		            sum->plotOn(mesframe[j],Components(expIbd),LineStyle(kDashed),LineColor(kGreen)) ;
		            sum->plotOn(mesframe[j],Components(RooArgSet(expLi9,expHe8)),LineStyle(kDashed),LineColor(kRed)) ;
		            xTitle[j]+="        time since last muon (s)";
		            mesframe[j]->GetXaxis()->SetTitle(xTitle[j]);
		            mesframe[j]->GetYaxis()->SetTitle("Entries");
		            c->cd(j+1);
		            mesframe[j]->Draw();
		            //gPad->SetLogy();
	            }
			}
	        
			n98total=(n98fit[0]+n98fit[1]+n98fit[2]+(n98fit[3]+n98fit[4]+n98fit[5])*exp(-1/0.257))/0.678;
			ifitn98total=sqrt(in98fit[0]*in98fit[0]+in98fit[1]*in98fit[1]+in98fit[2]*in98fit[2]+((in98fit[3]+in98fit[4]+in98fit[5])*exp(-1/0.257)*(in98fit[3]+in98fit[4]+in98fit[5])*exp(-1/0.257)))/0.678;
	        if( ihist==0 )

float ComputeTestStat(TString wsfile, double mu_susy_sig_val) {

  gROOT->Reset();

  TFile* wstf = new TFile( wsfile ) ;

  RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );
  ws->Print() ;
  
  ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
  
  modelConfig->Print() ;

  RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;
  
  rds->Print() ;
  rds->printMultiline(cout, 1, kTRUE, "") ;
  
  RooAbsPdf* likelihood = modelConfig->GetPdf() ;
  
  RooRealVar* rrv_mu_susy_sig = ws->var("mu_susy_all0lep") ;
  if ( rrv_mu_susy_sig == 0x0 ) {
    printf("\n\n\n *** can't find mu_susy_all0lep in workspace.  Quitting.\n\n\n") ;
    return ;
  } else {
    printf(" current value is : %8.3f\n", rrv_mu_susy_sig->getVal() ) ; cout << flush ;
    rrv_mu_susy_sig->setConstant(kFALSE) ;
  }

  /*
  // check the impact of varying the qcd normalization:

  RooRealVar *rrv_qcd_0lepLDP_ratioH1 = ws->var("qcd_0lepLDP_ratio_H1");
  RooRealVar *rrv_qcd_0lepLDP_ratioH2 = ws->var("qcd_0lepLDP_ratio_H2");
  RooRealVar *rrv_qcd_0lepLDP_ratioH3 = ws->var("qcd_0lepLDP_ratio_H3");
  
  rrv_qcd_0lepLDP_ratioH1->setVal(0.3);
  rrv_qcd_0lepLDP_ratioH2->setVal(0.3);
  rrv_qcd_0lepLDP_ratioH3->setVal(0.3);
  
  rrv_qcd_0lepLDP_ratioH1->setConstant(kTRUE);
  rrv_qcd_0lepLDP_ratioH2->setConstant(kTRUE);
  rrv_qcd_0lepLDP_ratioH3->setConstant(kTRUE);
  */
  
  printf("\n\n\n  ===== Doing a fit with SUSY component floating ====================\n\n") ;

  RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
  double logLikelihoodSusyFloat = fitResult->minNll() ;
  
  double logLikelihoodSusyFixed(0.) ;
  double testStatVal(-1.) ;
  if ( mu_susy_sig_val >= 0. ) {
    printf("\n\n\n  ===== Doing a fit with SUSY fixed ====================\n\n") ;
    printf(" fixing mu_susy_sig to %8.2f.\n", mu_susy_sig_val ) ;
    rrv_mu_susy_sig->setVal( mu_susy_sig_val ) ;
    rrv_mu_susy_sig->setConstant(kTRUE) ;
    
    fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
    logLikelihoodSusyFixed = fitResult->minNll() ;
    testStatVal = 2.*(logLikelihoodSusyFixed - logLikelihoodSusyFloat) ;
    printf("\n\n\n ======= test statistic : -2 * ln (L_fixed / ln L_max) = %8.3f\n\n\n", testStatVal ) ;
  }


  return testStatVal ;

}

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



       printf("\n\n\n === Model values for observables\n\n") ;

       printObservables() ;



      //--- save the actual values of the observables.

       saveObservables() ;











       //--- evaluate the test stat on the data: fit with susy floating.

       rrv_mu_susy_sig->setVal( poiVal ) ;
       rrv_mu_susy_sig->setConstant( kTRUE ) ;

       printf("\n\n\n ====== Fitting the data with susy fixed.\n\n") ;

       RooFitResult* dataFitResultSusyFixed = likelihood->fitTo(*rds, Save(true));
       int dataSusyFixedFitCovQual = dataFitResultSusyFixed->covQual() ;
       if ( dataSusyFixedFitCovQual != 3 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFixedFitCovQual ) ; return ; }
       double dataFitSusyFixedNll = dataFitResultSusyFixed->minNll() ;


       rrv_mu_susy_sig->setVal( 0.0 ) ;
       rrv_mu_susy_sig->setConstant( kFALSE ) ;

       printf("\n\n\n ====== Fitting the data with susy floating.\n\n") ;

       RooFitResult* dataFitResultSusyFloat = likelihood->fitTo(*rds, Save(true));
       int dataSusyFloatFitCovQual = dataFitResultSusyFloat->covQual() ;
       if ( dataSusyFloatFitCovQual != 3 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFloatFitCovQual ) ; return ; }
       double dataFitSusyFloatNll = dataFitResultSusyFloat->minNll() ;

       double dataTestStat = 2.*( dataFitSusyFixedNll - dataFitSusyFloatNll) ;

       printf("\n\n\n Data value of test stat : %8.2f\n", dataTestStat ) ;












       printf("\n\n\n === Nuisance parameters\n\n") ;

       {
          int npi(0) ;
          TIterator* npIter = nuisanceParameters->createIterator() ;