Skip to content
Snippets Groups Projects
Commit 75db82aa authored by Jan Schulze's avatar Jan Schulze
Browse files

initial import of MinuitFit

parent e1aff891
No related branches found
No related tags found
No related merge requests found
build-project DfuncClebschG ;
build-project MinuitFit ;
project :
;
lib MinuitFit :
[ glob *.cc : *App.cc ]
$(TOP)/qft++//qft++
$(TOP)/ErrLogger//ErrLogger
: <use>$(TOP)//Minuit2 <use>$(TOP)//Geneva
:
: <library>$(TOP)//Minuit2 <library>$(TOP)//Geneva ;
exe MinuitFitApp : MinuitFitApp.cc MinuitFit : ;
#include <getopt.h>
#include <fstream>
#include <sstream>
#include <string>
#include "Examples/Tutorial/MinuitFit/MinuitFit.hh"
#include "Examples/Tutorial/MinuitFit/MinuitFitFcn.hh"
#include "TFile.h"
#include "TMath.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TGraphErrors.h"
#include "ErrLogger/ErrLogger.hh"
MinuitFit::MinuitFit(double p0, double p1, double p2, double p3, double sigma) :
_theTFile(0),
_sigma(sigma)
{
// Display parameters for test distribution
cout << endl;
Info <<"Set p0 as "<< p0 << endmsg;
Info <<"Set p1 as "<< p1 << endmsg;
Info <<"Set p2 as "<< p2 << endmsg;
Info <<"Set p3 as "<< p3 << endmsg;
Info <<"Set sigma as "<< sigma << endmsg;
// Generate test distribution and smear them with a gaussian
TRandom randomNumber;
for(int i=0; i<1000; i++) {
double tmpXvalue=static_cast<double>((rand() % 10000 + 1)) / 100;
double tmpYvalue=randomNumber.Gaus(p0 + p1 * tmpXvalue + p2 * tmpXvalue * tmpXvalue + p3 * tmpXvalue * tmpXvalue * tmpXvalue, _sigma);
_xValue.push_back(tmpXvalue);
_yValue.push_back(tmpYvalue);
}
}
double MinuitFit::calcChiSqr(const std::vector<double>& minPar){
// Calculate chi^2 for current set of fit parameters
double result=0.;
for (unsigned int i=0; i<_xValue.size(); i++){
double yValFit=minPar[0]+minPar[1]*_xValue[i]+minPar[2]*_xValue[i]*_xValue[i]+minPar[3]*_xValue[i]*_xValue[i]*_xValue[i];
double yValExp=_yValue[i];
double tmpChi=((yValExp-yValFit)*(yValExp-yValFit))/(_sigma*_sigma);
result+=tmpChi;
}
return result;
}
void MinuitFit::drawGraph(double a, double b, double c, double d){
// Create arrays (which will be needed for feeding the TGraph)
const unsigned int dataEntries=_xValue.size();
double x[dataEntries];
double y[dataEntries];
double xerr[dataEntries];
double yerr[dataEntries];
// Convert vectors to arrays
for(unsigned int i=0; i<dataEntries; i++) {
x[i]=_xValue[i];
y[i]=_yValue[i];
xerr[i] = 0;
yerr[i] = _sigma;
}
// Create ROOT file
_theTFile = new TFile("myFit.root","RECREATE");
// Create TGraph and write to ROOT file
TCanvas* c1=new TCanvas("c1","c1",1200,800);
c1->cd();
TGraphErrors* linDist = new TGraphErrors(100, x, y, xerr, yerr);
linDist->SetMarkerStyle(6);
linDist->Draw("AP");
linDist->Write("myGraph");
// Draw fit with calculated parameters
TF1* fit1 = new TF1("fit1","pol3",0.,100.);
fit1->SetParameters(a, b, c, d);
fit1->SetLineColor(46);
fit1->SetLineWidth(3);
fit1->Draw();
fit1->Write();
}
MinuitFit::~MinuitFit()
{
_theTFile->Write();
_theTFile->Close();
}
#ifndef _MinuitFit_H
#define _MinuitFit_H
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <map>
#include <cassert>
#include <boost/shared_ptr.hpp>
#include "TROOT.h"
#include "qft++/topincludes/relativistic-quantum-mechanics.hh"
class TFile;
class TGraph;
class TCanvas;
class TRandom;
class MinuitFit {
public:
// create/copy/destroy:
///Constructor
MinuitFit(double p0, double p1, double p2, double p3, double sigma);
/** Destructor */
virtual ~MinuitFit();
double calcChiSqr(const std::vector<double>& minPar);
void drawGraph(double a, double b, double c, double d);
// Getters:
protected:
private:
TFile* _theTFile;
std::map <unsigned int, TGraph* > _myGraph;
std::vector< double > _xValue;
std::vector< double > _yValue;
double _sigma;
};
#endif
#include <iostream>
#include <cstring>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <boost/shared_ptr.hpp>
#include "Examples/Tutorial/MinuitFit/MinuitFit.hh"
#include "Examples/Tutorial/MinuitFit/MinuitFitFcn.hh"
#include "ErrLogger/ErrLogger.hh"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnStrategy.h"
using namespace ROOT::Minuit2;
int main(int __argc,char *__argv[]){
ErrLogger::instance()->setLevel(log4cpp::Priority::INFO);
if( __argc>1 && ( strcmp( __argv[1], "-help" ) == 0
|| strcmp( __argv[1], "--help" ) == 0 ) ){
Info << "\nThis is a test application for histogramming the Argand plot, the phase shift and the mass shape of the relativistic Breit-Wigner function with Blatt-Weisskopf barrier factors\n"
<< "The switches are:\n\n"
<< "-p0 (parameter p0 of test distribution; default -10.)\n\n"
<< "-p1 (parameter p1 of test distribution; default 10.)\n\n"
<< "-p2 (parameter p2 of test distribution; default 1.)\n\n"
<< "-p3 (parameter p3 of test distribution; default -0.01)\n\n"
<< "-sigma (sigma of gaussian smearing; default 3)\n\n"
<< endmsg;
return 0;
}
cout << endl;
cout << "___Starting..." << endl;
// Set default values
std::string p0Str="-10.";
std::string p1Str="10.";
std::string p2Str="1.";
std::string p3Str="-0.01";
std::string sigmaStr="3";
// Read arguments and replace default values
while ((optind < (__argc-1) ) && (__argv[optind][0]=='-')) {
bool found=false;
std::string sw = __argv[optind];
if (sw=="-p0"){
optind++;
p0Str = __argv[optind];
found=true;
}
if (sw=="-p1"){
optind++;
p1Str = __argv[optind];
found=true;
}
if (sw=="-p2"){
optind++;
p2Str = __argv[optind];
found=true;
}
if (sw=="-p3"){
optind++;
p3Str = __argv[optind];
found=true;
}
if (sw=="-sigma"){
optind++;
sigmaStr = __argv[optind];
found=true;
}
if (!found){
Warning << "Unknown switch: "
<< __argv[optind] << endmsg;
optind++;
}
while ( (optind < __argc ) && __argv[optind][0]!='-' ) optind++;
}
// Convert argument strings to doubles
std::stringstream p0StrStr(p0Str);
double p0=0.;
p0StrStr >> p0;
std::stringstream p1StrStr(p1Str);
double p1=0.;
p1StrStr >> p1;
std::stringstream p2StrStr(p2Str);
double p2=0.;
p2StrStr >> p2;
std::stringstream p3StrStr(p3Str);
double p3=0.;
p3StrStr >> p3;
std::stringstream sigmaStrStr(sigmaStr);
double sigma=0.;
sigmaStrStr >> sigma;
// Generate data distribution
boost::shared_ptr<MinuitFit> minuitFit(new MinuitFit(p0, p1, p2, p3, sigma));
MinuitFitFcn minuitFitFcn(minuitFit);
// Set user parameters for MinuitFitFcn
MnUserParameters upar;
upar.Add("a", -11, 3, 0., -20.);
upar.Add("b", 9.8, 2., 15., 5.);
upar.Add("c", 1.1, 0.3, 1.5, 0.5);
upar.Add("d", -0.008, 0.005, 0., -0.02);
MnMigrad migrad(minuitFitFcn, upar);
cout << endl;
Info <<"Start Migrad "<< endmsg;
FunctionMinimum min = migrad();
if(!min.IsValid()) {
// Try with higher strategy
Info <<"FM is invalid, try with strategy = 2."<< endmsg;
MnMigrad migrad2(minuitFitFcn, min.UserState(), MnStrategy(2));
min = migrad2();
}
Info << "Function value of minimization: " << min.Fval() << endmsg;
// Save final fit parameters and their errors in variables
double final_a = min.UserState().Value("a");
double final_b = min.UserState().Value("b");
double final_c = min.UserState().Value("c");
double final_d = min.UserState().Value("d");
double ea = min.UserState().Error("a");
double eb = min.UserState().Error("b");
double ec = min.UserState().Error("c");
double ed = min.UserState().Error("d");
cout << endl;
Info << "final a:\t" << final_a << " +- " << ea << endmsg;
Info << "final b:\t" << final_b << " +- " << eb << endmsg;
Info << "final c:\t" << final_c << " +- " << ec << endmsg;
Info << "final d:\t" << final_d << " +- " << ed << endmsg;
// Draw the graph and the fit in one canvas
minuitFit->drawGraph(final_a, final_b, final_c, final_d);
cout << endl;
cout << "___Finished" << endl;
cout << endl;
return 0;
}
#include "Examples/Tutorial/MinuitFit/MinuitFitFcn.hh"
#include "Examples/Tutorial/MinuitFit/MinuitFit.hh"
#include "ErrLogger/ErrLogger.hh"
#include <cassert>
using namespace ROOT::Minuit2;
MinuitFitFcn::MinuitFitFcn(boost::shared_ptr<MinuitFit> minuitFit) :
_minFitPtr(minuitFit)
{
if (0==_minFitPtr) {
Alert << "MinuitFit pointer is 0 !!!!" << endmsg;
exit(1);
}
}
MinuitFitFcn::~MinuitFitFcn()
{
}
double MinuitFitFcn::operator()(const std::vector<double>& par) const
{
double result=_minFitPtr->calcChiSqr(par);
DebugMsg << "current chi^2:\t"<< result << endmsg;
return result;
}
double MinuitFitFcn::Up() const
{
return 1.;
}
#ifndef _MinuitFitFcn_H
#define _MinuitFitFcn_H
#include <iostream>
#include <fstream>
//#include <string>
#include <vector>
#include <boost/shared_ptr.hpp>
//#include <cassert>
#include "Minuit2/FCNBase.h"
class MinuitFit;
namespace ROOT {
namespace Minuit2 {
class MinuitFitFcn : public FCNBase {
public:
MinuitFitFcn(boost::shared_ptr<MinuitFit>);
virtual ~MinuitFitFcn();
double operator()(const std::vector<double>& par) const;
double Up() const;
private:
boost::shared_ptr<MinuitFit> _minFitPtr;
};
} // namespace Minuit2
} // namespace ROOT
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment