Utilities/test/testExpressionFit.cpp

0001 #include "PhysicsTools/Utilities/interface/HistoChiSquare.h"
0002 #include "PhysicsTools/Utilities/interface/RootMinuitCommands.h"
0003 #include "PhysicsTools/Utilities/interface/RootMinuit.h"
0004 #include "PhysicsTools/Utilities/interface/Parameter.h"
0005 #include "PhysicsTools/Utilities/interface/rootTf1.h"
0006 #include "PhysicsTools/Utilities/interface/rootPlot.h"
0007 #include "PhysicsTools/Utilities/interface/Function.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 int main() {
0018   typedef fit::HistoChiSquare<funct::Function<funct::X> > ChiSquared;
0019   gROOT->SetStyle("Plain");
0020   try {
0021     fit::RootMinuitCommands<ChiSquared> commands("PhysicsTools/Utilities/test/testExpressionFit.txt");
0022
0023     const char* kYield = "Yield";
0024     const char* kMean = "Mean";
0025     const char* kSigma = "Sigma";
0026
0027     funct::Parameter yield(kYield, commands.par(kYield));
0028     funct::Parameter mean(kMean, commands.par(kMean));
0029     funct::Parameter sigma(kSigma, commands.par(kSigma));
0030     funct::Parameter c("C", 1. / sqrt(2 * M_PI));
0031     funct::X x;
0032     funct::Numerical<2> _2;
0033
0034     const double min = -5, max = 5;
0035
0036     funct::Function<funct::X> f = yield * c * exp(-(((x - mean) / sigma) ^ _2) / _2);
0037     TF1 startFun = root::tf1("startFun", f, min, max, yield, mean, sigma);
0038     TH1D histo("histo", "gaussian", 100, min, max);
0039     histo.FillRandom("startFun", yield);
0040     TCanvas canvas;
0041     startFun.Draw();
0042     canvas.SaveAs("gaussian.eps");
0043     histo.Draw();
0044     canvas.SaveAs("gaussianHisto.eps");
0045     startFun.Draw("same");
0046     canvas.SaveAs("gaussianHistoFun.eps");
0047     histo.Draw("e");
0048     startFun.Draw("same");
0049
0050     ChiSquared chi2(f, &histo, min, max);
0051     int fullBins = chi2.numberOfBins();
0052     std::cout << "N. deg. of freedom: " << fullBins << std::endl;
0053     fit::RootMinuit<ChiSquared> minuit(chi2, true);
0054     commands.add(minuit, yield);
0055     commands.add(minuit, mean);
0056     commands.add(minuit, sigma);
0057     commands.run(minuit);
0058     root::plot<funct::Function<funct::X> >("gaussianHistoFunFit.eps", histo, f, min, max, yield, mean, sigma);
0059   } catch (std::exception& err) {
0060     std::cerr << "Exception caught:\n" << err.what() << std::endl;
0061     return 1;
0062   }
0063
0064   return 0;
0065 }