// C o n s t r u c t A s i m o v d a t a s e t // --------------------------------------------- RooStats::ModelConfig mc("mc", combWS); mc.SetPdf("combModel"); mc.SetObservables(*(combWS->set("obs"))); mc.SetParametersOfInterest(poiSet); mc.SetNuisanceParameters(nuisSet); mc.SetGlobalObservables(globSet); mc.SetConstraintParameters(globSet); combWS->var("muZH_240GeV")->setVal(1); mc.SetSnapshot(*(combWS->var("muZH_240GeV"))); cout << combWS->var("muZH_240GeV")->getVal() << endl; mc.Print(); combWS->import(mc); RooArgSet *params2 = combModel->getParameters(*mrec); params2->Print("V"); RooArgSet combSet = RooArgSet(poiSet, nuisSet, "combSet"); const RooArgSet * genPoiValues = 0; RooArgSet genObs; RooAbsData * asimov = RooStats::AsymptoticCalculator::MakeAsimovData(*acombData, mc, poiSet, globSet, genPoiValues); RooCmdArg cmdList2; cmdList2.addArg(Range("range")); cmdList2.addArg(SplitRange(kTRUE)); cmdList2.addArg(Extended(kTRUE)); cmdList2.addArg(Offset(kTRUE)); cmdList2.addArg(NumCPU(8, 0)); cmdList2.addArg(GlobalObservables(globSet)); //cmdList2.addArg(SumCoefRange("range")); cmdList2.addArg(ExternalConstraints(smearSet)); cmdList2.addArg(Save(kTRUE)); cmdList2.addArg(Minimizer("Minuit2", "migrad")); cmdList2.addArg(Strategy(2)); //cmdList2.addArg(Minos(muZH_240GeV)); RooLinkedList linkedCmdList2 = cmdList2.subArgs(); RooFitResult * r98 = combModel->fitTo(*asimov, linkedCmdList2); r98->Print(); std::cout << "Asimov designed for mu = 1 gave mu-hat = " << combWS->var("muZH_240GeV")->getVal() << std::endl; combWS->import(*asimov, Rename("combined_Asimov")); RooArgSet *asimov_parameters = (RooArgSet*)combModel->getParameters(*mrec); combWS->saveSnapshot("asimovFit", *asimov_parameters, kTRUE); ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////// N O S Y S T E M A T I C S /////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// (combWS->var("Lmess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("Lmess_240GeV_eeH"))->setConstant(kTRUE); (combWS->var("COMmess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("COMmess_240GeV_eeH"))->setConstant(kTRUE); (combWS->var("BRmess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("BRmess_240GeV_eeH"))->setConstant(kTRUE); (combWS->var("Emess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("Emess_240GeV_eeH"))->setConstant(kTRUE); //(combWS->var("N_eeH"))->setConstant(kFALSE); //(combWS->var("N_mumuH"))->setConstant(kFALSE); (combWS->var("BRmess_240GeV_mumuH"))->setVal(1); (combWS->var("BRmess_240GeV_eeH"))->setVal(1); (combWS->var("Emess_240GeV_mumuH"))->setVal(1); (combWS->var("Emess_240GeV_eeH"))->setVal(1); (combWS->var("Lmess_240GeV_mumuH"))->setVal(1); (combWS->var("Lmess_240GeV_eeH"))->setVal(1); (combWS->var("COMmess_240GeV_mumuH"))->setVal(1); (combWS->var("COMmess_240GeV_eeH"))->setVal(1); RooCmdArg cmdList3; cmdList3.addArg(Range("range")); cmdList3.addArg(SplitRange(kTRUE)); cmdList3.addArg(Extended(kTRUE)); cmdList3.addArg(Offset(kTRUE)); cmdList3.addArg(NumCPU(8, 0)); //cmdList3.addArg(GlobalObservables(globSet)); //cmdList3.addArg(SumCoefRange("range")); //cmdList3.addArg(ExternalConstraints(smearSet)); cmdList3.addArg(Save(kTRUE)); cmdList3.addArg(Minimizer("Minuit2", "migrad")); cmdList3.addArg(Strategy(2)); //cmdList3.addArg(Minos(muZH_240GeV)); RooLinkedList linkedCmdList3 = cmdList3.subArgs(); RooFitResult* r5 = combModel_bare->fitTo(*acombDataBare, linkedCmdList3); r5->Print(); combWS->import(*combModel_bare); RooStats::ModelConfig mcB("mcB", combWS); mcB.SetPdf("combModel_bare"); mcB.SetObservables(*(combWS->set("obs"))); mcB.SetParametersOfInterest(poiSet); mcB.SetNuisanceParameters(nuisSet); combWS->var("muZH_240GeV")->setVal(1); mcB.SetSnapshot(*(combWS->var("muZH_240GeV"))); combWS->import(mcB); RooAbsData * combAsimov_bare = RooStats::AsymptoticCalculator::MakeAsimovData(*acombDataBare, mcB, poiSet, genObs, genPoiValues); RooFitResult * r6 = combWS->pdf("combModel_bare")->fitTo(*combAsimov_bare, linkedCmdList3); std::cout << "Asimov designed for mu = 1 gave mu-hat = " << combWS->var("muZH_240GeV")->getVal() << std::endl; r6->Print(); combWS->import(*combAsimov_bare, Rename("llH_Asimov_bare")); (combWS->var("Lmess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("Lmess_240GeV_eeH"))->setConstant(kFALSE); (combWS->var("COMmess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("COMmess_240GeV_eeH"))->setConstant(kFALSE); (combWS->var("BRmess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("BRmess_240GeV_eeH"))->setConstant(kFALSE); (combWS->var("Emess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("Emess_240GeV_eeH"))->setConstant(kFALSE); combWS->writeToFile("WScombinedModelZll_muZH.root"); // C o n s t r u c t p l a i n l i k e l i h o o d // --------------------------------------------------- TCanvas* c4 = new TCanvas ("c4", "c4", 100, 200, 700, 500); TCanvas* c5 = new TCanvas ("c5", "c5", 100, 200, 700, 500); // Construct unbinned likelihood + minimize with same settings as fitTo RooAbsReal * nll = combModel->createNLL(*asimov, Save(), Extended(kTRUE), ExternalConstraints(smearSet), Offset(kTRUE), NumCPU(8, 0), PrintEvalErrors(-1), GlobalObservables(globSet)); RooMinimizer minim(*nll); minim.setMinimizerType("Minuit2"); minim.optimizeConst(kTRUE); minim.setOffsetting(kTRUE); minim.setStrategy(2); minim.migrad(); minim.hesse(); // Create profile likelihood of my POI RooAbsReal *pll_muZH = nll->createProfile(muZH_240GeV); // Make sure that the combModel_bare is unconstrained (combWS->var("BRmess_240GeV_mumuH"))->setVal(1); (combWS->var("BRmess_240GeV_eeH"))->setVal(1); (combWS->var("Emess_240GeV_mumuH"))->setVal(1); (combWS->var("Emess_240GeV_eeH"))->setVal(1); (combWS->var("Lmess_240GeV_mumuH"))->setVal(1); (combWS->var("Lmess_240GeV_eeH"))->setVal(1); (combWS->var("COMmess_240GeV_mumuH"))->setVal(1); (combWS->var("COMmess_240GeV_eeH"))->setVal(1); (combWS->var("Lmess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("Lmess_240GeV_eeH"))->setConstant(kTRUE); (combWS->var("COMmess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("COMmess_240GeV_eeH"))->setConstant(kTRUE); (combWS->var("BRmess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("BRmess_240GeV_eeH"))->setConstant(kTRUE); (combWS->var("Emess_240GeV_mumuH"))->setConstant(kTRUE); (combWS->var("Emess_240GeV_eeH"))->setConstant(kTRUE); RooNLLVar nll_bare("nll_bare", "nll_bare", *combModel_bare, *combAsimov_bare, Extended(kTRUE), Offset(kTRUE), NumCPU(8, 0), Range("range"), SplitRange(kTRUE)); RooMinimizer minim2(nll_bare); minim2.setMinimizerType("Minuit2"); minim2.optimizeConst(kTRUE); minim2.setOffsetting(kTRUE); minim2.setStrategy(2); minim2.migrad(); minim2.hesse(); RooAbsReal *pll_muZH_bare = nll_bare.createProfile(muZH_240GeV); pll_muZH_bare->Print(); //Here, both PLL are OK (in the sense, that they have the "best fit value from the fit") pll_muZH->Print(); // Undo the unconstraining (combWS->var("Lmess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("Lmess_240GeV_eeH"))->setConstant(kFALSE); (combWS->var("COMmess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("COMmess_240GeV_eeH"))->setConstant(kFALSE); (combWS->var("BRmess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("BRmess_240GeV_eeH"))->setConstant(kFALSE); (combWS->var("Emess_240GeV_mumuH"))->setConstant(kFALSE); (combWS->var("Emess_240GeV_eeH"))->setConstant(kFALSE); //pll_muZH_bare->Print(); //pll_muZH->Print(); // C o n s t r u c t p r o f i l e l i k e l i h o o d i n m u Z H // ----------------------------------------------------------------------- // The profile likelihood estimator on nll for Nsig will minimize nll w.r.t // all floating parameters except Nsig for each evaluation c5 -> cd(); RooPlot *frame98 = muZH_240GeV.frame(Range(0.7, 1.3), Title("NLL in #mu ZH "), Bins(120)); frame98 -> GetXaxis() ->SetRangeUser(0.9, 1.1); //pll_muZH_bare->Print(); //pll_muZH->Print(); // Here, the "best fit value from the fit" gets disturbed pll_muZH-> plotOn(frame98, LineColor(kRed),LineStyle(kSolid)); pll_muZH_bare-> plotOn(frame98, LineColor(kBlack),LineStyle(kDashed)); frame98->GetYaxis()->SetTitleOffset(1.4); frame98->SetMinimum(0.); frame98->SetMaximum(4.); frame98->Draw(); c5->Draw(); c5->SaveAs("NLL_prf_muZH.pdf");