Context Navigation

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

Visit:

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: