Context Navigation

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

Visit:

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

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

Download in other formats: