C++ RooArgList类代码示例

RooAddPdf *buildSumOfGaussians(string name, RooRealVar *mass, RooRealVar *MH, int nGaussians){
    
  RooArgList *gaussians = new RooArgList();
  RooArgList *coeffs = new RooArgList();
  for (int g=0; g<nGaussians; g++){
    RooRealVar *dm = new RooRealVar(Form("dm_g%d",g),Form("dm_g%d",g),0.1,-5.*(1.+0.5*g),5.*(1.+0.5*g));
    RooAbsReal *mean = new RooFormulaVar(Form("mean_g%d",g),Form("mean_g%d",g),"@[email protected]",RooArgList(*MH,*dm));
    RooRealVar *sigma = new RooRealVar(Form("sigma_g%d",g),Form("sigma_g%d",g),2.,0.7,5.*(1.+0.5*g));
    RooGaussian *gaus = new RooGaussian(Form("gaus_g%d",g),Form("gaus_g%d",g),*mass,*mean,*sigma);
    //tempFitParams.insert(pair<string,RooRealVar*>(string(dm->GetName()),dm));
    //tempFitParams.insert(pair<string,RooRealVar*>(string(sigma->GetName()),sigma));
    //tempFitUtils.insert(pair<string,RooAbsReal*>(string(mean->GetName()),mean));
    //tempGaussians.insert(pair<string,RooGaussian*>(string(gaus->GetName()),gaus));
    gaussians->add(*gaus);
    if (g<nGaussians-1) {
      RooRealVar *frac = new RooRealVar(Form("frac_g%d",g),Form("frac_g%d",g),0.1,0.01,0.99);
      //tempFitParams.insert(pair<string,RooRealVar*>(string(frac->GetName()),frac));
      coeffs->add(*frac);
    }
    if (g==nGaussians-1 && forceFracUnity_){
      string formula="1.";
      for (int i=0; i<nGaussians-1; i++) formula += Form("[email protected]%d",i);
      RooAbsReal *recFrac = new RooFormulaVar(Form("frac_g%d",g),Form("frac_g%d",g),formula.c_str(),*coeffs);
      //tempFitUtils.insert(pair<string,RooAbsReal*>(string(recFrac->GetName()),recFrac));
      coeffs->add(*recFrac);
    }
  }
  RooAddPdf *tempSumOfGaussians = new RooAddPdf(name.c_str(),name.c_str(),*gaussians,*coeffs,recursive_);
  return tempSumOfGaussians;
}

///
/// Fills vector with floating pars names
///
void Utils::getParameters(const RooFitResult &result, std::vector<TString> &names){
	RooArgList pars           = result.floatParsFinal();
	TIterator * it            = pars.createIterator();
	while(RooRealVar* p = (RooRealVar*) it->Next()){
		names.push_back(TString(p->GetName()));
	}
};

void Fitter::MakePlots(const string& name, const string& title, RooPlot*& pPlot1, RooPlot*& pPlot2, RooAbsPdf* pPdf, RooDataSet* pData)
{
	bool plotbkg = false;
	RooPlot* pFrame = pPlot1? pPlot1 : pFitVar->frame(RooFit::Bins(50), RooFit::Name(name.c_str()), RooFit::Title(title.c_str()));
	RooPlot* pFramePull = pPlot2? pPlot2 : pFitVar->frame(RooFit::Bins(50), RooFit::Name((name+"_Pull").c_str()), RooFit::Title((title+" Pull").c_str()));
	RooHist* pHisPull = NULL;
	double pullMax = 0.;
	
	if(!pPdf && !pData) plotbkg = true;
	if(!pPdf) pPdf = pTotPdf;
	if(!pData) pData = this->pData;
	
	pData->plotOn(pFrame, RooFit::MarkerSize(0.9));
	
	if(plotbkg)
	{
		RooArgList bkgComponent;
		for(auto& pdf : bkgPdfs) bkgComponent.add(pdf.GetPdf());
		pTotPdf->plotOn(pFrame, RooFit::Components(bkgComponent), RooFit::LineColor(kGreen) );
	}
	
	pPdf->plotOn(pFrame, RooFit::LineColor(kRed));	
	pPdf->paramOn(pFrame, RooFit::Format("NEU", RooFit::AutoPrecision(1)), RooFit::Layout(0.65, 0.99, 0.93));
	pFrame->getAttText()->SetTextSize(0.75*pFrame->GetYaxis()->GetLabelSize());
	
	pHisPull = pFrame->pullHist();
	pFramePull->addPlotable(pHisPull, "P");
	pullMax = max(fabs(pFramePull->GetMinimum()), pFramePull->GetMaximum());
	pFramePull->GetYaxis()->SetRangeUser(-pullMax,pullMax);
	
	pPlot1 = pFrame;
	pPlot2 = pFramePull;
}

