Utilities/test/testZMassFitLikelihood.cpp

0001 #include "PhysicsTools/Utilities/interface/Likelihood.h"
0002 #include "PhysicsTools/Utilities/interface/BreitWigner.h"
0003 #include "PhysicsTools/Utilities/interface/RootMinuitCommands.h"
0004 #include "PhysicsTools/Utilities/interface/RootMinuit.h"
0005 #include "PhysicsTools/Utilities/interface/Parameter.h"
0006 #include "PhysicsTools/Utilities/interface/rootTf1.h"
0007 #include "PhysicsTools/Utilities/interface/rootPlot.h"
0008 #include "TFile.h"
0009 #include "TH1.h"
0010 #include "TF1.h"
0011 #include "TCanvas.h"
0012 #include "TROOT.h"
0013
0014 #include <iostream>
0015 #include "PhysicsTools/Utilities/interface/Operations.h"
0016 //using namespace std;
0017 //using namespace boost;
0018
0019 int main() {
0020   gROOT->SetStyle("Plain");
0021   typedef funct::BreitWigner PDF;
0022   typedef std::vector<double> Sample;
0023   typedef fit::Likelihood<Sample, PDF> Likelihood;
0024   typedef funct::Product<funct::Parameter, PDF>::type PlotFunction;
0025   try {
0026     fit::RootMinuitCommands<Likelihood> commands("PhysicsTools/Utilities/test/testZMassFitLikelihood.txt");
0027
0028     const char* kYield = "Yield";
0029     const char* kMass = "Mass";
0030     const char* kGamma = "Gamma";
0031
0032     funct::Parameter yield(kYield, commands.par(kYield));
0033     funct::Parameter mass(kMass, commands.par(kMass));
0034     funct::Parameter gamma(kGamma, commands.par(kGamma));
0035     funct::BreitWigner bw(mass, gamma);
0036
0037     PDF pdf = bw;
0038     PlotFunction f = yield * pdf;
0039     TF1 startFun = root::tf1("startFun", f, 0, 200, yield, mass, gamma);
0040     TH1D histo("histo", "Z mass (GeV/c)", 200, 0, 200);
0041     Sample sample;
0042     sample.reserve(yield);
0043     for (unsigned int i = 0; i < yield; ++i) {
0044       double m = startFun.GetRandom();
0045       histo.Fill(m);
0046       sample.push_back(m);
0047     }
0048     TCanvas canvas;
0049     startFun.Draw();
0050     canvas.SaveAs("breitWigner.eps");
0051     histo.Draw();
0052     canvas.SaveAs("breitWignerHisto.eps");
0053     startFun.Draw("same");
0054     canvas.SaveAs("breitWignerHistoFun.eps");
0055     histo.Draw("e");
0056     startFun.Draw("same");
0057
0058     Likelihood like(sample, pdf);
0059     fit::RootMinuit<Likelihood> minuit(like, true);
0060     commands.add(minuit, mass);
0061     commands.add(minuit, gamma);
0062     commands.run(minuit);
0063     ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > err;
0064     minuit.getErrorMatrix(err);
0065     std::cout << "error matrix:" << std::endl;
0066     for (size_t i = 0; i < 2; ++i) {
0067       for (size_t j = 0; j < 2; ++j) {
0068         std::cout << err(i, j) << "\t";
0069       }
0070       std::cout << std::endl;
0071     }
0072     root::plot<PlotFunction>("breitWignerHistoFunFit.eps", histo, f, 80, 120, yield, mass, gamma);
0073   } catch (std::exception& err) {
0074     std::cerr << "Exception caught:\n" << err.what() << std::endl;
0075     return 1;
0076   }
0077
0078   return 0;
0079 }