#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace RooFit; using namespace RooStats; using namespace std; void rf_simultaneous_runA_runB() { // ============================================================ // 1. Define nominal (global) yields and per–run nominal sub-yields. // ============================================================ // Nominal numbers for each run: const double nSigA_nom = 240; const double nBkgA_nom = 10000; const double nSigB_nom = 510; const double nBkgB_nom = 20000; // Nominal total events in each run: double eventsRunA = nSigA_nom + nBkgA_nom; double eventsRunB = nSigB_nom + nBkgB_nom; // Nominal bgk double nBkgTot = nBkgA_nom + nBkgB_nom; double nSigTot = nSigA_nom + nSigB_nom; RooRealVar x("x", "Observable", 50, 200); x.setRange("fitRange", 110, 150); x.setRange("fullRange", 50, 200); x.setRange("leftRange", 50, 110); x.setRange("rightRange", 150, 200); x.setBins(25); // Combined signal and single background yields and the fraction of signal, this is what we fit: RooRealVar nSig("nSig", "Signal Yield", nSigTot = 0 ? 10 : nSigTot, 0, 0.2 * nBkgTot); RooRealVar nBkgA("nBkgA", "Background Yield - runA", nBkgA_nom, 0.7 * nBkgA_nom, 1.3 * nBkgA_nom); RooRealVar nBkgB("nBkgB", "Background Yield - runB", nBkgB_nom, 0.7 * nBkgB_nom, 1.3 * nBkgB_nom); double fracA_val = nSigTot == 0 ? 0.5 : (nSigA_nom / nSigTot); RooRealVar fracA("fracA", "Fraction of signal runA events", fracA_val, 0, 1); // Roo formula var for Signal Yield RooFormulaVar nSigA("nSigA", "nSig*fracA", RooArgSet(nSig, fracA)); RooFormulaVar nSigB("nSigB", "nSig*(1-fracA)", RooArgSet(nSig, fracA)); // ============================================================ // 3. Build individual models for runA and runB. // ============================================================ // Run A model: RooRealVar meanA("meanA", "Signal Mean (runA)", 125); RooRealVar sigmaA("sigmaA", "Signal Width (runA)", 2); meanA.setConstant(kTRUE); sigmaA.setConstant(kTRUE); RooGaussian gaussA("gaussA", "Signal Gaussian (runA)", x, meanA, sigmaA); RooRealVar lambdaA("lambdaA", "Background Slope (runA)", -0.02); lambdaA.setConstant(kTRUE); RooExponential expA("expA", "Background Exponential (runA)", x, lambdaA); // The yields in the extended PDF come from our reparameterized formulas: RooAddPdf modelA_temp("modelA_temp", "RunA Model Temp", RooArgList(gaussA, expA), RooArgList(nSigA, nBkgA)); // Run B model: RooRealVar meanB("meanB", "Signal Mean (runB)", 125); RooRealVar sigmaB("sigmaB", "Signal Width (runB)", 4); meanB.setConstant(kTRUE); sigmaB.setConstant(kTRUE); RooGaussian gaussB("gaussB", "Signal Gaussian (runB)", x, meanB, sigmaB); RooRealVar lambdaB("lambdaB", "Background Slope (runB)", -0.01); lambdaB.setConstant(kTRUE); RooExponential expB("expB", "Background Exponential (runB)", x, lambdaB); RooAddPdf modelB_temp("modelB_temp", "RunB Model Temp", RooArgList(gaussB, expB), RooArgList(nSigB, nBkgB)); // ============================================================ // Fit background only to get the background yield nSig.setVal(0); nSig.setConstant(kTRUE); fracA.setVal(0.5); fracA.setConstant(kTRUE); RooDataSet *dataA = modelA_temp.generate(RooArgSet(x), Extended()); RooDataSet *dataB = modelB_temp.generate(RooArgSet(x), Extended()); RooFitResult *fitResultA = modelA_temp.fitTo(*dataA, Save(), Range("leftRange,rightRange"), PrintLevel(-1)); RooFitResult *fitResultB = modelB_temp.fitTo(*dataB, Save(), Range("leftRange,rightRange"), PrintLevel(-1)); // create a gaussian constraint for the background -- fit gives the full range estimates RooRealVar nBkgA_mean("nBkgA_mean", "Mean of nBkgA", nBkgA.getVal()); RooRealVar nBkgA_sigma("nBkgA_sigma", "Sigma of nBkgA", nBkgA.getError()); RooRealVar nBkgB_mean("nBkgB_mean", "Mean of nBkgB", nBkgB.getVal()); RooRealVar nBkgB_sigma("nBkgB_sigma", "Sigma of nBkgB", nBkgB.getError()); nBkgA_mean.setConstant(kTRUE); nBkgA_sigma.setConstant(kTRUE); nBkgB_mean.setConstant(kTRUE); nBkgB_sigma.setConstant(kTRUE); RooGaussian gaussBkgA("gaussBkgA", "Background Gaussian (runA)", nBkgA, nBkgA_mean, nBkgA_sigma); RooGaussian gaussBkgB("gaussBkgB", "Background Gaussian (runB)", nBkgB, nBkgB_mean, nBkgB_sigma); // Create the model multiplied by background constraints RooProdPdf modelA("modelA", "RunA Model with Bkg Constraint", RooArgList(modelA_temp, gaussBkgA)); RooProdPdf modelB("modelB", "RunB Model with Bkg Constraint", RooArgList(modelB_temp, gaussBkgB)); // list of the constraint RooArgList constraints; constraints.add(nBkgA); constraints.add(nBkgB); // list of global observables RooArgList globalObs; globalObs.add(nBkgA_mean); globalObs.add(nBkgB_mean); // 4. Build the simultaneous model. RooCategory runCat("runCat", "Run Category"); runCat.defineType("RunA"); runCat.defineType("RunB"); // ranges for categories runCat.setRange("fitRange", "RunA,RunB"); runCat.setRange("fullRange", "RunA,RunB"); RooSimultaneous simPdf("simPdf", "Simultaneous PDF", runCat); simPdf.addPdf(modelA, "RunA"); simPdf.addPdf(modelB, "RunB"); // ============================================================ // simPdf.Print("t"); // ============================================================ // reset values nSig.setVal(nSigTot); nSig.setConstant(kFALSE); fracA.setVal(fracA_val); fracA.setConstant(kFALSE); nBkgA.setVal(nBkgA_nom); nBkgA.setConstant(kFALSE); nBkgB.setVal(nBkgB_nom); nBkgB.setConstant(kFALSE); // data RooDataSet *combData = simPdf.generate(RooArgSet(x, runCat), Extended()); // fit the model to the data // set nSig away from its value to see if it converges nSig.setVal(nSigTot = 0 ? 0.01*nBkgTot : 1.2*nSigTot); RooFitResult *results = simPdf.fitTo(*combData, Save(), Range("fitRange"), //jonas github SumCoefRange("fullRange"), //jonas github SplitRange(), //jonas github EvalBackend("legacy"), //jonas github "cpu" or "legacy" does not change result Constrain(constraints), GlobalObservables(globalObs)); combData->table(runCat)->Print("V"); fitResultA->Print("v"); fitResultB->Print("v"); results->Print("v"); // --- RunA --- TCanvas *c1 = new TCanvas("c1", "Simultaneous Fit per Run", 1200, 400); c1->Divide(3, 1); // Plot for RunA c1->cd(1); RooPlot *xframeA = x.frame(Title("RunA")); combData->plotOn(xframeA, Cut("runCat==runCat::RunA")); // Draw signal+background (gaussA+expA) for RunA simPdf.plotOn(xframeA, Slice(runCat, "RunA"), ProjWData(runCat, *combData), Components("gaussA,expA"), LineStyle(kDashed), LineColor(kMagenta), Range("small")); // Draw signal only for RunA simPdf.plotOn(xframeA, Slice(runCat, "RunA"), ProjWData(runCat, *combData), Components("gaussA"), LineStyle(kDashed), LineColor(kGreen+2), Range("small")); // Draw background only for RunA simPdf.plotOn(xframeA, Slice(runCat, "RunA"), ProjWData(runCat, *combData), Components("expA"), LineStyle(kDotted), LineColor(kBlue), Range("small")); // Draw full model for RunA simPdf.plotOn(xframeA, Slice(runCat, "RunA"), ProjWData(runCat, *combData), LineColor(kRed), Range("small")); combData->plotOn(xframeA, Cut("runCat==runCat::RunA")); xframeA->Draw(); // --- RunB --- c1->cd(2); RooPlot *xframeB = x.frame(Title("RunB")); combData->plotOn(xframeB, Cut("runCat==runCat::RunB")); // Draw signal+background (gaussB+expB) for RunB simPdf.plotOn(xframeB, Slice(runCat, "RunB"), ProjWData(runCat, *combData), Components("gaussB,expB"), LineStyle(kDashed), LineColor(kMagenta), Range("small")); // Draw signal only for RunB simPdf.plotOn(xframeB, Slice(runCat, "RunB"), ProjWData(runCat, *combData), Components("gaussB"), LineStyle(kDashed), LineColor(kGreen+2), Range("small")); // Draw background only for RunB simPdf.plotOn(xframeB, Slice(runCat, "RunB"), ProjWData(runCat, *combData), Components("expB"), LineStyle(kDotted), LineColor(kRed), Range("small")); // Draw full model for RunB simPdf.plotOn(xframeB, Slice(runCat, "RunB"), ProjWData(runCat, *combData), LineColor(kBlue), Range("small")); combData->plotOn(xframeB, Cut("runCat==runCat::RunB")); xframeB->Draw(); // --- Combined --- c1->cd(3); RooPlot *xframeC = x.frame(Title("Combined")); combData->plotOn(xframeC); // Draw full model for combined simPdf.plotOn(xframeC, ProjWData(runCat, *combData), LineColor(kMagenta), Range("small")); // Draw signal only for combined simPdf.plotOn(xframeC, ProjWData(runCat, *combData), Components("expA,gaussA"), LineStyle(kDashed), LineColor(kBlue), Range("small")); // Draw background only for combined simPdf.plotOn(xframeC, ProjWData(runCat, *combData), Components("expB,gaussB"), LineStyle(kDotted), LineColor(kRed), Range("small")); combData->plotOn(xframeC); xframeC->Draw(); return; // test the model // Perform toy study with internal constraint on nBkgA and nBkgB RooMCStudy mcs(simPdf, RooArgSet(x, runCat), Constrain(constraints), Silence(), FitOptions(Range("small"), PrintLevel(-1))); // Run 100 toys, each with (nBkgA_small + nBkgB_small + nSigTot) events mcs.generateAndFit(1000, nBkgA_nom + nBkgB_nom + nSigA_nom + nSigB_nom); RooPlot *xframeMC = mcs.plotParam(nSig); auto can = new TCanvas(); can->cd(); xframeMC->Draw(); return; } /* // ============================================================ */ /* // fraction A double frac_blind_A = modelA.createIntegral(RooArgSet(x), NormSet(RooArgSet(x)), Range("small"))->getVal(); double frac_blind_B = modelB.createIntegral(RooArgSet(x), NormSet(RooArgSet(x)), Range("small"))->getVal(); // create a full model as sum of the two models RooAddPdf fullModel("fullModel", "Full Model", RooArgList(modelA, modelB)); // generate full data set RooDataSet *data = fullModel.generate(RooArgSet(x), Extended()); RooDataSet *data0 = (RooDataSet *)data->reduce(CutRange("small")); // fit in the full range and fit range RooFitResult *fitResult = fullModel.fitTo(*data, Save(), Range("full"), PrintLevel(1)); fullModel.removeStringAttribute("fitrange"); fitResult->Print("v"); // reset nuisance parameters nBkgA.setVal(nBkgA_nom*frac_blind_A); nBkgA.setRange(0, 3*nBkgA_nom*frac_blind_A); nBkgB.setVal(nBkgB_nom*frac_blind_B); nBkgB.setRange(0, 3*nBkgB_nom*frac_blind_B); nSig.setVal(nSigTot); fracA.setVal(0.5); fitResult = fullModel.fitTo(*data0, Save(), Range("small"), PrintLevel(1)); fitResult->Print("v"); */