Context Navigation

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

Visit:

Last change on this file since 2640 was 2574, checked in by wittek, 21 years ago
* empty log message *
File size: 33.2 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	fUseOrigDistribution = kFALSE;
250
251	//---------------------------
252	// camera geometry is needed for conversion mm ==> degree
253	//fCam = new MGeomCamCT1Daniel;
254	fCam = new MGeomCamCT1;
255
256	// matrices to contain the training/test samples
257	fMatrixTrain = new MHMatrix("MatrixTrain");
258	fMatrixTest = new MHMatrix("MatrixTest");
259
260	// objects of MCT1SupercutsCalc to which these matrices are attached
261	fCalcHadTrain = new MCT1SupercutsCalc("SupercutsCalcTrain");
262	fCalcHadTest = new MCT1SupercutsCalc("SupercutsCalcTest");
263
264	// Define columns of matrices
265	fCalcHadTrain->InitMapping(fMatrixTrain);
266	fCalcHadTest->InitMapping(fMatrixTest);
267	}
268
269	// --------------------------------------------------------------------------
270	//
271	// Default destructor.
272	//
273	MCT1FindSupercuts::~MCT1FindSupercuts()
274	{
275	delete fCam;
276	delete fMatrixTrain;
277	delete fMatrixTest;
278	delete fCalcHadTrain;
279	delete fCalcHadTest;
280	}
281
282
283	// --------------------------------------------------------------------------
284	//
285	//
286	//
287	void MCT1FindSupercuts::SetUseOrigDistribution(Bool_t b)
288	{
289	fUseOrigDistribution = b;
290
291	*fLog << "MCT1FindSupercuts : when defining training and test matrices use the original distribution"
292	<< endl;
293	}
294
295	// --------------------------------------------------------------------------
296	//
297	// Define the matrix 'fMatrixTrain' for the training sample
298	//
299	// alltogether 'howmanytrain' events are read from file 'nametrain';
300	// the events are selected according to a target distribution 'hreftrain'
301	//
302	//
303	Bool_t MCT1FindSupercuts::DefineTrainMatrix(
304	const TString &nametrain, MH3 &hreftrain,
305	const Int_t howmanytrain, const TString &filetrain)
306	{
307	if (nametrain.IsNull() \|\| howmanytrain <= 0)
308	return kFALSE;
309
310	*fLog << "=============================================" << endl;
311	*fLog << "fill training matrix from file '" << nametrain
312	<< "', select " << howmanytrain
313	<< " events " << endl;
314	if (!fUseOrigDistribution)
315	{
316	*fLog << " according to a distribution given by the MH3 object '"
317	<< hreftrain.GetName() << "'" << endl;
318	}
319	else
320	{
321	*fLog << " randomly" << endl;
322	}
323
324
325	MParList plist;
326	MTaskList tlist;
327
328	MReadMarsFile read("Events", nametrain);
329	read.DisableAutoScheme();
330
331	MFEventSelector2 seltrain(hreftrain);
332	seltrain.SetNumMax(howmanytrain);
333	seltrain.SetName("selectTrain");
334	if (fUseOrigDistribution)
335	seltrain.SetUseOrigDistribution(kTRUE);
336
337	MFillH filltrain(fMatrixTrain);
338	filltrain.SetFilter(&seltrain);
339	filltrain.SetName("fillMatrixTrain");
340
341	//******************************
342	// entries in MParList
343
344	plist.AddToList(&tlist);
345	plist.AddToList(fCam);
346	plist.AddToList(fMatrixTrain);
347
348	//******************************
349	// entries in MTaskList
350
351	tlist.AddToList(&read);
352	tlist.AddToList(&seltrain);
353	tlist.AddToList(&filltrain);
354
355	//******************************
356
357	MProgressBar bar;
358	MEvtLoop evtloop;
359	evtloop.SetParList(&plist);
360	evtloop.SetName("EvtLoopMatrixTrain");
361	evtloop.SetProgressBar(&bar);
362
363	if (!evtloop.Eventloop())
364	return kFALSE;
365
366	tlist.PrintStatistics(0, kTRUE);
367
368	fMatrixTrain->Print("SizeCols");
369	Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
370	if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
371	{
372	*fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events ("
373	<< howmanygenerated
374	<< ") is incompatible with the no.of requested events ("
375	<< howmanytrain << ")" << endl;
376	}
377
378	*fLog << "training matrix was filled" << endl;
379	*fLog << "=============================================" << endl;
380
381	//--------------------------
382	// write out training matrix
383
384	if (filetrain != "")
385	{
386	TFile filetr(filetrain, "RECREATE");
387	fMatrixTrain->Write();
388	filetr.Close();
389
390	*fLog << "MCT1FindSupercuts::DefineTrainMatrix; Training matrix was written onto file '"
391	<< filetrain << "'" << endl;
392	}
393
394
395	return kTRUE;
396	}
397
398	// --------------------------------------------------------------------------
399	//
400	// Define the matrix for the test sample
401	//
402	// alltogether 'howmanytest' events are read from file 'nametest'
403	//
404	// the events are selected according to a target distribution 'hreftest'
405	//
406	//
407	Bool_t MCT1FindSupercuts::DefineTestMatrix(
408	const TString &nametest, MH3 &hreftest,
409	const Int_t howmanytest, const TString &filetest)
410	{
411	if (nametest.IsNull() \|\| howmanytest<=0)
412	return kFALSE;
413
414	*fLog << "=============================================" << endl;
415	*fLog << "fill test matrix from file '" << nametest
416	<< "', select " << howmanytest
417	<< " events " << endl;
418	if (!fUseOrigDistribution)
419	{
420	*fLog << " according to a distribution given by the MH3 object '"
421	<< hreftest.GetName() << "'" << endl;
422	}
423	else
424	{
425	*fLog << " randomly" << endl;
426	}
427
428
429	MParList plist;
430	MTaskList tlist;
431
432	MReadMarsFile read("Events", nametest);
433	read.DisableAutoScheme();
434
435	MFEventSelector2 seltest(hreftest);
436	seltest.SetNumMax(howmanytest);
437	seltest.SetName("selectTest");
438	if (fUseOrigDistribution)
439	seltest.SetUseOrigDistribution(kTRUE);
440
441	MFillH filltest(fMatrixTest);
442	filltest.SetFilter(&seltest);
443
444	//******************************
445	// entries in MParList
446
447	plist.AddToList(&tlist);
448	plist.AddToList(fCam);
449	plist.AddToList(fMatrixTest);
450
451	//******************************
452	// entries in MTaskList
453
454	tlist.AddToList(&read);
455	tlist.AddToList(&seltest);
456	tlist.AddToList(&filltest);
457
458	//******************************
459
460	MProgressBar bar;
461	MEvtLoop evtloop;
462	evtloop.SetParList(&plist);
463	evtloop.SetName("EvtLoopMatrixTest");
464	evtloop.SetProgressBar(&bar);
465
466	if (!evtloop.Eventloop())
467	return kFALSE;
468
469	tlist.PrintStatistics(0, kTRUE);
470
471	fMatrixTest->Print("SizeCols");
472	const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
473	if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
474	{
475	*fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events ("
476	<< howmanygenerated
477	<< ") is incompatible with the no.of requested events ("
478	<< howmanytest << ")" << endl;
479	}
480
481	*fLog << "test matrix was filled" << endl;
482	*fLog << "=============================================" << endl;
483
484	//--------------------------
485	// write out test matrix
486
487	if (filetest != "")
488	{
489	TFile filete(filetest, "RECREATE", "");
490	fMatrixTest->Write();
491	filete.Close();
492
493	*fLog << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '"
494	<< filetest << "'" << endl;
495	}
496
497
498	return kTRUE;
499	}
500
501	// --------------------------------------------------------------------------
502	//
503	// Define the matrices for the training and test sample respectively
504	//
505	//
506	//
507	Bool_t MCT1FindSupercuts::DefineTrainTestMatrix(
508	const TString &name, MH3 &href,
509	const Int_t howmanytrain, const Int_t howmanytest,
510	const TString &filetrain, const TString &filetest)
511	{
512	*fLog << "=============================================" << endl;
513	*fLog << "fill training and test matrix from file '" << name
514	<< "', select " << howmanytrain
515	<< " training and " << howmanytest << " test events " << endl;
516	if (!fUseOrigDistribution)
517	{
518	*fLog << " according to a distribution given by the MH3 object '"
519	<< href.GetName() << "'" << endl;
520	}
521	else
522	{
523	*fLog << " randomly" << endl;
524	}
525
526
527	MParList plist;
528	MTaskList tlist;
529
530	MReadMarsFile read("Events", name);
531	read.DisableAutoScheme();
532
533	MFEventSelector2 selector(href);
534	selector.SetNumMax(howmanytrain+howmanytest);
535	selector.SetName("selectTrainTest");
536	selector.SetInverted();
537	if (fUseOrigDistribution)
538	selector.SetUseOrigDistribution(kTRUE);
539
540	MContinue cont(&selector);
541	cont.SetName("ContTrainTest");
542
543	Double_t prob = ( (Double_t) howmanytrain )
544	/ ( (Double_t)(howmanytrain+howmanytest) );
545	MFRandomSplit split(prob);
546
547	MFillH filltrain(fMatrixTrain);
548	filltrain.SetFilter(&split);
549	filltrain.SetName("fillMatrixTrain");
550
551
552	// consider this event as candidate for a test event
553	// only if event was not accepted as a training event
554
555	MContinue conttrain(&split);
556	conttrain.SetName("ContTrain");
557
558	MFillH filltest(fMatrixTest);
559	filltest.SetName("fillMatrixTest");
560
561
562	//******************************
563	// entries in MParList
564
565	plist.AddToList(&tlist);
566	plist.AddToList(fCam);
567	plist.AddToList(fMatrixTrain);
568	plist.AddToList(fMatrixTest);
569
570	//******************************
571	// entries in MTaskList
572
573	tlist.AddToList(&read);
574	tlist.AddToList(&cont);
575
576	tlist.AddToList(&split);
577	tlist.AddToList(&filltrain);
578	tlist.AddToList(&conttrain);
579
580	tlist.AddToList(&filltest);
581
582	//******************************
583
584	MProgressBar bar;
585	MEvtLoop evtloop;
586	evtloop.SetParList(&plist);
587	evtloop.SetName("EvtLoopMatrixTrainTest");
588	evtloop.SetProgressBar(&bar);
589
590	Int_t maxev = -1;
591	if (!evtloop.Eventloop(maxev))
592	return kFALSE;
593
594	tlist.PrintStatistics(0, kTRUE);
595
596	fMatrixTrain->Print("SizeCols");
597	const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
598	if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
599	{
600	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
601	<< generatedtrain
602	<< ") is incompatible with the no.of requested events ("
603	<< howmanytrain << ")" << endl;
604	}
605
606	fMatrixTest->Print("SizeCols");
607	const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
608	if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
609	{
610	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
611	<< generatedtest
612	<< ") is incompatible with the no.of requested events ("
613	<< howmanytest << ")" << endl;
614	}
615
616
617	*fLog << "training and test matrix were filled" << endl;
618	*fLog << "=============================================" << endl;
619
620
621	//--------------------------
622	// write out training matrix
623
624	if (filetrain != "")
625	{
626	TFile filetr(filetrain, "RECREATE");
627	fMatrixTrain->Write();
628	filetr.Close();
629
630	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
631	<< filetrain << "'" << endl;
632	}
633
634	//--------------------------
635	// write out test matrix
636
637	if (filetest != "")
638	{
639	TFile filete(filetest, "RECREATE", "");
640	fMatrixTest->Write();
641	filete.Close();
642
643	*fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
644	<< filetest << "'" << endl;
645	}
646
647	return kTRUE;
648	}
649
650	// --------------------------------------------------------------------------
651	//
652	// Read training and test matrices from files
653	//
654	//
655
656	Bool_t MCT1FindSupercuts::ReadMatrix(const TString &filetrain, const TString &filetest)
657	{
658	//--------------------------
659	// read in training matrix
660
661	TFile filetr(filetrain);
662	fMatrixTrain->Read("MatrixTrain");
663	fMatrixTrain->Print("SizeCols");
664
665	*fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '"
666	<< filetrain << "'" << endl;
667	filetr.Close();
668
669
670	//--------------------------
671	// read in test matrix
672
673	TFile filete(filetest);
674	fMatrixTest->Read("MatrixTest");
675	fMatrixTest->Print("SizeCols");
676
677	*fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '"
678	<< filetest << "'" << endl;
679	filete.Close();
680
681	return kTRUE;
682	}
683
684
685	//------------------------------------------------------------------------
686	//
687	// Steering program for optimizing the supercuts
688	// ---------------------------------------------
689	//
690	// the criterion for the 'optimum' is
691	//
692	// - a maximum significance of the gamma signal in the alpha plot,
693	// in which the supercuts have been applied
694	//
695	// The various steps are :
696	//
697	// - setup the event loop to be executed for each call to fcnSupercuts
698	// - call TMinuit to do the minimization of (-significance) :
699	// the fcnSupercuts function calculates the significance
700	// for the current cut values
701	// for this - the alpha plot is produced in the event loop,
702	// in which the cuts have been applied
703	// - the significance of the gamma signal in the alpha plot
704	// is calculated
705	//
706	// Needed as input : (to be set by the Set functions)
707	//
708	// - fFilenameParam name of file to which optimum values of the
709	// parameters are written
710	//
711	// - fHadronnessName name of container where MCT1SupercutsCalc will
712	// put the supercuts hadronness
713	//
714	// - for the minimization, the starting values of the parameters are taken
715	// - from file parSCinit (if it is != "")
716	// - or from the arrays params and/or steps
717	//
718	//----------------------------------------------------------------------
719	Bool_t MCT1FindSupercuts::FindParams(TString parSCinit,
720	TArrayD &params, TArrayD &steps)
721	{
722	// Setup the event loop which will be executed in the
723	// fcnSupercuts function of MINUIT
724	//
725	// parSCinit is the name of the file containing the initial values
726	// of the parameters;
727	// if parSCinit = "" 'params' and 'steps' are taken as initial values
728
729	*fLog << "=============================================" << endl;
730	*fLog << "Setup event loop for fcnSupercuts" << endl;
731
732	if (fHadronnessName.IsNull())
733	{
734	*fLog << "MCT1FindSupercuts::FindParams; hadronness name is not defined... aborting"
735	<< endl;
736	return kFALSE;
737	}
738
739	if (fMatrixTrain == NULL)
740	{
741	*fLog << "MCT1FindSupercuts::FindParams; training matrix is not defined... aborting"
742	<< endl;
743	return kFALSE;
744	}
745
746	if (fMatrixTrain->GetM().GetNrows() <= 0)
747	{
748	*fLog << "MCT1FindSupercuts::FindParams; training matrix has no entries"
749	<< endl;
750	return kFALSE;
751	}
752
753	MParList parlistfcn;
754	MTaskList tasklistfcn;
755
756	// loop over rows of matrix
757	MMatrixLoop loop(fMatrixTrain);
758
759	//--------------------------------
760	// create container for the supercut parameters
761	// and set them to their initial values
762	MCT1Supercuts super;
763
764	// take initial values from file parSCinit
765	if (parSCinit != "")
766	{
767	TFile inparam(parSCinit);
768	super.Read("MCT1Supercuts");
769	inparam.Close();
770	*fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file "
771	<< parSCinit << endl;
772	}
773
774	// take initial values from 'params' and/or 'steps'
775	else if (params.GetSize() != 0 \|\| steps.GetSize() != 0 )
776	{
777	if (params.GetSize() != 0)
778	{
779	*fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'"
780	<< endl;
781	super.SetParameters(params);
782	}
783	if (steps.GetSize() != 0)
784	{
785	*fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'"
786	<< endl;
787	super.SetStepsizes(steps);
788	}
789	}
790	else
791	{
792	*fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor"
793	<< endl;
794	}
795
796
797	//--------------------------------
798	// calculate supercuts hadronness
799	fCalcHadTrain->SetHadronnessName(fHadronnessName);
800
801	// apply the supercuts
802	MF scfilter(fHadronnessName+".fHadronness>0.5");
803	MContinue supercuts(&scfilter);
804
805	// plot \|alpha\|
806	const TString mh3Name = "AlphaFcn";
807	MBinning binsalpha("Binning"+mh3Name);
808	binsalpha.SetEdges(54, -12.0, 96.0);
809
810	*fLog << warn << "WARNING------------>ALPHA IS ASSUMED TO BE IN COLUMN 7!!!!!!!!" << endl;
811
812	// \|alpha\| is assumed to be in column 7 of the matrix
813	MH3 alpha("MatrixTrain[7]");
814	alpha.SetName(mh3Name);
815
816	MFillH fillalpha(&alpha);
817
818	// book histograms to be filled during the optimization
819	// ! not in the event loop !
820	MHCT1Supercuts plotsuper;
821	plotsuper.SetupFill(&parlistfcn);
822
823	//******************************
824	// entries in MParList (extension of old MParList)
825
826	parlistfcn.AddToList(&tasklistfcn);
827	parlistfcn.AddToList(&super);
828	parlistfcn.AddToList(fCam);
829	parlistfcn.AddToList(fMatrixTrain);
830
831	parlistfcn.AddToList(&binsalpha);
832	parlistfcn.AddToList(&alpha);
833
834	parlistfcn.AddToList(&plotsuper);
835
836	//******************************
837	// entries in MTaskList
838
839	tasklistfcn.AddToList(&loop);
840	tasklistfcn.AddToList(fCalcHadTrain);
841	tasklistfcn.AddToList(&supercuts);
842	tasklistfcn.AddToList(&fillalpha);
843
844
845	//******************************
846
847	MEvtLoop evtloopfcn("EvtLoopFCN");
848	evtloopfcn.SetParList(&parlistfcn);
849	*fLog << "Event loop for fcnSupercuts has been setup" << endl;
850
851	// address of evtloopfcn is used in CallMinuit()
852
853
854	//-----------------------------------------------------------------------
855	//
856	//---------- Start of minimization part --------------------
857	//
858	// Do the minimization with MINUIT
859	//
860	// Be careful: This is not thread safe
861	//
862	*fLog << "========================================================" << endl;
863	*fLog << "Start minimization for supercuts" << endl;
864
865
866	// -------------------------------------------
867	// prepare call to MINUIT
868	//
869
870	// get initial values of parameters
871	fVinit = super.GetParameters();
872	fStep = super.GetStepsizes();
873
874	TString name[fVinit.GetSize()];
875	fStep.Set(fVinit.GetSize());
876	fLimlo.Set(fVinit.GetSize());
877	fLimup.Set(fVinit.GetSize());
878	fFix.Set(fVinit.GetSize());
879
880	fNpar = fVinit.GetSize();
881
882	for (UInt_t i=0; i<fNpar; i++)
883	{
884	name[i] = "p";
885	name[i] += i+1;
886	//fStep[i] = TMath::Abs(fVinit[i]/10.0);
887	fLimlo[i] = -100.0;
888	fLimup[i] = 100.0;
889	fFix[i] = 0;
890	}
891
892	// these parameters make no sense, fix them at 0.0
893	fVinit[33] = 0.0;
894	fStep[33] = 0.0;
895	fFix[33] = 1;
896
897	fVinit[36] = 0.0;
898	fStep[36] = 0.0;
899	fFix[36] = 1;
900
901	fVinit[39] = 0.0;
902	fStep[39] = 0.0;
903	fFix[39] = 1;
904
905	fVinit[41] = 0.0;
906	fStep[41] = 0.0;
907	fFix[41] = 1;
908
909	fVinit[44] = 0.0;
910	fStep[44] = 0.0;
911	fFix[44] = 1;
912
913	fVinit[47] = 0.0;
914	fStep[47] = 0.0;
915	fFix[47] = 1;
916
917	// vary only first 48 parameters
918	//for (UInt_t i=0; i<fNpar; i++)
919	//{
920	// if (i >= 48)
921	// {
922	// fStep[i] = 0.0;
923	// fFix[i] = 1;
924	// }
925	//}
926
927
928	// -------------------------------------------
929	// call MINUIT
930
931	TStopwatch clock;
932	clock.Start();
933
934	*fLog << "before calling CallMinuit" << endl;
935
936	MMinuitInterface inter;
937	Bool_t rc = inter.CallMinuit(fcnSupercuts, name,
938	fVinit, fStep, fLimlo, fLimup, fFix,
939	&evtloopfcn, "SIMPLEX", kFALSE);
940
941	*fLog << "after calling CallMinuit" << endl;
942
943	*fLog << "Time spent for the minimization in MINUIT : " << endl;;
944	clock.Stop();
945	clock.Print();
946
947	plotsuper.DrawClone();
948
949	if (!rc)
950	return kFALSE;
951
952	*fLog << "Minimization for supercuts finished" << endl;
953	*fLog << "========================================================" << endl;
954
955
956	// -----------------------------------------------------------------
957	// in 'fcnSupercuts' (IFLAG=3) the optimum parameter values were put
958	// into MCT1Supercuts
959
960	// write optimum parameter values onto file filenameParam
961
962	TFile outparam(fFilenameParam, "RECREATE");
963	super.Write();
964	outparam.Close();
965
966	*fLog << "Optimum parameter values for supercuts were written onto file '"
967	<< fFilenameParam << "' :" << endl;
968
969	const TArrayD &check = super.GetParameters();
970	for (Int_t i=0; i<check.GetSize(); i++)
971	*fLog << check[i] << ", ";
972	*fLog << endl;
973
974
975
976	*fLog << "End of supercuts optimization part" << endl;
977	*fLog << "======================================================" << endl;
978
979	return kTRUE;
980	}
981
982
983	// -----------------------------------------------------------------------
984	//
985	// Test the supercuts on the test sample
986	//
987
988	Bool_t MCT1FindSupercuts::TestParams()
989	{
990	if (fMatrixTest->GetM().GetNrows() <= 0)
991	{
992	*fLog << "MCT1FindSupercuts::TestParams; test matrix has no entries"
993	<< endl;
994	return kFALSE;
995	}
996
997	// ------------- TEST the supercuts ------------------------------
998	//
999	*fLog << "Test the supercuts on the test sample" << endl;
1000
1001	// -----------------------------------------------------------------
1002	// read optimum parameter values from file filenameParam
1003	// into array 'supercutsPar'
1004
1005	TFile inparam(fFilenameParam);
1006	MCT1Supercuts scin;
1007	scin.Read("MCT1Supercuts");
1008	inparam.Close();
1009
1010	*fLog << "Optimum parameter values for supercuts were read from file '";
1011	*fLog << fFilenameParam << "' :" << endl;
1012
1013	const TArrayD &supercutsPar = scin.GetParameters();
1014	for (Int_t i=0; i<supercutsPar.GetSize(); i++)
1015	*fLog << supercutsPar[i] << ", ";
1016	*fLog << endl;
1017	//---------------------------
1018
1019
1020	// -----------------------------------------------------------------
1021	if (fHadronnessName.IsNull())
1022	{
1023	*fLog << "MCT1FindSupercuts::TestParams; hadronness name is not defined... aborting";
1024	*fLog << endl;
1025	return kFALSE;
1026	}
1027
1028	MParList parlist2;
1029	MTaskList tasklist2;
1030
1031	MCT1Supercuts supercuts;
1032	supercuts.SetParameters(supercutsPar);
1033
1034	fCalcHadTest->SetHadronnessName(fHadronnessName);
1035
1036
1037	//-------------------------------------------
1038
1039	MMatrixLoop loop(fMatrixTest);
1040
1041	// plot alpha before applying the supercuts
1042	const TString mh3NameB = "AlphaBefore";
1043	MBinning binsalphabef("Binning"+mh3NameB);
1044	binsalphabef.SetEdges(54, -12.0, 96.0);
1045
1046	// \|alpha\| is assumed to be in column 7 of the matrix
1047	MH3 alphabefore("MatrixTest[7]");
1048	alphabefore.SetName(mh3NameB);
1049
1050	TH1 &alphahistb = alphabefore.GetHist();
1051	alphahistb.SetName("alphaBefore-TestParams");
1052
1053	MFillH fillalphabefore(&alphabefore);
1054	fillalphabefore.SetName("FillAlphaBefore");
1055
1056	// apply the supercuts
1057	MF scfilter(fHadronnessName+".fHadronness>0.5");
1058	MContinue applysupercuts(&scfilter);
1059
1060	// plot alpha after applying the supercuts
1061	const TString mh3NameA = "AlphaAfter";
1062	MBinning binsalphaaft("Binning"+mh3NameA);
1063	binsalphaaft.SetEdges(54, -12.0, 96.0);
1064
1065	MH3 alphaafter("MatrixTest[7]");
1066	alphaafter.SetName(mh3NameA);
1067
1068	TH1 &alphahista = alphaafter.GetHist();
1069	alphahista.SetName("alphaAfter-TestParams");
1070
1071	MFillH fillalphaafter(&alphaafter);
1072	fillalphaafter.SetName("FillAlphaAfter");
1073
1074	//******************************
1075	// entries in MParList
1076
1077	parlist2.AddToList(&tasklist2);
1078
1079	parlist2.AddToList(&supercuts);
1080
1081	parlist2.AddToList(fCam);
1082	parlist2.AddToList(fMatrixTest);
1083
1084	parlist2.AddToList(&binsalphabef);
1085	parlist2.AddToList(&alphabefore);
1086
1087	parlist2.AddToList(&binsalphaaft);
1088	parlist2.AddToList(&alphaafter);
1089
1090	//******************************
1091	// entries in MTaskList
1092
1093	tasklist2.AddToList(&loop);
1094	tasklist2.AddToList(&fillalphabefore);
1095
1096	tasklist2.AddToList(fCalcHadTest);
1097	tasklist2.AddToList(&applysupercuts);
1098
1099	tasklist2.AddToList(&fillalphaafter);
1100
1101	//******************************
1102
1103	MProgressBar bar2;
1104	MEvtLoop evtloop2;
1105	evtloop2.SetParList(&parlist2);
1106	evtloop2.SetName("EvtLoopTestParams");
1107	evtloop2.SetProgressBar(&bar2);
1108	Int_t maxevents2 = -1;
1109	if (!evtloop2.Eventloop(maxevents2))
1110	return kFALSE;
1111
1112	tasklist2.PrintStatistics(0, kTRUE);
1113
1114
1115	//-------------------------------------------
1116	// draw the alpha plots
1117
1118	MH3* alphaBefore = (MH3*)parlist2.FindObject(mh3NameB, "MH3");
1119	TH1 &alphaHist1 = alphaBefore->GetHist();
1120	alphaHist1.SetXTitle("\|\\alpha\| [\\circ]");
1121
1122	MH3* alphaAfter = (MH3*)parlist2.FindObject(mh3NameA, "MH3");
1123	TH1 &alphaHist2 = alphaAfter->GetHist();
1124	alphaHist2.SetXTitle("\|\\alpha\| [\\circ]");
1125	alphaHist2.SetName("alpha-TestParams");
1126
1127	TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
1128	c->Divide(2,1);
1129
1130	gROOT->SetSelectedPad(NULL);
1131
1132	c->cd(1);
1133	alphaHist1.DrawCopy();
1134
1135	c->cd(2);
1136	alphaHist2.DrawCopy();
1137
1138
1139
1140	//-------------------------------------------
1141	// fit alpha distribution to get the number of excess events and
1142	// calculate significance of gamma signal in the alpha plot
1143
1144	const Double_t alphasig = 20.0;
1145	const Double_t alphamin = 30.0;
1146	const Double_t alphamax = 90.0;
1147	const Int_t degree = 2;
1148	const Bool_t drawpoly = kTRUE;
1149	const Bool_t fitgauss = kTRUE;
1150	const Bool_t print = kTRUE;
1151
1152	MHFindSignificance findsig;
1153	findsig.SetRebin(kTRUE);
1154	findsig.SetReduceDegree(kFALSE);
1155
1156	findsig.FindSigma(&alphaHist2, alphamin, alphamax, degree,
1157	alphasig, drawpoly, fitgauss, print);
1158
1159	const Double_t significance = findsig.GetSignificance();
1160	const Double_t alphasi = findsig.GetAlphasi();
1161
1162	*fLog << "Significance of gamma signal after supercuts in test sample : "
1163	<< significance << " (for \|alpha\| < " << alphasi << " degrees)"
1164	<< endl;
1165	//-------------------------------------------
1166
1167
1168	*fLog << "Test of supercuts part finished" << endl;
1169	*fLog << "======================================================" << endl;
1170
1171	return kTRUE;
1172	}
1173

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

Download in other formats: