Macros/RooFit/GaussianFitOnData_Psi2S.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 #include "RooExtendPdf.h"
0024 using namespace RooFit;
0025
0026 void GaussianFitOnData_Psi2S()
0027 {
0028   TFile inputFile("0_MuScleFit.root", "READ");
0029   TH1F * inputHisto = (TH1F*)inputFile.Get("hRecBestResAllEvents_Mass");
0030   inputHisto->Rebin(10);
0031
0032   double xMin = 3.3;
0033   double xMax = 4.4;
0034   RooRealVar x("x", "Mass (GeV)", xMin, xMax);
0035
0036   RooDataHist histo("dh", "dh", x, Import(*inputHisto));
0037
0038   // Build the Gaussian
0039   RooRealVar peak("peak","peak", 3.695, 3.6, 3.8);
0040   RooRealVar sigma("sigma","sigma", 0.06, 0.001, 0.1);
0041   // RooRealVar sigma("sigma","sigma", 0.044);
0042   RooGaussian signal("signal", "signal", x, peak, sigma);
0043
0044   // Build the exponential background pdf
0045   RooRealVar expCoeff("expCoeff", "exponential coefficient", -1., -10., 10.);
0046   RooExponential background("exponential", "exponential", x, expCoeff);
0047
0048   TFile outputFile("gaussianFitOnData_Psi2S.root", "RECREATE");
0049
0050   // // Fit the background in the sidebands
0051   // x.setRange("sb_lo", 2., 2.6);
0052   // x.setRange("sb_hi", 3.8, 4.4);
0053   // background.fitTo(histo, Range("sb_lo,sb_hi"));
0054
0055   // TCanvas canvasBackground("canvasBackgroundFit", "canvasBackgroundFit", 1000, 800);
0056   // canvasBackground.cd();
0057   // canvasBackground.Draw();
0058   // RooPlot* frameBackground = x.frame();
0059   // histo.plotOn(frameBackground);
0060   // background.plotOn(frameBackground);
0061   // background.paramOn(frameBackground);
0062
0063   // frameBackground->Draw();
0064   // frameBackground->Write();
0065   // canvasBackground.Write();
0066
0067   // canvasBackground.Write();
0068
0069   // // Fix the background parameter
0070   // expCoeff.setConstant(kTRUE);
0071
0072
0073   // // Prefit to determine peak and sigma
0074   // RooRealVar fSig("fSig", "signal fraction", 0.4, 0., 1.);
0075   // RooAddPdf prefitModel("prefitModel", "prefitModel", RooArgList(signal, background), fSig);
0076   // prefitModel.fitTo(histo, "hs");
0077   // prefitModel.fitTo(histo, "hs");
0078   // std::cout << "prefit peak = " << peak.getVal() << "+-" << peak.getError() << std::endl;
0079   // std::cout << "prefit sigma = " << sigma.getVal() << "+-" << sigma.getError() << std::endl;
0080
0081   // Now we have peak and sigma to the actual values, build an extended model to compute
0082   // also the yelds of signal and background in the peak+-3sigma window
0083
0084   int binXmin = inputHisto->FindBin( xMin );
0085   int binXmax = inputHisto->FindBin( xMax );
0086   Double_t fullIntegral = inputHisto->Integral(binXmin, binXmax);
0087   std::cout << "Events in fit interval = " << fullIntegral << std::endl;
0088
0089   // Build the model adding the exponential background
0090   RooRealVar nSig("nSig", "signal fraction", 30, 0., fullIntegral);
0091   RooRealVar nBkg("nBkg", "background fraction", 20, 0., fullIntegral);
0092
0093   Double_t sbSigma = 3.;
0094   // Do this twice: the first time to find the peak and sigma and the second time
0095   // to get the correct yelds in peak+-3sigma
0096   for( int i=0; i<2; ++i ) {
0097     // Compute the integral and write the signal and background events in +-sbSigma
0098     Double_t sbMin = peak.getVal() - sbSigma*sigma.getVal();
0099     Double_t sbMax = peak.getVal() + sbSigma*sigma.getVal();
0100     x.setRange("small", sbMin, sbMax);
0101     RooExtendPdf eSig("eSig", "eSig", signal, nSig, "small");
0102     RooExtendPdf eBkg("eBkg", "eBkg", background, nBkg, "small");
0103
0104     // RooAddPdf model("model", "model", RooArgList(signal, background), RooArgList(nSig, nBkg));
0105     RooAddPdf model("model", "model", RooArgList(eSig, eBkg));
0106     // model.fitTo(histo, "hs");
0107     model.fitTo(histo);
0108
0109     std::cout << "peak("<<i<<") = " << peak.getVal() << "+-" << peak.getError() << std::endl;
0110     std::cout << "sigma("<<i<<") = " << sigma.getVal() << "+-" << sigma.getError() << std::endl;
0111     if( i==1 ) {
0112       TCanvas canvas("canvas", "canvas", 1000, 800);
0113       canvas.cd();
0114       canvas.Draw();
0115       RooPlot* frame = x.frame();
0116       histo.plotOn(frame);
0117       model.plotOn(frame);
0118       model.paramOn(frame);
0119
0120       model.plotOn(frame, Components("exponential"), LineStyle(kDashed));
0121
0122       frame->Draw();
0123       frame->Write();
0124       canvas.Write();
0125       outputFile.Write();
0126     }
0127   }
0128
0129
0130   // RooAbsReal * backgroundFullIntegral = background.createIntegral(x);
0131   // Double_t backgroundFullIntegralValue = backgroundFullIntegral->getVal();
0132   // std::cout << "backgroundFullIntegralValue = " << backgroundFullIntegralValue << std::endl;
0133
0134   // int minBin = inputHisto->FindBin( sbMin );
0135   // int maxBin = inputHisto->FindBin( sbMax );
0136   // std::cout << "minBin("<<peak.getVal() - sbSigma*sigma.getVal()<<") = " << minBin << std::endl;
0137   // std::cout << "maxBin("<<peak.getVal() + sbSigma*sigma.getVal()<<") = " << maxBin << std::endl;
0138
0139   // Double_t integral = inputHisto->Integral(minBin, maxBin);
0140
0141
0142
0143   // RooAbsReal* fracSigRange = sig.createIntegral(x, x, "small");
0144   // Double_t nSigSmall = nSig.getVal()*fracSigRange->getVal();
0145
0146   // std::cout << "esig = " << esig.getVal() << "+-" << std::endl;
0147
0148
0149   // RooAbsReal * backgroundSmallIntegral = background.createIntegral(x, "small");
0150   // Double_t backgroundSmallIntegralValue = backgroundSmallIntegral->getVal();
0151   // std::cout << "backgroundSmallIntegralValue = " << backgroundSmallIntegralValue << std::endl;
0152
0153
0154   // double fSigSmall = 1 - fullIntegral*(1-fSig.getVal())/integral*backgroundSmallIntegralValue/backgroundFullIntegralValue;
0155
0156   // std::cout << "Events in ["<< xMin << ","<<xMax<<"] = " << fullIntegral << std::endl;
0157   // std::cout << "Signal events = " << (fSig.getVal())*fullIntegral << std::endl;
0158   // std::cout << "Background events = " << (1-fSig.getVal())*fullIntegral << std::endl;
0159   // std::cout << "S/B = " << fSig.getVal()/(1-fSig.getVal()) << std::endl;
0160   // std::cout << std::endl;
0161   // std::cout << "Events in peak +- "<<sbSigma<<"sigma = " << integral << std::endl;
0162   // std::cout << "Signal events = " << fSigSmall*integral << std::endl;
0163   // std::cout << "Background events = " << (1-fSigSmall)*integral << std::endl;
0164   // std::cout << "S/B = " << fSigSmall/(1-fSigSmall) << std::endl;
0165
0166 };