Changeset 2313


Ignore:
Timestamp:
08/22/03 10:33:25 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/manalysis
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc

    r2300 r2313  
    3636#include <math.h>            // fabs
    3737
     38#include <TArrayD.h>
     39#include <TArrayI.h>
    3840#include <TMinuit.h>
    3941#include <TStopwatch.h>
    … …  
    4244#include "MLogManip.h"
    4345#include "MParContainer.h"
    44 
    4546
    4647ClassImp(MMinuitInterface);
    … …  
    5758    fName  = name  ? name  : "MMinuitInterface";
    5859    fTitle = title ? title : "Interface for Minuit";
    59 
    6060}
    61 
    62 // --------------------------------------------------------------------------
    63 //
    64 // Default destructor.
    65 //
    66 MMinuitInterface::~MMinuitInterface()
    67 {
    68 }
    69 
    7061
    7162// -----------------------------------------------------------------------
    … …  
    7667Bool_t MMinuitInterface::CallMinuit(
    7768             void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
    78              UInt_t npar, char name[80][100],
    79              Double_t vinit[80], Double_t step[80],
    80              Double_t limlo[80], Double_t limup[80], Int_t fix[80],
    81              TObject *objectfit , TString method, Bool_t nulloutput)
     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)
    8272{
    83         TString str(method);
    84 
    85         //
    86         // Be carefull: This is not thread safe
    87         //
    88         UInt_t maxpar = 100;
    89 
    90         if (npar > maxpar)
     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]))
    91130        {
    92             *fLog << "CallMinuit : too many parameters,  npar = " << npar
    93                  << ",   maxpar = " << maxpar << endl;
     131            *fLog << "CallMinuit: Error in defining parameter "
     132                << name[i][0] << endl;
    94133            return kFALSE;
    95134        }
    96 
    97 
    98 
    99         //..............................................
    100         // Set the maximum number of parameters
    101         TMinuit &minuit = *(new TMinuit(maxpar));
    102         minuit.SetName("MinuitSupercuts");
    103         //TMinuit minuit(maxpar);
    104         //gMinuit = new TMinuit(maxpar);
    105         //TMinuit &minuit = *gMinuit;
    106 
    107 
    108         //..............................................
    109         // Set the print level
    110         // -1   no output except SHOW comands
    111         //  0   minimum output
    112         //  1   normal output (default)
    113         //  2   additional ouput giving intermediate results
    114         //  3   maximum output, showing progress of minimizations
    115         //
    116         //Int_t printLevel = -1;
    117         Int_t printLevel = 1;
    118         minuit.SetPrintLevel(printLevel);
    119 
    120         //..............................................       
    121         // Printout for warnings
    122         //    SET WAR      print warnings
    123         //    SET NOW      suppress warnings
    124         Int_t errWarn;
    125         Double_t tmpwar = 0;
    126         minuit.mnexcm("SET NOW", &tmpwar, 0, errWarn);
    127         //minuit.mnexcm("SET WAR", &tmpwar, 0, errWarn);
    128 
    129         //..............................................
    130         // Set the address of the minimization function
    131         minuit.SetFCN(fcn);
    132 
    133 
    134         //..............................................
    135         // Store address of object to be used in fcn
    136         minuit.SetObjectFit(objectfit );
    137 
    138 
    139         //..............................................
    140         // Set starting values and step sizes for parameters
    141         for (UInt_t i=0; i<npar; i++)
     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)
    142165        {
    143             if (minuit.DefineParameter(i, &name[i][0], vinit[i], step[i],
    144                                        limlo[i], limup[i]))
    145             {
    146                 *fLog << "CallMinuit: Error in defining parameter "
    147                      << name[i][0] << endl;
    148                 return kFALSE;
    149             }
     166            minuit.FixParameter(i);
     167            step[i] = 0.0;
    150168        }
    151 
    152         //..............................................
    153         //Int_t NumPars;
    154         //NumPars = minuit.GetNumPars();
    155         //*fLog << "CallMinuit :  number of free parameters = "
    156         //     << NumPars << endl;
    157 
    158 
    159         //..............................................
    160         // Error definition :
    161         //
    162         //    for chisquare function :
    163         //      up = 1.0   means calculate 1-standard deviation error
    164         //         = 4.0   means calculate 2-standard deviation error
    165         //
    166         //    for log(likelihood) function :
    167         //      up = 0.5   means calculate 1-standard deviation error
    168         //         = 2.0   means calculate 2-standard deviation error
    169         Double_t up = 1.0;
    170         minuit.SetErrorDef(up);
    171 
    172 
    173 
    174         // Int_t errMigrad;
    175         // Double_t tmp = 0;
    176         // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);
    177 
    178 
    179         //..............................................
    180         // fix a parameter
    181         for (UInt_t i=0; i<npar; i++)
    182         {
    183             if (fix[i] > 0)
    184             {
    185                 Int_t parNo = i;
    186                 minuit.FixParameter(parNo);
    187                 step[parNo] = 0.0;
    188             }
    189         }
    190 
    191         //..............................................
    192         //NumPars = minuit.GetNumPars();
    193         //*fLog << "CallMinuit :  number of free parameters = "
    194         //     << NumPars << endl;
    195 
    196         //..............................................
    197         // Set maximum number of iterations (default = 500)
    198         //Int_t maxiter = 100000;
    199         Int_t maxiter = 10;
    200         minuit.SetMaxIterations(maxiter);
    201 
    202         //..............................................
    203         // minimization by the method of Migrad
    204         if (str.Contains("Migrad", TString::kIgnoreCase))       
    205         {
    206           //*fLog << "call MIGRAD" << endl;
    207           if (nulloutput)
     169    }
     170
     171    //..............................................
     172    //NumPars = minuit.GetNumPars();
     173    //*fLog << "CallMinuit :  number of free parameters = "
     174    //     << NumPars << endl;
     175
     176    //..............................................
     177    // Set maximum number of iterations (default = 500)
     178    //Int_t maxiter = 100000;
     179    minuit.SetMaxIterations(10);
     180
     181    //..............................................
     182    // minimization by the method of Migrad
     183    if (method.Contains("Migrad", TString::kIgnoreCase))
     184    {
     185        //*fLog << "call MIGRAD" << endl;
     186        if (nulloutput)
    208187            fLog->SetNullOutput(kTRUE);
    209           Double_t tmp = 0;
    210           minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
    211           if (nulloutput)
     188        Double_t tmp = 0;
     189        minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
     190        if (nulloutput)
    212191            fLog->SetNullOutput(kFALSE);
    213           //*fLog << "return from MIGRAD" << endl;
    214         }
    215 
    216         //..............................................       
    217         // same minimization as by Migrad
    218         // but switches to the SIMPLEX method if MIGRAD fails to converge
    219         if (str.Contains("Minimize", TString::kIgnoreCase))       
    220         {
    221           *fLog << "call MINIMIZE" << endl;
    222           Double_t tmp = 0;
    223           minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
    224           *fLog << "return from MINIMIZE" << endl;
    225         }
    226 
    227         //..............................................       
    228         // minimization by the SIMPLEX method
    229         if (str.Contains("Simplex", TString::kIgnoreCase))       
    230         {
    231           *fLog << "call SIMPLEX" << endl;
    232           if (nulloutput)
     192        //*fLog << "return from MIGRAD" << endl;
     193    }
     194
     195    //..............................................
     196    // same minimization as by Migrad
     197    // but switches to the SIMPLEX method if MIGRAD fails to converge
     198    if (method.Contains("Minimize", TString::kIgnoreCase))
     199    {
     200        *fLog << "call MINIMIZE" << endl;
     201        Double_t tmp = 0;
     202        minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
     203        *fLog << "return from MINIMIZE" << endl;
     204    }
     205
     206    //..............................................
     207    // minimization by the SIMPLEX method
     208    if (method.Contains("Simplex", TString::kIgnoreCase))
     209    {
     210        *fLog << "call SIMPLEX" << endl;
     211        if (nulloutput)
    233212            fLog->SetNullOutput(kTRUE);
    234           Double_t tmp = 0;
    235           minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize);
    236           if (nulloutput)
     213        Double_t tmp = 0;
     214        minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize);
     215        if (nulloutput)
    237216            fLog->SetNullOutput(kFALSE);
    238           //*fLog << "return from SIMPLEX" << endl;
    239         }
    240 
    241         //..............................................       
    242         // check quality of minimization
    243         // istat = 0   covariance matrix not calculated
    244         //         1   diagonal approximation only (not accurate)
    245         //         2   full matrix, but forced positive-definite
    246         //         3   full accurate covariance matrix
    247         //             (indication of normal convergence)
    248         minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
    249 
    250         //if (fErrMinimize != 0  ||  fIstat < 3)
    251         if (fErrMinimize != 0)
    252         {
    253             *fLog << "CallMinuit : Minimization failed" << endl;
    254             *fLog << "       fMin = " << fMin   << ",   fEdm = "  << fEdm
    255                  << ",   fErrdef = "  << fErrdef << ",   fIstat = " << fIstat
    256                  << ",   fErrMinimize = " << fErrMinimize << endl;
    257             return kFALSE;
    258         }
    259 
    260         //*fLog << "CallMinuit : Minimization was successful" << endl;
    261         //*fLog << "       fMin = " << fMin   << ",   fEdm = "  << fEdm
    262         //     << ",   fErrdef = "  << fErrdef << ",   fIstat = " << fIstat
    263         //     << ",   fErrMinimize = " << fErrMinimize << endl;
    264 
    265 
    266         //..............................................
    267         // minimization by the method of Migrad
    268         if (str.Contains("Hesse", TString::kIgnoreCase))       
    269         {
    270           //*fLog << "call HESSE" << endl;
    271           Double_t tmp = 0;
    272           minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize);
    273           //*fLog << "return from HESSE" << endl;
    274         }
    275 
    276         //..............................................
    277         // Minos error analysis
    278         if (str.Contains("Minos", TString::kIgnoreCase))       
    279         {
    280           //*fLog << "call MINOS" << endl;
    281           Double_t tmp = 0;
    282           minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize);
    283           //*fLog << "return from MINOS" << endl;
    284         }
    285 
    286         //..............................................       
    287         // Print current status of minimization
    288         // if nkode = 0    only function value
    289         //            1    parameter values, errors, limits
    290         //            2    values, errors, step sizes, internal values
    291         //            3    values, errors, step sizes, 1st derivatives
    292         //            4    values, paraboloc errors, MINOS errors
    293  
    294         //Int_t nkode = 4;
    295         //minuit.mnprin(nkode, fmin);
    296 
    297         //..............................................       
    298         // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
    299         // iflag = 1   initial calculations only
    300         //         2   calculate 1st derivatives and function
    301         //         3   calculate function only
    302         //         4   calculate function + final calculations
    303         const char *command = "CALL";
    304         Double_t iflag = 3;
    305         Int_t errfcn3;
    306         minuit.mnexcm(command, &iflag, 1, errfcn3);
    307 
    308         //*fLog << "Exit CallMinuit" << endl;
    309 
    310         return kTRUE;
     217        //*fLog << "return from SIMPLEX" << endl;
     218    }
     219
     220    //..............................................
     221    // check quality of minimization
     222    // istat = 0   covariance matrix not calculated
     223    //         1   diagonal approximation only (not accurate)
     224    //         2   full matrix, but forced positive-definite
     225    //         3   full accurate covariance matrix
     226    //             (indication of normal convergence)
     227    minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
     228
     229    //if (fErrMinimize != 0  ||  fIstat < 3)
     230    if (fErrMinimize != 0)
     231    {
     232        *fLog << "CallMinuit : Minimization failed" << endl;
     233        *fLog << "       fMin = " << fMin   << ",   fEdm = "  << fEdm
     234            << ",   fErrdef = "  << fErrdef << ",   fIstat = " << fIstat
     235            << ",   fErrMinimize = " << fErrMinimize << endl;
     236        return kFALSE;
     237    }
     238
     239    //*fLog << "CallMinuit : Minimization was successful" << endl;
     240    //*fLog << "       fMin = " << fMin   << ",   fEdm = "  << fEdm
     241    //     << ",   fErrdef = "  << fErrdef << ",   fIstat = " << fIstat
     242    //     << ",   fErrMinimize = " << fErrMinimize << endl;
     243
     244
     245    //..............................................
     246    // minimization by the method of Migrad
     247    if (method.Contains("Hesse", TString::kIgnoreCase))
     248    {
     249        //*fLog << "call HESSE" << endl;
     250        Double_t tmp = 0;
     251        minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize);
     252        //*fLog << "return from HESSE" << endl;
     253    }
     254
     255    //..............................................
     256    // Minos error analysis
     257    if (method.Contains("Minos", TString::kIgnoreCase))
     258    {
     259        //*fLog << "call MINOS" << endl;
     260        Double_t tmp = 0;
     261        minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize);
     262        //*fLog << "return from MINOS" << endl;
     263    }
     264
     265    //..............................................
     266    // Print current status of minimization
     267    // if nkode = 0    only function value
     268    //            1    parameter values, errors, limits
     269    //            2    values, errors, step sizes, internal values
     270    //            3    values, errors, step sizes, 1st derivatives
     271    //            4    values, paraboloc errors, MINOS errors
     272
     273    //Int_t nkode = 4;
     274    //minuit.mnprin(nkode, fmin);
     275
     276    //..............................................
     277    // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
     278    // iflag = 1   initial calculations only
     279    //         2   calculate 1st derivatives and function
     280    //         3   calculate function only
     281    //         4   calculate function + final calculations
     282    Double_t iflag = 3;
     283    Int_t errfcn3;
     284    minuit.mnexcm("CALL", &iflag, 1, errfcn3);
     285
     286    delete &minuit;
     287    gMinuit = save;
     288
     289    return kTRUE;
    311290}
    312 
    313 
    314 //===========================================================================
    315 
    316 
    317 
    318 
    319 
    320 
    321 
    322 
    323 
    324 
    325 

  • trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h

    r2300 r2313  
    22#define MARS_MMinuitInterface
    33
    4 #ifndef MARS_MTask
    5 #include "MTask.h"
     4#ifndef MARS_MParContainer
     5#include "MParContainer.h"
    66#endif
    77
    8 #ifndef ROOT_TArrayD
    9 #include <TArrayD.h>
    10 #endif
     8class TArrayD;
     9class TArrayI;
    1110
    12 
    13 
    14 class MMinuitInterface : public MTask
     11class MMinuitInterface : public MParContainer
    1512{
    1613private:
    17 
    18   UInt_t   fNpar;
    19   Double_t fMin,   fEdm,   fErrdef;
    20   Int_t    fNpari, fNparx, fIstat;
    21   Int_t    fErrMinimize;
    22 
     14    UInt_t   fNpar;
     15    Double_t fMin,   fEdm,   fErrdef;
     16    Int_t    fNpari, fNparx, fIstat;
     17    Int_t    fErrMinimize;
    2318
    2419public:
    25   MMinuitInterface(const char *name=NULL, const char *title=NULL);
    26   ~MMinuitInterface();
     20    MMinuitInterface(const char *name=NULL, const char *title=NULL);
    2721
    28 Bool_t MMinuitInterface::CallMinuit(
    29               void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
    30               UInt_t npar, char name[80][100],
    31               Double_t vinit[80], Double_t step[80],
    32               Double_t limlo[80], Double_t limup[80], Int_t fix[80],
    33               TObject *fObjectFit, TString method, Bool_t nulloutput);
     22    Bool_t CallMinuit(
     23                      void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
     24                      const TString *name, const TArrayD &vinit, const TArrayD &step,
     25                      const TArrayD &limlo, const TArrayD &limup, const TArrayI &fix,
     26                      TObject *fObjectFit, const TString &method, Bool_t nulloutput);
    3427
    35   ClassDef(MMinuitInterface, 1) // Class for interfacing with Minuit
     28    ClassDef(MMinuitInterface, 0) // Class for interfacing with Minuit
    3629};
    3730
Note: See TracChangeset for help on using the changeset viewer.