source: trunk/MagicSoft/Mars/mjobs/MJOptimize.cc@ 6463

Last change on this file since 6463 was 6282, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 27.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 9/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJOptimize
28//
29// Class for otimizing the parameters of the supercuts
30//
31// Minimization Control
32// ====================
33//
34// To choose the minimization algorithm use:
35// void SetOptimizer(Optimizer_t o);
36//
37// Allowed options are:
38// enum Optimizer_t
39// {
40// kMigrad, // Minimize by the method of Migrad
41// kSimplex, // Minimize by the method of Simplex
42// kMinimize, // Migrad + Simplex (if Migrad fails)
43// kMinos, // Minos error determination
44// kImprove, // Local minimum search
45// kSeek, // Minimize by the method of Monte Carlo
46// kNone // Skip optimization
47// };
48//
49// For more details on the methods see TMinuit.
50//
51//
52// You can change the behaviour of the minimization using
53//
54// void SetNumMaxCalls(UInt_t num=0);
55// void SetTolerance(Float_t tol=0);
56//
57// While NumMaxCalls is the first, Tolerance the second arguement.
58// For more details start root and type
59//
60// gMinuit->mnhelp("command")
61//
62// while command can be
63// * MIGRAD
64// * SIMPLEX
65// * MINIMIZE
66// * MINOS
67// * IMPROVE
68// * SEEK
69//
70// The default (num==0 and tol==0) should always give you the
71// corresponding defaults used in Minuit.
72//
73//
74// FIXME: Implement changing cut in hadronness...
75// FIXME: Show MHSignificance on MStatusDisplay during filling...
76// FIXME: Choose step-size percentage as static data membewr
77// FIXME: Choose minimization method
78//
79/////////////////////////////////////////////////////////////////////////////
80#include "MJOptimize.h"
81
82#include <TMinuit.h>
83#include <TVirtualFitter.h>
84
85#include <TStopwatch.h>
86
87#include <TCanvas.h>
88
89#include <TGraph.h>
90#include <TMultiGraph.h>
91
92#include "MParList.h"
93#include "MTaskList.h"
94#include "MGeomCamCT1.h"
95#include "MFEventSelector.h"
96#include "MReadTree.h"
97#include "MHMatrix.h"
98#include "MEnergyEstimate.h"
99#include "MMatrixLoop.h"
100#include "MChisqEval.h"
101#include "MEvtLoop.h"
102#include "MDataElement.h"
103#include "MDataMember.h"
104#include "MLog.h"
105#include "MLogManip.h"
106#include "MParameters.h"
107#include "MFillH.h"
108#include "MF.h"
109#include "MFilterList.h"
110#include "../mfilter/MFSupercuts.h"
111#include "MContinue.h"
112#include "MGeomCamMagic.h"
113#include "../mimage/MHillasSrcCalc.h"
114#include "MHMatrix.h"
115#include "MMatrixLoop.h"
116#include "../mhflux/MHAlpha.h"
117#include "MStatusDisplay.h"
118#include "../mhflux/MHEnergyEst.h"
119#include "MDirIter.h"
120#include "../mjobs/MSequence.h"
121
122using namespace std;
123
124//------------------------------------------------------------------------
125//
126// fcn calculates the function to be minimized (using TMinuit::Migrad)
127//
128void MJOptimize::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
129{
130 MJOptimize *optim = (MJOptimize*)gMinuit->GetObjectFit();
131
132 // WORKAROUND --- FOR WHAT?
133 if (gMinuit->fEpsi<1e-2)
134 {
135 *optim->fLog << warn << "WARNING - For unknown reasons: fEspi<1e-100... resetting to 0.01." << endl;
136 gMinuit->fEpsi = 0.01;
137 }
138
139 TMinuit *minuit = gMinuit;
140 f = optim->Fcn(TArrayD(TMath::Min(gMinuit->fMaxpar, optim->fParameters.GetSize()), par), minuit);
141 gMinuit = minuit;
142
143}
144
145Double_t MJOptimize::Fcn(const TArrayD &par, TMinuit *minuit)
146{
147 if (fEvtLoop->GetDisplay()!=fDisplay)
148 return 0;
149 /*
150 switch(iflag)
151 {
152 case 1: // first call
153 case 2: // calc derivative
154 break;
155 case 3:
156 // last call
157 MStatusDisplay *d = new MStatusDisplay;
158 fEvtLoop->SetDisplay(d);
159 break;
160 }
161 */
162 MParList *plist = fEvtLoop->GetParList();
163
164 MParameterD *eval = (MParameterD*)plist->FindObject("MinimizationValue", "MParameterD");
165 MParContainer *pars = (MParContainer*)plist->FindObject("MParameters", "MParContainer");
166
167 MRead *read = (MRead*)plist->FindObject("MTaskList")->FindObject("MRead");
168 if (read)
169 read->Rewind();
170
171 if (fDebug>=0)
172 {
173 *fLog << inf << "New Set: ";
174 for (Int_t i=0; i<fParameters.GetSize(); i++)
175 *fLog << par[i] << " ";
176 *fLog << endl;
177 }
178
179 pars->SetVariables(par);
180 eval->SetVal(0);
181
182 if (fDebug<3)
183 gLog.SetNullOutput(kTRUE);
184
185 TStopwatch clock;
186 clock.Start();
187 fEvtLoop->Eventloop(fNumEvents);
188 clock.Stop();
189
190 if (fDebug<3)
191 gLog.SetNullOutput(kFALSE);
192
193 const Double_t f = eval->GetVal();
194
195 if (fDebug>=0)
196 *fLog << inf << "Result F=" << f << endl;
197
198 if (fDebug>=1 && minuit)
199 {
200 Double_t fmin, fedm, errdef;
201 Int_t n1, n2, istat;
202 minuit->mnstat(fmin, fedm, errdef, n1, n2, istat);
203 *fLog << inf << underline << "Minimization Status so far:" << endl;
204 *fLog << " Calls: " << minuit->fNfcn << " (max=" << gMinuit->fMaxIterations << ")" << endl;
205 *fLog << " Parameters: fixed=" << gMinuit->fNpfix << ", free=" << gMinuit->fNpar << endl;
206 *fLog << " Func min: " << fmin << " (Epsi=" << gMinuit->fEpsi << ", Apsi=" << gMinuit->fApsi << ")" << endl;
207 *fLog << " Found edm: " << fedm << endl;
208 *fLog << " ErrDef: " << errdef << endl;
209 *fLog << " Status: ";
210
211 switch (istat)
212 {
213 case 0: *fLog << "n/a" << endl; break;
214 case 1: *fLog << "approximation only, not accurate" << endl; break;
215 case 2: *fLog << "full matrix, but forced positive-definite" << endl; break;
216 case 3: *fLog << "full accurate covariance matrix" << endl; break;
217 default: *fLog << "undefined" << endl; break;
218 }
219 }
220
221 if (fDebug>=1)
222 {
223 clock.Print();
224 fEvtLoop->GetTaskList()->PrintStatistics();
225 }
226
227 return f;
228}
229
230MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0)
231{
232 fRules.SetOwner();
233 fFilter.SetOwner();
234
235 fNamesOn.SetOwner();
236 fNamesOff.SetOwner();
237}
238
239//------------------------------------------------------------------------
240//
241// Add seqeunces from list to reader
242//
243Bool_t MJOptimize::AddSequences(MRead &read, TList &list) const
244{
245 MDirIter files;
246
247 TIter Next(&list);
248 TObject *o=0;
249 while ((o=Next()))
250 {
251 MSequence seq(o->GetName());
252 if (!seq.IsValid())
253 return kFALSE;
254
255 seq.SetupDatRuns(files, o->GetTitle(), "I");
256 }
257
258 return read.AddFiles(files)>0;
259}
260
261//------------------------------------------------------------------------
262//
263// Add on-sequences:
264// - fname: sequence file name (with path)
265// - dir: directory were image files are stored
266//
267void MJOptimize::AddSequenceOn(const char *fname, const char *dir)
268{
269 fNamesOn.Add(new TNamed(fname, dir));
270}
271
272//------------------------------------------------------------------------
273//
274// Add off-sequences:
275// - fname: sequence file name (with path)
276// - dir: directory were image files are stored
277//
278void MJOptimize::AddSequenceOff(const char *fname, const char *dir)
279{
280 fNamesOff.Add(new TNamed(fname, dir));
281}
282
283//------------------------------------------------------------------------
284//
285// Empty list of on- and off-sequences
286//
287void MJOptimize::ResetSequences()
288{
289 fNamesOn.Delete();
290 fNamesOff.Delete();
291}
292
293//------------------------------------------------------------------------
294//
295// Add a parameter used in your filters (see AddFilter) The parameter
296// index is returned,
297//
298// Int_t idx = AddParameter("log10(MHillas.fSize)");
299//
300// The indices area starting with 0 always.
301//
302Int_t MJOptimize::AddParameter(const char *rule)
303{
304 fRules.Add(new TNamed(rule, ""));
305 return fRules.GetSize()-1;
306}
307
308//------------------------------------------------------------------------
309//
310// Add a filter which can be applied in the optimization (for deatils
311// see correspodning Run function) You can use the indices you got by
312// AddParameter, eg
313//
314// AddFilter("M[0] < 3.2");
315//
316// if used in optimization you can do
317//
318// AddFilter("M[0] < [0]");
319//
320// for more details, see SetParameter and FixParameter
321//
322void MJOptimize::AddFilter(const char *rule)
323{
324 fFilter.Add(new MF(rule));
325}
326
327//------------------------------------------------------------------------
328//
329// Add a cut which is used to fill the matrix, eg "MMcEvt.fOartId<1.5"
330// (The rule is applied, nit inverted: The matrix is filled with
331// the events fullfilling the condition)
332//
333void MJOptimize::AddPreCut(const char *rule)
334{
335 fPreCuts.Add(new MF(rule));
336}
337
338//------------------------------------------------------------------------
339//
340// Set the fParameters Array accoring to par.
341//
342void MJOptimize::SetParameters(const TArrayD &par)
343{
344 fParameters = par;
345}
346
347//------------------------------------------------------------------------
348//
349// Set the number of events processed by the eventloop. (Be carfull,
350// if you are doing on-off analysis and you only process the first
351// 1000 events which are on-events only the optimization may not work)
352//
353void MJOptimize::SetNumEvents(UInt_t n)
354{
355 fNumEvents=n;
356}
357
358//------------------------------------------------------------------------
359//
360// Set a debug level, which tells the optimization how much information
361// is displayed about and in the running eventloop.
362//
363void MJOptimize::SetDebug(UInt_t n)
364{
365 fDebug=n;
366}
367
368//------------------------------------------------------------------------
369//
370// Set a optimization algorithm to be used. For more information see
371// TMinuit.
372//
373// Available Algorithms are:
374// kMigrad, // Minimize by the method of Migrad
375// kSimplex, // Minimize by the method of Simplex
376// kSeek // Minimize by the method of Monte Carlo
377//
378void MJOptimize::SetOptimizer(Optimizer_t o)
379{
380 fType = o;
381}
382
383//------------------------------------------------------------------------
384//
385// If a status didplay is set, search for tab "Optimizer".
386// If not found, create it.
387// In the tab search for TMultiGraph "Parameters".
388// If not found create it.
389// If empty create TGraphs.
390// Check number of graphs vs. number of parameters.
391// return TList with graphs.
392//
393TList *MJOptimize::GetPlots() const
394{
395 if (!fDisplay)
396 return NULL;
397
398 TCanvas *c=fDisplay->GetCanvas("Optimizer");
399 if (!c)
400 c = &fDisplay->AddTab("Optimizer");
401
402 TMultiGraph *mg = dynamic_cast<TMultiGraph*>(c->FindObject("Parameters"));
403 if (!mg)
404 mg = new TMultiGraph("Parameters", "Parameters of optimization");
405
406 TList *l = mg->GetListOfGraphs();
407 if (!l)
408 {
409 const Int_t n = fParameters.GetSize();
410 for (int i=0; i<n+1; i++)
411 {
412 TGraph *g = new TGraph;
413 if (i==n)
414 g->SetLineColor(kBlue);
415 mg->Add(g, "");
416 AddPoint(mg->GetListOfGraphs(), i, i==n?1:fParameters[i]);
417 }
418 mg->SetBit(kCanDelete);
419 mg->Draw("al*");
420
421 l = mg->GetListOfGraphs();
422 }
423
424 return l->GetSize() == fParameters.GetSize()+1 ? l : NULL;
425}
426
427//------------------------------------------------------------------------
428//
429// Add a point with y=val as last point in idx-th Tgraph of list l.
430//
431void MJOptimize::AddPoint(TList *l, Int_t idx, Float_t val) const
432{
433 if (!l)
434 return;
435
436 TGraph *gr = (TGraph*)l->At(idx);
437 gr->SetPoint(gr->GetN(), gr->GetN(), val);
438}
439
440Int_t MJOptimize::Minuit(TMinuit &minuit, const char *cmd) const
441{
442 Int_t err;
443 Double_t tmp[2] = { fNumMaxCalls, fTolerance };
444 minuit.mnexcm(cmd, tmp, 2, err);
445
446 switch (err)
447 {
448 case 0:
449 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm excuted normally." << endl;
450 break;
451 case 1:
452 *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command is blank... ignored." << endl;
453 break;
454 case 2:
455 *fLog << warn << GetDescriptor() << " TMinuit::mnexcm command-line syntax error... ignored." << endl;
456 break;
457 case 3:
458 *fLog << warn << GetDescriptor() << " TMinuit::mnexcm unknown command... ignored." << endl;
459 break;
460 case 4:
461 *fLog << warn << GetDescriptor() << " TMinuit::mnexcm - Abnormal termination (eg Migrad not converged)" << endl;
462 break;
463 /*
464 case 5:
465 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - Parameters requested." << endl;
466 break;
467 case 6:
468 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET INPUT returned." << endl;
469 break;
470 case 7:
471 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET TITLE returned." << endl;
472 break;
473 case 8:
474 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - SET COVAR returned." << endl;
475 break;
476 case 9:
477 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - reserved." << endl;
478 break;
479 case 10:
480 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - END returned." << endl;
481 break;
482 case 11:
483 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - EXIT or STOP returned." << endl;
484 break;
485 case 12:
486 *fLog << inf << GetDescriptor() << " TMinuit::mnexcm - RETURN returned." << endl;
487 break;*/
488 }
489
490 return err;
491}
492
493Bool_t MJOptimize::Optimize(MEvtLoop &evtloop)
494{
495 if (fType==kNone)
496 return kTRUE;
497
498 gMinuit = new TMinuit(fParameters.GetSize());
499 //gMinuit->SetPrintLevel(-1);
500
501 gMinuit->SetFCN(fcn);
502 gMinuit->SetObjectFit(this);
503
504 //
505 // Set starting values and step sizes for parameters
506 //
507 for (Int_t i=0; i<fParameters.GetSize(); i++)
508 {
509 TString name = "par[";
510 name += i;
511 name += "]";
512 Double_t vinit = fParameters[i];
513 Double_t step = fStep[i];
514
515 Double_t limlo = fLimLo[i];
516 Double_t limup = fLimUp[i];
517
518 Bool_t rc = gMinuit->DefineParameter(i, name, vinit, step, limlo, limup);
519 if (rc)
520 {
521 *fLog << err << dbginf << "Error in defining parameter #" << i << endl;
522 return kFALSE;
523 }
524
525 if (step==0)
526 gMinuit->FixParameter(i);
527 }
528
529 fEvtLoop = &evtloop;
530
531 TList *g=GetPlots();
532
533 // Now ready for minimization step:
534 TStopwatch clock;
535 clock.Start();
536 switch (fType)
537 {
538 case kSimplex:
539 Simplex(*gMinuit);
540 break;
541 case kMigrad:
542 Migrad(*gMinuit);
543 break;
544 case kMinimize:
545 Minimize(*gMinuit);
546 break;
547 case kMinos:
548 Minos(*gMinuit);
549 break;
550 case kImprove:
551 Improve(*gMinuit);
552 break;
553 case kSeek:
554 Seek(*gMinuit);
555 break;
556 case kNone: // Should never happen
557 return kFALSE;
558 }
559 clock.Stop();
560 clock.Print();
561
562 if (evtloop.GetDisplay()!=fDisplay)
563 {
564 *fLog << inf << "Optimization aborted by user." << endl;
565 fDisplay = 0;
566 return kFALSE;
567 }
568
569 *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl;
570
571 //
572 // Update values of fA, fB:
573 //
574 for (Int_t i=0; i<fParameters.GetSize(); i++)
575 {
576 Double_t x1, x2;
577 gMinuit->GetParameter(i,x1,x2);
578 fParameters[i] = x1;
579 cout << i << ": " << fParameters[i] << endl;
580
581 AddPoint(g, i, x1);
582 }
583 AddPoint(g, fParameters.GetSize(), gMinuit->fAmin);
584
585 delete gMinuit;
586
587 return kTRUE;
588}
589
590//------------------------------------------------------------------------
591//
592// Optimize allows to use the optimizing by an eventloop based on
593// some requirements.
594//
595// 1) The tasklist to be executed must have the name MTaskList and
596// be an entry in the parameterlist.
597//
598// 2) The reading task (MReadMarsFile, MMatrixLoop) must have the name
599// "MRead". If it derives from MRead Rewind() must be implemented,
600// otherwise it must start reading from scratch each time its
601// PreProcess is called.
602//
603// 3) The parameters to be optimized must be accessed through (currently)
604// a single parameter container (MParContainer) called MParameters.
605// The parameters are set through SetVariables.
606//
607// 4) The result of a single function call for minimization (eg. chisquare)
608// must be available after the eventloop in a container of type
609// MParameterD with the name "MinimizationResult".
610//
611// 5) The parameters to start with must have been set using
612// MJOptimize::SetParameter or MJOptimize::SetParameters and
613// MJOptimize::FixParameter
614//
615// The behaviour of the optimization can be changed using:
616// void SetNumEvents(UInt_t n);
617// void SetDebug(UInt_t n);
618// void SetOptimizer(Optimizer_t o);
619//
620// After optimization the resulting parameters are set and another eventloop
621// with a MStatusDisplay is set is called. The resulting parameters can be
622// accessed using: GetParameters()
623//
624// To be fixed:
625// - MStatusDisplay should show status while optimization is running
626// - write result into MStatusDisplay
627// - save result
628//
629Bool_t MJOptimize::Optimize(MParList &parlist)
630{
631 // Checks to make sure, that fcn doesn't crash
632 if (!parlist.FindCreateObj("MParameterD", "MinimizationValue"))
633 return kFALSE;
634
635 if (!parlist.FindObject("MParameters", "MParContainer"))
636 {
637 *fLog << err << "MParameters [MParContainer] not found... abort." << endl;
638 return kFALSE;
639 }
640
641 TString txt("Starting ");
642 switch (fType)
643 {
644 case kMigrad: txt += "Migrad"; break;
645 case kMinimize: txt += "Minimize"; break;
646 case kMinos: txt += "Minot"; break;
647 case kImprove: txt += "Improve"; break;
648 case kSimplex: txt += "Simplex"; break;
649 case kSeek: txt += "Seek"; break;
650 case kNone: txt += "no"; break;
651 }
652 txt += " optimization";
653
654 fLog->Separator(txt);
655
656 // Setup eventloop
657 MEvtLoop evtloop;
658 evtloop.SetParList(&parlist);
659 evtloop.SetDisplay(fDisplay); // set display for evtloop and all childs
660 parlist.SetDisplay(0); // reset display for all contents of parlist and tasklist
661 evtloop.SetPrivateDisplay(); // prevent display from being cascaded again in PreProcess
662
663 *fLog << inf << "Number of Parameters: " << fParameters.GetSize() << endl;
664
665 if (!Optimize(evtloop))
666 return kFALSE;
667
668 gMinuit = 0;
669
670 fEvtLoop->SetDisplay(fDisplay);
671 return Fcn(fParameters);
672 // Done already in Fcn
673 // list.SetVariables(fParameters);
674}
675
676void MJOptimize::AddRulesToMatrix(MHMatrix &m) const
677{
678 TIter Next1(&fRules);
679 TObject *o1=0;
680 while ((o1=Next1()))
681 m.AddColumn(o1->GetName());
682}
683
684//------------------------------------------------------------------------
685//
686// Fill matrix m by using read. Use rules as a filter if userules.
687//
688Bool_t MJOptimize::FillMatrix(MReadTree &read, MParList &parlist, Bool_t userules)
689{
690 MHMatrix *m = (MHMatrix*)parlist.FindObject("M", "MHMatrix");
691 if (!m)
692 {
693 *fLog << err << "MJOptimize::FillMatrix - ERROR: M [MHMatrix] not found in parlist... abort." << endl;
694 return kFALSE;
695 }
696
697 m->Print("cols");
698
699 //MParList parlist;
700
701 // MGeomCamMagic cam;
702 // parlist.AddToList(&cam);
703
704 MTaskList tlist;
705 parlist.Replace(&tlist);
706
707 MFillH fillh(m);
708
709 tlist.AddToList(&read);
710
711 MFilterList list;
712 if (!list.AddToList(fPreCuts))
713 *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
714 if (userules)
715 SetupFilters(list);
716 list.SetName("PreCuts"); // reset Name set by SetupFilters
717 list.SetInverted(kFALSE); // reset inversion set by SetupFilters
718 fillh.SetFilter(&list);
719 tlist.AddToList(&list);
720
721 tlist.AddToList(&fillh);
722
723 MEvtLoop fillloop;
724 fillloop.SetParList(&parlist);
725 fillloop.SetDisplay(fDisplay);
726 if (!fillloop.Eventloop(fNumEvents))
727 {
728 *fLog << err << "Filling matrix failed..." << endl;
729 return kFALSE;
730 }
731 tlist.PrintStatistics();
732
733 *fLog << inf << "Read events from file '" << read.GetFileName() << "'" << endl;
734
735 if (fillloop.GetDisplay()!=fDisplay)
736 {
737 fDisplay = 0;
738 *fLog << inf << "Optimization aborted by user." << endl;
739 return kFALSE;
740 }
741
742 m->Print("size");
743
744 return kTRUE;
745}
746
747//------------------------------------------------------------------------
748//
749// Adds all filters to MFilterList
750//
751void MJOptimize::SetupFilters(MFilterList &list, MFilter *filter) const
752{
753 list.SetName("MParameters");
754 list.SetInverted();
755
756 if (filter)
757 {
758 if (fFilter.GetSize()>0)
759 {
760 *fLog << inf;
761 *fLog << "INFORMATION - You are using an external filter and internal filters." << endl;
762 *fLog << " Please make sure that all parameters '[i]' are starting" << endl;
763 *fLog << " behind the number of parameters of the external filter." << endl;
764 }
765 list.AddToList(filter);
766 }
767
768 if (!list.AddToList(fFilter))
769 *fLog << err << "ERROR - Calling MFilterList::AddToList fFilter failed!" << endl;
770
771 *fLog << inf << "Filter: ";
772 list.Print();
773 *fLog << endl;
774}
775
776//------------------------------------------------------------------------
777//
778Bool_t MJOptimize::Run(const char *fname, MFilter *filter, MAlphaFitter *fit)
779{
780 if (fParameters.GetSize()==0)
781 {
782 *fLog << err << "Sorry, no parameters defined." << endl;
783 return kFALSE;
784 }
785
786 fLog->Separator("Preparing On-only-optimization");
787
788 MParList parlist;
789
790 MGeomCamMagic geom; // For GetConvMm2Deg
791 parlist.AddToList(&geom);
792
793 MHMatrix m("M");
794 AddRulesToMatrix(m);
795 parlist.AddToList(&m);
796
797 MHAlpha hist;
798 hist.SkipHistTime();
799 hist.SkipHistTheta();
800 hist.SkipHistEnergy();
801 hist.InitMapping(&m);
802
803 if (filter && filter->InheritsFrom(MFSupercuts::Class()))
804 ((MFSupercuts*)filter)->InitMapping(&m);
805
806 MReadTree read("Events");
807 read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster!
808 if (fname)
809 read.AddFile(fname);
810 else
811 AddSequences(read, fNamesOn);
812 if (!FillMatrix(read, parlist))
813 return kFALSE;
814
815 MTaskList tasklist;
816 parlist.Replace(&tasklist);
817 if (fit)
818 parlist.AddToList(fit);
819
820 MFilterList list;
821 SetupFilters(list, filter);
822
823 MContinue contin(&list);
824 parlist.AddToList(&list);
825
826 MFillH fill(&hist);
827
828 MMatrixLoop loop(&m);
829
830 tasklist.AddToList(&loop);
831 tasklist.AddToList(&list);
832 tasklist.AddToList(&contin);
833 tasklist.AddToList(&fill);
834
835 if (!Optimize(parlist))
836 return kFALSE;
837
838 TObjArray cont;
839 cont.Add(&contin);
840 return WriteContainer(cont, fNameOut);
841}
842
843//------------------------------------------------------------------------
844//
845// Make sure, that filter->GetDataMember is correctly implemented!!!!
846//
847Bool_t MJOptimize::RunOnOff(const char *fname, MFilter *filter, MAlphaFitter *fit, const char *tree)
848{
849 if (fParameters.GetSize()==0)
850 {
851 *fLog << err << "Sorry, no parameters defined." << endl;
852 return kFALSE;
853 }
854
855 fLog->Separator("Preparing On/Off-optimization");
856
857 MParList parlist;
858
859 MGeomCamMagic geom; // For GetConvMm2Deg
860 parlist.AddToList(&geom);
861
862 MHMatrix m("M");
863 AddRulesToMatrix(m);
864 parlist.AddToList(&m);
865
866 const Int_t idxdatatype = m.AddColumn("DataType.fVal");
867
868 MHAlpha histon, histof("MHAlphaOff");
869 histon.SkipHistTime();
870 histon.SkipHistTheta();
871 //histon.SkipHistEnergy();
872 histof.SkipHistTime();
873 histof.SkipHistTheta();
874 //histof.SkipHistEnergy();
875 // FIXME: MHillasSrc.fAlpha is added twice!
876 histon.InitMapping(&m, 1);
877 histof.InitMapping(&m, 1);
878
879 if (filter && filter->InheritsFrom(MFSupercuts::Class()))
880 ((MFSupercuts*)filter)->InitMapping(&m);
881
882 parlist.AddToList(&histon);
883 parlist.AddToList(&histof);
884
885 if (fname)
886 {
887 MReadTree read(tree);
888 read.DisableAutoScheme(); // AutoScheme doesn't seem to be faster!
889 read.AddFile(fname);
890 if (!FillMatrix(read, parlist))
891 return kFALSE;
892 }
893 else
894 {
895 MParameterI par("DataType");
896 parlist.AddToList(&par);
897
898 gLog.Separator("Reading On-Data");
899 par.SetVal(1);
900 MReadTree readon(tree);
901 readon.DisableAutoScheme(); // AutoScheme doesn't seem to be faster!
902 AddSequences(readon, fNamesOn);
903 if (!FillMatrix(readon, parlist))
904 return kFALSE;
905
906 gLog.Separator("Reading Off-Data");
907 par.SetVal(0);
908 MReadTree readoff(tree);
909 readoff.DisableAutoScheme(); // AutoScheme doesn't seem to be faster!
910 AddSequences(readoff, fNamesOff);
911 if (!FillMatrix(readoff, parlist))
912 return kFALSE;
913 }
914
915 MTaskList tasklist;
916 parlist.Replace(&tasklist);
917 if (fit)
918 parlist.AddToList(fit);
919
920 MFilterList list;
921 SetupFilters(list, filter);
922
923 MContinue contin(&list);
924 parlist.AddToList(&list);
925
926 MFillH fillon(&histon, "", "FillHistOn");
927 MFillH fillof(&histof, "", "FillHistOff");
928
929 MF f0(Form("M[%d]<0.5", idxdatatype), "FilterOffData");
930 MF f1(Form("M[%d]>0.5", idxdatatype), "FilterOnData");
931
932 fillof.SetFilter(&f0);
933 fillon.SetFilter(&f1);
934
935 MMatrixLoop loop(&m);
936
937 tasklist.AddToList(&loop);
938 tasklist.AddToList(&list);
939 tasklist.AddToList(&contin);
940 tasklist.AddToList(&f0);
941 tasklist.AddToList(&f1);
942 tasklist.AddToList(&fillof);
943 tasklist.AddToList(&fillon);
944
945 if (!Optimize(parlist))
946 return kFALSE;
947
948 TObjArray cont;
949 cont.Add(&contin);
950 return WriteContainer(cont, fNameOut);
951}
952
953//------------------------------------------------------------------------
954//
955// Read all events from file which do match rules and optimize
956// energy estimator.
957//
958Bool_t MJOptimize::RunEnergy(const char *fname, const char *rule)
959{
960 if (fParameters.GetSize()==0)
961 {
962 *fLog << err << "Sorry, no parameters defined." << endl;
963 return kFALSE;
964 }
965
966 fLog->Separator("Preparing Energy optimization");
967
968 MParList parlist;
969
970 MGeomCamMagic geom; // For GetConvMm2Deg
971 parlist.AddToList(&geom);
972
973 MHMatrix m("M");
974 AddRulesToMatrix(m);
975 parlist.AddToList(&m);
976
977 MEnergyEstimate est("MParameters");
978 est.SetRule(rule);
979 parlist.AddToList(&est);
980
981 MHEnergyEst hist;
982 hist.InitMapping(&m);
983
984 MReadTree read("Events");
985 //read.DisableAutoScheme();
986 if (fname)
987 read.AddFile(fname);
988 else
989 AddSequences(read, fNamesOn);
990 if (!FillMatrix(read, parlist, kTRUE))
991 return kFALSE;
992
993 MTaskList tasklist;
994 parlist.Replace(&tasklist);
995
996 MFillH fill(&hist);
997
998 MMatrixLoop loop(&m);
999
1000 tasklist.AddToList(&loop);
1001 tasklist.AddToList(&est);
1002 tasklist.AddToList(&fill);
1003
1004 if (!Optimize(parlist))
1005 return kFALSE;
1006
1007 TObjArray cont;
1008 cont.Add(&est);
1009 return WriteContainer(cont, fNameOut);
1010}
Note: See TracBrowser for help on using the repository browser.