| 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 | !
|
|---|
| 20 | ! Copyright: MAGIC Software Development, 2000-2003
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 |
|
|---|
| 25 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 26 | // //
|
|---|
| 27 | // MMinuitInterface //
|
|---|
| 28 | // //
|
|---|
| 29 | // Class for interfacing with Minuit //
|
|---|
| 30 | // //
|
|---|
| 31 | // //
|
|---|
| 32 | // //
|
|---|
| 33 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 34 | #include "MMinuitInterface.h"
|
|---|
| 35 |
|
|---|
| 36 | #include <math.h> // fabs
|
|---|
| 37 |
|
|---|
| 38 | #include <TArrayD.h>
|
|---|
| 39 | #include <TArrayI.h>
|
|---|
| 40 | #include <TMinuit.h>
|
|---|
| 41 | #include <TStopwatch.h>
|
|---|
| 42 |
|
|---|
| 43 | #include "MLog.h"
|
|---|
| 44 | #include "MLogManip.h"
|
|---|
| 45 | #include "MParContainer.h"
|
|---|
| 46 |
|
|---|
| 47 | ClassImp(MMinuitInterface);
|
|---|
| 48 |
|
|---|
| 49 | using namespace std;
|
|---|
| 50 |
|
|---|
| 51 |
|
|---|
| 52 | // --------------------------------------------------------------------------
|
|---|
| 53 | //
|
|---|
| 54 | // Default constructor.
|
|---|
| 55 | //
|
|---|
| 56 | MMinuitInterface::MMinuitInterface(const char *name, const char *title)
|
|---|
| 57 | {
|
|---|
| 58 | fName = name ? name : "MMinuitInterface";
|
|---|
| 59 | fTitle = title ? title : "Interface for Minuit";
|
|---|
| 60 | }
|
|---|
| 61 |
|
|---|
| 62 | // -----------------------------------------------------------------------
|
|---|
| 63 | //
|
|---|
| 64 | // Interface to MINUIT
|
|---|
| 65 | //
|
|---|
| 66 | //
|
|---|
| 67 | Bool_t MMinuitInterface::CallMinuit(
|
|---|
| 68 | void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
|
|---|
| 69 | const TString* name, const TArrayD &vinit, const TArrayD &step2,
|
|---|
| 70 | const TArrayD &limlo, const TArrayD &limup, const TArrayI &fix,
|
|---|
| 71 | TObject *objectfit , const TString &method, Bool_t nulloutput)
|
|---|
| 72 | {
|
|---|
| 73 | TArrayD step(step2);
|
|---|
| 74 |
|
|---|
| 75 | //
|
|---|
| 76 | // Be carefull: This is not thread safe
|
|---|
| 77 | //
|
|---|
| 78 | if (vinit.GetSize() != step.GetSize() ||
|
|---|
| 79 | vinit.GetSize() != limlo.GetSize() ||
|
|---|
| 80 | vinit.GetSize() != limup.GetSize() ||
|
|---|
| 81 | vinit.GetSize() != fix.GetSize())
|
|---|
| 82 | {
|
|---|
| 83 | *fLog << "CallMinuit: Parameter Size Mismatch" << endl;
|
|---|
| 84 | return kFALSE;
|
|---|
| 85 | }
|
|---|
| 86 |
|
|---|
| 87 |
|
|---|
| 88 |
|
|---|
| 89 | //..............................................
|
|---|
| 90 | // Set the maximum number of parameters
|
|---|
| 91 | //TMinuit *const save = gMinuit;
|
|---|
| 92 |
|
|---|
| 93 | TMinuit &minuit = *new TMinuit(vinit.GetSize());
|
|---|
| 94 | minuit.SetName("MinuitSupercuts");
|
|---|
| 95 |
|
|---|
| 96 |
|
|---|
| 97 | //..............................................
|
|---|
| 98 | // Set the print level
|
|---|
| 99 | // -1 no output except SHOW comands
|
|---|
| 100 | // 0 minimum output
|
|---|
| 101 | // 1 normal output (default)
|
|---|
| 102 | // 2 additional ouput giving intermediate results
|
|---|
| 103 | // 3 maximum output, showing progress of minimizations
|
|---|
| 104 | //
|
|---|
| 105 | minuit.SetPrintLevel(-1);
|
|---|
| 106 |
|
|---|
| 107 | //..............................................
|
|---|
| 108 | // Printout for warnings
|
|---|
| 109 | // SET WAR print warnings
|
|---|
| 110 | // SET NOW suppress warnings
|
|---|
| 111 | Int_t errWarn;
|
|---|
| 112 | Double_t tmpwar = 0;
|
|---|
| 113 | minuit.mnexcm("SET NOW", &tmpwar, 0, errWarn);
|
|---|
| 114 | //minuit.mnexcm("SET WAR", &tmpwar, 0, errWarn);
|
|---|
| 115 |
|
|---|
| 116 | //..............................................
|
|---|
| 117 | // Set the address of the minimization function
|
|---|
| 118 | minuit.SetFCN(fcn);
|
|---|
| 119 |
|
|---|
| 120 | //..............................................
|
|---|
| 121 | // Store address of object to be used in fcn
|
|---|
| 122 | minuit.SetObjectFit(objectfit );
|
|---|
| 123 |
|
|---|
| 124 | //..............................................
|
|---|
| 125 | // Set starting values and step sizes for parameters
|
|---|
| 126 | for (Int_t i=0; i<vinit.GetSize(); i++)
|
|---|
| 127 | {
|
|---|
| 128 | if (minuit.DefineParameter(i, name[i], vinit[i], step[i],
|
|---|
| 129 | limlo[i], limup[i]))
|
|---|
| 130 | {
|
|---|
| 131 | *fLog << "CallMinuit: Error in defining parameter "
|
|---|
| 132 | << name[i][0] << endl;
|
|---|
| 133 | return kFALSE;
|
|---|
| 134 | }
|
|---|
| 135 | }
|
|---|
| 136 |
|
|---|
| 137 | //..............................................
|
|---|
| 138 | //Int_t NumPars;
|
|---|
| 139 | //NumPars = minuit.GetNumPars();
|
|---|
| 140 | //*fLog << "CallMinuit : number of free parameters = "
|
|---|
| 141 | // << NumPars << endl;
|
|---|
| 142 |
|
|---|
| 143 |
|
|---|
| 144 | //..............................................
|
|---|
| 145 | // Error definition :
|
|---|
| 146 | //
|
|---|
| 147 | // for chisquare function :
|
|---|
| 148 | // up = 1.0 means calculate 1-standard deviation error
|
|---|
| 149 | // = 4.0 means calculate 2-standard deviation error
|
|---|
| 150 | //
|
|---|
| 151 | // for log(likelihood) function :
|
|---|
| 152 | // up = 0.5 means calculate 1-standard deviation error
|
|---|
| 153 | // = 2.0 means calculate 2-standard deviation error
|
|---|
| 154 | minuit.SetErrorDef(1.0);
|
|---|
| 155 |
|
|---|
| 156 | // Int_t errMigrad;
|
|---|
| 157 | // Double_t tmp = 0;
|
|---|
| 158 | // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);
|
|---|
| 159 |
|
|---|
| 160 | //..............................................
|
|---|
| 161 | // fix a parameter
|
|---|
| 162 | for (Int_t i=0; i<vinit.GetSize(); i++)
|
|---|
| 163 | {
|
|---|
| 164 | if (fix[i] > 0)
|
|---|
| 165 | {
|
|---|
| 166 | minuit.FixParameter(i);
|
|---|
| 167 | step[i] = 0.0;
|
|---|
| 168 | }
|
|---|
| 169 | }
|
|---|
| 170 |
|
|---|
| 171 | //..............................................
|
|---|
| 172 | //NumPars = minuit.GetNumPars();
|
|---|
| 173 | //*fLog << "CallMinuit : number of free parameters = "
|
|---|
| 174 | // << NumPars << endl;
|
|---|
| 175 |
|
|---|
| 176 | //..............................................
|
|---|
| 177 | // This doesn't seem to have any effect
|
|---|
| 178 | // Set maximum number of iterations (default = 500)
|
|---|
| 179 | //Int_t maxiter = 100000;
|
|---|
| 180 | //minuit.SetMaxIterations(maxiter);
|
|---|
| 181 |
|
|---|
| 182 | //..............................................
|
|---|
| 183 | // minimization by the method of Migrad
|
|---|
| 184 | if (method.Contains("Migrad", TString::kIgnoreCase))
|
|---|
| 185 | {
|
|---|
| 186 | //*fLog << "call MIGRAD" << endl;
|
|---|
| 187 | if (nulloutput)
|
|---|
| 188 | fLog->SetNullOutput(kTRUE);
|
|---|
| 189 | Double_t tmp = 0;
|
|---|
| 190 | minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
|
|---|
| 191 | if (nulloutput)
|
|---|
| 192 | fLog->SetNullOutput(kFALSE);
|
|---|
| 193 | //*fLog << "return from MIGRAD" << endl;
|
|---|
| 194 | }
|
|---|
| 195 |
|
|---|
| 196 | //..............................................
|
|---|
| 197 | // same minimization as by Migrad
|
|---|
| 198 | // but switches to the SIMPLEX method if MIGRAD fails to converge
|
|---|
| 199 | if (method.Contains("Minimize", TString::kIgnoreCase))
|
|---|
| 200 | {
|
|---|
| 201 | *fLog << "call MINIMIZE" << endl;
|
|---|
| 202 | Double_t tmp = 0;
|
|---|
| 203 | minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
|
|---|
| 204 | *fLog << "return from MINIMIZE" << endl;
|
|---|
| 205 | }
|
|---|
| 206 |
|
|---|
| 207 | //..............................................
|
|---|
| 208 | // minimization by the SIMPLEX method
|
|---|
| 209 | if (method.Contains("Simplex", TString::kIgnoreCase))
|
|---|
| 210 | {
|
|---|
| 211 | *fLog << "call SIMPLEX" << endl;
|
|---|
| 212 | if (nulloutput)
|
|---|
| 213 | fLog->SetNullOutput(kTRUE);
|
|---|
| 214 | Int_t maxcalls = 3000;
|
|---|
| 215 | Double_t tolerance = 0.1;
|
|---|
| 216 | Double_t tmp[2];
|
|---|
| 217 | tmp[0] = maxcalls;
|
|---|
| 218 | tmp[1] = tolerance;
|
|---|
| 219 | minuit.mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize);
|
|---|
| 220 | if (nulloutput)
|
|---|
| 221 | fLog->SetNullOutput(kFALSE);
|
|---|
| 222 | *fLog << "return from SIMPLEX" << endl;
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 | //..............................................
|
|---|
| 226 | // check quality of minimization
|
|---|
| 227 | // istat = 0 covariance matrix not calculated
|
|---|
| 228 | // 1 diagonal approximation only (not accurate)
|
|---|
| 229 | // 2 full matrix, but forced positive-definite
|
|---|
| 230 | // 3 full accurate covariance matrix
|
|---|
| 231 | // (indication of normal convergence)
|
|---|
| 232 | minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
|
|---|
| 233 |
|
|---|
| 234 | //if (fErrMinimize != 0 || fIstat < 3)
|
|---|
| 235 | if (fErrMinimize != 0)
|
|---|
| 236 | {
|
|---|
| 237 | *fLog << "CallMinuit : Minimization failed" << endl;
|
|---|
| 238 | *fLog << " fMin = " << fMin << ", fEdm = " << fEdm
|
|---|
| 239 | << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat
|
|---|
| 240 | << ", fErrMinimize = " << fErrMinimize << endl;
|
|---|
| 241 | return kFALSE;
|
|---|
| 242 | }
|
|---|
| 243 |
|
|---|
| 244 | //*fLog << "CallMinuit : Minimization was successful" << endl;
|
|---|
| 245 | //*fLog << " fMin = " << fMin << ", fEdm = " << fEdm
|
|---|
| 246 | // << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat
|
|---|
| 247 | // << ", fErrMinimize = " << fErrMinimize << endl;
|
|---|
| 248 |
|
|---|
| 249 |
|
|---|
| 250 | //..............................................
|
|---|
| 251 | // minimization by the method of Migrad
|
|---|
| 252 | if (method.Contains("Hesse", TString::kIgnoreCase))
|
|---|
| 253 | {
|
|---|
| 254 | //*fLog << "call HESSE" << endl;
|
|---|
| 255 | Double_t tmp = 0;
|
|---|
| 256 | minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize);
|
|---|
| 257 | //*fLog << "return from HESSE" << endl;
|
|---|
| 258 | }
|
|---|
| 259 |
|
|---|
| 260 | //..............................................
|
|---|
| 261 | // Minos error analysis
|
|---|
| 262 | if (method.Contains("Minos", TString::kIgnoreCase))
|
|---|
| 263 | {
|
|---|
| 264 | //*fLog << "call MINOS" << endl;
|
|---|
| 265 | Double_t tmp = 0;
|
|---|
| 266 | minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize);
|
|---|
| 267 | //*fLog << "return from MINOS" << endl;
|
|---|
| 268 | }
|
|---|
| 269 |
|
|---|
| 270 | //..............................................
|
|---|
| 271 | // Print current status of minimization
|
|---|
| 272 | // if nkode = 0 only function value
|
|---|
| 273 | // 1 parameter values, errors, limits
|
|---|
| 274 | // 2 values, errors, step sizes, internal values
|
|---|
| 275 | // 3 values, errors, step sizes, 1st derivatives
|
|---|
| 276 | // 4 values, parabolic errors, MINOS errors
|
|---|
| 277 |
|
|---|
| 278 | //Int_t nkode = 4;
|
|---|
| 279 | //minuit.mnprin(nkode, fmin);
|
|---|
| 280 |
|
|---|
| 281 | //..............................................
|
|---|
| 282 | // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
|
|---|
| 283 | // iflag = 1 initial calculations only
|
|---|
| 284 | // 2 calculate 1st derivatives and function
|
|---|
| 285 | // 3 calculate function only
|
|---|
| 286 | // 4 calculate function + final calculations
|
|---|
| 287 | Double_t iflag = 3;
|
|---|
| 288 | Int_t errfcn3;
|
|---|
| 289 | minuit.mnexcm("CALL", &iflag, 1, errfcn3);
|
|---|
| 290 |
|
|---|
| 291 | // WW : the following statements were commented out because the
|
|---|
| 292 | // Minuit object will still be used;
|
|---|
| 293 | // this may be changed in the future
|
|---|
| 294 | //delete &minuit;
|
|---|
| 295 | //gMinuit = save;
|
|---|
| 296 |
|
|---|
| 297 | return kTRUE;
|
|---|
| 298 | }
|
|---|
| 299 |
|
|---|
| 300 |
|
|---|
| 301 |
|
|---|
| 302 |
|
|---|
| 303 |
|
|---|
| 304 |
|
|---|
| 305 |
|
|---|
| 306 |
|
|---|
| 307 |
|
|---|
