Macros/RooFit/CrystalBallFitOnData_JPsi.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
0025 void CrystalBallFitOnData_JPsi()
0026 {
0027   TFile inputFile("0_MuScleFit.root", "READ");
0028   TH1F * inputHisto = (TH1F*)inputFile.Get("hRecBestResAllEvents_Mass");
0029   inputHisto->Rebin(4);
0030
0031   double xMin = 2.5;
0032   double xMax = 4.;
0033   RooRealVar x("x", "Mass (GeV)", xMin, xMax);
0034
0035   RooDataHist histo("dh", "dh", x, Import(*inputHisto));
0036
0037   // Build the CB
0038   RooRealVar peak("peak","peak", 3.095, 2.5, 4.);
0039   RooRealVar sigma("sigma","sigma", 0.04, 0.001, 0.1);
0040   RooRealVar alpha("alpha", "alpha", 1., 0., 10.);
0041   RooRealVar n("n", "n", 1., 0., 100.);
0042   RooCBShape signal("signal", "signal", x, peak, sigma, alpha, n);
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("crystalBallFitOnData_JPsi.root", "RECREATE");
0049
0050   // Compute the total number of events in the fit window
0051   int binXmin = inputHisto->FindBin( xMin );
0052   int binXmax = inputHisto->FindBin( xMax );
0053   Double_t fullIntegral = inputHisto->Integral(binXmin, binXmax);
0054
0055   // Build the model adding the exponential background
0056   RooRealVar nSig("nSig", "signal fraction", 30, 0., fullIntegral);
0057   RooRealVar nBkg("nBkg", "background fraction", 20, 0., fullIntegral);
0058
0059   Double_t sbSigma = 3.;
0060   // Do this twice: the first time to find the peak and sigma and the second time
0061   // to get the correct yelds in peak+-3sigma
0062   for( int i=0; i<2; ++i ) {
0063     // Compute the integral and write the signal and background events in +-sbSigma
0064     Double_t sbMin = peak.getVal() - sbSigma*sigma.getVal();
0065     Double_t sbMax = peak.getVal() + sbSigma*sigma.getVal();
0066     x.setRange("small", sbMin, sbMax);
0067     RooExtendPdf eSig("eSig", "eSig", signal, nSig, "small");
0068     RooExtendPdf eBkg("eBkg", "eBkg", background, nBkg, "small");
0069
0070     // RooAddPdf model("model", "model", RooArgList(signal, background), RooArgList(nSig, nBkg));
0071     RooAddPdf model("model", "model", RooArgList(eSig, eBkg));
0072     // model.fitTo(histo, "hs");
0073     model.fitTo(histo);
0074
0075     std::cout << "peak("<<i<<") = " << peak.getVal() << "+-" << peak.getError() << std::endl;
0076     std::cout << "sigma("<<i<<") = " << sigma.getVal() << "+-" << sigma.getError() << std::endl;
0077     if( i==1 ) {
0078       TCanvas canvas("canvas", "canvas", 1000, 800);
0079       canvas.cd();
0080       canvas.Draw();
0081       RooPlot* frame = x.frame();
0082       histo.plotOn(frame);
0083       model.plotOn(frame);
0084       model.paramOn(frame);
0085
0086       model.plotOn(frame, Components("exponential"), LineStyle(kDashed));
0087
0088       frame->Draw();
0089       frame->Write();
0090       canvas.Write();
0091       outputFile.Write();
0092     }
0093   }
0094
0095
0096   // RooAbsReal * backgroundFullIntegral = background.createIntegral(x);
0097   // Double_t backgroundFullIntegralValue = backgroundFullIntegral->getVal();
0098   // std::cout << "backgroundFullIntegralValue = " << backgroundFullIntegralValue << std::endl;
0099
0100   // int minBin = inputHisto->FindBin( sbMin );
0101   // int maxBin = inputHisto->FindBin( sbMax );
0102   // std::cout << "minBin("<<peak.getVal() - sbSigma*sigma.getVal()<<") = " << minBin << std::endl;
0103   // std::cout << "maxBin("<<peak.getVal() + sbSigma*sigma.getVal()<<") = " << maxBin << std::endl;
0104
0105   // Double_t integral = inputHisto->Integral(minBin, maxBin);
0106
0107
0108
0109   // RooAbsReal* fracSigRange = sig.createIntegral(x, x, "small");
0110   // Double_t nSigSmall = nSig.getVal()*fracSigRange->getVal();
0111
0112   // std::cout << "esig = " << esig.getVal() << "+-" << std::endl;
0113
0114
0115   // RooAbsReal * backgroundSmallIntegral = background.createIntegral(x, "small");
0116   // Double_t backgroundSmallIntegralValue = backgroundSmallIntegral->getVal();
0117   // std::cout << "backgroundSmallIntegralValue = " << backgroundSmallIntegralValue << std::endl;
0118
0119
0120   // double fSigSmall = 1 - fullIntegral*(1-fSig.getVal())/integral*backgroundSmallIntegralValue/backgroundFullIntegralValue;
0121
0122   // std::cout << "Events in ["<< xMin << ","<<xMax<<"] = " << fullIntegral << std::endl;
0123   // std::cout << "Signal events = " << (fSig.getVal())*fullIntegral << std::endl;
0124   // std::cout << "Background events = " << (1-fSig.getVal())*fullIntegral << std::endl;
0125   // std::cout << "S/B = " << fSig.getVal()/(1-fSig.getVal()) << std::endl;
0126   // std::cout << std::endl;
0127   // std::cout << "Events in peak +- "<<sbSigma<<"sigma = " << integral << std::endl;
0128   // std::cout << "Signal events = " << fSigSmall*integral << std::endl;
0129   // std::cout << "Background events = " << (1-fSigSmall)*integral << std::endl;
0130   // std::cout << "S/B = " << fSigSmall/(1-fSigSmall) << std::endl;
0131
0132 };