using nlohmann::json; // Run with ROOT 6.28 (can be activated on lxplus with 'lb-run Moore/latest bash' ) // Argument of 'toyModel_simple.C' : // bool makeToyHistos ==> set 'true' to make histograms of data, signal and background // bool makeToyHistos ==> set 'false' to run fit // Run toymodel: // (1) make fitting histograms // root 'toyModel_simple.C(true)' // (2) run template fit // root 'toyModel_simple.C(false)' //============================================================================ // Create 1D histogram templates //============================================================================ void createToyHistos_1D(){ cout << "\n>>>>>> Create histogram templates..\n" << endl; // Save histograms in this root file TFile* histoFile = TFile::Open("histoFile.root", "RECREATE"); // Set number of events in signal and background (N_Sig = N_BG = N_DATA/2) int N_evt = 10000; // Set min/max range of observable x double x_min = 3000; double x_max = 6000; // Number of bins int N_bins = 15; // Functions used to fill histograms TF1* func1 = new TF1("func1","1 + 0.01 * x",3000, 6000); // background TF1* func2 = new TF1("func2","0.2 * exp(-( (x-4000.)*(x-4000.)/(2.0*200.*200.) ) )", 3000, 6000); // signal // Histograms TH1D* histoBG = new TH1D("histoBG","background histo", 15, 3000, 6000); TH1D* histoSig = new TH1D("histoSig","signal histo", 15, 3000, 6000); TH1D* histoData = new TH1D("histoData","data histo", 15, 3000, 6000); // define observable x double x_bg; // background double x_sig; // signal // Fill background and signal histograms for( int i = 0; i < N_evt ; ++i){ x_bg = func1->GetRandom(); x_sig = func2->GetRandom(); histoBG->Fill(x_bg); histoSig->Fill(x_sig); } // Fill data histogram for( int i = 0; i < N_evt ; ++i){ x_bg = func1->GetRandom(); histoData->Fill( x_bg ); x_sig = func2->GetRandom(); histoData->Fill( x_sig ); } cout << "Initial histogram yields : " << endl; cout << "N_Data = " << histoData->Integral() << endl; cout << "N_Sig = " << histoSig->Integral() << endl; cout << "N_BG = " << histoBG->Integral() << "\n" << endl; histoFile->Write(); histoFile->Close(); } //============================================================================ // Load histograms //============================================================================ std::vector getHistos(const std::string histoFileName){ cout << "\n>>>>>> Loading histograms fromroot file..\n" << endl; // Open root file TFile* histoFile = TFile::Open(histoFileName.c_str(), "READ"); // Define vector with histograms std::vector histos; // Get signal histogram and add to vector TH1D* histoSig = (TH1D*)histoFile->Get("histoSig")->Clone(); histos.push_back( histoSig ); // Get background histogram and add to vector TH1D* histoBG = (TH1D*)histoFile->Get("histoBG")->Clone(); histos.push_back( histoBG ); return histos; } //============================================================================ // Perform toymodel fit //============================================================================ void toyModel_simple(bool makeToyHistos = false){ cout << "\n>>>>>> Create histogram templates..\n" << endl; // -------------------------------------- // Make histograms // -------------------------------------- // make histograms if( makeToyHistos ){ TCanvas* c = new TCanvas(); createToyHistos_1D(); return; } // -------------------------------------- // Load histograms // -------------------------------------- // file with histograms std::string histoFileName = "histoFile.root"; // load histograms std::vector templateHistos = getHistos(histoFileName); // -------------------------------------- // Setup measurement & channel // -------------------------------------- // -- Define the measurement RooStats::HistFactory::Measurement B2D0MuNuMeas("B2D0MuNu", "B2D0MuNu fit"); B2D0MuNuMeas.SetExportOnly(true); // Tells histfactory to not run the fit B2D0MuNuMeas.SetPOI("num_histoSig"); // set to Bogus parma. of interest B2D0MuNuMeas.SetLumi(1.0); B2D0MuNuMeas.SetLumiRelErr(0.1); // -- Define the channel (in our case, there is only one) RooStats::HistFactory::Channel B2D0MuNu("BMCorr"); B2D0MuNu.SetData( "histoData", "histoFile.root" ); // -- Add the histograms to the channel std::cout << " " << std::endl; for( TH1D* histo : templateHistos ){ RooStats::HistFactory::Sample sample(histo->GetName(),histo->GetName(),histoFileName); sample.SetNormalizeByTheory(false); sample.AddNormFactor( "norm_"+(std::string)histo->GetName(), 1/histo->Integral(), 1/histo->Integral(), 1/histo->Integral() ); // fix normalisation sample.AddNormFactor( "num_"+(std::string)histo->GetName(), histo->Integral(), 1.0, histo->Integral()*2 ); // free yield // add histograms to channel B2D0MuNu.AddSample( sample ); } // add channel to measurements B2D0MuNuMeas.AddChannel( B2D0MuNu ); B2D0MuNuMeas.CollectHistograms(); // -------------------------------------- // Setting up workspace // -------------------------------------- // -- Define workspace RooWorkspace* workSpace = RooStats::HistFactory::MakeModelAndMeasurementFast(B2D0MuNuMeas, {.binnedFitOptimization=false}); // Get model manually RooStats::ModelConfig* modelConf = (RooStats::ModelConfig*)workSpace->obj("ModelConfig"); RooSimultaneous* model = (RooSimultaneous*)modelConf->GetPdf(); ((RooRealVar*)(modelConf->GetNuisanceParameters()->find("Lumi")))->setConstant(true); //fix the Lumi to not have it also fitted ((RooRealVar*)(modelConf->GetNuisanceParameters()->find("norm_histoSig")))->setConstant(true); //fix norm ((RooRealVar*)(modelConf->GetNuisanceParameters()->find("norm_histoBG")))->setConstant(true); //fix norm RooArgSet* obs = (RooArgSet*)modelConf->GetObservables(); RooCategory *idx = (RooCategory*)obs->find("channelCat"); RooAbsData* data = (RooAbsData*)workSpace->data("obsData"); RooRealVar* obs_x = (RooRealVar*)obs->find("obs_x_BMCorr"); // -------------------------------------- // Perform fit // -------------------------------------- RooFitResult *fitResult; RooAbsReal* nll = model->createNLL( *data , RooFit::Offset(kTRUE) );// RooMinuit* minuitHF = new RooMinuit( *nll ) ; minuitHF->setStrategy(2); minuitHF->setErrorLevel(0.5); minuitHF->optimizeConst(1); minuitHF->fit("smh"); fitResult = minuitHF->save("Result","Result"); fitResult->Print("v"); // -------------------------------------- // Plotting fit result // -------------------------------------- // plot data, total pdf and individual components RooPlot* frame = obs_x->frame(); data->plotOn(frame, RooFit::DataError(RooAbsData::Poisson), RooFit::MarkerSize(0.4), RooFit::DrawOption("ZP"), RooFit::Name("pdf_Data")); model->plotOn(frame, RooFit::ProjWData(*idx, *data), RooFit::Name("pdf_model"), RooFit::Components("*histoSig*,*histoBG*"), RooFit::DrawOption("L"), RooFit::LineColor(kGreen)); model->plotOn(frame, RooFit::ProjWData(*idx, *data), RooFit::Name("pdf_sig"), RooFit::Components("*histoSig*"), RooFit::DrawOption("L"), RooFit::LineColor(kBlue)); model->plotOn(frame, RooFit::ProjWData(*idx, *data), RooFit::Name("pdf_bg"), RooFit::Components("*histoBG*"), RooFit::DrawOption("L"), RooFit::LineColor(kRed)); // set axis frame->GetYaxis()->SetRangeUser(0.0, 5.5e3); frame->GetXaxis()->SetTitle(" "); frame->SetTitle(" "); // - - - - - - legend TLegend* legend = new TLegend(0.70, 0.4, 0.90, 0.96); legend->SetBorderSize(0); legend->SetFillColorAlpha(kWhite, 0.0); legend->AddEntry(frame->findObject("pdf_Data"), "Data", "lep"); legend->AddEntry(frame->findObject("pdf_model"), "Fit", "l"); legend->AddEntry(frame->findObject("pdf_sig"), "Signal", "l"); legend->AddEntry(frame->findObject("pdf_bg"), "Background", "l"); legend->Draw(); // - - - - - - combined plot RooPlot* frame_tot = obs_x->frame(); data->plotOn(frame_tot, RooFit::DataError(RooAbsData::Poisson), RooFit::MarkerSize(0.4), RooFit::DrawOption("ZP"), RooFit::Name("pdf_Data")); model->plotOn(frame_tot, RooFit::ProjWData(*idx, *data), RooFit::Name("total"), RooFit::Components("*histoSig*,*histoBG*"), RooFit::DrawOption("L"), RooFit::LineColor(kGreen)); // - - - - - - pull plot // pull mine/max for lines double min = 3000; double max = 6000; RooPlot* pull = obs_x->frame(); RooHist* h1 = frame_tot->pullHist(); pull->addPlotable(h1,"E1"); pull->GetXaxis()->SetLabelSize(0.13); pull->GetYaxis()->SetLabelSize(0.13); pull->GetYaxis()->SetNdivisions(5); pull->SetTitle(" "); pull->GetXaxis()->SetTitle(" "); pull->SetMaximum(4.5); pull->SetMinimum(-4.5); TLine l_up(min, 2., max, 2.); TLine l_m(min, 0., max, 0.); TLine l_down(min, -2., max, -2.); l_up.SetLineColor(kRed); l_down.SetLineColor(kRed); l_m.SetLineStyle(3); // - - - - - - plot fit and pull plots together TCanvas c("Fit_and_pull", "Fit and pull", 900, 600); c.Divide(1,2); c.cd(1); gPad->SetPad(0.0, 0.20, 1.0, 1.0); gPad->SetRightMargin(0.025); gPad->SetTopMargin(0.02); frame->Draw(); legend->Draw("same"); c.cd(2); gPad->SetPad(0.0, 0.00, 1.0, 0.25); pull->GetXaxis()->SetTitle("obs_x_BMCorr"); pull->GetXaxis()->SetTitleOffset(0.9); pull->GetXaxis()->SetTitleSize(0.17); gPad->SetRightMargin(0.025); gPad->SetBottomMargin(0.40); pull->Draw("e1"); l_up.Draw("same"); l_down.Draw("same"); l_m.Draw("same"); c.Update(); c.SaveAs( "obs_x_BMCorr.pdf" ); }