testWrapper(){
	gSystem->Load("libTMVAWrapper.so");
	RooWorkspace w("w");
	w.factory("x[-5,5]");
	w.factory("alpha[-5,5]");
	x = w.var("x");
	alpha = w.var("alpha");
	//	vars = w.allVars();
	RooArgList vars;
	vars.add(*x);
	vars.add(*alpha);
	vars.Print();

	TMVAWrapper wrap("wrap","wrap", vars, "TMVARegression_alphavary.root_MLP.weights.xml");
	wrap.Print("v");
	cout << "wrap = " << wrap.getVal() <<endl;
	x->setVal(1);
	alpha->setVal(2.);
	cout << "wrap = " << wrap.getVal() <<endl;

//	w.import(wrap);
//	w.Print();
//	return;
	//	cout << "\n\nvars are: " << wrap.getListOfVars() << endl;
	frame = x->frame();
	wrap.plotOn(frame);
	frame->Draw();
}

void THSRooFit::SetDataWeight(){
  if(!fInWeights) return;
  if(fInWeights->GetSpeciesID(fWeightName)<0) return;
  RooArgList setWeights;
  RooRealVar wvar(fWeightName,fWeightName,0);
  setWeights.add(wvar);
  setWeights.Print();
  RooDataSet wData("weights","weights",setWeights) ;

  //Loop over the weights find the entry for the event ID
  for(Int_t ev=0;ev<fData->numEntries();ev++){
    const RooArgSet* vars=fData->get(ev);
    fInWeights->GetEntryBinarySearch((Long64_t)vars->getRealValue(fIDBranchName));
    // wvar.setVal(fInWeights->GetWeight(fWeightName));
    wvar.setVal(fInWeights->GetWeight(fWeightName));
    wData.add(setWeights);
  }
  wData.Print();
  ((RooDataSet*)fData)->merge(&wData);
  fData->Print();
  RooDataSet* tmpData=0;
  tmpData=new RooDataSet(TString("Weight_")+fWeightName+fData->GetName(),fData->GetTitle(),(RooDataSet*)fData,*fData->get(),0,fWeightName);
  fData->SetName("Empty");fData->reset();fData=0;
  fData=tmpData;
  fData->Print();
  if(ftoWS) fWS->import(*fData); //import if told to
  else if(!(fDataBins))fWS->import(*fData); //or no databins to confuse it 
}

