Context Navigation

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

Visit:

Last change on this file since 3033 was 2742, checked in by tbretz, 21 years ago
* empty log message *
File size: 32.7 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 "MFEventSelector.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	MFEventSelector split;
525	split.SetSelectionRatio(prob);
526
527	MFillH filltrain(fMatrixTrain);
528	filltrain.SetFilter(&split);
529	filltrain.SetName("fillMatrixTrain");
530
531
532	// consider this event as candidate for a test event
533	// only if event was not accepted as a training event
534
535	MContinue conttrain(&split);
536	conttrain.SetName("ContTrain");
537
538	MFillH filltest(fMatrixTest);
539	filltest.SetName("fillMatrixTest");
540
541
542	//******************************
543	// entries in MParList
544
545	plist.AddToList(&tlist);
546	plist.AddToList(fCam);
547	plist.AddToList(fMatrixTrain);
548	plist.AddToList(fMatrixTest);
549
550	//******************************
551	// entries in MTaskList
552
553	tlist.AddToList(&read);
554	tlist.AddToList(&cont);
555
556	tlist.AddToList(&split);
557	tlist.AddToList(&filltrain);
558	tlist.AddToList(&conttrain);
559
560	tlist.AddToList(&filltest);
561
562	//******************************
563
564	MProgressBar bar;
565	MEvtLoop evtloop;
566	evtloop.SetParList(&plist);
567	evtloop.SetName("EvtLoopMatrixTrainTest");
568	evtloop.SetProgressBar(&bar);
569
570	Int_t maxev = -1;
571	if (!evtloop.Eventloop(maxev))
572	return kFALSE;
573
574	tlist.PrintStatistics(0, kTRUE);
575
576	fMatrixTrain->Print("SizeCols");
577	const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
578	if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
579	{
580	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
581	<< generatedtrain
582	<< ") is incompatible with the no.of requested events ("
583	<< howmanytrain << ")" << endl;
584	}
585
586	fMatrixTest->Print("SizeCols");
587	const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
588	if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
589	{
590	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
591	<< generatedtest
592	<< ") is incompatible with the no.of requested events ("
593	<< howmanytest << ")" << endl;
594	}
595
596
597	*fLog << "training and test matrix were filled" << endl;
598	*fLog << "=============================================" << endl;
599
600
601	//--------------------------
602	// write out training matrix
603
604	if (filetrain != "")
605	{
606	TFile filetr(filetrain, "RECREATE");
607	fMatrixTrain->Write();
608	filetr.Close();
609
610	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
611	<< filetrain << "'" << endl;
612	}
613
614	//--------------------------
615	// write out test matrix
616
617	if (filetest != "")
618	{
619	TFile filete(filetest, "RECREATE", "");
620	fMatrixTest->Write();
621	filete.Close();
622
623	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
624	<< filetest << "'" << endl;
625	}
626
627	return kTRUE;
628	}
629
630	// --------------------------------------------------------------------------
631	//
632	// Read training and test matrices from files
633	//
634	//
635
636	Bool_t MCT1FindSupercuts::ReadMatrix(const TString &filetrain, const TString &filetest)
637	{
638	//--------------------------
639	// read in training matrix
640
641	TFile filetr(filetrain);
642	fMatrixTrain->Read("MatrixTrain");
643	fMatrixTrain->Print("SizeCols");
644
645	*fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '"
646	<< filetrain << "'" << endl;
647	filetr.Close();
648
649
650	//--------------------------
651	// read in test matrix
652
653	TFile filete(filetest);
654	fMatrixTest->Read("MatrixTest");
655	fMatrixTest->Print("SizeCols");
656
657	*fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '"
658	<< filetest << "'" << endl;
659	filete.Close();
660
661	return kTRUE;
662	}
663
664
665	//------------------------------------------------------------------------
666	//
667	// Steering program for optimizing the supercuts
668	// ---------------------------------------------
669	//
670	// the criterion for the 'optimum' is
671	//
672	// - a maximum significance of the gamma signal in the alpha plot,
673	// in which the supercuts have been applied
674	//
675	// The various steps are :
676	//
677	// - setup the event loop to be executed for each call to fcnSupercuts
678	// - call TMinuit to do the minimization of (-significance) :
679	// the fcnSupercuts function calculates the significance
680	// for the current cut values
681	// for this - the alpha plot is produced in the event loop,
682	// in which the cuts have been applied
683	// - the significance of the gamma signal in the alpha plot
684	// is calculated
685	//
686	// Needed as input : (to be set by the Set functions)
687	//
688	// - fFilenameParam name of file to which optimum values of the
689	// parameters are written
690	//
691	// - fHadronnessName name of container where MCT1SupercutsCalc will
692	// put the supercuts hadronness
693	//
694	// - for the minimization, the starting values of the parameters are taken
695	// - from file parSCinit (if it is != "")
696	// - or from the arrays params and/or steps
697	//
698	//----------------------------------------------------------------------
699	Bool_t MCT1FindSupercuts::FindParams(TString parSCinit,
700	TArrayD &params, TArrayD &steps)
701	{
702	// Setup the event loop which will be executed in the
703	// fcnSupercuts function of MINUIT
704	//
705	// parSCinit is the name of the file containing the initial values
706	// of the parameters;
707	// if parSCinit = "" 'params' and 'steps' are taken as initial values
708
709	*fLog << "=============================================" << endl;
710	*fLog << "Setup event loop for fcnSupercuts" << endl;
711
712	if (fHadronnessName.IsNull())
713	{
714	*fLog << "MCT1FindSupercuts::FindParams; hadronness name is not defined... aborting"
715	<< endl;
716	return kFALSE;
717	}
718
719	if (fMatrixTrain == NULL)
720	{
721	*fLog << "MCT1FindSupercuts::FindParams; training matrix is not defined... aborting"
722	<< endl;
723	return kFALSE;
724	}
725
726	if (fMatrixTrain->GetM().GetNrows() <= 0)
727	{
728	*fLog << "MCT1FindSupercuts::FindParams; training matrix has no entries"
729	<< endl;
730	return kFALSE;
731	}
732
733	MParList parlistfcn;
734	MTaskList tasklistfcn;
735
736	// loop over rows of matrix
737	MMatrixLoop loop(fMatrixTrain);
738
739	//--------------------------------
740	// create container for the supercut parameters
741	// and set them to their initial values
742	MCT1Supercuts super;
743
744	// take initial values from file parSCinit
745	if (parSCinit != "")
746	{
747	TFile inparam(parSCinit);
748	super.Read("MCT1Supercuts");
749	inparam.Close();
750	*fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file "
751	<< parSCinit << endl;
752	}
753
754	// take initial values from 'params' and/or 'steps'
755	else if (params.GetSize() != 0 \|\| steps.GetSize() != 0 )
756	{
757	if (params.GetSize() != 0)
758	{
759	*fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'"
760	<< endl;
761	super.SetParameters(params);
762	}
763	if (steps.GetSize() != 0)
764	{
765	*fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'"
766	<< endl;
767	super.SetStepsizes(steps);
768	}
769	}
770	else
771	{
772	*fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor"
773	<< endl;
774	}
775
776
777	//--------------------------------
778	// calculate supercuts hadronness
779	fCalcHadTrain->SetHadronnessName(fHadronnessName);
780
781	// apply the supercuts
782	MF scfilter(fHadronnessName+".fHadronness>0.5");
783	MContinue supercuts(&scfilter);
784
785	// plot \|alpha\|
786	const TString mh3Name = "AlphaFcn";
787	MBinning binsalpha("Binning"+mh3Name);
788	binsalpha.SetEdges(54, -12.0, 96.0);
789
790	*fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;
791
792	// \|alpha\| is assumed to be in column 7 of the matrix
793	MH3 alpha("MatrixTrain[7]");
794	alpha.SetName(mh3Name);
795
796	MFillH fillalpha(&alpha);
797
798	// book histograms to be filled during the optimization
799	// ! not in the event loop !
800	MHCT1Supercuts plotsuper;
801	plotsuper.SetupFill(&parlistfcn);
802
803	//******************************
804	// entries in MParList (extension of old MParList)
805
806	parlistfcn.AddToList(&tasklistfcn);
807	parlistfcn.AddToList(&super);
808	parlistfcn.AddToList(fCam);
809	parlistfcn.AddToList(fMatrixTrain);
810
811	parlistfcn.AddToList(&binsalpha);
812	parlistfcn.AddToList(&alpha);
813
814	parlistfcn.AddToList(&plotsuper);
815
816	//******************************
817	// entries in MTaskList
818
819	tasklistfcn.AddToList(&loop);
820	tasklistfcn.AddToList(fCalcHadTrain);
821	tasklistfcn.AddToList(&supercuts);
822	tasklistfcn.AddToList(&fillalpha);
823
824
825	//******************************
826
827	MEvtLoop evtloopfcn("EvtLoopFCN");
828	evtloopfcn.SetParList(&parlistfcn);
829	*fLog << "Event loop for fcnSupercuts has been setup" << endl;
830
831	// address of evtloopfcn is used in CallMinuit()
832
833
834	//-----------------------------------------------------------------------
835	//
836	//---------- Start of minimization part --------------------
837	//
838	// Do the minimization with MINUIT
839	//
840	// Be careful: This is not thread safe
841	//
842	*fLog << "========================================================" << endl;
843	*fLog << "Start minimization for supercuts" << endl;
844
845
846	// -------------------------------------------
847	// prepare call to MINUIT
848	//
849
850	// get initial values of parameters
851	fVinit = super.GetParameters();
852	fStep = super.GetStepsizes();
853
854	TString name[fVinit.GetSize()];
855	fStep.Set(fVinit.GetSize());
856	fLimlo.Set(fVinit.GetSize());
857	fLimup.Set(fVinit.GetSize());
858	fFix.Set(fVinit.GetSize());
859
860	fNpar = fVinit.GetSize();
861
862	for (UInt_t i=0; i<fNpar; i++)
863	{
864	name[i] = "p";
865	name[i] += i+1;
866	//fStep[i] = TMath::Abs(fVinit[i]/10.0);
867	fLimlo[i] = -100.0;
868	fLimup[i] = 100.0;
869	fFix[i] = 0;
870	}
871
872	// these parameters make no sense, fix them at 0.0
873	fVinit[33] = 0.0;
874	fStep[33] = 0.0;
875	fFix[33] = 1;
876
877	fVinit[36] = 0.0;
878	fStep[36] = 0.0;
879	fFix[36] = 1;
880
881	fVinit[39] = 0.0;
882	fStep[39] = 0.0;
883	fFix[39] = 1;
884
885	fVinit[41] = 0.0;
886	fStep[41] = 0.0;
887	fFix[41] = 1;
888
889	fVinit[44] = 0.0;
890	fStep[44] = 0.0;
891	fFix[44] = 1;
892
893	fVinit[47] = 0.0;
894	fStep[47] = 0.0;
895	fFix[47] = 1;
896
897	// vary only first 48 parameters
898	//for (UInt_t i=0; i<fNpar; i++)
899	//{
900	// if (i >= 48)
901	// {
902	// fStep[i] = 0.0;
903	// fFix[i] = 1;
904	// }
905	//}
906
907
908	// -------------------------------------------
909	// call MINUIT
910
911	TStopwatch clock;
912	clock.Start();
913
914	*fLog << "before calling CallMinuit" << endl;
915
916	MMinuitInterface inter;
917	Bool_t rc = inter.CallMinuit(fcnSupercuts, name,
918	fVinit, fStep, fLimlo, fLimup, fFix,
919	&evtloopfcn, "SIMPLEX", kFALSE);
920
921	*fLog << "after calling CallMinuit" << endl;
922
923	*fLog << "Time spent for the minimization in MINUIT : " << endl;;
924	clock.Stop();
925	clock.Print();
926
927	plotsuper.DrawClone();
928
929	if (!rc)
930	return kFALSE;
931
932	*fLog << "Minimization for supercuts finished" << endl;
933	*fLog << "========================================================" << endl;
934
935
936	// -----------------------------------------------------------------
937	// in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put
938	// into MCT1Supercuts
939
940	// write optimum parameter values onto file filenameParam
941
942	TFile outparam(fFilenameParam, "RECREATE");
943	super.Write();
944	outparam.Close();
945
946	*fLog << "Optimum parameter values for supercuts were written onto file '"
947	<< fFilenameParam << "' :" << endl;
948
949	const TArrayD &check = super.GetParameters();
950	for (Int_t i=0; i<check.GetSize(); i++)
951	*fLog << check[i] << ", ";
952	*fLog << endl;
953
954
955
956	*fLog << "End of supercuts optimization part" << endl;
957	*fLog << "======================================================" << endl;
958
959	return kTRUE;
960	}
961
962
963	// -----------------------------------------------------------------------
964	//
965	// Test the supercuts on the test sample
966	//
967
968	Bool_t MCT1FindSupercuts::TestParams()
969	{
970	if (fMatrixTest->GetM().GetNrows() <= 0)
971	{
972	*fLog << "MCT1FindSupercuts::TestParams; test matrix has no entries"
973	<< endl;
974	return kFALSE;
975	}
976
977	// ------------- TEST the supercuts ------------------------------
978	//
979	*fLog << "Test the supercuts on the test sample" << endl;
980
981	// -----------------------------------------------------------------
982	// read optimum parameter values from file filenameParam
983	// into array 'supercutsPar'
984
985	TFile inparam(fFilenameParam);
986	MCT1Supercuts scin;
987	scin.Read("MCT1Supercuts");
988	inparam.Close();
989
990	*fLog << "Optimum parameter values for supercuts were read from file '";
991	*fLog << fFilenameParam << "' :" << endl;
992
993	const TArrayD &supercutsPar = scin.GetParameters();
994	for (Int_t i=0; i<supercutsPar.GetSize(); i++)
995	*fLog << supercutsPar[i] << ", ";
996	*fLog << endl;
997	//---------------------------
998
999
1000	// -----------------------------------------------------------------
1001	if (fHadronnessName.IsNull())
1002	{
1003	*fLog << "MCT1FindSupercuts::TestParams; hadronness name is not defined... aborting";
1004	*fLog << endl;
1005	return kFALSE;
1006	}
1007
1008	MParList parlist2;
1009	MTaskList tasklist2;
1010
1011	MCT1Supercuts supercuts;
1012	supercuts.SetParameters(supercutsPar);
1013
1014	fCalcHadTest->SetHadronnessName(fHadronnessName);
1015
1016
1017	//-------------------------------------------
1018
1019	MMatrixLoop loop(fMatrixTest);
1020
1021	// plot alpha before applying the supercuts
1022	const TString mh3NameB = "AlphaBefore";
1023	MBinning binsalphabef("Binning"+mh3NameB);
1024	binsalphabef.SetEdges(54, -12.0, 96.0);
1025
1026	// \|alpha\| is assumed to be in column 7 of the matrix
1027	MH3 alphabefore("MatrixTest[7]");
1028	alphabefore.SetName(mh3NameB);
1029
1030	TH1 &alphahistb = alphabefore.GetHist();
1031	alphahistb.SetName("alphaBefore-TestParams");
1032
1033	MFillH fillalphabefore(&alphabefore);
1034	fillalphabefore.SetName("FillAlphaBefore");
1035
1036	// apply the supercuts
1037	MF scfilter(fHadronnessName+".fHadronness>0.5");
1038	MContinue applysupercuts(&scfilter);
1039
1040	// plot alpha after applying the supercuts
1041	const TString mh3NameA = "AlphaAfter";
1042	MBinning binsalphaaft("Binning"+mh3NameA);
1043	binsalphaaft.SetEdges(54, -12.0, 96.0);
1044
1045	MH3 alphaafter("MatrixTest[7]");
1046	alphaafter.SetName(mh3NameA);
1047
1048	TH1 &alphahista = alphaafter.GetHist();
1049	alphahista.SetName("alphaAfter-TestParams");
1050
1051	MFillH fillalphaafter(&alphaafter);
1052	fillalphaafter.SetName("FillAlphaAfter");
1053
1054	//******************************
1055	// entries in MParList
1056
1057	parlist2.AddToList(&tasklist2);
1058
1059	parlist2.AddToList(&supercuts);
1060
1061	parlist2.AddToList(fCam);
1062	parlist2.AddToList(fMatrixTest);
1063
1064	parlist2.AddToList(&binsalphabef);
1065	parlist2.AddToList(&alphabefore);
1066
1067	parlist2.AddToList(&binsalphaaft);
1068	parlist2.AddToList(&alphaafter);
1069
1070	//******************************
1071	// entries in MTaskList
1072
1073	tasklist2.AddToList(&loop);
1074	tasklist2.AddToList(&fillalphabefore);
1075
1076	tasklist2.AddToList(fCalcHadTest);
1077	tasklist2.AddToList(&applysupercuts);
1078
1079	tasklist2.AddToList(&fillalphaafter);
1080
1081	//******************************
1082
1083	MProgressBar bar2;
1084	MEvtLoop evtloop2;
1085	evtloop2.SetParList(&parlist2);
1086	evtloop2.SetName("EvtLoopTestParams");
1087	evtloop2.SetProgressBar(&bar2);
1088	Int_t maxevents2 = -1;
1089	if (!evtloop2.Eventloop(maxevents2))
1090	return kFALSE;
1091
1092	tasklist2.PrintStatistics(0, kTRUE);
1093
1094
1095	//-------------------------------------------
1096	// draw the alpha plots
1097
1098	MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
1099	TH1 &alphaHist1 = alphaBefore->GetHist();
1100	alphaHist1.SetXTitle("\|\\alpha\| [\\circ]");
1101
1102	MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
1103	TH1 &alphaHist2 = alphaAfter->GetHist();
1104	alphaHist2.SetXTitle("\|\\alpha\| [\\circ]");
1105	alphaHist2.SetName("alpha-TestParams");
1106
1107	TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
1108	c->Divide(2,1);
1109
1110	gROOT->SetSelectedPad(NULL);
1111
1112	c->cd(1);
1113	alphaHist1.DrawCopy();
1114
1115	c->cd(2);
1116	alphaHist2.DrawCopy();
1117
1118
1119
1120	//-------------------------------------------
1121	// fit alpha distribution to get the number of excess events and
1122	// calculate significance of gamma signal in the alpha plot
1123
1124	const Double_t alphasig = 20.0;
1125	const Double_t alphamin = 30.0;
1126	const Double_t alphamax = 90.0;
1127	const Int_t degree = 2;
1128	const Bool_t drawpoly = kTRUE;
1129	const Bool_t fitgauss = kTRUE;
1130	const Bool_t print = kTRUE;
1131
1132	MHFindSignificance findsig;
1133	findsig.SetRebin(kTRUE);
1134	findsig.SetReduceDegree(kFALSE);
1135
1136	findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree,
1137	alphasig, drawpoly, fitgauss, print);
1138
1139	const Double_t significance = findsig.GetSignificance();
1140	const Double_t alphasi = findsig.GetAlphasi();
1141
1142	*fLog << "Significance of gamma signal after supercuts in test sample : "
1143	<< significance << " (for \|alpha\| < " << alphasi << " degrees)"
1144	<< endl;
1145	//-------------------------------------------
1146
1147
1148	*fLog << "Test of supercuts part finished" << endl;
1149	*fLog << "======================================================" << endl;
1150
1151	return kTRUE;
1152	}
1153

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

Download in other formats: