Context Navigation

source: tags/Mars-V0.8.2/manalysis/MCT1FindSupercuts.cc

Visit:

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

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

Download in other formats: