Context Navigation

source: trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc@ 2734

Visit:

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

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

Download in other formats: