Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7151) +++ trunk/MagicSoft/Mars/Changelog (revision 7152) @@ -33,4 +33,54 @@ - improved output - added plotting other spectras + + * Makefile: + - added new directory mjoptim + + * manalysis/MEnergyEstimate.[h,cc]: + - derives now from a new more genral class MParameterCalc + + * manalysis/Makefile, manalysis/AnalysisLinkDef.h: + - added new class MParamterCalc + + * manalysis/MParameterCalc.[h,cc]: + - added new class MParamterCalc + + * mbase/MLog.cc: + - do not display 'underline' if fIsNull output + + * mhbase/MBinning.cc: + - removed some obsolete comments fixed old stylish SetBinningCos + + * mhflux/MAlphaFitter.[h,cc]: + - background fit can now be switched off + - removed some obsolete comments + - adapted Print output + - added new option to change signal function type + - adde new Strategy option: kGaussSigma + + * mhflux/MHAlpha.[h,cc]: + - output the gauss sigma to the parameter list + + * mhflux/MHThetaSq.h: + - mase InitMapping public (used in MJOptimizeDisp) + + * mjobs/JobsLinkDef.h, mjobs/Makefile: + - removed MJOptimize* + + * mjobs/MJOptimize.[h,cc]: + - allow minimization of something else than MinimizationValue + + * mjobs/MJOptimizeEnergy.[h,cc]: + - added class description with example + - added weighting task to RunEnergy (untested) + + * mtools/MChisqEval.[h,cc]: + - added comments + - print result in Postprocess + - allow to give a different name to the result + + * mjobs/MJOptimize*: + - finally moved to new dir mjoptim + - new class MJOptimizeDisp Index: trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 7151) +++ trunk/MagicSoft/Mars/mhflux/MAlphaFitter.cc (revision 7152) @@ -380,5 +380,5 @@ return rc; } - +/* Bool_t MAlphaFitter::FitTime(const TH3D &hon, UInt_t bin, Bool_t paint) { @@ -395,5 +395,5 @@ return rc; } - +*/ Bool_t MAlphaFitter::FitAlpha(const TH3D &hon, Bool_t paint) { Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 7151) +++ (revision ) @@ -1,808 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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, 9/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MJOptimize -// -// Class for otimizing the parameters of the supercuts -// -// Minimization Control -// ==================== -// -// To choose the minimization algorithm use: -// void SetOptimizer(Optimizer_t o); -// -// Allowed options are: -// enum Optimizer_t -// { -// kMigrad, // Minimize by the method of Migrad -// kSimplex, // Minimize by the method of Simplex -// kMinimize, // Migrad + Simplex (if Migrad fails) -// kMinos, // Minos error determination -// kImprove, // Local minimum search -// kSeek, // Minimize by the method of Monte Carlo -// kNone // Skip optimization -// }; -// -// For more details on the methods see TMinuit. -// -// -// You can change the behaviour of the minimization using -// -// void SetNumMaxCalls(UInt_t num=0); -// void SetTolerance(Float_t tol=0); -// -// While NumMaxCalls is the first, Tolerance the second arguement. -// For more details start root and type -// -// gMinuit->mnhelp("command") -// -// while command can be -// * MIGRAD -// * SIMPLEX -// * MINIMIZE -// * MINOS -// * IMPROVE -// * SEEK -// -// The default (num==0 and tol==0) should always give you the -// corresponding defaults used in Minuit. -// -// -// FIXME: Implement changing cut in hadronness... -// FIXME: Show MHSignificance on MStatusDisplay during filling... -// FIXME: Choose step-size percentage as static data membewr -// FIXME: Choose minimization method -// -///////////////////////////////////////////////////////////////////////////// -#include "MJOptimize.h" - -#include -#include - -#include - -#include - -#include -#include - -#include "MHMatrix.h" - -// environment -#include "MLog.h" -#include "MLogManip.h" - -#include "MDirIter.h" -#include "MStatusDisplay.h" - -// eventloop -#include "MParList.h" -#include "MTaskList.h" -#include "MEvtLoop.h" - -// parameters -#include "MParameters.h" - -// tasks -#include "MReadTree.h" -#include "MMatrixLoop.h" -#include "MFillH.h" - -// filters -#include "MF.h" -#include "MFilterList.h" - -using namespace std; - -//------------------------------------------------------------------------ -// -// fcn calculates the function to be minimized (using TMinuit::Migrad) -// -void MJOptimize::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - MJOptimize *optim = (MJOptimize*)gMinuit->GetObjectFit(); - - // WORKAROUND --- FOR WHAT? - if (gMinuit->fEpsi<1e-2) - { - *optim->fLog << warn << "WARNING - For unknown reasons: fEspi<1e-100... resetting to 0.01." << endl; - gMinuit->fEpsi = 0.01; - } - - TMinuit *minuit = gMinuit; - f = optim->Fcn(TArrayD(TMath::Min(gMinuit->fMaxpar, optim->fParameters.GetSize()), par), minuit); - gMinuit = minuit; - -} - -Double_t MJOptimize::Fcn(const TArrayD &par, TMinuit *minuit) -{ - if (fEvtLoop->GetDisplay()!=fDisplay) - return 0; - /* - switch(iflag) - { - case 1: // first call - case 2: // calc derivative - break; - case 3: - // last call - MStatusDisplay *d = new MStatusDisplay; - fEvtLoop->SetDisplay(d); - break; - } - */ - MParList *plist = fEvtLoop->GetParList(); - - MParameterD *eval = (MParameterD*)plist->FindObject(fNameMinimizationValue, "MParameterD"); - MParContainer *pars = (MParContainer*)plist->FindObject("MParameters", "MParContainer"); - - MRead *read = (MRead*)plist->FindObject("MTaskList")->FindObject("MRead"); - if (read) - read->Rewind(); - - if (fDebug>=0) - { - *fLog << inf << "New Set: "; - for (Int_t i=0; iSetVariables(par); - eval->SetVal(0); - - const Bool_t isnull = gLog.IsNullOutput(); - if (fDebug<3) - gLog.SetNullOutput(kTRUE); - - TStopwatch clock; - clock.Start(); - fEvtLoop->Eventloop(fNumEvents, MEvtLoop::kNoStatistics); - clock.Stop(); - - if (fDebug<3) - gLog.SetNullOutput(isnull); - - const Double_t f = eval->GetVal(); - - if (fDebug>=0) - *fLog << inf << "Result F=" << f << endl; - - if (fDebug>=1 && minuit) - { - Double_t fmin, fedm, errdef; - Int_t n1, n2, istat; - minuit->mnstat(fmin, fedm, errdef, n1, n2, istat); - *fLog << inf << underline << "Minimization Status so far:" << endl; - *fLog << " Calls: " << minuit->fNfcn << " (max=" << gMinuit->fMaxIterations << ")" << endl; - *fLog << " Parameters: fixed=" << gMinuit->fNpfix << ", free=" << gMinuit->fNpar << endl; - *fLog << " Func min: " << fmin << " (Epsi=" << gMinuit->fEpsi << ", Apsi=" << gMinuit->fApsi << ")" << 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>=1) - { - clock.Print(); - fEvtLoop->GetTaskList()->PrintStatistics(); - } - - return f; -} - -MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0), fTestTrain(0), fNameMinimizationValue("MinimizationValue") -{ - fRules.SetOwner(); - fFilter.SetOwner(); - - fNamesOn.SetOwner(); - fNamesOff.SetOwner(); -} - -//------------------------------------------------------------------------ -// -// Add seqeunces from list to reader -// -Bool_t MJOptimize::AddSequences(MRead &read, TList &list) const -{ - MDirIter files; - - TIter Next(&list); - TObject *o=0; - while ((o=Next())) - { - MSequence seq(o->GetName()); - if (!seq.IsValid()) - return kFALSE; - - seq.SetupDatRuns(files, o->GetTitle(), "I"); - } - - return read.AddFiles(files)>0; -} - -//------------------------------------------------------------------------ -// -// Add on-sequences: -// - fname: sequence file name (with path) -// - dir: directory were image files are stored -// -void MJOptimize::AddSequenceOn(const char *fname, const char *dir) -{ - fNamesOn.Add(new TNamed(fname, dir)); -} - -//------------------------------------------------------------------------ -// -// Add off-sequences: -// - fname: sequence file name (with path) -// - dir: directory were image files are stored -// -void MJOptimize::AddSequenceOff(const char *fname, const char *dir) -{ - fNamesOff.Add(new TNamed(fname, dir)); -} - -//------------------------------------------------------------------------ -// -// Empty list of on- and off-sequences -// -void MJOptimize::ResetSequences() -{ - fNamesOn.Delete(); - fNamesOff.Delete(); -} - -//------------------------------------------------------------------------ -// -// Add a parameter used in your filters (see AddFilter) The parameter -// index is returned, -// -// Int_t idx = AddParameter("log10(MHillas.fSize)"); -// -// The indices area starting with 0 always. -// -Int_t MJOptimize::AddParameter(const char *rule) -{ - fRules.Add(new TNamed(rule, "")); - return fRules.GetSize()-1; -} - -//------------------------------------------------------------------------ -// -// Add a filter which can be applied in the optimization (for deatils -// see correspodning Run function) You can use the indices you got by -// AddParameter, eg -// -// AddFilter("M[0] < 3.2"); -// -// if used in optimization you can do -// -// AddFilter("M[0] < [0]"); -// -// for more details, see SetParameter and FixParameter -// -void MJOptimize::AddFilter(const char *rule) -{ - fFilter.Add(new MF(rule)); -} - -//------------------------------------------------------------------------ -// -// Add a cut which is used to fill the matrix, eg "MMcEvt.fOartId<1.5" -// (The rule is applied, nit inverted: The matrix is filled with -// the events fullfilling the condition) -// -void MJOptimize::AddPreCut(const char *rule) -{ - MFilter *f = new MF(rule); - f->SetBit(kCanDelete); - AddPreCut(f); -} - -//------------------------------------------------------------------------ -// -// Add a cut which is used to fill the matrix. If kCanDelete is set -// MJOptimize takes the ownership. -// -void MJOptimize::AddPreCut(MFilter *f) -{ - fPreCuts.Add(f); -} - -//------------------------------------------------------------------------ -// -// Set the fParameters Array accoring to par. -// -void MJOptimize::SetParameters(const TArrayD &par) -{ - fParameters = par; -} - -//------------------------------------------------------------------------ -// -// Set the number of events processed by the eventloop. (Be carfull, -// if you are doing on-off analysis and you only process the first -// 1000 events which are on-events only the optimization may not work) -// -void MJOptimize::SetNumEvents(UInt_t n) -{ - fNumEvents=n; -} - -//------------------------------------------------------------------------ -// -// Set a debug level, which tells the optimization how much information -// is displayed about and in the running eventloop. -// -void MJOptimize::SetDebug(UInt_t n) -{ - fDebug=n; -} - -//------------------------------------------------------------------------ -// -// Set a optimization algorithm to be used. For more information see -// TMinuit. -// -// Available Algorithms are: -// kMigrad, // Minimize by the method of Migrad -// kSimplex, // Minimize by the method of Simplex -// kSeek // Minimize by the method of Monte Carlo -// -void MJOptimize::SetOptimizer(Optimizer_t o) -{ - fType = o; -} - -//------------------------------------------------------------------------ -// -// If a status didplay is set, search for tab "Optimizer". -// If not found, create it. -// In the tab search for TMultiGraph "Parameters". -// If not found create it. -// If empty create TGraphs. -// Check number of graphs vs. number of parameters. -// return TList with graphs. -// -TList *MJOptimize::GetPlots() const -{ - if (!fDisplay) - return NULL; - - TCanvas *c=fDisplay->GetCanvas("Optimizer"); - if (!c) - c = &fDisplay->AddTab("Optimizer"); - - TMultiGraph *mg = dynamic_cast(c->FindObject("Parameters")); - if (!mg) - mg = new TMultiGraph("Parameters", "Parameters of optimization"); - - TList *l = mg->GetListOfGraphs(); - if (!l) - { - const Int_t n = fParameters.GetSize(); - for (int i=0; iSetLineColor(kBlue); - mg->Add(g, ""); - AddPoint(mg->GetListOfGraphs(), i, i==n?1:fParameters[i]); - } - mg->SetBit(kCanDelete); - mg->Draw("al*"); - - l = mg->GetListOfGraphs(); - } - - return l->GetSize() == fParameters.GetSize()+1 ? l : NULL; -} - -//------------------------------------------------------------------------ -// -// Add a point with y=val as last point in idx-th Tgraph of list l. -// -void MJOptimize::AddPoint(TList *l, Int_t idx, Float_t val) const -{ - if (!l) - return; - - TGraph *gr = (TGraph*)l->At(idx); - gr->SetPoint(gr->GetN(), gr->GetN(), val); -} - -Int_t MJOptimize::Minuit(TMinuit &minuit, const char *cmd) const -{ - Int_t err; - Double_t tmp[2] = { fNumMaxCalls, fTolerance }; - minuit.mnexcm(cmd, tmp, 2, err); - - switch (err) - { - case 0: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm excuted normally." << endl; - break; - case 1: - *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command is blank... ignored." << endl; - break; - case 2: - *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command-line syntax error... ignored." << endl; - break; - case 3: - *fLog << warn << GetDescriptor() << " TMinuit::mnexcm unknown command... ignored." << endl; - break; - case 4: - *fLog << warn << GetDescriptor() << " TMinuit::mnexcm - Abnormal termination (eg Migrad not converged)" << endl; - break; - /* - case 5: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - Parameters requested." << endl; - break; - case 6: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET INPUT returned." << endl; - break; - case 7: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET TITLE returned." << endl; - break; - case 8: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET COVAR returned." << endl; - break; - case 9: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - reserved." << endl; - break; - case 10: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - END returned." << endl; - break; - case 11: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - EXIT or STOP returned." << endl; - break; - case 12: - *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - RETURN returned." << endl; - break;*/ - } - - return err; -} - -Bool_t MJOptimize::Optimize(MEvtLoop &evtloop) -{ - if (fParameters.GetSize()==0) - { - *fLog << err << GetDescriptor() << "::Optimize: ERROR - Sorry, no parameters defined." << endl; - return kFALSE; - } - - if (fType==kNone) - return kTRUE; - - gMinuit = new TMinuit(fParameters.GetSize()); - - gMinuit->SetFCN(fcn); - gMinuit->SetObjectFit(this); - gMinuit->SetPrintLevel(-1); // Don't print when DefineParameter - - // - // Set starting values and step sizes for parameters - // - for (Int_t i=0; iDefineParameter(i, name, vinit, step, limlo, limup); - if (rc) - { - *fLog << err << dbginf << "Error in defining parameter #" << i << endl; - return kFALSE; - } - - if (step==0) - gMinuit->FixParameter(i); - } - - gMinuit->SetPrintLevel(1); // Switch on pritning again - gMinuit->mnprin(1,0); // Print all parameters - - fEvtLoop = &evtloop; - - TList *g=GetPlots(); - - // Now ready for minimization step: - TStopwatch clock; - clock.Start(); - switch (fType) - { - case kSimplex: - Simplex(*gMinuit); - break; - case kMigrad: - Migrad(*gMinuit); - break; - case kMinimize: - Minimize(*gMinuit); - break; - case kMinos: - Minos(*gMinuit); - break; - case kImprove: - Improve(*gMinuit); - break; - case kSeek: - Seek(*gMinuit); - break; - case kNone: // Should never happen - return kFALSE; - } - clock.Stop(); - clock.Print(); - - if (evtloop.GetDisplay()!=fDisplay) - { - *fLog << inf << "Optimization aborted by user." << endl; - fDisplay = 0; - return kFALSE; - } - - *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl; - - // - // Update values of fA, fB: - // - for (Int_t i=0; iGetParameter(i,x1,x2); - fParameters[i] = x1; - cout << i << ": " << fParameters[i] << endl; - - AddPoint(g, i, x1); - } - AddPoint(g, fParameters.GetSize(), gMinuit->fAmin); - - delete gMinuit; - - return kTRUE; -} - -//------------------------------------------------------------------------ -// -// Optimize allows to use the optimizing by an eventloop based on -// some requirements. -// -// 1) The tasklist to be executed must have the name MTaskList and -// be an entry in the parameterlist. -// -// 2) The reading task (MReadMarsFile, MMatrixLoop) must have the name -// "MRead". If it derives from MRead Rewind() must be implemented, -// otherwise it must start reading from scratch each time its -// PreProcess is called. -// -// 3) The parameters to be optimized must be accessed through (currently) -// a single parameter container (MParContainer) called MParameters. -// The parameters are set through SetVariables. -// -// 4) The result of a single function call for minimization (eg. chisquare) -// must be available after the eventloop in a container of type -// MParameterD with the name "MinimizationResult". -// -// 5) The parameters to start with must have been set using -// MJOptimize::SetParameter or MJOptimize::SetParameters and -// MJOptimize::FixParameter -// -// The behaviour of the optimization can be changed using: -// void SetNumEvents(UInt_t n); -// void SetDebug(UInt_t n); -// void SetOptimizer(Optimizer_t o); -// -// After optimization the resulting parameters are set and another eventloop -// with a MStatusDisplay is set is called. The resulting parameters can be -// accessed using: GetParameters() -// -// To be fixed: -// - MStatusDisplay should show status while optimization is running -// - write result into MStatusDisplay -// - save result -// -Bool_t MJOptimize::Optimize(MParList &parlist) -{ - // Checks to make sure, that fcn doesn't crash - if (!parlist.FindCreateObj("MParameterD", fNameMinimizationValue)) - return kFALSE; - - if (!parlist.FindObject("MParameters", "MParContainer")) - { - *fLog << err << "MParameters [MParContainer] not found... abort." << endl; - return kFALSE; - } - - MMatrixLoop *loop = dynamic_cast(parlist.FindTask("MRead")); - - TString txt("Starting "); - switch (fType) - { - case kMigrad: txt += "Migrad"; break; - case kMinimize: txt += "Minimize"; break; - case kMinos: txt += "Minos"; break; - case kImprove: txt += "Improve"; break; - case kSimplex: txt += "Simplex"; break; - case kSeek: txt += "Seek"; break; - case kNone: txt += "no"; break; - } - txt += " optimization"; - - fLog->Separator(txt); - - // Setup eventloop - MEvtLoop evtloop; - evtloop.SetParList(&parlist); - evtloop.SetDisplay(fDisplay); // set display for evtloop and all childs - parlist.SetDisplay(0); // reset display for all contents of parlist and tasklist - evtloop.SetPrivateDisplay(); // prevent display from being cascaded again in PreProcess - - *fLog << inf << "Number of Parameters: " << fParameters.GetSize() << endl; - - // In case the reader is the matrix loop and testrain is enabled - // switch on even mode... - if (loop && TMath::Abs(fTestTrain)>0) - loop->SetOperationMode(fTestTrain>0?MMatrixLoop::kEven:MMatrixLoop::kOdd); - - if (!Optimize(evtloop)) - return kFALSE; - - gMinuit = 0; - - fEvtLoop->SetDisplay(fDisplay); - if (!Fcn(fParameters)) - return kFALSE; - - // In case the reader is the matrix loop and testrain is enabled - // switch on odd mode... - if (!loop || fTestTrain==0) - return kTRUE; - - loop->SetOperationMode(fTestTrain<0?MMatrixLoop::kEven:MMatrixLoop::kOdd); - - // Done already in Fcn - // list.SetVariables(fParameters); - return Fcn(fParameters); -} - -void MJOptimize::AddRulesToMatrix(MHMatrix &m) const -{ - TIter Next1(&fRules); - TObject *o1=0; - while ((o1=Next1())) - m.AddColumn(o1->GetName()); -} - -//------------------------------------------------------------------------ -// -// Fill matrix m by using read. Use rules as a filter if userules. -// -Bool_t MJOptimize::FillMatrix(MReadTree &read, MParList &parlist, Bool_t userules) -{ - MHMatrix *m = (MHMatrix*)parlist.FindObject("M", "MHMatrix"); - if (!m) - { - *fLog << err << "MJOptimize::FillMatrix - ERROR: M [MHMatrix] not found in parlist... abort." << endl; - return kFALSE; - } - - m->Print("cols"); - - //MParList parlist; - - // MGeomCamMagic cam; - // parlist.AddToList(&cam); - - MTaskList tlist; - parlist.Replace(&tlist); - - MFillH fillh(m); - - tlist.AddToList(&read); - - MFilterList list; - if (!list.AddToList(fPreCuts)) - *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl; - if (userules) - SetupFilters(list); - list.SetName("PreCuts"); // reset Name set by SetupFilters - list.SetInverted(kFALSE); // reset inversion set by SetupFilters - fillh.SetFilter(&list); - tlist.AddToList(&list); - - tlist.AddToList(&fillh); - - MEvtLoop fillloop; - fillloop.SetParList(&parlist); - fillloop.SetDisplay(fDisplay); - if (!fillloop.Eventloop(fNumEvents)) - { - *fLog << err << "Filling matrix failed..." << endl; - return kFALSE; - } - - *fLog << inf << "Read events from file '" << read.GetFileName() << "'" << endl; - - if (fillloop.GetDisplay()!=fDisplay) - { - fDisplay = 0; - *fLog << inf << "Optimization aborted by user." << endl; - return kFALSE; - } - - m->Print("size"); - - return kTRUE; -} - -//------------------------------------------------------------------------ -// -// Adds all filters to MFilterList -// -void MJOptimize::SetupFilters(MFilterList &list, MFilter *filter) const -{ - list.SetName("MParameters"); - list.SetInverted(); - - if (filter) - { - if (fFilter.GetSize()>0) - { - *fLog << inf; - *fLog << "INFORMATION - You are using an external filter and internal filters." << endl; - *fLog << " Please make sure that all parameters '[i]' are starting" << endl; - *fLog << " behind the number of parameters of the external filter." << endl; - } - list.AddToList(filter); - } - - if (!list.AddToList(fFilter)) - *fLog << err << "ERROR - Calling MFilterList::AddToList fFilter failed!" << endl; - - *fLog << inf << "Filter: "; - list.Print(); - *fLog << endl; -} Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.h (revision 7151) +++ (revision ) @@ -1,161 +1,0 @@ -#ifndef MARS_MJOptimize -#define MARS_MJOptimize - -#ifndef MARS_MJob -#include "MJob.h" -#endif - -#ifndef ROOT_TArrayD -#include -#endif - -class TMinuit; - -class MAlphaFitter; - -class MEvtLoop; -class MParList; -class MFilter; -class MFilterList; -class MFitParameters; -class MFitParametersCalc; - -class MHMatrix; -class MGeomCam; -class MRead; -class MReadTree; - -class MJOptimize : public MJob -{ -public: - enum Optimizer_t - { - kMigrad, // Minimize by the method of Migrad - kSimplex, // Minimize by the method of Simplex - kMinimize, // Migrad + Simplex (if Migrad fails) - kMinos, // Minos error determination - kImprove, // Local minimum search - kSeek, // Minimize by the method of Monte Carlo - kNone // Skip optimization - }; - -private: - Int_t fDebug; // -1 no output, 0 MJOptimize output, 1 PrintStatistics output - Int_t fNumEvents; - - TList fRules; - TList fFilter; - - void AddPoint(TList *l, Int_t idx, Float_t val) const; - TList *GetPlots() const; - - MEvtLoop *fEvtLoop; //! - - // Minuit Interface - static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); - Double_t Fcn(const TArrayD &par, TMinuit *minuit=0); - - Int_t Minuit(TMinuit &minuit, const char *cmd) const; - Int_t Migrad(TMinuit &minuit) const { return Minuit(minuit, "MIGRAD"); } - Int_t Simplex(TMinuit &minuit) const { return Minuit(minuit, "SIMPLEX"); } - Int_t Minimize(TMinuit &minuit) const { return Minuit(minuit, "MINIMIZE"); } - Int_t Seek(TMinuit &minuit) const { return Minuit(minuit, "SEEK"); } - Int_t Improve(TMinuit &minuit) const { return Minuit(minuit, "IMPROVE"); } - Int_t Minos(TMinuit &minuit) const { return Minuit(minuit, "MINOS"); } - - TArrayD fParameters; //! - TArrayD fLimLo; //! - TArrayD fLimUp; //! - TArrayD fStep; //! - - Optimizer_t fType; - UInt_t fNumMaxCalls; - Float_t fTolerance; - Int_t fTestTrain; - TString fNameMinimizationValue; - - Bool_t Optimize(MEvtLoop &evtloop); - -protected: - TList fPreCuts; - - TList fNamesOn; - TList fNamesOff; - - TString fNameOut; - - void AddRulesToMatrix(MHMatrix &m) const; - void SetupFilters(MFilterList &list, MFilter *filter=0) const; - Bool_t AddSequences(MRead &read, TList &list) const; - Bool_t FillMatrix(MReadTree &read, MParList &l, Bool_t userules=kFALSE); - -public: - MJOptimize(); - - // I/O - void AddSequenceOn(const char *fname, const char *dir=""); - void AddSequenceOff(const char *fname, const char *dir=""); - - void ResetSequences(); - - // Interface for filter cuts - Int_t AddParameter(const char *rule); - void AddFilter(const char *rule); - void AddPreCut(const char *rule); - void AddPreCut(MFilter *f); - - // Steering of optimization - void SetNumEvents(UInt_t n); - void SetDebug(UInt_t n); - void SetNameOut(const char *name="") { fNameOut = name; } - void SetOptimizer(Optimizer_t o); - void SetNumMaxCalls(UInt_t num=0) { fNumMaxCalls=num; } - void SetTolerance(Float_t tol=0) { fTolerance=tol; } - void EnableTestTrain(Int_t b=1) { fTestTrain=b; } // Use 1 and -1 - void SetNameMinimizationValue(const char *name="MinimizationValue") { fNameMinimizationValue = name; } - - // Parameter access - void SetParameters(const TArrayD &par); - void SetParameter(Int_t idx, Double_t start, Double_t lo=0, Double_t up=0, Double_t step=-1) - { - if (fParameters.GetSize()<=idx) - { - fParameters.Set(idx+1); - fLimLo.Set(idx+1); - fLimUp.Set(idx+1); - fStep.Set(idx+1); - } - - fParameters[idx] = start; - fLimLo[idx] = lo; - fLimUp[idx] = up; - if (step<=0) - fStep[idx] = start==0 ? 0.1 : TMath::Abs(start*0.15); - else - fStep[idx] = step; - } - void FixParameter(Int_t idx, Double_t start, Double_t lo=0, Double_t up=0, Double_t step=-1) - { - if (fParameters.GetSize()<=idx) - { - fParameters.Set(idx+1); - fLimLo.Set(idx+1); - fLimUp.Set(idx+1); - fStep.Set(idx+1); - } - - fParameters[idx] = start; - fLimLo[idx] = 0; - fLimUp[idx] = 0; - fStep[idx] = 0; - } - - const TArrayD &GetParameters() const { return fParameters; } - - // Generalized optimizing routines - Bool_t Optimize(MParList &list); - - ClassDef(MJOptimize, 0) // Class for optimization of the Supercuts -}; - -#endif Index: trunk/MagicSoft/Mars/mjobs/MJOptimizeCuts.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimizeCuts.cc (revision 7151) +++ (revision ) @@ -1,344 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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, 9/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MJOptimizeCuts -// -// Class for otimizing the parameters of the supercuts -// -// Minimization Control -// ==================== -// -// To choose the minimization algorithm use: -// void SetOptimizer(Optimizer_t o); -// -// Allowed options are: -// enum Optimizer_t -// { -// kMigrad, // Minimize by the method of Migrad -// kSimplex, // Minimize by the method of Simplex -// kMinimize, // Migrad + Simplex (if Migrad fails) -// kMinos, // Minos error determination -// kImprove, // Local minimum search -// kSeek, // Minimize by the method of Monte Carlo -// kNone // Skip optimization -// }; -// -// For more details on the methods see TMinuit. -// -// -// You can change the behaviour of the minimization using -// -// void SetNumMaxCalls(UInt_t num=0); -// void SetTolerance(Float_t tol=0); -// -// While NumMaxCalls is the first, Tolerance the second arguement. -// For more details start root and type -// -// gMinuit->mnhelp("command") -// -// while command can be -// * MIGRAD -// * SIMPLEX -// * MINIMIZE -// * MINOS -// * IMPROVE -// * SEEK -// -// The default (num==0 and tol==0) should always give you the -// corresponding defaults used in Minuit. -// -// -// FIXME: Implement changing cut in hadronness... -// FIXME: Show MHSignificance on MStatusDisplay during filling... -// FIXME: Choose step-size percentage as static data membewr -// FIXME: Choose minimization method -// -///////////////////////////////////////////////////////////////////////////// -#include "MJOptimizeCuts.h" - -#include "MHMatrix.h" - -// Environment -#include "MLog.h" -#include "MLogManip.h" - -// Eventloop -#include "MReadTree.h" -#include "MParList.h" -#include "MTaskList.h" - -// Parameter container -#include "MGeomCamMagic.h" -#include "MParameters.h" -#include "MFilterList.h" - -// histograms -#include "../mhflux/MHAlpha.h" - -// Tasks -#include "MF.h" -#include "MFillH.h" -#include "MContinue.h" -#include "MMatrixLoop.h" - -#include "MFMagicCuts.h" - -ClassImp(MJOptimizeCuts); - -using namespace std; - -//------------------------------------------------------------------------ -// -MHAlpha *MJOptimizeCuts::CreateNewHist(const char *name) const -{ - TClass *cls = gROOT->GetClass(fNameHist); - if (!cls) - { - *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl; - return NULL; - } - if (!cls->InheritsFrom(MHAlpha::Class())) - { - *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl; - return NULL; - } - - MHAlpha *h = (MHAlpha*)cls->New(); - if (h && name) - h->SetName(name); - - return h; -} - -//------------------------------------------------------------------------ -// -Bool_t MJOptimizeCuts::RunOnOffCore(MHAlpha &histon, MHAlpha &histof, const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree) -{ - fLog->Separator("Preparing On/Off-optimization"); - - MParList parlist; - - MGeomCamMagic geom; // For GetConvMm2Deg - parlist.AddToList(&geom); - - MHMatrix m("M"); - AddRulesToMatrix(m); - parlist.AddToList(&m); - - const Int_t idxdatatype = m.AddColumn("DataType.fVal"); - - histon.SkipHistTime(); - histon.SkipHistTheta(); - //histon.SkipHistEnergy(); - histof.SkipHistTime(); - histof.SkipHistTheta(); - //histof.SkipHistEnergy(); - histon.ForceUsingSize(); - histof.ForceUsingSize(); - histon.InitMapping(&m, 1); - histof.InitMapping(&m, 1); - - if (filter && filter->InheritsFrom(MFMagicCuts::Class())) - ((MFMagicCuts*)filter)->InitMapping(&m); - - parlist.AddToList(&histon); - parlist.AddToList(&histof); - - if (fname) - { - MReadTree read(tree); - read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! - read.AddFile(fname); - if (!FillMatrix(read, parlist)) - return kFALSE; - } - else - { - MParameterI par("DataType"); - parlist.AddToList(&par); - - gLog.Separator("Reading On-Data"); - par.SetVal(1); - MReadTree readon(tree); - readon.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! - AddSequences(readon, fNamesOn); - if (!FillMatrix(readon, parlist)) - return kFALSE; - - gLog.Separator("Reading Off-Data"); - par.SetVal(0); - MReadTree readoff(tree); - readoff.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! - AddSequences(readoff, fNamesOff); - if (!FillMatrix(readoff, parlist)) - return kFALSE; - } - - MTaskList tasklist; - parlist.Replace(&tasklist); - if (fit) - parlist.AddToList(fit); - - MFilterList list; - SetupFilters(list, filter); - - MContinue contin(&list); - parlist.AddToList(&list); - - MFillH fillof(&histof, "", "FillHistOff"); - MFillH fillon(&histon, "", "FillHistOn"); - - MF f0(Form("M[%d]<0.5", idxdatatype), "FilterOffData"); - MF f1(Form("M[%d]>0.5", idxdatatype), "FilterOnData"); - - fillof.SetFilter(&f0); - fillon.SetFilter(&f1); - - MMatrixLoop loop(&m); - - tasklist.AddToList(&loop); - tasklist.AddToList(&list); - tasklist.AddToList(&contin); - tasklist.AddToList(&f0); - tasklist.AddToList(&f1); - tasklist.AddToList(&fillof); - tasklist.AddToList(&fillon); - - // Optimize with the tasklist in this parameterlist - if (!Optimize(parlist)) - return kFALSE; - - // Copy the result back to be accessible by the user - if (fit) - histon.GetAlphaFitter().Copy(*fit); - - // Print the result - histon.GetAlphaFitter().Print("result"); - - // Store result if requested - TObjArray cont; - cont.Add(&contin); - return WriteContainer(cont, fNameOut); -} - -Bool_t MJOptimizeCuts::RunOnCore(MHAlpha &hist, const char *fname, MFilter *filter, MAlphaFitter *fit) -{ - fLog->Separator("Preparing On-only-optimization"); - - MParList parlist; - - MGeomCamMagic geom; // For GetConvMm2Deg - parlist.AddToList(&geom); - - MHMatrix m("M"); - AddRulesToMatrix(m); - parlist.AddToList(&m); - - hist.SkipHistTime(); - hist.SkipHistTheta(); - hist.SkipHistEnergy(); - hist.InitMapping(&m); - - if (filter && filter->InheritsFrom(MFMagicCuts::Class())) - ((MFMagicCuts*)filter)->InitMapping(&m); - - MReadTree read("Events"); - read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! - if (fname) - read.AddFile(fname); - else - AddSequences(read, fNamesOn); - if (!FillMatrix(read, parlist)) - return kFALSE; - - MTaskList tasklist; - parlist.Replace(&tasklist); - if (fit) - parlist.AddToList(fit); - - MFilterList list; - SetupFilters(list, filter); - - MContinue contin(&list); - parlist.AddToList(&list); - - MFillH fill(&hist); - - MMatrixLoop loop(&m); - - tasklist.AddToList(&loop); - tasklist.AddToList(&list); - tasklist.AddToList(&contin); - tasklist.AddToList(&fill); - - // Optimize with the tasklist in this parameterlist - if (!Optimize(parlist)) - return kFALSE; - - // Copy the result back to be accessible by the user - if (fit) - hist.GetAlphaFitter().Copy(*fit); - // Print the result - hist.GetAlphaFitter().Print("result"); - - // Store result if requested - TObjArray cont; - cont.Add(&contin); - return WriteContainer(cont, fNameOut); -} - -//------------------------------------------------------------------------ -// -Bool_t MJOptimizeCuts::RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree) -{ - MHAlpha *histon = CreateNewHist("Hist"); - MHAlpha *histof = CreateNewHist("HistOff"); - - if (!histon || !histof) - return kFALSE; - - const Bool_t rc = RunOnOffCore(*histon, *histof, fname, filter, fit, tree); - - delete histon; - delete histof; - - return rc; -} - -//------------------------------------------------------------------------ -// -Bool_t MJOptimizeCuts::RunOn(const char *fname, MFilter *filter, MAlphaFitter *fit) -{ - MHAlpha *histon = CreateNewHist(); - - if (!histon) - return kFALSE; - - const Bool_t rc = RunOnCore(*histon, fname, filter, fit); - - delete histon; - return rc; -} Index: trunk/MagicSoft/Mars/mjobs/MJOptimizeCuts.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimizeCuts.h (revision 7151) +++ (revision ) @@ -1,57 +1,0 @@ -#ifndef MARS_MJOptimizeCuts -#define MARS_MJOptimizeCuts - -#ifndef MARS_MJOptimize -#include "MJOptimize.h" -#endif - -class MHAlpha; - -class MJOptimizeCuts : public MJOptimize -{ -private: - TString fNameHist; - - MHAlpha *CreateNewHist(const char *name=0) const; - - Bool_t RunOnCore(MHAlpha &hist, const char *fname, MFilter *filter, MAlphaFitter *fit); - Bool_t RunOnOffCore(MHAlpha &histon, MHAlpha &histoff, const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree); - -public: - MJOptimizeCuts(const char *name="MHAlpha") : MJOptimize(), fNameHist(name) { } - - // Optimization On-Off/Wobble - Bool_t RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit=0, const char *tree="Events"); - Bool_t RunOnOff(const char *fname, MAlphaFitter *fit=0, const char *tree="Events") - { - return RunOnOff(fname, 0, fit, tree); - } - Bool_t RunOnOff(MFilter *filter, MAlphaFitter *fit=0, const char *tree="Events") - { - return RunOnOff(0, filter, fit, tree); - } - Bool_t RunOnOff(MAlphaFitter *fit=0, const char *tree="Events") - { - return RunOnOff(fit, tree); - } - - // Optimization On-only - Bool_t RunOn(const char *fname, MFilter *filter, MAlphaFitter *fit=0); - Bool_t RunOn(const char *fname, MAlphaFitter *fit=0) - { - return RunOn(fname, 0, fit); - } - - Bool_t RunOn(MFilter *filter, MAlphaFitter *fit=0) - { - return RunOn(0, filter, fit); - } - Bool_t RunOn(MAlphaFitter *fit=0) - { - return RunOn(0, 0, fit); - } - - ClassDef(MJOptimizeCuts, 0) // Class for optimization of the Supercuts -}; - -#endif Index: trunk/MagicSoft/Mars/mjobs/MJOptimizeEnergy.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimizeEnergy.cc (revision 7151) +++ (revision ) @@ -1,151 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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, 9/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MJOptimize -// -// Class for otimizing a rule to estimate the energy. For more details see -// MJOptimize. -// -// Example: -// -------- -// -// MJOptimizeEnergy opt; -// opt.SetDebug(2); -// opt.SetOptimizer(MJOptimize::kMigrad); -// opt.SetNumEvents(20000); -// opt.EnableTestTrain(); -// opt.AddParameter("MHillas.fSize"); -// opt.SetParameter(0, 1, 0, 2); -// char *r = "[0]*M[0]"; //Rule to calculate estimated energy -// MStatusDisplay *d = new MStatusDisplay; -// opt.SetDisplay(d); -// opt.RunDisp("ganymed-summary.root", r); -// -///////////////////////////////////////////////////////////////////////////// -#include "MJOptimizeEnergy.h" - -#include "MHMatrix.h" - -// environment -#include "MLog.h" -#include "MLogManip.h" - -// eventloop -#include "MParList.h" -#include "MTaskList.h" -#include "MEvtLoop.h" - -// parameters -#include "MParameters.h" -#include "MGeomCamMagic.h" - -// histograms -#include "../mhflux/MHEnergyEst.h" - -// tasks -#include "MReadTree.h" -#include "MMatrixLoop.h" -#include "MEnergyEstimate.h" -#include "MFillH.h" - -// filters -#include "MFDataMember.h" - -using namespace std; - -ClassImp(MJOptimizeEnergy); - -//------------------------------------------------------------------------ -// -// Read all events from file which do match rules and optimize -// energy estimator. -// -Bool_t MJOptimizeEnergy::RunEnergy(const char *fname, const char *rule, MTask *weights) -{ - fLog->Separator("Preparing Energy optimization"); - - MParList parlist; - - MParameterI par("DataType"); - par.SetVal(1); - parlist.AddToList(&par); - - MFDataMember filter("DataType.fVal", '>', 0.5); - fPreCuts.Add(&filter); - - MGeomCamMagic geom; // For GetConvMm2Deg - parlist.AddToList(&geom); - - MHMatrix m("M"); - AddRulesToMatrix(m); - parlist.AddToList(&m); - - MHEnergyEst hist; - hist.InitMapping(&m); - - MParameterCalc est(rule, "MParameters"); - est.SetNameParameter("MEnergyEst"); - parlist.AddToList(&est); - - MReadTree read("Events"); - // NECESSARY BECAUSE OF MDataFormula GetRules missing - read.DisableAutoScheme(); - if (fname) - read.AddFile(fname); - else - AddSequences(read, fNamesOn); - if (!FillMatrix(read, parlist, kTRUE)) - return kFALSE; - - fPreCuts.Remove(&filter); - - MTaskList tasklist; - parlist.Replace(&tasklist); - - MFillH fill(&hist); - if (weights) - fill.SetWeight(); - - MMatrixLoop loop(&m); - - tasklist.AddToList(&loop); - tasklist.AddToList(&est); - if (weights) - tasklist.AddToList(weights); - tasklist.AddToList(&fill); - - // Optimize with the tasklist in this parameterlist - if (!Optimize(parlist)) - return kFALSE; - - // Print the result - hist.Print(); - - // Store result if requested - TObjArray cont; - cont.Add(&est); - return WriteContainer(cont, fNameOut); -} Index: trunk/MagicSoft/Mars/mjobs/MJOptimizeEnergy.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimizeEnergy.h (revision 7151) +++ (revision ) @@ -1,25 +1,0 @@ -#ifndef MARS_MJOptimizeEnergy -#define MARS_MJOptimizeEnergy - -#ifndef MARS_MJOptimize -#include "MJOptimize.h" -#endif - -class MTask; - -class MJOptimizeEnergy : public MJOptimize -{ -public: - MJOptimizeEnergy() : MJOptimize() { } - - // Special optimizing routines - Bool_t RunEnergy(const char *fname, const char *rule, MTask *weights=0); - Bool_t RunEnergy(const char *rule, MTask *weights=0) - { - return RunEnergy(0, rule, weights); - } - - ClassDef(MJOptimizeEnergy, 0) // Class for optimization of the Supercuts -}; - -#endif Index: trunk/MagicSoft/Mars/mjobs/Makefile =================================================================== --- trunk/MagicSoft/Mars/mjobs/Makefile (revision 7151) +++ trunk/MagicSoft/Mars/mjobs/Makefile (revision 7152) @@ -35,7 +35,4 @@ MJStar.cc \ MJCut.cc \ - MJOptimize.cc \ - MJOptimizeCuts.cc \ - MJOptimizeEnergy.cc \ MJSpectrum.cc Index: trunk/MagicSoft/Mars/mjoptim/MJOptimize.cc =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimize.cc (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimize.cc (revision 7152) @@ -0,0 +1,808 @@ +/* ======================================================================== *\ +! +! * +! * 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, 9/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJOptimize +// +// Class for otimizing the parameters of the supercuts +// +// Minimization Control +// ==================== +// +// To choose the minimization algorithm use: +// void SetOptimizer(Optimizer_t o); +// +// Allowed options are: +// enum Optimizer_t +// { +// kMigrad, // Minimize by the method of Migrad +// kSimplex, // Minimize by the method of Simplex +// kMinimize, // Migrad + Simplex (if Migrad fails) +// kMinos, // Minos error determination +// kImprove, // Local minimum search +// kSeek, // Minimize by the method of Monte Carlo +// kNone // Skip optimization +// }; +// +// For more details on the methods see TMinuit. +// +// +// You can change the behaviour of the minimization using +// +// void SetNumMaxCalls(UInt_t num=0); +// void SetTolerance(Float_t tol=0); +// +// While NumMaxCalls is the first, Tolerance the second arguement. +// For more details start root and type +// +// gMinuit->mnhelp("command") +// +// while command can be +// * MIGRAD +// * SIMPLEX +// * MINIMIZE +// * MINOS +// * IMPROVE +// * SEEK +// +// The default (num==0 and tol==0) should always give you the +// corresponding defaults used in Minuit. +// +// +// FIXME: Implement changing cut in hadronness... +// FIXME: Show MHSignificance on MStatusDisplay during filling... +// FIXME: Choose step-size percentage as static data membewr +// FIXME: Choose minimization method +// +///////////////////////////////////////////////////////////////////////////// +#include "MJOptimize.h" + +#include +#include + +#include + +#include + +#include +#include + +#include "MHMatrix.h" + +// environment +#include "MLog.h" +#include "MLogManip.h" + +#include "MDirIter.h" +#include "MStatusDisplay.h" + +// eventloop +#include "MParList.h" +#include "MTaskList.h" +#include "MEvtLoop.h" + +// parameters +#include "MParameters.h" + +// tasks +#include "MReadTree.h" +#include "MMatrixLoop.h" +#include "MFillH.h" + +// filters +#include "MF.h" +#include "MFilterList.h" + +using namespace std; + +//------------------------------------------------------------------------ +// +// fcn calculates the function to be minimized (using TMinuit::Migrad) +// +void MJOptimize::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + MJOptimize *optim = (MJOptimize*)gMinuit->GetObjectFit(); + + // WORKAROUND --- FOR WHAT? + if (gMinuit->fEpsi<1e-2) + { + *optim->fLog << warn << "WARNING - For unknown reasons: fEspi<1e-100... resetting to 0.01." << endl; + gMinuit->fEpsi = 0.01; + } + + TMinuit *minuit = gMinuit; + f = optim->Fcn(TArrayD(TMath::Min(gMinuit->fMaxpar, optim->fParameters.GetSize()), par), minuit); + gMinuit = minuit; + +} + +Double_t MJOptimize::Fcn(const TArrayD &par, TMinuit *minuit) +{ + if (fEvtLoop->GetDisplay()!=fDisplay) + return 0; + /* + switch(iflag) + { + case 1: // first call + case 2: // calc derivative + break; + case 3: + // last call + MStatusDisplay *d = new MStatusDisplay; + fEvtLoop->SetDisplay(d); + break; + } + */ + MParList *plist = fEvtLoop->GetParList(); + + MParameterD *eval = (MParameterD*)plist->FindObject(fNameMinimizationValue, "MParameterD"); + MParContainer *pars = (MParContainer*)plist->FindObject("MParameters", "MParContainer"); + + MRead *read = (MRead*)plist->FindObject("MTaskList")->FindObject("MRead"); + if (read) + read->Rewind(); + + if (fDebug>=0) + { + *fLog << inf << "New Set: "; + for (Int_t i=0; iSetVariables(par); + eval->SetVal(0); + + const Bool_t isnull = gLog.IsNullOutput(); + if (fDebug<3) + gLog.SetNullOutput(kTRUE); + + TStopwatch clock; + clock.Start(); + fEvtLoop->Eventloop(fNumEvents, MEvtLoop::kNoStatistics); + clock.Stop(); + + if (fDebug<3) + gLog.SetNullOutput(isnull); + + const Double_t f = eval->GetVal(); + + if (fDebug>=0) + *fLog << inf << "Result F=" << f << endl; + + if (fDebug>=1 && minuit) + { + Double_t fmin, fedm, errdef; + Int_t n1, n2, istat; + minuit->mnstat(fmin, fedm, errdef, n1, n2, istat); + *fLog << inf << underline << "Minimization Status so far:" << endl; + *fLog << " Calls: " << minuit->fNfcn << " (max=" << gMinuit->fMaxIterations << ")" << endl; + *fLog << " Parameters: fixed=" << gMinuit->fNpfix << ", free=" << gMinuit->fNpar << endl; + *fLog << " Func min: " << fmin << " (Epsi=" << gMinuit->fEpsi << ", Apsi=" << gMinuit->fApsi << ")" << 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>=1) + { + clock.Print(); + fEvtLoop->GetTaskList()->PrintStatistics(); + } + + return f; +} + +MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0), fTestTrain(0), fNameMinimizationValue("MinimizationValue") +{ + fRules.SetOwner(); + fFilter.SetOwner(); + + fNamesOn.SetOwner(); + fNamesOff.SetOwner(); +} + +//------------------------------------------------------------------------ +// +// Add seqeunces from list to reader +// +Bool_t MJOptimize::AddSequences(MRead &read, TList &list) const +{ + MDirIter files; + + TIter Next(&list); + TObject *o=0; + while ((o=Next())) + { + MSequence seq(o->GetName()); + if (!seq.IsValid()) + return kFALSE; + + seq.SetupDatRuns(files, o->GetTitle(), "I"); + } + + return read.AddFiles(files)>0; +} + +//------------------------------------------------------------------------ +// +// Add on-sequences: +// - fname: sequence file name (with path) +// - dir: directory were image files are stored +// +void MJOptimize::AddSequenceOn(const char *fname, const char *dir) +{ + fNamesOn.Add(new TNamed(fname, dir)); +} + +//------------------------------------------------------------------------ +// +// Add off-sequences: +// - fname: sequence file name (with path) +// - dir: directory were image files are stored +// +void MJOptimize::AddSequenceOff(const char *fname, const char *dir) +{ + fNamesOff.Add(new TNamed(fname, dir)); +} + +//------------------------------------------------------------------------ +// +// Empty list of on- and off-sequences +// +void MJOptimize::ResetSequences() +{ + fNamesOn.Delete(); + fNamesOff.Delete(); +} + +//------------------------------------------------------------------------ +// +// Add a parameter used in your filters (see AddFilter) The parameter +// index is returned, +// +// Int_t idx = AddParameter("log10(MHillas.fSize)"); +// +// The indices area starting with 0 always. +// +Int_t MJOptimize::AddParameter(const char *rule) +{ + fRules.Add(new TNamed(rule, "")); + return fRules.GetSize()-1; +} + +//------------------------------------------------------------------------ +// +// Add a filter which can be applied in the optimization (for deatils +// see correspodning Run function) You can use the indices you got by +// AddParameter, eg +// +// AddFilter("M[0] < 3.2"); +// +// if used in optimization you can do +// +// AddFilter("M[0] < [0]"); +// +// for more details, see SetParameter and FixParameter +// +void MJOptimize::AddFilter(const char *rule) +{ + fFilter.Add(new MF(rule)); +} + +//------------------------------------------------------------------------ +// +// Add a cut which is used to fill the matrix, eg "MMcEvt.fOartId<1.5" +// (The rule is applied, nit inverted: The matrix is filled with +// the events fullfilling the condition) +// +void MJOptimize::AddPreCut(const char *rule) +{ + MFilter *f = new MF(rule); + f->SetBit(kCanDelete); + AddPreCut(f); +} + +//------------------------------------------------------------------------ +// +// Add a cut which is used to fill the matrix. If kCanDelete is set +// MJOptimize takes the ownership. +// +void MJOptimize::AddPreCut(MFilter *f) +{ + fPreCuts.Add(f); +} + +//------------------------------------------------------------------------ +// +// Set the fParameters Array accoring to par. +// +void MJOptimize::SetParameters(const TArrayD &par) +{ + fParameters = par; +} + +//------------------------------------------------------------------------ +// +// Set the number of events processed by the eventloop. (Be carfull, +// if you are doing on-off analysis and you only process the first +// 1000 events which are on-events only the optimization may not work) +// +void MJOptimize::SetNumEvents(UInt_t n) +{ + fNumEvents=n; +} + +//------------------------------------------------------------------------ +// +// Set a debug level, which tells the optimization how much information +// is displayed about and in the running eventloop. +// +void MJOptimize::SetDebug(UInt_t n) +{ + fDebug=n; +} + +//------------------------------------------------------------------------ +// +// Set a optimization algorithm to be used. For more information see +// TMinuit. +// +// Available Algorithms are: +// kMigrad, // Minimize by the method of Migrad +// kSimplex, // Minimize by the method of Simplex +// kSeek // Minimize by the method of Monte Carlo +// +void MJOptimize::SetOptimizer(Optimizer_t o) +{ + fType = o; +} + +//------------------------------------------------------------------------ +// +// If a status didplay is set, search for tab "Optimizer". +// If not found, create it. +// In the tab search for TMultiGraph "Parameters". +// If not found create it. +// If empty create TGraphs. +// Check number of graphs vs. number of parameters. +// return TList with graphs. +// +TList *MJOptimize::GetPlots() const +{ + if (!fDisplay) + return NULL; + + TCanvas *c=fDisplay->GetCanvas("Optimizer"); + if (!c) + c = &fDisplay->AddTab("Optimizer"); + + TMultiGraph *mg = dynamic_cast(c->FindObject("Parameters")); + if (!mg) + mg = new TMultiGraph("Parameters", "Parameters of optimization"); + + TList *l = mg->GetListOfGraphs(); + if (!l) + { + const Int_t n = fParameters.GetSize(); + for (int i=0; iSetLineColor(kBlue); + mg->Add(g, ""); + AddPoint(mg->GetListOfGraphs(), i, i==n?1:fParameters[i]); + } + mg->SetBit(kCanDelete); + mg->Draw("al*"); + + l = mg->GetListOfGraphs(); + } + + return l->GetSize() == fParameters.GetSize()+1 ? l : NULL; +} + +//------------------------------------------------------------------------ +// +// Add a point with y=val as last point in idx-th Tgraph of list l. +// +void MJOptimize::AddPoint(TList *l, Int_t idx, Float_t val) const +{ + if (!l) + return; + + TGraph *gr = (TGraph*)l->At(idx); + gr->SetPoint(gr->GetN(), gr->GetN(), val); +} + +Int_t MJOptimize::Minuit(TMinuit &minuit, const char *cmd) const +{ + Int_t err; + Double_t tmp[2] = { fNumMaxCalls, fTolerance }; + minuit.mnexcm(cmd, tmp, 2, err); + + switch (err) + { + case 0: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm excuted normally." << endl; + break; + case 1: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command is blank... ignored." << endl; + break; + case 2: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command-line syntax error... ignored." << endl; + break; + case 3: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm unknown command... ignored." << endl; + break; + case 4: + *fLog << warn << GetDescriptor() << " TMinuit::mnexcm - Abnormal termination (eg Migrad not converged)" << endl; + break; + /* + case 5: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - Parameters requested." << endl; + break; + case 6: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET INPUT returned." << endl; + break; + case 7: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET TITLE returned." << endl; + break; + case 8: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET COVAR returned." << endl; + break; + case 9: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - reserved." << endl; + break; + case 10: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - END returned." << endl; + break; + case 11: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - EXIT or STOP returned." << endl; + break; + case 12: + *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - RETURN returned." << endl; + break;*/ + } + + return err; +} + +Bool_t MJOptimize::Optimize(MEvtLoop &evtloop) +{ + if (fParameters.GetSize()==0) + { + *fLog << err << GetDescriptor() << "::Optimize: ERROR - Sorry, no parameters defined." << endl; + return kFALSE; + } + + if (fType==kNone) + return kTRUE; + + gMinuit = new TMinuit(fParameters.GetSize()); + + gMinuit->SetFCN(fcn); + gMinuit->SetObjectFit(this); + gMinuit->SetPrintLevel(-1); // Don't print when DefineParameter + + // + // Set starting values and step sizes for parameters + // + for (Int_t i=0; iDefineParameter(i, name, vinit, step, limlo, limup); + if (rc) + { + *fLog << err << dbginf << "Error in defining parameter #" << i << endl; + return kFALSE; + } + + if (step==0) + gMinuit->FixParameter(i); + } + + gMinuit->SetPrintLevel(1); // Switch on pritning again + gMinuit->mnprin(1,0); // Print all parameters + + fEvtLoop = &evtloop; + + TList *g=GetPlots(); + + // Now ready for minimization step: + TStopwatch clock; + clock.Start(); + switch (fType) + { + case kSimplex: + Simplex(*gMinuit); + break; + case kMigrad: + Migrad(*gMinuit); + break; + case kMinimize: + Minimize(*gMinuit); + break; + case kMinos: + Minos(*gMinuit); + break; + case kImprove: + Improve(*gMinuit); + break; + case kSeek: + Seek(*gMinuit); + break; + case kNone: // Should never happen + return kFALSE; + } + clock.Stop(); + clock.Print(); + + if (evtloop.GetDisplay()!=fDisplay) + { + *fLog << inf << "Optimization aborted by user." << endl; + fDisplay = 0; + return kFALSE; + } + + *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl; + + // + // Update values of fA, fB: + // + for (Int_t i=0; iGetParameter(i,x1,x2); + fParameters[i] = x1; + cout << i << ": " << fParameters[i] << endl; + + AddPoint(g, i, x1); + } + AddPoint(g, fParameters.GetSize(), gMinuit->fAmin); + + delete gMinuit; + + return kTRUE; +} + +//------------------------------------------------------------------------ +// +// Optimize allows to use the optimizing by an eventloop based on +// some requirements. +// +// 1) The tasklist to be executed must have the name MTaskList and +// be an entry in the parameterlist. +// +// 2) The reading task (MReadMarsFile, MMatrixLoop) must have the name +// "MRead". If it derives from MRead Rewind() must be implemented, +// otherwise it must start reading from scratch each time its +// PreProcess is called. +// +// 3) The parameters to be optimized must be accessed through (currently) +// a single parameter container (MParContainer) called MParameters. +// The parameters are set through SetVariables. +// +// 4) The result of a single function call for minimization (eg. chisquare) +// must be available after the eventloop in a container of type +// MParameterD with the name "MinimizationResult". +// +// 5) The parameters to start with must have been set using +// MJOptimize::SetParameter or MJOptimize::SetParameters and +// MJOptimize::FixParameter +// +// The behaviour of the optimization can be changed using: +// void SetNumEvents(UInt_t n); +// void SetDebug(UInt_t n); +// void SetOptimizer(Optimizer_t o); +// +// After optimization the resulting parameters are set and another eventloop +// with a MStatusDisplay is set is called. The resulting parameters can be +// accessed using: GetParameters() +// +// To be fixed: +// - MStatusDisplay should show status while optimization is running +// - write result into MStatusDisplay +// - save result +// +Bool_t MJOptimize::Optimize(MParList &parlist) +{ + // Checks to make sure, that fcn doesn't crash + if (!parlist.FindCreateObj("MParameterD", fNameMinimizationValue)) + return kFALSE; + + if (!parlist.FindObject("MParameters", "MParContainer")) + { + *fLog << err << "MParameters [MParContainer] not found... abort." << endl; + return kFALSE; + } + + MMatrixLoop *loop = dynamic_cast(parlist.FindTask("MRead")); + + TString txt("Starting "); + switch (fType) + { + case kMigrad: txt += "Migrad"; break; + case kMinimize: txt += "Minimize"; break; + case kMinos: txt += "Minos"; break; + case kImprove: txt += "Improve"; break; + case kSimplex: txt += "Simplex"; break; + case kSeek: txt += "Seek"; break; + case kNone: txt += "no"; break; + } + txt += " optimization"; + + fLog->Separator(txt); + + // Setup eventloop + MEvtLoop evtloop; + evtloop.SetParList(&parlist); + evtloop.SetDisplay(fDisplay); // set display for evtloop and all childs + parlist.SetDisplay(0); // reset display for all contents of parlist and tasklist + evtloop.SetPrivateDisplay(); // prevent display from being cascaded again in PreProcess + + *fLog << inf << "Number of Parameters: " << fParameters.GetSize() << endl; + + // In case the reader is the matrix loop and testrain is enabled + // switch on even mode... + if (loop && TMath::Abs(fTestTrain)>0) + loop->SetOperationMode(fTestTrain>0?MMatrixLoop::kEven:MMatrixLoop::kOdd); + + if (!Optimize(evtloop)) + return kFALSE; + + gMinuit = 0; + + fEvtLoop->SetDisplay(fDisplay); + if (!Fcn(fParameters)) + return kFALSE; + + // In case the reader is the matrix loop and testrain is enabled + // switch on odd mode... + if (!loop || fTestTrain==0) + return kTRUE; + + loop->SetOperationMode(fTestTrain<0?MMatrixLoop::kEven:MMatrixLoop::kOdd); + + // Done already in Fcn + // list.SetVariables(fParameters); + return Fcn(fParameters); +} + +void MJOptimize::AddRulesToMatrix(MHMatrix &m) const +{ + TIter Next1(&fRules); + TObject *o1=0; + while ((o1=Next1())) + m.AddColumn(o1->GetName()); +} + +//------------------------------------------------------------------------ +// +// Fill matrix m by using read. Use rules as a filter if userules. +// +Bool_t MJOptimize::FillMatrix(MReadTree &read, MParList &parlist, Bool_t userules) +{ + MHMatrix *m = (MHMatrix*)parlist.FindObject("M", "MHMatrix"); + if (!m) + { + *fLog << err << "MJOptimize::FillMatrix - ERROR: M [MHMatrix] not found in parlist... abort." << endl; + return kFALSE; + } + + m->Print("cols"); + + //MParList parlist; + + // MGeomCamMagic cam; + // parlist.AddToList(&cam); + + MTaskList tlist; + parlist.Replace(&tlist); + + MFillH fillh(m); + + tlist.AddToList(&read); + + MFilterList list; + if (!list.AddToList(fPreCuts)) + *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl; + if (userules) + SetupFilters(list); + list.SetName("PreCuts"); // reset Name set by SetupFilters + list.SetInverted(kFALSE); // reset inversion set by SetupFilters + fillh.SetFilter(&list); + tlist.AddToList(&list); + + tlist.AddToList(&fillh); + + MEvtLoop fillloop; + fillloop.SetParList(&parlist); + fillloop.SetDisplay(fDisplay); + if (!fillloop.Eventloop(fNumEvents)) + { + *fLog << err << "Filling matrix failed..." << endl; + return kFALSE; + } + + *fLog << inf << "Read events from file '" << read.GetFileName() << "'" << endl; + + if (fillloop.GetDisplay()!=fDisplay) + { + fDisplay = 0; + *fLog << inf << "Optimization aborted by user." << endl; + return kFALSE; + } + + m->Print("size"); + + return kTRUE; +} + +//------------------------------------------------------------------------ +// +// Adds all filters to MFilterList +// +void MJOptimize::SetupFilters(MFilterList &list, MFilter *filter) const +{ + list.SetName("MParameters"); + list.SetInverted(); + + if (filter) + { + if (fFilter.GetSize()>0) + { + *fLog << inf; + *fLog << "INFORMATION - You are using an external filter and internal filters." << endl; + *fLog << " Please make sure that all parameters '[i]' are starting" << endl; + *fLog << " behind the number of parameters of the external filter." << endl; + } + list.AddToList(filter); + } + + if (!list.AddToList(fFilter)) + *fLog << err << "ERROR - Calling MFilterList::AddToList fFilter failed!" << endl; + + *fLog << inf << "Filter: "; + list.Print(); + *fLog << endl; +} Index: trunk/MagicSoft/Mars/mjoptim/MJOptimize.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimize.h (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimize.h (revision 7152) @@ -0,0 +1,161 @@ +#ifndef MARS_MJOptimize +#define MARS_MJOptimize + +#ifndef MARS_MJob +#include "MJob.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class TMinuit; + +class MAlphaFitter; + +class MEvtLoop; +class MParList; +class MFilter; +class MFilterList; +class MFitParameters; +class MFitParametersCalc; + +class MHMatrix; +class MGeomCam; +class MRead; +class MReadTree; + +class MJOptimize : public MJob +{ +public: + enum Optimizer_t + { + kMigrad, // Minimize by the method of Migrad + kSimplex, // Minimize by the method of Simplex + kMinimize, // Migrad + Simplex (if Migrad fails) + kMinos, // Minos error determination + kImprove, // Local minimum search + kSeek, // Minimize by the method of Monte Carlo + kNone // Skip optimization + }; + +private: + Int_t fDebug; // -1 no output, 0 MJOptimize output, 1 PrintStatistics output + Int_t fNumEvents; + + TList fRules; + TList fFilter; + + void AddPoint(TList *l, Int_t idx, Float_t val) const; + TList *GetPlots() const; + + MEvtLoop *fEvtLoop; //! + + // Minuit Interface + static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag); + Double_t Fcn(const TArrayD &par, TMinuit *minuit=0); + + Int_t Minuit(TMinuit &minuit, const char *cmd) const; + Int_t Migrad(TMinuit &minuit) const { return Minuit(minuit, "MIGRAD"); } + Int_t Simplex(TMinuit &minuit) const { return Minuit(minuit, "SIMPLEX"); } + Int_t Minimize(TMinuit &minuit) const { return Minuit(minuit, "MINIMIZE"); } + Int_t Seek(TMinuit &minuit) const { return Minuit(minuit, "SEEK"); } + Int_t Improve(TMinuit &minuit) const { return Minuit(minuit, "IMPROVE"); } + Int_t Minos(TMinuit &minuit) const { return Minuit(minuit, "MINOS"); } + + TArrayD fParameters; //! + TArrayD fLimLo; //! + TArrayD fLimUp; //! + TArrayD fStep; //! + + Optimizer_t fType; + UInt_t fNumMaxCalls; + Float_t fTolerance; + Int_t fTestTrain; + TString fNameMinimizationValue; + + Bool_t Optimize(MEvtLoop &evtloop); + +protected: + TList fPreCuts; + + TList fNamesOn; + TList fNamesOff; + + TString fNameOut; + + void AddRulesToMatrix(MHMatrix &m) const; + void SetupFilters(MFilterList &list, MFilter *filter=0) const; + Bool_t AddSequences(MRead &read, TList &list) const; + Bool_t FillMatrix(MReadTree &read, MParList &l, Bool_t userules=kFALSE); + +public: + MJOptimize(); + + // I/O + void AddSequenceOn(const char *fname, const char *dir=""); + void AddSequenceOff(const char *fname, const char *dir=""); + + void ResetSequences(); + + // Interface for filter cuts + Int_t AddParameter(const char *rule); + void AddFilter(const char *rule); + void AddPreCut(const char *rule); + void AddPreCut(MFilter *f); + + // Steering of optimization + void SetNumEvents(UInt_t n); + void SetDebug(UInt_t n); + void SetNameOut(const char *name="") { fNameOut = name; } + void SetOptimizer(Optimizer_t o); + void SetNumMaxCalls(UInt_t num=0) { fNumMaxCalls=num; } + void SetTolerance(Float_t tol=0) { fTolerance=tol; } + void EnableTestTrain(Int_t b=1) { fTestTrain=b; } // Use 1 and -1 + void SetNameMinimizationValue(const char *name="MinimizationValue") { fNameMinimizationValue = name; } + + // Parameter access + void SetParameters(const TArrayD &par); + void SetParameter(Int_t idx, Double_t start, Double_t lo=0, Double_t up=0, Double_t step=-1) + { + if (fParameters.GetSize()<=idx) + { + fParameters.Set(idx+1); + fLimLo.Set(idx+1); + fLimUp.Set(idx+1); + fStep.Set(idx+1); + } + + fParameters[idx] = start; + fLimLo[idx] = lo; + fLimUp[idx] = up; + if (step<=0) + fStep[idx] = start==0 ? 0.1 : TMath::Abs(start*0.15); + else + fStep[idx] = step; + } + void FixParameter(Int_t idx, Double_t start, Double_t lo=0, Double_t up=0, Double_t step=-1) + { + if (fParameters.GetSize()<=idx) + { + fParameters.Set(idx+1); + fLimLo.Set(idx+1); + fLimUp.Set(idx+1); + fStep.Set(idx+1); + } + + fParameters[idx] = start; + fLimLo[idx] = 0; + fLimUp[idx] = 0; + fStep[idx] = 0; + } + + const TArrayD &GetParameters() const { return fParameters; } + + // Generalized optimizing routines + Bool_t Optimize(MParList &list); + + ClassDef(MJOptimize, 0) // Class for optimization of the Supercuts +}; + +#endif Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeCuts.cc =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeCuts.cc (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeCuts.cc (revision 7152) @@ -0,0 +1,344 @@ +/* ======================================================================== *\ +! +! * +! * 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, 9/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJOptimizeCuts +// +// Class for otimizing the parameters of the supercuts +// +// Minimization Control +// ==================== +// +// To choose the minimization algorithm use: +// void SetOptimizer(Optimizer_t o); +// +// Allowed options are: +// enum Optimizer_t +// { +// kMigrad, // Minimize by the method of Migrad +// kSimplex, // Minimize by the method of Simplex +// kMinimize, // Migrad + Simplex (if Migrad fails) +// kMinos, // Minos error determination +// kImprove, // Local minimum search +// kSeek, // Minimize by the method of Monte Carlo +// kNone // Skip optimization +// }; +// +// For more details on the methods see TMinuit. +// +// +// You can change the behaviour of the minimization using +// +// void SetNumMaxCalls(UInt_t num=0); +// void SetTolerance(Float_t tol=0); +// +// While NumMaxCalls is the first, Tolerance the second arguement. +// For more details start root and type +// +// gMinuit->mnhelp("command") +// +// while command can be +// * MIGRAD +// * SIMPLEX +// * MINIMIZE +// * MINOS +// * IMPROVE +// * SEEK +// +// The default (num==0 and tol==0) should always give you the +// corresponding defaults used in Minuit. +// +// +// FIXME: Implement changing cut in hadronness... +// FIXME: Show MHSignificance on MStatusDisplay during filling... +// FIXME: Choose step-size percentage as static data membewr +// FIXME: Choose minimization method +// +///////////////////////////////////////////////////////////////////////////// +#include "MJOptimizeCuts.h" + +#include "MHMatrix.h" + +// Environment +#include "MLog.h" +#include "MLogManip.h" + +// Eventloop +#include "MReadTree.h" +#include "MParList.h" +#include "MTaskList.h" + +// Parameter container +#include "MGeomCamMagic.h" +#include "MParameters.h" +#include "MFilterList.h" + +// histograms +#include "../mhflux/MHAlpha.h" + +// Tasks +#include "MF.h" +#include "MFillH.h" +#include "MContinue.h" +#include "MMatrixLoop.h" + +#include "MFMagicCuts.h" + +ClassImp(MJOptimizeCuts); + +using namespace std; + +//------------------------------------------------------------------------ +// +MHAlpha *MJOptimizeCuts::CreateNewHist(const char *name) const +{ + TClass *cls = gROOT->GetClass(fNameHist); + if (!cls) + { + *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl; + return NULL; + } + if (!cls->InheritsFrom(MHAlpha::Class())) + { + *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl; + return NULL; + } + + MHAlpha *h = (MHAlpha*)cls->New(); + if (h && name) + h->SetName(name); + + return h; +} + +//------------------------------------------------------------------------ +// +Bool_t MJOptimizeCuts::RunOnOffCore(MHAlpha &histon, MHAlpha &histof, const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree) +{ + fLog->Separator("Preparing On/Off-optimization"); + + MParList parlist; + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + const Int_t idxdatatype = m.AddColumn("DataType.fVal"); + + histon.SkipHistTime(); + histon.SkipHistTheta(); + //histon.SkipHistEnergy(); + histof.SkipHistTime(); + histof.SkipHistTheta(); + //histof.SkipHistEnergy(); + histon.ForceUsingSize(); + histof.ForceUsingSize(); + histon.InitMapping(&m, 1); + histof.InitMapping(&m, 1); + + if (filter && filter->InheritsFrom(MFMagicCuts::Class())) + ((MFMagicCuts*)filter)->InitMapping(&m); + + parlist.AddToList(&histon); + parlist.AddToList(&histof); + + if (fname) + { + MReadTree read(tree); + read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + read.AddFile(fname); + if (!FillMatrix(read, parlist)) + return kFALSE; + } + else + { + MParameterI par("DataType"); + parlist.AddToList(&par); + + gLog.Separator("Reading On-Data"); + par.SetVal(1); + MReadTree readon(tree); + readon.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + AddSequences(readon, fNamesOn); + if (!FillMatrix(readon, parlist)) + return kFALSE; + + gLog.Separator("Reading Off-Data"); + par.SetVal(0); + MReadTree readoff(tree); + readoff.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + AddSequences(readoff, fNamesOff); + if (!FillMatrix(readoff, parlist)) + return kFALSE; + } + + MTaskList tasklist; + parlist.Replace(&tasklist); + if (fit) + parlist.AddToList(fit); + + MFilterList list; + SetupFilters(list, filter); + + MContinue contin(&list); + parlist.AddToList(&list); + + MFillH fillof(&histof, "", "FillHistOff"); + MFillH fillon(&histon, "", "FillHistOn"); + + MF f0(Form("M[%d]<0.5", idxdatatype), "FilterOffData"); + MF f1(Form("M[%d]>0.5", idxdatatype), "FilterOnData"); + + fillof.SetFilter(&f0); + fillon.SetFilter(&f1); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&list); + tasklist.AddToList(&contin); + tasklist.AddToList(&f0); + tasklist.AddToList(&f1); + tasklist.AddToList(&fillof); + tasklist.AddToList(&fillon); + + // Optimize with the tasklist in this parameterlist + if (!Optimize(parlist)) + return kFALSE; + + // Copy the result back to be accessible by the user + if (fit) + histon.GetAlphaFitter().Copy(*fit); + + // Print the result + histon.GetAlphaFitter().Print("result"); + + // Store result if requested + TObjArray cont; + cont.Add(&contin); + return WriteContainer(cont, fNameOut); +} + +Bool_t MJOptimizeCuts::RunOnCore(MHAlpha &hist, const char *fname, MFilter *filter, MAlphaFitter *fit) +{ + fLog->Separator("Preparing On-only-optimization"); + + MParList parlist; + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + hist.SkipHistTime(); + hist.SkipHistTheta(); + hist.SkipHistEnergy(); + hist.InitMapping(&m); + + if (filter && filter->InheritsFrom(MFMagicCuts::Class())) + ((MFMagicCuts*)filter)->InitMapping(&m); + + MReadTree read("Events"); + read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster! + if (fname) + read.AddFile(fname); + else + AddSequences(read, fNamesOn); + if (!FillMatrix(read, parlist)) + return kFALSE; + + MTaskList tasklist; + parlist.Replace(&tasklist); + if (fit) + parlist.AddToList(fit); + + MFilterList list; + SetupFilters(list, filter); + + MContinue contin(&list); + parlist.AddToList(&list); + + MFillH fill(&hist); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&list); + tasklist.AddToList(&contin); + tasklist.AddToList(&fill); + + // Optimize with the tasklist in this parameterlist + if (!Optimize(parlist)) + return kFALSE; + + // Copy the result back to be accessible by the user + if (fit) + hist.GetAlphaFitter().Copy(*fit); + // Print the result + hist.GetAlphaFitter().Print("result"); + + // Store result if requested + TObjArray cont; + cont.Add(&contin); + return WriteContainer(cont, fNameOut); +} + +//------------------------------------------------------------------------ +// +Bool_t MJOptimizeCuts::RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree) +{ + MHAlpha *histon = CreateNewHist("Hist"); + MHAlpha *histof = CreateNewHist("HistOff"); + + if (!histon || !histof) + return kFALSE; + + const Bool_t rc = RunOnOffCore(*histon, *histof, fname, filter, fit, tree); + + delete histon; + delete histof; + + return rc; +} + +//------------------------------------------------------------------------ +// +Bool_t MJOptimizeCuts::RunOn(const char *fname, MFilter *filter, MAlphaFitter *fit) +{ + MHAlpha *histon = CreateNewHist(); + + if (!histon) + return kFALSE; + + const Bool_t rc = RunOnCore(*histon, fname, filter, fit); + + delete histon; + return rc; +} Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeCuts.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeCuts.h (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeCuts.h (revision 7152) @@ -0,0 +1,57 @@ +#ifndef MARS_MJOptimizeCuts +#define MARS_MJOptimizeCuts + +#ifndef MARS_MJOptimize +#include "MJOptimize.h" +#endif + +class MHAlpha; + +class MJOptimizeCuts : public MJOptimize +{ +private: + TString fNameHist; + + MHAlpha *CreateNewHist(const char *name=0) const; + + Bool_t RunOnCore(MHAlpha &hist, const char *fname, MFilter *filter, MAlphaFitter *fit); + Bool_t RunOnOffCore(MHAlpha &histon, MHAlpha &histoff, const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree); + +public: + MJOptimizeCuts(const char *name="MHAlpha") : MJOptimize(), fNameHist(name) { } + + // Optimization On-Off/Wobble + Bool_t RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit=0, const char *tree="Events"); + Bool_t RunOnOff(const char *fname, MAlphaFitter *fit=0, const char *tree="Events") + { + return RunOnOff(fname, 0, fit, tree); + } + Bool_t RunOnOff(MFilter *filter, MAlphaFitter *fit=0, const char *tree="Events") + { + return RunOnOff(0, filter, fit, tree); + } + Bool_t RunOnOff(MAlphaFitter *fit=0, const char *tree="Events") + { + return RunOnOff(fit, tree); + } + + // Optimization On-only + Bool_t RunOn(const char *fname, MFilter *filter, MAlphaFitter *fit=0); + Bool_t RunOn(const char *fname, MAlphaFitter *fit=0) + { + return RunOn(fname, 0, fit); + } + + Bool_t RunOn(MFilter *filter, MAlphaFitter *fit=0) + { + return RunOn(0, filter, fit); + } + Bool_t RunOn(MAlphaFitter *fit=0) + { + return RunOn(0, 0, fit); + } + + ClassDef(MJOptimizeCuts, 0) // Class for optimization of the Supercuts +}; + +#endif Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.cc (revision 7152) @@ -0,0 +1,192 @@ +/* ======================================================================== *\ +! +! * +! * 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, 6/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ +///////////////////////////////////////////////////////////////////////////// +// +// MJOptimizeDisp +// +// Class for otimizing the disp parameters. For more details see +// MJOptimize. +// +// MJOptimizeDisp minimzes the width of the Theta^2 distribution by +// minimizing (the RMS^2 of Theta^2), while Theta^2 is +// calculated as: +// Theta^2 = d^2 + p^2 - 2*d*p*cos(a) +// with: +// d: MHillasSrc.fDist [deg] +// p: Disp as calculated by the given rule, eg: [0] (a constant) +// a: MHillasSrc.fAlpha [rad] +// +// Example: +// -------- +// MJOptimizeDisp opt; +// opt.SetDebug(2); +// opt.SetOptimizer(MJOptimize::kMigrad); +// opt.SetNumEvents(20000); +// opt.EnableTestTrain(); +// opt.AddParameter("1-(MHillas.fWidth/MHillas.fLength)"); +// opt.SetParameter(0, 1, 0, 2); +// char *r = "[0]*M[0]"; //Rule to calculate disp +// MStatusDisplay *d = new MStatusDisplay; +// opt.SetDisplay(d); +// opt.RunDisp("ganymed-summary.root", r); +// +///////////////////////////////////////////////////////////////////////////// +#include "MJOptimizeDisp.h" + +#include "MHMatrix.h" + +// environment +#include "MLog.h" +#include "MLogManip.h" + +// eventloop +#include "MParList.h" +#include "MTaskList.h" +#include "MEvtLoop.h" + +// parameters +#include "MParameters.h" +#include "MGeomCamMagic.h" + +// histograms +#include "../mhflux/MAlphaFitter.h" +#include "../mhflux/MHThetaSq.h" +#include "../mtools/MChisqEval.h" + +// tasks +#include "MReadTree.h" +#include "MMatrixLoop.h" +#include "MParameterCalc.h" +#include "MFillH.h" + +// filters +#include "MFDataMember.h" + +using namespace std; + +ClassImp(MJOptimizeDisp); + +//------------------------------------------------------------------------ +// +// Read all events from file which do match rules and optimize +// energy estimator. +// +Bool_t MJOptimizeDisp::RunDisp(const char *fname, const char *rule, MTask *weights) +{ + fLog->Separator("Preparing Energy optimization"); + + MParList parlist; + + MParameterI par1("DataType"); + par1.SetVal(1); + parlist.AddToList(&par1); + + MParameterD par2("ThetaSquaredCut"); + par2.SetVal(0.2); + parlist.AddToList(&par2); + + MAlphaFitter fit; + fit.SetPolynomOrder(0); + fit.SetSignalFitMax(0.8); + fit.EnableBackgroundFit(kFALSE); + fit.SetSignalFunction(MAlphaFitter::kThetaSq); + fit.SetMinimizationStrategy(MAlphaFitter::kGaussSigma); + parlist.AddToList(&fit); + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + const Int_t num1 = m.AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg"); + const Int_t num2 = m.AddColumn("MHillasSrc.fAlpha*kDegToRad"); + + MHThetaSq hist; + hist.SkipHistTime(); + hist.SkipHistTheta(); + hist.SkipHistEnergy(); + hist.InitMapping(&m); + + MFDataMember filter("DataType.fVal", '>', 0.5); + fPreCuts.Add(&filter); + + MParameterCalc calc1(rule, "MParameters"); + calc1.SetNameParameter("Disp"); + parlist.AddToList(&calc1); + + const char *disp = "Disp.fVal"; + const char *n1 = Form("M[%d]", num1); + + const char *rule2 = Form("(%s*%s) + (%s*%s) - (2*%s*%s*cos(M[%d]))", + n1, n1, disp, disp, n1, disp, num2); + + MParameterCalc calc2(rule2, "MThetaSqCalc"); + calc2.SetNameParameter("ThetaSquared"); + + MReadTree read("Events"); + // NECESSARY BECAUSE OF MDataFormula GetRules missing + read.DisableAutoScheme(); + if (fname) + read.AddFile(fname); + else + AddSequences(read, fNamesOn); + if (!FillMatrix(read, parlist, kTRUE)) + return kFALSE; + + fPreCuts.Remove(&filter); + + MTaskList tasklist; + parlist.Replace(&tasklist); + + MFillH fill(&hist); + if (weights) + fill.SetWeight(); + + MChisqEval eval; + eval.SetY1("ThetaSquared.fVal"); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&calc1); + tasklist.AddToList(&calc2); + if (weights) + tasklist.AddToList(weights); + tasklist.AddToList(&fill); + tasklist.AddToList(&eval); + + // Optimize with the tasklist in this parameterlist + if (!Optimize(parlist)) + return kFALSE; + + // Print the result + hist.GetAlphaFitter().Print("result"); + + // Store result if requested + TObjArray cont; + cont.Add(&calc1); + return WriteContainer(cont, fNameOut); +} Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.h (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeDisp.h (revision 7152) @@ -0,0 +1,25 @@ +#ifndef MARS_MJOptimizeDisp +#define MARS_MJOptimizeDisp + +#ifndef MARS_MJOptimize +#include "MJOptimize.h" +#endif + +class MTask; + +class MJOptimizeDisp : public MJOptimize +{ +public: + MJOptimizeDisp() : MJOptimize() { } + + // Special optimizing routines + Bool_t RunDisp(const char *fname, const char *rule, MTask *weights=0); + Bool_t RunDisp(const char *rule, MTask *weights=0) + { + return RunDisp(0, rule, weights); + } + + ClassDef(MJOptimizeDisp, 0) // Class for optimization of the Supercuts +}; + +#endif Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.cc =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.cc (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.cc (revision 7152) @@ -0,0 +1,151 @@ +/* ======================================================================== *\ +! +! * +! * 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, 9/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MJOptimize +// +// Class for otimizing a rule to estimate the energy. For more details see +// MJOptimize. +// +// Example: +// -------- +// +// MJOptimizeEnergy opt; +// opt.SetDebug(2); +// opt.SetOptimizer(MJOptimize::kMigrad); +// opt.SetNumEvents(20000); +// opt.EnableTestTrain(); +// opt.AddParameter("MHillas.fSize"); +// opt.SetParameter(0, 1, 0, 2); +// char *r = "[0]*M[0]"; //Rule to calculate estimated energy +// MStatusDisplay *d = new MStatusDisplay; +// opt.SetDisplay(d); +// opt.RunDisp("ganymed-summary.root", r); +// +///////////////////////////////////////////////////////////////////////////// +#include "MJOptimizeEnergy.h" + +#include "MHMatrix.h" + +// environment +#include "MLog.h" +#include "MLogManip.h" + +// eventloop +#include "MParList.h" +#include "MTaskList.h" +#include "MEvtLoop.h" + +// parameters +#include "MParameters.h" +#include "MGeomCamMagic.h" + +// histograms +#include "../mhflux/MHEnergyEst.h" + +// tasks +#include "MReadTree.h" +#include "MMatrixLoop.h" +#include "MEnergyEstimate.h" +#include "MFillH.h" + +// filters +#include "MFDataMember.h" + +using namespace std; + +ClassImp(MJOptimizeEnergy); + +//------------------------------------------------------------------------ +// +// Read all events from file which do match rules and optimize +// energy estimator. +// +Bool_t MJOptimizeEnergy::RunEnergy(const char *fname, const char *rule, MTask *weights) +{ + fLog->Separator("Preparing Energy optimization"); + + MParList parlist; + + MParameterI par("DataType"); + par.SetVal(1); + parlist.AddToList(&par); + + MFDataMember filter("DataType.fVal", '>', 0.5); + fPreCuts.Add(&filter); + + MGeomCamMagic geom; // For GetConvMm2Deg + parlist.AddToList(&geom); + + MHMatrix m("M"); + AddRulesToMatrix(m); + parlist.AddToList(&m); + + MHEnergyEst hist; + hist.InitMapping(&m); + + MParameterCalc est(rule, "MParameters"); + est.SetNameParameter("MEnergyEst"); + parlist.AddToList(&est); + + MReadTree read("Events"); + // NECESSARY BECAUSE OF MDataFormula GetRules missing + read.DisableAutoScheme(); + if (fname) + read.AddFile(fname); + else + AddSequences(read, fNamesOn); + if (!FillMatrix(read, parlist, kTRUE)) + return kFALSE; + + fPreCuts.Remove(&filter); + + MTaskList tasklist; + parlist.Replace(&tasklist); + + MFillH fill(&hist); + if (weights) + fill.SetWeight(); + + MMatrixLoop loop(&m); + + tasklist.AddToList(&loop); + tasklist.AddToList(&est); + if (weights) + tasklist.AddToList(weights); + tasklist.AddToList(&fill); + + // Optimize with the tasklist in this parameterlist + if (!Optimize(parlist)) + return kFALSE; + + // Print the result + hist.Print(); + + // Store result if requested + TObjArray cont; + cont.Add(&est); + return WriteContainer(cont, fNameOut); +} Index: trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.h (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/MJOptimizeEnergy.h (revision 7152) @@ -0,0 +1,25 @@ +#ifndef MARS_MJOptimizeEnergy +#define MARS_MJOptimizeEnergy + +#ifndef MARS_MJOptimize +#include "MJOptimize.h" +#endif + +class MTask; + +class MJOptimizeEnergy : public MJOptimize +{ +public: + MJOptimizeEnergy() : MJOptimize() { } + + // Special optimizing routines + Bool_t RunEnergy(const char *fname, const char *rule, MTask *weights=0); + Bool_t RunEnergy(const char *rule, MTask *weights=0) + { + return RunEnergy(0, rule, weights); + } + + ClassDef(MJOptimizeEnergy, 0) // Class for optimization of the Supercuts +}; + +#endif Index: trunk/MagicSoft/Mars/mjoptim/Makefile =================================================================== --- trunk/MagicSoft/Mars/mjoptim/Makefile (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/Makefile (revision 7152) @@ -0,0 +1,42 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = Optim + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhbase -I../mfbase -I../mdata -I../mfileio \ + -I../manalysis -I../mgeom -I../mfilter -I../mjobs +# mdata: MDataMember +# mfbase: MF +# mfilter: MFMagicCuts +# mhbase: MHMatrix +# manalysis: MMatrixLoop +# mjobs: MJob +# mgeom: MGeomCam + +SRCFILES = MJOptimize.cc \ + MJOptimizeCuts.cc \ + MJOptimizeDisp.cc \ + MJOptimizeEnergy.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak Index: trunk/MagicSoft/Mars/mjoptim/OptimIncl.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/OptimIncl.h (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/OptimIncl.h (revision 7152) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/mjoptim/OptimLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mjoptim/OptimLinkDef.h (revision 7152) +++ trunk/MagicSoft/Mars/mjoptim/OptimLinkDef.h (revision 7152) @@ -0,0 +1,12 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MJOptimize+; +#pragma link C++ class MJOptimizeCuts+; +#pragma link C++ class MJOptimizeDisp+; +#pragma link C++ class MJOptimizeEnergy+; + +#endif Index: trunk/MagicSoft/Mars/mtools/MChisqEval.cc =================================================================== --- trunk/MagicSoft/Mars/mtools/MChisqEval.cc (revision 7151) +++ trunk/MagicSoft/Mars/mtools/MChisqEval.cc (revision 7152) @@ -25,4 +25,7 @@ ///////////////////////////////////////////////////////////////////////////// // +// MChisqEval +// +// Evaluates a chisq from one or two MParameterD calculated in an eventloop. // ///////////////////////////////////////////////////////////////////////////// @@ -30,4 +33,7 @@ #include + +#include "MLog.h" +#include "MLogManip.h" #include "MDataChain.h" @@ -40,14 +46,8 @@ using namespace std; -void MChisqEval::StreamPrimitive(ofstream &out) const -{ - out << " MChisqEval " << GetUniqueName() << ";"; - if (fData0) - out << " " << GetUniqueName() << ".SetY1(\"" << fData0->GetRule() << "\");" << endl; - if (fData1) - out << " " << GetUniqueName() << ".SetY1(\"" << fData1->GetRule() << "\");" << endl; -} +const TString MChisqEval::gsDefName = "MChisqEval"; +const TString MChisqEval::gsDefTitle = "Evaluate a chisq"; -MChisqEval::MChisqEval(const char *name, const char *title) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale) +MChisqEval::MChisqEval(const char *name, const char *title) : fData0(NULL), fData1(NULL), fNameResult("MinimizationValue") { fName = name ? name : gsDefName.Data(); @@ -55,5 +55,5 @@ } -MChisqEval::MChisqEval(MData *y1, const char *name, const char *title) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale) +MChisqEval::MChisqEval(MData *y1, const char *name, const char *title) : fData0(NULL), fData1(NULL), fNameResult("MinimizationValue") { fName = name ? name : gsDefName.Data(); @@ -62,5 +62,5 @@ } -MChisqEval::MChisqEval(MData *y1, MData *y2, const char *name, const char *title) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale) +MChisqEval::MChisqEval(MData *y1, MData *y2, const char *name, const char *title) : fData0(NULL), fData1(NULL), fNameResult("MinimizationValue") { fName = name ? name : gsDefName.Data(); @@ -123,5 +123,5 @@ return kFALSE; - fResult = (MParameterD*)plist->FindCreateObj("MParameterD", "MFitResult"); + fResult = (MParameterD*)plist->FindCreateObj("MParameterD", fNameResult); if (!fResult) return kFALSE; @@ -136,5 +136,5 @@ const Double_t dy = y2-y1; - const Double_t err = fData1 ? y1*y1 : 1; + const Double_t err = fData1 ? y1*y2 : 1; fChisq += dy*dy/err; @@ -144,12 +144,21 @@ Int_t MChisqEval::PostProcess() { - fChisq /= GetNumExecutions(); + if (GetNumExecutions()>0) + fChisq /= GetNumExecutions(); fResult->SetVal(fChisq); + + *fLog << inf << GetDescriptor() << ": Result=" << fChisq << endl; return kTRUE; } -const TString MChisqEval::gsDefName = "MChisqEval"; -const TString MChisqEval::gsDefTitle = "Evaluate a chisq"; +void MChisqEval::StreamPrimitive(ofstream &out) const +{ + out << " MChisqEval " << GetUniqueName() << ";"; + if (fData0) + out << " " << GetUniqueName() << ".SetY1(\"" << fData0->GetRule() << "\");" << endl; + if (fData1) + out << " " << GetUniqueName() << ".SetY1(\"" << fData1->GetRule() << "\");" << endl; +} Index: trunk/MagicSoft/Mars/mtools/MChisqEval.h =================================================================== --- trunk/MagicSoft/Mars/mtools/MChisqEval.h (revision 7151) +++ trunk/MagicSoft/Mars/mtools/MChisqEval.h (revision 7152) @@ -21,10 +21,6 @@ MData *fData1; // Data Member two (measured data) - // -------------------------------------------------------------------------- - // - // Implementation of SavePrimitive. Used to write the call to a constructor - // to a macro. In the original root implementation it is used to write - // gui elements to a macro-file. - // + TString fNameResult; + void StreamPrimitive(ofstream &out) const;