C++ RooRealVar::GetName方法代码示例

void Fitter::addGausFitModel(RooArgList *pdfs, RooArgList *coeffs, bool isLast)
{
    for( int i=0;i< nGaussians ;++i)
    {
        RooRealVar *m = new RooRealVar( Form("fitmodel_gaus_mean_%d",i),Form("fitmodel_gaus_mean_%d",i), startMean_[i] );
        RooRealVar *s = new RooRealVar( Form("fitmodel_gaus_sigma_%d",i), Form("fitmodel_gaus_sigma_%d",i), startSigma_[i] );
        RooGaussian *gaus = new RooGaussian(Form("fitmodel_gaus_g%d",i),Form("fitmodel_gaus_g%d",i),*x_, *m, *s);

        vars_[ m->GetName() ] = m;
        vars_[ s->GetName() ] = s;

        pdfs -> add( *gaus );

        if ( (i != nGaussians-1 )  or not isLast)
        {
            RooRealVar *f = new RooRealVar( Form("fitmodel_gaus_frac_%d",i) ,Form("fitmodel_gaus_frac_%d",i) , startFraction_[i], 0.01,1.0 );
            vars_ [ f->GetName() ] = f;
            coeffs -> add(   *f  ) ;
        }
    }

}

 void GenericModel::ImportModelData(Int_t parameter_num, RooArgList* plist)
 {
    // Import all model datasets corresponding to the defined parameters, ranges & binnings
    RooArgList PLIST;
    if (plist) PLIST.addClone(*plist);
    RooRealVar* par = dynamic_cast<RooRealVar*>(fParameters.at(parameter_num));
    RooAbsBinning& bins = par->getBinning();
    Int_t N = bins.numBins();
    RooRealVar* par_in_list = (RooRealVar*)PLIST.find(par->GetName());
    if (!par_in_list) {
       PLIST.addClone(RooRealVar(par->GetName(), par->GetTitle(), 0.));
       par_in_list = (RooRealVar*)PLIST.find(par->GetName());
    }
    for (int i = 0; i < N; i++) {
       par_in_list->setMax(bins.binHigh(i));
       par_in_list->setMin(bins.binLow(i));
       par_in_list->setVal(bins.binCenter(i));
       if ((parameter_num + 1) < GetNumberOfParameters())
          ImportModelData(parameter_num + 1, &PLIST);
       else {
          AddModelData(PLIST, GetModelDataSet(PLIST));
       }
    }
 }

double TwoBody::printMcmcUpperLimit( double peak, ModelConfig &_mc, std::string filename ){
  //
  // print out the upper limit on the first Parameter of Interest
  //

  std::string _legend = "[TwoBody::printMcmcUpperLimit]: ";

  char buf[1024];

  double _limit = numeric_limits<double>::max();

  if (mcInt){
    //mc.SetWorkspace(*ws);
    RooRealVar * firstPOI = (RooRealVar*) _mc.GetParametersOfInterest()->first();
    _limit = mcInt->UpperLimit(*firstPOI);
    std::cout << "\n95% upper limit on " << firstPOI->GetName() << " is : " <<
      _limit << endl;
    
    if (bMcmcConverged){
      sprintf(buf, "%7.1f   %7.6f", peak, _limit);
    }
    else{
      sprintf(buf, "# %7.1f   %7.6f  # MCMC did not converge", peak, _limit);
    }

  }
  else{
    sprintf(buf, "# MCMC did not converge");
  }

  if (filename.size()!=0){
    
    std::ofstream aFile;
    
    // append to file if exists
    aFile.open(filename.c_str(), std::ios_base::app);
    
    aFile << buf << std::endl;
    
    // close outfile here so it is safe even if subsequent iterations crash
    aFile.close();
    
  }
    
  return _limit;
}

///
/// Set all observables to 'toy' values drawn from the
/// PDF using the current parameter values. A certain number
/// of toys is pregenerated to speed up when doing mulitple toy fits.
///
void PDF_Abs::setObservablesToy()
{
	obsValSource = "toy";
	if( !pdf ){ cout<< "PDF_Abs::setObservables(): ERROR: pdf not initialized."<<endl; exit(1); }
	if ( toyObservables==0 || iToyObs==nToyObs )
	{
		RooRandom::randomGenerator()->SetSeed(0);
		if ( iToyObs==nToyObs ) delete toyObservables;
		toyObservables = pdf->generate(*(RooArgSet*)observables, nToyObs);
		iToyObs=0;
	}
	for ( int i=0; i<nObs; i++ )
	{
		RooRealVar* pObs = (RooRealVar*)((RooArgList*)observables)->at(i);
		pObs->setVal(((RooRealVar*)toyObservables->get(iToyObs)->find(pObs->GetName()))->getVal());
	}
	iToyObs+=1;
}

TTree *dataset2tree(RooDataSet *dataset){

  const RooArgSet *args = dataset->get();
  if(args==NULL) return NULL;
  RooArgList argList(*args);

  Double_t variables[50];
  Long64_t nEntries= dataset->numEntries();
  //nEntries=1;
  TTree *tree = new TTree("tree","tree");
  tree->SetDirectory(0);
  TIterator *it1=NULL; 
  it1 = argList.createIterator();
  int index1=0;
  for(RooRealVar *var = (RooRealVar *) it1->Next(); var!=NULL;
      var = (RooRealVar *) it1->Next(),index1++){
    TString name(var->GetName());
    name.ReplaceAll("-","_");
    tree->Branch(name, &(variables[index1]), name+"/D");
  }

  //  tree->Print();

  for(Long64_t jentry=0; jentry<nEntries; jentry++){
    
    TIterator *it=NULL; 
    RooArgList argList1(*(dataset->get(jentry)));
    it = argList1.createIterator();
    //(dataset->get(jentry))->Print();
    int index=0;
    for(RooRealVar *var = (RooRealVar *) it->Next(); var!=NULL;
	var = (RooRealVar *) it->Next(), index++){
      variables[index]=var->getVal();
      //var->Print();
    }
   
    delete it;
    tree->Fill();
  }
  tree->ResetBranchAddresses();
  //  tree->Print();
  return tree;
}

double getChisq(RooAbsData &dat, RooAbsPdf &pdf, RooRealVar &var, bool prt=false) {

    // Find total number of events
    double nEvt;
    double nTot=0.0;

    for(int j=0; j<dat.numEntries(); j++) {
        dat.get(j);
        nEvt=dat.weight();
        nTot+=nEvt;
    }

    // Find chi-squared equivalent 2NLL
    //RooRealVar *var=(RooRealVar*)(pdf.getParameters(*dat)->find("CMS_hgg_mass"));
    double totNLL=0.0;
    double prbSum=0.0;

    for(int j=0; j<dat.numEntries(); j++) {
        double m=dat.get(j)->getRealValue(var.GetName());
        if ( m < var.getMin() || m > var.getMax())  continue;
        // Find probability density and hence probability
        var.setVal(m);
        double prb = var.getBinWidth(0)*pdf.getVal(var);
        prbSum+=prb;

        dat.get(j);
        nEvt=dat.weight();

        double mubin=nTot*prb;
        double contrib(0.);
        if (nEvt < 1) contrib = mubin;
        else contrib=mubin-nEvt+nEvt*log(nEvt/mubin);
        totNLL+=contrib;

        if(prt) cout << "Bin " << j << " prob = " << prb << " nEvt = " << nEvt << ", mu = " << mubin << " contribution " << contrib << endl;
    }

    totNLL*=2.0;
    if(prt) cout << pdf.GetName() << " nTot = " << nTot << " 2NLL constant = " << totNLL << endl;

    return totNLL;
}

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


  // Now get the interval
  PointSetInterval* interval = fc.GetInterval();
  ConfidenceBelt* belt = fc.GetConfidenceBelt();
 
  // print out the iterval on the first Parameter of Interest
  cout << "\n95% interval on " <<firstPOI->GetName()<<" is : ["<<
    interval->LowerLimit(*firstPOI) << ", "<<
    interval->UpperLimit(*firstPOI) <<"] "<<endl;

  // get observed UL and value of test statistic evaluated there
  RooArgSet tmpPOI(*firstPOI);
  double observedUL = interval->UpperLimit(*firstPOI);
  firstPOI->setVal(observedUL);
  double obsTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*data,tmpPOI);


  // Ask the calculator which points were scanned
  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();

void THSEventsPDF::adjustBinning(Int_t* offset1) const
{
   RooRealVar* xvar = fx_off ;
  if (!dynamic_cast<RooRealVar*>(xvar)) {
    coutE(InputArguments) << "RooDataHist::adjustBinning(" << GetName() << ") ERROR: dimension " << xvar->GetName() << " must be real" << endl ;
    assert(0) ;
  }
  Double_t xlo = xvar->getMin() ;
  Double_t xhi = xvar->getMax() ;
  //adjust bin range limits with new scale parameter
  //cout<<scale<<" "<<fMean<<" "<<xlo<<" "<<xhi<<endl;
  xlo=(xlo-fMean)/scale+fMean;
  xhi=(xhi-fMean)/scale+fMean;
  if(xvar->getBinning().lowBound()==xlo&&xvar->getBinning().highBound()==xhi) return;
  xvar->setRange(xlo,xhi) ;
  // Int_t xmin(0) ;
  // cout<<"THSEventsPDF::adjustBinning( "<<xlo <<" "<<xhi<<endl;
  //now adjust fitting range to bin limits??Possibly not
  if (fRHist->GetXaxis()->GetXbins()->GetArray()) {

    RooBinning xbins(fRHist->GetNbinsX(),fRHist->GetXaxis()->GetXbins()->GetArray()) ;

    Double_t tolerance = 1e-6*xbins.averageBinWidth() ;
    
    // Adjust xlo/xhi to nearest boundary
    Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+tolerance)) ;
    Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-tolerance)) ;
    xbins.setRange(xloAdj,xhiAdj) ;

    xvar->setBinning(xbins) ;
    if (fabs(xloAdj-xlo)>tolerance||fabs(xhiAdj-xhi)<tolerance) {
      coutI(DataHandling) << "RooDataHist::adjustBinning(" << GetName() << "): fit range of variable " << xvar->GetName() << " expanded to nearest bin boundaries: [" 
			  << xlo << "," << xhi << "] --> [" << xloAdj << "," << xhiAdj << "]" << endl ;
    }


  } else {

    RooBinning xbins(fRHist->GetXaxis()->GetXmin(),fRHist->GetXaxis()->GetXmax()) ;
    xbins.addUniform(fRHist->GetNbinsX(),fRHist->GetXaxis()->GetXmin(),fRHist->GetXaxis()->GetXmax()) ;

    Double_t tolerance = 1e-6*xbins.averageBinWidth() ;

    // Adjust xlo/xhi to nearest boundary
    Double_t xloAdj = xbins.binLow(xbins.binNumber(xlo+tolerance)) ;
    Double_t xhiAdj = xbins.binHigh(xbins.binNumber(xhi-tolerance)) ;
    xbins.setRange(xloAdj,xhiAdj) ;
    xvar->setRange(xloAdj,xhiAdj) ;
    //xvar->setRange(xlo,xhi) ;
 
  }
  return;
}

