Logo ROOT   6.08/07
Reference Guide
ReneTest.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 // $Id$
11 #ifdef _WIN32
12  #define _USE_MATH_DEFINES
13 #endif
15 #include "Minuit2/MnMigrad.h"
16 #include "Minuit2/MnMinos.h"
18 #include "Minuit2/MnPrint.h"
19 #include "Minuit2/FCNBase.h"
20 #include "Minuit2/MnScan.h"
21 #include "Minuit2/MnPlot.h"
22 
23 using namespace ROOT::Minuit2;
24 
25 class ReneFcn : public FCNBase {
26 
27 
28 public:
29 
30  ReneFcn(const std::vector<double>& meas) : fMeasurements(meas) {}
31 
32  virtual ~ReneFcn() {}
33 
34  virtual double operator()(const std::vector<double>& par) const {
35  double a = par[2];
36  double b = par[1];
37  double c = par[0];
38  double p0 = par[3];
39  double p1 = par[4];
40  double p2 = par[5];
41  double fval = 0.;
42 // double nbin = 3./double(fMeasurements.size());
43  for(unsigned int i = 0; i < fMeasurements.size(); i++) {
44  double ni = fMeasurements[i];
45  if(ni < 1.e-10) continue;
46 // double xi = (i + 0.5)*nbin; //xi=0-3
47 // double xi = (i+0.5)/20.; //xi=0-3
48 // double xi = (i+0.5)/40.; //xi=0-3
49  double xi = (i+1.)/40. - 1./80.; //xi=0-3
50  double ei = ni;
51 // double nexp1 = a*xi*xi + b*xi + c;
52 // double nexp2 = 0.5*p0*p1/M_PI;
53 // double nexp3 = std::max(1.e-10, (xi-p2)*(xi-p2) + 0.25*p1*p1);
54 // double nexp = nexp1 + nexp2/nexp3;
55  double nexp = a*xi*xi + b*xi + c + (0.5*p0*p1/M_PI)/std::max(1.e-10, (xi-p2)*(xi-p2) + 0.25*p1*p1);
56  fval += (ni-nexp)*(ni-nexp)/ei;
57  }
58  return fval;
59  }
60 
61  virtual double Up() const {return 1.;}
62 
63 private:
64  std::vector<double> fMeasurements;
65 };
66 
67 
68 /*
69 extern "C" void fcnr_(int&, double[], double&, double[], int&);
70 extern "C" void stand_() {}
71 
72 class ReneFcn : public FCNBase {
73 
74 public:
75 
76  ReneFcn(const std::vector<double>& meas) : fMeasurements(meas) {}
77 
78  virtual ~ReneFcn() {}
79 
80  virtual double operator()(const std::vector<double>& par) const {
81  double mypar[6];
82  for(std::vector<double>::const_iterator ipar = par.begin();
83  ipar != par.end(); ipar++)
84  mypar[ipar-par.begin()] = par[ipar-par.begin()];
85  double fval = 0.;
86  int iflag = 4;
87  int npar = par.size();
88  fcnr_(npar, 0, fval, mypar, iflag);
89 
90  return fval;
91  }
92 
93  virtual double Up() const {return 1.;}
94 
95 private:
96  std::vector<double> fMeasurements;
97 };
98 
99 */
100 
101 
102 int main() {
103  /*
104  double tmp[60] = {6., 1.,10.,12., 6.,13.,23.,22.,15.,21.,
105  23.,26.,36.,25.,27.,35.,40.,44.,66.,81.,
106  75.,57.,48.,45.,46.,41.,35.,36.,53.,32.,
107  40.,37.,38.,31.,36.,44.,42.,37.,32.,32.,
108  43.,44.,35.,33.,33.,39.,29.,41.,32.,44.,
109  26.,39.,29.,35.,32.,21.,21.,15.,25.,15.};
110  std::vector<double> measurements(tmp, tmp+60);
111  */
112  /*
113  double tmp[120] = {2.,1.,1.,0.,1.,1.,0.,1.,3.,0.,0.,1.,0.,
114  1.,1.,0.,0.,1.,0.,0.,0.,0.,2.,1.,1.,2.,
115  2.,0.,2.,4.,2.,6.,2.,1.,4.,0.,3.,6.,16.,
116  30.,34.,18.,8.,2.,3.,4.,4.,5.,6.,3.,5.,
117  0.,1.,1.,7.,3.,2.,5.,1.,3.,5.,3.,2.,3.,
118  2.,2.,1.,1.,5.,2.,3.,7.,2.,7.,6.,5.,1.,
119  4.,5.,0.,6.,3.,4.,3.,3.,6.,8.,8.,3.,4.,
120  4.,8.,9.,7.,3.,4.,6.,2.,5.,10.,7.,6.,4.,
121  4.,7.,7.,5.,4.,12.,4.,6.,3.,7.,4.,3.,4.,
122  3,10.,8.,7.};
123  */
124  double tmp[120] = {38.,36.,46.,52.,54.,52.,61.,52.,64.,77.,
125  60.,56.,78.,71.,81.,83.,89.,96.,118.,96.,
126  109.,111.,107.,107.,135.,156.,196.,137.,
127  160.,153.,185.,222.,251.,270.,329.,422.,
128  543.,832.,1390.,2835.,3462.,2030.,1130.,
129  657.,469.,411.,375.,295.,281.,281.,289.,
130  273.,297.,256.,274.,287.,280.,274.,286.,
131  279.,293.,314.,285.,322.,307.,313.,324.,
132  351.,314.,314.,301.,361.,332.,342.,338.,
133  396.,356.,344.,395.,416.,406.,411.,422.,
134  393.,393.,409.,455.,427.,448.,459.,403.,
135  441.,510.,501.,502.,482.,487.,506.,506.,
136  526.,517.,534.,509.,482.,591.,569.,518.,
137  609.,569.,598.,627.,617.,610.,662.,666.,
138  652.,671.,647.,650.,701.};
139 
140  std::vector<double> measurements(tmp, tmp+120);
141 
142  ReneFcn fFCN(measurements);
143 
144  MnUserParameters upar;
145  upar.Add("p0", 100., 10.);
146  upar.Add("p1", 100., 10.);
147  upar.Add("p2", 100., 10.);
148  upar.Add("p3", 100., 10.);
149  upar.Add("p4", 1., 0.3);
150  upar.Add("p5", 1., 0.3);
151  /*
152 # ext. || Name || type || Value || Error +/-
153 
154  0 || p0 || free || 32.04 || 9.611
155  1 || p1 || free || 98.11 || 29.43
156  2 || p2 || free || 39.15 || 11.75
157  3 || p3 || free || 362.4 || 108.7
158  4 || p4 || free || 0.06707 || 0.02012
159  5 || p5 || free || 1.006 || 0.3019
160 
161  upar.Add(0, "p0", 32.04, 9.611);
162  upar.Add(1, "p1", 98.11, 29.43);
163  upar.Add(2, "p2", 39.15, 11.75);
164  upar.Add(3, "p3", 362.4, 108.7);
165  upar.Add(4, "p4", 0.06707, 0.02012);
166  upar.Add(5, "p5", 1.006, 0.3019);
167  */
168 
169  std::cout<<"initial parameters: "<<upar<<std::endl;
170 
171  std::cout<<"start migrad "<<std::endl;
172  MnMigrad migrad(fFCN, upar);
173  FunctionMinimum min = migrad();
174  if(!min.IsValid()) {
175  //try with higher strategy
176  std::cout<<"FM is invalid, try with strategy = 2."<<std::endl;
177  MnMigrad migrad2(fFCN, min.UserState(), MnStrategy(2));
178  min = migrad2();
179  }
180  std::cout<<"minimum: "<<min<<std::endl;
181  /*
182  std::cout<<"start Minos"<<std::endl;
183  MnMinos Minos(migrad, min);
184  AsymmetricError e0 = Minos(0);
185  AsymmetricError e1 = Minos(1);
186  AsymmetricError e2 = Minos(2);
187 
188  std::cout<<"par0: "<<e0.Value()<<" + "<<e0.Upper()<<e0.Lower()<<std::endl;
189  std::cout<<"par1: "<<e1.Value()<<" + "<<e1.Upper()<<e1.Lower()<<std::endl;
190  std::cout<<"par2: "<<e2.Value()<<" + "<<e2.Upper()<<e2.Lower()<<std::endl;
191  */
192 
193  {
194  std::vector<double> params(6, 1.);
195  std::vector<double> Error(6, 1.);
196  MnScan scan(fFCN, params, Error);
197  std::cout<<"scan parameters: "<<scan.Parameters()<<std::endl;
198  MnPlot plot;
199  for(unsigned int i = 0; i < upar.VariableParameters(); i++) {
200  std::vector<std::pair<double, double> > xy = scan.Scan(i);
201 // std::vector<std::pair<double, double> > xy = scan.scan(0);
202  plot(xy);
203  }
204  std::cout<<scan.Parameters()<<std::endl;
205  }
206 
207  {
208  std::vector<double> params(6, 1.);
209  std::vector<double> Error(6, 1.);
210  MnScan scan(fFCN, params, Error);
211  std::cout<<"scan parameters: "<<scan.Parameters()<<std::endl;
212  FunctionMinimum min2 = scan();
213 // std::cout<<min2<<std::endl;
214  std::cout<<scan.Parameters()<<std::endl;
215  }
216 
217  return 0;
218 }
219 
double par[1]
Definition: unuranDistr.cxx:38
API class for minimization using Variable Metric technology ("MIGRAD"); allows for user interaction: ...
Definition: MnMigrad.h:31
const MnUserParameters & Parameters() const
Definition: MnApplication.h:65
return c
MnPlot produces a text-screen graphical output of (x,y) points, e.g.
Definition: MnPlot.h:26
TArc * a
Definition: textangle.C:12
bool Add(const std::string &, double, double)
Add free Parameter Name, Value, Error.
static double p2(double t, double a, double b, double c)
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:47
int main()
Definition: ReneTest.cxx:102
#define M_PI
Definition: Rotated.cxx:105
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
static double p1(double t, double a, double b)
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
API class for minimization using a scan method to find the minimum; allows for user interaction: set/...
Definition: MnScan.h:31
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
std::vector< std::pair< double, double > > Scan(unsigned int par, unsigned int maxsteps=41, double low=0., double high=0.)
Definition: MnScan.cxx:18
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
unsigned int VariableParameters() const
const MnUserParameterState & UserState() const
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
void Error(ErrorHandler_t func, int code, const char *va_(fmt),...)
Write error message and call a handler, if required.