rf510_wsnamedsets.C

Go to the documentation of this file.
00001 /////////////////////////////////////////////////////////////////////////
00002 //
00003 // 'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #510
00004 // 
00005 //  Working with named parameter sets and parameter snapshots in 
00006 //  workspaces
00007 //
00008 //  04/2009 - Wouter Verkerke 
00009 //
00010 /////////////////////////////////////////////////////////////////////////
00011 
00012 #ifndef __CINT__
00013 #include "RooGlobalFunc.h"
00014 #endif
00015 #include "RooRealVar.h"
00016 #include "RooDataSet.h"
00017 #include "RooGaussian.h"
00018 #include "RooConstVar.h"
00019 #include "RooChebychev.h"
00020 #include "RooAddPdf.h"
00021 #include "RooWorkspace.h"
00022 #include "RooPlot.h"
00023 #include "TCanvas.h"
00024 #include "TAxis.h"
00025 #include "TFile.h"
00026 #include "TH1.h"
00027 using namespace RooFit ;
00028 
00029 
00030 void fillWorkspace(RooWorkspace& w) ;
00031 
00032 void rf510_wsnamedsets()
00033 {
00034   // C r e a t e   m o d e l   a n d   d a t a s e t
00035   // -----------------------------------------------
00036 
00037   RooWorkspace* w = new RooWorkspace("w") ;
00038   fillWorkspace(*w) ;
00039 
00040   // Exploit convention encoded in named set "parameters" and "observables"
00041   // to use workspace contents w/o need for introspected
00042   RooAbsPdf* model = w->pdf("model") ;
00043 
00044   // Generate data from p.d.f. in given observables
00045   RooDataSet* data = model->generate(*w->set("observables"),1000) ;
00046 
00047   // Fit model to data
00048   model->fitTo(*data) ;
00049   
00050   // Plot fitted model and data on frame of first (only) observable
00051   RooPlot* frame = ((RooRealVar*)w->set("observables")->first())->frame() ;
00052   data->plotOn(frame) ;
00053   model->plotOn(frame) ;
00054 
00055   // Overlay plot with model with reference parameters as stored in snapshots
00056   w->loadSnapshot("reference_fit") ;
00057   model->plotOn(frame,LineColor(kRed)) ;
00058   w->loadSnapshot("reference_fit_bkgonly") ;
00059   model->plotOn(frame,LineColor(kRed),LineStyle(kDashed)) ;
00060 
00061 
00062   // Draw the frame on the canvas
00063   new TCanvas("rf510_wsnamedsets","rf503_wsnamedsets",600,600) ;
00064   gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame->Draw() ;
00065 
00066 
00067   // Print workspace contents
00068   w->Print() ;
00069 
00070 
00071   // Workspace will remain in memory after macro finishes
00072   gDirectory->Add(w) ;
00073 
00074 }
00075 
00076 
00077 
00078 void fillWorkspace(RooWorkspace& w) 
00079 {
00080   // C r e a t e   m o d e l
00081   // -----------------------
00082 
00083   // Declare observable x
00084   RooRealVar x("x","x",0,10) ;
00085 
00086   // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their paramaters
00087   RooRealVar mean("mean","mean of gaussians",5,0,10) ;
00088   RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
00089   RooRealVar sigma2("sigma2","width of gaussians",1) ;
00090 
00091   RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;  
00092   RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;  
00093   
00094   // Build Chebychev polynomial p.d.f.  
00095   RooRealVar a0("a0","a0",0.5,0.,1.) ;
00096   RooRealVar a1("a1","a1",-0.2,0.,1.) ;
00097   RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;
00098 
00099   // Sum the signal components into a composite signal p.d.f.
00100   RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
00101   RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;
00102 
00103   // Sum the composite signal and background 
00104   RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
00105   RooAddPdf  model("model","g1+g2+a",RooArgList(bkg,sig),bkgfrac) ;
00106 
00107   // Import model into p.d.f.
00108   w.import(model) ;
00109 
00110 
00111   // E n c o d e   d e f i n i t i o n   o f   p a r a m e t e r s   i n   w o r k s p a c e
00112   // ---------------------------------------------------------------------------------------
00113 
00114 
00115   // Define named sets "parameters" and "observables", which list which variables should be considered
00116   // parameters and observables by the users convention
00117   // 
00118   // Variables appearing in sets _must_ live in the workspace already, or the autoImport flag
00119   // of defineSet must be set to import them on the fly. Named sets contain only references
00120   // to the original variables, therefore the value of observables in named sets already
00121   // reflect their 'current' value
00122   RooArgSet* params = (RooArgSet*) model.getParameters(x) ;
00123   w.defineSet("parameters",*params) ;
00124   w.defineSet("observables",x) ;
00125 
00126 
00127   // E n c o d e   r e f e r e n c e   v a l u e   f o r   p a r a m e t e r s   i n   w o r k s p a c e
00128   // ---------------------------------------------------------------------------------------------------
00129 
00130 
00131   // Define a parameter 'snapshot' in the p.d.f.
00132   // Unlike a named set, a parameter snapshot stores an independent set of values for
00133   // a given set of variables in the workspace. The values can be stored and reloaded
00134   // into the workspace variable objects using the loadSnapshot() and saveSnapshot()
00135   // methods. A snapshot saves the value of each variable, any errors that are stored
00136   // with it as well as the 'Constant' flag that is used in fits to determine if a 
00137   // parameter is kept fixed or not.
00138 
00139   // Do a dummy fit to a (supposedly) reference dataset here and store the results
00140   // of that fit into a snapshot
00141   RooDataSet* refData = model.generate(x,10000) ;
00142   model.fitTo(*refData,PrintLevel(-1)) ;
00143   
00144   // The kTRUE flag imports the values of the objects in (*params) into the workspace
00145   // If not set, the present values of the workspace parameters objects are stored
00146   w.saveSnapshot("reference_fit",*params,kTRUE) ;
00147 
00148   // Make another fit with the signal componentforced to zero
00149   // and save those parameters too
00150 
00151   bkgfrac.setVal(1) ;
00152   bkgfrac.setConstant(kTRUE) ;
00153   bkgfrac.removeError() ;
00154   model.fitTo(*refData,PrintLevel(-1)) ;  
00155   
00156   w.saveSnapshot("reference_fit_bkgonly",*params,kTRUE) ;
00157   
00158 
00159 }

Generated on Tue Jul 5 15:45:08 2011 for ROOT_528-00b_version by  doxygen 1.5.1