Logo ROOT   6.08/07
Reference Guide
PaulTest.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 #ifdef USE_SEALBASE
17 #include "SealBase/Filename.h"
18 #include "SealBase/ShellEnvironment.h"
19 #endif
20 
21 
22 #include <iostream>
23 #include <fstream>
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  int nmeas = 0;
34 
35 #ifdef USE_SEALBASE
36  seal::Filename inputFile (seal::Filename ("$SEAL/src/MathLibs/Minuit/tests/MnSim/paul.txt").substitute (seal::ShellEnvironment ()));
37  std::ifstream in(inputFile.Name() );
38 #else
39  std::ifstream in("paul.txt");
40 #endif
41  if (!in) {
42  std::cerr << "Error opening input data file" << std::endl;
43  return 1;
44  }
45 
46  // read input data
47  {
48  double x = 0., weight = 0., width = 0., err = 0.;
49  while(in>>x>>weight>>width>>err) {
50  positions.push_back(x);
51  double ni = weight*width;
52  measurements.push_back(ni);
53  var.push_back(ni);
54  nmeas += int(ni);
55  }
56  std::cout<<"size= "<<var.size()<<std::endl;
57  assert(var.size() > 0);
58  std::cout<<"nmeas: "<<nmeas<<std::endl;
59  }
60 
61  // create FCN function
62  GaussFcn fFCN(measurements, positions, var);
63 
64  std::vector<double> meas = fFCN.Measurements();
65  std::vector<double> pos = fFCN.Positions();
66 
67  // create initial starting values for parameters
68  double x = 0.;
69  double x2 = 0.;
70  double norm = 0.;
71  double dx = pos[1]-pos[0];
72  double area = 0.;
73  for(unsigned int i = 0; i < meas.size(); i++) {
74  norm += meas[i];
75  x += (meas[i]*pos[i]);
76  x2 += (meas[i]*pos[i]*pos[i]);
77  area += dx*meas[i];
78  }
79  double mean = x/norm;
80  double rms2 = x2/norm - mean*mean;
81 
82  std::cout<<"initial mean: "<<mean<<std::endl;
83  std::cout<<"initial sigma: "<<sqrt(rms2)<<std::endl;
84  std::cout<<"initial area: "<<area<<std::endl;
85 
86  MnUserParameters upar;
87  upar.Add("mean", mean, 0.1);
88  upar.Add("sigma", sqrt(rms2), 0.1);
89  upar.Add("area", area, 0.1);
90 
91  MnMigrad migrad(fFCN, upar);
92  std::cout<<"start migrad "<<std::endl;
93  FunctionMinimum min = migrad();
94  std::cout<<"minimum: "<<min<<std::endl;
95 
96  std::cout<<"start Minos"<<std::endl;
97  MnMinos Minos(fFCN, min);
98  std::pair<double,double> e0 = Minos(0);
99  std::pair<double,double> e1 = Minos(1);
100  std::pair<double,double> e2 = Minos(2);
101 
102  std::cout<<"par0: "<<min.UserState().Value("mean")<<" "<<e0.first<<" "<<e0.second<<std::endl;
103  std::cout<<"par1: "<<min.UserState().Value("sigma")<<" "<<e1.first<<" "<<e1.second<<std::endl;
104  std::cout<<"par2: "<<min.UserState().Value("area")<<" "<<e2.first<<" "<<e2.second<<std::endl;
105 
106  return 0;
107 }
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
RooCmdArg Minos(Bool_t flag=kTRUE)
std::vector< double > Measurements() const
Definition: GaussFcn.h:38
int main()
Definition: PaulTest.cxx:28
std::vector< double > Positions() const
Definition: GaussFcn.h:39
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
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
const MnUserParameterState & UserState() const
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40