map<string,RooAddPdf*> getMITPdfs(RooWorkspace *work, int ncats, bool is2011){

  map<string,RooAddPdf*> result;

  RooArgList pdfs;
  RooArgList norms;

  string sqrts;
  string catstring;
  if (is2011){
    sqrts="7TeV";
    catstring="cat";
  }
  else {
    sqrts="8TeV";
    catstring="mvacat";
  }

  for (int cat=0; cat<ncats; cat++){
    RooAddPdf *ggh = (RooAddPdf*)work->pdf(Form("hggpdfsmrel_%s%d_%s_ggh",catstring.c_str(),cat,sqrts.c_str())); 
    RooAddPdf *vbf = (RooAddPdf*)work->pdf(Form("hggpdfsmrel_%s%d_%s_vbf",catstring.c_str(),cat,sqrts.c_str())); 
    RooAddPdf *wzh = (RooAddPdf*)work->pdf(Form("hggpdfsmrel_%s%d_%s_wzh",catstring.c_str(),cat,sqrts.c_str())); 
    RooAddPdf *tth = (RooAddPdf*)work->pdf(Form("hggpdfsmrel_%s%d_%s_tth",catstring.c_str(),cat,sqrts.c_str()));
    RooFormulaVar *ggh_norm = (RooFormulaVar*)work->function(Form("hggpdfsmrel_%s%d_%s_ggh_norm",catstring.c_str(),cat,sqrts.c_str()));
    RooFormulaVar *vbf_norm = (RooFormulaVar*)work->function(Form("hggpdfsmrel_%s%d_%s_vbf_norm",catstring.c_str(),cat,sqrts.c_str()));
    RooFormulaVar *wzh_norm = (RooFormulaVar*)work->function(Form("hggpdfsmrel_%s%d_%s_wzh_norm",catstring.c_str(),cat,sqrts.c_str()));
    RooFormulaVar *tth_norm = (RooFormulaVar*)work->function(Form("hggpdfsmrel_%s%d_%s_tth_norm",catstring.c_str(),cat,sqrts.c_str()));
    pdfs.add(*ggh); pdfs.add(*vbf); pdfs.add(*wzh); pdfs.add(*tth);
    norms.add(*ggh_norm); norms.add(*vbf_norm); norms.add(*wzh_norm); norms.add(*tth_norm);
    result.insert(pair<string,RooAddPdf*>(Form("cat%d",cat),new RooAddPdf(Form("sig_cat%d",cat),"sig",RooArgList(*ggh,*vbf,*wzh,*tth),RooArgList(*ggh_norm,*vbf_norm,*wzh_norm,*tth_norm))));
  }
  result.insert(pair<string,RooAddPdf*>("all",new RooAddPdf("sig_all","sig",pdfs,norms)));
  return result;

}

RooArgList Fitter::GetYieldsList()
{
	RooArgList yldList;
	for(auto& pVec_pdfs : {&sigPdfs, &bkgPdfs})
	for(auto& pdf : *pVec_pdfs)
		yldList.add(pdf.GetYld());
		
	return yldList;
}

void makeSystPlot(RooWorkspace *tempWS, RooWorkspace *WS,  TString systT, int Bin, int toMassNo, int fromMassNo)  
{

  std::cout<< "SystName--> " << systT.Data() << std::endl;
  RooArgList  * hobs = new RooArgList("hobs");

  TString xname;

  if(IFILE.Contains ("Zmm") | IFILE.Contains ("Zee"))
  { 
    xname = "CMS_vhbb_BDT_Zll_8TeV";
  }
  if(IFILE.Contains("We") | IFILE.Contains("Wm")) 
  {
    xname = "CMS_vhbb_BDT_Wln_8TeV";
  }
  if(IFILE.Contains ("Zn") && (Bin == 0))
  {
    xname = "CMS_vhbb_BDT_ZnunuLowPt_8TeV";
  }

  if(IFILE.Contains ("Zn") && (Bin == 1)) 
  {
    xname = "CMS_vhbb_BDT_ZnunuMedPt_8TeV";
  }

  if(IFILE.Contains ("Zn") && (Bin == 2)) 
  {
    xname = "CMS_vhbb_BDT_ZnunuHighPt_8TeV";
  }
  if(IFILE.Contains("Wt")) 
  {
    xname = "BDT";
  }

  std::cout<< "BDT Name--> " << xname.Data() << std::endl;
  hobs->add(*WS->var(xname.Data()));
  RooRealVar BDT(xname.Data(),xname.Data(),-1,1); 

  tempRooDataHistNom = (RooDataHist*)  tempWS->data(systT.Data());
  tempHistNom = tempRooDataHistNom->createHistogram(systT.Data(),BDT,RooFit::Binning(bins));
  
  if(systT.Contains("WH"))
  {
     tempHistNom->Scale(xSec8WH[toMassNo]/xSec8WH[fromMassNo]);
  }  
 
  if(systT.Contains("ZH"))
  {
     tempHistNom->Scale(xSec8ZH[toMassNo]/xSec8ZH[fromMassNo]);
  }  
 
  RooDataHist *DHnom;
  DHnom = new RooDataHist(systT.Data(),"",*hobs,tempHistNom);
  WS->import(*(new RooHistPdf(systT.Data(),"",*hobs,*DHnom))); 
}

