source: trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc@ 5741

Last change on this file since 5741 was 5671, checked in by domingo, 21 years ago
*** empty log message ***
File size: 31.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): Eva Domingo, 12/2004 <mailto:domingo@ifae.es>
19! Wolfgang Wittek, 12/2004 <mailto:wittek@mppmu.mpg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2005
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27// //
28// MFindDisp //
29// //
30// Class for otimizing the estimation of Disp //
31// //
32// //
33// //
34/////////////////////////////////////////////////////////////////////////////
35#include "MFindDisp.h"
36
37#include <math.h> // fabs
38
39#include <TArrayD.h>
40#include <TCanvas.h>
41#include <TFile.h>
42#include <TH1.h>
43#include <TH2.h>
44#include <TMinuit.h>
45#include <TStopwatch.h>
46#include <TVirtualFitter.h>
47
48#include "MBinning.h"
49#include "MContinue.h"
50#include "MDisp.h"
51#include "MDispCalc.h"
52#include "MDataElement.h"
53#include "MDataMember.h"
54
55#include "MEvtLoop.h"
56#include "MFSelFinal.h"
57#include "MF.h"
58#include "MFDisp.h"
59#include "MFEventSelector.h"
60#include "MFEventSelector2.h"
61#include "MFillH.h"
62#include "MFEventSelector.h"
63#include "MGeomCamMagic.h"
64#include "MH3.h"
65#include "MHDisp.h"
66#include "MHFindSignificance.h"
67#include "MHMatrix.h"
68#include "MHOnSubtraction.h"
69
70#include "MLog.h"
71#include "MLogManip.h"
72#include "MMatrixLoop.h"
73#include "MMinuitInterface.h"
74#include "MParList.h"
75#include "MProgressBar.h"
76#include "MReadMarsFile.h"
77#include "MReadTree.h"
78#include "MStatusDisplay.h"
79#include "MTaskList.h"
80
81
82ClassImp(MFindDisp);
83
84using namespace std;
85
86
87//------------------------------------------------------------------------
88// fcnDisp
89//------------------------------------------------------------------------
90//
91// - calculates the quantity to be minimized (by TMinuit)
92//
93// - the quantity to be minimized is defined in
94// MHDisp::Fill() and MHDisp::Finalize()
95//
96// - the parameters to be varied in the minimization are the parameters
97// appearing in the parametrization of Disp (MDisp::Calc())
98//
99//------------------------------------------------------------------------
100
101static void fcnDisp(Int_t &npar, Double_t *gin, Double_t &f,
102 Double_t *par, Int_t iflag)
103{
104 cout << "entry fcnDisp" << endl;
105
106 // save pointer to the MINUIT oject (for optimizing Disp) for the case
107 // that it is overwritten by the pointer to another Minuit object
108 // TMinuit *savePointer = gMinuit;
109 //-------------------------------------------------------------
110
111
112 MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();
113
114 MParList *plistfcn = (MParList*)evtloopfcn->GetParList();
115 MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList");
116
117 // get needed Disp containers from the existing parList
118 MDisp *disp = (MDisp*)plistfcn->FindObject("MDisp");
119 if (!disp)
120 {
121 gLog << "fcnDisp : MDisp object '" << "MDisp"
122 << "' not found... aborting" << endl;
123 return;
124 }
125
126 MHDisp *hdisp = (MHDisp*)plistfcn->FindObject("MHDisp");
127 if (!hdisp)
128 {
129 gLog << "fcnDisp : MHDisp object '" << "MHDisp"
130 << "' not found... aborting" << endl;
131 return;
132 }
133
134 //
135 // transfer current parameter values to MDisp
136 //
137 // Attention : npar is the number of variable parameters
138 // not the total number of parameters
139 //
140 Double_t fMin, fEdm, fErrdef;
141 Int_t fNpari, fNparx, fIstat;
142 gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
143
144 disp->SetParameters(TArrayD(fNparx, par));
145
146 // checking that parameters have been properly set
147 TArrayD checkparameters = disp->GetParameters();
148 gLog << "fcnDisp : fNpari, fNparx =" << fNpari << ", "
149 << fNparx << endl;
150 gLog << "fcnDisp : i, par, checkparameters =" << endl;
151 for (Int_t i=0; i<fNparx; i++)
152 {
153 gLog << i << ", " << par[i] << ", " << checkparameters[i] << endl;
154 }
155
156 //
157 // loop over the events in the training matrix to compute the Chi^2
158 // for the current parameter values
159 evtloopfcn->Eventloop();
160
161 tasklistfcn->PrintStatistics(0, kTRUE);
162
163 //-------------------------------------------
164 // get the Chi^2 value for the current values of the parameters
165 f = hdisp->GetChi2();
166}
167
168
169// --------------------------------------------------------------------------
170//
171// Default constructor.
172//
173MFindDisp::MFindDisp(MFDisp *fdisp, const char *name, const char *title)
174 : fHowManyTrain(10000), fHowManyTest(10000)
175{
176 fName = name ? name : "MFindDisp";
177 fTitle = title ? title : "Optimizer of Disp";
178
179 //---------------------------
180 // camera geometry is needed for conversion mm ==> degree
181 fCam = new MGeomCamMagic;
182
183 // initialize MFDisp filter cuts to default (no cutting) values
184 // (done at MFDisp constructor)
185 fDispFilter = fdisp;
186 *fLog << inf << "MFindDisp::MFindDisp; address of MFDisp object = "
187 << fDispFilter << endl;
188
189 // matrices to contain the training/test samples
190 fMatrixTrain = new MHMatrix("MatrixTrain");
191 fMatrixTest = new MHMatrix("MatrixTest");
192
193 // objects of MDispCalc to which these matrices are attached
194 fDispCalcTrain = new MDispCalc("DispTrain");
195 fDispCalcTest = new MDispCalc("DispTest");
196
197 // define columns of matrices
198 fDispCalcTrain->InitMapping(fMatrixTrain);
199 fDispCalcTest->InitMapping(fMatrixTest);
200}
201
202// --------------------------------------------------------------------------
203//
204// Default destructor.
205//
206MFindDisp::~MFindDisp()
207{
208 delete fCam;
209 delete fMatrixTrain;
210 delete fMatrixTest;
211 delete fDispCalcTrain;
212 delete fDispCalcTest;
213}
214
215
216// --------------------------------------------------------------------------
217//
218// Define the matrix 'fMatrixTrain' for the TRAINING sample
219//
220// alltogether 'howmanytrain' events are read from file 'nametrain';
221// the events are selected according to a target distribution 'hreftrain'
222//
223//
224Bool_t MFindDisp::DefineTrainMatrix(
225 const TString &nametrain, MH3 &hreftrain,
226 const Int_t howmanytrain, const TString &filetrain,
227 Int_t iflag)
228{
229 if (nametrain.IsNull() || howmanytrain <= 0)
230 return kFALSE;
231
232 *fLog << "=============================================" << endl;
233 *fLog << "Fill TRAINING Matrix from file '" << nametrain << endl;
234 *fLog << " select " << howmanytrain << " events " << endl;
235 if (!hreftrain.GetHist().GetEntries()==0)
236 {
237 *fLog << " according to a distribution given by the MH3 object '"
238 << hreftrain.GetName() << "'" << endl;
239 }
240 else
241 {
242 *fLog << " randomly" << endl;
243 }
244 *fLog << "=============================================" << endl;
245
246
247 MParList plist;
248 MTaskList tlist;
249
250 // initialize display to check the reference distribution
251 // for the event selection (just if iflag==1)
252 MStatusDisplay *display = NULL;
253 if (iflag)
254 display = new MStatusDisplay;
255
256 MReadMarsFile read("Events", nametrain);
257 read.DisableAutoScheme();
258
259 // apply cuts for event selection
260 MContinue contdisp(fDispFilter);
261 contdisp.SetName("ContFilterSelector2");
262
263 // choose a reference distribution
264 MFEventSelector2 seltrain(hreftrain);
265 seltrain.SetNumMax(howmanytrain);
266 seltrain.SetName("selectTrain");
267
268 MFillH filltrain(fMatrixTrain);
269 filltrain.SetFilter(&seltrain);
270 filltrain.SetName("fillMatrixTrain");
271
272 //******************************
273 // entries in MParList
274
275 plist.AddToList(&tlist);
276 plist.AddToList(fCam);
277 plist.AddToList(fMatrixTrain);
278
279 //******************************
280 // entries in MTaskList
281
282 tlist.AddToList(&read);
283 if (fDispFilter != NULL)
284 tlist.AddToList(&contdisp);
285 tlist.AddToList(&seltrain);
286 tlist.AddToList(&filltrain);
287
288 //******************************
289
290 MProgressBar bar;
291 MEvtLoop evtloop;
292 if (display != NULL)
293 evtloop.SetDisplay(display);
294 evtloop.SetParList(&plist);
295 evtloop.SetName("EvtLoopMatrixTrain");
296 evtloop.SetProgressBar(&bar);
297
298 if (!evtloop.Eventloop())
299 return kFALSE;
300
301 tlist.PrintStatistics(0, kTRUE);
302
303
304 // print the filled Training Matrix
305 fMatrixTrain->Print("SizeCols");
306
307 // check that number of generated events is compatible with the resquested
308 Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
309 if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
310 {
311 *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events ("
312 << howmanygenerated
313 << ") is incompatible with the no.of requested events ("
314 << howmanytrain << ")" << endl;
315 }
316
317 *fLog << "TRAINING Matrix was filled" << endl;
318 *fLog << "=============================================" << endl;
319
320
321 //--------------------------
322 // write out training matrix
323
324 if (filetrain != "")
325 {
326 TFile filetr(filetrain, "RECREATE");
327 fMatrixTrain->Write();
328 filetr.Close();
329
330 *fLog << "MFindDisp::DefineTrainMatrix; Training matrix was written onto file '"
331 << filetrain << "'" << endl;
332 }
333
334 if (display != NULL)
335 delete display;
336
337 return kTRUE;
338}
339
340
341// --------------------------------------------------------------------------
342//
343// Define the matrix 'fMatrixTest' for the TEST sample
344//
345// alltogether 'howmanytest' events are read from file 'nametest'
346// the events are selected according to a target distribution 'hreftest'
347//
348//
349Bool_t MFindDisp::DefineTestMatrix(
350 const TString &nametest, MH3 &hreftest,
351 const Int_t howmanytest, const TString &filetest,
352 Int_t iflag)
353{
354 if (nametest.IsNull() || howmanytest<=0)
355 return kFALSE;
356
357 *fLog << "=============================================" << endl;
358 *fLog << "Fill TEST Matrix from file '" << nametest << endl;
359 *fLog << " select " << howmanytest << " events " << endl;
360 if (!hreftest.GetHist().GetEntries()==0)
361 {
362 *fLog << " according to a distribution given by the MH3 object '"
363 << hreftest.GetName() << "'" << endl;
364 }
365 else
366 {
367 *fLog << " randomly" << endl;
368 }
369
370
371 MParList plist;
372 MTaskList tlist;
373
374 // initialize display to check the reference distribution
375 // for the event selection (just if iflag==1)
376 MStatusDisplay *display = NULL;
377 if (iflag)
378 display = new MStatusDisplay;
379
380 MReadMarsFile read("Events", nametest);
381 read.DisableAutoScheme();
382
383 // apply cuts for event selection
384 MContinue contdisp(fDispFilter);
385 contdisp.SetName("ContFilterSelector2");
386
387 // choose a reference distribution
388 MFEventSelector2 seltest(hreftest);
389 seltest.SetNumMax(howmanytest);
390 seltest.SetName("selectTest");
391
392 MFillH filltest(fMatrixTest);
393 filltest.SetFilter(&seltest);
394 filltest.SetName("fillMatrixTest");
395
396 //******************************
397 // entries in MParList
398
399 plist.AddToList(&tlist);
400 plist.AddToList(fCam);
401 plist.AddToList(fMatrixTest);
402
403 //******************************
404 // entries in MTaskList
405
406 tlist.AddToList(&read);
407 if (fDispFilter != NULL)
408 tlist.AddToList(&contdisp);
409 tlist.AddToList(&seltest);
410 tlist.AddToList(&filltest);
411
412 //******************************
413
414 MProgressBar bar;
415 MEvtLoop evtloop;
416 if (display != NULL)
417 evtloop.SetDisplay(display);
418 evtloop.SetParList(&plist);
419 evtloop.SetName("EvtLoopMatrixTest");
420 evtloop.SetProgressBar(&bar);
421
422 if (!evtloop.Eventloop())
423 return kFALSE;
424
425 tlist.PrintStatistics(0, kTRUE);
426
427
428 // print the filled Test Matrix
429 fMatrixTest->Print("SizeCols");
430
431 // check that number of generated events is compatible with the resquested
432 const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
433 if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
434 {
435 *fLog << "MFindDisp::DefineTestMatrix; no.of generated events ("
436 << howmanygenerated
437 << ") is incompatible with the no.of requested events ("
438 << howmanytest << ")" << endl;
439 }
440
441 *fLog << "TEST Matrix was filled" << endl;
442 *fLog << "=============================================" << endl;
443
444
445 //--------------------------
446 // write out test matrix
447
448 if (filetest != "")
449 {
450 TFile filete(filetest, "RECREATE");
451 fMatrixTest->Write();
452 filete.Close();
453
454 *fLog << "MFindDisp::DefineTestMatrix; Test matrix was written onto file '"
455 << filetest << "'" << endl;
456 }
457
458 if (display != NULL)
459 delete display;
460
461 return kTRUE;
462}
463
464
465// --------------------------------------------------------------------------
466//
467// Define the matrices for the TRAINING and TEST sample respectively
468//
469// - this is for the case that training and test samples are generated
470// from the same input file (which has the name 'name')
471//
472Bool_t MFindDisp::DefineTrainTestMatrix(
473 const TString &name, MH3 &href,
474 const Int_t howmanytrain, const Int_t howmanytest,
475 const TString &filetrain, const TString &filetest,
476 Int_t iflag)
477{
478 *fLog << "=============================================" << endl;
479 *fLog << "Fill TRAINING and TEST Matrix from file '" << name << endl;
480 *fLog << " select " << howmanytrain
481 << " training and " << howmanytest << " test events " << endl;
482 if (!href.GetHist().GetEntries()==0)
483 {
484 *fLog << " according to a distribution given by the MH3 object '"
485 << href.GetName() << "'" << endl;
486 }
487 else
488 {
489 *fLog << " randomly" << endl;
490 }
491
492
493 MParList plist;
494 MTaskList tlist;
495
496 // initialize display to check the reference distribution
497 // for the event selection (just if iflag==1)
498 MStatusDisplay *display = NULL;
499 if (iflag)
500 display = new MStatusDisplay;
501
502 MReadMarsFile read("Events", name);
503 read.DisableAutoScheme();
504
505 // apply cuts for event selection
506 MContinue contdisp(fDispFilter);
507 contdisp.SetName("ContFilterSelector2");
508
509 // choose a reference distribution
510 MFEventSelector2 selector(href);
511 selector.SetNumMax(howmanytrain+howmanytest);
512 selector.SetName("selectTrainTest");
513 selector.SetInverted();
514 MContinue cont(&selector);
515 cont.SetName("ContTrainTest");
516
517 // choose randomly the events to fill the Training Matrix
518 Double_t prob = ( (Double_t) howmanytrain )
519 / ( (Double_t)(howmanytrain+howmanytest) );
520 MFEventSelector split;
521 split.SetSelectionRatio(prob);
522
523 MFillH filltrain(fMatrixTrain);
524 filltrain.SetFilter(&split);
525 filltrain.SetName("fillMatrixTrain");
526
527 // consider this event as candidate for a test event
528 // only if event was not accepted as a training event
529 MContinue conttrain(&split);
530 conttrain.SetName("ContTrain");
531
532 MFillH filltest(fMatrixTest);
533 filltest.SetName("fillMatrixTest");
534
535
536 //******************************
537 // entries in MParList
538
539 plist.AddToList(&tlist);
540 plist.AddToList(fCam);
541 plist.AddToList(fMatrixTrain);
542 plist.AddToList(fMatrixTest);
543
544 //******************************
545 // entries in MTaskList
546
547 tlist.AddToList(&read);
548 if (fDispFilter != NULL)
549 tlist.AddToList(&contdisp);
550 tlist.AddToList(&cont);
551 tlist.AddToList(&filltrain);
552 tlist.AddToList(&conttrain);
553 tlist.AddToList(&filltest);
554
555 //******************************
556
557 MProgressBar bar;
558 MEvtLoop evtloop;
559 if (display != NULL)
560 evtloop.SetDisplay(display);
561 evtloop.SetParList(&plist);
562 evtloop.SetName("EvtLoopMatrixTrainTest");
563 evtloop.SetProgressBar(&bar);
564
565 Int_t maxev = -1;
566 if (!evtloop.Eventloop(maxev))
567 return kFALSE;
568
569 tlist.PrintStatistics(0, kTRUE);
570
571
572 // print the filled Train Matrix
573 fMatrixTrain->Print("SizeCols");
574
575 // check that number of generated events is compatible with the resquested
576 const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
577 if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
578 {
579 *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
580 << generatedtrain
581 << ") is incompatible with the no.of requested events ("
582 << howmanytrain << ")" << endl;
583 }
584
585
586 // print the filled Test Matrix
587 fMatrixTest->Print("SizeCols");
588
589 // check that number of generated events is compatible with the resquested
590 const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
591 if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
592 {
593 *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
594 << generatedtest
595 << ") is incompatible with the no.of requested events ("
596 << howmanytest << ")" << endl;
597 }
598
599
600 *fLog << "TRAINING and TEST Matrix were filled" << endl;
601 *fLog << "=============================================" << endl;
602
603
604 //--------------------------
605 // write out training matrix
606
607 if (filetrain != "")
608 {
609 TFile filetr(filetrain, "RECREATE");
610 fMatrixTrain->Write();
611 filetr.Close();
612
613 *fLog << "MFindDisp::DefineTrainTestMatrix; Training matrix was written onto file '"
614 << filetrain << "'" << endl;
615 }
616
617 //--------------------------
618 // write out test matrix
619
620 if (filetest != "")
621 {
622 TFile filete(filetest, "RECREATE");
623 fMatrixTest->Write();
624 filete.Close();
625
626 *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrix was written onto file '"
627 << filetest << "'" << endl;
628 }
629
630 if (display != NULL)
631 delete display;
632
633 return kTRUE;
634}
635
636
637// --------------------------------------------------------------------------
638//
639// Read training and test matrices from files
640//
641//
642Bool_t MFindDisp::ReadMatrix(const TString &filetrain, const TString &filetest)
643{
644 //--------------------------
645 // read in training matrix
646
647 TFile filetr(filetrain);
648 fMatrixTrain->Read("MatrixTrain");
649 fMatrixTrain->Print("SizeCols");
650
651 *fLog << "MFindDisp::ReadMatrix; Training matrix was read in from file '"
652 << filetrain << "'" << endl;
653 filetr.Close();
654
655
656 //--------------------------
657 // read in test matrix
658
659 TFile filete(filetest);
660 fMatrixTest->Read("MatrixTest");
661 fMatrixTest->Print("SizeCols");
662
663 *fLog << "MFindDisp::ReadMatrix; Test matrix was read in from file '"
664 << filetest << "'" << endl;
665 filete.Close();
666
667 return kTRUE;
668}
669
670
671//------------------------------------------------------------------------
672//
673// Steering program for optimizing Disp
674// ------------------------------------
675//
676// the criterion for the 'optimum' is defined in
677// MHDisp::Fill() and MHDisp::Finalize();
678// for example : smallest sum (over all events) of d^2, where d is the
679// distance between the estimated source position
680// (calculated using the current value of Disp) and
681// the true source position
682//
683// The various steps are :
684//
685// - setup the event loop to be executed for each call to fcnDisp
686//
687// - call TMinuit to do the minimization :
688// the fcnDisp function calculates the Chi^2 for the current
689// parameter values;
690// for this - Disp is calculated in the event loop by calling
691// MDispCalc::Process() ==> MDisp::Calc()
692// - the Chi^2 contributions are summed up in the event loop
693// by calling MHDisp::Fill()
694// - after the event loop the Chi^2 is calculated by calling
695// MHDisp::Finalize()
696//
697// Needed as input : (to be set by the Set functions)
698//
699// - fFilenameParam name of file to which optimum values of the
700// parameters are written
701//
702// - for the minimization, the starting values of the parameters are taken
703// - from the file parDispInit (if it is != "")
704// - or from the arrays params and/or steps
705// - or from the Disp constructor
706//
707//----------------------------------------------------------------------
708Bool_t MFindDisp::FindParams(TString parDispInit,
709 TArrayD &params, TArrayD &steps)
710{
711 // Setup the event loop which will be executed in the
712 // fcnDisp function of MINUIT
713 //
714 // parDispInit is the name of the file containing the initial values
715 // of the parameters;
716 // if parDispInit = "" 'params' and 'steps' are taken as initial values
717 //
718
719 *fLog << "=============================================" << endl;
720 *fLog << "Setup event loop for fcnDisp" << endl;
721
722
723 if (fMatrixTrain == NULL)
724 {
725 *fLog << "MFindDisp::FindParams; training matrix is not defined... aborting"
726 << endl;
727 return kFALSE;
728 }
729
730 if (fMatrixTrain->GetM().GetNrows() <= 0)
731 {
732 *fLog << "MFindDisp::FindParams; training matrix has no entries"
733 << endl;
734 return kFALSE;
735 }
736
737 //---------------------------------------------------------
738 MParList parlistfcn;
739 MTaskList tasklistfcn;
740
741 // loop over rows of matrix
742 MMatrixLoop loop(fMatrixTrain);
743
744 //--------------------------------
745 // create container for the Disp parameters
746 // and set them to their initial values
747 MDisp disp;
748
749 // take initial values from file parDispInit
750 if (parDispInit != "")
751 {
752 TFile inparam(parDispInit);
753 disp.Read("MDisp");
754 inparam.Close();
755 *fLog << "MFindDisp::FindParams; initial values of parameters are taken from file "
756 << parDispInit << endl;
757 }
758
759 // take initial values from 'params' and/or 'steps'
760 else if (params.GetSize() != 0 || steps.GetSize() != 0 )
761 {
762 if (params.GetSize() != 0)
763 {
764 *fLog << "MFindDisp::FindParams; initial values of parameters are taken from 'params'"
765 << endl;
766 disp.SetParameters(params);
767 }
768 if (steps.GetSize() != 0)
769 {
770 *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'"
771 << endl;
772 disp.SetStepsizes(steps);
773 }
774 }
775 else
776 {
777 *fLog << "MFindDisp::FindParams; initial values and step sizes are taken "
778 << "from the MDisp constructor" << endl;
779 }
780
781 // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object
782 TArrayI map;
783 MHMatrix *tmpmatrix = NULL;
784 tmpmatrix = fDispCalcTrain->GetMatrixMap(map);
785
786 // attention: argument of MHDisp is the name of MImageParDisp container, that should
787 // be the same than the name given to it when creating MDispCalc object at the MFindDisp
788 // constructor: fDispCalcTrain = new MDispCalc("DispTrain");
789 MHDisp hdisp("DispTrain");
790 hdisp.SetMatrixMap(tmpmatrix,map);
791 // fill the plots for Disp and sum up the Chi^2 contributions
792 MFillH filldispplots("MHDisp", "");
793
794 //******************************
795 // entries in MParList
796
797 parlistfcn.AddToList(&tasklistfcn);
798 parlistfcn.AddToList(&disp);
799 parlistfcn.AddToList(&hdisp);
800 parlistfcn.AddToList(fCam);
801 parlistfcn.AddToList(fMatrixTrain);
802
803 //******************************
804 // entries in MTaskList
805
806 tasklistfcn.AddToList(&loop);
807 tasklistfcn.AddToList(fDispCalcTrain);
808 tasklistfcn.AddToList(&filldispplots);
809
810 //******************************
811
812 MEvtLoop evtloopfcn("EvtLoopFCN");
813 evtloopfcn.SetParList(&parlistfcn);
814 *fLog << "Event loop for fcnDisp has been setup" << endl;
815
816 // address of evtloopfcn is used in CallMinuit()
817
818
819 //-----------------------------------------------------------------------
820 //
821 //---------- Start of minimization part --------------------
822 //
823 // Do the minimization with MINUIT
824 //
825 // Be careful: This is not thread safe
826 //
827 *fLog << "========================================================" << endl;
828 *fLog << "Start minimization for Disp" << endl;
829
830
831 // -------------------------------------------
832 // prepare call to MINUIT
833 //
834
835 // get initial values of parameters
836 fVinit = disp.GetParameters();
837 fStep = disp.GetStepsizes();
838
839 TString name[fVinit.GetSize()];
840 fStep.Set(fVinit.GetSize());
841 fLimlo.Set(fVinit.GetSize());
842 fLimup.Set(fVinit.GetSize());
843 fFix.Set(fVinit.GetSize());
844
845 fNpar = fVinit.GetSize();
846
847 // define names, step sizes, lower and upper limits of the parameters
848 // fFix[] = 0; means the parameter is not fixed
849 for (UInt_t i=0; i<fNpar; i++)
850 {
851 name[i] = "p";
852 name[i] += i+1;
853 //fStep[i] = TMath::Abs(fVinit[i]/10.0);
854 fLimlo[i] = -100.0;
855 fLimup[i] = 100.0;
856 fFix[i] = 0;
857 }
858
859 // fix some parameters
860 //for (UInt_t i=0; i<fNpar; i++)
861 //{
862 // if (i == 1 || i==3)
863 // {
864 // fStep[i] = 0.0;
865 // fFix[i] = 1;
866 // }
867 //}
868
869
870 // -------------------------------------------
871 // call MINUIT
872
873 TStopwatch clock;
874 clock.Start();
875
876 *fLog << "before calling CallMinuit" << endl;
877
878 MMinuitInterface inter;
879 Bool_t rc = inter.CallMinuit(fcnDisp, name,
880 fVinit, fStep, fLimlo, fLimup, fFix,
881 &evtloopfcn, "MIGRAD", kFALSE);
882
883 *fLog << "after calling CallMinuit" << endl;
884
885 *fLog << "Time spent for the minimization in MINUIT : " << endl;;
886 clock.Stop();
887 clock.Print();
888
889
890 if (!rc)
891 return kFALSE;
892
893 *fLog << "Minimization for Disp finished" << endl;
894 *fLog << "========================================================" << endl;
895
896
897 // -----------------------------------------------------------------
898 // in 'fcnDisp' the optimum parameter values were put into MDisp
899
900 // write optimum parameter values onto file fFilenameParam
901
902 TFile outparam(fFilenameParam, "RECREATE");
903 disp.Write();
904 outparam.Close();
905
906 *fLog << "Optimum parameter values for Disp were written onto file '"
907 << fFilenameParam << "' :" << endl;
908
909 const TArrayD &check = disp.GetParameters();
910 for (Int_t i=0; i<check.GetSize(); i++)
911 *fLog << check[i] << ", ";
912 *fLog << endl;
913
914
915
916 //-------------------------------------------
917 // draw the plots
918 hdisp.Draw();
919
920 *fLog << "End of Disp optimization part" << endl;
921 *fLog << "======================================================" << endl;
922
923 return kTRUE;
924}
925
926
927// -----------------------------------------------------------------------
928//
929// Test the result of the Disp optimization on the test sample
930//
931
932Bool_t MFindDisp::TestParams()
933{
934 if (fMatrixTest == NULL)
935 {
936 *fLog << "MFindDisp::TestParams; test matrix is not defined... aborting"
937 << endl;
938 return kFALSE;
939 }
940
941 if (fMatrixTest->GetM().GetNrows() <= 0)
942 {
943 *fLog << "MFindDisp::TestParams; test matrix has no entries"
944 << endl;
945 return kFALSE;
946 }
947
948 // ------------- TEST the Disp optimization -------------------
949 //
950 *fLog << "Test the Disp optimization on the test sample" << endl;
951
952 // -----------------------------------------------------------------
953 // read optimum parameter values from file fFilenameParam
954 // into array 'dispPar'
955
956 TFile inparam(fFilenameParam);
957 MDisp dispin;
958 dispin.Read("MDisp");
959 inparam.Close();
960
961 *fLog << "Optimum parameter values for Disp were read from file '";
962 *fLog << fFilenameParam << "' :" << endl;
963
964 const TArrayD &dispPar = dispin.GetParameters();
965 for (Int_t i=0; i<dispPar.GetSize(); i++)
966 *fLog << dispPar[i] << ", ";
967 *fLog << endl;
968 //---------------------------
969
970
971 // -----------------------------------------------------------------
972 MParList parlist2;
973 MTaskList tasklist2;
974
975 MDisp disp;
976 disp.SetParameters(dispPar);
977
978 MMatrixLoop loop(fMatrixTest);
979
980 // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object
981 TArrayI map;
982 MHMatrix *tmpmatrix = NULL;
983 tmpmatrix = fDispCalcTest->GetMatrixMap(map);
984
985 // attention: argument of MHDisp is the name of MImageParDisp container, that should
986 // be the same than the name given to it when creating MDispCalc object at the MFindDisp
987 // constructor: fDispCalcTrain = new MDispCalc("DispTest");
988 // fill the plots for Disp and sum up the Chi^2 contributions
989 MHDisp hdisp1("DispTest");
990 hdisp1.SetName("MHDispCorr");
991 hdisp1.SetSelectedPos(1);
992 hdisp1.SetMatrixMap(tmpmatrix,map);
993 MFillH filldispplots1("MHDispCorr[MHDisp]", "");
994
995 MHDisp hdisp2("DispTest");
996 hdisp2.SetName("MHDispWrong");
997 hdisp2.SetSelectedPos(2);
998 hdisp2.SetMatrixMap(tmpmatrix,map);
999 MFillH filldispplots2("MHDispWrong[MHDisp]", "");
1000
1001 MHDisp hdisp3("DispTest");
1002 hdisp3.SetName("MHDispM3Long");
1003 hdisp3.SetSelectedPos(3);
1004 hdisp3.SetMatrixMap(tmpmatrix,map);
1005 MFillH filldispplots3("MHDispM3Long[MHDisp]", "");
1006
1007 MHDisp hdisp4("DispTest");
1008 hdisp4.SetName("MHDispAsym");
1009 hdisp4.SetSelectedPos(4);
1010 hdisp4.SetMatrixMap(tmpmatrix,map);
1011 MFillH filldispplots4("MHDispAsym[MHDisp]", "");
1012
1013
1014 //******************************
1015 // entries in MParList
1016
1017 parlist2.AddToList(&tasklist2);
1018 parlist2.AddToList(&disp);
1019 parlist2.AddToList(&hdisp1);
1020 parlist2.AddToList(&hdisp2);
1021 parlist2.AddToList(&hdisp3);
1022 parlist2.AddToList(&hdisp4);
1023 parlist2.AddToList(fCam);
1024 parlist2.AddToList(fMatrixTest);
1025
1026 //******************************
1027 // entries in MTaskList
1028
1029 tasklist2.AddToList(&loop);
1030 tasklist2.AddToList(fDispCalcTest);
1031 tasklist2.AddToList(&filldispplots1);
1032 tasklist2.AddToList(&filldispplots2);
1033 tasklist2.AddToList(&filldispplots3);
1034 tasklist2.AddToList(&filldispplots4);
1035
1036 //******************************
1037
1038 MProgressBar bar2;
1039 MEvtLoop evtloop2;
1040 evtloop2.SetParList(&parlist2);
1041 evtloop2.SetName("EvtLoopTestParams");
1042 evtloop2.SetProgressBar(&bar2);
1043 Int_t maxevents2 = -1;
1044 if (!evtloop2.Eventloop(maxevents2))
1045 return kFALSE;
1046
1047 tasklist2.PrintStatistics(0, kTRUE);
1048
1049
1050 //-------------------------------------------
1051 // draw the plots
1052
1053 hdisp1.Draw();
1054 hdisp2.Draw();
1055 hdisp3.Draw();
1056 hdisp4.Draw();
1057
1058 //-------------------------------------------
1059
1060
1061 *fLog << "Test of Disp otimization finished" << endl;
1062 *fLog << "======================================================" << endl;
1063
1064 return kTRUE;
1065}
1066
1067
1068
1069
1070
1071
1072
1073
Note: See TracBrowser for help on using the repository browser.