ROOT  6.06/09
Reference Guide
PaulTest2.cxx
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
10 #include "GaussFcn.h"
12 #include "Minuit2/MnMigrad.h"
13 #include "Minuit2/MnMinos.h"
14 #include "Minuit2/MnPrint.h"
15 
16 #include <iostream>
17 #include <fstream>
18 
19 
20 #ifdef USE_SEALBASE
21 #include "SealBase/Filename.h"
22 #include "SealBase/ShellEnvironment.h"
23 #endif
24 
25 using namespace ROOT::Minuit2;
26 
27 
28 int main() {
29 
30  std::vector<double> positions;
31  std::vector<double> measurements;
32  std::vector<double> var;
33  double nmeas = 0;
34 
35 #ifdef USE_SEALBASE
36  seal::Filename inputFile (seal::Filename ("$SEAL/src/MathLibs/Minuit/tests/MnSim/paul2.txt").substitute (seal::ShellEnvironment ()));
37  std::ifstream in(inputFile.Name() );
38 #else
39  std::ifstream in("paul2.txt");
40 #endif
41  if (!in) {
42  std::cerr << "Error opening input data file" << std::endl;
43  return 1;
44  }
45 
46 
47  // read input data
48  {
49  double x = 0., y = 0., width = 0., err = 0., un1 = 0., un2 = 0.;
50  while(in>>x>>y>>width>>err>>un1>>un2) {
51  if(err < 1.e-8) continue;
52  positions.push_back(x);
53  measurements.push_back(y);
54  var.push_back(err*err);
55  nmeas += y;
56  }
57  std::cout<<"size= "<<var.size()<<std::endl;
58  assert(var.size() > 0);
59  std::cout<<"nmeas: "<<nmeas<<std::endl;
60  }
61 
62  // create FCN function
63  GaussFcn fFCN(measurements, positions, var);
64 
65  std::vector<double> meas = fFCN.Measurements();
66  std::vector<double> pos = fFCN.Positions();
67 
68  // create initial starting values for parameters
69  double x = 0.;
70  double x2 = 0.;
71  double norm = 0.;
72  double area = 0.;
73  double dx = pos[1]-pos[0];
74  for(unsigned int i = 0; i < meas.size(); i++) {
75  norm += meas[i];
76  x += (meas[i]*pos[i]);
77  x2 += (meas[i]*pos[i]*pos[i]);
78  area += dx*meas[i];
79  }
80  double mean = x/norm;
81  double rms2 = x2/norm - mean*mean;
82 
83  std::cout<<"initial mean: "<<mean<<std::endl;
84  std::cout<<"initial sigma: "<<sqrt(rms2)<<std::endl;
85  std::cout<<"initial area: "<<area<<std::endl;
86  std::vector<double> init_val(3);
87  init_val[0] = mean;
88  init_val[1] = sqrt(rms2);
89  init_val[2] = area;
90  std::cout<<"initial fval: "<<fFCN(init_val)<<std::endl;
91 
92  MnUserParameters upar;
93  upar.Add("mean", mean, 1.);
94  upar.Add("sigma", sqrt(rms2), 1.);
95  upar.Add("area", area, 10.);
96 
97  MnMigrad migrad(fFCN, upar);
98  std::cout<<"start migrad "<<std::endl;
99  FunctionMinimum min = migrad();
100  std::cout<<"minimum: "<<min<<std::endl;
101 
102  std::cout<<"start Minos"<<std::endl;
103  MnMinos Minos(fFCN, min);
104  std::pair<double,double> e0 = Minos(0);
105  std::pair<double,double> e1 = Minos(1);
106  std::pair<double,double> e2 = Minos(2);
107 
108  std::cout<<"par0: "<<min.UserState().Value("mean")<<" "<<e0.first<<" "<<e0.second<<std::endl;
109  std::cout<<"par1: "<<min.UserState().Value(1)<<" "<<e1.first<<" "<<e1.second<<std::endl;
110  std::cout<<"par2: "<<min.UserState().Value(2)<<" "<<e2.first<<" "<<e2.second<<std::endl;
111 
112  return 0;
113 }
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
RooCmdArg Minos(Bool_t flag=kTRUE)
#define assert(cond)
Definition: unittest.h:542
const MnUserParameterState & UserState() const
bool Add(const std::string &, double, double)
Add free Parameter Name, Value, Error.
double sqrt(double)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
API class for Minos Error analysis (asymmetric errors); minimization has to be done before and Minimu...
Definition: MnMinos.h:34
std::vector< double > Positions() const
Definition: GaussFcn.h:39
int main()
Definition: PaulTest2.cxx:28
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
Double_t y[n]
Definition: legend1.C:17
std::vector< double > Measurements() const
Definition: GaussFcn.h:38
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40