///
/// Find the global minimum in a more thorough way.
/// First fit with external start parameters, then
/// for each parameter that starts with "d" or "r" (typically angles and ratios):
///   - at upper scan range, rest at start parameters
///   - at lower scan range, rest at start parameters
/// This amounts to a maximum of 1+2^n fits, where n is the number
/// of parameters to be varied.
///
/// \param w Workspace holding the pdf.
/// \param name Name of the pdf without leading "pdf_".
/// \param forceVariables Apply the force method for these variables only. Format
/// "var1,var2,var3," (list must end with comma). Default is to apply for all angles,
/// all ratios except rD_k3pi and rD_kpi, and the k3pi coherence factor.
///
RooFitResult* Utils::fitToMinForce(RooWorkspace *w, TString name, TString forceVariables)
{
	bool debug = true;

	TString parsName = "par_"+name;
	TString obsName  = "obs_"+name;
	TString pdfName  = "pdf_"+name;
	RooFitResult *r = 0;
	int printlevel = -1;
	RooMsgService::instance().setGlobalKillBelow(ERROR);

	// save start parameters
	if ( !w->set(parsName) ){
		cout << "MethodProbScan::scan2d() : ERROR : parsName not found: " << parsName << endl;
		exit(1);
	}
	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
//.........这里部分代码省略.........

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

      double confLevel = optPL.confLevel;
      double nScanPoints = optPL.nScanPoints;
      bool plotAsTF1 = optPL.plotAsTF1;
      double poiXMin = optPL.poiMinPlot;
      double poiXMax = optPL.poiMaxPlot;
      bool doHypoTest = optPL.doHypoTest;
      double nullParamValue = optPL.nullValue;

   // -------------------------------------------------------
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
      const char* filename = "";
      if (!strcmp(infile,"")) {
         filename = "results/example_combined_GaussExample_model.root";
         bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
         // if file does not exists generate with histfactory
         if (!fileExist) {
#ifdef _WIN32
            cout << "HistFactory file cannot be generated on Windows - exit" << endl;
            return;
#endif
            // Normally this would be run on the command line
            cout <<"will run standard hist2workspace example"<<endl;
            gROOT->ProcessLine(".! prepareHistFactory .");
            gROOT->ProcessLine(".! hist2workspace config/example.xml");
            cout <<"\n\n---------------------"<<endl;
            cout <<"Done creating example input"<<endl;
            cout <<"---------------------\n\n"<<endl;
         }

      }
      else
         filename = infile;

      // Try to open the file
      TFile *file = TFile::Open(filename);

      // if input file was specified byt not found, quit
      if(!file ){
         cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
         return;
      }

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

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

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

   // get the modelConfig out of the file
   RooAbsData* data = w->data(dataName);

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

   // ---------------------------------------------
   // create and use the ProfileLikelihoodCalculator
   // to find and plot the 95% confidence interval
   // on the parameter of interest as specified
   // in the model config
   ProfileLikelihoodCalculator pl(*data,*mc);
   pl.SetConfidenceLevel(confLevel); // 95% interval
   LikelihoodInterval* interval = pl.GetInterval();

   // print out the interval on the first Parameter of Interest
   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   cout << "\n>>>> RESULT : " << confLevel*100 << "% interval on " <<firstPOI->GetName()<<" is : ["<<
      interval->LowerLimit(*firstPOI) << ", "<<
      interval->UpperLimit(*firstPOI) <<"]\n "<<endl;

   // make a plot

   cout << "making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the TF1 drawing option)\n";
   LikelihoodIntervalPlot plot(interval);
   plot.SetNPoints(nScanPoints);  // do not use too many points, it could become very slow for some models
   if (poiXMin < poiXMax) plot.SetRange(poiXMin, poiXMax);
   TString opt;
   if (plotAsTF1) opt += TString("tf1");
   plot.Draw(opt);  // use option TF1 if too slow (plot.Draw("tf1")


   // if requested perform also an hypothesis test for the significance
   if (doHypoTest) {
//.........这里部分代码省略.........

void rf403_weightedevts()
{
  // C r e a t e   o b s e r v a b l e   a n d   u n w e i g h t e d   d a t a s e t 
  // -------------------------------------------------------------------------------

  // Declare observable
  RooRealVar x("x","x",-10,10) ;
  x.setBins(40) ;

  // Construction a uniform pdf
  RooPolynomial p0("px","px",x) ;

  // Sample 1000 events from pdf
  RooDataSet* data = p0.generate(x,1000) ;

 

  // C a l c u l a t e   w e i g h t   a n d   m a k e   d a t a s e t   w e i g h t e d 
  // -----------------------------------------------------------------------------------

  // Construct formula to calculate (fake) weight for events
  RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;

  // Add column with variable w to previously generated dataset
  RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;

  // Dataset d is now a dataset with two observable (x,w) with 1000 entries
  data->Print() ;

  // Instruct dataset wdata in interpret w as event weight rather than as observable
  RooDataSet wdata(data->GetName(),data->GetTitle(),data,*data->get(),0,w->GetName()) ;

  // Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430K
  wdata.Print() ;



  // U n b i n n e d   M L   f i t   t o   w e i g h t e d   d a t a 
  // ---------------------------------------------------------------

  // Construction quadratic polynomial pdf for fitting
  RooRealVar a0("a0","a0",1) ;
  RooRealVar a1("a1","a1",0,-1,1) ;
  RooRealVar a2("a2","a2",1,0,10) ;
  RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;

  // Fit quadratic polynomial to weighted data

  // NOTE: A plain Maximum likelihood fit to weighted data does in general 
  //       NOT result in correct error estimates, unless individual
  //       event weights represent Poisson statistics themselves.
  //       
  // Fit with 'wrong' errors
  RooFitResult* r_ml_wgt = p2.fitTo(wdata,Save()) ;
  
  // A first order correction to estimated parameter errors in an 
  // (unbinned) ML fit can be obtained by calculating the
  // covariance matrix as
  //
  //    V' = V C-1 V
  //
  // where V is the covariance matrix calculated from a fit
  // to -logL = - sum [ w_i log f(x_i) ] and C is the covariance
  // matrix calculated from -logL' = -sum [ w_i^2 log f(x_i) ] 
  // (i.e. the weights are applied squared)
  //
  // A fit in this mode can be performed as follows:

  RooFitResult* r_ml_wgt_corr = p2.fitTo(wdata,Save(),SumW2Error(kTRUE)) ;



  // P l o t   w e i g h e d   d a t a   a n d   f i t   r e s u l t 
  // ---------------------------------------------------------------

  // Construct plot frame
  RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;

  // Plot data using sum-of-weights-squared error rather than Poisson errors
  wdata.plotOn(frame,DataError(RooAbsData::SumW2)) ;

  // Overlay result of 2nd order polynomial fit to weighted data
  p2.plotOn(frame) ;



  // M L  F i t   o f   p d f   t o   e q u i v a l e n t  u n w e i g h t e d   d a t a s e t
  // -----------------------------------------------------------------------------------------
  
  // Construct a pdf with the same shape as p0 after weighting
  RooGenericPdf genPdf("genPdf","x*x+10",x) ;

  // Sample a dataset with the same number of events as data
  RooDataSet* data2 = genPdf.generate(x,1000) ;

  // Sample a dataset with the same number of weights as data
  RooDataSet* data3 = genPdf.generate(x,43000) ;

  // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
  RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  int i=0;//color index
  TLegend *leg = new TLegend(0.2, 0.02, .4, .42);  
  leg->SetTextSize(0.06);
  leg->AddEntry(framex->findObject("curvetot"),"Total PDF","l");
  while(vartPC){//loop over compotents of totalPdf
    TString vartPCtitle = vartPC->GetTitle();
    TIterator* itercompPars;//forward declare so it persists outside the if statement
    RooRealVar* varcompPars;
    if(!(vartPCtitle.Contains(":")||vartPCtitle.Contains("@"))){//only for non-sub-shapes
      while(i==0||i==10||i==4||i==1||i==5||(i>=10&&i<=27))i++;//avoid white and blue and black and yellow and horribleness
      RooArgSet* compPars = vartPC->getParameters(data);//set of the parameters of the component the loop is on
      itercompPars = compPars->createIterator();
      varcompPars = (RooRealVar*) itercompPars->Next();
    
      while(varcompPars){//write and print mean, sig, etc. of sub-shapes
        TString vartitle = varcompPars->GetTitle();
        double varval = varcompPars->getVal();
        TString varvalstring = Form("%f",varval);
        double hi = varcompPars->getErrorHi();
        
        TString varerrorstring = "[exact]";
        if(hi!=-1){
          double lo = varcompPars->getErrorLo();
          double varerror = TMath::Max(fabs(lo),hi);
          varerrorstring = Form("%E",varerror);
        }
        
        outputtext = vartitle+" = "+varvalstring+" +/- "+varerrorstring;
        textfile<<outputtext<<endl;
        cout<<outputtext<<endl;
        
        varcompPars = (RooRealVar*) itercompPars->Next(); 
      }
      totalPdf.plotOn(framex,Name(vartPC->GetName()),LineStyle(kDashed),LineColor(i),Components(vartPC->GetName()));
      leg->AddEntry(framex->findObject(vartPC->GetName()),vartPCtitle,"l");
    
      i++;
    }
    vartPC = (RooAddPdf*) itertPC->Next();
    itercompPars->Reset();//make sure it's ready for the next vartPC
  }
  
  // Calculate chi2/ndf
  RooArgSet *floatpar = totalPdf.getParameters(data);
  int floatpars = (floatpar->selectByAttrib("Constant",kFALSE))->getSize();
  Double_t chi2 = framex->chiSquare("curvetot","Hist",floatpars);
  TString chi2string = Form("%f",chi2);
  //create text box to list important parameters on the plot
  // TPaveText* txt = new TPaveText(0.1,0.5,0.7,0.9,"NBNDC");
  // txt->SetTextSize(0.06);
  // txt->SetTextColor(kBlack);
  // txt->SetBorderSize(0);
  // txt->SetFillColor(0);
  // txt->SetFillStyle(0);
  outputtext = "#chi^{2}/N_{DoF} = "+chi2string;
  cout<<outputtext<<endl;
  textfile<<outputtext<<endl;
  // txt->AddText(outputtext);
  
  // Print stuff
  TIterator* iteryields =  yields.createIterator();
  RooRealVar* varyields = (RooRealVar*) iteryields->Next();//only inherits things from TObject unless class specified
  vector<double> Y, E;//holds yields and associated errors
  vector<TString> YS, ES;//holds strings of the corresponding yields
  int j=0;//count vector position
  int jS=0, jL=0;//these hold the position of the S and L results;initialized in case there is no nsigS or nsigL

void makeAndImportDataSets(TFile *fin, RooWorkspace *wspace, RooRealVar &x){

   // Get TTrees from input files and fill RooDataSets!
   //const char *name = x.GetName();
   
   RooRealVar weight("weight","weight",0,100);
   weight.removeRange();

   RooArgSet treevariables(x,weight);
   // List all branches in the tree and add them as variables 
   TObjArray *branches = (TObjArray*) ((TTree*)fin->Get("data_obs"))->GetListOfBranches();
   TIter next(branches);
   TBranch *br;
   while ( (br = (TBranch*)next()) ){
	const char *brname = br->GetName();
	std::cout << brname << std::endl;
	if ( *brname==*weight.GetName() ) continue; 
	if ( *brname==*x.GetName() ) continue;
	RooRealVar *vartmp = new RooRealVar(brname,brname,0,1); vartmp->removeRange();
	treevariables.add(*vartmp);
   }

 //  RooRealVar cvar(cutvar.c_str(),cutvar.c_str(),0,100);
 //  cvar.removeRange();
    
   std::cout << "Building RooDataSets " << std::endl; 
   treevariables.Print("V");
   //assert(0);
   // Zvv Signal region MC sample 
   RooDataSet zvv("DY","DY",treevariables,RooFit::Import(*((TTree*)fin->Get("DY"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));

   // Zvv Control sample 
   // Data ... (double mu)  Zvv_control
   RooDataSet zvvcontrol("Zvv_control","Zvv_control",treevariables,RooFit::Import(*((TTree*)fin->Get("Zvv_control"))),RooFit::Cut(cutstring.c_str()));

   // Backgrounds ...	    Zvv_control_bkg_mc
   RooDataSet zvvcontrolbkgmc("Zvv_control_bkg_mc","Zvv_control_bkg_mc",treevariables,RooFit::Import(*((TTree*)fin->Get("T_control_bkg_mc"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet zvvcontrolbkgmc_1("1","1",treevariables,RooFit::Import(*(TTree*)fin->Get("TT_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet zvvcontrolbkgmc_2("2","2",treevariables,RooFit::Import(*(TTree*)fin->Get("WW_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet zvvcontrolbkgmc_3("3","3",treevariables,RooFit::Import(*(TTree*)fin->Get("WZ_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet zvvcontrolbkgmc_4("4","4",treevariables,RooFit::Import(*(TTree*)fin->Get("ZZ_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   zvvcontrolbkgmc.append(zvvcontrolbkgmc_1);
   zvvcontrolbkgmc.append(zvvcontrolbkgmc_2);
   zvvcontrolbkgmc.append(zvvcontrolbkgmc_3);
   zvvcontrolbkgmc.append(zvvcontrolbkgmc_4);

   // MC ...		    Zvv_control_mc
   RooDataSet zvvcontrolmc("Zvv_control_mc","Zvv_control_mc",treevariables,RooFit::Import(*((TTree*)fin->Get("Zvv_control_mc"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   std::cout << " Weighs of datasets "  << std::endl;
   std::cout << " Zvv_control " << zvvcontrol.sumEntries() << std::endl;
   std::cout << " Zvv_control_mc " << zvvcontrolmc.sumEntries() << std::endl;

   // Wlv Control sample (single mu)
   // Data ... (double mu)
   // Backgrounds ...
   //.append()
   // MC ...

   // Store these to the output workspace  
   wspace->import(zvvcontrol);
   wspace->import(zvvcontrolmc);
   wspace->import(zvvcontrolbkgmc);
   wspace->import(zvv);


   // Wlv Signal region MC sample 
   RooDataSet wlv("W","W",treevariables,RooFit::Import(*((TTree*)fin->Get("W"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlv_ht("WHT","WHT",treevariables,RooFit::Import(*((TTree*)fin->Get("WHT"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   wlv.append(wlv_ht);

   // Wlv Control sample 
   // Data ... (single mu)  Wlv_control
   RooDataSet wlvcontrol("Wlv_control","Wlv_control",treevariables,RooFit::Import(*((TTree*)fin->Get("Wlv_control"))),RooFit::Cut(cutstring.c_str()));

   // Backgrounds ...	    Wlv_control_bkg_mc
   RooDataSet wlvcontrolbkgmc("Wlv_control_bkg_mc","Wlv_control_bkg_mc",treevariables,RooFit::Import(*((TTree*)fin->Get("T_sl_control_bkg_mc"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlvcontrolbkgmc_1("1","1",treevariables,RooFit::Import(*(TTree*)fin->Get("TT_sl_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlvcontrolbkgmc_2("2","2",treevariables,RooFit::Import(*(TTree*)fin->Get("DY_sl_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlvcontrolbkgmc_3("3","3",treevariables,RooFit::Import(*(TTree*)fin->Get("WZ_sl_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlvcontrolbkgmc_4("4","4",treevariables,RooFit::Import(*(TTree*)fin->Get("ZZ_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlvcontrolbkgmc_5("5","5",treevariables,RooFit::Import(*(TTree*)fin->Get("WW_control_bkg_mc")),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   wlvcontrolbkgmc.append(wlvcontrolbkgmc_1);
   wlvcontrolbkgmc.append(wlvcontrolbkgmc_2);
   wlvcontrolbkgmc.append(wlvcontrolbkgmc_3);
   wlvcontrolbkgmc.append(wlvcontrolbkgmc_4);
   wlvcontrolbkgmc.append(wlvcontrolbkgmc_5);

   // MC ...		    Zvv_control_mc
   RooDataSet wlvcontrolmc("Wlv_control_mc","Wlv_control_mc",treevariables,RooFit::Import(*((TTree*)fin->Get("Wlv_control_mc_1"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   RooDataSet wlvcontrolmc_2("Wlv_control_mc_2","Wlv_control_mc_2",treevariables,RooFit::Import(*((TTree*)fin->Get("Wlv_control_mc_2"))),RooFit::Cut(cutstring.c_str()),RooFit::WeightVar("weight"));
   wlvcontrolmc.append(wlvcontrolmc_2);
   std::cout << " Weighs of datasets "  << std::endl;
   std::cout << " Wlv_control "    << wlvcontrol.sumEntries() << std::endl;
   std::cout << " Wlv_control_mc " << wlvcontrolmc.sumEntries() << std::endl;

   // Wlv Control sample (single mu)
   // Data ... (double mu)
   // Backgrounds ...
   //.append()
   // MC ...
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    cout <<"will run standard hist2workspace example"<<endl;
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;

    // now try to access the file again
    file = TFile::Open(filename);
  }

  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

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

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

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

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

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

  /////////////////////////////////////////////
  // create and use the BayesianCalculator
  // to find and plot the 95% credible interval
  // on the parameter of interest as specified
  // in the model config
  
  // before we do that, we must specify our prior
  // it belongs in the model config, but it may not have
  // been specified
  RooUniform prior("prior","",*mc->GetParametersOfInterest());
  w->import(prior);
  mc->SetPriorPdf(*w->pdf("prior"));

  // do without systematics
  //mc->SetNuisanceParameters(RooArgSet() );

  
  BayesianCalculator bayesianCalc(*data,*mc);
  bayesianCalc.SetConfidenceLevel(0.95); // 95% interval

  // default of the calculator is central interval.  here use shortest , central or upper limit depending on input
  // doing a shortest interval might require a longer time since it requires a scan of the posterior function
  if (intervalType == 0)  bayesianCalc.SetShortestInterval(); // for shortest interval
  if (intervalType == 1)  bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval
  if (intervalType == 2)  bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit

  if (!integrationType.IsNull() ) { 
     bayesianCalc.SetIntegrationType(integrationType); // set integrationType
     bayesianCalc.SetNumIters(nToys); // set number of ietrations (i.e. number of toys for MC integrations)
  }

  // compute interval by scanning the posterior function
  if (scanPosterior)   
     bayesianCalc.SetScanOfPosterior(nScanPoints);


  SimpleInterval* interval = bayesianCalc.GetInterval();

  // print out the iterval on the first Parameter of Interest
  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  cout << "\n95% interval on " <<firstPOI->GetName()<<" is : ["<<
    interval->LowerLimit() << ", "<<
    interval->UpperLimit() <<"] "<<endl;


  // make a plot 
  // since plotting may take a long time (it requires evaluating 
  // the posterior in many points) this command will speed up 
  // by reducing the number of points to plot - do 50

  cout << "\nDrawing plot of posterior function....." << endl;

  bayesianCalc.SetScanOfPosterior(nScanPoints);

  RooPlot * plot = bayesianCalc.GetPosteriorPlot();
  plot->Draw();  

}

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


  // the number of toy MC used to generate the distribution
  int nToyMC = 1000;
  // The parameter below is needed for asymptotic distribution to be chi-square,
  // but set to false if your model is not numerically stable if mu<0
  bool allowNegativeMu=true;


  /////////////////////////////////////////////////////////////
  // First part is just to access a user-defined file
  // or create the standard example file if it doesn't exist
  ////////////////////////////////////////////////////////////
   const char* filename = "";
   if (!strcmp(infile,"")) {
      filename = "results/example_combined_GaussExample_model.root";
      bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
      // if file does not exists generate with histfactory
      if (!fileExist) {
#ifdef _WIN32
         cout << "HistFactory file cannot be generated on Windows - exit" << endl;
         return;
#endif
         // Normally this would be run on the command line
         cout <<"will run standard hist2workspace example"<<endl;
         gROOT->ProcessLine(".! prepareHistFactory .");
         gROOT->ProcessLine(".! hist2workspace config/example.xml");
         cout <<"\n\n---------------------"<<endl;
         cout <<"Done creating example input"<<endl;
         cout <<"---------------------\n\n"<<endl;
      }

   }
   else
      filename = infile;

   // Try to open the file
   TFile *file = TFile::Open(filename);

   // if input file was specified byt not found, quit
   if(!file ){
      cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
      return;
   }


  /////////////////////////////////////////////////////////////
  // Now get the data and workspace
  ////////////////////////////////////////////////////////////

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

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

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

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

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

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

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

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