C++ RooRealVar类代码示例

void simplePrintResults() {
    vector<string> filenames;
    filenames.push_back("Output/Test2/result/FIT_DATA_Psi2SJpsi_PPPrompt_Bkg_SecondOrderChebychev_pt65300_rap016_cent0200_262620_263757.root");

    const char* parname = "N_Jpsi_PP";

    vector<string>::iterator it = filenames.begin();
    for (it; it<filenames.end(); it++) {
        TFile *f = new TFile(it->c_str());
        if (!f) {
            cout << "Error, " << *it << " not found" << endl;
            continue;
        }
        RooWorkspace *ws = (RooWorkspace*) f->Get("workspace");
        if (!ws) {
            cout << "Error, workspace not found in " << *it << endl;
            continue;
        }

        RooRealVar *var = ws->var(parname);
        if (!ws) {
            cout << "Error, variable " << parname << " not found in " << *it << endl;
            continue;
        }
        cout << *it << " " << var->getVal() << " +- " << var->getError() << endl;
    }
}

RooStats::ModelConfig * Tprime::SetBModel( void ) {
    //
    // Define model config and parameter snapshot to describe the b model.
    // Import to workspace.
    //

    std::string legend = "[Tprime::SetBModel]: ";

    // full signal+background model
    //RooStats::ModelConfig * pSBModel = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");

    // let's make the b model (bg-only) from the alt model (s+b) with xsec=0
    //RooStats::ModelConfig * pBModel =
    //  new RooStats::ModelConfig(*(RooStats::ModelConfig *)pWs->genobj("ModelConfig"));
    RooStats::ModelConfig * _sbModel = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");
    RooStats::ModelConfig * pBModel = _sbModel->Clone("BModel");

    //pBModel->SetName("BModel");

    pBModel->SetWorkspace(*pWs);

    //pBModel->SetParametersOfInterest(RooArgSet());
    pWs->import(*pBModel);

    // set POI to the b model value and take snapshot
    RooRealVar * pPoi = (RooRealVar *)pBModel->GetParametersOfInterest()->first();
    pPoi->setVal(0.0);
    pBModel->SetSnapshot(*pPoi);

    pBModel->Print();

    return pBModel;
}

