From 9aec74381afe5ca54c156b678a65807271deba27 Mon Sep 17 00:00:00 2001 From: arabella Date: Mon, 9 Mar 2020 09:35:40 +0100 Subject: [PATCH 1/3] AxE and simultaneous fits --- AxE_MC_withFits_root.C | 232 +++++++++++++++++ README.md | 54 +++- makePlot_fitPeak_unbinned.py | 5 +- makeSimultaneousFits.C | 467 ++++++++++++++++++++++++++++++++++ models.cc => unused/models.cc | 0 models.h => unused/models.h | 0 6 files changed, 755 insertions(+), 3 deletions(-) create mode 100644 AxE_MC_withFits_root.C create mode 100644 makeSimultaneousFits.C rename models.cc => unused/models.cc (100%) rename models.h => unused/models.h (100%) diff --git a/AxE_MC_withFits_root.C b/AxE_MC_withFits_root.C new file mode 100644 index 0000000..7b0665f --- /dev/null +++ b/AxE_MC_withFits_root.C @@ -0,0 +1,232 @@ +//root AxE_MC_withFits_root.C'(1, 1, 1)' + +//look at MC to compute AxE = nEvents surviving a selection in a q2 bin / original nEvents in MC +//compute AxE by counting (events in q2 bin) and by fitting +//results are AxE and signal model +//save results in workspace to be loaded in the final fits + +//roofit +#include "RooRealVar.h" +#include "RooDataSet.h" +#include "RooDataHist.h" +#include "RooGaussian.h" +#include "RooChebychev.h" +#include "RooCBShape.h" +#include "RooArgusBG.h" +#include "RooFFTConvPdf.h" +#include "RooAddPdf.h" +#include "RooWorkspace.h" +#include "RooFitResult.h" +#include "RooFormulaVar.h" +#include "TText.h" +#include "TCanvas.h" +#include "TAxis.h" +#include "TH1F.h" +#include "TH1D.h" +#include "RooPlot.h" + + #include + #include + #include + #include + + +using namespace RooFit; +void AxE_MC_withFits_root(int isMC=1, int isEE=1, int isResonant=0, float mvaCut = 12.68){ + + gSystem->Load("../HiggsAnalysis/CombinedLimit/lib/libHiggsAnalysisCombinedLimit.so"); + + // float LumiNonReso = 1.53E+04; // fb-1 + // float LumiReso = 45.16919828; // fb-1 + + float NOriginalNonReso = 7587169; + float NOriginalReso = 2122456; + + + int nBins = 3; + std::vector> mll_bins_Range; + mll_bins_Range.push_back(std::pair(1.1, 2.4)); + mll_bins_Range.push_back(std::pair(2.4, 4.)); + mll_bins_Range.push_back(std::pair(1.1, 4.)); + + + //current files from Otto, need to adapt to other's outputs + std::string inputFileList = "/afs/cern.ch/user/k/klau/myWorkspace/public/ForArabella/03Mar2020_pf/"; + if(!isResonant && isEE) + inputFileList += "RootTree_2020Jan16_BuToKee_BToKEEAnalyzer_2020Feb18_fullq2_pf_isoPFMVADphiptImb_weighted_pauc02_mva.root"; + else if(isResonant && isEE) + inputFileList += "RootTree_2020Jan16_BuToKJpsi_Toee_BToKEEAnalyzer_2020Feb18_fullq2_pf_isoPFMVADphiptImb_weighted_pauc02_mva.root"; + //need to add for muons + + + std::cout << " isEE = " << isEE << " isResonant = " << isResonant << std::endl; + std::cout << "inputFileList = " << inputFileList << std::endl; + + + TChain* tree = new TChain("tree"); + tree->Add(inputFileList.c_str()); + + int nEntriesTot = tree->GetEntries(); + std::cout << " loaded nEvents = " << nEntriesTot << std::endl; + + float B_mll_fullfit; + float B_fit_mass; + float B_bdt_weight; + + tree->SetBranchStatus("*", 0); + if(isEE){ + tree->SetBranchStatus("BToKEE_mll_fullfit", 1); tree->SetBranchAddress("BToKEE_mll_fullfit", &B_mll_fullfit); + tree->SetBranchStatus("BToKEE_fit_mass", 1); tree->SetBranchAddress("BToKEE_fit_mass", &B_fit_mass); + tree->SetBranchStatus("BToKEE_xgb", 1); tree->SetBranchAddress("BToKEE_xgb", &B_bdt_weight); + } + else{ + tree->SetBranchStatus("BToKMuMu_mll_fullfit", 1); tree->SetBranchAddress("BToKMuMu_mll_fullfit", &B_mll_fullfit); + tree->SetBranchStatus("BToKMuMu_fit_mass", 1); tree->SetBranchAddress("BToKMuMu_fit_mass", &B_fit_mass); + tree->SetBranchStatus("BToKMuMu_xgb", 1); tree->SetBranchAddress("BToKMuMu_xgb", &B_bdt_weight); + } + + //loop over tree load hostos and RooDataSet + TH1F* hHisto[3]; + RooWorkspace w("signalModel"); + RooRealVar x("x", "", 4.4, 6.); + RooDataSet* data[3]; + + int ijC = 0; + for(auto ij: mll_bins_Range){ + hHisto[ijC] = new TH1F(Form("hHisto_mll_%.1f-%.1f", ij.first, ij.second), "", 500, 4.4, 6.); + + data[ijC] = new RooDataSet(Form("data_mll_%.1f-%.1f", ij.first, ij.second), + Form("data_mll_%.1f-%.1f", ij.first, ij.second), RooArgSet(x)); + ++ijC; + } + for(int ij=0; ijGetEntry(ij); + + if(B_bdt_weight <= mvaCut) continue; + int ijC = 0; + for(auto ijB: mll_bins_Range){ + if(B_mll_fullfit >= ijB.first && B_mll_fullfit < ijB.second) { + x = B_fit_mass; + data[ijC]->add(RooArgSet(x)); + hHisto[ijC]->Fill(B_fit_mass); + } + ++ijC; + } + } + w.import(x); + + //not really needed. Just to take a look at histos + TFile out(Form("AxE_isEE%d_isResonant%d.root", isEE, isResonant), "recreate"); + out.cd(); + for(auto ij=0; ijfitTo(*(data[ij]), Minimizer("Minuit2"),Save(true)); + + RooPlot * plot = w.var("x")->frame(); + if(isEE){ + if(isResonant) plot->SetXTitle("K(JPsi)ee mass (GeV)"); + else plot->SetXTitle("Kee mass (GeV)"); + } + plot->SetTitle(""); + data[ij]->plotOn(plot, Binning(200)); + model->plotOn(plot); + model->paramOn(plot, RooFit::Layout(0.2,0.4,0.9),RooFit::Format("NEA",AutoPrecision(1))); + plot->getAttLine()->SetLineColorAlpha(kWhite, 0.2); + plot->getAttText()->SetTextSize(0.03); + plot->getAttText()->SetTextFont(42); + + float chi2 = plot->chiSquare(6); + + + TCanvas * cc = new TCanvas(); + cc->SetLogy(0); + plot->Draw(); + + TLatex tL; + tL.SetNDC(); + tL.SetTextSize(0.04); + tL.SetTextFont(42); + tL.DrawLatex(0.65,0.8, Form("chi2 = %.1f ",chi2)); + TLatex tL1; + tL1.SetNDC(); + tL1.SetTextSize(0.05); + tL1.SetTextFont(42); + tL1.DrawLatex(0.65,0.7, Form("mll in %.1f-%.1f", mll_bins_Range[ij].first, mll_bins_Range[ij].second)); + std::string outName = "plots/Bmass_MC/"; + if(isEE) outName += "Kee"; + else outName += "Kmm"; + outName += Form("_mll_%.1f-%.1f", mll_bins_Range[ij].first, mll_bins_Range[ij].second); + cc->Print((outName+".png").c_str(), "png"); + cc->Print((outName+".pdf").c_str(), "pdf"); + cc->Print((outName+".root").c_str(), "root"); + + RooRealVar* parS = (RooRealVar*) r->floatParsFinal().find(("nSig"+suffixName).c_str()); + nEv_postFit[ij] = parS->getValV(); + nEvError_postFit[ij] = parS->getError(); + RooRealVar* parMean = (RooRealVar*) r->floatParsFinal().find(("mu"+suffixName).c_str()); + RooRealVar* parSigma = (RooRealVar*) r->floatParsFinal().find(("sigma"+suffixName).c_str()); + meanVal_postFit[ij] = parMean->getValV(); + sigmaVal_postFit[ij] = parSigma->getValV(); + + } + + // print AxE results + for(auto ij=0; ijLoad("../HiggsAnalysis/CombinedLimit/lib/libHiggsAnalysisCombinedLimit.so"); + gROOT->Macro("/afs/cern.ch/user/a/amartell/public/setStyle.C"); + + RooWorkspace w("w"); + RooRealVar x("x", "", 4.5, 6.); + + + //read workspace for signal model + TFile* inF1 = TFile::Open("signalModel_MC_isResonant1_isEE1_workspace.root"); + RooWorkspace* wSignal = (RooWorkspace*)inF1->Get("signalModel"); //->Clone("wSignal"); + + //non resonant for AxE + TFile* inF1nnR = TFile::Open("signalModel_MC_isResonant0_isEE1_workspace.root"); + RooWorkspace* wSignalnnR = (RooWorkspace*)inF1nnR->Get("signalModel"); //->Clone("wSignal"); + + //read workspace for partially reco bkg + TFile* inF2 = TFile::Open("part_workspace.root"); + RooWorkspace* wPartial = (RooWorkspace*)inF2->Get("myPartialWorkSpace"); //->Clone("wPartial"); + + wSignal->Print(); + wPartial->Print(); + + + float BmassFit = wSignal->var("mu_mll_2.4-4.0")->getVal(); + float BsigmaFit = wSignal->var("sigma_mll_2.4-4.0")->getVal(); + + float Blow_sideband = BmassFit - 3.*BsigmaFit; + float Bhigh_sideband = BmassFit + 3.*BsigmaFit; + + float Blow_blindmin = x.getBinning().lowBound(); + float Bhigh_blindmax = x.getBinning().highBound(); + + std::cout << " Blow_blindmin = " << Blow_blindmin << " Bhigh_blindmax = " << Bhigh_blindmax << std::endl; + + //load data + std::string inputFileList = "/afs/cern.ch/user/k/klau/myWorkspace/public/ForArabella/03Mar2020_pf/"; + if(isEE) + inputFileList += "RootTree_2020Jan16_Run2018ABCD_BToKEEAnalyzer_2020Feb18_fullq2_pf_isoPFMVADphiptImb_weighted_pauc02_mvaCut02.root"; + + std::cout << "inputFileList = " << inputFileList << std::endl; + + TChain* tree = new TChain("tree"); + tree->Add(inputFileList.c_str()); + + int nEntriesTot = tree->GetEntries(); + std::cout << " loaded nEvents = " << nEntriesTot << std::endl; + + float B_mll_fullfit; + float B_fit_mass; + float B_bdt_weight; + + tree->SetBranchStatus("*", 0); + if(isEE){ + tree->SetBranchStatus("BToKEE_mll_fullfit", 1); tree->SetBranchAddress("BToKEE_mll_fullfit", &B_mll_fullfit); + tree->SetBranchStatus("BToKEE_fit_mass", 1); tree->SetBranchAddress("BToKEE_fit_mass", &B_fit_mass); + tree->SetBranchStatus("BToKEE_xgb", 1); tree->SetBranchAddress("BToKEE_xgb", &B_bdt_weight); + } + else{ + tree->SetBranchStatus("BToKMuMu_mll_fullfit", 1); tree->SetBranchAddress("BToKMuMu_mll_fullfit", &B_mll_fullfit); + tree->SetBranchStatus("BToKMuMu_fit_mass", 1); tree->SetBranchAddress("BToKMuMu_fit_mass", &B_fit_mass); + tree->SetBranchStatus("BToKMuMu_xgb", 1); tree->SetBranchAddress("BToKMuMu_xgb", &B_bdt_weight); + } + + //for the moment JPsi and lowq2 + std::vector> mll_bins_Range; + mll_bins_Range.push_back(std::pair(1.1, 2.4)); + mll_bins_Range.push_back(std::pair(2.4, 4.)); + + RooDataSet* data[2]; + int ijC = 0; + for(auto ij: mll_bins_Range){ + data[ijC] = new RooDataSet(Form("data_mll_%.1f-%.1f", ij.first, ij.second), + Form("data_mll_%.1f-%.1f", ij.first, ij.second), RooArgSet(x)); + ++ijC; + } + for(int ij=0; ijGetEntry(ij); + + if(B_bdt_weight <= mvaCut) continue; + int ijC = 0; + for(auto ijB: mll_bins_Range){ + if(B_mll_fullfit >= ijB.first && B_mll_fullfit < ijB.second) { + //blind lowq2 + if(ijC == 0 && + B_fit_mass >= Blow_sideband && B_fit_mass <= Bhigh_sideband) continue; + x = B_fit_mass; + data[ijC]->add(RooArgSet(x)); + } + ++ijC; + } + } + w.import(x); + + + //Resonant category: background model + signal model + w.factory("nsigR[5.e3, 0., 1.e4]"); + w.factory("nbkgR[1.e4, 0., 1.e5]"); + + w.factory("Exponential::combBkgR(x, tauR[-3.0, -100.0, -1.0e-5])"); + w.import("part_workspace.root:myPartialWorkSpace:partial", RenameAllNodes("Bkg")); + + RooRealVar smodelR_mu("smodelR_mu", "", (wSignal->var("mu_mll_2.4-4.0"))->getVal()); + RooRealVar smodelR_sigma("smodelR_sigma", "", (wSignal->var("sigma_mll_2.4-4.0"))->getVal()); + RooRealVar smodelR_aL("smodelR_aL", "", (wSignal->var("aL_mll_2.4-4.0"))->getVal()); + RooRealVar smodelR_nL("smodelR_nL", "", (wSignal->var("nL_mll_2.4-4.0"))->getVal()); + RooRealVar smodelR_aR("smodelR_aR", "", (wSignal->var("aR_mll_2.4-4.0"))->getVal()); + RooRealVar smodelR_nR("smodelR_nR", "", (wSignal->var("nR_mll_2.4-4.0"))->getVal()); + w.import(smodelR_mu); + w.import(smodelR_sigma); + w.import(smodelR_aL); + w.import(smodelR_nL); + w.import(smodelR_aR); + w.import(smodelR_nR); + w.factory("RooDoubleCBFast::smodelR(x, smodelR_mu, smodelR_sigma, smodelR_aL, smodelR_nL, smodelR_aR, smodelR_nR)"); + + w.factory("SUM::modelBkg(f1[0.5, 0., 10.] * partial_Bkg, combBkgR)"); + w.factory("SUM::modelResonant(nsigR * smodelR, nbkgR * modelBkg)"); + + RooAbsPdf * modelResonant = w.pdf("modelResonant"); + + + //fit Bmass in resonant category + RooFitResult* rR = modelResonant->fitTo(*(data[1]), Extended(true), Minimizer("Minuit2"),Save(true)); + + RooRealVar* parNSignalResonant = (RooRealVar*)rR->floatParsFinal().find("nsigR"); + + + w.var("x")->setRange("signalRegion", Blow_sideband, Bhigh_sideband); + RooAbsReal* backgroundFraction = w.pdf("modelBkg")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + RooAbsReal* signalFraction = w.pdf("smodelR")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + + float nBkgInsignalRegion = backgroundFraction->getVal() * w.var("nbkgR")->getVal(); + float nSignalInsignalRegion = signalFraction->getVal() * w.var("nsigR")->getVal(); + + std::cout << " resonant q2: nSignal = " << nSignalInsignalRegion << " nBkg = " << nBkgInsignalRegion + << " S/sqrt(S+B) = " << nSignalInsignalRegion / sqrt(nSignalInsignalRegion+nBkgInsignalRegion) << std::endl; + std::cout << "\n\n " << std::endl; + + + RooPlot* plotR = w.var("x")->frame(); + plotR->SetXTitle("K(JPsi)ee mass (GeV)"); + plotR->SetTitle(""); + data[1]->plotOn(plotR, Binning(50)); + modelResonant->plotOn(plotR); + modelResonant->plotOn(plotR, Name("fitBkg"), Components("modelBkg"),LineStyle(kDashed), LineColor(kBlue)); + modelResonant->plotOn(plotR, Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), DrawOption("F"), MoveToBack()); + modelResonant->plotOn(plotR, Name("fitPartialBkg"), Components("partial_Bkg"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); + modelResonant->plotOn(plotR, Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), DrawOption("F"), MoveToBack()); + modelResonant->plotOn(plotR, Name("fitSig"), Components("smodelR"),LineColor(kRed+1)); + data[1]->plotOn(plotR, Binning(50)); + + TCanvas* ccR = new TCanvas(); + ccR->SetLogy(0); + plotR->Draw(); + TPaveText pt(0.7,0.70,0.9,0.90,"nbNDC"); + pt.SetFillColor(0); + pt.SetLineWidth(0); + pt.SetBorderSize(1); + pt.SetTextFont(42); + pt.SetTextSize(0.04); + pt.SetTextAlign(12); + pt.AddText(Form("MVA cut = %.2f", mvaCut)); + pt.AddText(Form("S = %.0f #pm %.f", nSignalInsignalRegion, w.var("nsigR")->getError())); + pt.AddText(Form("B = %.0f #pm %.f", nBkgInsignalRegion, w.var("nbkgR")->getError())); + pt.AddText(Form("S/#sqrt{S+B} = %.1f", nSignalInsignalRegion / sqrt(nSignalInsignalRegion+nBkgInsignalRegion))); + pt.Draw("same"); + ccR->Update(); + + TLegend tl(0.70,0.5,0.90,0.70); + tl.AddEntry(plotR->findObject("fitPartialBkg"), "partial","l"); + tl.AddEntry(plotR->findObject("fitBkg"), "all bkg","l"); + tl.AddEntry(plotR->findObject("fitSig"), "signal","l"); + tl.SetTextFont(42); + tl.SetTextSize(0.04); + tl.Draw("same"); + + ccR->Print("modelResonant.png"); + ccR->Print("modelResonant.pdf"); + + //Lowq2 + w.factory("nbkgLowq2[500, 0., 1.e5]"); + + w.factory("Exponential::combBkgL(x, tauL[-3.0, -100.0, -1.0e-5])"); + w.import("part_workspace.root:myPartialWorkSpace:partial", RenameAllNodes("BkgL")); + + w.factory("SUM::modelBkgL(f2[0.5,0.,1.] * partial_BkgL, combBkgL)"); + w.factory("SUM::modelBkgLowq2(nbkgLowq2 * modelBkgL)"); + + RooAbsPdf * modelBkgLowq2 = w.pdf("modelBkgLowq2"); + + + //fit Bmass in lowq2 category + w.var("x")->setRange("Slow", Blow_blindmin, Blow_sideband); + w.var("x")->setRange("Shigh", Bhigh_sideband, Bhigh_blindmax); + + RooFitResult* rL = modelBkgLowq2->fitTo(*(data[0]), Extended(true), Minimizer("Minuit2"), Save(true), Range("Slow,Shigh")); + + RooAbsReal* backgroundFractionL = w.pdf("modelBkgLowq2")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + float nBkgInsignalRegionL = backgroundFractionL->getVal() * w.var("nbkgLowq2")->getVal(); + + RooPlot* plotL = w.var("x")->frame(); + plotL->SetXTitle("Kee mass (GeV)"); + plotL->SetTitle(""); + data[0]->plotOn(plotL, Binning(50)); + modelBkgLowq2->plotOn(plotL); + modelBkgLowq2->plotOn(plotL, Name("fitBkg"), Components("modelBkgLowq2"),LineStyle(kDashed), LineColor(kBlue)); + modelBkgLowq2->plotOn(plotL, Name("fitPartialBkgL"), Components("partial_BkgL"), FillColor(40), LineColor(40), DrawOption("F")); + data[0]->plotOn(plotL, Binning(50)); + + + TCanvas* ccL = new TCanvas(); + ccL->SetLogy(0); + plotL->Draw(); + TPaveText ptL(0.7,0.70,0.9,0.90,"NBNDC"); + ptL.SetFillColor(0); + ptL.SetLineWidth(0); + ptL.SetBorderSize(1); + ptL.SetTextFont(42); + ptL.SetTextSize(0.04); + ptL.SetTextAlign(12); + ptL.AddText(Form("MVA cut = %.2f", mvaCut)); + ptL.AddText(Form("B = %.0f #pm %.f", nBkgInsignalRegionL, w.var("nbkgLowq2")->getError())); + ptL.Draw("same"); + ccL->Update(); + + TLegend tlL(0.70,0.5,0.90,0.70); + tlL.AddEntry(plotL->findObject("fitPartialBkgL"), "partial","l"); + tlL.AddEntry(plotL->findObject("fitBkg"), "all bkg","l"); + tlL.SetTextFont(42); + tlL.SetTextSize(0.04); + tlL.Draw("same"); + + ccL->Print("modelBkgLowq2.png"); + ccL->Print("modelBkgLowq2.pdf"); + + + //now build the full lowq2 model with toys: dataset from combined signal(reascaled from resonant) + bkg + RooRealVar BR_Kll ("BR_Kll", "", 4.51e-7); + RooRealVar BR_KJPsill ("BR_KJPsill", "", 1.01e-3 * 0.0597); + + RooRealVar AxE_KJPsill_mll_2p4_4 ("AxE_KJPsill_mll_2p4_4", "", wSignal->var("axeVal_mll_2.4-4.0")->getValV()); + RooRealVar AxE_KJPsill_mll_1p1_4 ("AxE_KJPsill_mll_1p1_4", "", wSignal->var("axeVal_mll_1.1-4.0")->getValV()); + RooRealVar AxE_Kll_mll_1p1_2p4 ("AxE_Kll_mll_1p1_2p4", "", wSignalnnR->var("axeVal_mll_1.1-2.4")->getValV()); + RooRealVar AxE_Kll_mll_1p1_4 ("AxE_Kll_mll_1p1_4", "", wSignalnnR->var("axeVal_mll_1.1-4.0")->getValV()); + + // RooFormulaVar nsigLowq2("nsigLowq2", "", "@0 * @1/@2 * @3/@4 / (@5/@6) ", RooArgList(*parNSignalResonant, BR_Kll, BR_KJPsill, AxE_Kll_mll_1p1_2p4, AxE_Kll_mll_1p1_4, AxE_KJPsill_mll_2p4_4, AxE_KJPsill_mll_1p1_4)); + + RooFormulaVar nsigLowq2("nsigLowq2", "", "@0 * @1/@2 * @3/@4 ", RooArgList(*parNSignalResonant, BR_Kll, BR_KJPsill, AxE_Kll_mll_1p1_2p4, AxE_KJPsill_mll_2p4_4)); + + std::cout << " nResonant = " << parNSignalResonant->getVal() << " AxE = " << AxE_Kll_mll_1p1_2p4.getVal()/AxE_KJPsill_mll_2p4_4.getVal() + << " BR ratio = " << BR_Kll.getVal() / BR_KJPsill.getVal() << " resonant to nn resonant = " << nsigLowq2.getVal() << std::endl; + std::cout << "\n " << std::endl; + w.import(nsigLowq2); + + + w.factory("RooDoubleCBFast::smodelL(x, smodelR_mu, smodelR_sigma, smodelR_aL, smodelR_nL, smodelR_aR, smodelR_nR)"); + w.factory("SUM::modelToy(nsigLowq2 * smodelL, nbkgLowq2 * modelBkgLowq2)"); + RooAbsPdf* modelToy = w.pdf("modelToy"); + + //generate toy roodataset + RooDataSet* toyDataRaw = modelToy->generate(RooArgSet(*(w.var("x"))), 100 * w.var("nbkgLowq2")->getVal()); + //weight for events(want same number of bkg events but reducing the stat fluctuations) + RooFormulaVar evtWeight("evtWeight", "", " 0.01", *(w.var("x"))); + RooRealVar* evtWeightVal = (RooRealVar*) toyDataRaw->addColumn(evtWeight) ; + RooDataSet toyData(toyDataRaw->GetName(), toyDataRaw->GetTitle(), toyDataRaw, *toyDataRaw->get(), 0, evtWeightVal->GetName()); + + RooPlot* plotT = w.var("x")->frame(); + plotT->SetXTitle("Kee mass (GeV) - toy"); + plotT->SetTitle(""); + toyData.plotOn(plotT, Binning(50)); + + modelToy->fitTo(toyData, Extended(true), Minimizer("Minuit2"), Save(true), SumW2Error(kTRUE)); + modelToy->plotOn(plotT); + + + TCanvas* ccT = new TCanvas(); + ccT->SetLogy(0); + plotT->Draw(); + ccT->Print("toyData.png"); + ccT->Print("toyData.pdf"); + + + std::cout << " fit toy nbkgLowq2 = " << w.var("nbkgLowq2")->getVal() << " +/- " << w.var("nbkgLowq2")->getError() << std::endl; + std::cout << "\n\n " << std::endl; + + //fit simultaneously resonant (data) and non resonant (toy) to have correct estimate of errors + //Define categories + w.factory("nsigL[5.e3, 0., 1.e4]"); + w.factory("nbkgL[1.e4, 0., 1.e5]"); + w.factory("SUM::modelLowq2(nsigL * smodelL, nbkgL * modelBkgL)"); + RooAbsPdf* modelLowq2 = w.pdf("modelLowq2"); + + + RooCategory sample("sample", "sample"); + sample.defineType("lowq2"); + sample.defineType("resonant"); + + //Construct combined dataset in (x, sample) with weights + RooRealVar combW("combW", "", 0, 1); + + RooArgSet varsPlusCat(*(w.var("x"))); varsPlusCat.add(sample); + RooArgSet varsPlusWeight(varsPlusCat); varsPlusWeight.add(combW); + RooDataSet combData("combData", "", varsPlusWeight, WeightVar(combW)); + + sample.setLabel("lowq2"); + for(int i=0, n=toyData.numEntries(); inumEntries(); iget(i); + combData.add(varsPlusCat, data[1]->weight()); + } + + //combined dataset without weights + //RooDataSet combData("combData", "", *(w.var("x")), RooFit::Index(sample), RooFit::Import("lowq2", toyDataW), RooFit::Import("resonant", resonantData)); + + //Construct a simultaneous pdf in (x, sample) + //with category sample as index + RooSimultaneous simPdf("simPdf", "", sample); + + //Associate each model to the corresponding sample + simPdf.addPdf(*(w.pdf("modelLowq2")), "lowq2"); + simPdf.addPdf(*(w.pdf("modelResonant")), "resonant"); + + RooFitResult * r = simPdf.fitTo(combData, Extended(true), Minimizer("Minuit2"),Save(true), SumW2Error(kTRUE)); + + RooAbsReal* comb_backgroundFractionR = w.pdf("modelBkg")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + RooAbsReal* comb_signalFractionR = w.pdf("smodelR")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + RooAbsReal* comb_backgroundFractionL = w.pdf("modelBkgL")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + RooAbsReal* comb_signalFractionL = w.pdf("smodelL")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); + + float comb_nBkgInsignalRegionR = comb_backgroundFractionR->getVal() * w.var("nbkgR")->getVal(); + float comb_nSignalInsignalRegionR = comb_signalFractionR->getVal() * w.var("nsigR")->getVal(); + float comb_nBkgInsignalRegionL = comb_backgroundFractionR->getVal() * w.var("nbkgL")->getVal(); + float comb_nSignalInsignalRegionL = comb_signalFractionR->getVal() * w.var("nsigL")->getVal(); + + std::cout << " resonant q2: nSignal = " << comb_nSignalInsignalRegionR << " nBkg = " << comb_nBkgInsignalRegionR + << " S/sqrt(S+B) = " << comb_nSignalInsignalRegionR / sqrt(comb_nSignalInsignalRegionR+comb_nBkgInsignalRegionR) << std::endl; + std::cout << " low q2: nSignal = " << comb_nSignalInsignalRegionL << " nBkg = " << comb_nBkgInsignalRegionL + << " S/sqrt(S+B) = " << comb_nSignalInsignalRegionL / sqrt(comb_nSignalInsignalRegionL+comb_nBkgInsignalRegionL) << std::endl; + std::cout << "\n\n " << std::endl; + + + + RooPlot* frame1 = w.var("x")->frame(); + frame1->SetXTitle("Kee mass (GeV) - toy"); + frame1->SetTitle(""); + RooPlot* frame2 = w.var("x")->frame(); + frame2->SetXTitle("K(JPsi)ee mass (GeV)"); + frame2->SetTitle(""); + + combData.plotOn(frame1, Binning(50), RooFit::Cut("sample==sample::lowq2")); + combData.plotOn(frame2, Binning(50), RooFit::Cut("sample==sample::resonant")); + + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData)); + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitBkg"), Components("modelBkgL"), LineStyle(kDashed), LineColor(kBlue)); + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgL"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitPartialBkg"), Components("partial_BkgL"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgL"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitSig"), Components("smodelL"), LineColor(kRed+1)); + combData.plotOn(frame1, Binning(50), RooFit::Cut("sample==sample::lowq2")); + + + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData)); + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitBkg"), Components("modelBkg"), LineStyle(kDashed), LineColor(kBlue)); + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitPartialBkg"), Components("partial_Bkg"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitSig"), Components("smodelR"), LineColor(kRed+1)); + combData.plotOn(frame2, Binning(50), RooFit::Cut("sample==sample::resonant")); + + + + TCanvas* tc = new TCanvas("tc", "", 1300, 600); + tc->Divide(2); + tc->cd(1); + //gPad->SetLeftMargin(0.15); + // frame1->GetYaxis()->SetTitleOffset(1.4); + frame1->Draw(); + + TPaveText ptCombL(0.65,0.70,0.85,0.90,"nbNDC"); + ptCombL.SetFillColor(0); + ptCombL.SetLineWidth(0); + ptCombL.SetBorderSize(1); + ptCombL.SetTextFont(42); + ptCombL.SetTextSize(0.04); + ptCombL.SetTextAlign(12); + ptCombL.AddText(Form("MVA cut = %.2f", mvaCut)); + ptCombL.AddText(Form("S = %.0f #pm %.f", comb_nSignalInsignalRegionL, w.var("nsigL")->getError())); + ptCombL.AddText(Form("B = %.0f #pm %.f", comb_nBkgInsignalRegionL, w.var("nbkgL")->getError())); + ptCombL.AddText(Form("S/#sqrt{S+B} = %.1f", comb_nSignalInsignalRegionL / sqrt(comb_nSignalInsignalRegionL+comb_nBkgInsignalRegionL))); + ptCombL.Draw("same"); + tc->Update(); + TLegend tlCombL(0.70,0.5,0.90,0.70); + tlCombL.AddEntry(frame1->findObject("combBkg"), "combinatorial","l"); + tlCombL.AddEntry(frame1->findObject("fitPartialBkg"), "partial","l"); + tlCombL.AddEntry(frame1->findObject("fitBkg"), "all bkg","l"); + tlCombL.AddEntry(frame1->findObject("fitSig"), "signal","l"); + tlCombL.SetTextFont(42); + tlCombL.SetTextSize(0.04); + tlCombL.Draw("same"); + + + tc->cd(2); + //gPad->SetLeftMargin(0.15); + // frame2->GetYaxis()->SetTitleOffset(1.4); + frame2->Draw(); + + TPaveText ptCombR(0.65,0.70,0.85,0.90,"nbNDC"); + ptCombR.SetFillColor(0); + ptCombR.SetLineWidth(0); + ptCombR.SetBorderSize(1); + ptCombR.SetTextFont(42); + ptCombR.SetTextSize(0.04); + ptCombR.SetTextAlign(12); + ptCombR.AddText(Form("MVA cut = %.2f", mvaCut)); + ptCombR.AddText(Form("S = %.0f #pm %.f", comb_nSignalInsignalRegionR, w.var("nsigR")->getError())); + ptCombR.AddText(Form("B = %.0f #pm %.f", comb_nBkgInsignalRegionR, w.var("nbkgR")->getError())); + ptCombR.AddText(Form("S/#sqrt{S+B} = %.1f", comb_nSignalInsignalRegionR / sqrt(comb_nSignalInsignalRegionR+comb_nBkgInsignalRegionR))); + ptCombR.Draw("same"); + tc->Update(); + TLegend tlCombR(0.70,0.5,0.90,0.70); + tlCombR.AddEntry(frame2->findObject("combBkg"), "combinatorial","l"); + tlCombR.AddEntry(frame2->findObject("fitPartialBkg"), "partial","l"); + tlCombR.AddEntry(frame2->findObject("fitBkg"), "all bkg","l"); + tlCombR.AddEntry(frame2->findObject("fitSig"), "signal","l"); + tlCombR.SetTextFont(42); + tlCombR.SetTextSize(0.04); + tlCombR.Draw("same"); + + tc->Print("combined.png"); + tc->Print("combined.pdf"); + + +} diff --git a/models.cc b/unused/models.cc similarity index 100% rename from models.cc rename to unused/models.cc diff --git a/models.h b/unused/models.h similarity index 100% rename from models.h rename to unused/models.h From 1e67ecc374c0f572c4e17d98fc3f5e7807bcfcd1 Mon Sep 17 00:00:00 2001 From: arabella Date: Mon, 9 Mar 2020 11:36:13 +0100 Subject: [PATCH 2/3] fix typo in signalRegion --- makeSimultaneousFits.C | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/makeSimultaneousFits.C b/makeSimultaneousFits.C index 7db9fc9..b701430 100644 --- a/makeSimultaneousFits.C +++ b/makeSimultaneousFits.C @@ -152,10 +152,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ //fit Bmass in resonant category RooFitResult* rR = modelResonant->fitTo(*(data[1]), Extended(true), Minimizer("Minuit2"),Save(true)); - RooRealVar* parNSignalResonant = (RooRealVar*)rR->floatParsFinal().find("nsigR"); - - - w.var("x")->setRange("signalRegion", Blow_sideband, Bhigh_sideband); + w.var("x")->setRange("signalRegioncut", Blow_sideband, Bhigh_sideband); RooAbsReal* backgroundFraction = w.pdf("modelBkg")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); RooAbsReal* signalFraction = w.pdf("smodelR")->createIntegral(*(w.var("x")), RooFit::NormSet(*(w.var("x"))), Range("signalRegioncut")); @@ -265,6 +262,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ //now build the full lowq2 model with toys: dataset from combined signal(reascaled from resonant) + bkg + RooRealVar nResonantEvents ("nResonantEvents", "", nSignalInsignalRegion); RooRealVar BR_Kll ("BR_Kll", "", 4.51e-7); RooRealVar BR_KJPsill ("BR_KJPsill", "", 1.01e-3 * 0.0597); @@ -275,9 +273,9 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ // RooFormulaVar nsigLowq2("nsigLowq2", "", "@0 * @1/@2 * @3/@4 / (@5/@6) ", RooArgList(*parNSignalResonant, BR_Kll, BR_KJPsill, AxE_Kll_mll_1p1_2p4, AxE_Kll_mll_1p1_4, AxE_KJPsill_mll_2p4_4, AxE_KJPsill_mll_1p1_4)); - RooFormulaVar nsigLowq2("nsigLowq2", "", "@0 * @1/@2 * @3/@4 ", RooArgList(*parNSignalResonant, BR_Kll, BR_KJPsill, AxE_Kll_mll_1p1_2p4, AxE_KJPsill_mll_2p4_4)); + RooFormulaVar nsigLowq2("nsigLowq2", "", "@0 * @1/@2 * @3/@4 ", RooArgList(nResonantEvents, BR_Kll, BR_KJPsill, AxE_Kll_mll_1p1_2p4, AxE_KJPsill_mll_2p4_4)); - std::cout << " nResonant = " << parNSignalResonant->getVal() << " AxE = " << AxE_Kll_mll_1p1_2p4.getVal()/AxE_KJPsill_mll_2p4_4.getVal() + std::cout << " nResonant = " << nResonantEvents.getVal() << " AxE = " << AxE_Kll_mll_1p1_2p4.getVal()/AxE_KJPsill_mll_2p4_4.getVal() << " BR ratio = " << BR_Kll.getVal() / BR_KJPsill.getVal() << " resonant to nn resonant = " << nsigLowq2.getVal() << std::endl; std::cout << "\n " << std::endl; w.import(nsigLowq2); From 4fa5d5b13db03db47c49671b628ae34ec8537f9a Mon Sep 17 00:00:00 2001 From: arabella Date: Tue, 10 Mar 2020 18:17:07 +0100 Subject: [PATCH 3/3] partial reco bkg WP dependent plus fixes --- AxE_MC_withFits_root.C | 4 +- README.md | 6 ++ makeSimultaneousFits.C | 27 ++----- modelFrom_partialRecoBkg.C | 157 +++++++++++++++++++++++++++++++++++++ 4 files changed, 173 insertions(+), 21 deletions(-) create mode 100644 modelFrom_partialRecoBkg.C diff --git a/AxE_MC_withFits_root.C b/AxE_MC_withFits_root.C index 7b0665f..70b364b 100644 --- a/AxE_MC_withFits_root.C +++ b/AxE_MC_withFits_root.C @@ -88,12 +88,12 @@ void AxE_MC_withFits_root(int isMC=1, int isEE=1, int isResonant=0, float mvaCut //loop over tree load hostos and RooDataSet TH1F* hHisto[3]; RooWorkspace w("signalModel"); - RooRealVar x("x", "", 4.4, 6.); + RooRealVar x("x", "", 4.5, 6.); RooDataSet* data[3]; int ijC = 0; for(auto ij: mll_bins_Range){ - hHisto[ijC] = new TH1F(Form("hHisto_mll_%.1f-%.1f", ij.first, ij.second), "", 500, 4.4, 6.); + hHisto[ijC] = new TH1F(Form("hHisto_mll_%.1f-%.1f", ij.first, ij.second), "", 500, 4.5, 6.); data[ijC] = new RooDataSet(Form("data_mll_%.1f-%.1f", ij.first, ij.second), Form("data_mll_%.1f-%.1f", ij.first, ij.second), RooArgSet(x)); diff --git a/README.md b/README.md index b067ffe..cc7e3b0 100644 --- a/README.md +++ b/README.md @@ -39,6 +39,12 @@ Need to adapt format for input files to other's ntuples ```shell python makePlot_fitPeak_unbinned.py -p -i /afs/cern.ch/user/k/klau/myWorkspace/public/ForArabella/03Mar2020_pf/RootTree_2020Jan16_BdToKstarJpsi_ToKPiee_BToKEEAnalyzer_2020Feb18_fullq2_pf_isoPFMVADphiptImb_weighted_pauc02_mva.root -o part_workspace.root ``` +Model WP dependent, q2 bin dependent and for both resonant and non-resonant +```shell +root modelFrom_partialRecoBkg.C +``` + + * run the simultaneous fit on the resonant (data) and non-resonant (toy MC) to extract expected S/sqrt(S+B) with error - import models for signal (MC resonant) and partially reconstructed background. Build resonant_PDF (signal + partialRecoBkg + expo for combinatorial bkg) diff --git a/makeSimultaneousFits.C b/makeSimultaneousFits.C index b701430..85003c7 100644 --- a/makeSimultaneousFits.C +++ b/makeSimultaneousFits.C @@ -43,12 +43,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ TFile* inF1nnR = TFile::Open("signalModel_MC_isResonant0_isEE1_workspace.root"); RooWorkspace* wSignalnnR = (RooWorkspace*)inF1nnR->Get("signalModel"); //->Clone("wSignal"); - //read workspace for partially reco bkg - TFile* inF2 = TFile::Open("part_workspace.root"); - RooWorkspace* wPartial = (RooWorkspace*)inF2->Get("myPartialWorkSpace"); //->Clone("wPartial"); - wSignal->Print(); - wPartial->Print(); float BmassFit = wSignal->var("mu_mll_2.4-4.0")->getVal(); @@ -127,7 +122,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ w.factory("nbkgR[1.e4, 0., 1.e5]"); w.factory("Exponential::combBkgR(x, tauR[-3.0, -100.0, -1.0e-5])"); - w.import("part_workspace.root:myPartialWorkSpace:partial", RenameAllNodes("Bkg")); + w.import("partiallyRecoBkgModel_MC_isResonant1_isEE1_workspace.root:partialModel:partial_mll_2.4-4.0", RenameAllNodes("BkgR")); RooRealVar smodelR_mu("smodelR_mu", "", (wSignal->var("mu_mll_2.4-4.0"))->getVal()); RooRealVar smodelR_sigma("smodelR_sigma", "", (wSignal->var("sigma_mll_2.4-4.0"))->getVal()); @@ -143,7 +138,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ w.import(smodelR_nR); w.factory("RooDoubleCBFast::smodelR(x, smodelR_mu, smodelR_sigma, smodelR_aL, smodelR_nL, smodelR_aR, smodelR_nR)"); - w.factory("SUM::modelBkg(f1[0.5, 0., 10.] * partial_Bkg, combBkgR)"); + w.factory("SUM::modelBkg(f1[0.5, 0., 10.] * partial_mll_2.4-4.0_BkgR, combBkgR)"); w.factory("SUM::modelResonant(nsigR * smodelR, nbkgR * modelBkg)"); RooAbsPdf * modelResonant = w.pdf("modelResonant"); @@ -171,7 +166,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ modelResonant->plotOn(plotR); modelResonant->plotOn(plotR, Name("fitBkg"), Components("modelBkg"),LineStyle(kDashed), LineColor(kBlue)); modelResonant->plotOn(plotR, Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), DrawOption("F"), MoveToBack()); - modelResonant->plotOn(plotR, Name("fitPartialBkg"), Components("partial_Bkg"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); + modelResonant->plotOn(plotR, Name("fitPartialBkg"), Components("partial_mll_2.4-4.0_BkgR"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); modelResonant->plotOn(plotR, Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), DrawOption("F"), MoveToBack()); modelResonant->plotOn(plotR, Name("fitSig"), Components("smodelR"),LineColor(kRed+1)); data[1]->plotOn(plotR, Binning(50)); @@ -208,9 +203,9 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ w.factory("nbkgLowq2[500, 0., 1.e5]"); w.factory("Exponential::combBkgL(x, tauL[-3.0, -100.0, -1.0e-5])"); - w.import("part_workspace.root:myPartialWorkSpace:partial", RenameAllNodes("BkgL")); + w.import("partiallyRecoBkgModel_MC_isResonant0_isEE1_workspace.root:partialModel:partial_mll_1.1-2.4", RenameAllNodes("BkgL")); - w.factory("SUM::modelBkgL(f2[0.5,0.,1.] * partial_BkgL, combBkgL)"); + w.factory("SUM::modelBkgL(f2[0.5,0.,1.] * partial_mll_1.1-2.4_BkgL, combBkgL)"); w.factory("SUM::modelBkgLowq2(nbkgLowq2 * modelBkgL)"); RooAbsPdf * modelBkgLowq2 = w.pdf("modelBkgLowq2"); @@ -230,8 +225,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ plotL->SetTitle(""); data[0]->plotOn(plotL, Binning(50)); modelBkgLowq2->plotOn(plotL); - modelBkgLowq2->plotOn(plotL, Name("fitBkg"), Components("modelBkgLowq2"),LineStyle(kDashed), LineColor(kBlue)); - modelBkgLowq2->plotOn(plotL, Name("fitPartialBkgL"), Components("partial_BkgL"), FillColor(40), LineColor(40), DrawOption("F")); + modelBkgLowq2->plotOn(plotL, Name("fitBkg"), Components("modelBkgLowq2"),LineStyle(kDashed), LineColor(kBlue), Range("Slow,Shigh")); data[0]->plotOn(plotL, Binning(50)); @@ -251,7 +245,6 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ ccL->Update(); TLegend tlL(0.70,0.5,0.90,0.70); - tlL.AddEntry(plotL->findObject("fitPartialBkgL"), "partial","l"); tlL.AddEntry(plotL->findObject("fitBkg"), "all bkg","l"); tlL.SetTextFont(42); tlL.SetTextSize(0.04); @@ -385,7 +378,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData)); simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitBkg"), Components("modelBkgL"), LineStyle(kDashed), LineColor(kBlue)); simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgL"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); - simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitPartialBkg"), Components("partial_BkgL"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); + simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitPartialBkg"), Components("partial_mll_1.1-2.4_BkgL"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgL"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); simPdf.plotOn(frame1, RooFit::Slice(sample, "lowq2"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitSig"), Components("smodelL"), LineColor(kRed+1)); combData.plotOn(frame1, Binning(50), RooFit::Cut("sample==sample::lowq2")); @@ -394,7 +387,7 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData)); simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitBkg"), Components("modelBkg"), LineStyle(kDashed), LineColor(kBlue)); simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); - simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitPartialBkg"), Components("partial_Bkg"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); + simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitPartialBkg"), Components("partial_mll_2.4-4.0_BkgR"), FillColor(40), LineColor(40), DrawOption("F"), AddTo("combBkg")); simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("combBkg"), Components("combBkgR"), FillColor(42), LineColor(42), MoveToBack(), DrawOption("F")); simPdf.plotOn(frame2, RooFit::Slice(sample, "resonant"), RooFit::ProjWData(RooArgSet(sample), combData), Name("fitSig"), Components("smodelR"), LineColor(kRed+1)); combData.plotOn(frame2, Binning(50), RooFit::Cut("sample==sample::resonant")); @@ -404,8 +397,6 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ TCanvas* tc = new TCanvas("tc", "", 1300, 600); tc->Divide(2); tc->cd(1); - //gPad->SetLeftMargin(0.15); - // frame1->GetYaxis()->SetTitleOffset(1.4); frame1->Draw(); TPaveText ptCombL(0.65,0.70,0.85,0.90,"nbNDC"); @@ -432,8 +423,6 @@ void makeSimultaneousFits(int isEE=1, float mvaCut = 12.68){ tc->cd(2); - //gPad->SetLeftMargin(0.15); - // frame2->GetYaxis()->SetTitleOffset(1.4); frame2->Draw(); TPaveText ptCombR(0.65,0.70,0.85,0.90,"nbNDC"); diff --git a/modelFrom_partialRecoBkg.C b/modelFrom_partialRecoBkg.C new file mode 100644 index 0000000..db4feb5 --- /dev/null +++ b/modelFrom_partialRecoBkg.C @@ -0,0 +1,157 @@ +//root modelFrom_partialRecoBkg.C'(1, 0)' + + +//roofit +#include "RooRealVar.h" +#include "RooDataSet.h" +#include "RooDataHist.h" +#include "RooGaussian.h" +#include "RooChebychev.h" +#include "RooCBShape.h" +#include "RooArgusBG.h" +#include "RooFFTConvPdf.h" +#include "RooAddPdf.h" +#include "RooWorkspace.h" +#include "RooFitResult.h" +#include "RooFormulaVar.h" +#include "TText.h" +#include "TCanvas.h" +#include "TAxis.h" +#include "TH1F.h" +#include "TH1D.h" +#include "RooPlot.h" + + #include + #include + #include + #include + + +using namespace RooFit; +void modelFrom_partialRecoBkg(int isEE=1, int isResonant=0, float mvaCut = 12.68){ + + gSystem->Load("../HiggsAnalysis/CombinedLimit/lib/libHiggsAnalysisCombinedLimit.so"); + + + int nBins = 3; + std::vector> mll_bins_Range; + mll_bins_Range.push_back(std::pair(1.1, 2.4)); + mll_bins_Range.push_back(std::pair(2.4, 4.)); + mll_bins_Range.push_back(std::pair(1.1, 4.)); + + + //current files from Otto, need to adapt to other's outputs + std::string inputFileList = "/afs/cern.ch/user/k/klau/myWorkspace/public/ForArabella/03Mar2020_pf/"; + if(!isResonant && isEE) + inputFileList += "RootTree_2020Jan16_BdToKstaree_BToKEEAnalyzer_2020Feb18_fullq2_pf_isoPFMVADphiptImb_weighted_pauc02_mva.root"; + else if(isResonant && isEE) + inputFileList += "RootTree_2020Jan16_BdToKstarJpsi_ToKPiee_BToKEEAnalyzer_2020Feb18_fullq2_pf_isoPFMVADphiptImb_weighted_pauc02_mva.root"; + //need to add for muons + + + std::cout << " isEE = " << isEE << " isResonant = " << isResonant << std::endl; + std::cout << "inputFileList = " << inputFileList << std::endl; + + + TChain* tree = new TChain("tree"); + tree->Add(inputFileList.c_str()); + + int nEntriesTot = tree->GetEntries(); + std::cout << " loaded nEvents = " << nEntriesTot << std::endl; + + float B_mll_fullfit; + float B_fit_mass; + float B_bdt_weight; + + tree->SetBranchStatus("*", 0); + if(isEE){ + tree->SetBranchStatus("BToKEE_mll_fullfit", 1); tree->SetBranchAddress("BToKEE_mll_fullfit", &B_mll_fullfit); + tree->SetBranchStatus("BToKEE_fit_mass", 1); tree->SetBranchAddress("BToKEE_fit_mass", &B_fit_mass); + tree->SetBranchStatus("BToKEE_xgb", 1); tree->SetBranchAddress("BToKEE_xgb", &B_bdt_weight); + } + else{ + tree->SetBranchStatus("BToKMuMu_mll_fullfit", 1); tree->SetBranchAddress("BToKMuMu_mll_fullfit", &B_mll_fullfit); + tree->SetBranchStatus("BToKMuMu_fit_mass", 1); tree->SetBranchAddress("BToKMuMu_fit_mass", &B_fit_mass); + tree->SetBranchStatus("BToKMuMu_xgb", 1); tree->SetBranchAddress("BToKMuMu_xgb", &B_bdt_weight); + } + + //loop over tree load RooDataSet + RooWorkspace w("partialModel"); + RooRealVar x("x", "", 4.5, 6.); + RooDataSet* data[3]; + + int ijC = 0; + for(auto ij: mll_bins_Range){ + data[ijC] = new RooDataSet(Form("data_mll_%.1f-%.1f", ij.first, ij.second), + Form("data_mll_%.1f-%.1f", ij.first, ij.second), RooArgSet(x)); + ++ijC; + } + for(int ij=0; ijGetEntry(ij); + + if(B_bdt_weight <= mvaCut) continue; + int ijC = 0; + for(auto ijB: mll_bins_Range){ + if(B_mll_fullfit >= ijB.first && B_mll_fullfit < ijB.second) { + x = B_fit_mass; + data[ijC]->add(RooArgSet(x)); + } + ++ijC; + } + } + w.import(x); + + + for(auto ij=0; ijframe(); + if(isEE){ + if(isResonant) plot->SetXTitle("K*(JPsi)ee mass (GeV)"); + else plot->SetXTitle("Kee mass (GeV)"); + } + plot->SetTitle(""); + data[ij]->plotOn(plot, Binning(50)); + partialModel->plotOn(plot, LineColor(kRed)); //MoveToBack()); + + float chi2 = plot->chiSquare(6); + + + TCanvas * cc = new TCanvas(); + cc->SetLogy(0); + plot->Draw(); + + TLatex tL; + tL.SetNDC(); + tL.SetTextSize(0.04); + tL.SetTextFont(42); + tL.DrawLatex(0.65,0.8, Form("chi2 = %.1f ",chi2)); + TLatex tL1; + tL1.SetNDC(); + tL1.SetTextSize(0.05); + tL1.SetTextFont(42); + tL1.DrawLatex(0.65,0.7, Form("mll in %.1f-%.1f", mll_bins_Range[ij].first, mll_bins_Range[ij].second)); + std::string outName = "plots/PartiallyReco_MC/"; + if(isEE) outName += "Kee"; + else outName += "Kmm"; + outName += Form("_mll_%.1f-%.1f", mll_bins_Range[ij].first, mll_bins_Range[ij].second); + cc->Print((outName+".png").c_str(), "png"); + cc->Print((outName+".pdf").c_str(), "pdf"); + cc->Print((outName+".root").c_str(), "root"); + } + + w.Print(); + TFile outWS(Form("partiallyRecoBkgModel_MC_isResonant%d_isEE%d_workspace.root", isResonant, isEE), "recreate"); + outWS.cd(); + w.Write(); + outWS.Close(); +}