////////////////////////////////////////////////////////////////////////// // // Modified based on RooFit tutorial macro #501 rf501_simultaneouspdf.C // // My question is: // 1. I use Least-Chi2 method to fit data with model, fit data_ctl with model_ctl, // and obtain the minimum FCN_1min and FCN_2min, respectively; // 2. I use Least-Chi2 method to simultaneously fit data with model and data_ctl with model_ctl, // and obtain the FCN_min of the simultaneous fit; (in the simultaneous fit, the parameters // of model and model_ctl are free); // 3. I find FCN_min = FCN_1min + FCN_2min, which is same as mathematical prediction; // 4. in mathematics, a simultaneously fit means to search the minmum value of FCN = FCN_1 + FCN_2; // it should be FCN_min = FCN_1min + FCN_2min; because FCN_1min and FCN_2min are the minimum // value of each part; the minimum value of FCN shoule be FCN_1min + FCN_2min; // 5. However, when I use real data, for Least-Chi2 method I find FCN_min != FCN_1min + FCN_2min, // so I am very confused now; // // wencheng@mail.ustc.edu.cn; yanwc@ihep.ac.cn; // ///////////////////////////////////////////////////////////////////////// #ifndef __CINT__ #include "RooGlobalFunc.h" #endif #include "RooRealVar.h" #include "RooDataSet.h" #include "RooGaussian.h" #include "RooConstVar.h" #include "RooChebychev.h" #include "RooAddPdf.h" #include "RooSimultaneous.h" #include "RooCategory.h" #include "TCanvas.h" #include "TAxis.h" #include "RooPlot.h" gSystem->Load("libRooFit"); using namespace RooFit; using namespace std; void rf501_simultaneouspdf_chi2() { // C r e a t e m o d e l f o r p h y s i c s s a m p l e // ------------------------------------------------------------- // Create observables RooRealVar x("x","x",-8,8) ; // Construct signal pdf RooRealVar mean("mean","mean",0,-8,8) ; RooRealVar sigma("sigma","sigma",0.3,0.1,10) ; RooGaussian gx("gx","gx",x,mean,sigma) ; // Construct background pdf RooRealVar a0("a0","a0",-0.1,-1,1) ; RooRealVar a1("a1","a1",0.004,-1,1) ; RooChebychev px("px","px",x,RooArgSet(a0,a1)) ; // Construct Extended composite pdf RooRealVar Ngx("Ngx","Ngx", 500, 0., 20000.); RooRealVar Npx("Npx","Npx", 500, 0., 20000.); RooAddPdf model("model","model",RooArgList(gx,px),RooArgList(Ngx,Npx)) ; // C r e a t e m o d e l f o r c o n t r o l s a m p l e // -------------------------------------------------------------- // Construct signal pdf. // NOTE that sigma is shared with the signal sample model RooRealVar mean_ctl("mean_ctl","mean_ctl",-3,-8,8) ; RooRealVar sigma_ctl("sigma_ctl","sigma",0.5,0.1,10) ; RooGaussian gx_ctl("gx_ctl","gx_ctl",x,mean_ctl,sigma_ctl) ; // Construct the background pdf RooRealVar a0_ctl("a0_ctl","a0_ctl",-0.1,-1,1) ; RooRealVar a1_ctl("a1_ctl","a1_ctl", 0.5,-0.1,1) ; RooChebychev px_ctl("px_ctl","px_ctl",x,RooArgSet(a0_ctl,a1_ctl)) ; // Construct the composite model RooRealVar Ngx_ctl("Ngx_ctl","Ngx_ctl", 1000, 0., 20000.); RooRealVar Npx_ctl("Npx_ctl","Npx_ctl", 1000, 0., 20000.); RooAddPdf model_ctl("model_ctl","model_ctl",RooArgList(gx_ctl,px_ctl), RooArgList(Ngx_ctl,Npx_ctl)) ; // G e n e r a t e e v e n t s f o r b o t h s a m p l e s // --------------------------------------------------------------- // Generate binned data from 1000 events in x and y from model // RooDataSet *data = model.generate(RooArgSet(x),1000) ; // RooDataSet *data_ctl = model_ctl.generate(RooArgSet(x),2000) ; RooDataHist *data = model.generateBinned(RooArgSet(x),1000) ; RooDataHist *data_ctl = model_ctl.generateBinned(RooArgSet(x),2000) ; //RooDataHist *data_ctl = px_ctl.generateBinned(RooArgSet(x),1000) ; // C r e a t e i n d e x c a t e g o r y a n d j o i n s a m p l e s // --------------------------------------------------------------------------- // Define category to distinguish physics and control samples events RooCategory sample("sample","sample") ; sample.defineType("physics") ; sample.defineType("control") ; // Construct combined binned dataset in (x,sample) //RooDataSet combData("combData","combined data",x,Index(sample),Import("physics",*data),Import("control",*data_ctl)) ; RooDataHist combData("combData","combined data",x,Index(sample),Import("physics",*data),Import("control",*data_ctl)) ; // C o n s t r u c t a s i m u l t a n e o u s p d f i n ( x , s a m p l e ) // ----------------------------------------------------------------------------------- // Construct a simultaneous pdf using category sample as index RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ; // Associate model with the physics state and model_ctl with the control state simPdf.addPdf(model,"physics") ; simPdf.addPdf(model_ctl,"control") ; // P e r f o r m a s i m u l t a n e o u s f i t // --------------------------------------------------- // Perform Chi2 fit of model to data and model_ctl to data_ctl and simultaneous fit RooFitResult* result_data = model.chi2FitTo(*data, Extended(kTRUE),Save(kTRUE)); RooFitResult* result_data_ctl = model_ctl.chi2FitTo(*data_ctl, Extended(kTRUE),Save(kTRUE)); RooFitResult* result_comData = simPdf.chi2FitTo(combData, Extended(kTRUE),Save(kTRUE)); // P l o t m o d e l s l i c e s o n d a t a s l i c e s // ---------------------------------------------------------------- // Make a frame for the physics sample RooPlot* frame1 = x.frame(Bins(30),Title("Physics sample")) ; // Plot all data tagged as physics sample combData.plotOn(frame1,Cut("sample==sample::physics")) ; // Plot "physics" slice of simultaneous pdf. // NBL You _must_ project the sample index category with data using ProjWData // as a RooSimultaneous makes no prediction on the shape in the index category // and can thus not be integrated simPdf.plotOn(frame1,Slice(sample,"physics"),ProjWData(sample,combData)) ; simPdf.plotOn(frame1,Slice(sample,"physics"),Components("px"),ProjWData(sample,combData),LineStyle(kDashed)) ; // The same plot for the control sample slice RooPlot* frame2 = x.frame(Bins(30),Title("Control sample")) ; combData.plotOn(frame2,Cut("sample==sample::control")) ; simPdf.plotOn(frame2,Slice(sample,"control"),ProjWData(sample,combData)) ; simPdf.plotOn(frame2,Slice(sample,"control"),Components("px_ctl"),ProjWData(sample,combData),LineStyle(kDashed)) ; TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf",800,400) ; c->Divide(2) ; c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ; c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ; cout<<"FCN of data fit="<minNll()<