void Utils::setParameters(RooWorkspace* w, TString parname, RooSlimFitResult* r, bool constAndFloat)
{
	if ( constAndFloat ){
		RooArgList list = r->floatParsFinal();
		list.add(r->constPars());
		setParameters(w, parname, &list);
		return;
	}
	setParameters(w, parname, &(r->floatParsFinal()));
}

//#include "/uscms_data/d3/cvernier/DiH_13TeV/CMSSW_7_1_5/src/HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
//#include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
void makeRooMultiPdfWorkspace(){

   // Load the combine Library 
   gSystem->Load("libHiggsAnalysisCombinedLimit.so");

   // Open the dummy H->gg workspace 
   TFile *f_hgg = TFile::Open("w_background_Bern.root");
   RooWorkspace *w_hgg = (RooWorkspace*)f_hgg->Get("HbbHbb");
   // The observable (CMS_hgg_mass in the workspace)
   RooRealVar *mass =  w_hgg->var("x");

   // Get three of the functions inside, exponential, linear polynomial, power law
   RooAbsPdf *pdf_exp = w_hgg->pdf("bg_exp");
   RooAbsPdf *pdf_pol = w_hgg->pdf("bg");


   // Fit the functions to the data to set the "prefit" state (note this can and should be redone with combine when doing 
   // bias studies as one typically throws toys from the "best-fit"
   RooAbsData *data = w_hgg->data("data_obs");
   pdf_exp->fitTo(*data);  // index 0
   pdf_pol->fitTo(*data);   // index 2

   // Make a plot (data is a toy dataset)
   RooPlot *plot = mass->frame();   data->plotOn(plot);
   pdf_exp->plotOn(plot,RooFit::LineColor(kBlue));
   pdf_pol->plotOn(plot,RooFit::LineColor(kRed));
   plot->SetTitle("PDF fits to toy data");
   plot->Draw();

   // Make a RooCategory object. This will control which of the pdfs is "active"
   RooCategory cat("pdf_index","Index of Pdf which is active");

   // Make a RooMultiPdf object. The order of the pdfs will be the order of their index, ie for below 
   // 0 == exponential
   // 1 == linear function
   // 2 == powerlaw
   RooArgList mypdfs;
   mypdfs.add(*pdf_exp);
   mypdfs.add(*pdf_pol);
   
   RooMultiPdf multipdf("roomultipdf","All Pdfs",cat,mypdfs);
   
   // As usual make an extended term for the background with _norm for freely floating yield
   RooRealVar norm("roomultipdf_norm","Number of background events",0,10000);
   
   // Save to a new workspace
   TFile *fout = new TFile("background_pdfs.root","RECREATE");
   RooWorkspace wout("backgrounds","backgrounds");
   wout.import(cat);
   wout.import(norm);
   wout.import(multipdf);
   wout.Print();
   wout.Write();

}

   void GenericModel::AddModelData(RooArgList& params, RooDataSet* data)
   {
      // Add a set of data calculated with the model for the given set of parameter values.

      if (!data) return;

      fDataSets.Add(data);
      RooArgList* ral = new RooArgList;
      ral->addClone(params);
      fDataSetParams.Add(ral);//keep a copy of the parameter values
      ++fNDataSets;
   }

