Index: trunk/MagicSoft/Mars/mbase/MMath.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MMath.h (revision 3581) +++ trunk/MagicSoft/Mars/mbase/MMath.h (revision 3582) @@ -9,6 +9,6 @@ { public: - Double_t Significance(Double_t s, Double_t b); - Double_t SignificanceLiMa(Double_t s, Double_t b, Double_t alpha=1); + static Double_t Significance(Double_t s, Double_t b); + static Double_t SignificanceLiMa(Double_t s, Double_t b, Double_t alpha=1); ClassDef(MMath, 0) Index: trunk/MagicSoft/Mars/mfit/MTFitLoop.cc =================================================================== --- trunk/MagicSoft/Mars/mfit/MTFitLoop.cc (revision 3582) +++ trunk/MagicSoft/Mars/mfit/MTFitLoop.cc (revision 3582) @@ -0,0 +1,230 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 3/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MTFitLoop +// +// This class optimized parameters which are calculated in an eventloop. +// For this it is unimportant whether the loop reads from a file a +// matrix or from somewhere else... +// +// The parameters which are optimized must be stored in a MParContainer +// which overwrites SetVariables(). Eg. a MDataChain or MF overwrites +// SetVariables(). In a MF all arguments given as [0], [1] are +// set by SetVariables (in MDataValue). +// eg: In you loop you have a cut like this: +// MF filter("MHillas.fWidth<[0]"); +// filter.SetName("MyParameters"); +// +// Now for each time the eventloop is executed +// (MEvtLoop::Eventloop(fNumEvents)) the parameters from TMinuit are +// passed to MF::SetVariables and changed. +// +// +///////////////////////////////////////////////////////////////////////////// +#include "MTFitLoop.h" + +#include +#include + +#include "MParList.h" +#include "MRead.h" +#include "MEvtLoop.h" + +#include "MLog.h" +#include "MLogManip.h" + +//------------------------------------------------------------------------ +// +// fcn calculates the function to be minimized (using TMinuit::Migrad) +// +void MTFitLoop::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + MTFitLoop *optim = (MTFitLoop*)gMinuit->GetObjectFit(); + + TMinuit *minuit = gMinuit; + f = optim->Fcn(npar, gin, par, iflag); + gMinuit = minuit; + +} +Double_t MTFitLoop::Fcn(Int_t &npar, Double_t *gin, Double_t *par, Int_t iflag) +{ + MParList *plist = fEvtLoop->GetParList(); + + MParameterD *eval = (MParameterD*)plist->FindCreateObj("MParameterD", fFitParameter); + MParContainer *pars = plist->FindObject(fParametersName); + + MRead *read = (MRead*)plist->FindObject("MTaskList")->FindObject("MRead"); + if (read) + read->Rewind(); + + if (fDebug>=1) + { + Double_t fmin, fedm, errdef; + Int_t n1, n2, istat; + gMinuit->mnstat(fmin, fedm, errdef, n1, n2, istat); + *fLog << inf << underline << "Minimization Status so far:" << endl; + *fLog << " Calls: " << gMinuit->fNfcn << endl; + *fLog << " Func min: " << fmin << endl; + *fLog << " Found edm: " << fedm << endl; + *fLog << " ErrDef: " << errdef << endl; + *fLog << " Status: "; + switch (istat) + { + case 0: *fLog << "n/a" << endl; break; + case 1: *fLog << "approximation only, not accurate" << endl; break; + case 2: *fLog << "full matrix, but forced positive-definite" << endl; break; + case 3: *fLog << "full accurate covariance matrix" << endl; break; + default: *fLog << "undefined" << endl; break; + } + } + + if (fDebug>=0) + { + *fLog << inf << "Set: "; + for (Int_t i=0; i<7; i++) + *fLog << par[i] << " "; + *fLog << endl; + } + + pars->SetVariables(TArrayD(gMinuit->fMaxpar, par)); + + if (fDebug<3) + gLog.SetNullOutput(kTRUE); + fEvtLoop->Eventloop(fNumEvents); + if (fDebug<3) + gLog.SetNullOutput(kFALSE); + + const Double_t f = eval->GetVal(); + + if (fDebug>=0) + *fLog << inf << "F=" << f << endl; + + if (fDebug>=1) + fEvtLoop->GetTaskList()->PrintStatistics(); + + return f; +} + +MTFitLoop::MTFitLoop() +{ + fDebug = -1; + fNumEvents = -1; +} + +void MTFitLoop::Optimize(MEvtLoop &loop) +{ + *fLog << inf << "Event loop was setup" << endl; + MParList *parlist = loop->GetParList(); + if (!parlist) + return; + + MParContainer *pars = plist->FindObject(fParametersName); + if (!pars) + return; + + fEvtLoop = evtloop; + + MParContainer ¶meters = *pars; + + TMinuit *minsave = gMinuit; + + gMinuit = new TMinuit(parameters.GetSize()); + gMinuit->SetPrintLevel(-1); + + gMinuit->SetFCN(fcn); + gMinuit->SetObjectFit(this); + + // For chisq fits call this // seems to be something like %) + // + // The default tolerance is 0.1, and the minimization will stop"); + // when the estimated vertical distance to the minimum (EDM) is"); + // less than 0.001*[tolerance]*UP (see [SET ERRordef])."); + // + if (gMinuit->SetErrorDef(1000)) + { + *fLog << err << dbginf << "SetErrorDef failed." << endl; + return; + } + + // + // Set starting values and step sizes for parameters + // + for (Int_t i=0; iDefineParameter(i, name, vinit, step, 0, limup); + if (!rc) + continue; + + *fLog << err << dbginf << "Error in defining parameter #" << i << endl; + return; + } + + for (int i=0; iFixParameter(i); + else + gMinuit->Release(i); + + // Now ready for minimization step: + + TStopwatch clock; + clock.Start(); + const Bool_t rc = gMinuit->Migrad(); + clock.Stop(); + clock.Print(); + + if (rc) + { + *fLog << err << dbginf << "Migrad failed." << endl; + return; + } + + *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl; + + // + // Update values of fA, fB: + // + for (Int_t i=0; iGetParameter(i,x1,x2); + parameters[i] = x1; + cout << i << ": " << parameters[i] << endl; + } + + list.SetVariables(parameters); + + gMinuit = minsave; +} Index: trunk/MagicSoft/Mars/mfit/MTFitLoop.h =================================================================== --- trunk/MagicSoft/Mars/mfit/MTFitLoop.h (revision 3582) +++ trunk/MagicSoft/Mars/mfit/MTFitLoop.h (revision 3582) @@ -0,0 +1,46 @@ +#ifndef MARS_MTFitLoop +#define MARS_MTFitLoop + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TArrayC +#include +#endif + +class MEvtLoop; + +class MTFitLoop : public MParContainer +{ +private: + Int_t fDebug; // -1 no output, 0 MTFitLoop output, 1 PrintStatistics output + Int_t fNumEvents; + + static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); + virtual Double_t Fcn(Int_t &npar, Double_t *gin, Double_t *par, Int_t iflag); + + MEvtLoop *fEvtLoop; + + TString fParametersName; + TString fFitParameter; + + TArrayC fFixedParams; + +public: + MTFitLoop(); + + void Optimize(MEvtLoop &loop); + + void SetNameParameters(const char *parm) { fParametersName = parm; } + void SetFitParameter(const char *parm) { fFitParameter = parm; } + + void SetFixedParameters(const TArrayC &c) { fFixedParams = c; } + + void SetDebug(Int_t n) { fDebug = n; } + void SetNumEvents(Int_t n) { fNumEvents = n; } + + ClassDef(MTFitLoop, 0) // Class which can optimize a value (chi^2, significance, etc) calculated in an eventloop +}; + +#endif