Utilities/test/testSideBandSubtraction.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <TCanvas.h>
0004 #include <TFile.h>
0005 #include <TH1D.h>
0006 #include <TH2D.h>
0007 #include <TLorentzVector.h>
0008 #include <exception>
0009 #include "PhysicsTools/Utilities/interface/SideBandSubtraction.h"
0010 #include "RooGenericPdf.h"
0011 #include "RooPolynomial.h"
0012 #include "RooGaussian.h"
0013 #include "RooProdPdf.h"
0014 #include "RooPlot.h"
0015 #include "RooAddPdf.h"
0016 #include "RooDataSet.h"
0017
0018 using std::cout;
0019 using std::endl;
0020 using std::vector;
0021
0022 using namespace RooFit;
0023 using std::string;
0024 using std::vector;
0025
0026 void print_plot(RooDataSet *dataSet, RooRealVar printVar, string outname, string cut) {
0027   RooPlot *genericFrame = printVar.frame();
0028   dataSet->plotOn(genericFrame, Cut(cut.c_str()));
0029   TCanvas genericCanvas;
0030   genericFrame->Draw();
0031   outname = outname + ".eps";
0032   genericCanvas.SaveAs(outname.c_str());
0033   outname.replace(outname.size() - 3, outname.size(), "gif");
0034   genericCanvas.SaveAs(outname.c_str());
0035 }
0036 int main(void) {
0037   //variables to generate data for
0038   RooRealVar gfrac("gfrac", "fraction of gaussian", 0.5);
0039   RooRealVar mass("mass", "Separation Variable (Mass)", 2, 4, "GeV/c^2");  //separation variable
0040   RooRealVar pt("pt", "pt", 2, 4, "GeV");
0041
0042   //signal peak
0043   RooRealVar mean("mean", "mean of signal", 3.2, 2, 4, "GeV");
0044   RooRealVar sigma("sigma", "width of signal", 0.2, "GeV");
0045   RooGaussian signal("signal", "signal gaussian", pt, mean, sigma);
0046
0047   //constant background
0048   RooPolynomial background("background", "Constant Background", mass, RooArgList());
0049
0050   //interest variable gaussian
0051   RooRealVar pt_sig_mean("pt_sig_mean", "mean of pt signal", 3.1, 2, 4, "GeV");
0052   RooRealVar pt_sig_sigma("pt_sig_sigma", "width of pt signal", 0.2, "GeV");
0053   RooGaussian pt_sig_gauss("pt_sig_gauss", "pt signal gaussian", mass, pt_sig_mean, pt_sig_sigma);
0054
0055   //interest variable background
0056   RooRealVar pt_bkg_mean("pt_bkg_mean", "mean of pt background", 2.5, 2, 4, "GeV/c");
0057   RooRealVar pt_bkg_sigma("pt_bkg_sigma", "width of pt background", 0.2, "GeV");
0058   RooGaussian pt_bkg_gauss("pt_bkg_gauss", "Gaussian of pt background", pt, pt_bkg_mean, pt_bkg_sigma);
0059
0060   //background and model pdf for SBS
0061   RooProdPdf bg_pdf("signal pdf", "background*mass", RooArgSet(background, signal));  //model_pdf
0062   RooProdPdf signal_pdf("signal_pdf", "pt_sig_gauss*pt_bkg_gauss", RooArgSet(pt_sig_gauss, pt_bkg_gauss));
0063
0064   //Add the two
0065   RooAddPdf sum_pdf("sum_pdf", "bg_pdf+signal_pdf", RooArgList(bg_pdf, signal_pdf), gfrac);
0066
0067   RooDataSet *datasum = sum_pdf.generate(RooArgSet(mass, pt), 25000);  //data
0068   print_plot(datasum, pt, "pt_no_cut", "");
0069
0070   //build the sbs object
0071   vector<TH1F *> basehistos(0);
0072   TH1F ptHist("pt", "Pt Hist", 100, 2, 4);
0073   basehistos.push_back(&ptHist);
0074
0075   SideBandSubtract sbs(&sum_pdf, &bg_pdf, datasum, &mass, basehistos, true);
0076   sbs.addSideBandRegion(2, 2.8);
0077   sbs.addSignalRegion(2.8, 3.4);
0078   sbs.addSideBandRegion(3.4, 4.0);
0079   sbs.doGlobalFit();
0080   sbs.printResults();
0081   return 0;
0082 }