JetTagBin::JetTagBin(TFile & file, string jtBinName, vector<Sample> theSamples):
  jtBinName_(jtBinName)
{
  jtBinVar_= new RooRealVar(jtBinName.c_str(), jtBinName.c_str(), 50, 600);
  int nbins_def = 0;
  for(unsigned int i=0;i<theSamples.size();++i) {
    string name = theSamples[i].name()+jtBinName;
    string histname = "";
    if ( theSamples[i].sampName()=="" ) histname = name;
    else histname = theSamples[i].sampName()+jtBinName;
    if ( verbose ) cout << "Opening histogram " << histname << endl;
    TH1F * tempHisto = (TH1F*)file.Get(histname.c_str());
    TH1F * tempHisto2;
    if(tempHisto==0) {
      cout << "Template " << name << " does not exist" << endl;
      assert(0); // maybe we don't assert but rather break/continue so we can easily work with different jtbin names
    }
    else {
      tempHisto2 = (TH1F*) tempHisto->Clone();
      //      delete tempHisto;
    }

    // We set the binning, assuming all samples have the same binning
    int nbins = tempHisto2->GetNbinsX();
    if(i==0) {
      nbins_def = nbins;
      int xlo   = tempHisto2->GetBinLowEdge(1);
      int xhi   = tempHisto2->GetBinLowEdge(nbins + 1);
      RooUniformBinning jtBinning(xlo,xhi,nbins_def);
      jtBinVar_->setBinning(jtBinning);
    }
    else {
      if (nbins!=nbins_def) {
        cout << "something here " << endl;
        assert(0);
      }
    }
    
    Template * temp = new Template(*tempHisto2, name, *jtBinVar_, theSamples[i]);
    jtTemplates_.push_back( temp );
    delete tempHisto2;

  }
  if ( verbose ) cout << "Making combined pdf" << endl;
  RooArgList TemplateList;
  for(unsigned int i=0;i<jtTemplates_.size();++i) {
    TemplateList.add(*(jtTemplates_[i]->epdf()));
  }
  jtPdf_ = new RooAddPdf((jtBinName_+"_pdf").c_str(),(jtBinName_+"_pdf").c_str(),TemplateList);
  jtData_ = 0;
  jtNll_ = 0;
}

//--------------------------------------------------------------------------------------------------
void printCorrelations(ostream& os, RooFitResult *res)
{
  ios_base::fmtflags flags = os.flags();
  const RooArgList parlist = res->floatParsFinal();
  
  os << "  Correlation Matrix" << endl;
  os << " --------------------" << endl;
  for(Int_t i=0; i<parlist.getSize(); i++) {
    for(Int_t j=0; j<parlist.getSize(); j++) 
      os << "  " << setw(7) << setprecision(4) << fixed << res->correlationMatrix()(i,j);    
    os << endl;
  }
  os.flags(flags);
}

///
/// searches for unphysical values in a fit result, fills problematic var names in vectors
/// returns true if everything is within boundaries
///
bool Utils::checkBoundary(const RooSlimFitResult& r, std::vector<TString> lowProb, std::vector<TString> highProb){
	// make sure vectors are empty
	lowProb.clear(); highProb.clear();
	RooArgList floats = r.floatParsFinal();
	TIterator* floatIt = floats.createIterator();
	while(RooRealVar* var = (RooRealVar*)floatIt->Next()){
		//check lower bound
		if(var->getVal() < var->getMin("phys")){
			lowProb.push_back(var->GetName());
		}
		// check higher bound
		if( var->getVal() > var->getMax("phys")){
			highProb.push_back(var->GetName());
		}
	}
	return (lowProb.size() != 0 || highProb.size() != 0) ? false : true;
}

void ParameterCache::printFitResultToOutStream(ofstream &out, RooSlimFitResult *slimFitRes) {

	out << "### FCN: " << slimFitRes->minNll() << ", EDM: " << slimFitRes->edm() << endl;
	out << "### COV quality: " << slimFitRes->covQual() << ", status: " << slimFitRes->status()
		<< ", confirmed: " << (slimFitRes->isConfirmed()?"yes":"no") << endl;
	RooArgList argList = slimFitRes->floatParsFinal();
	argList.add(slimFitRes->constPars());
	argList.sort();
	TIterator *iter = argList.createIterator();
	while ( RooRealVar *arg=(RooRealVar*)iter->Next() ) {
		if ( TString(arg->GetName()).Contains("obs") ) continue;
		out << Form("%-25s",arg->GetName()) << " " << Form("%12.6f",arg->getVal())
			<< " " << Form("%12.6f",arg->getErrorLo())
			<< " " << Form("%12.6f",arg->getErrorHi()) << endl;
	}
	delete iter;
}