Utilities/test/testZMassFitExtLikelihood.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 funct::Product<funct::Parameter, PDF>::type FitFunction;
0024   typedef fit::Likelihood<Sample, FitFunction, funct::Parameter> Likelihood;
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     FitFunction 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, f, yield);
0059     fit::RootMinuit<Likelihood> minuit(like, true);
0060     commands.add(minuit, yield);
0061     commands.add(minuit, mass);
0062     commands.add(minuit, gamma);
0063     commands.run(minuit);
0064     ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > err;
0065     minuit.getErrorMatrix(err);
0066     std::cout << "error matrix:" << std::endl;
0067     for (size_t i = 0; i < 3; ++i) {
0068       for (size_t j = 0; j < 3; ++j) {
0069         std::cout << err(i, j) << "\t";
0070       }
0071       std::cout << std::endl;
0072     }
0073     root::plot<FitFunction>("breitWignerHistoFunFit.eps", histo, f, 80, 120, yield, mass, gamma);
0074   } catch (std::exception& err) {
0075     std::cerr << "Exception caught:\n" << err.what() << std::endl;
0076     return 1;
0077   }
0078
0079   return 0;
0080 }