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