本文整理汇总了C++中RooPlot类的典型用法代码示例。如果您正苦于以下问题:C++ RooPlot类的具体用法?C++ RooPlot怎么用?C++ RooPlot使用的例子?那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了RooPlot类的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: getFitResults
FitStruct getFitResults( const std::string& outputdir, TTree* tree, const std::string& name, float xMin, float xMax, std::string& cut, std::string& whatToProject ) {
DrawTools::setStyle();
std::string histoName(Form("h1_%s", name.c_str() ));
std::string histoName2(Form("h2_%s", name.c_str() ));
TH1D* h1;
TF1* f1;
gStyle->SetOptFit(1);
h1 = new TH1D(histoName.c_str(), "", 5000, xMin, xMax);
tree->Project( histoName.c_str(), whatToProject.c_str(), cut.c_str() );
f1 = new TF1( Form("gaus_%s", name.c_str()), "gaus", xMin, xMax);
h1->Fit( f1, "RQN" );
f1->SetLineColor(kRed);
double peakpos = h1->GetMean();
double sigma = h1->GetRMS();
// double peakpos = f1->GetParameter(1);
// double sigma = f1->GetParameter(2);
double fitmin;
double fitmax;
if( (peakpos-8*sigma) < 4000 ) {
fitmin = 4000;
} else {
fitmin = peakpos-8*sigma;
}
if(((peakpos>0) && (sigma > 2*peakpos)) ) {
fitmax = 2*peakpos;
} else {
fitmax = peakpos+5*sigma;
}
if (sigma < 500) sigma = 500;
TH1D* h2 = new TH1D("h2", "", 200, fitmin,fitmax);
tree->Project( "h2", whatToProject.c_str(), cut.c_str() );
RooRealVar x("x","ADC Channels", fitmin, fitmax);
RooDataHist data("data","dataset with x",x,Import(*h2) );
RooPlot* frame;
RooPlot* xframe = x.frame();
frame = x.frame("Title");
data.plotOn(frame); //this will show histogram data points on canvas
// data.statOn(frame); //this will display hist stat on canvas
RooRealVar meanr("meanr","Mean",peakpos,peakpos-sigma, peakpos+ sigma);
RooRealVar width("width","#sigma",sigma *0.3, 150.0, 5.*sigma);
RooRealVar A("A","Dist",2., 0.0, 7.0);
RooRealVar N("N","Deg",5, 0.0, 10);
meanr.setRange( 30000. , 1000000.);
width.setRange(500, 22000);
RooCBShape fit_fct("fit_fct","fit_fct",x,meanr,width,A,N);
int ndf = 4;
fit_fct.fitTo(data);
fit_fct.plotOn(frame,LineColor(4));//this will show fit overlay on canvas
// fit_fct.paramOn(frame); //this will display the fit parameters on canvas
double mean = meanr.getVal();
double meanErr = meanr.getError();
double rms = width.getVal();
double rmsErr = width.getError();
double reso = 100.* rms/mean; //in percent
double resoErr = 100.* getRatioError( rms, mean, meanErr, rmsErr );
TCanvas* cans = new TCanvas("cans", "un canvas", 600,600);
cans->cd();
frame->Draw();
TLegend* lego = new TLegend(0.6, 0.7, 0.9, 0.92);
lego->SetTextSize(0.038);
lego->AddEntry( (TObject*)0 ,Form("#mu = %.0f #pm %.0f", meanr.getVal(), meanr.getError() ), "");
lego->AddEntry( (TObject*)0 ,Form("#sigma = %.0f #pm %.0f ", width.getVal(), width.getError() ), "");
lego->AddEntry( (TObject*)0 ,Form("#chi^{2} = %.2f / %d ", frame->chiSquare(ndf) , ndf ), "");
lego->AddEntry( (TObject*)0 ,Form("#sigma/#mu = %.1f #pm %.1f %s ", reso , resoErr ,"%"), "");
lego->SetFillColor(0);
lego->Draw("same");
cans->SaveAs( Form( "%s/CBFit_%s.png", outputdir.c_str(), name.c_str() ) );
FitStruct fs;
//.........这里部分代码省略.........
示例2: fitDilMass
void fitDilMass(){
using namespace RooFit;
gROOT->SetStyle("Plain");
//declare variables
RooRealVar mll("mll", "dilepton mass", 0,100 ,"GeV");
RooRealVar edge("edge", "dilepton mass edge", 0,100 ,"GeV");
RooRealVar min("min", "min dilepton mass", 0,100 ,"GeV");
RooRealVar k("k", "constant", 0,100 ,"GeV");
RooRealVar mean("mean", "mean of mll resolution", 0,100 ,"GeV");
RooRealVar sigma("sigma", "sigma of mll resolution", 0,10 ,"GeV");
//set variables
mll.setBins(100);
edge.setVal(50.);
mean.setVal(0.);
sigma.setVal(2.);
min.setVal(10.);
min.setConstant();
k.setVal(1.);
mean.setConstant();
sigma.setConstant();
mll.setMin(0);
mll.setMax(100);
//get data
TFile* file = TFile::Open("LM0_dilmass.root");
TH1D* hist = static_cast<TH1D *>((file)->Get("dilmass"));
RooDataHist data("data","data histogram",RooArgSet(mll),hist);
//declare and make PDFs
RooGenericPdf tri("tri","( (mll < min && mll > 0) * k + (mll > min && mll < edge ) * (mll - min) + (mll > edge && mll < 200) * k)", RooArgList(mll , edge, min, k) );
RooGaussian gaus("gaus","Detector Gaussian",mll,mean,sigma);
RooNumConvPdf trigaus("trigaus","triangle + gaussian function",mll,tri,gaus);
trigaus.setConvolutionWindow(mean,sigma,25);
//do fit
RooFitResult *result = trigaus.fitTo( data , Save() , Minos(kFALSE) , SumW2Error(kFALSE));
//make plot
TCanvas *c1=new TCanvas("c1","",800,600);
c1->cd();
RooPlot* frame = mll.frame();
frame->SetXTitle("dilepton mass (GeV)");
frame->SetYTitle("");
frame->SetTitle("");
data.plotOn(frame);
tri.plotOn(frame , LineStyle(kDashed));
trigaus.plotOn( frame );
trigaus.paramOn( frame );
frame->Draw();
//frame->SetMinimum(-10);
//frame->SetMaximum(100);
int nfloatpars = result->floatParsFinal().getSize();
int ndf = 20 - nfloatpars;
float chi2 = frame->chiSquare(nfloatpars)*ndf;
float prob = TMath::Prob(chi2,ndf);
cout << endl << endl;
cout << "Fit Results-------------------------------------------" << endl;
cout << "chi2/ndf = " << chi2 << " / " << ndf << endl;
cout << "prob = " << prob << endl;
cout << "edge = " << edge.getVal() << " +/- " << edge.getError() << endl;
cout << "k = " << k.getVal() << " +/- " << k.getError() << endl;
cout << "------------------------------------------------------" << endl;
cout << endl << endl;
}
示例3: deFit
//.........这里部分代码省略.........
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);
const double nrho = intRho->getVal()*Nrho.getVal();
const double nrho_err = intRho->getVal()*Nrho.getError();
const double nrho_err_npq = TMath::Sqrt(nrho*(Nrho.getVal()-nrho)/Nrho.getVal());
const double nrho_err_total = TMath::Sqrt(nrho_err*nrho_err+nrho_err_npq*nrho_err_npq);
const double ncmb = intCmb->getVal()*Ncmb.getVal();
const double ncmb_err = intCmb->getVal()*Ncmb.getError();
const double ncmb_err_npq = TMath::Sqrt(ncmb*(Ncmb.getVal()-ncmb)/Ncmb.getVal());
const double ncmb_err_total = TMath::Sqrt(ncmb_err*ncmb_err+ncmb_err_npq*ncmb_err_npq);
const double purity = nsig/(nsig+nrho+ncmb);
const double purity_err = nsig_err_total/(nsig+nrho+ncmb);
double sig_frac;
double pdf_sig_val;
double pdf_DK_val;
double pdf_smooth_val;
fstream ofile("de_sig_fraction.txt",fstream::out);
for(int i=0; i<1000; i++){
const double dde = 0.2/1000;
de.setVal(-0.1+(i+0.5)*dde);
pdf_sig_val = Nsig.getVal()*pdf_sig.getVal(de);
pdf_DK_val = Nrho.getVal()*gDK.getVal(de);
pdf_smooth_val = Ncmb.getVal()*pdf_comb.getVal(de);
sig_frac = pdf_sig_val/(pdf_sig_val+pdf_DK_val+pdf_smooth_val);
ofile << de.getVal() << " " << sig_frac << endl;
}
ofile.close();
/////////////
// Plots //
/////////////
// de //
RooPlot* deFrame = de.frame();
ds.plotOn(deFrame,DataError(RooAbsData::SumW2),MarkerSize(1));
pdf.plotOn(deFrame,Components(gDK),LineStyle(kDashed));
pdf.plotOn(deFrame,Components(pdf_sig),LineStyle(kDashed));
pdf.plotOn(deFrame,Components(pdf_comb),LineStyle(kDashed));
pdf.plotOn(deFrame,LineWidth(2));
RooHist* hdepull = deFrame->pullHist();
RooPlot* dePull = de.frame(Title("#Delta E pull distribution"));
dePull->addPlotable(hdepull,"P");
dePull->GetYaxis()->SetRangeUser(-5,5);
TCanvas* cm = new TCanvas("Delta E","Delta E",600,700);
cm->cd();
TPad *pad3 = new TPad("pad3","pad3",0.01,0.20,0.99,0.99);
TPad *pad4 = new TPad("pad4","pad4",0.01,0.01,0.99,0.20);
pad3->Draw();
pad4->Draw();
pad3->cd();
pad3->SetLeftMargin(0.15);
pad3->SetFillColor(0);
pad3->SetGrid();
deFrame->GetXaxis()->SetTitleSize(0.05);
deFrame->GetXaxis()->SetTitleOffset(0.85);
deFrame->GetXaxis()->SetLabelSize(0.04);
deFrame->GetYaxis()->SetTitleOffset(1.6);
deFrame->Draw();
const int height = svd == 2 ? 500 : 120;
TLine *de_line_RIGHT = new TLine(de_max,0,de_max,height);
示例4: fitM3
//.........这里部分代码省略.........
RooMCStudy *mc_wjets = new RooMCStudy(ext_hpdf_wjets,mass,Binned(kTRUE),Silence(kTRUE));
mc_wjets->generate(Nsamples,0,kFALSE,"data/toymc_wjets_%04d.dat");
RooMCStudy *mc_zjets = new RooMCStudy(ext_hpdf_zjets,mass,Binned(kTRUE),Silence(kTRUE));
mc_zjets->generate(Nsamples,0,kFALSE,"data/toymc_zjets_%04d.dat");
RooMCStudy *mc_qcd = new RooMCStudy(ext_hpdf_qcd,mass,Binned(kTRUE),Silence(kTRUE));
mc_qcd->generate(Nsamples,0,kFALSE,"data/toymc_qcd_%04d.dat");
RooMCStudy *mc_ST = new RooMCStudy(ext_hpdf_ST,mass,Binned(kTRUE),Silence(kTRUE),FitOptions(ExternalConstraints(STgaussConstraint)));
mc_ST->generate(Nsamples,0,kFALSE,"data/toymc_ST_%04d.dat");
return;
*/
RooMCStudy *mcstudy = new RooMCStudy(model_M3, mass, FitModel(model_histpdf),Binned(kTRUE),Silence(kTRUE), Extended() ,
//FitOptions(Save(kTRUE),Minos(kTRUE),Extended()) );
FitOptions(Save(kTRUE),Minos(kTRUE),Extended(),ExternalConstraints(STgaussConstraint)));//RooArgList(fconstext,fconstext2)) )); //gaussian constraint
//mcstudyM3->generate(Nsamples,0,kFALSE,"toymc.dat");
//mcstudyM3->generateAndFit(Nsamples,0,kFALSE,"toymc.dat");
//TList dataList;
//for (int isample=0; isample<Nsamples; ++isample) dataList.Add( mcstudyM3->genData(isample));
// Fit
mcstudy->generateAndFit(Nsamples,0,kTRUE);
//mcstudy->fit(Nsamples, "data/toymc_%04d.dat");
gDirectory->Add(mcstudy) ;
// E x p l o r e r e s u l t s o f s t u d y
// ------------------------------------------------
// Make plots of the distributions of mean, the error on mean and the pull of mean
RooPlot* frame1 = mcstudy->plotParam(Ntt,Bins(40));
RooPlot* frame2 = mcstudy->plotError(Ntt,Bins(40)) ;
RooPlot* frame3 = mcstudy->plotPull(Ntt,Bins(40),FitGauss(kTRUE)) ;
RooPlot* frame1w = mcstudy->plotParam(Nbkg,Bins(40)) ;
RooPlot* frame2w = mcstudy->plotError(Nbkg,Bins(40)) ;
RooPlot* frame3w = mcstudy->plotPull(Nbkg,Bins(40),FitGauss(kTRUE)) ;
RooPlot* frame1st = mcstudy->plotParam(NST,Bins(40)) ;
RooPlot* frame2st = mcstudy->plotError(NST,Bins(40)) ;
//RooPlot* frame3st = mcstudy->plotPull(NST,Bins(40),FitGauss(kTRUE)) ;
// Plot distribution of minimized likelihood
RooPlot* frame4 = mcstudy->plotNLL(Bins(40)) ;
// Make some histograms from the parameter dataset
TH1* hh_cor_ttbar_w = mcstudy->fitParDataSet().createHistogram("hh",Ntt,YVar(Nbkg)) ;
// Access some of the saved fit results from individual toys
//TH2* corrHist000 = mcstudy->fitResult(0)->correlationHist("c000") ;
//TH2* corrHist127 = mcstudy->fitResult(127)->correlationHist("c127") ;
//TH2* corrHist953 = mcstudy->fitResult(953)->correlationHist("c953") ;
// Draw all plots on a canvas
gStyle->SetPalette(1) ;
gStyle->SetOptStat(0) ;
TCanvas* cv = new TCanvas("cv","cv",600,600) ;
hM3->SetFillColor(kRed);
hWjets->SetFillColor(kGreen);
hM3->Draw();
hWjets->Draw("same");
gPad->RedrawAxis();
示例5: TCanvas
void MakeSpinPlots::DrawBlindFit(TString tag, TString mcName,TString cosThetaBin){
TString fitTag="FULLFIT";
TString cat = "evtcat";
if(cosThetaBin!=""){
tag+="_"+cosThetaBin;
fitTag="FULLCOSTFIT";
cat = "evtcat_cosT";
}
TString dataTag = "_Combined";
if(cosThetaBin!="") dataTag+="_CosTBin";
TCanvas *cv = new TCanvas(Form("%s_%s",mcName.Data(),tag.Data()));
RooRealVar* mass = ws->var("mass");
mass->setBins( (mass->getMax() - mass->getMin())/1.5 ); //enfore 1.5GeV bin width
RooPlot* frame = mass->frame();
double Nb = ws->var(Form("Data_BKGFIT_%s_Nbkg",tag.Data()))->getVal();
cout << Nb << endl;
RooDataSet *blind = (RooDataSet*)ws->data("Data"+dataTag)->reduce(TString("((mass<119) || (mass>135.5)) && ")+cat+"=="+cat+"::"+tag);
blind->plotOn(frame);
tPair lbl(mcName,tag);
double nBkg = ws->data("Data"+dataTag)->sumEntries(cat+"=="+cat+"::"+tag);
ws->pdf( Form("Data_BKGFIT_%s_bkgModel",tag.Data()) )->plotOn(frame,RooFit::Range("all"),RooFit::Normalization(nBkg/blind->sumEntries()),
RooFit::LineColor(kRed));
//TLatex *prelim = new TLatex(250,x->GetXmax()-40.,"CMS Preliminary");
TLatex *prelim = new TLatex(0.12,0.96,"CMS Preliminary");
TLatex *lum = new TLatex(0.7,0.96,Form("#sqrt{s}=8 TeV L = %0.1f fb^{-1}",lumi));
prelim->SetNDC();
lum->SetNDC();
prelim->SetTextSize(0.045);
prelim->SetTextColor(kBlack);
lum->SetTextSize(0.045);
lum->SetTextColor(kBlack);
TLatex *owner = new TLatex(0.6,0.88,"Caltech-CMS Preliminary");
owner->SetNDC();
owner->SetTextSize(0.045);
owner->SetTextColor(kBlack);
TLatex *Nbkg = new TLatex(0.7,0.8,Form("N_{bkg}= %0.1f #pm %0.1f",nBackground[lbl].first,nBackground[lbl].second));
Nbkg->SetNDC();
Nbkg->SetTextSize(0.045);
TLatex *sig = new TLatex(0.7,0.72,Form("#sigma_{eff} = %0.1f #pm %0.2f",fitSigEff[lbl].first,fitSigEff[lbl].second));
sig->SetNDC();
sig->SetTextSize(0.045);
TLatex *expBkg = new TLatex(0.7,0.64,Form("B @ 125 = %0.1f",fitBkg1Sigma[lbl].first));
expBkg->SetNDC();
expBkg->SetTextSize(0.045);
frame->addObject(prelim);
frame->addObject(lum);
//frame->addObject(owner);
frame->addObject(Nbkg);
frame->addObject(sig);
frame->addObject(expBkg);
frame->Draw();
cv->SaveAs( basePath+Form("/mgg-%s-%s-%s_BLIND.png",outputTag.Data(),mcName.Data(),tag.Data()) );
cv->SaveAs( basePath+Form("/C/mgg-%s-%s-%s_BLIND.C",outputTag.Data(),mcName.Data(),tag.Data()) );
cv->SaveAs( basePath+Form("/mgg-%s-%s-%s_BLIND.pdf",outputTag.Data(),mcName.Data(),tag.Data()) );
delete cv;
}
示例6: drawCtauFrom2DPlot
void drawCtauFrom2DPlot(RooWorkspace& myws, // Local workspace
string outputDir, // Output directory
struct InputOpt opt, // Variable with run information (kept for legacy purpose)
struct KinCuts cut, // Variable with current kinematic cuts
map<string, string> parIni, // Variable containing all initial parameters
string plotLabel, // The label used to define the output file name
// Select the type of datasets to fit
string DSTAG, // Specifies the type of datasets: i.e, DATA, MCJPSINP, ...
bool isPbPb, // Define if it is PbPb (True) or PP (False)
// Select the type of object to fit
bool incJpsi, // Includes Jpsi model
bool incPsi2S, // Includes Psi(2S) model
bool incBkg, // Includes Background model
// Select the fitting options
// Select the drawing options
bool setLogScale, // Draw plot with log scale
bool incSS, // Include Same Sign data
double binWidth // Bin width
)
{
RooMsgService::instance().getStream(0).removeTopic(Caching);
RooMsgService::instance().getStream(1).removeTopic(Caching);
RooMsgService::instance().getStream(0).removeTopic(Plotting);
RooMsgService::instance().getStream(1).removeTopic(Plotting);
RooMsgService::instance().getStream(0).removeTopic(Integration);
RooMsgService::instance().getStream(1).removeTopic(Integration);
RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ;
if (DSTAG.find("_")!=std::string::npos) DSTAG.erase(DSTAG.find("_"));
string pdfTotName = Form("pdfCTAUMASS_Tot_%s", (isPbPb?"PbPb":"PP"));
string dsOSName = Form("dOS_%s_%s", DSTAG.c_str(), (isPbPb?"PbPb":"PP"));
string dsOSNameCut = dsOSName+"_CTAUCUT";
string hOSName = Form("dhCTAUERRTot_Tot_%s", (isPbPb?"PbPb":"PP"));
string hOSNameBkg = Form("dhCTAUERR_Bkg_%s", (isPbPb?"PbPb":"PP"));
string hOSNameJpsi = Form("dhCTAUERR_Jpsi_%s", (isPbPb?"PbPb":"PP"));
string hOSNamePsi2S = Form("dhCTAUERR_Psi2S_%s", (isPbPb?"PbPb":"PP"));
string dsSSName = Form("dSS_%s_%s", DSTAG.c_str(), (isPbPb?"PbPb":"PP"));
bool isWeighted = myws.data(dsOSName.c_str())->isWeighted();
vector<double> range; range.push_back(cut.dMuon.ctau.Min); range.push_back(cut.dMuon.ctau.Max);
double minRange = -4.0;
double maxRange = 7.0;
Double_t outTot = myws.data(dsOSName.c_str())->numEntries();
Double_t outErr = myws.data(dsOSName.c_str())->reduce(Form("(ctau>%.6f || ctau<%.6f)", range[1], range[0]))->numEntries();
int nBins = min(int( round((maxRange - minRange)/binWidth) ), 1000);
double normDSTot = 1.0; if (myws.data(dsOSNameCut.c_str())) { normDSTot = myws.data(dsOSName.c_str())->sumEntries()/myws.data(dsOSNameCut.c_str())->sumEntries(); }
double normJpsi = 1.0; if (myws.data(hOSNameJpsi.c_str())) { normJpsi = myws.data(dsOSName.c_str())->sumEntries()*normDSTot/myws.data(hOSNameJpsi.c_str())->sumEntries(); }
double normPsi2S = 1.0; if (myws.data(hOSNamePsi2S.c_str())) { normPsi2S = myws.data(dsOSName.c_str())->sumEntries()*normDSTot/myws.data(hOSNamePsi2S.c_str())->sumEntries(); }
double normBkg = 1.0; if (myws.data(hOSNameBkg.c_str())) { normBkg = myws.data(dsOSName.c_str())->sumEntries()*normDSTot/myws.data(hOSNameBkg.c_str())->sumEntries(); }
double normTot = 1.0; if (myws.data(hOSName.c_str())) { normTot = myws.data(dsOSName.c_str())->sumEntries()*normDSTot/myws.data(hOSName.c_str())->sumEntries(); }
// Create the main plot of the fit
RooPlot* frame = myws.var("ctau")->frame(Bins(nBins), Range(minRange, maxRange));
frame->updateNormVars(RooArgSet(*myws.var("invMass"), *myws.var("ctau"), *myws.var("ctauErr"))) ;
myws.data(dsOSName.c_str())->plotOn(frame, Name("dOS"), DataError(RooAbsData::SumW2), XErrorSize(0), MarkerColor(kBlack), LineColor(kBlack), MarkerSize(1.2));
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("PDF"),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSNameCut.c_str()), kTRUE),
Normalization(normDSTot, RooAbsReal::NumEvent),
FillStyle(1001), FillColor(kViolet+6), VLines(), DrawOption("LF"), NumCPU(32), LineColor(kBlack)
);
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("BKG"),Components(RooArgSet( *myws.pdf(Form("pdfCTAUMASS_Bkg_%s", (isPbPb?"PbPb":"PP"))) )),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSName.c_str()), kTRUE),
Normalization(normDSTot, RooAbsReal::NumEvent),
FillStyle(1001), FillColor(kAzure-9), VLines(), DrawOption("LF"), NumCPU(32)
);
if (incJpsi) {
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("JPSIPR"),Components(RooArgSet( *myws.pdf(Form("pdfCTAUMASS_JpsiPR_%s", (isPbPb?"PbPb":"PP"))) )),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSName.c_str()), kTRUE),
Normalization(normDSTot, RooAbsReal::NumEvent),
LineColor(kRed+3), Precision(1e-5), NumCPU(32)
);
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("JPSINOPR"),Components(RooArgSet( *myws.pdf(Form("pdfCTAUMASS_JpsiNoPR_%s", (isPbPb?"PbPb":"PP"))) )),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSName.c_str()), kTRUE),
Normalization(normDSTot, RooAbsReal::NumEvent),
LineColor(kGreen+3), Precision(1e-5), NumCPU(32)
);
}
if (incPsi2S) {
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("PSI2SPR"),Components(RooArgSet( *myws.pdf(Form("pdfCTAUMASS_Psi2SPR_%s", (isPbPb?"PbPb":"PP"))) )),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSName.c_str()), kTRUE),
Normalization(normDSTot, RooAbsReal::NumEvent),
LineColor(kRed+3), Precision(1e-5), NumCPU(32)
);
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("PSI2SNOPR"),Components(RooArgSet( *myws.pdf(Form("pdfCTAUMASS_Psi2SNo_%s", (isPbPb?"PbPb":"PP"))) )),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSName.c_str()), kTRUE),
Normalization(normDSTot, RooAbsReal::NumEvent),
LineColor(kGreen+3), Precision(1e-5), NumCPU(32)
);
}
if (incSS) {
myws.data(dsSSName.c_str())->plotOn(frame, Name("dSS"), MarkerColor(kRed), LineColor(kRed), MarkerSize(1.2));
}
myws.data(dsOSName.c_str())->plotOn(frame, Name("dOS"), DataError(RooAbsData::SumW2), XErrorSize(0), MarkerColor(kBlack), LineColor(kBlack), MarkerSize(1.2));
myws.pdf(pdfTotName.c_str())->plotOn(frame,Name("PDFLINE"),
ProjWData(RooArgSet(*myws.var("ctauErr")), *myws.data(dsOSNameCut.c_str()), kTRUE),
//.........这里部分代码省略.........
示例7: eregtesting_13TeV_Pi0
//.........这里部分代码省略.........
RooAbsReal *signlim = ws->function("signlim");
RooAbsReal *sign2lim = ws->function("sign2lim");
RooAbsReal *sigalphalim = ws->function("sigalphalim");
RooAbsReal *sigalpha2lim = ws->function("sigalpha2lim");
//formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
ecorvar->setRange(0.,2.);
ecorvar->setBins(800);
//formula for raw energy/true energy (1.0/(etrue/eraw))
RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
rawvar->setRange(0.,2.);
rawvar->setBins(800);
//clone data and add regression outputs for plotting
RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
RooRealVar *meanvar = (RooRealVar*)hdataclone->addColumn(*sigmeanlim);
RooRealVar *widthvar = (RooRealVar*)hdataclone->addColumn(*sigwidthlim);
RooRealVar *nvar = (RooRealVar*)hdataclone->addColumn(*signlim);
RooRealVar *n2var = (RooRealVar*)hdataclone->addColumn(*sign2lim);
RooRealVar *alphavar = (RooRealVar*)hdataclone->addColumn(*sigalphalim);
RooRealVar *alpha2var = (RooRealVar*)hdataclone->addColumn(*sigalpha2lim);
//plot target variable and weighted regression prediction (using numerical integration over reduced testing dataset)
TCanvas *craw = new TCanvas;
//RooPlot *plot = tgtvar->frame(0.6,1.2,100);
RooPlot *plot = tgtvar->frame(0.6,2.0,100);
hdata->plotOn(plot);
sigpdf->plotOn(plot,ProjWData(*hdatasmall));
plot->Draw();
craw->SaveAs("RawE.pdf");
craw->SaveAs("RawE.png");
craw->SetLogy();
plot->SetMinimum(0.1);
craw->SaveAs("RawElog.pdf");
craw->SaveAs("RawElog.png");
//plot distribution of regressed functions over testing dataset
TCanvas *cmean = new TCanvas;
RooPlot *plotmean = meanvar->frame(0.8,2.0,100);
hdataclone->plotOn(plotmean);
plotmean->Draw();
cmean->SaveAs("mean.pdf");
cmean->SaveAs("mean.png");
TCanvas *cwidth = new TCanvas;
RooPlot *plotwidth = widthvar->frame(0.,0.05,100);
hdataclone->plotOn(plotwidth);
plotwidth->Draw();
cwidth->SaveAs("width.pdf");
cwidth->SaveAs("width.png");
TCanvas *cn = new TCanvas;
RooPlot *plotn = nvar->frame(0.,111.,200);
hdataclone->plotOn(plotn);
plotn->Draw();
cn->SaveAs("n.pdf");
cn->SaveAs("n.png");
示例8: fitlambda
void fitlambda(bool realdata=1){
gSystem->Load("libRooFit");
using namespace RooFit;
TFile *data;
if (realdata==1) data=new TFile("./RecoRoutines_Z-selection_data36pb.rew2.corr1.root","read");
if (realdata==0) data=new TFile("./RecoRoutines_Z-selection_ZJets_TuneZ2_7TeV_alpgen_tauola.rew2.corr1.root","read");
TString dir="MuPtScale0toinf/";
TString title="dimuonstree";
TTree *tree = (TTree*)data->Get((dir+title).Data());
RooRealVar RMpPt("RMpPt","RMpPt",0,400);
RooRealVar RMpEta("RMpEta","RMpEta",-10,10);
RooRealVar RMpPhi("RMpPhi","RMpPhi",-10,10);
RooRealVar RMmPt("RMmPt","RMmPt",0,400);
RooRealVar RMmEta("RMmEta","RMmEta",-10,10);
RooRealVar RMmPhi("RMmPhi","RMmPhi",-10,10);
RooRealVar RZmass("RZmass","RZmass",81.2,101.2);
RooDataSet dataset("dataset","dataset",RooArgSet(RMpPt,RMpEta,RMpPhi,RMmPt,RMmEta,RMmPhi,RZmass),StoreError(RooArgSet(RZmass)));
Double_t MpPt ;
Double_t MpEta;
Double_t MpPhi;
Double_t MmPt;
Double_t MmEta;
Double_t MmPhi;
Double_t Zmass;
Double_t EvWeight;
tree->SetBranchAddress("MpPt",&MpPt);
tree->SetBranchAddress("MpEta",&MpEta);
tree->SetBranchAddress("MpPhi",&MpPhi);
tree->SetBranchAddress("MmPt",&MmPt);
tree->SetBranchAddress("MmEta",&MmEta);
tree->SetBranchAddress("MmPhi",&MmPhi);
tree->SetBranchAddress("Zmass",&Zmass);
tree->SetBranchAddress("EvWeight",&EvWeight);
Long64_t nentries = tree->GetEntries();
for (Long64_t i=0;i<nentries;i++) {
tree->GetEntry(i);
RMpPt=MpPt;
RMpEta=MpEta;
RMpPhi=MpPhi;
RMmPt=MmPt;
RMmEta=MmEta;
RMmPhi=MmPhi;
RZmass=Zmass;
RZmass.setError(3);
dataset.add(RooArgSet(RMpPt,RMpEta,RMpPhi,RMmPt,RMmEta,RMmPhi,RZmass),EvWeight);
}
RooRealVar mean("mean","mean",91.2,85.0,95.0);
RooRealVar cs01("cs01","cs01",0,-0.00,0.00);
RooRealVar cs02("cs02","cs02",0,-0.00,0.00);
RooRealVar ca01("ca01","ca01",0,-0.00,0.00);
RooRealVar ca02("ca02","ca02",0,-0.05,0.05);
RooRealVar cas1("cas1","cas1",0,-0.00,0.00);
RooRealVar cas2("cas2","cas2",0,-0.05,0.05);
RooRealVar phase("phase","phase",0,-3.14,3.14);
RooRealVar cae1("cae1","cae1",0,-0.00,0.00);
RooRealVar caeB2("caeB2","caeB2",0,-0.05,0.05);
RooRealVar caeE2("caeE2","caeE2",0,-0.05,0.05);
RooArgSet deps(RMpPt,RMmPt,RMpPhi,RMmPhi,RMpEta,RMmEta);
deps.add(mean);
deps.add(cs01);
deps.add(cs02);
deps.add(ca01);
deps.add(ca02);
deps.add(cas1);
deps.add(cas2);
deps.add(phase);
deps.add(cae1);
deps.add(caeB2);
deps.add(caeE2);
// corrfunction:
// dp/p = CS01+CS02*p+charge*(CA01+CA02*p)+charge*(CAS1+CAS2*p)*sin(phi+PHASE)+charge*(CAE1+CAE2*p)*eta
// 0.5*(caeB2+caeE2+TMath::Sign(caeB2,1-abs(eta))+TMath::Sign(caeE2,abs(eta)-1))
RooFormulaVar model("model","mean+mean/2*(cs01+cs02*RMpPt/100+(ca01+ca02*RMpPt/100)+(cas1+cas2*RMpPt/100)*sin(RMpPhi+phase)+\
(cae1+0.5*(caeB2+caeE2+TMath::Sign(caeB2,0.5-abs(RMpEta))+TMath::Sign(caeE2,abs(RMpEta)-0.5))*RMpPt/100)*RMpEta)\
+mean/2*(cs01+cs02*RMmPt/100-(ca01+ca02*RMmPt/100)-(cas1+cas2*RMmPt/100)*sin(RMmPhi+phase)-\
(cae1+0.5*(caeB2+caeE2+TMath::Sign(caeB2,0.5-abs(RMmEta))+TMath::Sign(caeE2,abs(RMmEta)-0.5))*RMmPt/100)*RMmEta)",deps);
RooFitResult *r = model.chi2FitTo(dataset,YVar(RZmass),Save());
r->Print();
RooPlot* frame = RMpPhi.frame() ;
dataset.plotOnXY(frame,YVar(RZmass));
//.........这里部分代码省略.........
示例9: TCanvas
//.........这里部分代码省略.........
h->SetTitle(titbuf);
sample.defineType(TString(sName));
hmap[sName.Data()] = h;
}
f->Close();
delete c;
// divide the binned data in categories according to the generated top quark mass
RooRealVar mass("m","Mass", 100, 500);
RooDataHist combData("combData", "combined data",mass, sample, hmap );
//the parameters to fit and the variable
RooRealVar g_mean_slope("#mu_{G}(slope)","g_mean_slope",0.01,0.,1.);
RooRealVar g_mean_shift("#mu_{G}(intercept)","g_mean_shift",162,100,180);
RooRealVar g_sigma_slope("#sigma_{G}(slope)","g_sigma_slope",0.01,0.,1.);
RooRealVar g_sigma_shift("#sigma_{G}(intercept)","g_sigma_shift",10,0.,25);
RooRealVar l_mean_slope("mpv_{L}(slope)","l_mean_slope",0.,0.,1.);//1,0,10);
RooRealVar l_mean_shift("mpv_{L}(intercept)","l_mean_shift",212,150,250);
RooRealVar l_sigma_slope("#sigma_{L}(slope)","l_sigma_slope",0.,0.,1.);//1,0,10);
RooRealVar l_sigma_shift("#sigma_{L}(intercept)","l_sigma_shift",10,0,25);
RooRealVar massfrac_slope("#alpha(slope)","massfrac_slope",0,0,0.01);
RooRealVar massfrac_shift("#alpha(intercept)","massfrac_shift",0.38,0.,1.);
//build the prototype pdf
RooRealVar topmass( "mtop","mtop",100,300);
RooFormulaVar g_mean( "g_mean", "(@0-172)*@[email protected]", RooArgSet(topmass,g_mean_slope,g_mean_shift));
RooFormulaVar g_sigma( "g_sigma", "(@0-172)*@[email protected]", RooArgSet(topmass,g_sigma_slope,g_sigma_shift));
RooGaussian gaus("gaus", "Mass component 1", mass, g_mean, g_sigma);
RooFormulaVar l_mean( "l_mean", "(@0-172)*@[email protected]", RooArgSet(topmass,l_mean_slope,l_mean_shift));
RooFormulaVar l_sigma( "l_sigma", "(@0-172)*@[email protected]", RooArgSet(topmass,l_sigma_slope,l_sigma_shift));
RooLandau lan("lan", "Mass component 2", mass, l_mean, l_sigma);
RooFormulaVar massfrac( "#alpha", "(@0-172)*@[email protected]", RooArgSet(topmass,massfrac_slope,massfrac_shift));
RooAddPdf massmodel("model","Model",RooArgList(lan,gaus),massfrac);
//RooNumConvPdf massmodel("model","Model",topmass,lan,gaus);
//now split per categories
RooSimPdfBuilder builder(massmodel) ;
RooArgSet* config = builder.createProtoBuildConfig() ;
config->setStringValue("physModels","model"); // Name of the PDF we are going to work with
config->setStringValue("splitCats","signal"); // Category used to differentiate sub-datasets
config->setStringValue("model","signal : mtop"); // Prescription to taylor PDF parameters mtop for each subset in signal
RooSimultaneous* simPdf = builder.buildPdf(*config,&combData) ;
config = simPdf->getParameters(combData);
for(size_t ipt=0; ipt<MassPointCollection.size(); ipt++)
{
TString sName("m"); sName+=(ipt+1);
Float_t imass=MassPointCollection[ipt].second;
(((RooRealVar &)(*config)["mtop_"+sName])).setRange(imass,imass);
(((RooRealVar &)(*config)["mtop_"+sName])).setVal(imass);
}
//fit to data
simPdf->fitTo(combData,Range(100.,400.));
//display
for(size_t ipt=0; ipt<MassPointCollection.size(); ipt++)
{
if(ipt%5==0)
{
TString name("SignalPDFs_"); name+=ipt;
c = new TCanvas(name,name);
c->SetBorderSize(0);
c->SetFillStyle(0);
c->SetFillColor(0);
c->SetWindowSize(1750,350);
c->Clear();
c->Divide(5,1);
}
TPad *p = (TPad *)c->cd(ipt%5+1);
p->SetGridx();
p->SetGridy();
TString procName("m"); procName += (ipt+1);
char buf[100];
sprintf(buf,"m_{t}=%3.1lf GeV/c^{2}",MassPointCollection[ipt].second);
RooPlot* frame = mass.frame(Title(buf));
RooDataSet* dataslice = (RooDataSet *)combData.reduce("signal==signal::"+procName);
dataslice->plotOn(frame,DataError(RooAbsData::SumW2));
RooCategory newCat(procName,procName);
simPdf->plotOn(frame,Slice(newCat),ProjWData(mass,*dataslice));
frame->GetYaxis()->SetTitleOffset(1.0);
frame->GetYaxis()->SetTitle("Events");
frame->GetXaxis()->SetTitleOffset(0.8);
frame->GetXaxis()->SetTitle("Reconstructed Mass [GeV/c^{2}]");
frame->Draw();
TPaveText *pt = new TPaveText(0.75,0.85,0.97,0.95,"brNDC");
pt->SetBorderSize(0);
pt->SetFillColor(0);
pt->SetFillStyle(0);
char buf2[50];
sprintf(buf2,"%3.1lf GeV/c^{2}",MassPointCollection[ipt].second);
pt->AddText(buf2);
pt->Draw();
}
return simPdf;
}
示例10: fit_and_weights_norm
//.........这里部分代码省略.........
alpha1.setConstant( true );
alpha2.setConstant( true );
frac2.setConstant( true );
n1.setConstant( true );
n2.setConstant( true );
sigma1.setConstant( true );
sigma2.setConstant( true );
// -- bg, mass shape
RooRealVar a1("a1","a1", -0.1, -0.5, 0.5);
RooChebychev comb("comb","comb", Lambda_b0_DTF_MASS_constr1, a1);
RooRealVar mean3("mean3","mean3", 5560., 5500., 5600.);
RooRealVar sigma3("sigma3","sigma3", 5., 1., 10.);
RooRealVar frac3("frac3","frac", 0.2, 0.0, 0.3);
RooGaussian gauss3("gauss3","gauss3", Lambda_b0_DTF_MASS_constr1, mean3, sigma3);
RooAddPdf bg("bg","bg", RooArgList(gauss3, comb), RooArgList(frac3));
// -- add signal & bg
RooAddPdf pdf("pdf", "pdf", RooArgList(sig, comb), RooArgList( sigYield, bgYield));
RooArgSet obs;
obs.add(Lambda_b0_DTF_MASS_constr1);
obs.add(Jpsi_M);
//obs.add(chi_c_M);
//obs.add(proton_ProbNNp);
//obs.add(proton_ProbNNk);
//obs.add(kaon_ProbNNp);
//obs.add(kaon_ProbNNk);
RooDataSet ds("ds","ds", obs, RooFit::Import(*tree));
RooPlot* plot = Lambda_b0_DTF_MASS_constr1.frame();
RooFitResult * result = pdf.fitTo( ds, RooFit::Extended() );
ds.plotOn( plot );
pdf.plotOn( plot );
RooPlot* plotPullMass = Lambda_b0_DTF_MASS_constr1.frame();
plotPullMass->addPlotable( plot->pullHist() );
//plotPullMass->SetMinimum();
//plotPullMass->SetMaximum();
TCanvas* c = new TCanvas();
TPad* pad1 = new TPad("pad1","pad1", 0, 0.3, 1, 1.0);
pad1->SetBottomMargin(0.0);
pad1->SetTopMargin(0.01);
pad1->Draw();
c->cd();
TPad* pad2 = new TPad("pad2","pad2", 0, 0., 1, 0.3);
pad2->SetBottomMargin(0.0);
pad2->SetTopMargin(0.0);
pad2->Draw();
pdf.plotOn( plot, RooFit::Components( sig ), RooFit::LineColor( kTeal ), RooFit::LineStyle(kDashed) );
pdf.plotOn( plot, RooFit::Components( comb ), RooFit::LineColor( kOrange ), RooFit::LineStyle(kDashed) );
pdf.plotOn( plot, RooFit::Components( gauss3 ), RooFit::LineColor( kViolet ), RooFit::LineStyle(kDashed) );
pad1->cd();
plot->Draw();
示例11: makefit
RooWorkspace* makefit(char* filename, char* FitTitle, char* Outfile, double minMass, double maxMass, double mean_bw, double gamma_bw, double cutoff_cb, const char* plotOpt, const int nbins) {
gROOT->ProcessLine(".L tdrstyle.C");
setTDRStyle();
gStyle->SetPadRightMargin(0.05);
//Create Data Set
RooRealVar mass("mass","m(EE)",minMass,maxMass,"GeV/c^{2}");
RooDataSet *data = RooDataSet::read(filename,RooArgSet(mass));
//====================== Parameters===========================
//Crystal Ball parameters
RooRealVar cbBias ("#Deltam_{CB}", "CB Bias", -.01, -10, 10, "GeV/c^{2}");
RooRealVar cbSigma("sigma_{CB}", "CB Width", 1.7, 0.02, 5.0, "GeV/c^{2}");
RooRealVar cbCut ("a_{CB}","CB Cut", 1.05, 0.1, 3.0);
RooRealVar cbPower("n_{CB}","CB Order", 2.45, 0.1, 20.0);
cbCut.setVal(cutoff_cb);
//Breit_Wigner parameters
RooRealVar bwMean("m_{Z}","BW Mean", 91.1876, "GeV/c^{2}");
bwMean.setVal(mean_bw);
RooRealVar bwWidth("#Gamma_{Z}", "BW Width", 2.4952, "GeV/c^{2}");
bwWidth.setVal(gamma_bw);
// Fix the Breit-Wigner parameters to PDG values
bwMean.setConstant(kTRUE);
bwWidth.setConstant(kTRUE);
// Exponential Background parameters
RooRealVar expRate("#lambda_{exp}", "Exponential Rate", -0.064, -1, 1);
RooRealVar c0("c_{0}", "c0", 1., 0., 50.);
//Number of Signal and Background events
RooRealVar nsig("N_{S}", "# signal events", 524, 0.1, 10000000000.);
RooRealVar nbkg("N_{B}", "# background events", 43, 1., 10000000.);
//============================ P.D.F.s=============================
// Mass signal for two decay electrons p.d.f.
RooBreitWigner bw("bw", "bw", mass, bwMean, bwWidth);
RooCBShape cball("cball", "Crystal Ball", mass, cbBias, cbSigma, cbCut, cbPower);
RooFFTConvPdf BWxCB("BWxCB", "bw X crystal ball", mass, bw, cball);
// Mass background p.d.f.
RooExponential bg("bg", "exp. background", mass, expRate);
// Mass model for signal electrons p.d.f.
RooAddPdf model("model", "signal", RooArgList(BWxCB), RooArgList(nsig));
TStopwatch t ;
t.Start() ;
model.fitTo(*data,FitOptions("mh"),Optimize(0),Timer(1));
t.Print() ;
TCanvas* c = new TCanvas("c","Unbinned Invariant Mass Fit", 0,0,800,600);
//========================== Plotting ============================
//Create a frame
RooPlot* plot = mass.frame(Range(minMass,maxMass),Bins(nbins));
data->plotOn(plot);
model.plotOn(plot);
model.paramOn(plot, Format(plotOpt, AutoPrecision(1)), Parameters(RooArgSet(cbBias, cbSigma, cbCut, cbPower, bwMean, bwWidth, expRate, nsig, nbkg)), Layout(0.66,0.63));
plot->getAttText()->SetTextSize(.03);
plot->Draw();
// Print Fit Values
TLatex *tex = new TLatex();
tex->SetNDC();
tex->SetTextSize(.04);
tex->SetTextFont(2);
tex->DrawLatex(0.195,0.875, "CMS ECAL, 2012");
tex->Draw();
tex->SetTextSize(0.022);
tex->DrawLatex(0.195, 0.81, FitTitle);
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->Print(Outfile);
RooWorkspace *w = new RooWorkspace("zfit");
w->import(model);
w->import(*data);
return w;
}
示例12: StandardHistFactoryPlotsWithCategories
//.........这里部分代码省略.........
RooRealVar * firstPOI = dynamic_cast<RooRealVar*>(mc->GetParametersOfInterest()->first());
firstPOI->setVal(muVal);
// firstPOI->setConstant();
if(doFit){
mc->GetPdf()->fitTo(*data);
}
// -------------------------------------------------------
mc->GetNuisanceParameters()->Print("v");
int nPlotsMax = 1000;
cout <<" check expectedData by category"<<endl;
RooDataSet* simData=NULL;
RooSimultaneous* simPdf = NULL;
if(strcmp(mc->GetPdf()->ClassName(),"RooSimultaneous")==0){
cout <<"Is a simultaneous PDF"<<endl;
simPdf = (RooSimultaneous *)(mc->GetPdf());
} else {
cout <<"Is not a simultaneous PDF"<<endl;
}
if(doFit) {
RooCategory* channelCat = (RooCategory*) (&simPdf->indexCat());
TIterator* iter = channelCat->typeIterator() ;
RooCatType* tt = NULL;
tt=(RooCatType*) iter->Next();
RooAbsPdf* pdftmp = ((RooSimultaneous*)mc->GetPdf())->getPdf(tt->GetName()) ;
RooArgSet* obstmp = pdftmp->getObservables(*mc->GetObservables()) ;
obs = ((RooRealVar*)obstmp->first());
RooPlot* frame = obs->frame();
cout <<Form("%s==%s::%s",channelCat->GetName(),channelCat->GetName(),tt->GetName())<<endl;
cout << tt->GetName() << " " << channelCat->getLabel() <<endl;
data->plotOn(frame,MarkerSize(1),Cut(Form("%s==%s::%s",channelCat->GetName(),channelCat->GetName(),tt->GetName())),DataError(RooAbsData::None));
Double_t normCount = data->sumEntries(Form("%s==%s::%s",channelCat->GetName(),channelCat->GetName(),tt->GetName())) ;
pdftmp->plotOn(frame,LineWidth(2.),Normalization(normCount,RooAbsReal::NumEvent)) ;
frame->Draw();
cout <<"expected events = " << mc->GetPdf()->expectedEvents(*data->get()) <<endl;
return;
}
int nPlots=0;
if(!simPdf){
TIterator* it = mc->GetNuisanceParameters()->createIterator();
RooRealVar* var = NULL;
while( (var = (RooRealVar*) it->Next()) != NULL){
RooPlot* frame = obs->frame();
frame->SetYTitle(var->GetName());
data->plotOn(frame,MarkerSize(1));
var->setVal(0);
mc->GetPdf()->plotOn(frame,LineWidth(1.));
var->setVal(1);
mc->GetPdf()->plotOn(frame,LineColor(kRed),LineStyle(kDashed),LineWidth(1));
var->setVal(-1);
mc->GetPdf()->plotOn(frame,LineColor(kGreen),LineStyle(kDashed),LineWidth(1));
list->Add(frame);
var->setVal(0);
}
示例13: Tbroomfit
//.........这里部分代码省略.........
}
RooAbsPdf* model = (RooAbsPdf*) cust.build(kTRUE) ; //build a clone...comment on changes...
// model->Print("t") ;
//delete model ; // eventualy delete the model...
/*
* DISPLAY RESULTS >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
*/
TPad *orig_gpad=(TPad*)gPad;
TCanvas *c;
c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("fitresult");
if (c==NULL){
printf("making new canvas\n%s","");
c=new TCanvas("fitresult",h2->GetName(),1000,700);
}else{
printf("using old canvas\n%s","");
c->SetTitle( h2->GetName() );
}
c->Clear();
printf(" canvas cleared\n%s","");
c->Divide(1,2) ;
printf(" canvas divided\n%s","");
c->Modified();c->Update();
RooDataHist datah("datah","datah with x",x,h2);
RooPlot* xframe = x.frame();
datah.plotOn(xframe, DrawOption("logy") );
// return;
if (TPRegexp("chi2").Match(command)!=0){//----------------------CHI2
//from lorenzo moneta
// TH1 * h1 = datah.createHistogram(x);
// TF1 * f = model->asTF(RooArgList(x) , parameters ); //???
// h2->Fit(f);
//It will work but you need to create a THNSparse and fit it
//or use directly the ROOT::Fit::BinData class to create a ROOT::Fit::Chi2Function to minimize.
// THIS CANNOT DO ZERO BINS
fitresult = model->chi2FitTo( datah , Save() );
}else{
// FIT FIT FIT FIT FIT FIT FIT FIT FIT FIT
fitresult = model->fitTo( datah , Save() );
}
fitresult->SetTitle( h2->GetName() ); // I PUT histogram name to global fitresult
// will be done by printResult ... fitresult->Print("v") ;
//duplicite fitresult->floatParsFinal().Print("s") ;
// later - after parsfinale .... : printResult();
// model->Print(); // not interesting........
model->plotOn(xframe, LineColor(kRed), DrawOption("l0z") );
//,Minos(kFALSE)
/*
* Posledni nakreslena vec je vychodiskem pro xframe->resid...?
* NA PORADI ZALEZI....
*/
示例14: higgsMassFit
//.........这里部分代码省略.........
//---------------------------------
// 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();
RooPlot* rooBKGPlot = x.frame();
rooBKGTotPdf -> plotOn(rooBKGPlot,LineColor(kBlack));
enum EColor color = (enum EColor)(BKGColors[0]);
rooBKGTotPdf -> plotOn(rooBKGPlot,Components(("rooBKGPdf_"+BKGShortNames[0]).c_str()),LineColor(color));
color = (enum EColor)(BKGColors[1]);
rooBKGTotPdf -> plotOn(rooBKGPlot,Components(("rooBKGPdf_"+BKGShortNames[1]).c_str()),LineColor(color));
rooSIGPdf -> plotOn(rooBKGPlot,LineColor(kBlack),LineStyle(1),LineWidth(1));
rooBKGPlot->Draw();
TH1F* BKGShapeHistoNorm = (TH1F*) BKGTotShapeHisto -> Clone();
BKGShapeHistoNorm -> Scale(1./BKGTotShapeHisto->Integral()/nBKG);
BKGShapeHistoNorm -> Draw("HIST,same");
char pngFileName[50];
sprintf(pngFileName,"BKGShapeNorm_H%d_%s.png",mH,mode.c_str());
c2 -> Print(pngFileName,"png");
}
//------------------------
// generate toy experiment
std::cout << "***********************************************************************" << std::endl;
std::cout << ">>> 1st toy experiment - " << mode << " mode" << std::endl;
int NBKGToy = int(BKGTotShapeHisto->Integral());
int NSIGToy = 0;
if( mode == "discovery" ) NSIGToy = int(SIGShapeHisto->Integral());
示例15: Gaussian_fit
//before: gROOT->ProcessLine(".include /afs/cern.ch/cms/slc5_ia32_gcc434/lcg/roofit/5.26.00-cms5/include")
//5_3_6: gROOT->ProcessLine(".include /afs/cern.ch/cms/slc5_amd64_gcc462/lcg/roofit/5.32.03-cms9/include/")
//Usage: .x Gaussian_fit.C+("Compare2012DMINECCnot","2012D_EtaRingMY_Vs_MYnoCC.root")
void Gaussian_fit( string Dir, string File ){
string PathL = Dir + "/" + File;
TCanvas* myc1 = new TCanvas("myc1", "CMS", 600, 600);
TString outname = Dir + "/Gaus_Fit.root";
//Files
cout<<"Opening: "<<PathL<<endl;
TFile* fin = TFile::Open(PathL.c_str());
//Histos
string List[4]={"Ratio1D_EB","Ratio1D_EEm","Ratio1D_EEp","Ratio1D_EB_eta1"};
for(int nH=0; nH<4; nH++){
TH1F *h_toFit = (TH1F*) fin->Get( List[nH].c_str() );
float hmean = h_toFit->GetMean();
float hrms = h_toFit->GetRMS();
TString outName = Dir + "/" + "Ratio1D_EB.png";
if(nH==1) outName = Dir + "/" + "Ratio1D_EEm.png";
if(nH==2) outName = Dir + "/" + "Ratio1D_EEp.png";
if(nH==3) outName = Dir + "/" + "Ratio1D_EB_eta1.png";
//Fit Method
RooRealVar x("x","IC distribution",hmean-1.6*hrms, hmean+1.6*hrms, "");
RooDataHist dh("dh","#gamma#gamma invariant mass",RooArgList(x),h_toFit);
RooRealVar mean("mean","mean",hmean, hmean-1.5*hrms,hmean+1.5*hrms,"");
RooRealVar sigma("sigma","#sigma",hrms, hrms-hrms/40.,hrms+hrms/40.,"");
mean.setRange(hmean-hmean/2.,hmean+hmean/2.);
sigma.setRange(0.,hrms+hrms/2.);
RooGaussian gaus("gaus","Core Gaussian",x, mean,sigma);
RooRealVar Nsig("Nsig","#pi^{0} yield",1000.,0.,1.e7);
Nsig.setVal( h_toFit->Integral()-h_toFit->Integral()/1000. );
Nsig.setRange(h_toFit->Integral()/50.,h_toFit->Integral());
RooRealVar cb0("cb0","cb0", 0.2, -1.,1.);
RooRealVar cb1("cb1","cb1",-0.1, -1.,1.);
RooRealVar cb2("cb2","cb2", 0.1, -1.,1.);
RooArgList cbpars(cb0,cb1,cb2);
RooChebychev bkg("bkg","bkg model", x, cbpars );
RooRealVar Nbkg("Nbkg","background yield",h_toFit->Integral()/100.,0.,h_toFit->Integral());
Nbkg.setVal( h_toFit->Integral()/100. );
RooAddPdf model1("model","only_gaus",RooArgList(gaus),RooArgList(Nsig));
RooAddPdf model2("model","sig+bkg",RooArgList(gaus,bkg),RooArgList(Nsig,Nbkg));
RooAbsPdf* model=0; model = &model2;
RooNLLVar nll("nll","log likelihood var",*model,dh);//,Extended());
RooMinuit m(nll);
m.setVerbose(kFALSE);
m.migrad();
//RooFitResult* res = m.save();
RooPlot* xframe = x.frame(h_toFit->GetNbinsX());
xframe->SetTitle("");
dh.plotOn(xframe);
model->plotOn(xframe);
//model->plotOn(xframe,Components(bkg),LineStyle(kDashed), LineColor(kRed));
h_toFit->Draw();
xframe->Draw("same");
//myc1->SaveAs(outName.Data());
float sigma_plot = sigma.getVal();
TLatex lat;
char line[300];
lat.SetNDC();
lat.SetTextSize(0.030);
lat.SetTextColor(1);
sprintf(line,"SIGMA: %.6f ", sigma_plot );
float xmin(0.55), yhi(0.80);// ypass(0.05);
lat.DrawLatex(xmin,yhi, line);
myc1->SaveAs(outName.Data());
}
//rms_EB->Write();
//rms_EEp->Write();
//rms_EEm->Write();
//output->Close();
}