Context Navigation

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

Visit:

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

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

Download in other formats: