Macros/RooFit/GaussianFitOnData.cc

0001 #ifndef __CINT__
0002 #include "RooGlobalFunc.h"
0003 #endif
0004 #include "RooRealVar.h"
0005 #include "RooDataSet.h"
0006 #include "RooGaussModel.h"
0007 #include "RooAddModel.h"
0008 #include "RooTruthModel.h"
0009 #include "RooDecay.h"
0010 #include "RooPlot.h"
0011 #include "RooLandau.h"
0012 #include "RooGaussian.h"
0013 #include "RooNumConvPdf.h"
0014 #include "RooExponential.h"
0015 #include "RooCBShape.h"
0016 #include "TCanvas.h"
0017 #include "TH1.h"
0018 #include "RooGenericPdf.h"
0019 #include "RooAddPdf.h"
0020 #include <sstream>
0021 #include "TFile.h"
0022 #include "RooDataHist.h"
0023 using namespace RooFit;
0024
0025 void GaussianFitOnData()
0026 {
0027   TFile inputFile("0_MuScleFit.root", "READ");
0028   TH1F * inputHisto = (TH1F*)inputFile.Get("hRecBestResAllEvents_Mass");
0029   inputHisto->Rebin(10);
0030
0031   // RooRealVar x("x", "Mass (GeV)", 2., 4.);
0032   // double xMin = 3.3;
0033   // double xMax = 4.4;
0034   double xMin = 2.;
0035   double xMax = 4.;
0036   RooRealVar x("x", "Mass (GeV)", xMin, xMax);
0037
0038   RooDataHist histo("dh", "dh", x, Import(*inputHisto));
0039
0040   // Build the CB
0041   RooRealVar peak("peak","peak", 3.1, 2.8, 3.8);
0042   // RooRealVar peak("peak","peak", 3.695, 3.6, 3.8);
0043   RooRealVar sigma("sigma","sigma", 0.06, 0.001, 0.1);
0044   RooGaussModel signal("signal", "signal", x, peak, sigma);
0045
0046   // Build the exponential background pdf
0047   RooRealVar expCoeff("expCoeff", "exponential coefficient", -1., -10., 10.);
0048   RooExponential background("exponential", "exponential", x, expCoeff);
0049
0050   TFile outputFile("gaussianFitOnData.root", "RECREATE");
0051
0052   // // Fit the background in the sidebands
0053   // x.setRange("sb_lo", 2., 2.6);
0054   // x.setRange("sb_hi", 3.8, 4.4);
0055   // background.fitTo(histo, Range("sb_lo,sb_hi"));
0056
0057   // TCanvas canvasBackground("canvasBackgroundFit", "canvasBackgroundFit", 1000, 800);
0058   // canvasBackground.cd();
0059   // canvasBackground.Draw();
0060   // RooPlot* frameBackground = x.frame();
0061   // histo.plotOn(frameBackground);
0062   // background.plotOn(frameBackground);
0063   // background.paramOn(frameBackground);
0064
0065   // frameBackground->Draw();
0066   // frameBackground->Write();
0067   // canvasBackground.Write();
0068
0069   // canvasBackground.Write();
0070
0071   // // Fix the background parameter
0072   // expCoeff.setConstant(kTRUE);
0073
0074   // Build the model adding the exponential background
0075   RooRealVar fSig("fSig", "signal fraction", 0.4, 0., 1.);
0076   RooAddPdf model("model", "model", RooArgList(signal, background), fSig);
0077
0078   model.fitTo(histo, "hs");
0079
0080   TCanvas canvas("canvas", "canvas", 1000, 800);
0081   canvas.cd();
0082   canvas.Draw();
0083   RooPlot* frame = x.frame();
0084   histo.plotOn(frame);
0085   model.plotOn(frame);
0086   model.paramOn(frame);
0087
0088   model.plotOn(frame, Components("exponential"), LineStyle(kDashed));
0089
0090
0091   int binXmin = inputHisto->FindBin( xMin );
0092   int binXmax = inputHisto->FindBin( xMax );
0093   Double_t fullIntegral = inputHisto->Integral(binXmin, binXmax);
0094   std::cout << "Events in fit interval = " << fullIntegral << std::endl;
0095
0096   RooAbsReal * backgroundFullIntegral = background.createIntegral(x);
0097   Double_t backgroundFullIntegralValue = backgroundFullIntegral->getVal();
0098   std::cout << "backgroundFullIntegralValue = " << backgroundFullIntegralValue << std::endl;
0099
0100   // Compute the integral and write the signal and background events in +-2.5sigma
0101   Double_t sbMin = peak.getVal() - 2.5*sigma.getVal();
0102   Double_t sbMax = peak.getVal() + 2.5*sigma.getVal();
0103   int minBin = inputHisto->FindBin( sbMin );
0104   int maxBin = inputHisto->FindBin( sbMax );
0105   std::cout << "minBin("<<peak.getVal() - 2.5*sigma.getVal()<<") = " << minBin << std::endl;
0106   std::cout << "maxBin("<<peak.getVal() + 2.5*sigma.getVal()<<") = " << maxBin << std::endl;
0107
0108   Double_t integral = inputHisto->Integral(minBin, maxBin);
0109
0110   x.setRange("small", sbMin, sbMax);
0111   RooAbsReal * backgroundSmallIntegral = background.createIntegral(x, "small");
0112   Double_t backgroundSmallIntegralValue = backgroundSmallIntegral->getVal();
0113   std::cout << "backgroundSmallIntegralValue = " << backgroundSmallIntegralValue << std::endl;
0114
0115   double fSigSmall = 1 - fullIntegral*(1-fSig.getVal())/integral*backgroundSmallIntegralValue/backgroundFullIntegralValue;
0116
0117
0118   std::cout << "Events in ["<< xMin << ","<<xMax<<"]:" << std::endl;
0119   std::cout << "Signal events = " << (fSig.getVal())*fullIntegral << std::endl;
0120   std::cout << "Background events = " << (1-fSig.getVal())*fullIntegral << std::endl;
0121
0122   std::cout << "Events in peak +- 2.5sigma:" << std::endl;
0123   std::cout << "Signal events = " << fSigSmall*integral << std::endl;
0124   std::cout << "Background events = " << (1-fSigSmall)*integral << std::endl;
0125
0126   frame->Draw();
0127   frame->Write();
0128   canvas.Write();
0129   outputFile.Write();
0130 };