source: trunk/MagicSoft/Mars/mtools/MTMinuit.cc@ 3055

Last change on this file since 3055 was 2729, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 9.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Wolfgang Wittek, 7/2003 <mailto:wittek@mppmu.mpg.de>
19! Author(s): Thomas Bretz, 12/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MTMinuit
29//
30// Class for interfacing with Minuit
31//
32/////////////////////////////////////////////////////////////////////////////
33#include "MTMinuit.h"
34
35#include <TMinuit.h>
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40ClassImp(MTMinuit);
41
42using namespace std;
43
44// --------------------------------------------------------------------------
45//
46// Default constructor.
47//
48MTMinuit::MTMinuit(const char *name, const char *title)
49 : fFcn(NULL), fNames(NULL), fMinuit(NULL)
50{
51 fName = name ? name : "MTMinuit";
52 fTitle = title ? title : "Interface for Minuit";
53}
54
55MTMinuit::~MTMinuit()
56{
57 if (fMinuit)
58 delete fMinuit;
59}
60
61TMinuit *MTMinuit::GetMinuit() const
62{
63 return fMinuit;
64}
65
66void MTMinuit::SetFcn(void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t))
67{
68 fFcn = fcn;
69}
70
71void MTMinuit::InitParameters(const TArrayD &vinit, const TArrayD *step, const TString *name)
72{
73 const int n = fVinit.GetSize();
74
75 fVinit = vinit;
76
77 fStep.Set(n);
78 memset(fStep.GetArray(), 0, n*sizeof(Double_t));
79
80 if (step && step->GetSize()!=n)
81 {
82 *fLog << warn << "MTMinuit::SetParameters: number of steps doesn't match number of parameters... ignored." << endl;
83 step = NULL;
84 }
85
86 for (int i=0; i<n; i++)
87 fStep[i] = step ? (*step)[i] : TMath::Abs(fVinit[i]/10.0);
88
89 if (fNames)
90 delete fNames;
91
92 fNames = new TString[n];
93
94 // FIXME: Sanity check for TString array size!
95 for (int i=0; i<n ; i++)
96 fNames[i] = name ? (const char*)name[i] : Form("par%02d", i);
97
98 // ---------------------------------------------------------------
99
100 fLimLo.Set(n);
101 fLimUp.Set(n);
102 memset(fLimLo.GetArray(), 0, n*sizeof(Double_t));
103 memset(fLimUp.GetArray(), 0, n*sizeof(Double_t));
104
105 fFix.Set(n);
106 memset(fFix.GetArray(), 0, n*sizeof(Char_t));
107}
108
109void MTMinuit::SetLimits(const TArrayD &limlo, const TArrayD &limup)
110{
111 if (fVinit.GetSize()==limlo.GetSize())
112 fLimLo = limlo;
113 else
114 *fLog << warn << "MTMinuit::SetLimits: size of limlo doesn't match number of parameters... ignored." << endl;
115
116 if (fVinit.GetSize()==limup.GetSize())
117 fLimUp = limup;
118 else
119 *fLog << warn << "MTMinuit::SetLimits: size of limup doesn't match number of parameters... ignored." << endl;
120}
121
122void MTMinuit::SetFixedParameters(const TArrayC &fix)
123{
124 if (fVinit.GetSize()==fix.GetSize())
125 fFix = fix;
126 else
127 *fLog << warn << "MTMinuit::SetFixedParameters: size of fix doesn't match number of parameters... ignored." << endl;
128}
129
130// -----------------------------------------------------------------------
131//
132// Interface to MINUIT
133//
134Bool_t MTMinuit::CallMinuit(TObject *objectfit , const TString &method, Bool_t nulloutput)
135{
136 if (!fFcn(NULL))
137 {
138 *fLog << err << "CallMinuit: ERROR - fFcn not set... abort." << endl;
139 return kFALSE;
140 }
141 if (!fNames)
142 {
143 *fLog << err << "CallMinuit: ERROR - Parameters not set... abort." << endl;
144 return kFALSE;
145 }
146
147 // Make a copy of fStep
148 TArrayD step(fStep);
149
150 // Save gMinuit
151 TMinuit *const save = gMinuit;
152
153 // Set the maximum number of parameters
154 if (fMinuit)
155 delete fMinuit;
156 fMinuit = new TMinuit(fVinit.GetSize());
157
158 //
159 // Set the print level
160 // -1 no output except SHOW comands
161 // 0 minimum output
162 // 1 normal output (default)
163 // 2 additional ouput giving intermediate results
164 // 3 maximum output, showing progress of minimizations
165 //
166 fMinuit->SetPrintLevel(-1);
167
168 //..............................................
169 // Printout for warnings
170 // SET WAR print warnings
171 // SET NOW suppress warnings
172 Int_t errWarn;
173 Double_t tmpwar = 0;
174 fMinuit->mnexcm("SET NOW", &tmpwar, 0, errWarn);
175 //fMinuit->mnexcm("SET WAR", &tmpwar, 0, errWarn);
176
177 //..............................................
178 // Set the address of the minimization function
179 fMinuit->SetFCN(fFcn);
180
181 //..............................................
182 // Store address of object to be used in fcn
183 fMinuit->SetObjectFit(objectfit);
184
185 //..............................................
186 // Set starting values and step sizes for parameters
187 for (Int_t i=0; i<fVinit.GetSize(); i++)
188 {
189 if (!fMinuit->DefineParameter(i, fNames[i], fVinit[i], fStep[i], fLimLo[i], fLimUp[i]))
190 continue;
191
192 *fLog << err << "CallMinuit: Error in defining parameter " << fNames[i] << endl;
193 return kFALSE;
194 }
195
196 //
197 // Error definition :
198 //
199 // for chisquare function :
200 // up = 1.0 means calculate 1-standard deviation error
201 // = 4.0 means calculate 2-standard deviation error
202 //
203 // for log(likelihood) function :
204 // up = 0.5 means calculate 1-standard deviation error
205 // = 2.0 means calculate 2-standard deviation error
206 //
207 fMinuit->SetErrorDef(1.0);
208
209 // Int_t errMigrad;
210 // Double_t tmp = 0;
211 // fMinuit->mnexcm("MIGRAD", &tmp, 0, errMigrad);
212
213 // fix a parameter
214 for (Int_t i=0; i<fVinit.GetSize(); i++)
215 {
216 if (!fFix[i])
217 continue;
218
219 fMinuit->FixParameter(i);
220 }
221
222 //..............................................
223 // This doesn't seem to have any effect
224 // Set maximum number of iterations (default = 500)
225 //Int_t maxiter = 100000;
226 //fMinuit->SetMaxIterations(maxiter);
227
228 // minimization by the method of Migrad
229 Int_t fErrMinimize=0;
230 if (method.Contains("Migrad", TString::kIgnoreCase))
231 {
232 if (nulloutput)
233 fLog->SetNullOutput(kTRUE);
234 Double_t tmp = 0;
235 fMinuit->mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
236 if (nulloutput)
237 fLog->SetNullOutput(kFALSE);
238 }
239
240 //..............................................
241 // same minimization as by Migrad
242 // but switches to the SIMPLEX method if MIGRAD fails to converge
243 if (method.Contains("Minimize", TString::kIgnoreCase))
244 {
245 Double_t tmp = 0;
246 fMinuit->mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
247 }
248
249 //..............................................
250 // minimization by the SIMPLEX method
251 if (method.Contains("Simplex", TString::kIgnoreCase))
252 {
253 Double_t tmp[2];
254 tmp[0] = 3000; // maxcalls;
255 tmp[1] = 0.1; // tolerance;
256
257 if (nulloutput)
258 fLog->SetNullOutput(kTRUE);
259 fMinuit->mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize);
260 if (nulloutput)
261 fLog->SetNullOutput(kFALSE);
262 }
263
264 //..............................................
265 // check quality of minimization
266 // istat = 0 covariance matrix not calculated
267 // 1 diagonal approximation only (not accurate)
268 // 2 full matrix, but forced positive-definite
269 // 3 full accurate covariance matrix
270 // (indication of normal convergence)
271 Double_t fMin, fEdm, fErrdef;
272 Int_t fNpari, fNparx, fIstat;
273 fMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
274
275 if (fErrMinimize != 0)
276 {
277 *fLog << err << "CallMinuit: Minimization failed with:" << endl;
278 *fLog << " fMin = " << fMin << endl;
279 *fLog << " fEdm = " << fEdm << endl;
280 *fLog << " fErrdef = " << fErrdef << endl;
281 *fLog << " fIstat = " << fIstat << endl;
282 *fLog << " fErrMinimize = " << fErrMinimize << endl;
283 return kFALSE;
284 }
285
286 //..............................................
287 // minimization by the method of Migrad
288 if (method.Contains("Hesse", TString::kIgnoreCase))
289 {
290 Double_t tmp = 0;
291 fMinuit->mnexcm("HESSE", &tmp, 0, fErrMinimize);
292 }
293
294 //..............................................
295 // Minos error analysis
296 if (method.Contains("Minos", TString::kIgnoreCase))
297 {
298 Double_t tmp = 0;
299 fMinuit->mnexcm("MINOS", &tmp, 0, fErrMinimize);
300 }
301
302 //..............................................
303 // Print current status of minimization
304 // if nkode = 0 only function value
305 // 1 parameter values, errors, limits
306 // 2 values, errors, step sizes, internal values
307 // 3 values, errors, step sizes, 1st derivatives
308 // 4 values, parabolic errors, MINOS errors
309
310 //Int_t nkode = 4;
311 //fMinuit->mnprin(nkode, fmin);
312
313 //..............................................
314 // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
315 // iflag = 1 initial calculations only
316 // 2 calculate 1st derivatives and function
317 // 3 calculate function only
318 // 4 calculate function + final calculations
319 Double_t iflag = 3;
320 Int_t errfcn3;
321 fMinuit->mnexcm("CALL", &iflag, 1, errfcn3);
322
323 // WW : the following statements were commented out because the
324 // Minuit object will still be used;
325 // this may be changed in the future
326 gMinuit = save;
327
328 return kTRUE;
329}
Note: See TracBrowser for help on using the repository browser.