void PEs(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	 RooRealVar &iScale,RooRealVar &iRes) { 
  double iM0 = iMean.getVal(); double iS0 = iSig.getVal();
  //iScale.setVal(iMeanScale); iRes.setVal(iSigScale);
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  TNtuple * lDN= new TNtuple( "xxx","xxx","ntot:m_r:m:merr:sig_r:sig:sigerr");
  for(int i0=0;i0<iN;i0++){
    if(i0 % 10 == 0) cout << "+++++++++++++++++++++++++++ running ======> " << i0 << endl;
    int lN    = lRand->Poisson(iNEvents);
    RooDataSet * lSignal  = iGen->generate(iVar,lN);
    iMean.setVal(iM0); iSig.setVal(iS0);
    iFit->fitTo(*lSignal,Strategy(1));//,Save(kTRUE),PrintLevel(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
    Float_t values[]={
      (Float_t) lN,
      (Float_t) 90.78/iScale.getVal(),
      (Float_t) iMean.getVal(),
      (Float_t) iMean.getError(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getError()
    };
    lDN->Fill(values);
  }
  TFile *lF = new TFile("XXX.root","RECREATE");
  lDN->Write();
  lF->Close();
}

TGraph *graphLH(std::string nuisname, double err ){

	w->loadSnapshot("bestfitall"); // SetTo BestFit values as start

	// Get The parameter we want 
	RooRealVar *nuis =(RooRealVar*) w->var(nuisname.c_str());
	double bf = nuis->getVal();
	double nll_0=nll->getVal();


	TGraph *gr = new TGraph(2*npoints+1);
	for (int i=-1*npoints;i<=npoints;i++){
		nuis->setVal(bf+err*( ((float)i)*nsigma/npoints));
		double nll_v = nll->getVal();
		gr->SetPoint(i+npoints,nuis->getVal(),nll_v-nll_0);
	}

	gr->SetTitle("");
	gr->GetYaxis()->SetTitle("NLL - obs data");
	gr->GetYaxis()->SetTitleOffset(1.1);
	gr->GetXaxis()->SetTitleSize(0.05);
	gr->GetYaxis()->SetTitleSize(0.05);
	gr->GetXaxis()->SetTitle(nuisname.c_str());
	gr->SetLineColor(4);
	gr->SetLineWidth(2);
	gr->SetMarkerStyle(21);
	gr->SetMarkerSize(0.6);
	
	return gr;
	
}

RooAbsPdf *MakeModelNoSignal(RooDataHist *data, RooRealVar *mww, char *name) {
    char *modelName = (char*)calloc(50, sizeof(char));
    strcat(modelName, name);
    strcat(modelName, "model");
    char *signalName = (char*)calloc(50, sizeof(char));
    char *backgroundName = (char*)calloc(50, sizeof(char));
    strcat(signalName, name);
    strcat(signalName, "signal");
    strcat(backgroundName, name);
    strcat(backgroundName, "background");

    char *meanName = (char*)calloc(50, sizeof(char));
    strcat(meanName, name);
    strcat(meanName, "mean");

    char *sigmaName = (char*)calloc(50, sizeof(char));
    strcat(sigmaName, name);
    strcat(sigmaName, "sigma");

    char *decayName = (char*)calloc(50, sizeof(char));
    strcat(decayName, name);
    strcat(decayName, "decay");

    char *weightName = (char*)calloc(50, sizeof(char));
    strcat(weightName, name);
    strcat(weightName, "weight");

    RooRealVar *decay = new RooRealVar(decayName, "decay", -10, 10, "GeV");
    RooExponential *background = new RooExponential(backgroundName, "background", *mww, *decay);

    background->fitTo(*data);
    decay->setConstant();
    return background;
}

void Plot(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	  RooRealVar &iScale,RooRealVar &iRes,RooDataSet *iData=0) { 
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  RooDataSet * lSignal  = iGen->generate(iVar,iNEvents);
  if(iData == 0) {
    iFit->fitTo(*lSignal,Strategy(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
  } else {
    iFit->fitTo(*iData,Strategy(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*iData,Strategy(2));
  }    
  iVar.setBins(30);
  RooPlot *lFrame1 = iVar.frame(RooFit::Title("XXX")) ;
  if(iData == 0) lSignal->plotOn(lFrame1);
  if(iData != 0) iData->plotOn(lFrame1);
  iFit->plotOn(lFrame1);
  TCanvas *iC =new TCanvas("A","A",800,600);
  iC->cd(); lFrame1->Draw();
  iC->SaveAs("Crap.png");
  if(iData != 0) { 
    RooPlot *lFrame2 = iVar.frame(RooFit::Title("XXX")) ;
    iData->plotOn(lFrame2);
    iGen->plotOn(lFrame2);
    TCanvas *iC1 =new TCanvas("B","B",800,600);
    iC1->cd(); lFrame2->Draw();
    iC1->SaveAs("Crap.png");
  }
}

void Fitter::addBernFitModel(RooArgList *pdfs,RooArgList*coeffs, bool isLast)
{
    if (nBernstein <=0 ) return;

    RooArgList * args = new RooArgList();
    for(int b=0;b<nBernstein; b++)
    {
        string name = Form("fitmodel_bern_param_%d",b);
        RooRealVar *tmp = new RooRealVar(name.c_str(), name.c_str(), startBern_[b],0.01,100.);
        vars_ [ name ] = tmp;
        args -> add( *tmp );
    }
    string name = Form("fitmodel_bern_%d",nBernstein);
    RooAbsPdf *bern; 
    if (nBernstein ==1 ) bern= new RooBernsteinFast<1>(name.c_str(),name.c_str(),*x_, *args);
    if (nBernstein ==2 ) bern= new RooBernsteinFast<2>(name.c_str(),name.c_str(),*x_, *args);
    if (nBernstein ==3 ) bern= new RooBernsteinFast<3>(name.c_str(),name.c_str(),*x_, *args);
    if (nBernstein ==4 ) bern= new RooBernsteinFast<4>(name.c_str(),name.c_str(),*x_, *args);
    if (nBernstein ==5 ) bern= new RooBernsteinFast<5>(name.c_str(),name.c_str(),*x_, *args);
    if (nBernstein ==6 ) bern= new RooBernsteinFast<6>(name.c_str(),name.c_str(),*x_, *args);
    if (nBernstein ==7 ) bern= new RooBernsteinFast<7>(name.c_str(),name.c_str(),*x_, *args);

    pdfs->add( *bern );

    if ( not isLast ){ // if I have other pdfs add fraction
        RooRealVar *f = new RooRealVar( Form("fitmodel_bern_frac") ,Form("fitmodel_bern_frac") , 0.3, 0.01,1.0 );
        vars_[ f->GetName() ] = f;
        coeffs -> add ( *f );
    }
} 

void FitterUtils::PlotShape2D(RooDataSet& originDataSet, RooDataSet& genDataSet, RooAbsPdf& shape, string plotsfile, string canvName, RooRealVar& B_plus_M, RooRealVar& misPT)
{
   //**************Prepare TFile to save the plots

   TFile f2(plotsfile.c_str(), "UPDATE");

   //**************Plot Signal Zero Gamma

   TH2F* th2fKey = (TH2F*)shape.createHistogram("th2Shape", B_plus_M, Binning(20), YVar(misPT, Binning(20)));
   cout<<genDataSet.sumEntries()<<endl;
   TH2F* th2fGen = (TH2F*)genDataSet.createHistogram("th2fGen", B_plus_M, Binning(20), YVar(misPT, Binning(20)));

   RooPlot* plotM = B_plus_M.frame();
   originDataSet.plotOn(plotM);
   shape.plotOn(plotM);

   RooPlot* plotMisPT = misPT.frame();
   originDataSet.plotOn(plotMisPT);
   shape.plotOn(plotMisPT);

   TCanvas canv(canvName.c_str(), canvName.c_str(), 800, 800);
   canv.Divide(2,2);
   canv.cd(1); th2fGen->Draw("lego");
   canv.cd(2); th2fKey->Draw("surf");
   canv.cd(3); plotM->Draw();
   canv.cd(4); plotMisPT->Draw();

   canv.Write();

   f2.Close();
}

void chi2(int xmin = 0, int xmax = 200, TString filename="../DsubMC/met2j0bIso2ewkScale_0_200.root", int nparam = 2){
RooAbsData::ErrorType errorType = RooAbsData::SumW2;

file = new TFile(filename);
RooRealVar* h = new RooRealVar("h","h",xmin,xmax); 

//Get Data
TH1D* hData = file->Get("dataih");
hData->SetName("hData");
//hData->Draw();
RooDataHist* data = new RooDataHist("data","data",*h,hData);

//Get Summed MC
TH1D* hMC = file->Get("hh");
hMC->SetName("hMC");
hMC->Draw();
RooDataHist rdhMC("MC","MC",*h,hMC);
RooHistPdf pdfMC("MCpdf","MCpdf",*h,rdhMC);

//make (plot) the curves
RooPlot* hFrame = h->frame(Name("hFrame"));
data->plotOn(hFrame,RooFit::DataError(errorType));
pdfMC.plotOn(hFrame,ProjWData(*data),Components(pdfMC),Name("h_total"));

//Determine Chi^2
double chi2fit = hFrame->chiSquare("h_total", "h_data", nparam);
cout<<"Chi 2 / dof: "<<chi2fit<<endl;

//Determine K-S
double ks = hMC->KolmogorovTest(hData);
cout<<"Kolmogorov-Smirnov: "<<ks<<endl;
}

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

    RooStats::ModelConfig * _mc = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");

    RooRealVar * firstPOI = (RooRealVar*) _mc->GetParametersOfInterest()->first();
    double _limit = mcInt->UpperLimit(*firstPOI);
    cout << "\n95% upper limit on " <<firstPOI->GetName()<<" is : "<<
         _limit <<endl;

    if (filename.size()!=0) {

        std::ofstream aFile;

        // append to file if exists
        aFile.open(filename.c_str(), std::ios_base::app);

        char buf[1024];
        sprintf(buf, "%7.1f   %7.6f", peak, _limit);

        aFile << buf << std::endl;

        // close outfile here so it is safe even if subsequent iterations crash
        aFile.close();

    }

    return _limit;
}

void PDF_GLWADS_Dpi_K3pi::setUncertainties(config c)
{
    switch(c)
    {
    case lumi1fb:
    {
        obsErrSource = "1fb-1, ExpNll/sept2012K3PIResult.root";
        TString File = this->dir+"/ExpNll/sept2012K3PIResult.root";
        TFile *fr = TFile::Open(File);
        RooFitResult *r = (RooFitResult*)fr->Get("fitresult_model_reducedData_binned");
        assert(r);
        for ( int i=0; i<nObs; i++ )
        {
            RooRealVar* pObs = (RooRealVar*)((RooArgList*)observables)->at(i);
            RooRealVar* pRes = (RooRealVar*)r->floatParsFinal().find(obsTmkToMalcolm(pObs->GetName()));
            assert(pRes);
            StatErr[i] = pRes->getError();
        }
        SystErr[0] = 0.010;   // afav_dpi_obs
        SystErr[1] = 0.00011; // rp_dpi_obs
        SystErr[2] = 0.00011; // rm_dpi_obs
        fr->Close();
        delete r;
        delete fr;
        break;
    }
    default:
        cout << "PDF_GLWADS_Dpi_K3pi::setUncertainties() : ERROR : config "+ConfigToTString(c)+" not found." << endl;
        exit(1);
    }
}

int main(){
  
  RooMsgService::instance().setGlobalKillBelow(ERROR);
  
  TFile *bkgFile = TFile::Open("comb_svn/hgg.inputbkgdata_8TeV_MVA.root");
  RooWorkspace *bkgWS = (RooWorkspace*)bkgFile->Get("cms_hgg_workspace");
  RooRealVar *mass = (RooRealVar*)bkgWS->var("CMS_hgg_mass");
  RooDataSet *data = (RooDataSet*)bkgWS->data("data_mass_cat0");

  RooRealVar *p1 = new RooRealVar("p1","p1",-2.,-10.,0.);
  RooRealVar *p2 = new RooRealVar("p2","p2",-0.001,-0.01,0.01);
  RooRealVar *p3 = new RooRealVar("p3","p3",-0.0001,-0.01,0.01);

  //RooPowerLawSum *pow1 = new RooPowerLawSum("pow","pow",*mass,RooArgList(*p1));
  //RooPowerLawSum *pow2 = new RooPowerLawSum("pow","pow",*mass,RooArgList(*p1,*p2));
  //RooPowerLawSum *pow3 = new RooPowerLawSum("pow","pow",*mass,RooArgList(*p1,*p2,*p3));

  RooExponentialSum *pow1 = new RooExponentialSum("pow1","pow1",*mass,RooArgList(*p1));
  RooExponentialSum *pow2 = new RooExponentialSum("pow2","pow2",*mass,RooArgList(*p1,*p2));
  RooExponentialSum *pow3 = new RooExponentialSum("pow3","pow3",*mass,RooArgList(*p1,*p2,*p3));

  TCanvas *canv = new TCanvas();
  RooPlot *frame = mass->frame();
  data->plotOn(frame);
  
  cout << "bus" << endl;
  pow1->fitTo(*data,PrintEvalErrors(-1),PrintLevel(-1),Warnings(-1),Verbose(-1));
  cout << "bus" << endl;
  pow1->plotOn(frame);
  pow2->fitTo(*data);
  pow2->plotOn(frame,LineColor(kRed),LineStyle(kDashed));
  pow3->fitTo(*data);
  pow3->plotOn(frame,LineColor(kGreen),LineStyle(7));
  frame->Draw();
  canv->Print("test.pdf");

  TFile *outFile = new TFile("test.root","RECREATE");

  TTree *tree = new TTree("tree","tree");

  vector<double> mu;
  mu.push_back(1.);
  mu.push_back(2.);
  mu.push_back(3.);
  mu.push_back(4.);
  vector<string> label;
  label.push_back("pol");
  label.push_back("pow");
  label.push_back("lau");
  label.push_back("exp");

  tree->Branch("mu",&mu);
  tree->Branch("label",&label);

  tree->Fill();
  outFile->cd();
  tree->Write();
  outFile->Close();
  return 0;
}

void PDF_GLWADS_DKDpi_K3pi::setObservables(config c)
{
	switch(c)
	{
		case truth:{
					   setObservablesTruth();
					   break;
				   }
		case toy:{
					 setObservablesToy();
					 break;
				 }
		case lumi1fb:{
						 obsValSource = "1fb-1, ExpNll/sept2012K3PIResult.root";
						 TString File = this->dir+"/ExpNll/sept2012K3PIResult.root";
						 TFile *fr = TFile::Open(File);
						 RooFitResult *r = (RooFitResult*)fr->Get("fitresult_model_reducedData_binned");
						 assert(r);
						 TIterator* it = observables->createIterator();
						 while ( RooRealVar* pObs = (RooRealVar*)it->Next() )
						 {
							 RooRealVar* pRes = (RooRealVar*)r->floatParsFinal().find(obsTmkToMalcolm(pObs->GetName()));
							 pObs->setVal(pRes->getVal());
						 }
						 fr->Close();
						 delete r;
						 delete fr;
						 break;
					 }
		case lumi3fb:{
						 obsValSource = "3fb-1 ANA v7 unblind"; // https://twiki.cern.ch/twiki/pub/LHCbPhysics/B2D0K/LHCb-ANA-2014-071-v7.pdf (see Vavas email 04/08/15)

						 // these get transformed over from the new inputs using ExpNll/transportGLWADS_new_to_old.py
						 // in the case of the DK only (robust) combination some of the observables don't exist
						 // usemap as the temp store
						 std::map< TString, double > vals;
						 vals["rkp_k3pi_obs"]      =  0.0793;
						 vals["afav_dk_k3pi_obs"]  =  -0.0004;
						 vals["afav_dpi_k3pi_obs"] =  0.0;
						 vals["rp_dk_k3pi_obs"]    =  0.018369;
						 vals["rm_dk_k3pi_obs"]    =  0.009611;
						 vals["rp_dpi_k3pi_obs"]   =  0.003683;
						 vals["rm_dpi_k3pi_obs"]   =  0.003857;

						 // now can loop the observables and set the values
						 TIterator* it = observables->createIterator();
						 while ( RooRealVar* pObs = (RooRealVar*)it->Next() ){
							 pObs->setVal(vals[pObs->GetName()]);
						 }

						 vals.clear();
						 break;
				}
		default:{
					cout << "PDF_GLWADS_DKDpi_K3pi::setObservables() : ERROR : config "+ConfigToTString(c)+" not found." << endl;
					exit(1);
				}
	}
}

void checkBestFitPoint(std::string workspace, std::string fitFile, bool splusb){
	
	// Open the ws file...
	TFile *fd_=0;
	TFile *fw_=0;
		
	gSystem->Load("$CMSSW_BASE/lib/$SCRAM_ARCH/libHiggsAnalysisCombinedLimit.so");
	gROOT->SetBatch(true);
	gStyle->SetOptFit(0);
	gStyle->SetOptStat(0);
	gStyle->SetPalette(1,0);

	fw_ =  TFile::Open(workspace.c_str());
	w   = (RooWorkspace*) fw_->Get("w");
  w->Print();
	RooDataSet *data = (RooDataSet*) w->data("data_obs");
  if (splusb) {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig");
  } else {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig_bonly");
  }
	std::cout << "make nll"<<std::endl;
	nll = mc_s->GetPdf()->createNLL(
		*data,RooFit::Constrain(*mc_s->GetNuisanceParameters())
		,RooFit::Extended(mc_s->GetPdf()->canBeExtended()));
	
	// Now get the best fit result
	fd_ =  TFile::Open(fitFile.c_str());
	RooFitResult *fit;
  if (splusb) {
    fit =(RooFitResult*)fd_->Get("fit_s");
  } else {
    fit =(RooFitResult*)fd_->Get("fit_b");
  }
	RooArgSet fitargs = fit->floatParsFinal();
	
	std::cout << "Got the best fit values" <<std::endl;		
	w->saveSnapshot("bestfitall",fitargs,true);
	
  TString filename;
  if (splusb) {
    filename = "minimum_s.pdf";
  } else {
    filename = "minimum_b.pdf";
  }
	// Now make the plots!	
	TCanvas *c = new TCanvas("c","",600,600);
	c->SaveAs((filename+"["));

	TIterator* iter(fitargs->createIterator());
        for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
                 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);      
                 std::string name = rrv->GetName();
		 TGraph *gr = graphLH(name,rrv->getError());
		 gr->Draw("ALP");
		 c->SaveAs((filename+"["));
	}
	c->SaveAs((filename+"]"));
}

//
// set value and range for a variable in the workspace
//
void setValRange (RooWorkspace* workspace, const char* name, double val, double vmin, double vmax)
{
  RooRealVar* var = workspace->var(name);
  if ( var ) {
    if ( vmax>vmin )  var->setRange(vmin,vmax);
    var->setVal(val);
  }
}