File indexing completed on 2023-03-17 11:14:48
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
0032
0033
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
0041 RooRealVar peak("peak","peak", 3.1, 2.8, 3.8);
0042
0043 RooRealVar sigma("sigma","sigma", 0.06, 0.001, 0.1);
0044 RooGaussModel signal("signal", "signal", x, peak, sigma);
0045
0046
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
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
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
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 };