File indexing completed on 2024-04-06 12:24:23
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
0017
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 }