Utilities/test/testZMassFit.cpp

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 //using namespace std;
0018 //using namespace boost;
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 }