本文整理汇总了C++中RooAddPdf类的典型用法代码示例。如果您正苦于以下问题：C++ RooAddPdf类的具体用法？C++ RooAddPdf怎么用？C++ RooAddPdf使用的例子？那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了RooAddPdf类的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: mes
RooRealVar mes("mes","m_{ES} (GeV)",5.20,5.30) ;
// --- Parameters ---
RooRealVar sigmean("sigmean","B^{#pm} mass",5.28,5.20,5.30) ;
RooRealVar sigwidth("sigwidth","B^{#pm} width",0.0027,0.001,1.) ;
// --- Build Gaussian PDF ---
RooGaussian signal("signal","signal PDF",mes,sigmean,sigwidth) ;
// --- Build Argus background PDF ---
RooRealVar argpar("argpar","argus shape parameter",-20.0,-100.,-1.) ;
RooArgusBG background("background","Argus PDF",mes,RooConst(5.291),argpar) ;
// --- Construct signal+background PDF ---
RooRealVar nsig("nsig","#signal events",100,0.,10000) ;
RooRealVar nbkg("nbkg","#background events",100,0.,10000) ;
RooAddPdf model("model","g+a",RooArgList(signal,background),RooArgList(nsig,nbkg)) ;
// --- Generate a toyMC sample from composite PDF ---
RooDataSet *data = model.generate(mes,2000) ;
TFile* f= new TFile("rooFit.root","RECREATE");
TCanvas* cc = new TCanvas("cc","cc",800,600) ;
//TH1F* h=(TH1F*)f->Get("histo__mes");
//h->Rebin(2);
//RooDataHist* data=new RooDataHist("data", "data", mes, h);
//TFile* f= new TFile("rooFit.root","RECREATE");
//TH1F* h=data->createHistogram("histo",mes,Binning(100,5.2,5.3));
//h->Draw();
//h->Write();
//f->Close();示例2: addNuisanceWithToys
//.........这里部分代码省略.........
//Perform the tail fit with the alternative fit function (once initially, before allowing tail fit to float in toy fit).
//=============================================================================================================================================
RooFitResult *lRFit1 = 0;
//lRFit1=lFit1->fitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(iFirst,iLast),RooFit::Strategy(0));
lRFit1=lFit1->fitTo(*pH0,RooFit::Save(kTRUE),RooFit::Range(200,1500),RooFit::Strategy(0));
//Generate the background pdf corresponding to the result of the alternative tail fit
RooGenericPdf *lFit1Final = 0; lFit1Final = new RooGenericPdf("genPdf","exp(-m/(a1+b1*m))",RooArgList(lM,lA1,lB1));
if(iFitModel1 == 1) lFit1Final = new RooGenericPdf("genPdf","exp(-a1*pow(m,b1))",RooArgList(lM,lA1,lB1));
if(iFitModel1 == 2) lFit1Final = new RooGenericPdf("genPdf","a1*exp(b1*m)",RooArgList(lM,lA1,lB1));
if(iFitModel1 == 3) lFit1Final = new RooGenericPdf("genPdf","a1/pow(m,b1)",RooArgList(lM,lA1,lB1));
// lA1.removeRange();
// lB1.removeRange();
//=============================================================================================================================================
//Define RooRealVar for the normalization of the signal and background, starting from the initial integral of the input histograms
lM.setRange(300,1500);
RooRealVar lNB("nb","nb",lNB0,0,10000);
RooRealVar lNSig("nsig","nsig",lNSig0,-1000,1000);
//Define a PDF for the signal histogram lSig
RooDataHist *pS = new RooDataHist("sigH","sigH",RooArgList(lM),lSig);
RooHistPdf *lSPdf = new RooHistPdf ("sigPdf","sigPdf",lM,*pS);
//Define generator and fit functions for the RooMCStudy
RooAddPdf *lGenMod = new RooAddPdf ("genmod","genmod",RooArgList(*lFitFinal ,*lSPdf),RooArgList(lNB,lNSig));
RooAddPdf *lFitMod = new RooAddPdf ("fitmod","fitmod",RooArgList(*lFit1Final,*lSPdf),RooArgList(lNB,lNSig));
//Generate plot of the signal and background models going into the toy generation
RooPlot* plot=lM.frame();
lGenMod->plotOn(plot);
lGenMod->plotOn(plot,RooFit::Components(*lSPdf),RooFit::LineColor(2));
TCanvas* lC11 = new TCanvas("pdf","pdf",600,600) ;
lC11->cd();
plot->Draw();
lC11->SaveAs(("SBModel_"+iBkg+"_" + iDir + "_" + iEnergy+".pdf").c_str());
std::cout << "===================================================================================================================================================" <<std::endl;
std::cout << "FIT PARAMETERS BEFORE ROOMCSTUDY: lA: " << lA.getVal() << " lB: " << lB.getVal() << " lA1: " << lA1.getVal() << " lB1: " << lB1.getVal() << std::endl;
std::cout << "===================================================================================================================================================" <<std::endl;
RooMCStudy *lToy = new RooMCStudy(*lGenMod,lM,RooFit::FitModel(*lFitMod),RooFit::Binned(kTRUE),RooFit::Silence(),RooFit::Extended(kTRUE),RooFit::Verbose(kTRUE),RooFit::FitOptions(RooFit::Save(kTRUE),RooFit::Strategy(0)));
// Generate and fit iNToys toy samples
std::cout << "Number of background events: " << lNB0 << " Number of signal events: " << lNSig0 << " Sum: " << lNB0+lNSig0 << std::endl;
//=============================================================================================================================================
// Generate and fit toys
//=============================================================================================================================================
lToy->generateAndFit(iNToys,lNB0+lNSig0,kTRUE);
std::cout << "===================================================================================================================================================" <<std::endl;
std::cout << "FIT PARAMETERS AFTER ROOMCSTUDY: lA: " << lA.getVal() << " lB: " << lB.getVal() << " lA1: " << lA1.getVal() << " lB1: " << lB1.getVal() << std::endl;
std::cout << "===================================================================================================================================================" <<std::endl;示例3: higgsMassFit
//.........这里部分代码省略.........
//-----------------------------------
// draw the background + signal stack
if( drawPlots )
{
TCanvas* c1 = new TCanvas("BKGShapeStack","BKGShapeStack");
c1 -> SetGridx();
c1 -> SetGridy();
BKGShapeStack -> Draw("HIST");
SIGShapeHisto -> Draw("HIST,same");
char pngFileName[50];
sprintf(pngFileName,"BKGShapeStack_H%d_%s.png",mH,mode.c_str());
c1 -> Print(pngFileName,"png");
}
//---------------------------------
// define the bkg shape with roofit
std::cout << ">>> Define the background pdf" << std::endl;
RooKeysPdf** rooBKGPdf = new RooKeysPdf*[nBKG];
RooRealVar** rooNBKG = new RooRealVar*[nBKG];
RooRealVar* rooNBKGTot = new RooRealVar("rooNBKGTot","",BKGTotShapeHisto->Integral(),0.,1000000.);
for(int i = 0; i < nBKG; ++i)
{
rooBKGPdf[i] = new RooKeysPdf(("rooBKGPdf_"+BKGShortNames[i]).c_str(),"",x,*rooBKGDataSet[i]);
rooNBKG[i] = new RooRealVar(("rooNBKG_"+BKGShortNames[i]).c_str(),"",BKGShapeHisto[i]->Integral(),BKGShapeHisto[i]->Integral()),BKGShapeHisto[i]->Integral();
}
RooAddPdf* rooBKGTotPdf = new RooAddPdf("rooBKGTotPdf","",RooArgList(*rooBKGPdf[0],*rooBKGPdf[1]),RooArgList(*rooNBKG[0],*rooNBKG[1]));
//---------------------------------
// define the sig shape with roofit
std::cout << ">>> Define the signal pdf" << std::endl;
RooKeysPdf* rooSIGPdf = new RooKeysPdf("rooSIGPdf","",x,*rooSIGDataSet);
RooRealVar* rooNSIG = new RooRealVar("rooNSIG","",1.,-1000000.,1000000.);
//---------------------------------
// define the tot shape with roofit
std::cout << ">>> Define the total pdf" << std::endl;
RooAddPdf* rooTotPdf = NULL;
if( mode == "exclusion") rooTotPdf = new RooAddPdf("rooTotPdf","",RooArgList(*rooBKGTotPdf),RooArgList(*rooNBKGTot));
if( mode == "discovery") rooTotPdf = new RooAddPdf("rooTotPdf","",RooArgList(*rooBKGTotPdf,*rooSIGPdf),RooArgList(*rooNBKGTot,*rooNSIG));
//----
// plot
if( drawPlots )
{
TCanvas* c2 = new TCanvas("rooTotPdf","rooTotPdf");
c2 -> SetGridx();
c2 -> SetGridy();示例4: fitTopMass
void fitTopMass(const char * inFileName)
{
gROOT->SetBatch(true);
gSystem->Load("libRooFit");
TFile* inFile = TFile::Open(inFileName);
RooRealVar x("mass","M_{top}",100,200);
using namespace std;
vector<TH1F*> hists;
vector<int> nTotal;
std::vector<std::pair<float, float> > range;
const char* types[] = { "top1","top2","top"};
for( auto type : types ) {
hists.push_back ( (TH1F*)inFile->Get(TString::Format("%s_mass_GenJet",type)) ) ;
hists.push_back ( (TH1F*)inFile->Get(TString::Format("%s_mass_Jet",type)) ) ;
hists.push_back ( (TH1F*)inFile->Get(TString::Format("%s_mass_btagged_Jet",type)) ) ;
hists.push_back ( (TH1F*)inFile->Get(TString::Format("%s_mass_charged_Jet",type)) ) ;
hists.push_back ( (TH1F*)inFile->Get(TString::Format("%s_mass_btagged_charged_Jet",type)) ) ;
}
for( int i= 0 ; i<hists.size() ; i++) {
if ( hists[i] == nullptr) std::cout<<"Error! Hist is nullptr."<<std::endl;
}
for( int i = 0 ; i < hists.size() ; i++) {
nTotal.push_back(hists[i]->Integral());
}
// For top1,
for(int i= 0 ; i < 5 ; i++) {
range.push_back( make_pair(150,180) );
}
// For top2,
for(int i= 0 ; i < 5 ; i++) {
range.push_back( make_pair(150,180) );
}
// For all top,
for(int i= 0 ; i < 5 ; i++) {
range.push_back( make_pair(150,180) );
}
for( int i= 0 ; i< hists.size() ; i++) {
std::cout<<"hist["<<i<<"]"<<std::endl;
RooDataHist* dh = new RooDataHist("dh","data histogram",x,hists[i]);
auto gaus_mean = RooRealVar("gaus_mean","gaus_mean",172.5, 100,200);
auto gaus_sigma = RooRealVar("gaus_sigma","gaus_sigma",1, 0,100);
auto nGaus = RooRealVar("nGaus","nGaus",nTotal[i]*0.2, 0., nTotal[i]);
auto gaus_pdf = RooGaussian("gaus_pdf","Gaussian p.d.f",x,gaus_mean, gaus_sigma);
auto expo_tau = RooRealVar("exp_const","exp_const",0,-10000,10000);
auto nExp = RooRealVar("nExp","nExp",1,0,nTotal[i]);
auto expo_pdf = RooExponential("bkg","bkg p.d.f",x,expo_tau);
auto CB_mean = RooRealVar("CB_mean","CB_mean",172.5,100,200);
auto CB_sigma = RooRealVar("CB_sigma","CB_sigma",5,0,100);
auto CB_alpha = RooRealVar("CB_alpha","CB_alpha",1,0,100);
auto CB_n = RooRealVar("CB_n","CB_n",1,0,100);
auto nCB = RooRealVar("nCB","nCB",nTotal[i]*0.8,0.,nTotal[i]);
auto CB_pdf = RooCBShape("sig","signal p.d.f",x,CB_mean, CB_sigma, CB_alpha, CB_n );
auto BW_mean = RooRealVar("BW_mean","BW_mean",172.5,100,200);
auto BW_sigma = RooRealVar("BW_sigma","BW_sigma",10,0,100);
auto nBW = RooRealVar("nBW","nBW",nTotal[i]*0.7, 0, nTotal[i]);
auto BW_pdf = RooBreitWigner("sig","signal p.d.f",x,BW_mean, BW_sigma );
RooFormulaVar mirrorX("neg_x","-mass", RooArgSet(x));
//RooRealVar landau_mean("lan_mean","Landau_mean",150.5,100,200);
RooRealVar landau_mean("lan_mean","Landau_mean",-160,-200,-100);
RooRealVar landau_sigma("lan_sigma","Landau_sigma",2.5, 0,100);
RooRealVar nLandau("nLan","nlan",nTotal[i]*0.3, 0, nTotal[i]);
RooLandau landau_pdf("lx","lx",mirrorX,landau_mean, landau_sigma);
//TF1* f1 = new TF1("landau","nLan*ROOT::Math::landau_pdf(-x,lan_sigma,-lan_mean)",100,300);
//RooAbsPdf* invlandau_pdf = RooFit::bindPdf(f1,x,nLandau, landau_mean,landau_sigma);
//RooGenericPdf invlandau_pdf("invLandau","ROOT::Math::landau_pdf(-x,lan_sigma,-lan_mean)",RooArgSet(x,landau_sigma,landau_mean));
//sig_pdf = cb_pdf;
//bkg_pdf = bw_pdf;
//nsig = nCB;
//nbkg = nBW;
//RooFFTConvPdf model("lxg","CB (x) gauss",x,CB_pdf,gaus_pdf);
//RooFFTConvPdf model("lxg","CB (x) BW",x,CB_pdf,BW_pdf);
//auto model = RooAddPdf("model","model",RooArgList(landau),RooArgList(nlandau));
//auto model = RooAddPdf("model","model",RooArgList(sig_pdf, bkg_pdf),RooArgList(nsig, nbkg));
//auto model = RooAddPdf("model","model",RooArgList(CB_pdf, gaus_pdf),RooArgList(nCB, nGaus));
RooAddPdf* model;
RooFitResult * fitResult;
TPaveText* t1 = new TPaveText(0.15,0.5,0.45,0.8,"brNDC");
//auto model = RooAddPdf("model","model",RooArgList(CB_pdf, expo_pdf),RooArgList(nCB, nExp));
//auto model = RooAddPdf("model","model",RooArgList(CB_pdf, BW_pdf),RooArgList(nCB, nBW));
//auto model = RooAddPdf("model","model",RooArgList(CB_pdf),RooArgList(nCB));
//auto model = RooAddPdf("model","model",RooArgList(cb_pdf, gaus_pdf),RooArgList(nCB, nGaus));
//auto model = sig_pdf;
//.........这里部分代码省略.........
示例5: fit_mass
void fit_mass(TString fileN="") {//suffix added before file extension, e.g., '.pdf'
TString placeholder;//to add strings before using them, e.g., for saving text files
gROOT->SetBatch(kTRUE);
gROOT->ProcessLine(".x /afs/cern.ch/user/m/mwilkins/cmtuser/src/lhcbStyle.C");
// gStyle->SetPadTickX(1);
// gStyle->SetPadTickY(1);
// gStyle->SetPadLeftMargin(0.15);
// gStyle->SetTextSize(0.3);
// //open file and get histogram
// TFile *inHistos = new TFile("/afs/cern.ch/work/m/mwilkins/Lb2JpsiLtr/data/histos_data.root", "READ");
// TH1F * h100 = (TH1F*)inHistos->Get("h70");
// cout<<"data histogram gotten"<<endl;
//unbinned
TFile *hastree = new TFile("/afs/cern.ch/work/m/mwilkins/Lb2JpsiLtr/data/cutfile_Optimized.root", "READ");
TTree * h100 = (TTree*)hastree->Get("mytree");
cout<<"tree gotten"<<endl;
TFile *SMChistos= new TFile("/afs/cern.ch/work/m/mwilkins/Lb2JpsiLtr/MC/withKScut/histos_SMCfile_fullMC.root", "READ");
cout<<"SMC file opened"<<endl;
TH1F *SMCh = (TH1F*)SMChistos->Get("h00");
cout<<"SMC hist gotten"<<endl;
RooRealVar *mass = new RooRealVar("Bs_LOKI_MASS_JpsiConstr","m(J/#psi #Lambda)",4100,6100,"MeV");
mass->setRange("bkg1",4300,4800);
mass->setRange("bkg2",5700,5950);
mass->setRange("bkg3",4300,5500);
mass->setRange("bkg4",5100,5500);
mass->setRange("L",5350,5950);
mass->setRange("tot",4300,5950);
cout<<"mass declared"<<endl;
// RooDataHist *data = new RooDataHist("data","1D",RooArgList(*mass),h100);
//unbinned
RooDataSet *data = new RooDataSet("data","1D",h100,*mass);
cout<<"data declared"<<endl;
RooDataHist *SMC = new RooDataHist("SMC","1D",RooArgList(*mass),SMCh);
cout<<"SMC hist assigned to RooDataHist"<<endl;
// Construct Pdf Model
// /\0
//gaussian
RooRealVar mean1L("mean1L","/\\ gaus 1: mean",5621.103095,5525,5700);
RooRealVar sig1L("sig1L","/\\ gaus 1: sigma",6.898126,0,100);
RooGaussian gau1L("gau1L","#Lambda signal: gaussian 1",*mass,mean1L,sig1L);
RooFormulaVar mean2L("mean2L","@0",mean1L);
RooRealVar sig2L("sig2L","/\\ gaus 2: sigma",14.693117,0,100);
RooGaussian gau2L("gau2L","#Lambda signal: gaussian 2",*mass,mean2L,sig2L);
RooRealVar f1L("f1L","/\\ signal: fraction gaussian 1",0.748776,0,1);
RooAddPdf sigL("sigL","#Lambda signal",RooArgList(gau1L,gau2L),RooArgList(f1L));
// //CB
// RooRealVar mean3L("mean3L","/\\ CB: mean",5621.001,5525,5700);
// RooRealVar sig3L("sig3L","/\\ CB: sigma",5.161,0,100);
// RooRealVar alphaL3("alphaL3","/\\ CB: alpha",2.077,0,1000);
// RooRealVar nL3("nL1","/\\ CB: n",0.286,0,1000);
// RooCBShape CBL("CBL","#Lambda signal: CB",*mass,mean3L,sig3L,alphaL3,nL3);
// RooRealVar mean4L("mean4L","/\\ gaus: mean",5621.804,5525,5700);
// RooRealVar sig4L("sig4L","/\\ gaus: sigma",10.819,0,100);
// RooGaussian gauL("gauL","#Lambda signal: gaussian",*mass,mean4L,sig4L);
// RooRealVar f1L("f1L","/\\ signal: fraction CB",0.578,0,1);
// RooAddPdf sigL("sigL","#Lambda signal",RooArgList(CBL,gauL),RooArgList(f1L));
// sigma0
//using RooHistPdf from MC--no need to build pdf here
RooHistPdf sigS = makeroohistpdf(SMC,mass,"sigS","#Sigma^{0} signal (RooHistPdf)");
// /\*
cout<<"Lst stuff"<<endl;
RooRealVar meanLst1("meanLst1","/\\*(misc.): mean1",5011.031237,4900,5100);
RooRealVar sigLst1("sigLst1","/\\*(misc.): sigma1",70.522092,0,100);
RooRealVar meanLst2("mean5Lst2","/\\*(1405): mean2",5245.261703,5100,5350);
RooRealVar sigLst2("sigLst2","/\\*(1405): sigma2",64.564763,0,100);
RooRealVar alphaLst2("alphaLst2","/\\*(1405): alpha2",29.150301);
RooRealVar nLst2("nLst2","/\\*(1405): n2",4.615817,0,50);
RooGaussian gauLst1("gauLst1","#Lambda*(misc.), gaus",*mass,meanLst1,sigLst1);
RooCBShape gauLst2("gauLst2","#Lambda*(1405), CB",*mass,meanLst2,sigLst2,alphaLst2,nLst2);
// RooRealVar fLst1("fLst1","/\\* bkg: fraction gaus 1",0.743,0,1);
// RooAddPdf bkgLst("bkgLst","#Lambda* signal",RooArgList(gauLst1,gauLst2),RooArgList(fLst1));
//Poly func BKG mass
// RooRealVar b0("b0","Background: Chebychev b0",-1.071,-10000,10000);
RooRealVar b1("b1","Background: Chebychev b1",-1.323004,-10,-0.00000000000000000000001);
RooRealVar b2("b2","Background: Chebychev b2",0.145494,0,10);
RooRealVar b3("b3","Background: Chebychev b3",-0.316,-10000,10000);
RooRealVar b4("b4","Background: Chebychev b4",0.102,-10000,10000);
RooRealVar b5("b5","Background: Chebychev b5",0.014,-10000,10000);
RooRealVar b6("b6","Background: Chebychev b6",-0.015,-10000,10000);
RooRealVar b7("b7","Background: Chebychev b7",0.012,-10000,10000);
RooArgList bList(b1,b2);
RooChebychev bkg("bkg","Background", *mass, bList);
// TF1 *ep = new TF1("ep","[2]*exp([0]*x+[1]*x*x)",4300,5950);
// ep->SetParameter(0,1);
// ep->SetParameter(1,-1);
// ep->SetParameter(2,2000);
// ep->SetParName(0,"a");
// ep->SetParName(1,"b");
// ep->SetParName(2,"c");
// RooRealVar a("a","Background: Coefficent of x",1,-10000,10000);
// RooRealVar b("b","Background: Coefficent of x*x",-1,-10000,10000);
// RooRealVar c("c","Background: Coefficent of exp()",2000,-10000,10000);
// RooTFnPdfBinding bkg("ep","ep",ep,RooArgList(*mass,a,b));
//.........这里部分代码省略.........
示例6: main
int main (int argc, char **argv) {
TFile *tf = TFile::Open("tmp/DataSets.root");
RooWorkspace *w = (RooWorkspace*)tf->Get("w");
RooDataSet *Data = (RooDataSet*)w->data("Data2011")->Clone("Data");
Data->append( *((RooDataSet*)w->data("Data2012")) );
RooDataSet *Bs2Kst0Kst0_MC = (RooDataSet*)w->data("Bs2Kst0Kst0_MC2011")->Clone("Bs2KstKst0_MC");
Bs2Kst0Kst0_MC->append( *((RooDataSet*)w->data("Bs2Kst0Kst0_MC2012")) );
RooDataSet *Bs2Kst0Kst01430_MC = (RooDataSet*)w->data("Bs2Kst0Kst01430_MC2011")->Clone("Bs2KstKst0_MC");
Bs2Kst0Kst01430_MC->append( *((RooDataSet*)w->data("Bs2Kst0Kst01430_MC2012")) );
RooDataSet *Bs2Kst01430Kst01430_MC = (RooDataSet*)w->data("Bs2Kst01430Kst01430_MC2011")->Clone("Bs2KstKst0_MC");
Bs2Kst01430Kst01430_MC->append( *((RooDataSet*)w->data("Bs2Kst01430Kst01430_MC2012")) );
RooDataSet *Bd2Kst0Kst0_MC = (RooDataSet*)w->data("Bd2Kst0Kst0_MC2011")->Clone("Bs2KstKst0_MC");
Bd2Kst0Kst0_MC->append( *((RooDataSet*)w->data("Bd2Kst0Kst0_MC2012")) );
RooDataSet *Bd2PhiKst0_MC = (RooDataSet*)w->data("Bd2PhiKst0_MC2011")->Clone("Bs2KstKst0_MC");
Bd2PhiKst0_MC->append( *((RooDataSet*)w->data("Bd2PhiKst0_MC2012")) );
RooDataSet *Bs2PhiKst0_MC = (RooDataSet*)w->data("Bs2PhiKst0_MC2011")->Clone("Bs2KstKst0_MC");
Bs2PhiKst0_MC->append( *((RooDataSet*)w->data("Bs2PhiKst0_MC2012")) );
RooDataSet *Bd2RhoKst0_MC = (RooDataSet*)w->data("Bd2RhoKst0_MC2011")->Clone("Bs2KstKst0_MC");
Bd2RhoKst0_MC->append( *((RooDataSet*)w->data("Bd2RhoKst0_MC2012")) );
RooDataSet *Lb2ppipipi_MC = (RooDataSet*)w->data("Lb2ppipipi_MC2011")->Clone("Bs2KstKst0_MC");
Lb2ppipipi_MC->append( *((RooDataSet*)w->data("Lb2ppipipi_MC2012")) );
RooDataSet *Lb2pKpipi_MC = (RooDataSet*)w->data("Lb2pKpipi_MC2011")->Clone("Bs2KstKst0_MC");
Lb2pKpipi_MC->append( *((RooDataSet*)w->data("Lb2pKpipi_MC2012")) );
w->import(*Data);
w->import(*Bs2Kst0Kst0_MC);
w->import(*Bs2Kst0Kst01430_MC);
w->import(*Bs2Kst01430Kst01430_MC);
w->import(*Bd2Kst0Kst0_MC);
w->import(*Bd2PhiKst0_MC);
w->import(*Bs2PhiKst0_MC);
w->import(*Bd2RhoKst0_MC);
w->import(*Lb2ppipipi_MC);
w->import(*Lb2pKpipi_MC);
RooRealVar *mass = (RooRealVar*)w->var("B_s0_DTF_B_s0_M");
fitIpatia( w, "bs2kstkst_mc", "Bs2KstKst0_MC");
// Make the PDF here
RooRealVar *p1 = new RooRealVar("p1","p1",-0.002,-0.004,0.);
RooExponential *exp = new RooExponential("exp","exp",*mass,*p1);
//RooRealVar *m1 = new RooRealVar("m1","m1",5320,5380);
//RooRealVar *s1 = new RooRealVar("s1","s1",1,20);
//RooGaussian *sig = new RooGaussian("sig","sig",*mass,*m1,*s1);
RooRealVar *m2 = new RooRealVar("m2","m2",5320,5380);
RooRealVar *s2 = new RooRealVar("s2","s2",1,20);
RooGaussian *sig_bd = new RooGaussian("sig_bd","sig_bd",*mass,*m2,*s2);
//
RooRealVar *bs2kstkst_l = new RooRealVar( "bs2kstkst_l" ,"", -5, -20, -1.);
RooConstVar *bs2kstkst_zeta = new RooConstVar( "bs2kstkst_zeta","",0. );
RooConstVar *bs2kstkst_fb = new RooConstVar( "bs2kstkst_fb" ,"",0. );
RooRealVar *bs2kstkst_sigma = new RooRealVar( "bs2kstkst_sigma","",15 ,10 ,20 );
RooRealVar *bs2kstkst_mu = new RooRealVar( "bs2kstkst_mu" ,"",5350 ,5380 );
RooRealVar *bs2kstkst_a = new RooRealVar( "bs2kstkst_a" ,"",2.5 , 0 ,10 );
RooRealVar *bs2kstkst_n = new RooRealVar( "bs2kstkst_n" ,"",2.5 , 0 ,10 );
RooRealVar *bs2kstkst_a2 = new RooRealVar( "bs2kstkst_a2" ,"",2.5 , 0 ,10 );
RooRealVar *bs2kstkst_n2 = new RooRealVar( "bs2kstkst_n2" ,"",2.5 , 0 ,10 );
RooIpatia2 *sig = new RooIpatia2("sig","sig",*mass,*bs2kstkst_l,*bs2kstkst_zeta,*bs2kstkst_fb,*bs2kstkst_sigma,*bs2kstkst_mu,*bs2kstkst_a,*bs2kstkst_n,*bs2kstkst_a2,*bs2kstkst_n2);
RooRealVar *bkg_y = new RooRealVar("bkg_y","bkg_y",10e3,10e5);
RooRealVar *sig_y = new RooRealVar("sig_y","sig_y",0,20e3);
RooRealVar *sig_bd_y = new RooRealVar("sig_bd_y","sig_bd_y",0,3000);
RooArgList *pdfs = new RooArgList();
RooArgList *yields = new RooArgList();
pdfs->add( *exp );
pdfs->add( *sig );
pdfs->add( *sig_bd );
yields->add( *bkg_y );
yields->add( *sig_y );
yields->add( *sig_bd_y );
RooAddPdf *pdf = new RooAddPdf("pdf","pdf",*pdfs,*yields);
pdf->fitTo(*Data, Extended() );
RooPlot *plot = mass->frame();
Data->plotOn(plot);
// set fit params constant;
pdf->plotOn(plot);
TCanvas *c = new TCanvas();
plot->Draw();
//.........这里部分代码省略.........
示例7: fitHist
double* fitHist(TCanvas *iC, bool iDoMu,int iPlot, std::string iName,TH1D *iData,TH1D *iW,TH1D *iEWK,TH1D *iAntiData,TH1D *iAntiW,TH1D *iAntiEWK,const Double_t METMAX,const Int_t NBINS,const Double_t lumi,const Int_t Ecm,int iAltQCD) {
//
// Declare fit parameters for signal and background yields
// Note: W signal and EWK+top PDFs are constrained to the ratio described in MC
//
RooRealVar nSig(("nSig"+iName).c_str(),("nSig"+iName).c_str(),0.7*(iData->Integral()),0,iData->Integral());
RooRealVar nQCD(("nQCD"+iName).c_str(),("nQCD"+iName).c_str(),0.3*(iData->Integral()),0,iData->Integral());
RooRealVar cewk(("cewk"+iName).c_str(),("cewk"+iName).c_str(),0.1,0,5) ;
cewk.setVal(iEWK->Integral()/iW->Integral());
cewk.setConstant(kTRUE);
RooFormulaVar nEWK(("nEWK"+iName).c_str(),("nEWK"+iName).c_str(),("cewk"+iName+"*nSig"+iName).c_str(),RooArgList(nSig,cewk));
RooRealVar nAntiSig(("nAntiSig"+iName).c_str(),("nAntiSig"+iName).c_str(),0.05*(iAntiData->Integral()),0,iAntiData->Integral());
RooRealVar nAntiQCD(("nAntiQCD"+iName).c_str(),("nAntiQCD"+iName).c_str(),0.9 *(iAntiData->Integral()),0,iAntiData->Integral());
RooRealVar dewk (("dewk" +iName).c_str(),("dewk" +iName).c_str(),0.1,0,5) ;
dewk.setVal(iAntiEWK->Integral()/iAntiW->Integral());
dewk.setConstant(kTRUE);
RooFormulaVar nAntiEWK(("nAntiEWK"+iName).c_str(),("nAntiEWK"+iName).c_str(),("dewk"+iName+"*nAntiSig"+iName).c_str(),RooArgList(nAntiSig,dewk));
//
// Construct PDFs for fitting
//
RooRealVar pfmet(("pfmet"+iName).c_str(),("pfmet"+iName).c_str(),0,METMAX);
pfmet.setBins(NBINS);
// Signal PDFs
RooDataHist wMet (("wMET" +iName).c_str(), ("wMET" +iName).c_str(), RooArgSet(pfmet),iW);
RooHistPdf pdfW(( "w"+iName).c_str(), ( "w"+iName).c_str(), pfmet, wMet, 1);
RooDataHist awMet (("awMET"+iName).c_str(), ("awMET"+iName).c_str(), RooArgSet(pfmet),iAntiW);
RooHistPdf apdfW(("aw"+iName).c_str(), ("aw"+iName).c_str(), pfmet,awMet, 1);
// EWK+top PDFs
RooDataHist ewkMet (("ewkMET" +iName).c_str(),( "ewkMET"+iName).c_str(), RooArgSet(pfmet),iEWK);
RooHistPdf pdfEWK (( "ewk"+iName).c_str(),( "ewk"+iName).c_str(), pfmet,ewkMet, 1);
RooDataHist aewkMet(("aewkMET"+iName).c_str(),("aewkMET"+iName).c_str(), RooArgSet(pfmet),iAntiEWK);
RooHistPdf apdfEWK (("aewk"+iName).c_str(),("aewk"+iName).c_str(), pfmet,aewkMet, 1);
// QCD Pdfs
CPepeModel0 qcd0 (("qcd0" +iName).c_str(),pfmet);
CPepeModel1 qcd1 (("qcd1" +iName).c_str(),pfmet);
CPepeModel2 qcd2 (("qcd2" +iName).c_str(),pfmet);
CPepeModel0 aqcd0(("aqcd0"+iName).c_str(),pfmet);
CPepeModel1 aqcd1(("aqcd1"+iName).c_str(),pfmet,qcd1.sigma);
CPepeModel2 aqcd2(("aqcd2"+iName).c_str(),pfmet);
RooGenericPdf *lQCD = qcd0.model;
RooGenericPdf *lAQCD = aqcd0.model;
if(iDoMu) lQCD = qcd1.model;
if(iDoMu) lAQCD = aqcd1.model;
if(iAltQCD == 0) lQCD = qcd0.model;
if(iAltQCD == 1) lQCD = qcd1.model;
if(iAltQCD == 2) lQCD = qcd2.model;
if(iAltQCD == 0) lAQCD = aqcd0.model;
if(iAltQCD == 1) lAQCD = aqcd1.model;
if(iAltQCD == 2) lAQCD = aqcd2.model;
// Signal + Background PDFs
RooAddPdf pdfMet (("pdfMet"+iName).c_str(), ("pdfMet" +iName).c_str(), RooArgList(pdfW ,pdfEWK ,*lQCD), RooArgList(nSig,nEWK,nQCD));
RooAddPdf apdfMet(("apdfMet"+iName).c_str(),("apdfMet"+iName).c_str(), RooArgList(apdfW,apdfEWK,*lAQCD), RooArgList(nAntiSig,nAntiEWK,nAntiQCD));
// PDF for simultaneous fit
RooCategory rooCat("rooCat","rooCat");
rooCat.defineType("Select");
rooCat.defineType("Anti");
RooSimultaneous pdfTotal("pdfTotal","pdfTotal",rooCat);
pdfTotal.addPdf(pdfMet, "Select");
if(iDoMu) pdfTotal.addPdf(apdfMet,"Anti");
// Perform fits
RooDataHist dataMet (("dataMet"+iName).c_str(),("dataMet"+iName).c_str(), RooArgSet(pfmet),iData);
RooDataHist antiMet (("antiMet"+iName).c_str(),("antiMet"+iName).c_str(), RooArgSet(pfmet),iAntiData);
RooDataHist dataTotal(("data" +iName).c_str(),("data" +iName).c_str(), RooArgList(pfmet), Index(rooCat),
Import("Select", dataMet),
Import("Anti", antiMet));
RooFitResult *fitRes = 0;
bool runMinos = kTRUE;
if(iPlot == 0 || iPlot == 3) runMinos = kFALSE; //Remove Minos when running toys (too damn slow)
if(!iDoMu) fitRes = pdfMet .fitTo(dataMet ,Extended(),Minos(runMinos),Save(kTRUE));
if( iDoMu) fitRes = pdfTotal.fitTo(dataTotal,Extended(),Minos(runMinos),Save(kTRUE));
double *lResults = new double[16];
lResults[0] = nSig.getVal();
lResults[1] = nEWK.getVal();
lResults[2] = nQCD.getVal();
lResults[3] = nAntiSig.getVal();
lResults[4] = nAntiEWK.getVal();
lResults[5] = nAntiQCD.getVal();
if(!iDoMu) lResults[6] = double(qcd0.sigma->getVal());
if( iDoMu) lResults[6] = double(qcd1.sigma->getVal());
lResults[7] = 0.;
if(!iDoMu) lResults[7] = qcd1.a1->getVal();
lResults[8] = nSig .getPropagatedError(*fitRes);
lResults[9] = nEWK .getPropagatedError(*fitRes);
lResults[10] = nQCD .getPropagatedError(*fitRes);
lResults[11] = nAntiSig.getPropagatedError(*fitRes);
lResults[12] = nAntiEWK.getPropagatedError(*fitRes);
lResults[13] = nAntiQCD.getPropagatedError(*fitRes);
if( iDoMu) lResults[14] = qcd0.sigma->getError();
//.........这里部分代码省略.........
示例8: RooArgList
void Fitter::finalModel(){
// -- Save final model
// final model uses splines to get the values of the multigaussians
// it will also include nuisances variables if necessary
//
for( int cat=0;cat<int(inputMasks.size()) ;++cat)
{
RooArgList *coeffs = new RooArgList();
RooArgList *pdfs = new RooArgList();
if(nBernstein >0 )
{
RooArgList *pars=new RooArgList();
for(int b=0;b<nBernstein; b++)
{
RooAbsReal *p = splines_ [ Form("spline_cat%d_par%d_bern",cat,b) ] ;
if (p == NULL ) cout<<"[Fitter]::[finalModel]::[ERROR] bern parameter"<<b<<"for cat"<<cat<<"is NULL "<<endl;
pars -> add( *p );
}
RooAbsPdf* pdf;
if (nBernstein==1) pdf=new RooBernsteinFast<1>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
if (nBernstein==2) pdf=new RooBernsteinFast<2>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
if (nBernstein==3) pdf=new RooBernsteinFast<3>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
if (nBernstein==4) pdf=new RooBernsteinFast<4>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
if (nBernstein==5) pdf=new RooBernsteinFast<5>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
if (nBernstein==6) pdf=new RooBernsteinFast<6>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
if (nBernstein==7) pdf=new RooBernsteinFast<7>(Form("bern_%d",nBernstein), Form("bern_%d",nBernstein), *x_, *pars);
pdfs ->add( *pdf );
if (nGaussians >0 ) coeffs -> add( *splines_[Form("spline_cat%d_bern_frac",cat)] );
}
for( int i=0;i< nGaussians ;++i)
{
RooAbsReal *m = splines_[ Form("spline_cat%d_g%d_mean",cat,i) ] ;
RooAbsReal *s = splines_[ Form("spline_cat%d_g%d_sigma",cat,i) ] ;
if( m== NULL ) cout<<"[Fitter]::[finalModel]::[ERROR] mean spline is NULL"<< Form("spline_cat%d_g%d_mean",cat,i)<<endl;
if( s== NULL ) cout<<"[Fitter]::[finalModel]::[ERROR] sigma spline is NULL"<< Form("spline_cat%d_g%d_sigma",cat,i)<<endl;
if (m== NULL or s==NULL) for(auto& p: splines_) cout<<"[Fitter]::[finalModel]::[DEBUG] splines_["<<p.first<<"] = "<<p.second<<endl;
RooGaussian *gaus = new RooGaussian(Form("gaus_g%d",i),Form("gaus_g%d",i),*x_, *m, *s);
pdfs -> add( *gaus );
if (i != nGaussians-1) // sum to 1
{
RooAbsReal *f = splines_ [ Form("spline_cat%d_g%d_frac",cat,i) ];
coeffs -> add( *f ) ;
}
}
string name = Form("pdf_sigmodel_cat%d",cat);
// ------------------------------------------------------------------------------\/recursive
RooAddPdf *sigModel = new RooAddPdf(name.c_str(),name.c_str(),*pdfs,*coeffs,true);
w_ -> import( *sigModel );
// -- Plot
if (plot ) {
TCanvas *c = new TCanvas();
RooPlot *p = x_ -> frame();
for(float mh=200;mh<1000;mh+= 50)
{
mh_ -> setVal(mh);
sigModel->plotOn(p);
}
p -> Draw();
c -> SaveAs( Form("%s/interpolation_cat%d.pdf",plotDir.c_str(),cat) );
c -> SaveAs( Form("%s/interpolation_cat%d.png",plotDir.c_str(),cat) );
delete p;
delete c;
} // end plot
} //end cat
}示例9: drawPlot
//.........这里部分代码省略.........
&& fabs( zeeTree->fEle2Eta) < 2.5 )) continue;
//*************************************************************************
//pt bins and eta bins
//*************************************************************************
Int_t Ele1PtBin = -1;
Int_t Ele1EtaBin = -1;
Int_t Ele2PtBin = -1;
Int_t Ele2EtaBin = -1;
if (ele1pt > 10 && ele1pt < 20) Ele1PtBin = 0;
else if (ele1pt < 30) Ele1PtBin = 1;
else if (ele1pt < 40) Ele1PtBin = 2;
else Ele1PtBin = 3;
if (ele2pt > 10 && ele2pt < 20) Ele2PtBin = 0;
else if (ele2pt < 30) Ele2PtBin = 1;
else if (ele2pt < 40) Ele2PtBin = 2;
else Ele2PtBin = 3;
if (fabs(zeeTree->fEle1SCEta) < 1.0) Ele1EtaBin = 0;
else if (fabs(zeeTree->fEle1SCEta) < 1.479) Ele1EtaBin = 1;
else Ele1EtaBin = 2;
if (fabs(zeeTree->fEle2SCEta) < 1.0) Ele2EtaBin = 0;
else if (fabs(zeeTree->fEle2SCEta) < 1.479) Ele2EtaBin = 1;
else Ele2EtaBin = 2;
if (CategoryBin == 0) {
if (!(Ele1EtaBin == 0 && Ele2EtaBin == 0)) continue;
}
else if (CategoryBin == 1) {
if (!(Ele1EtaBin == 1 && Ele2EtaBin == 1)) continue;
}
else if (CategoryBin == 2) {
if (!(Ele1EtaBin == 2 && Ele2EtaBin == 2)) continue;
}
//*************************************************************************
// restrict range of mass
//*************************************************************************
double zMass = (ele1FourVector+ele2FourVector).M();
if (zMass < minMass || zMass > maxMass) continue;
//*************************************************************************
//set mass variable
//*************************************************************************
zMassArgSet.setRealValue("mass", zMass);
data->add(zMassArgSet);
}
cout << "dataset size: " << data->numEntries() << endl;
RooRealVar *cbBias = (RooRealVar*)w->var("#Deltam_{CB}");
RooRealVar *cbSigma = (RooRealVar*)w->var("sigma_{CB}");
RooRealVar *cbCut = (RooRealVar*)w->var("a_{CB}");
RooRealVar *cbPower = (RooRealVar*)w->var("n_{CB}");
// // Now if it's a restricted fit, fix values of cbCut and cbPower to MC values.
// if (isRestricted) {
// cbCut.setConstant(kTRUE);
// cbPower.setConstant(kTRUE);
// }
// Mass model for signal electrons p.d.f.
RooAddPdf *model = (RooAddPdf*)w->pdf("model");
TCanvas* c = new TCanvas("c","c", 0,0,800,600);
//========================== Plotting ============================
//Create a frame
RooPlot* plot = mass->frame(Range(minMass,maxMass),Bins(nbins));
// Add data and model to canvas
plot->SetTitle("");
plot->GetYaxis()->SetTitleOffset(1.4);
plot->GetXaxis()->SetTitle("m_{ee} (GeV/c^{2})");
data->plotOn(plot);
model->plotOn(plot);
// model->paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(*cbBias, *cbSigma, *cbCut, *cbPower)), Layout(0.60,0.90,0.90));
model->paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(*cbBias, *cbSigma, *cbCut, *cbPower)), Layout(0.12,0.38,0.60));
plot->getAttText()->SetTextSize(.025);
plot->Draw();
// Print Fit Values
TLatex *tex = new TLatex();
tex->SetNDC();
tex->SetTextSize(.04);
tex->SetTextFont(2);
tex->DrawLatex(0.195,0.775, "Run 2012A/B");
tex->Draw();
// tex->SetTextSize(0.022);
// tex->DrawLatex(0.195, 0.75, "Z #rightarrow ee^{+}");
// tex->SetTextSize(0.024);
// tex->DrawLatex(0.645, 0.59, Form("BW Mean = %.2f GeV/c^{2}", bwMean.getVal()));
// tex->DrawLatex(0.645, 0.54, Form("BW #sigma = %.2f GeV/c^{2}", bwWidth.getVal()));
c->Update();
c->SaveAs((outFilename + ".gif").c_str());
}示例10: fitHisto
//.........这里部分代码省略.........
}
double showerTh[7] = {0.02, 0.5, 1.5, 2.5, 3.5, 4.5, 5.0};
TH1F* LiResult[5];
LiResult[0]= new TH1F("LiResult", "All", 6, showerTh);
LiResult[1]= new TH1F("LiResult", "AD1", 6, showerTh);
LiResult[2]= new TH1F("LiResult", "AD2", 6, showerTh);
LiResult[3]= new TH1F("LiResult", "AD3", 6, showerTh);
LiResult[4]= new TH1F("LiResult", "AD4", 6, showerTh);
double NumMuon[6]={0.};
double RateMuon[6]={0.};
double NumIbd[6]={0.};
RooRealVar x("x","x",0.001,40., "s");
RooRealVar tauLi9("tauLi9", "tauLi9", -0.257, "s");
RooRealVar tauHe8("tauHe8", "tauHe8", -0.172, "s");
RooRealVar rateMu("rateMu","rate of showermuon",-0.1,-10., 0.,"Hz");
RooRealVar Li9Frac("Li9Frac","Li9's Frac", 0.0);//R
RooRealVar N98("N98","total number of Li9 and He8",500.,0.,1.e5);
//RooRealVar N98("N98","total number of Li9 and He8",3e2, 1e1, 2.e5);
RooFormulaVar lambdaLi9("lambdaLi9","lambdaLi9","1/@0 + @1",RooArgList(tauLi9, rateMu));
RooFormulaVar lambdaHe8("lambdaHe8","lambdaHe8","1/@0 + @1",RooArgList(tauHe8, rateMu));
RooFormulaVar coeLi9("coeLi9","coeLi9","@0 * @1 ",RooArgList(N98,Li9Frac));
RooFormulaVar coeHe8("coeHe8","coeHe8","@0 * ( 1 - @1 )",RooArgList(N98,Li9Frac));
RooRealVar NIbd("NIbd","number of background",3e2, 1e1, 2.e6);
RooExponential expLi9("expLi9", "Li9 distribution", x, lambdaLi9);
RooExponential expHe8("expHe8", "He8 distribution", x, lambdaHe8);
double Ibdrate[3]={0.0076,0.0067,0.00082};
RooRealVar rateIbd("rateIbd","rateIbd",-Ibdrate[siteNum-1],-10,0.,"Hz");
RooFormulaVar lambdaIbd("lambdaIbd","lambdaIbd","@0 + @1",RooArgList(rateIbd, rateMu));
RooExponential expIbd("expIbd","Ibd distribution",x,rateMu);
//RooExponential expIbd("expIbd","Ibd distribution",x,lambdaIbd);
RooFitResult* fitres;
RooAddPdf* sum;
RooDataHist* data;
RooPlot* mesframe[6];
double n98[41]={0.};
double in98[41]={0.};
double nIbd[41]={0.};
double inIbd[41]={0.};
double minNLL[41]={0.};
double rateMuValue[41]={0.};
double rateMuErr[41]={0.};
//double minNl[6]={0.};
int minindex[6]={0};
double n98fit[6]={0.};
double in98fit[6]={0.};
double rateMufit[6]={0.};
double binwidth=0.;
TH1F* lh[5][6];
TString lname[5]={"","AD1","AD2","AD3","AD4"};
int lcolor[5]={4,3,2,6,5};
bool draw6slice=false;
TCanvas* c;
TString xTitle[6]={"slice 1 : 0.02~0.5GeV","slice 2 : 0.5~1.5GeV","slice 3 : 1.5~2.5GeV","slice 4 : 2.5~3.5GeV","slice 5 : 3.5~4.5GeV","slice 6 : >4.5GeV"};
if( draw6slice )
{
c = new TCanvas("c","c",2000,800) ;
c->Divide(3,2);
gStyle->SetEndErrorSize(5.0);
gStyle->SetMarkerSize(0.1);
//gStyle->SetHistLineWidth(1);
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]);示例11: sort
void Fitter::fitSignal(){
// sort mass points
sort( mIn.begin(), mIn.end() );
// -- Construct Model
RooArgList *pdfs = new RooArgList();
RooArgList *coeffs = new RooArgList();
//bernstein
addBernFitModel(pdfs,coeffs,nGaussians<=0);
// gaussians
addGausFitModel(pdfs,coeffs,true);
// construct Fitmodel
string name = "model";
RooAddPdf *fitModel = new RooAddPdf(name.c_str(),name.c_str(),*pdfs,*coeffs,true);
// -- Fit Model to data
for(int cat=0; cat < int(inputMasks.size()) ;++cat)
{
for( auto & m: mIn )
{
// ------------------------- FIT ----------------
string mass = Form( massMask_.c_str(), m );
string name = Form("cat_%d_mass_%s",cat,mass.c_str());
// switch on and off verbosity of roofit/minuit
int errlevel = -1;
int printlevel = -1;
int warnlevel = 0;
if (verbose)
{
errlevel = 1;
printlevel = 1;
warnlevel = 1;
}
// -- fit
fitModel->fitTo( *hist_[ name ] ,SumW2Error(kTRUE), PrintEvalErrors(errlevel),PrintLevel(printlevel),Warnings(warnlevel) );
fitModel->fitTo( *hist_[ name ] ,SumW2Error(kTRUE), PrintEvalErrors(errlevel),PrintLevel(printlevel),Warnings(warnlevel) );
// -- Plot
if (plot ) {
TCanvas *c = new TCanvas(Form("c_%s_cat%d",mass.c_str(),cat),"c");
RooPlot *p = x_ -> frame();
hist_[ name ] -> plotOn(p,DataError(RooAbsData::SumW2));
fitModel -> plotOn(p);
for(int i=0;i < nGaussians ;++i)
{
fitModel -> plotOn(p, Components(Form("fitmodel_gaus_g%d",i)),LineStyle(kDashed) );
}
if( nBernstein >0 )
{
fitModel -> plotOn(p, Components(Form("fitmodel_bern_%d",nBernstein)),LineStyle(kDotted) );
}
//c->SetLogy();
p -> Draw();
c -> Update() ;
c -> SaveAs( Form("%s/fit_mh_%s_cat%d.pdf",plotDir.c_str(),mass.c_str(),cat) );
c -> SaveAs( Form("%s/fit_mh_%s_cat%d.png",plotDir.c_str(),mass.c_str(),cat) );
//
delete p;
delete c;
} // end plots
// -- Save coefficients vs mh
saveCoefficientsBern(cat,mass, nGaussians<=0);
saveCoefficientsGaus(cat,mass,true);
} // end mass
// -- Interpolate coefficients
// -- bern
interpolateBern(cat,nGaussians<=0);
// -- gaus
interpolateGaus(cat,true);
}
return;
}示例12: RooProdPdf
// pdf //
/////////
RooProdPdf *pdf_rho = new RooProdPdf("pdf_rho","pdf_rho",pdf_de_rho,Conditional(pdf_mbc_rho,mbc));
// RooProdPdf pdf_rho("pdf_rho","pdf_rho",pdf_de_rho,pdf_mbc_rho);
}
//////////////////
// Complete PDF //
//////////////////
// RooRealVar fsig("fsig","fsig",0.1,0.,1.);// fsig.setConstant(kTRUE);
// RooRealVar frho("frho","frho",0.1,0.,1.);// frho.setConstant(kTRUE);
RooRealVar Nsig("Nsig","Nsig",1150,100.,10000.);
RooRealVar Nrho("Nrho","Nrho",9308,100,100000.);
RooRealVar Ncmb("Ncmb","Ncmb",21288,1000,100000);
// RooAddPdf pdf("pdf","pdf",RooArgList(pdf_sig,pdf_rho,pdf_comb),RooArgList(fsig,frho));
RooAddPdf pdf("pdf","pdf",RooArgList(pdf_sig,*pdf_rho,pdf_comb),RooArgList(Nsig,Nrho,Ncmb));
RooArgSet* params = pdf.getParameters(RooArgSet(de,mbc));
// RooArgset* initParams = (RooArgSet*) params->snapshot();
pdf.fitTo(ds,Verbose(),Timer(true));
params->printLatex(OutputFile("PurityFit.tex"));
RooAbsReal* intSig = pdf_sig.createIntegral(RooArgSet(de,mbc),NormSet(RooArgSet(de,mbc)),Range("Signal"));
RooAbsReal* intRho = pdf_rho->createIntegral(RooArgSet(de,mbc),NormSet(RooArgSet(de,mbc)),Range("Signal"));
RooAbsReal* intCmb = pdf_comb.createIntegral(RooArgSet(de,mbc),NormSet(RooArgSet(de,mbc)),Range("Signal"));
const double nsig = intSig->getVal()*Nsig.getVal();
const double nsig_err = intSig->getVal()*Nsig.getError();
const double nsig_err_npq = TMath::Sqrt(nsig*(Nsig.getVal()-nsig)/Nsig.getVal());
const double nsig_err_total = TMath::Sqrt(nsig_err*nsig_err+nsig_err_npq*nsig_err_npq);示例13: mass
RooRealVar mass("mass","The invariant mass",100.,150.);
//Signal PDF
RooRealVar mean1("mean1","Mean of first signal gaussian",125.,110.,140.);
RooRealVar sigma1("sigma1","Sigma of first signal gaussian",1.,0.001,4.);
RooGaussian SigGauss1("SigGauss1","First signal gaussian",mass,mean1,sigma1);
//Background PDF
RooRealVar a0("a0","a0",-1.,-3.,3.);
RooRealVar a1("a1","a1",0.5,-3.,3.);
RooChebychev BkgPDF("BkgPDF","BkgPDF",mass,RooArgSet(a0,a1));
//Total first PDF
RooRealVar frac1("frac1","Fraction for first PDF",0.6,0.,1.);
RooAddPdf totPDF1("totPDF1","Total PDF 1",RooArgList(SigGauss1,BkgPDF),RooArgList(frac1));
//Build the second signal PDF
RooRealVar sigma2("sigma2","Sigma of second signal gaussian",2.,0.001,4.);
RooGaussian SigGauss2("SigGauss2","Second signal gaussian",mass,mean1,sigma2);
//Total second PDF
RooRealVar frac2("frac2","Fraction for second PDF",0.4,0.,1.);
RooAddPdf totPDF2("totPDF2","Total PDF 2",RooArgList(SigGauss2,BkgPDF),RooArgList(frac2));
//Generate the two samples
RooDataSet *data1 = totPDF1.generate(RooArgSet(mass),1000);
RooDataSet *data2 = totPDF2.generate(RooArgSet(mass),2000);
RooCategory SigCat("SigCat","Signal categories");
SigCat.defineType("Signal1");示例14: runQuickRejig
string runQuickRejig(string oldfilename, int ncats){
string newfilename="TestWS.root";
TFile *oldFile = TFile::Open(oldfilename.c_str());
TFile *newFile = new TFile(newfilename.c_str(),"RECREATE");
RooWorkspace *oldWS = (RooWorkspace*)oldFile->Get("cms_hgg_workspace");
RooWorkspace *newWS = new RooWorkspace("wsig_8TeV");
RooRealVar *mass = (RooRealVar*)oldWS->var("CMS_hgg_mass");
RooDataSet *tot;
RooRealVar *normT1 = new RooRealVar("nT1","nT1",500,0,550);
RooRealVar *meanT1 = new RooRealVar("mT1","mT1",125,122,128);
RooRealVar *sigmaT1 = new RooRealVar("sT1","sT1",3.1,1.,10.);
RooGaussian *gausT1 = new RooGaussian("gT1","gT1",*mass,*meanT1,*sigmaT1);
RooRealVar *normT2 = new RooRealVar("nT2","nT2",2,0,500);
RooRealVar *meanT2 = new RooRealVar("mT2","mT2",125,122,128);
RooRealVar *sigmaT2 = new RooRealVar("sT2","sT2",1.7,0.,10.);
RooGaussian *gausT2 = new RooGaussian("gT2","gT2",*mass,*meanT2,*sigmaT2);
RooRealVar *normT3 = new RooRealVar("nT3","nT3",2,0,500);
RooRealVar *meanT3 = new RooRealVar("mT3","mT3",125,122,128);
RooRealVar *sigmaT3 = new RooRealVar("sT3","sT3",1.7,0.,10.);
RooGaussian *gausT3 = new RooGaussian("gT3","gT3",*mass,*meanT3,*sigmaT3);
RooAddPdf *gausT = new RooAddPdf("gT","gT",RooArgList(*gausT1,*gausT2,*gausT3),RooArgList(*normT1,*normT2,*normT3));
for (int cat=0; cat<ncats; cat++){
RooDataSet *thisData = (RooDataSet*)oldWS->data(Form("sig_ggh_mass_m125_cat%d",cat));
newWS->import(*thisData);
RooRealVar *norm1 = new RooRealVar(Form("n1%d",cat),Form("n1%d",cat),500,0,550);
RooRealVar *mean1 = new RooRealVar(Form("m1%d",cat),Form("m1%d",cat),125,122,128);
RooRealVar *sigma1 = new RooRealVar(Form("s1%d",cat),Form("s1%d",cat),3.1,1.,10.);
RooGaussian *gaus1 = new RooGaussian(Form("g1%d",cat),Form("g1%d",cat),*mass,*mean1,*sigma1);
RooRealVar *norm2 = new RooRealVar(Form("n2%d",cat),Form("n2%d",cat),2,0,500);
RooRealVar *mean2 = new RooRealVar(Form("m2%d",cat),Form("m2%d",cat),125,122,128);
RooRealVar *sigma2 = new RooRealVar(Form("s2%d",cat),Form("s2%d",cat),1.7,0.,10.);
RooGaussian *gaus2 = new RooGaussian(Form("g2%d",cat),Form("g2%d",cat),*mass,*mean2,*sigma2);
RooRealVar *norm3 = new RooRealVar(Form("n3%d",cat),Form("n3%d",cat),2,0,500);
RooRealVar *mean3 = new RooRealVar(Form("m3%d",cat),Form("m3%d",cat),125,122,128);
RooRealVar *sigma3 = new RooRealVar(Form("s3%d",cat),Form("s3%d",cat),1.7,0.,10.);
RooGaussian *gaus3 = new RooGaussian(Form("g3%d",cat),Form("g3%d",cat),*mass,*mean3,*sigma3);
RooAddPdf *gaus = new RooAddPdf(Form("g%d",cat),"g",RooArgList(*gaus1,*gaus2,*gaus3),RooArgList(*norm1,*norm2,*norm3));
gaus->fitTo(*thisData,SumW2Error(kTRUE));
newWS->import(*gaus);
if (cat==0) {
tot = thisData;
tot->SetName("sig_ggh_m125");
}
else tot->append(*thisData);
}
newWS->import(*tot);
gausT->fitTo(*tot,SumW2Error(kTRUE));
newWS->import(*gausT);
newWS->Write();
oldFile->Close();
newFile->Close();
delete newFile;
delete newWS;
return newfilename;
}示例15: PDF
PDF(TString chrg, double scale, bool addqcd, RooRealVar* x, RooRealVar* f1, RooRealVar *f2, bool redice){
TFile *data = new TFile("eWPol_Signal_Wjets.root");
TFile *fbkg= new TFile("eWPol_Signal_QCD.root");
TH1D *hsig, *hbkg;
if (chrg=="plus") {
mc1 = (TH1D*)data->Get(pHist1); mc2 = (TH1D*)data->Get(pHist2); mc3 = (TH1D*)data->Get(pHist3);
hsig = (TH1D*)data->Get(pHist_data); hsig->Rebin(rbin);
hbkg = (TH1D*)fbkg->Get(pHist_data); hbkg->Rebin(rbin);
} else {
mc1 = (TH1D*)data->Get(mHist1); mc2 = (TH1D*)data->Get(mHist2); mc3 = (TH1D*)data->Get(mHist3);
hsig = (TH1D*)data->Get(mHist_data); hsig->Rebin(rbin);
hbkg = (TH1D*)fbkg->Get(mHist_data); hbkg->Rebin(rbin);
}
mc1->Rebin(rbin); mc2->Rebin(rbin); mc3->Rebin(rbin);
TH1D *test = hsig->Clone();
// QCD background
Double_t nbkg=0;
ScaleErrors(hbkg,scale);
ScaleErrors(test,scale);
ScaleErrors(hsig,scale);
ScaleErrors(mc1,scale); ScaleErrors(mc2,scale); ScaleErrors(mc3,scale);
if(addqcd){
test->Add(hbkg); nbkg=hbkg->Integral();
}
stat=test->Integral();
Double_t f_sig=(hsig->Integral())/stat;
Double_t f_bkg=nbkg/stat;
// Import binned Data
data1 = new RooDataHist ("data1","dataset with ICvar",*x,mc1);
data2 = new RooDataHist ("data2","dataset with ICvar",*x,mc2);
data3 = new RooDataHist ("data3","dataset with ICvar",*x,mc3);
// Background template
data4 = new RooDataHist("data4","dataset with ICVar",*x,hbkg);
bkg = new RooHistPdf("bkg","bkg",*x,*data4);
// Helicity fractions
//3 = new RooRealVar("f_{0}","f3 fraction",100.,0.,100000.) ;
f3 = new RooRealVar("f_{0}","f3 fraction",0.3,0.,1.) ;
fs = new RooRealVar("f_{s}","f4 fraction",f_sig);//,0.,1.) ;
// Templates
h1 = new RooHistPdf("h1","h1",*x,*data1);
h2 = new RooHistPdf("h2","h2",*x,*data2);
h3 = new RooHistPdf("h3","h3",*x,*data3);
// Model which changes the interpretation of parameters
RooAddPdf *com = new RooAddPdf("cplus1","component 1",RooArgList(*h2,*h3),RooArgList(*f2)) ;
sig = new RooAddPdf("sig","component 1",RooArgList(*h1,*com),RooArgList(*f1)) ;
// Original Model (allows negative values of f0...)
//sig = new RooAddPdf("sig","component 1",RooArgList(*h1,*h2,*h3),RooArgList(*f1,*f2)) ;
// composite PDF
model = new RooAddPdf("model","model",RooArgList(*sig,*bkg),*fs);
// if (redice) {
test1=model->generate(*x,int(stat));
//} else {
rtest1 = new RooDataHist("data","dataset with LPvar",*x,test);
//}
}