File indexing completed on 2024-04-06 12:24:23
0001 #include "PhysicsTools/Utilities/interface/BreitWigner.h"
0002 #include "PhysicsTools/Utilities/interface/HistoChiSquare.h"
0003 #include "PhysicsTools/Utilities/interface/HistoPoissonLikelihoodRatio.h"
0004 #include "PhysicsTools/Utilities/interface/RootMinuitCommands.h"
0005 #include "PhysicsTools/Utilities/interface/RootMinuit.h"
0006 #include "PhysicsTools/Utilities/interface/Parameter.h"
0007 #include "PhysicsTools/Utilities/interface/rootTf1.h"
0008 #include "PhysicsTools/Utilities/interface/rootPlot.h"
0009 #include "TFile.h"
0010 #include "TH1.h"
0011 #include "TF1.h"
0012 #include "TCanvas.h"
0013 #include "TROOT.h"
0014
0015 #include <iostream>
0016 #include "PhysicsTools/Utilities/interface/Operations.h"
0017
0018
0019
0020 typedef funct::Product<funct::Parameter, funct::BreitWigner>::type FitFunction;
0021 typedef fit::HistoChiSquare<FitFunction> ChiSquared;
0022 typedef fit::HistoPoissonLikelihoodRatio<FitFunction> PoissonLR;
0023
0024 template <typename T>
0025 int main_t(const std::string tag) {
0026 try {
0027 fit::RootMinuitCommands<T> commands("PhysicsTools/Utilities/test/testZMassFit.txt");
0028
0029 const char* kYield = "Yield";
0030 const char* kMass = "Mass";
0031 const char* kGamma = "Gamma";
0032
0033 funct::Parameter yield(kYield, commands.par(kYield));
0034 funct::Parameter mass(kMass, commands.par(kMass));
0035 funct::Parameter gamma(kGamma, commands.par(kGamma));
0036 funct::BreitWigner bw(mass, gamma);
0037
0038 FitFunction f = yield * bw;
0039 TF1 startFun = root::tf1("startFun", f, 0, 200, yield, mass, gamma);
0040 TH1D histo("histo", "Z mass (GeV/c)", 200, 0, 200);
0041 histo.FillRandom("startFun", yield);
0042 TCanvas canvas;
0043 startFun.Draw();
0044 canvas.SaveAs((tag + "breitWigner.eps").c_str());
0045 histo.Draw();
0046 canvas.SaveAs((tag + "breitWignerHisto.eps").c_str());
0047 startFun.Draw("same");
0048 canvas.SaveAs((tag + "breitWignerHistoFun.eps").c_str());
0049 histo.Draw("e");
0050 startFun.Draw("same");
0051
0052 T chi2(f, &histo, 80, 120);
0053 int fullBins = chi2.numberOfBins();
0054 std::cout << "N. deg. of freedom: " << fullBins << std::endl;
0055 fit::RootMinuit<T> minuit(chi2, true);
0056 commands.add(minuit, yield);
0057 commands.add(minuit, mass);
0058 commands.add(minuit, gamma);
0059 commands.run(minuit);
0060 ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > err;
0061 minuit.getErrorMatrix(err);
0062 std::cout << "error matrix:" << std::endl;
0063 for (size_t i = 0; i < 3; ++i) {
0064 for (size_t j = 0; j < 3; ++j) {
0065 std::cout << err(i, j) << "\t";
0066 }
0067 std::cout << std::endl;
0068 }
0069 root::plot<FitFunction>((tag + "breitWignerHistoFunFit.eps").c_str(), histo, f, 80, 120, yield, mass, gamma);
0070 } catch (std::exception& err) {
0071 std::cerr << "Exception caught:\n" << err.what() << std::endl;
0072 return 1;
0073 }
0074
0075 return 0;
0076 }
0077
0078 int main() {
0079 gROOT->SetStyle("Plain");
0080 std::cout << "=== chi-2 fit ===" << std::endl;
0081 int ret1 = main_t<ChiSquared>("chi2_");
0082 if (ret1 != 0)
0083 return ret1;
0084 std::cout << "=== poisson LR fit ===" << std::endl;
0085 int ret2 = main_t<PoissonLR>("possLR_");
0086 if (ret2 != 0)
0087 return ret2;
0088 return 0;
0089 }